Merge branch 'occ_fix_maxh' into 'master'

Propagate OCC maxh settings correctly

See merge request ngsolve/netgen!673

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<const GeometryVertex*> GetFaceVertices(const GeometryFace& face) const { return Array<const GeometryVertex*>{}; }
 
     void Clear();

libsrc/occ/occ_utils.cpp

@@ -5,6 +5,7 @@
 #include <BRepTools.hxx>
 
 #include "occ_utils.hpp"
+#include "occgeom.hpp"
 
 namespace netgen
 {

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:

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))
+      auto & props = geom.GetSolid(dom).properties;
-        continue;
-      auto& props = OCCGeometry::GetProperties(s);
       maxhdom[dom] = min2(maxhdom[dom], props.maxh);
       maxlayer = max2(maxlayer, props.layer);
-      dom++;
     }
-    for(const auto& f : GetFaces(geom.GetShape()))
+
-      if(OCCGeometry::HaveProperties(f))
+    for(auto & f : geom.Faces())
-        maxlayer = max2(maxlayer, OCCGeometry::GetProperties(f).layer);
+      maxlayer = max2(maxlayer, f->properties.layer);
-    for(const auto& e : GetEdges(geom.GetShape()))
+
-      if(OCCGeometry::HaveProperties(e))
+    for(auto & e : geom.Edges())
-        maxlayer = max2(maxlayer, OCCGeometry::GetProperties(e).layer);
+      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 = v->properties;
-              {
+            if(props.maxh < 1e99)
-                auto& props = OCCGeometry::GetProperties(v);
+                mesh.GetLocalH(props.layer)->SetH(v->GetPoint(), props.maxh);
-                if(props.maxh < 1e99)
-                  mesh.GetLocalH(props.layer)->SetH(occ2ng(BRep_Tool::Pnt(TopoDS::Vertex(v))), 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;

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<OCCSolid>(s);
           if(HaveProperties(s))
             occ_solid->properties = GetProperties(s);
-          solids.Append(std::move(occ_solid));
-
           for(auto f : GetFaces(s))
           {
               auto & face = static_cast<OCCFace&>(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<OCCSolid&>(*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,

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<GeometryFace&>(as_const(*this).GetFace(shape));
     }
 
+    GeometrySolid & GetSolid(const TopoDS_Shape & shape)
+    {
+        return const_cast<GeometrySolid&>(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;

libsrc/occ/python_occ_shapes.cpp

@@ -827,6 +827,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
           return nullopt;
       }, [](const TopoDS_Shape & self, optional<string> 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;
+                    OCCGeometry::GetProperties(self).hpref = val;
-                    hpref = max2(val, hpref);
+                    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 : shapes)
-          {
+          OCCGeometry::GetProperties(s).maxh = maxh;
-            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;
-          }
       }, "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)
-          {
+          OCCGeometry::GetProperties(shape).hpref = hpref;
-            auto& val = OCCGeometry::GetProperties(shape).hpref;
-            val = max2(hpref, val);
-          }
       }, "set hpref for all elements of list")
     .def_property("quad_dominated", [](ListOfShapes& shapes)
                   {