From 44471fe649c75fb54901110bfe6db0e0655fed15 Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Tue, 14 Dec 2021 12:16:03 +0100 Subject: [PATCH 1/7] occ - identify with given trafo (allows identifying multiple faces at once) --- libsrc/occ/occ_utils.cpp | 16 ++++++++ libsrc/occ/occ_utils.hpp | 3 ++ libsrc/occ/occgeom.cpp | 66 +++++++++++++++++++++++--------- libsrc/occ/occgeom.hpp | 4 +- libsrc/occ/python_occ_shapes.cpp | 2 +- 5 files changed, 70 insertions(+), 21 deletions(-) diff --git a/libsrc/occ/occ_utils.cpp b/libsrc/occ/occ_utils.cpp index 8b747137..210e358e 100644 --- a/libsrc/occ/occ_utils.cpp +++ b/libsrc/occ/occ_utils.cpp @@ -11,6 +11,22 @@ namespace netgen return occ2ng( Handle(BRep_TVertex)::DownCast(shape)->Pnt() ); } + Transformation<3> occ2ng (const gp_Trsf & occ_trafo) + { + Transformation<3> trafo; + auto v = occ_trafo.TranslationPart(); + auto m = occ_trafo.VectorialPart(); + auto & tv = trafo.GetVector(); + auto & tm = trafo.GetMatrix(); + for(auto i : Range(3)) + { + tv[i] = v.Coord(i+1); + for(auto k : Range(3)) + tm(i,k) = m(i+1,k+1); + } + return trafo; + } + Box<3> GetBoundingBox( const TopoDS_Shape & shape ) { Bnd_Box bb; diff --git a/libsrc/occ/occ_utils.hpp b/libsrc/occ/occ_utils.hpp index 1e9fa9a1..e51a2b51 100644 --- a/libsrc/occ/occ_utils.hpp +++ b/libsrc/occ/occ_utils.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include "meshing.hpp" @@ -47,6 +48,8 @@ namespace netgen return occ2ng (BRep_Tool::Pnt (v)); } + DLL_HEADER Transformation<3> occ2ng (const gp_Trsf & t); + inline gp_Pnt ng2occ (const Point<3> & p) { return gp_Pnt(p(0), p(1), p(2)); diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 290c4ece..b013b073 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -1940,12 +1940,19 @@ namespace netgen return occ2ng( props.CentreOfMass() ); } - void OCCGeometry :: IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type) + void OCCGeometry :: IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo) { - auto cme = GetCenter(me); - auto cyou = GetCenter(you); - Transformation<3> trafo{cyou-cme}; - identifications[me.TShape()].push_back( {me.TShape(), you.TShape(), Transformation<3>(cyou - cme), name, type} ); + Transformation<3> trafo; + if(opt_trafo) + trafo = occ2ng(*opt_trafo); + else + { + auto cme = GetCenter(me); + auto cyou = GetCenter(you); + trafo = Transformation<3>{cyou-cme}; + } + + identifications[me.TShape()].push_back( {me.TShape(), you.TShape(), trafo, name, type} ); auto vme = GetVertices(me); auto vyou = GetVertices(you); @@ -1976,6 +1983,11 @@ namespace netgen std::map vmap; auto verts_me = GetVertices(me); + auto verts_you = GetVertices(you); + + if(verts_me.size() != verts_you.size()) + return false; + for (auto i : Range(verts_me.size())) { auto s = verts_me[i].TShape(); @@ -1987,7 +1999,7 @@ namespace netgen } bool all_verts_mapped = true; - for (auto vert : GetVertices(you)) + for (auto vert : verts_you) { auto s = vert.TShape(); auto p = occ2ng(s); @@ -2006,22 +2018,40 @@ namespace netgen return all_verts_mapped; } - void OCCGeometry :: IdentifyFaces(const TopoDS_Shape & solid, const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type) + void OCCGeometry :: IdentifyFaces(const TopoDS_Shape & solid, const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo) { - auto cme = GetCenter(me); - auto cyou = GetCenter(you); - Transformation<3> trafo(cyou-cme); + Transformation<3> trafo; + if(opt_trafo) + { + trafo = occ2ng(*opt_trafo); + } + else + { + auto cme = GetCenter(me); + auto cyou = GetCenter(you); + trafo = Transformation<3>(cyou-cme); + } - identifications[me.TShape()].push_back - (OCCIdentification { me.TShape(), you.TShape(), trafo, name, type }); + auto faces_me = GetFaces(me); + auto faces_you = GetFaces(you); - auto edges_me = GetEdges(me); - auto edges_you = GetEdges(you); + for(auto face_me : faces_me) + for(auto face_you : faces_you) + { + if(!IsMappedShape(trafo, face_me, face_you)) + continue; - for (auto e_me : edges_me) - for (auto e_you : edges_you) - if(IsMappedShape(trafo, e_me, e_you)) - IdentifyEdges(e_me, e_you, name, type); + identifications[face_me.TShape()].push_back + (OCCIdentification { face_me.TShape(), face_you.TShape(), trafo, name, type }); + + auto edges_me = GetEdges(me); + auto edges_you = GetEdges(you); + + for (auto e_me : edges_me) + for (auto e_you : edges_you) + if(IsMappedShape(trafo, e_me, e_you)) + IdentifyEdges(e_me, e_you, name, type, opt_trafo); + } } void OCCParameters :: Print(ostream & ost) const diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index 27eae5b7..790d9727 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -341,8 +341,8 @@ namespace netgen bool ErrorInSurfaceMeshing (); // void WriteOCC_STL(char * filename); - static void IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type); - static void IdentifyFaces(const TopoDS_Shape & solid,const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type); + static void IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo = nullopt); + static void IdentifyFaces(const TopoDS_Shape & solid,const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo = nullopt); private: //bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const; diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index e52d270b..bf4d90ab 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -1091,7 +1091,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) }, py::arg("other"), py::arg("name"), py::arg("type")=Identifications::PERIODIC, "Identify shapes for periodic meshing") .def("Identify", OCCGeometry::IdentifyFaces, "Identify faces", - py::arg("from"), py::arg("to"), py::arg("name"), py::arg("type")=Identifications::PERIODIC) + py::arg("from"), py::arg("to"), py::arg("name"), py::arg("type")=Identifications::PERIODIC, py::arg("trafo")=nullopt) .def("Distance", [](const TopoDS_Shape& self, const TopoDS_Shape& other) From de813df0c25c0e098c93d69b4fd3511686efb77d Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Tue, 14 Dec 2021 12:50:40 +0100 Subject: [PATCH 2/7] add prisms for between closesurface identifications explicitly (no attached faces/edges needed as in prism mesh rules) not active yet (still buggy for CSG) --- libsrc/meshing/meshfunc.cpp | 54 +++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 597a61fa..05f08b72 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -166,6 +166,58 @@ namespace netgen return ret; } + // Add between identified surface elements (only consider closesurface identifications) + void FillCloseSurface( MeshingData & md) + { + static Timer timer("FillCloseSurface"); RegionTimer rtimer(timer); + + auto & mesh = md.mesh; + auto & identifications = mesh->GetIdentifications(); + auto nmax = identifications.GetMaxNr(); + + bool have_closesurfaces = false; + for(auto i : Range(1,nmax+1)) + if(identifications.GetType(i) == Identifications::CLOSESURFACES) + have_closesurfaces = true; + if(!have_closesurfaces) + return; + + NgArray map; + for(auto identnr : Range(1,nmax+1)) + { + if(identifications.GetType(identnr) != Identifications::CLOSESURFACES) + continue; + + identifications.GetMap(identnr, map); + + for(auto & sel : mesh->SurfaceElements()) + { + bool is_mapped = true; + for(auto pi : sel.PNums()) + if(!PointIndex(map[pi]).IsValid()) + is_mapped = false; + + if(!is_mapped) + continue; + + // in case we have symmetric mapping (used in csg), only map in one direction + if(map[map[sel[0]]] == sel[0] && map[sel[0]] < sel[0]) + continue; + + // insert prism + auto np = sel.GetNP(); + Element el(2*np); + for(auto i : Range(np)) + { + el[i] = sel[i]; + el[i+np] = map[sel[i]]; + } + el.SetIndex(md.domain); + mesh->AddVolumeElement(el); + } + } + } + void CloseOpenQuads( MeshingData & md) { auto & mesh = *md.mesh; @@ -488,6 +540,8 @@ namespace netgen if (md[i].mesh->CheckOverlappingBoundary()) throw NgException ("Stop meshing since boundary mesh is overlapping"); + // TODO: FillCloseSurface is still not working with CSG closesurfaces + // FillCloseSurface( md[i] ); CloseOpenQuads( md[i] ); MeshDomain(md[i]); }); From bf261d533f1a0e5f79be556637602c1ab39c28f4 Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Tue, 14 Dec 2021 12:50:59 +0100 Subject: [PATCH 3/7] keep direction of identifications --- libsrc/meshing/meshfunc.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 05f08b72..248ae0c7 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -156,8 +156,6 @@ namespace netgen if(!pi0.IsValid() || !pi1.IsValid()) continue; - if(pi1 Date: Wed, 15 Dec 2021 20:05:18 +0100 Subject: [PATCH 4/7] fix PropagateIdentifications after Glue --- libsrc/meshing/basegeom.cpp | 82 ++++++++++++++++++++++++-- libsrc/meshing/basegeom.hpp | 16 +++++- libsrc/meshing/meshfunc.cpp | 1 + libsrc/occ/occ_edge.cpp | 22 ++----- libsrc/occ/occ_edge.hpp | 8 +-- libsrc/occ/occ_face.cpp | 1 - libsrc/occ/occgeom.cpp | 98 ++++++++++++++------------------ libsrc/occ/occgeom.hpp | 4 +- libsrc/occ/python_occ_shapes.cpp | 75 +++++++++++------------- 9 files changed, 178 insertions(+), 129 deletions(-) diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp index dd108592..4dcbedbc 100644 --- a/libsrc/meshing/basegeom.cpp +++ b/libsrc/meshing/basegeom.cpp @@ -56,18 +56,18 @@ namespace netgen return false; auto & e = *other_ptr; - if(tol < Dist(GetCenter(), e.GetCenter())) + if(tol < Dist(trafo(GetCenter()), e.GetCenter())) return false; - auto v0 = GetStartVertex().GetPoint(); - auto v1 = GetEndVertex().GetPoint(); + auto v0 = trafo(GetStartVertex().GetPoint()); + auto v1 = trafo(GetEndVertex().GetPoint()); auto w0 = e.GetStartVertex().GetPoint(); auto w1 = e.GetEndVertex().GetPoint(); // have two closed edges, use midpoints to compare if(Dist(v0,v1) < tol && Dist(w0,w1) < tol) { - v1 = GetPoint(0.5); + v1 = trafo(GetPoint(0.5)); w1 = other_ptr->GetPoint(0.5); } @@ -75,6 +75,38 @@ namespace netgen (Dist(v0, w1) < tol && Dist(v1, w0) < tol) ); } + bool GeometryFace :: IsMappedShape( const GeometryShape & other_, const Transformation<3> & trafo, double tol ) const + { + const auto other_ptr = dynamic_cast(&other_); + if(!other_ptr) + return false; + auto & f = *other_ptr; + + if(tol < Dist(GetCenter(), f.GetCenter())) + return false; + + // simple check: check if there is a bijective mapping of mapped edges + auto & other_edges = f.edges; + if(edges.Size() != other_edges.Size()) + return false; + + auto nedges = edges.Size(); + Array is_mapped(nedges); + is_mapped = false; + + for(auto e : edges) + { + int found_mapping = 0; + for(auto other_e : other_edges) + if(e->IsMappedShape(*other_e, trafo, tol)) + found_mapping++; + if(found_mapping != 1) + return false; + } + + return true; + } + void GeometryFace :: RestrictHTrig(Mesh& mesh, const PointGeomInfo& gi0, const PointGeomInfo& gi1, @@ -169,6 +201,15 @@ namespace netgen void NetgenGeometry :: ProcessIdentifications() { + for(auto i : Range(vertices)) + vertices[i]->nr = i; + for(auto i : Range(edges)) + edges[i]->nr = i; + for(auto i : Range(faces)) + faces[i]->nr = i; + for(auto i : Range(solids)) + solids[i]->nr = i; + auto mirror_identifications = [&] ( auto & shapes ) { for(auto i : Range(shapes)) @@ -181,11 +222,39 @@ namespace netgen } }; + auto tol = 1e-8 * bounding_box.Diam(); + for(auto & f : faces) + for(auto & ident: f->identifications) + for(auto e : static_cast(ident.from)->edges) + for(auto e_other : static_cast(ident.to)->edges) + if(e->IsMappedShape(*e_other, ident.trafo, tol)) + e->identifications.Append( {e, e_other, ident.trafo, ident.type, ident.name} ); + + for(auto & e : edges) + for(auto & ident: e->identifications) + { + auto & from = static_cast(*ident.from); + auto & to = static_cast(*ident.to); + + GeometryVertex * pfrom[] = { &from.GetStartVertex(), &from.GetEndVertex() }; + GeometryVertex * pto[] = { &to.GetStartVertex(), &to.GetEndVertex() }; + + // swap points of other edge if necessary + Point<3> p_from0 = ident.trafo(from.GetStartVertex().GetPoint()); + Point<3> p_from1 = ident.trafo(from.GetEndVertex().GetPoint()); + Point<3> p_to0 = ident.trafo(to.GetStartVertex().GetPoint()); + + if(Dist(p_from1, p_to0) < Dist(p_from1, p_to0)) + swap(pto[0], pto[1]); + + for(auto i : Range(2)) + pfrom[i]->identifications.Append( {pfrom[i], pto[i], ident.trafo, ident.type, ident.name} ); + } + mirror_identifications(vertices); mirror_identifications(edges); mirror_identifications(faces); - // todo: propagate identifications faces -> edges -> vertices auto find_primary = [&] (auto & shapes) { @@ -729,6 +798,9 @@ namespace netgen sel[1] = s[1]; sel[2] = tree.Find(trafo(mesh[s[1]])); sel[3] = tree.Find(trafo(mesh[s[0]])); + for(auto i : Range(4)) + sel.GeomInfo()[i] = face.Project(mesh[sel[i]]); + sel.SetIndex(face.nr+1); mesh.AddSurfaceElement(sel); } diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index e1e3016b..d7371e2b 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -70,11 +70,19 @@ namespace netgen class DLL_HEADER GeometryEdge : public GeometryShape { + protected: + GeometryVertex *start, *end; public: int domin=-1, domout=-1; - virtual const GeometryVertex& GetStartVertex() const = 0; - virtual const GeometryVertex& GetEndVertex() const = 0; + GeometryEdge( GeometryVertex &start_, GeometryVertex &end_ ) + : start(&start_), end(&end_) + {} + + virtual const GeometryVertex& GetStartVertex() const { return *start; } + virtual const GeometryVertex& GetEndVertex() const { return *end; } + virtual GeometryVertex& GetStartVertex() { return *start; } + virtual GeometryVertex& GetEndVertex() { return *end; } virtual double GetLength() const = 0; virtual Point<3> GetCenter() const = 0; virtual Point<3> GetPoint(double t) const = 0; @@ -100,11 +108,13 @@ namespace netgen class DLL_HEADER GeometryFace : public GeometryShape { public: + Array edges; int domin=-1, domout=-1; virtual Point<3> GetCenter() const = 0; virtual size_t GetNBoundaries() const = 0; virtual Array GetBoundary(const Mesh& mesh) const = 0; + virtual PointGeomInfo Project(Point<3>& p) const = 0; // Project point using geo info. Fast if point is close to // parametrization in geo info. @@ -143,6 +153,8 @@ namespace netgen newgi = Project(newp); } + virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override; + protected: void RestrictHTrig(Mesh& mesh, const PointGeomInfo& gi0, diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 248ae0c7..6714ccf8 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -210,6 +210,7 @@ namespace netgen el[i] = sel[i]; el[i+np] = map[sel[i]]; } + el.SetIndex(md.domain); mesh->AddVolumeElement(el); } diff --git a/libsrc/occ/occ_edge.cpp b/libsrc/occ/occ_edge.cpp index 901417bc..7f2c3d66 100644 --- a/libsrc/occ/occ_edge.cpp +++ b/libsrc/occ/occ_edge.cpp @@ -7,8 +7,9 @@ namespace netgen { - OCCEdge::OCCEdge(TopoDS_Shape edge_) - : tedge(edge_.TShape()), + OCCEdge::OCCEdge(TopoDS_Shape edge_, GeometryVertex & start_, GeometryVertex & end_) + : GeometryEdge(start_, end_), + tedge(edge_.TShape()), edge(TopoDS::Edge(edge_)) { curve = BRep_Tool::Curve(edge, s0, s1); @@ -18,26 +19,13 @@ namespace netgen if(verts.size() != 2) throw Exception("OCC edge does not have 2 vertices"); - start = OCCVertex(verts[0]); - end = OCCVertex(verts[1]); - // swap start/end if necessary - double d00 = Dist(GetPoint(0), start.GetPoint()); - double d01 = Dist(GetPoint(0), end.GetPoint()); + double d00 = Dist(GetPoint(0), start->GetPoint()); + double d01 = Dist(GetPoint(0), end->GetPoint()); if(d01 < d00) swap(start, end); } - const GeometryVertex& OCCEdge::GetStartVertex() const - { - return start; - } - - const GeometryVertex& OCCEdge::GetEndVertex() const - { - return end; - } - double OCCEdge::GetLength() const { return props.Mass(); diff --git a/libsrc/occ/occ_edge.hpp b/libsrc/occ/occ_edge.hpp index c1d95428..f34cee56 100644 --- a/libsrc/occ/occ_edge.hpp +++ b/libsrc/occ/occ_edge.hpp @@ -14,6 +14,7 @@ namespace netgen { class OCCEdge : public GeometryEdge { + public: T_Shape tedge; TopoDS_Edge edge; Handle(Geom_Curve) curve; @@ -21,16 +22,11 @@ namespace netgen GProp_GProps props; public: - OCCVertex start; - OCCVertex end; - - OCCEdge(TopoDS_Shape edge_); + OCCEdge(TopoDS_Shape edge_, GeometryVertex & start_, GeometryVertex & end_); auto Shape() const { return edge; } T_Shape TShape() const { return tedge; } - const GeometryVertex& GetStartVertex() const override; - const GeometryVertex& GetEndVertex() const override; double GetLength() const override; Point<3> GetCenter() const override; Point<3> GetPoint(double t) const override; diff --git a/libsrc/occ/occ_face.cpp b/libsrc/occ/occ_face.cpp index d4d1a7c4..7b985f06 100644 --- a/libsrc/occ/occ_face.cpp +++ b/libsrc/occ/occ_face.cpp @@ -98,7 +98,6 @@ namespace netgen // auto cof = curve_on_face[ORIENTATION][edgenr]; auto edge = edge_on_face[ORIENTATION][edgenr]; - OCCEdge gedge(edge); double s0, s1; auto cof = BRep_Tool::CurveOnSurface (edge, face, s0, s1); diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index b013b073..278b0aab 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -1113,8 +1113,10 @@ namespace netgen auto tshape = v.TShape(); if(vertex_map.count(tshape)!=0) continue; - vertex_map[tshape] = vertices.Size(); auto occ_vertex = make_unique(TopoDS::Vertex(v)); + occ_vertex->nr = vertices.Size(); + vertex_map[tshape] = occ_vertex->nr; + if(global_shape_properties.count(tshape)>0) occ_vertex->properties = global_shape_properties[tshape]; vertices.Append(std::move(occ_vertex)); @@ -1129,10 +1131,9 @@ namespace netgen if(BRep_Tool::Degenerated(edge)) continue; edge_map[tshape] = edges.Size(); - auto occ_edge = make_unique(edge); + auto verts = GetVertices(e); + auto occ_edge = make_unique(edge, *vertices[vertex_map[verts[0].TShape()]], *vertices[vertex_map[verts[1].TShape()]] ); occ_edge->properties = global_shape_properties[tshape]; - occ_edge->start.nr = vertex_map[occ_edge->start.TShape()]; - occ_edge->end.nr = vertex_map[occ_edge->end.TShape()]; edges.Append(std::move(occ_edge)); } @@ -1145,6 +1146,10 @@ namespace netgen auto k = faces.Size(); face_map[tshape] = k; auto occ_face = make_unique(f); + + for(auto e : GetEdges(f)) + occ_face->edges.Append( edges[edge_map[e.TShape()]].get() ); + if(global_shape_properties.count(tshape)>0) occ_face->properties = global_shape_properties[tshape]; faces.Append(std::move(occ_face)); @@ -1194,6 +1199,9 @@ namespace netgen if(identifications.count(tshape)) for(auto & ident : identifications[tshape]) { + if(shape_map.count(ident.from)==0 || shape_map.count(ident.to)==0) + continue; + ShapeIdentification si{ shapes[shape_map[ident.from]].get(), shapes[shape_map[ident.to]].get(), @@ -1208,9 +1216,8 @@ namespace netgen add_identifications( edges, edge_map ); add_identifications( faces, face_map ); - ProcessIdentifications(); - bounding_box = ::netgen::GetBoundingBox( shape ); + ProcessIdentifications(); } @@ -1940,32 +1947,6 @@ namespace netgen return occ2ng( props.CentreOfMass() ); } - void OCCGeometry :: IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo) - { - Transformation<3> trafo; - if(opt_trafo) - trafo = occ2ng(*opt_trafo); - else - { - auto cme = GetCenter(me); - auto cyou = GetCenter(you); - trafo = Transformation<3>{cyou-cme}; - } - - identifications[me.TShape()].push_back( {me.TShape(), you.TShape(), trafo, name, type} ); - - auto vme = GetVertices(me); - auto vyou = GetVertices(you); - Point<3> pme0 = trafo(occ2ng(vme[0])); - Point<3> pme1 = trafo(occ2ng(vme[1])); - Point<3> pyou = occ2ng(vyou[0]); - - bool do_swap = Dist(pme1, pyou) < Dist(pme0, pyou); - - for(auto i : Range(2)) - identifications[vme[i].TShape()].push_back( {vme[i].TShape(), vyou[do_swap ? 1-i : i].TShape(), trafo, name, type} ); - } - bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you) { if(me.ShapeType() != you.ShapeType()) return false; @@ -1998,7 +1979,6 @@ namespace netgen vmap[s] = nullptr; } - bool all_verts_mapped = true; for (auto vert : verts_you) { auto s = vert.TShape(); @@ -2010,16 +1990,18 @@ namespace netgen return true; }); if(!vert_mapped) - { - all_verts_mapped = false; - break; - } + return false; } - return all_verts_mapped; + return true; } - void OCCGeometry :: IdentifyFaces(const TopoDS_Shape & solid, const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo) + void OCCGeometry :: Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo) { + auto type_me = me.ShapeType(); + auto type_you = you.ShapeType(); + if(type_me != type_you) + throw NgException ("Identify: cannot identify different shape types"); + Transformation<3> trafo; if(opt_trafo) { @@ -2032,25 +2014,33 @@ namespace netgen trafo = Transformation<3>(cyou-cme); } - auto faces_me = GetFaces(me); - auto faces_you = GetFaces(you); + ListOfShapes id_me; + ListOfShapes id_you; - for(auto face_me : faces_me) - for(auto face_you : faces_you) + if(auto faces_me = GetFaces(me); faces_me.size()>0) + { + id_me = faces_me; + id_you = GetFaces(you); + } + else if(auto edges_me = GetEdges(me); edges_me.size()>0) + { + id_me = edges_me; + id_you = GetEdges(you); + } + else + { + id_me = GetVertices(me); + id_you = GetVertices(you); + } + + for(auto shape_me : id_me) + for(auto shape_you : id_you) { - if(!IsMappedShape(trafo, face_me, face_you)) + if(!IsMappedShape(trafo, shape_me, shape_you)) continue; - identifications[face_me.TShape()].push_back - (OCCIdentification { face_me.TShape(), face_you.TShape(), trafo, name, type }); - - auto edges_me = GetEdges(me); - auto edges_you = GetEdges(you); - - for (auto e_me : edges_me) - for (auto e_you : edges_you) - if(IsMappedShape(trafo, e_me, e_you)) - IdentifyEdges(e_me, e_you, name, type, opt_trafo); + identifications[shape_me.TShape()].push_back + (OCCIdentification { shape_me.TShape(), shape_you.TShape(), trafo, name, type }); } } diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index 790d9727..e8d6a7b7 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -44,6 +44,7 @@ namespace netgen #define OCCGEOMETRYVISUALIZATIONFULLCHANGE 1 // Compute transformation matrices and redraw #define OCCGEOMETRYVISUALIZATIONHALFCHANGE 2 // Redraw + bool IsMappedShape(const Transformation<3> & trafo, const TopoDS_Shape & me, const TopoDS_Shape & you); class EntityVisualizationCode { @@ -341,8 +342,7 @@ namespace netgen bool ErrorInSurfaceMeshing (); // void WriteOCC_STL(char * filename); - static void IdentifyEdges(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo = nullopt); - static void IdentifyFaces(const TopoDS_Shape & solid,const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo = nullopt); + static void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo = nullopt); private: //bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const; diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index bf4d90ab..46fc1b4f 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -199,14 +199,14 @@ py::object CastShape(const TopoDS_Shape & s) template void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape) { - std::map mod_map; + std::map> mod_map; std::map tshape_handled; for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }) for (TopExp_Explorer e(shape, typ); e.More(); e.Next()) { auto tshape = e.Current().TShape(); - mod_map[tshape] = tshape; + mod_map[tshape].insert(tshape); tshape_handled[tshape] = false; } @@ -215,11 +215,8 @@ void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape) { auto tshape = e.Current().TShape(); - auto & modified = builder.Modified(e.Current()); - if(modified.Size()!=1) - continue; - - mod_map[tshape] = modified.First().TShape(); + for (auto mods : builder.Modified(e.Current())) + mod_map[tshape].insert(mods.TShape()); } for (auto typ : { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX }) @@ -236,21 +233,33 @@ void PropagateIdentifications (TBuilder & builder, TopoDS_Shape shape) auto tshape_mapped = mod_map[tshape]; - for(auto & ident : OCCGeometry::identifications[tshape]) + for(auto ident : OCCGeometry::identifications[tshape]) { - // update existing identification - if(tshape == tshape_mapped) - { - ident.to = mod_map[ident.to]; - ident.from = mod_map[ident.from]; - } - else - { - OCCIdentification id_new = ident; - id_new.to = mod_map[id_new.to]; - id_new.from = mod_map[id_new.from]; - OCCGeometry::identifications[mod_map[tshape]].push_back(id_new); - } + // nothing happened + if(mod_map[ident.to].size()==1 && mod_map[ident.from].size() ==1) + continue; + + auto from = ident.from; + auto to = ident.to; + + for(auto from_mapped : mod_map[from]) + for(auto to_mapped : mod_map[to]) + { + if(from==from_mapped && to==to_mapped) + continue; + + TopoDS_Shape s_from; s_from.TShape(from_mapped); + TopoDS_Shape s_to; s_to.TShape(to_mapped); + + if(!IsMappedShape(ident.trafo, s_from, s_to)) + continue; + + OCCIdentification id_new = ident; + id_new.to = to_mapped; + id_new.from = from_mapped; + auto id_owner = from == tshape ? from_mapped : to_mapped; + OCCGeometry::identifications[id_owner].push_back(id_new); + } } } } @@ -1071,27 +1080,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) BRepMesh_IncrementalMesh (shape, deflection, true); }) - .def("Identify", [](const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE idtype) { - // only edges supported, by now - auto type = me.ShapeType(); - auto tyou = you.ShapeType(); - if(type != tyou) - throw NgException ("Identify: cannot identify different shape types"); - - switch(type) - { - case TopAbs_VERTEX: - case TopAbs_EDGE: - OCCGeometry::IdentifyEdges(me, you, name, idtype); - break; - default: - throw NgException ("Identify: unsupported shape type"); - break; - } - }, py::arg("other"), py::arg("name"), py::arg("type")=Identifications::PERIODIC, "Identify shapes for periodic meshing") - - .def("Identify", OCCGeometry::IdentifyFaces, "Identify faces", - py::arg("from"), py::arg("to"), py::arg("name"), py::arg("type")=Identifications::PERIODIC, py::arg("trafo")=nullopt) + .def("Identify", OCCGeometry::Identify, py::arg("other"), py::arg("name"), + py::arg("type")=Identifications::PERIODIC, py::arg("trafo")=nullopt, + "Identify shapes for periodic meshing") .def("Distance", [](const TopoDS_Shape& self, const TopoDS_Shape& other) From d467621edd306caa5f546a62d3b81690abda842b Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Thu, 16 Dec 2021 19:00:10 +0100 Subject: [PATCH 5/7] change interface for identifications --- libsrc/meshing/meshfunc.cpp | 5 ++-- libsrc/occ/occgeom.cpp | 40 ++++++++++++++++++-------------- libsrc/occ/occgeom.hpp | 4 +++- libsrc/occ/python_occ_shapes.cpp | 9 ++++++- 4 files changed, 36 insertions(+), 22 deletions(-) diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 6714ccf8..68dcd347 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -539,8 +539,9 @@ namespace netgen if (md[i].mesh->CheckOverlappingBoundary()) throw NgException ("Stop meshing since boundary mesh is overlapping"); - // TODO: FillCloseSurface is still not working with CSG closesurfaces - // FillCloseSurface( md[i] ); + // TODO: FillCloseSurface is not working with CSG closesurfaces + if(md[i].mesh->GetGeometry()->GetGeomType() == Mesh::GEOM_OCC) + FillCloseSurface( md[i] ); CloseOpenQuads( md[i] ); MeshDomain(md[i]); }); diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 278b0aab..5a81f981 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -1995,42 +1995,46 @@ namespace netgen return true; } - void OCCGeometry :: Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo) + void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo) { - auto type_me = me.ShapeType(); - auto type_you = you.ShapeType(); - if(type_me != type_you) - throw NgException ("Identify: cannot identify different shape types"); - - Transformation<3> trafo; + gp_Trsf trafo; if(opt_trafo) { - trafo = occ2ng(*opt_trafo); + trafo = *opt_trafo; } else { - auto cme = GetCenter(me); - auto cyou = GetCenter(you); - trafo = Transformation<3>(cyou-cme); + auto v = GetCenter(you) - GetCenter(me); + trafo.SetTranslation(gp_Vec(v[0], v[1], v[2])); } + ListOfShapes list_me, list_you; + list_me.push_back(me); + list_you.push_back(you); + Identify(list_me, list_you, name, type, trafo); + } + + void Identify(const ListOfShapes & me, const ListOfShapes & you, string name, Identifications::ID_TYPE type, gp_Trsf occ_trafo) + { + Transformation<3> trafo = occ2ng(occ_trafo); + ListOfShapes id_me; ListOfShapes id_you; - if(auto faces_me = GetFaces(me); faces_me.size()>0) + if(auto faces_me = me.Faces(); faces_me.size()>0) { id_me = faces_me; - id_you = GetFaces(you); + id_you = you.Faces(); } - else if(auto edges_me = GetEdges(me); edges_me.size()>0) + else if(auto edges_me = me.Edges(); edges_me.size()>0) { id_me = edges_me; - id_you = GetEdges(you); + id_you = you.Edges(); } else { - id_me = GetVertices(me); - id_you = GetVertices(you); + id_me = me.Vertices(); + id_you = you.Vertices(); } for(auto shape_me : id_me) @@ -2039,7 +2043,7 @@ namespace netgen if(!IsMappedShape(trafo, shape_me, shape_you)) continue; - identifications[shape_me.TShape()].push_back + OCCGeometry::identifications[shape_me.TShape()].push_back (OCCIdentification { shape_me.TShape(), shape_you.TShape(), trafo, name, type }); } } diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index e8d6a7b7..0619cb97 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -342,11 +342,13 @@ namespace netgen bool ErrorInSurfaceMeshing (); // void WriteOCC_STL(char * filename); - static void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo = nullopt); private: //bool FastProject (int surfi, Point<3> & ap, double& u, double& v) const; }; + + void Identify(const ListOfShapes & me, const ListOfShapes & you, string name, Identifications::ID_TYPE type, gp_Trsf occ_trafo); + void Identify(const TopoDS_Shape & me, const TopoDS_Shape & you, string name, Identifications::ID_TYPE type, std::optional opt_trafo); void PrintContents (OCCGeometry * geom); diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index 46fc1b4f..178dfbc5 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -1080,7 +1080,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) BRepMesh_IncrementalMesh (shape, deflection, true); }) - .def("Identify", OCCGeometry::Identify, py::arg("other"), py::arg("name"), + + .def("Identify", py::overload_cast>(&Identify), + py::arg("other"), py::arg("name"), py::arg("type")=Identifications::PERIODIC, py::arg("trafo")=nullopt, "Identify shapes for periodic meshing") @@ -1625,6 +1627,11 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) } }, "set hpref for all elements of list") + .def("Identify", py::overload_cast(&Identify), + py::arg("other"), py::arg("name"), + py::arg("type")=Identifications::PERIODIC, py::arg("trafo"), + "Identify shapes for periodic meshing") + ; From 2d3c0e7186c3462947708f8c1e3a4174170ad375 Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Thu, 16 Dec 2021 20:54:56 +0100 Subject: [PATCH 6/7] add tests for occ identifications --- tests/pytest/test_occ_identifications.py | 76 ++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 tests/pytest/test_occ_identifications.py diff --git a/tests/pytest/test_occ_identifications.py b/tests/pytest/test_occ_identifications.py new file mode 100644 index 00000000..d15abc5e --- /dev/null +++ b/tests/pytest/test_occ_identifications.py @@ -0,0 +1,76 @@ +import pytest + +from netgen.meshing import IdentificationType +idtype = IdentificationType.CLOSESURFACES + +def test_two_boxes(): + occ = pytest.importorskip("netgen.occ") + inner = occ.Box((0,0,0), (1,1,1)) + trafo = occ.gp_Trsf().Scale(inner.center, 1.1) + outer = trafo(inner) + + inner.Identify(outer, "", idtype, trafo) + shape = occ.Glue([outer-inner, inner]) + + geo = occ.OCCGeometry(shape) + mesh = geo.GenerateMesh(maxh=0.3) + have_prisms = False + + for el in mesh.Elements3D(): + if len(el.vertices)==6: + have_prisms = True + break + + assert have_prisms + +def test_two_circles(): + occ = pytest.importorskip("netgen.occ") + circ1 = occ.WorkPlane().Circle(1).Face() + trafo = occ.gp_Trsf().Scale(circ1.center, 1.1) + + circ2 = trafo(circ1) + circ1.edges[0].Identify(circ2.edges[0], "", idtype, trafo) + circ2 -= circ1 + shape = occ.Glue([circ1, circ2]) + + geo = occ.OCCGeometry(shape, 2) + mesh = geo.GenerateMesh(maxh=0.2) + have_quads = False + + for el in mesh.Elements2D(): + if len(el.vertices)==4: + have_quads = True + break + + assert have_quads + +def test_cut_identified_face(): + occ = pytest.importorskip("netgen.occ") + from netgen.occ import Z, Box, Cylinder, Glue, OCCGeometry + box = Box((-1,-1,0), (1,1,1)) + cyl = Cylinder( (0,0,0), Z, 0.5, 1 ) + + box.faces.Min(Z).Identify(box.faces.Max(Z), "", idtype) + shape = Glue([cyl, box]) + geo = OCCGeometry(shape) + mesh = geo.GenerateMesh(maxh=0.5) + + for el in mesh.Elements3D(): + assert len(el.vertices)==6 + +def test_identify_multiple_faces(): + occ = pytest.importorskip("netgen.occ") + from netgen.occ import Z, Box, Cylinder, Glue, OCCGeometry, gp_Trsf + box = Box((-1,-1,0), (1,1,1)) + cyl = Cylinder( (0,0,0), Z, 0.5, 1 ) + + shape = Glue([box, cyl]) + bot_faces = shape.faces[Z < 0.1] + top_faces = shape.faces[Z > 0.1] + bot_faces.Identify(top_faces, "", idtype, gp_Trsf.Translation((0,0,1))) + + geo = OCCGeometry(shape) + mesh = geo.GenerateMesh(maxh=0.3) + + for el in mesh.Elements3D(): + assert len(el.vertices)==6 From 1e86bc2c59832f3bfcf9381b618f8e4cf4f458e2 Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Fri, 17 Dec 2021 10:38:15 +0100 Subject: [PATCH 7/7] occ - consistent ordering of shapes --- libsrc/occ/occgeom.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 5a81f981..1b84b5ac 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -1108,8 +1108,9 @@ namespace netgen fsingular = esingular = vsingular = false; // Add shapes - for(auto v : GetVertices(shape)) + for(auto i1 : Range(1, vmap.Extent()+1)) { + auto v = vmap(i1); auto tshape = v.TShape(); if(vertex_map.count(tshape)!=0) continue; @@ -1122,8 +1123,9 @@ namespace netgen vertices.Append(std::move(occ_vertex)); } - for(auto e : GetEdges(shape)) + for(auto i1 : Range(1, emap.Extent()+1)) { + auto e = emap(i1); auto tshape = e.TShape(); auto edge = TopoDS::Edge(e); if(edge_map.count(tshape)!=0) @@ -1137,8 +1139,9 @@ namespace netgen edges.Append(std::move(occ_edge)); } - for(auto f : GetFaces(shape)) + for(auto i1 : Range(1, fmap.Extent()+1)) { + auto f = fmap(i1); auto tshape = f.TShape(); if(face_map.count(tshape)==0) { @@ -1167,8 +1170,9 @@ namespace netgen } - for(auto s : GetSolids(shape)) + for(auto i1 : Range(1, somap.Extent()+1)) { + auto s = somap(i1); auto tshape = s.TShape(); int k; if(solid_map.count(tshape)==0)