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)