How to apply OCCIdentification correctly?

Hi everyone,
I am playing around with netgen for a few days now. I am using it as a library together with OpenCascade in order to get surface and volume meshes from OCCT shapes. This works well for me, but now I got stucked as it comes to periodic meshing. I am using netgen 6.2.2305 and OCCT 7.6.0.

The geometry is a periodic segment of a rotationally symmetric solid. My goal is to get a periodic surface mesh for the segment.

int main(const int argc, const char** argv)
	TopoDS_Shape s;
	BRepTools::Read(s, argv[1], BRep_Builder());
	auto occgeo = make_shared<netgen::OCCGeometry>(s);

	// periodicity: face 1 -> face 26
	const auto angle = 2.0 * M_PI / 7;
	gp_Trsf transf;
	transf.SetRotation(gp::OZ(), angle);
	netgen::Transformation<3> transfNG;
	transfNG.SetAxisRotation(3, angle);

	const int faceIdxA = 1;
	const int faceIdxB = 26;
	const auto perA = occgeo->fmap(faceIdxA);
	const auto perB = occgeo->fmap(faceIdxB);
	netgen::OCCIdentification ident{ perA, perB, transfNG, "periodic", netgen::Identifications::PERIODIC };
	const auto identOK = occgeo->HaveIdentifications(perA);

	// meshing
	netgen::MeshingParameters mp;
	mp.perfstepsend = netgen::MESHCONST_MESHSURFACE;

	auto mesh = make_shared<netgen::Mesh>();

	occgeo->GenerateMesh(mesh, mp);

The generated mesh is not periodic. I check this by exporting the mesh parts related to face 1 and face 26. The transformed mesh of face 1 is not congruent with the one of face 26:

    // result
	for (int fdIdx = 1; fdIdx <= mesh->GetNFD(); fdIdx++)
		// get infos from face descriptor
		netgen::FaceDescriptor fd = mesh->GetFaceDescriptor(fdIdx);
		auto surfNr = fd.SurfNr(); // occ-face index

		if (!(surfNr == faceIdxA || surfNr == faceIdxB))
		// get mesh elements related to face descriptor
		netgen::Array<netgen::SurfaceElementIndex> se_face;
		mesh->GetSurfaceElementsOfFace(fdIdx, se_face);

		std::map<int, int> nodes; // global node index to local node index
		for (auto elemIdx : se_face)
			const netgen::Element2d& elem = mesh->SurfaceElement(elemIdx);
			for (int pt = 0; pt < 3; pt++)
				auto t = nodes.find(elem[pt]);
				if (t == nodes.end())
					nodes.emplace((int)elem[pt], (int)nodes.size() + 1); // one-based index

		// build poly-triangulation
		Handle(Poly_Triangulation) subMesh = new Poly_Triangulation((int)nodes.size(), (int)se_face.Size(), false);

		for (auto node : nodes)
			const netgen::MeshPoint& point = mesh->Point(netgen::PointIndex(node.first));
			subMesh->SetNode(node.second, gp_Pnt(point[0], point[1], point[2]));

		for (int i = 0; i < se_face.Size(); i++)
			const netgen::Element2d& elem = mesh->SurfaceElement(se_face[i]);
			Poly_Triangle tri(nodes[elem[0]], nodes[elem[1]], nodes[elem[2]]);
			subMesh->SetTriangle(1 + i, tri); // one-based index

		// write sub-mesh
		auto pathMesh = "C:\\temp\\mesh_" + std::to_string(surfNr) + ".stl";
		RWStl::WriteAscii(subMesh, OSD_Path(pathMesh.c_str()));

		if (surfNr != faceIdxA)

		auto pathMeshTransf = "C:\\temp\\mesh_" + std::to_string(surfNr) + "_transf.stl";
		for (int nodeIdx = 1; nodeIdx <= subMesh->NbNodes(); nodeIdx++)
			subMesh->SetNode(nodeIdx, subMesh->Node(nodeIdx).Transformed(transf));
		RWStl::WriteAscii(subMesh, OSD_Path(pathMeshTransf.c_str()));

Are my assumptions wrong or what is wrong with the code?

Thank you!

(as this is my first post I am obviously not allowed to attach any files)

