Merge branch 'occ_better_divide_edge' into 'master'

OCC better divide edge algorithm

See merge request ngsolve/netgen!551

libsrc/meshing/basegeom.cpp

@@ -479,61 +479,82 @@ namespace netgen
     mesh.LoadLocalMeshSize(mparam.meshsizefilename);
   }
 
-  void DivideEdge(GeometryEdge * edge, const MeshingParameters & mparam, const Mesh & mesh, Array<Point<3>> & points, Array<double> & params, int layer)
+  void GeometryEdge :: Divide(const MeshingParameters & mparam, const Mesh & mesh, Array<Point<3>> & points, Array<double> & params)
   {
-      static Timer tdivedgesections("Divide edge sections");
+    static Timer tdivedgesections("Divide edge sections");
-      static Timer tdivide("Divide Edges");
+    static Timer tdivide("Divide Edges");
-      RegionTimer rt(tdivide);
+    RegionTimer rt(tdivide);
-      // -------------------- DivideEdge -----------------
+    // -------------------- DivideEdge -----------------
-      static constexpr size_t divide_edge_sections = 10000;
+    tdivedgesections.Start();
-      double hvalue[divide_edge_sections+1];
+    auto layer = properties.layer;
-      hvalue[0] = 0;
+    double safety = 0.5*(1.-mparam.grading);
 
-      Point<3> old_pt = edge->GetPoint(0.);
+    double lam = 0.0;
-      // calc local h for edge
+    Point<3> p = GetPoint(0.0);
-      tdivedgesections.Start();
+    auto old_p = p;
-      for(auto i : Range(divide_edge_sections))
+    Array<double> hvalue, fine_params;
+    hvalue.Append(.0);
+
+    while (lam<1. && hvalue.Size() < 20000) {
+      fine_params.Append(lam);
+      auto h = mesh.GetH(old_p, layer);
+      auto step = safety * h/GetTangent(lam).Length();
+      lam += step;
+      lam = min2(lam, 1.0);
+      p = GetPoint(lam);
+      hvalue.Append((hvalue.Size()==0 ? 0.0 : hvalue.Last()) + 1./h * (p-old_p).Length());
+      old_p = p;
+    }
+
+    fine_params.Append(1.0);
+
+    if(hvalue.Size()==20000 && lam<1.0)
+      cout << "Warning: Could not divide Edge" << endl;
+
+    tdivedgesections.Stop();
+
+    auto n = hvalue.Size()-1;
+    int nsubedges = max2(1, int(floor(hvalue.Last()+0.5)));
+    points.SetSize(nsubedges-1);
+    params.SetSize(nsubedges+1);
+
+    int i1 = 0;
+    for(auto i : Range(1,nsubedges))
+    {
+      auto h_target = i*hvalue.Last()/nsubedges;
+      while(hvalue[i1]<h_target && i1<hvalue.Size())
+        i1++;
+
+      if(i1==hvalue.Size())
       {
-          auto pt = edge->GetPoint(double(i+1)/divide_edge_sections);
+        points.SetSize(i-1);
-          hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pt, layer) * (pt-old_pt).Length();
+        params.SetSize(i+1);
-          old_pt = pt;
+        cout << "divide edge: local h too small" << endl;
+        break;
       }
-      int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5)));
-      tdivedgesections.Stop();
-      points.SetSize(nsubedges-1);
-      params.SetSize(nsubedges+1);
 
-      int i = 1;
+      // interpolate lam between points
-      int i1 = 0;
+      auto lam0 = fine_params[i1-1];
-      do
+      auto lam1 = fine_params[i1];
-      {
+      auto h0 = hvalue[i1-1];
-          if (hvalue[i1]/hvalue[divide_edge_sections]*nsubedges >= i)
+      auto h1 = hvalue[i1];
-          {
-              params[i] = (double(i1)/divide_edge_sections);
-              points[i-1] = MeshPoint(edge->GetPoint(params[i]));
-              i++;
-          }
-          i1++;
-          if (i1 > divide_edge_sections)
-          {
-              nsubedges = i;
-              points.SetSize(nsubedges-1);
-              params.SetSize(nsubedges+1);
-              cout << "divide edge: local h too small" << endl;
-          }
 
-      } while(i < nsubedges);
+      auto fac = (h_target-h0)/(h1-h0);
+      auto lam = lam0 + fac*(lam1-lam0);
+      params[i] = lam;
+      points[i-1] = MeshPoint(GetPoint(params[i]));
+    }
 
