clean up geometry interface, fix if number of subdivided edges was not correct

libsrc/meshing/basegeom.cpp

@@ -10,14 +10,6 @@ namespace netgen
   GeometryRegister :: ~GeometryRegister()
   { ; }
 
-  Array<Point<3>> GeometryEdge :: GetEquidistantPointArray(size_t npoints) const
-  {
-    Array<Point<3>> pts(npoints);
-    for(auto i : Range(npoints))
-      pts[i] = GetPoint(double(i)/(npoints-1));
-    return pts;
-  }
-
   void GeometryFace :: RestrictHTrig(Mesh& mesh,
                                      const PointGeomInfo& gi0,
                                      const PointGeomInfo& gi1,
@@ -148,6 +140,8 @@ namespace netgen
     static Timer t1("MeshEdges"); RegionTimer regt(t1);
     static Timer tdivide("Divide Edges");
     static Timer tdivedgesections("Divide edge sections");
+    const char* savetask = multithread.task;
+    multithread.task = "Mesh Edges";
 
     // create face descriptors and set bc names
     mesh.SetNBCNames(faces.Size());
@@ -196,11 +190,15 @@ namespace netgen
                 double hvalue[divide_edge_sections+1];
                 hvalue[0] = 0;
 
+                Point<3> old_pt = edge->GetPoint(0.);
                 // calc local h for edge
                 tdivedgesections.Start();
-                auto edgepts = edge->GetEquidistantPointArray(divide_edge_sections+1);
+                for(auto i : Range(divide_edge_sections))
-                for(auto i : Range(edgepts.Size()-1))
+                  {
-                  hvalue[i+1] = hvalue[i] + 1./mesh.GetH(edgepts[i+1]) * (edgepts[i+1]-edgepts[i]).Length();
+                    auto pt = edge->GetPoint(double(i+1)/divide_edge_sections);
+                    hvalue[i+1] = hvalue[i] + 1./mesh.GetH(pt) * (pt-old_pt).Length();
+                    old_pt = pt;
+                  }
                 int nsubedges = max2(1, int(floor(hvalue[divide_edge_sections]+0.5)));
                 tdivedgesections.Stop();
                 mps.SetSize(nsubedges-1);
@@ -235,7 +233,7 @@ namespace netgen
                     cout << "CORRECTED" << endl;
                     mps.SetSize (nsubedges-2);
                     params.SetSize (nsubedges);
-                    params[nsubedges] = 1.;
+                    params[nsubedges-1] = 1.;
                   }
                 tdivide.Stop();
                 // ----------- Add Points to mesh and create segments -----
@@ -292,6 +290,8 @@ namespace netgen
               }
           }
       }
+    mesh.CalcSurfacesOfNode();
+    multithread.task = savetask;
   }
 
   void NetgenGeometry :: MeshSurface(Mesh& mesh,

libsrc/meshing/basegeom.hpp

@@ -32,7 +32,15 @@ namespace netgen
     virtual double CalcStep(double t, double sag) const = 0;
     virtual bool OrientedLikeGlobal() const = 0;
     virtual size_t GetHash() const = 0;
-    virtual Array<Point<3>> GetEquidistantPointArray(size_t npoints) const;
+    virtual void ProjectPoint(Point<3>& p, EdgePointGeomInfo& gi) const = 0;
+    virtual void PointBetween(const Point<3>& p1,
+                              const Point<3>& p2,
+                              double secpoint,
+                              const EdgePointGeomInfo& gi1,
+                              const EdgePointGeomInfo& gi2,
+                              Point<3>& newp,
+                              EdgePointGeomInfo& newgi) const = 0;
+    virtual Vec<3> GetTangent(const Point<3>& p) const = 0;
   };
 
   class GeometryFace
@@ -42,9 +50,11 @@ namespace netgen
     virtual size_t GetNBoundaries() const = 0;
     virtual Array<unique_ptr<GeometryEdge>> GetBoundary(size_t index) const = 0;
     virtual string GetName() const { return "default"; }