Any help from the community is really appreciated here.

I now uploaded the geometry on a file hoster, so anyone can reproduce the case:

I expected that the mesh nodes of the periodic faces being congruent, but they are not:


Hi, the angle needs to be in degree not radian. But then I get an exception in occ egde projection in your code.
Not sure if this is a bug on our or on occ side…

I don’t get any exceptions when running the code. Can you share the code line where the exception occurs?

I corrected the transformation definition using degrees instead of radians now to transfNG.SetAxisRotation(3, angle * 180 / M_PI), but this didn’t change the missing periodicity. The result is still the same.

What me confuses is, that const auto identOK = occgeo->HaveIdentifications(perA); returns true in my code, but when debugging auto & identifications = mesh.GetIdentifications(); is empty inside NetgenGeometry :: FindEdges and have_identifications is false inside NetgenGeometry :: MeshSurface. It seems my defined identification gets lost during processing.

yes netgen throws away non-matching identifications when meshing (as during geometry construction a lot of faces won’t exist any more,…)

from netgen.occ import *
from ngsolve import *

shape = OCCGeometry("s3.brep").shape

for i, face in enumerate(shape.faces): = "Face" + str(i+1)

s1 = shape.faces["Face1"]
s2 = shape.faces["Face26"]
s1.col = (1,0,0)
s2.col = (0,1,0)
trf = gp_Trsf.Rotation(Axis((0,0,0), Z), 360/7)
s1.Identify(s2, "periodic", IdentificationType.PERIODIC, trf)
geo = OCCGeometry(shape)
mesh = Mesh(geo.GenerateMesh(meshsize.moderate))

this is what I did with your brep. Then the face seems to be identified, but in edge identification (netgen also identifies subshapes hierarchically) a GeomAPI_ProjectPointOnCurve::NearestPoint fails…

Okay, since I’m not using Python I cannot reproduce this, at least not with the C++ code I provided.

What I see from your code is, that you initialize OCCGeometry twice, before and after defining the identification. In my code OCCGeometry is initialized only once, before defining the identification. What me confuses here is, that when initializing OCCGeometry the method OCCGeometry :: BuildFMap is called and identifications are processed inside. In my case the periodic identification is not yet defined there. So, do I initialize in a wrong order?

Ah, yes, you need to create a new occgeom object after identification.

Okay, I changed the code accordingly. After defining the identification the OCCGeometry is created again:

TopoDS_Shape s;
BRepTools::Read(s, argv[1], BRep_Builder());
auto occgeoPre = make_shared<netgen::OCCGeometry>(s);
// periodicity: face 1 -> face 26
const auto angle = 2.0 * M_PI / 7;
gp_Trsf transf;
transf.SetRotation(gp::OZ(), angle);
netgen::Transformation<3> transfNG;
transfNG.SetAxisRotation(3, angle * 180 / M_PI); // in degree!

const int faceIdxA = 1;
const int faceIdxB = 26;
const auto perA = occgeoPre->fmap(faceIdxA);
const auto perB = occgeoPre->fmap(faceIdxB);
netgen::OCCIdentification ident{ perA, perB, transfNG, "periodic", netgen::Identifications::PERIODIC };
// build OCCGeometry a second time in order to recognize OCCIdentification correctly
auto occgeo = make_shared<netgen::OCCGeometry>(occgeoPre->shape); // create new OCCGeometry
const auto identOK = occgeo->HaveIdentifications(perA);

That is not working out of the box. When hierarchically adding identifications inside NetgenGeometry :: ProcessIdentifications wrong edge pairs are found at e->identifications.Append( {e, e_other, ident.trafo, ident.type,} );. There is only found 1 edge pair (not as many as as boundary edges exist which I would suspect since periodicity was defined for 1 face). Furthermore, the related edges e and e_other are the same (e == e_other). This naturally continues for point pairs.

Maybe the prior defined identification does not reflect the geometry of the second OCCGeometry instance somehow?

Compare with edge.IsSame(other) not ==, == also checks for orientation (which we do not want here)

I guess the one that is exactly the same is the one in the origin

e and e_other inside NetgenGeometry :: ProcessIdentifications have no IsSame function, do they? They have the same address in memory when debugging.