-      params[0] = 0.;
+    params[0] = 0.;
-      params[nsubedges] = 1.;
+    params[nsubedges] = 1.;
 
-      if(params[nsubedges] <= params[nsubedges-1])
+    if(params[nsubedges] <= params[nsubedges-1])
-      {
+    {
-          cout << "CORRECTED" << endl;
+        cout << "CORRECTED" << endl;
-          points.SetSize (nsubedges-2);
+        points.SetSize (nsubedges-2);
-          params.SetSize (nsubedges);
+        params.SetSize (nsubedges);
-          params[nsubedges-1] = 1.;
+        params[nsubedges-1] = 1.;
-      }
+    }
   }
 
   void NetgenGeometry :: FindEdges(Mesh& mesh,
@@ -616,7 +637,7 @@ namespace netgen
             }
             else
             {
-                DivideEdge(edge, mparam, mesh, edge_points, params, edge->properties.layer);
+                edge->Divide(mparam, mesh, edge_points, params);
             }
         }
         else

libsrc/meshing/basegeom.hpp

@@ -110,6 +110,7 @@ namespace netgen
     }
     virtual Vec<3> GetTangent(double t) const = 0;
     virtual bool IsMappedShape( const GeometryShape & other, const Transformation<3> & trafo, double tolerance ) const override;
+    virtual void Divide(const MeshingParameters & mparam, const Mesh & mesh, Array<Point<3>> & points, Array<double> & params);
   };
 
   class DLL_HEADER GeometryFace : public GeometryShape

libsrc/occ/occ_edge.cpp

@@ -69,6 +69,7 @@ namespace netgen
         gp_Pnt p;
         gp_Vec v;
         curve->D1(t, p, v);
-        return occ2ng(v);
+        return occ2ng(v) * (s1-s0);
     }
+
 }

libsrc/occ/occ_face.cpp

@@ -1,6 +1,7 @@
 #include <BRepGProp.hxx>
 #include <BRep_Tool.hxx>
 #include <GeomAPI_ProjectPointOnCurve.hxx>
+#include <BRepLProp_SLProps.hxx>
 
 #include "occ_edge.hpp"
 #include "occ_face.hpp"
@@ -111,9 +112,16 @@ namespace netgen
 
                 for(auto i : Range(2))
                 {
+                    // take uv from CurveOnSurface as start value but project again for better accuracy
+                    // (cof->Value yields wrong values (outside of surface) for complicated faces
                     auto uv = cof->Value(s[i]);
-                    seg.epgeominfo[i].u = uv.X();
+                    PointGeomInfo gi;
-                    seg.epgeominfo[i].v = uv.Y();
+                    gi.u = uv.X();
+                    gi.v = uv.Y();
+                    Point<3> pproject = mesh[seg[i]];
+                    ProjectPointGI(pproject, gi);
+                    seg.epgeominfo[i].u = gi.u;
+                    seg.epgeominfo[i].v = gi.v;
                 }
 
                 bool do_swap = ORIENTATION == REVERSED;
@@ -234,7 +242,11 @@ namespace netgen
 
     double OCCFace::GetCurvature(const PointGeomInfo& gi) const
     {
-        throw Exception(ToString("not implemented") + __FILE__ + ":" + ToString(__LINE__));
+        BRepAdaptor_Surface sf(face, Standard_True);
+        BRepLProp_SLProps prop2(sf, 2, 1e-5);
+        prop2.SetParameters (gi.u, gi.v);
+        return max(fabs(prop2.MinCurvature()),
+                   fabs(prop2.MaxCurvature()));
     }
 
     void OCCFace::RestrictH(Mesh& mesh, const MeshingParameters& mparam) const