+    virtual void Project(Point<3>& p) const = 0;
     // Project point using geo info. Fast if point is close to
     // parametrization in geo info.
     virtual bool ProjectPointGI(Point<3>& p, PointGeomInfo& gi) const =0;
+    virtual bool CalcPointGeomInfo(const Point<3>& p, PointGeomInfo& gi) const = 0;
     virtual Point<3> GetPoint(const PointGeomInfo& gi) const = 0;
     virtual void CalcEdgePointGI(const GeometryEdge& edge,
                                  double t,
@@ -55,6 +65,21 @@ namespace netgen
     virtual double GetCurvature(const PointGeomInfo& gi) const = 0;
 
     virtual void RestrictH(Mesh& mesh, const MeshingParameters& mparam) const = 0;
+    virtual Vec<3> GetNormal(const Point<3>& p) const
+    {
+      return {0.,0.,1.};
+    }
+    virtual Vec<3> GetNormal(const Point<3>& p, const PointGeomInfo& gi) const
+    {
+      return GetNormal(p);
+    }
+    virtual void PointBetween(const Point<3>& p1,
+                              const Point<3>& p2,
+                              double secpoint,
+                              const PointGeomInfo& gi1,
+                              const PointGeomInfo& gi2,
+                              Point<3>& newp,
+                              PointGeomInfo& newgi) const = 0;
 
   protected:
     void RestrictHTrig(Mesh& mesh,
@@ -100,21 +125,32 @@ namespace netgen
     virtual void FinalizeMesh(Mesh& mesh) const {}
 
     virtual void ProjectPoint (int surfind, Point<3> & p) const
-    { }
+    {
-    virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const { }
+      faces[surfind-1]->Project(p);
-  virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo& gi) const
+    }
-  { ProjectPointEdge(surfind, surfind2, p); }
 
-    virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const {return false;}
+    virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p) const
+    {
+      throw Exception("In ProjectPointEdge of basegeometry");
+    }
+  virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo& gi) const
+  {
+    edges[gi.edgenr]->ProjectPoint(p, gi);
+  }
+
+    virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const
+    {
+      return faces[surfind-1]->CalcPointGeomInfo(p3, gi);
+    }
     virtual bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
     {
-      throw Exception("ProjectPointGI not overloaded in class" + Demangle(typeid(*this).name()));
+      return faces[surfind-1]->ProjectPointGI(p, gi);
     }
 
     virtual Vec<3> GetNormal(int surfind, const Point<3> & p) const
-    { return {0.,0.,1.}; }
+    { return faces[surfind-1]->GetNormal(p); }
     virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo & gi) const
-    { return GetNormal(surfind, p); }
+    { return faces[surfind-1]->GetNormal(p, gi); }
     [[deprecated]]
     void GetNormal(int surfind, const Point<3> & p, Vec<3> & n) const
     {
@@ -129,6 +165,11 @@ namespace netgen
                                Point<3> & newp,
                                PointGeomInfo & newgi) const
     {
+      if(faces.Size())
+        {
+          faces[surfi-1]->PointBetween(p1, p2, secpoint, gi1, gi2, newp, newgi);
+          return;
+        }
       newp = p1 + secpoint * (p2-p1);
     }
 
@@ -140,13 +181,21 @@ namespace netgen
                                   Point<3> & newp,
                                   EdgePointGeomInfo & newgi) const
     {
+      if(edges.Size())
+        {
+          edges[ap1.edgenr]->PointBetween(p1, p2, secpoint,
+                                          ap1, ap2, newp, newgi);
+          return;
+        }
       newp = p1+secpoint*(p2-p1);
     }
 
     virtual Vec<3> GetTangent(const Point<3> & p, int surfi1,
                               int surfi2,
                               const EdgePointGeomInfo & egi) const
-    { throw Exception("Call GetTangent of " + Demangle(typeid(*this).name())); }
+    {
+      return edges[egi.edgenr]->GetTangent(p);
+    }
 
     virtual size_t GetEdgeIndex(const GeometryEdge& edge) const
     {