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 };
occgeo->GetIdentifications(perA).push_back(ident);
const auto identOK = occgeo->HaveIdentifications(perA);
// meshing
netgen::MeshingParameters mp;
mp.perfstepsend = netgen::MESHCONST_MESHSURFACE;
auto mesh = make_shared<netgen::Mesh>();
mesh->SetGeometry(occgeo);
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))
continue;
// 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)
continue;
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!
Matthias
(as this is my first post I am obviously not allowed to attach any files)