libsrc/occ/occ_utils.hpp

@@ -276,13 +276,12 @@ namespace netgen
     }
   };
 
-
+  inline auto Properties (TopoDS_Shape shape)
-  inline gp_Pnt Center (TopoDS_Shape shape)
   {
     GProp_GProps props;
     double tol;
     switch (shape.ShapeType())
-      {
+    {
       case TopAbs_SOLID:
       case TopAbs_COMPOUND:
       case TopAbs_COMPSOLID:
@@ -298,8 +297,18 @@ namespace netgen
         BRepGProp::LinearProperties(shape, props, tol); break;
       default:
         BRepGProp::LinearProperties(shape, props);
-      }
+    }
-    return props.CentreOfMass();
+    return props;
+  }
+
+  inline gp_Pnt Center (TopoDS_Shape shape)
+  {
+    return Properties(shape).CentreOfMass();
+  }
+
+  inline double Mass (TopoDS_Shape shape)
+  {
+    return Properties(shape).Mass();
   }
 
 }

libsrc/occ/occgenmesh.cpp

@@ -537,9 +537,7 @@ namespace netgen
             multithread.percent = 100 * (i-1)/double(nedges);
             if (BRep_Tool::Degenerated(e)) continue;
 
-            GProp_GProps system;
+            double len = Mass(e);
-            BRepGProp::LinearProperties(e, system);
-            double len = system.Mass();
 
             if (len < mincurvelength)
               {
@@ -587,7 +585,7 @@ namespace netgen
             localh = min2(localh, OCCGeometry::GetProperties(e).maxh);
             maxedgelen = max (maxedgelen, len);
             minedgelen = min (minedgelen, len);
-            int maxj = max((int) ceil(len/localh), 2);
+            int maxj = max((int) ceil(len/localh)*2, 2);
 
             for (int j = 0; j <= maxj; j++)
               {

libsrc/occ/python_occ_shapes.cpp

@@ -715,23 +715,8 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
 
     .def("Properties", [] (const TopoDS_Shape & shape)
          {
-           GProp_GProps props;
+           auto props = Properties(shape);
-           switch (shape.ShapeType())
+           return tuple( py::cast(props.Mass()), py::cast(props.CentreOfMass()) );
-             {
-             case TopAbs_FACE:
-             case TopAbs_SHELL:
-               BRepGProp::SurfaceProperties (shape, props); break;
-	     case TopAbs_SOLID:
-	     case TopAbs_COMPOUND:
-	     case TopAbs_COMPSOLID:
-               BRepGProp::VolumeProperties (shape, props); break;
-             default:
-               BRepGProp::LinearProperties(shape, props);
-               // throw Exception("Properties implemented only for FACE");
-             }
-           double mass = props.Mass();
-           gp_Pnt center = props.CentreOfMass();
-           return tuple( py::cast(mass), py::cast(center) );
          }, "returns tuple of shape properties, currently ('mass', 'center'")
     
     .def_property_readonly("center", [](const TopoDS_Shape & shape) {
@@ -739,20 +724,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m)
       }, "returns center of gravity of shape")
     
     .def_property_readonly("mass", [](const TopoDS_Shape & shape) {
-           GProp_GProps props;
+           return Mass(shape);
-           switch (shape.ShapeType())
-             {
-             case TopAbs_FACE:
-             case TopAbs_SHELL:
-               BRepGProp::SurfaceProperties (shape, props); break;
-	     case TopAbs_SOLID:
-	     case TopAbs_COMPOUND:
-	     case TopAbs_COMPSOLID:
-               BRepGProp::VolumeProperties (shape, props); break;
-             default:
-               BRepGProp::LinearProperties(shape, props);
-             }
-           return props.Mass();
       }, "returns mass of shape, what is length, face, or volume")
 
     .def("Move", [](const TopoDS_Shape & shape, const gp_Vec v)