diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 8d5a0ee5..2eb7a203 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -212,11 +212,18 @@ namespace netgen size_t GetNVertices() const { return vertices.Size(); } size_t GetNEdges() const { return edges.Size(); } size_t GetNFaces() const { return faces.Size(); } + size_t GetNSolids() const { return solids.Size(); } + const GeometrySolid & GetSolid(int i) const { return *solids[i]; } const GeometryFace & GetFace(int i) const { return *faces[i]; } const GeometryEdge & GetEdge(int i) const { return *edges[i]; } const GeometryVertex & GetVertex(int i) const { return *vertices[i]; } + auto Solids() const { return FlatArray{solids}; } + auto Faces() const { return FlatArray{faces}; } + auto Edges() const { return FlatArray{edges}; } + auto Vertices() const { return FlatArray{vertices}; } + virtual Array GetFaceVertices(const GeometryFace& face) const { return Array{}; } void Clear(); diff --git a/libsrc/occ/occ_utils.cpp b/libsrc/occ/occ_utils.cpp index 243349fc..96a8a41e 100644 --- a/libsrc/occ/occ_utils.cpp +++ b/libsrc/occ/occ_utils.cpp @@ -5,6 +5,7 @@ #include #include "occ_utils.hpp" +#include "occgeom.hpp" namespace netgen { diff --git a/libsrc/occ/occ_utils.hpp b/libsrc/occ/occ_utils.hpp index 8696cb55..04973789 100644 --- a/libsrc/occ/occ_utils.hpp +++ b/libsrc/occ/occ_utils.hpp @@ -170,6 +170,17 @@ namespace netgen common.push_back(shape); return common; } + + ListOfShapes GetHighestDimShapes() const + { + for (auto type : {TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX}) + { + auto ret = SubShapes(type); + if (ret.size() > 0) + return ret; + } + return ListOfShapes(); + } }; inline ListOfShapes GetSolids(const TopoDS_Shape & shape) @@ -212,6 +223,16 @@ namespace netgen return sub; } + inline ListOfShapes GetHighestDimShapes(const TopoDS_Shape & shape) + { + auto ret = GetSolids(shape); if(ret.size() > 0) return ret; + ret = GetFaces(shape); if(ret.size() > 0) return ret; + ret = GetEdges(shape); if(ret.size() > 0) return ret; + ret = GetVertices(shape); if(ret.size() > 0) return ret; + return ListOfShapes(); + } + + class DirectionalInterval { public: diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp index 9a6c269e..1ce34c55 100644 --- a/libsrc/occ/occgenmesh.cpp +++ b/libsrc/occ/occgenmesh.cpp @@ -414,7 +414,7 @@ namespace netgen // Philippose - 15/01/2009 - auto& props = OCCGeometry::GetProperties(geom.fmap(k)); + auto& props = occface.properties; double maxh = min2(geom.face_maxh[k-1], props.maxh); //double maxh = mparam.maxh; // int noldpoints = mesh->GetNP(); @@ -485,21 +485,18 @@ namespace netgen int maxlayer = 1; int dom = 0; - for(const auto& s : GetSolids(geom.GetShape())) + for(auto dom : Range(geom.GetNSolids())) { - if(!OCCGeometry::HaveProperties(s)) - continue; - auto& props = OCCGeometry::GetProperties(s); + auto & props = geom.GetSolid(dom).properties; maxhdom[dom] = min2(maxhdom[dom], props.maxh); maxlayer = max2(maxlayer, props.layer); - dom++; } - for(const auto& f : GetFaces(geom.GetShape())) - if(OCCGeometry::HaveProperties(f)) - maxlayer = max2(maxlayer, OCCGeometry::GetProperties(f).layer); - for(const auto& e : GetEdges(geom.GetShape())) - if(OCCGeometry::HaveProperties(e)) - maxlayer = max2(maxlayer, OCCGeometry::GetProperties(e).layer); + + for(auto & f : geom.Faces()) + maxlayer = max2(maxlayer, f->properties.layer); + + for(auto & e : geom.Edges()) + maxlayer = max2(maxlayer, e->properties.layer); mesh.SetMaxHDomain (maxhdom); @@ -514,14 +511,11 @@ namespace netgen for(auto layer : Range(1, maxlayer+1)) mesh.SetLocalH (bb.PMin(), bb.PMax(), mparam.grading, layer); - for(const auto& v : GetVertices(geom.GetShape())) + for(auto& v : geom.Vertices()) { - if(OCCGeometry::HaveProperties(v)) - { - auto& props = OCCGeometry::GetProperties(v); - if(props.maxh < 1e99) - mesh.GetLocalH(props.layer)->SetH(occ2ng(BRep_Tool::Pnt(TopoDS::Vertex(v))), props.maxh); - } + auto& props = v->properties; + if(props.maxh < 1e99) + mesh.GetLocalH(props.layer)->SetH(v->GetPoint(), props.maxh); } int nedges = geom.emap.Extent(); @@ -538,19 +532,10 @@ namespace netgen multithread.task = "Setting local mesh size (elements per edge)"; - // Philippose - 23/01/2009 - // Find all the parent faces of a given edge - // and limit the mesh size of the edge based on the - // mesh size limit of the face - TopTools_IndexedDataMapOfShapeListOfShape edge_face_map; - edge_face_map.Clear(); - TopExp::MapShapesAndAncestors(geom.shape, TopAbs_EDGE, TopAbs_FACE, edge_face_map); - // setting elements per edge for (int i = 1; i <= nedges && !multithread.terminate; i++) { TopoDS_Edge e = TopoDS::Edge (geom.emap(i)); - int layer = OCCGeometry::GetProperties(e).layer; multithread.percent = 100 * (i-1)/double(nedges); if (BRep_Tool::Degenerated(e)) continue; @@ -583,23 +568,10 @@ namespace netgen double localh = len/mparam.segmentsperedge; double s0, s1; - const TopTools_ListOfShape& parent_faces = edge_face_map.FindFromKey(e); - - TopTools_ListIteratorOfListOfShape parent_face_list; - - for(parent_face_list.Initialize(parent_faces); parent_face_list.More(); parent_face_list.Next()) - { - TopoDS_Face parent_face = TopoDS::Face(parent_face_list.Value()); - - int face_index = geom.fmap.FindIndex(parent_face); - - if(face_index >= 1) localh = min(localh,geom.face_maxh[face_index - 1]); - localh = min2(localh, OCCGeometry::GetProperties(parent_face).maxh); - } - Handle(Geom_Curve) c = BRep_Tool::Curve(e, s0, s1); - localh = min2(localh, OCCGeometry::GetProperties(e).maxh); + const auto & props = gedge.properties; + localh = min2(localh, props.maxh); maxedgelen = max (maxedgelen, len); minedgelen = min (minedgelen, len); int maxj = max((int) ceil(len/localh)*2, 2); @@ -607,7 +579,7 @@ namespace netgen for (int j = 0; j <= maxj; j++) { gp_Pnt pnt = c->Value (s0+double(j)/maxj*(s1-s0)); - mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh, layer); + mesh.RestrictLocalH (Point3d(pnt.X(), pnt.Y(), pnt.Z()), localh, props.layer); } } @@ -622,12 +594,12 @@ namespace netgen double maxcur = 0; multithread.percent = 100 * (i-1)/double(nedges); TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); - int layer = OCCGeometry::GetProperties(edge).layer; if (BRep_Tool::Degenerated(edge)) continue; double s0, s1; Handle(Geom_Curve) c = BRep_Tool::Curve(edge, s0, s1); BRepAdaptor_Curve brepc(edge); BRepLProp_CLProps prop(brepc, 2, 1e-5); + auto layer = geom.GetEdge(edge).properties.layer; for (int j = 1; j <= nsections; j++) { @@ -658,7 +630,6 @@ namespace netgen { multithread.percent = 100 * (i-1)/double(nfaces); TopoDS_Face face = TopoDS::Face(geom.fmap(i)); - int layer = OCCGeometry::GetProperties(face).layer; TopLoc_Location loc; Handle(Geom_Surface) surf = BRep_Tool::Surface (face); Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); @@ -675,6 +646,7 @@ namespace netgen // one prop for evaluating and one for derivatives BRepLProp_SLProps prop(sf, 0, 1e-5); BRepLProp_SLProps prop2(sf, 2, 1e-5); + auto layer = geom.GetFace(face).properties.layer; int ntriangles = triangulation -> NbTriangles(); for (int j = 1; j <= ntriangles; j++) @@ -725,7 +697,7 @@ namespace netgen for (int i = 1; i <= nedges && !multithread.terminate; i++) { TopoDS_Edge edge = TopoDS::Edge (geom.emap(i)); - int layer = OCCGeometry::GetProperties(edge).layer; + int layer = geom.GetEdge(edge).properties.layer; if (BRep_Tool::Degenerated(edge)) continue; double s0, s1; diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index f1607bd7..13156783 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -200,6 +200,11 @@ namespace netgen return *faces[fmap.FindIndex(shape)-1]; } + const GeometrySolid & OCCGeometry :: GetSolid(const TopoDS_Shape & shape) const + { + return *solids[somap.FindIndex(shape)-1]; + } + string STEP_GetEntityName(const TopoDS_Shape & theShape, STEPCAFControl_Reader * aReader) { @@ -1239,11 +1244,11 @@ namespace netgen auto occ_solid = make_unique(s); if(HaveProperties(s)) occ_solid->properties = GetProperties(s); - solids.Append(std::move(occ_solid)); - for(auto f : GetFaces(s)) { auto & face = static_cast(GetFace(f)); + face.properties.maxh = min2(face.properties.maxh, occ_solid->properties.maxh); + if(face.domin==-1) face.domin = k; else @@ -1251,21 +1256,10 @@ namespace netgen if(face.Shape().Orientation() == TopAbs_INTERNAL) face.domout = k; } + solids.Append(std::move(occ_solid)); } // Propagate maxh to children - for(auto& solid : solids) - { - auto& shape = static_cast(*solid).GetShape(); - if(!OCCGeometry::HaveProperties(shape)) - continue; - for(auto& f : GetFaces(shape)) - { - auto& face = GetFace(f); - face.properties.maxh = min2(face.properties.maxh, - GetProperties(shape).maxh); - } - } for(auto& face : faces) for(auto& edge : face->edges) edge->properties.maxh = min2(edge->properties.maxh, diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index fa00163a..6d2abbc8 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -232,6 +232,7 @@ namespace netgen using NetgenGeometry::GetVertex; using NetgenGeometry::GetEdge; using NetgenGeometry::GetFace; + using NetgenGeometry::GetSolid; GeometryShape & GetShape(const TopoDS_Shape & shape) { @@ -252,10 +253,16 @@ namespace netgen return const_cast(as_const(*this).GetFace(shape)); } + GeometrySolid & GetSolid(const TopoDS_Shape & shape) + { + return const_cast(as_const(*this).GetSolid(shape)); + } + const GeometryShape & GetShape(const TopoDS_Shape & shape) const; const GeometryVertex & GetVertex(const TopoDS_Shape & shape) const; const GeometryEdge & GetEdge(const TopoDS_Shape & shape) const; const GeometryFace & GetFace(const TopoDS_Shape & shape) const; + const GeometrySolid & GetSolid(const TopoDS_Shape & shape) const; void Analyse(Mesh& mesh, const MeshingParameters& mparam) const override; diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index e2c819d1..62c8f387 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -827,6 +827,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) return nullopt; }, [](const TopoDS_Shape & self, optional name) { OCCGeometry::GetProperties(self).name = name; + for (auto & s : GetHighestDimShapes(self)) + OCCGeometry::GetProperties(s).name = name; }, "'name' of shape") .def_property("maxh", @@ -837,6 +839,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) [](TopoDS_Shape& self, double val) { OCCGeometry::GetProperties(self).maxh = val; + for(auto & s : GetHighestDimShapes(self)) + OCCGeometry::GetProperties(s).maxh = val; }, "maximal mesh-size for shape") .def_property("hpref", @@ -846,8 +850,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) }, [](TopoDS_Shape& self, double val) { - auto & hpref = OCCGeometry::GetProperties(self).hpref; - hpref = max2(val, hpref); + OCCGeometry::GetProperties(self).hpref = val; + for(auto & s : GetHighestDimShapes(self)) + OCCGeometry::GetProperties(s).hpref = val; }, "number of refinement levels for geometric refinement") .def_property("col", [](const TopoDS_Shape & self) -> py::object { @@ -860,6 +865,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) if(c.size() == 4) col[3] = c[3]; OCCGeometry::GetProperties(self).col = col; + for(auto & s : GetHighestDimShapes(self)) + OCCGeometry::GetProperties(s).col = col; }, "color of shape as RGB - tuple") .def_property("layer", [](const TopoDS_Shape& self) { if (!OCCGeometry::HaveProperties(self)) @@ -867,6 +874,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) return OCCGeometry::GetProperties(self).layer; }, [](const TopoDS_Shape& self, int layer) { OCCGeometry::GetProperties(self).layer = layer; + for(auto & s : GetHighestDimShapes(self)) + OCCGeometry::GetProperties(s).layer = layer; }, "layer of shape") .def("UnifySameDomain", [](const TopoDS_Shape& shape, bool edges, bool faces, @@ -1813,17 +1822,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) }, [](ListOfShapes& shapes, double maxh) { - for(auto& shape : shapes) - { - for(auto& s : GetSolids(shape)) - OCCGeometry::GetProperties(s).maxh = maxh; - for(auto& s : GetFaces(shape)) - OCCGeometry::GetProperties(s).maxh = maxh; - for(auto& s : GetEdges(shape)) - OCCGeometry::GetProperties(s).maxh = maxh; - for(auto& s : GetVertices(shape)) - OCCGeometry::GetProperties(s).maxh = maxh; - } + for(auto & s : shapes) + OCCGeometry::GetProperties(s).maxh = maxh; }, "set maxh for all elements of list") .def_property("hpref", [](ListOfShapes& shapes) { @@ -1832,10 +1832,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) [](ListOfShapes& shapes, double hpref) { for(auto& shape : shapes) - { - auto& val = OCCGeometry::GetProperties(shape).hpref; - val = max2(hpref, val); - } + OCCGeometry::GetProperties(shape).hpref = hpref; }, "set hpref for all elements of list") .def_property("quad_dominated", [](ListOfShapes& shapes) {