However, you are correct, they reflect the edge on the z-axis from (0, 0, 50) to (0, 0 ,60). That also means, that this edge geometrically exist in both faces defined as periodic. Maybe this disturbs internal processing of the identifications?

But again, for whatever reason no other edge pairs are found except this wrong one.

For this check do not compare the netgen-Geometryedges, they should differ.
downcast the GeometryEdge to a OCCEdge and then query the TopoDS_Edge inside. This should have a IsSame function

Ah you for sure need a ngsolve upgrade for this. “Cake pieces” up to the origin was something I fixed like a month ago (basically a problem with same edges like you said), it was not possible before, but should work in the newest release.
See here: fix cake identification, allow until revolve axis · NGSolve/netgen@fb211a5 · GitHub

But I still can not run your geometry with newest version, so there is something else…

Edges are found (for me) correctly, it seems, if I add

                 cout << "Edge " << e->nr << " is mapped to edge " << e_other->nr << endl;

in basegeom.cpp line 271 it outputs:

Edge 0 is mapped to edge 26
Edge 1 is mapped to edge 29
Edge 2 is mapped to edge 35
Edge 3 is mapped to edge 37
Edge 4 is mapped to edge 39
Edge 5 is mapped to edge 41
Edge 6 is mapped to edge 43
Edge 7 is mapped to edge 45
Edge 8 is mapped to edge 47
Edge 9 is mapped to edge 49
Edge 10 is mapped to edge 51
Edge 11 is mapped to edge 53
Edge 12 is mapped to edge 55
Edge 13 is mapped to edge 57
Edge 14 is mapped to edge 59
Edge 15 is mapped to edge 61
Edge 16 is mapped to edge 63
Edge 17 is mapped to edge 65
Edge 18 is mapped to edge 67
Edge 19 is mapped to edge 69
Edge 20 is mapped to edge 71
Edge 21 is mapped to edge 73
Edge 22 is mapped to edge 75
Edge 23 is mapped to edge 77
Edge 24 is mapped to edge 24

further inspection yields:

point (0, 0, 60.0352) mapped to (0, 0, 60.0352)
proj: (0, 0, 60.0352) onto edge 77
Traceback (most recent call last):
  File "<string>", line 20, in <module>
netgen.libngpy._NgOCC.OCCException: StdFail_NotDone: GeomAPI_ProjectPointOnCurve::NearestPoint

So it seems occ failes to project the axis point (0,0,60) onto the edge 77, but this point is the edge endpoint. So I’d rather say this is a occ-bug.

Okay, now I can confirm the results you get.

The exception occurs inside OCCEdge::ProjectPoint when trying to retrieve the projection result although the projection failed. Before retrieving result one should check for validity using proj.NbPoints() != 0.

I just returned in case of invalid projection, but this lead to other exceptions in further processing. So I increased the parametric range for projection computation GeomAPI_ProjectPointOnCurve proj(pnt, curve, s0, s1 * 1.1); as a quick test and the projection run successfully. It turned out that the projection result is slightly outside the defined parametric range of the curve (proj.LowerDistanceParameter() returns something like 1.00000003). Furthermore, further processing run successfully as well and I could achieve a periodic mesh.

So, how to correctly handle the projection inside OCCEdge::ProjectPoint? Or in other words: in case of invalid projection why just returning without updating Point<3>& p and EdgePointGeomInfo* gi does not work?

Furthermore, I think Edge 24 is mapped to edge 24 is a bug. It must be mapped to the (geometrically identical) edge inside the other face instead.

Side node: I had to correct the transformation for periodic identification again to transfNG.SetAxisRotation(3, -angle); // in radians and with orientation from perB to perA!.

We have edges only once in the netgen-occ geometry, so mapping 24 to 24 is correct.

I extended the parameter range for the projection now by some small tolerance. This seems to fix this issue, although I’m not sure if this is the occ-ish way to do this.

Fix will be on github soon (runs through our pipelines)

I’m not sure, too. BRepExtrema_DistShapeShape might be another option, it also returns the nearest point and parameter on the curve.

I think this might be more expensive though… I’m a bit hesitant to use this without benchmarking