Interpolate growth vectors on surfaces

also:
- clean up delaunay2d code (Use Point<2>, remove comments)
- implement CalcWeights() used to interpolate data from boundary points
  to surface points

libsrc/gprim/adtree.hpp

@@ -1184,6 +1184,8 @@ public:
         : DelaunayTree(box.PMin(), box.PMax())
       { }
 
+    double GetTolerance() { return tol; }
+
     size_t GetNLeaves()
       {
         return n_leaves;

libsrc/meshing/CMakeLists.txt

@@ -26,7 +26,7 @@ add_library(mesh ${NG_LIB_TYPE}
         ruler2.cpp ruler3.cpp secondorder.cpp smoothing2.5.cpp		    
         smoothing2.cpp smoothing3.cpp specials.cpp
         topology.cpp validate.cpp bcfunctions.cpp   
-        parallelmesh.cpp  paralleltop.cpp  paralleltop.hpp basegeom.cpp     
+        parallelmesh.cpp  paralleltop.cpp  basegeom.cpp
         python_mesh.cpp surfacegeom.cpp
         ../../ng/onetcl.cpp
         ${rules_sources}
@@ -53,6 +53,6 @@ install(FILES
   localh.hpp meshclass.hpp meshfunc.hpp meshing2.hpp meshing3.hpp
   meshing.hpp meshtool.hpp meshtype.hpp msghandler.hpp paralleltop.hpp
   ruler2.hpp ruler3.hpp specials.hpp topology.hpp validate.hpp
-  python_mesh.hpp surfacegeom.hpp
+  python_mesh.hpp surfacegeom.hpp delaunay2d.hpp
   DESTINATION ${NG_INSTALL_DIR_INCLUDE}/meshing COMPONENT netgen_devel
 )

libsrc/meshing/boundarylayer.cpp

@@ -1,6 +1,7 @@
 #include <mystdlib.h>
 #include "meshing.hpp"
 #include "meshing2.hpp"
+#include "delaunay2d.hpp"
 #include "global.hpp"
 #include "../geom2d/csg2d.hpp"
 
@@ -262,6 +263,8 @@ namespace netgen
           {
             for(auto pi : sel.PNums())
               {
+                if(mesh[pi].Type() >= SURFACEPOINT)
+                    continue;
                 auto & np = growthvectors[pi];
                 if(np.Length() == 0) { np = n; continue; }
                 auto npn = np * n;
@@ -352,6 +355,8 @@ namespace netgen
               for(auto i : Range(sel.PNums()))
                 {
                   auto pi = sel.PNums()[i];
+                  if(mesh[pi].Type() >= SURFACEPOINT)
+                    continue;
                   if(growthvectors[pi].Length2() == 0.)
                     continue;
                   auto next = sel.PNums()[(i+1)%sel.GetNV()];
@@ -393,6 +398,73 @@ namespace netgen
           }
       }
 
+    Array<Point<2>, PointIndex> delaunay_points(mesh.GetNP());
+    Array<int, PointIndex> p2face(mesh.GetNP());
+    p2face = 0;
+
+    // interpolate growth vectors at inner surface points from surrounding edge points
+    Array<SurfaceElementIndex> surface_els;
+    Array<PointIndex> edge_points;
+    Array<PointIndex> surface_points;
+    for(auto facei : Range(1, fd_old+1))
+      {
+        if(!blp.surfid.Contains(facei))
+            continue;
+
+        p2face = 0;
+
+        edge_points.SetSize(0);
+        surface_points.SetSize(0);
+        surface_els.SetSize(0);
+        mesh.GetSurfaceElementsOfFace (facei, surface_els);
+        Box<2> bbox ( Box<2>::EMPTY_BOX );
+        for(auto sei : surface_els)
+        {
+            const auto & sel = mesh[sei];
+            for (auto i : Range(sel.GetNP()))
+            {
+                auto pi = sel[i];
+                if(p2face[pi] != 0)
+                    continue;
+
+                p2face[pi] = facei;
+
+                if(mesh[pi].Type() <= EDGEPOINT)
+                    edge_points.Append(pi);
+                else
+                    surface_points.Append(pi);
+
+                auto & gi = sel.GeomInfo()[i];
+                // TODO: project to plane if u,v not available?
+                delaunay_points[pi] = {gi.u, gi.v};
+                bbox.Add(delaunay_points[pi]);
+            }
+        }
+
+        if(surface_points.Size()==0)
+            continue;
+
+        DelaunayMesh dmesh( delaunay_points, bbox  );
+
+        for(auto pi : edge_points)
+        {
+            p2face[pi] = 0;
+            dmesh.AddPoint(pi);
+        }
+
+        std::map<PointIndex, double> weights;
+        for(auto pi : surface_points)
+        {
+            dmesh.AddPoint(pi, &weights);
+            auto & v = growthvectors[pi];
+            for(auto & [pi_other, weight] : weights)
+                v += weight * growthvectors[pi_other];
+        }
+
+      }
+
+
+
     // insert new points
     for (PointIndex pi = 1; pi <= np; pi++)
       if (growthvectors[pi].Length2() != 0)

libsrc/meshing/delaunay2d.cpp

@@ -1,69 +1,22 @@
 #include <mystdlib.h>
-#include "meshing.hpp"
+#include "delaunay2d.hpp"
 
 #include <geom2d/csg2d.hpp>
 
 
-// not yet working ....
-
 namespace netgen
 {
   using ngcore::INT;
-  extern void Optimize2d (Mesh & mesh, MeshingParameters & mp);
 
-  static inline Point<2> P2( Point<3> p )
+
+  void DelaunayTrig::CalcCenter (FlatArray<Point<2>, PointIndex> points)
   {
-    return {p[0], p[1]};
+    Point<2> p1 = points[pnums[0]];
-  }
+    Point<2> p2 = points[pnums[1]];
-
+    Point<2> p3 = points[pnums[2]];
-  static inline Point<3> P3( Point<2> p )
-  {
-    return {p[0], p[1], 0};
-  }
-
-
-  class DelaunayTrig
-  {
-    PointIndex pnums[3];
-    Point<2> c;
-  public:    
-    double r;
-    double rad2;
-    DelaunayTrig () = default;
-    DelaunayTrig (int p1, int p2, int p3)
-    { 
-      pnums[0] = p1;
-      pnums[1] = p2;
-      pnums[2] = p3;
-    }
-
-    PointIndex & operator[] (int j) { return pnums[j]; }
-    const PointIndex & operator[] (int j) const { return pnums[j]; }
-
-    void CalcCenter (Mesh & mesh)
-    {
-      Point<2> p1 = P2(mesh[pnums[0]]);
-      Point<2> p2 = P2(mesh[pnums[1]]);
-      Point<2> p3 = P2(mesh[pnums[2]]);
 
     Vec<2> v1 = p2-p1;
     Vec<2> v2 = p3-p1;
-      /*
-      Mat<2,2> mat, inv;
-      mat(0,0) = v1*v1;
-      mat(0,1) = v1*v2;
-      mat(1,0) = v2*v1;
-      mat(1,1) = v2*v2;
-      Vec<2> rhs, sol;
-      rhs(0) = 0.5 * v1*v1;
-      rhs(1) = 0.5 * v2*v2;
-      CalcInverse (mat, inv);
-      sol = inv * rhs;
-      
-      c = p1 + sol(0) * v1 + sol(1) * v2;
-      rad2 = Dist2(c, p1);
-      r = sqrt(rad2);
-      */
 
     // without normal equation ...
     Mat<2,2> mat, inv;
@@ -80,25 +33,7 @@ namespace netgen
     r = sqrt(rad2);
   }
 
-    Point<2> Center() const { return c; }
+  int DelaunayMesh::GetNeighbour( int eli, int edge )
-    double Radius2() const { return rad2; }
-    Box<2> BoundingBox() const { return Box<2> (c-Vec<2>(r,r), c+Vec<2>(r,r)); }
-
-    mutable PointIndex visited_pi = -1;
-  };
-
-  class DelaunayMesh
-  {
-    ngcore::ClosedHashTable<INT<2>, INT<2>> edge_to_trig;
-    Array<DelaunayTrig> trigs;
-    unique_ptr<DelaunayTree<2>> tree;
-    Mesh & mesh;
-
-    Array<int> closeels;
-    Array<int> intersecting;
-    Array<INT<2>> edges;
-
-    int GetNeighbour( int eli, int edge )
   {
     auto p0 = trigs[eli][(edge+1)%3];
     auto p1 = trigs[eli][(edge+2)%3];
@@ -106,21 +41,14 @@ namespace netgen
       Swap(p0,p1);
 
     INT<2> hash = {p0,p1};
-      /*
-      if(!edge_to_trig.Used(hash))
-        return -1;
 
-      auto i2 = edge_to_trig.Get({p0,p1});
-
-      return i2[0] == eli ? i2[1] : i2[0];
-      */
     auto pos = edge_to_trig.Position(hash);
     if (pos == -1) return -1;
     auto i2 = edge_to_trig.GetData(pos);
     return i2[0] == eli ? i2[1] : i2[0];
   }
 
-    void SetNeighbour( int eli, int edge )
+  void DelaunayMesh::SetNeighbour( int eli, int edge )
   {
     auto p0 = trigs[eli][(edge+1)%3];
     auto p1 = trigs[eli][(edge+2)%3];
@@ -143,28 +71,9 @@ namespace netgen
           }
         edge_to_trig.SetData (pos, i2);
       }
-      /*
-      if(!edge_to_trig.Used(hash))
-        edge_to_trig[hash] = {eli, -1};
-      else
-      {
-
-        auto i2 = edge_to_trig.Get({p0,p1});
-
-        if(i2[0]==-1)
-          i2[0] = eli;
-        else
-        {
-          if(i2[1]==-1)
-            i2[1] = eli;
   }
 
-        edge_to_trig[hash] = i2;
+  void DelaunayMesh::UnsetNeighbours( int eli )
-      }
-      */
-    }
-
-    void UnsetNeighbours( int eli )
   {
     for(int edge : Range(3))
     {
@@ -174,7 +83,6 @@ namespace netgen
         Swap(p0,p1);
 
       INT<2> hash = {p0,p1};
-        // auto i2 = edge_to_trig.Get({p0,p1});
       auto pos = edge_to_trig.Position(hash);
       auto i2 = edge_to_trig.GetData(pos);
 
@@ -182,20 +90,19 @@ namespace netgen
         i2[0] = i2[1];
       i2[1] = -1;
 
-        // edge_to_trig[hash] = i2;
       edge_to_trig.SetData (pos, i2);
     }
   }
 
 
-    void AppendTrig( int pi0, int pi1, int pi2 )
+  void DelaunayMesh::AppendTrig( int pi0, int pi1, int pi2 )
   {
     DelaunayTrig el;
     el[0] = pi0;
     el[1] = pi1;
     el[2] = pi2;
 
-      el.CalcCenter(mesh);
+    el.CalcCenter(points);
 
     trigs.Append(el);
     int ti = trigs.Size()-1;
@@ -205,9 +112,8 @@ namespace netgen
       SetNeighbour(ti, i);
   }
 
-    public:
+  DelaunayMesh::DelaunayMesh( Array<Point<2>, PointIndex> & points_, Box<2> box  )
-    DelaunayMesh( Mesh & mesh_, Box<2> box  ) 
+    : points(points_)
-      : mesh(mesh_)
   {
     Vec<2> vdiag = box.PMax()-box.PMin();
 
@@ -224,16 +130,17 @@ namespace netgen
 
     tree = make_unique<DelaunayTree<2>>(box);
 
-      auto pi0 = mesh.AddPoint (P3(p0));
+    auto pi0 = points.Append (p0);
-      auto pi1 = mesh.AddPoint (P3(p1));
+    auto pi1 = points.Append (p1);
-      auto pi2 = mesh.AddPoint (P3(p2));
+    auto pi2 = points.Append (p2);
     AppendTrig(pi0, pi1, pi2);
   }
 
-    void AddPoint( PointIndex pi_new )
+  void DelaunayMesh::CalcIntersecting( PointIndex pi_new )
   {
-      static Timer t("AddPoint"); RegionTimer reg(t);
+    static Timer t("CalcIntersecting"); RegionTimer reg(t);
-      Point<2> newp = P2(mesh[pi_new]);
+
+    Point<2> newp = points[pi_new];
     intersecting.SetSize(0);
     edges.SetSize(0);
 
@@ -284,18 +191,34 @@ namespace netgen
       }
     }
 
-    // static Timer tvis("trig visited");
-    // tvis.Start();
-    // BitArray trig_visited(trigs.Size());
-    // trig_visited.Clear();
     if(definitive_overlapping_trig==-1)
+    {
+      Mesh m;
+      m.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0));
+      for(auto pi : points.Range())
+          m.AddPoint(P3(points[pi]));
+
+      for (DelaunayTrig & trig : trigs)
+      {
+        if (trig[0] < 0) continue;
+
+        Vec<3> n = Cross (P3(points[trig[1]])-P3(points[trig[0]]), 
+            P3(points[trig[2]])-P3(points[trig[0]]));
+        if (n(2) < 0) Swap (trig[1], trig[2]);
+
+        Element2d el(trig[0], trig[1], trig[2]);
+        el.SetIndex (1);
+        m.AddSurfaceElement (el);
+      }
+      m.Compress();
+      m.AddPoint(P3(points[pi_new]));
+      m.Save("error.vol.gz");
       throw Exception("point not in any circle "+ ToString(pi_new));
-    // tvis.Stop();
+    }
-    // static Timer t2("addpoint - rest"); RegionTimer r2(t2);
+
     Array<int> trigs_to_visit;
     trigs_to_visit.Append(definitive_overlapping_trig);
     intersecting.Append(definitive_overlapping_trig);
-    // trig_visited.SetBit(definitive_overlapping_trig);
     trigs[definitive_overlapping_trig].visited_pi = pi_new;
 
     while(trigs_to_visit.Size())
@@ -303,8 +226,6 @@ namespace netgen
       int ti = trigs_to_visit.Last();
       trigs_to_visit.DeleteLast();
 
-      // trig_visited.SetBit(ti);
-
       auto & trig = trigs[ti];
       trig.visited_pi = pi_new;
 
@@ -313,23 +234,20 @@ namespace netgen
         auto nb = GetNeighbour(ti, ei);
         if(nb==-1)
           continue;
-        // if(trig_visited.Test(nb))
-        // continue;
 
         const auto & trig_nb = trigs[nb];
         if (trig_nb.visited_pi == pi_new)
           continue;
 
-        // trig_visited.SetBit(nb);
         trig_nb.visited_pi = pi_new;
 
         bool is_intersecting = Dist2(newp, trig_nb.Center()) < trig_nb.Radius2()*(1+1e-12);
 
         if(!is_intersecting)
         {
-          const Point<2> p0 = P2(mesh[PointIndex (trig[(ei+1)%3])]);
+          const Point<2> p0 = points[PointIndex (trig[(ei+1)%3])];
-          const Point<2> p1 = P2(mesh[PointIndex (trig[(ei+2)%3])]);
+          const Point<2> p1 = points[PointIndex (trig[(ei+2)%3])];
-          const Point<2> p2 = P2(mesh[PointIndex (trig[ei])]);
+          const Point<2> p2 = points[PointIndex (trig[ei])];
           auto v = p1-p0;
 
           Vec<2> n = {-v[1], v[0]};
@@ -373,6 +291,41 @@ namespace netgen
           edges.Append (edge);
       }
     }
+  }
+
+  void DelaunayMesh::CalcWeights( PointIndex pi_new, std::map<PointIndex, double> & weights )
+  {
+      double eps = tree->GetTolerance();
+      weights.clear();
+      double sum = 0.0;
+      auto p = points[pi_new];
+      auto pi_last = *points.Range().end()-3;
+      for(auto edge : edges)
+      {
+          for(PointIndex pi : {edge[0], edge[1]})
+          {
+          if(pi>=pi_last)
+              continue;
+          if(weights.count(pi))
+              continue;
+          double weight = 1.0/(eps+Dist(p, points[pi]));
+          sum += weight;
+          weights[pi] = weight;
+          }
+      }
+      double isum = 1.0/sum;
+      for(auto & [pi, weight] : weights)
+          weight *= isum;
+  }
+
+  void DelaunayMesh::AddPoint( PointIndex pi_new, std::map<PointIndex, double> * weights )
+  {
+      static Timer t("AddPoint"); RegionTimer reg(t);
+
+      CalcIntersecting(pi_new);
+
+      if(weights)
+          CalcWeights(pi_new, *weights);
 
       for (int j : intersecting)
         {
@@ -387,33 +340,7 @@ namespace netgen
 
       for (int j : intersecting)
           tree->DeleteElement (j);
-
-    static int counter=0;
-    if(0)
-    {
-      Mesh m;
-      m = mesh;
-      m.ClearSurfaceElements();
-      for (DelaunayTrig & trig : trigs)
-      {
-        if (trig[0] < 0) continue;
-
-        Vec<3> n = Cross (mesh[trig[1]]-mesh[trig[0]], 
-            mesh[trig[2]]-mesh[trig[0]]);
-        if (n(2) < 0) Swap (trig[1], trig[2]);
-
-        Element2d el(trig[0], trig[1], trig[2]);
-        el.SetIndex (1);
-        m.AddSurfaceElement (el);
   }
-      m.Save("meshes/mesh_"+ToString(counter++)+".vol.gz");
-    }
-
-    }
-
-    Array<DelaunayTrig> & GetElements() { return trigs; }
-
-  };
 
   ostream & operator<< (ostream & ost, DelaunayTrig trig)
   {
@@ -699,26 +626,25 @@ namespace netgen
     Array<PointIndex, PointIndex> compress;
     Array<PointIndex, PointIndex> icompress(mesh.Points().Size());
 
-    /*
+    Array<Point<2>, PointIndex> temp_points;
-    for(auto pi : mesh.Points().Range())
-      if(add_point.Test(pi))
-    */
     for (PointIndex pi : addpoints)
       {
         icompress[pi] = tempmesh.AddPoint(mesh[pi]);
         compress.Append(pi);
+        temp_points.Append(P2(mesh[pi]));
       }
 
     for (PointIndex pi : Range(first_point_blockfill, last_point_blockfill))
       {
         icompress[pi] = tempmesh.AddPoint(mesh[pi]);
         compress.Append(pi);
+        temp_points.Append(P2(mesh[pi]));
       }
     t3.Stop();
     // DelaunayMesh adds surrounding trig (don't add the last 3 points to delaunay AGAIN!
-    auto tempmesh_points = tempmesh.Points().Range();
+    auto points_range = temp_points.Range();
 
-    DelaunayMesh dmesh(tempmesh, bbox);
+    DelaunayMesh dmesh(temp_points, bbox);
 
     timer_addpoints.Start();
 
@@ -736,7 +662,7 @@ namespace netgen
 //     for (PointIndex pi : old_points)
 //       mixed[pi] = PointIndex ( (prim * pi) % old_points.Size() + PointIndex::BASE );
     
-    for (auto pi : tempmesh_points)
+    for (auto pi : points_range)
       dmesh.AddPoint(pi);
 
     timer_addpoints.Stop();

libsrc/meshing/delaunay2d.hpp

@@ -0,0 +1,72 @@
+#include "meshing.hpp"
+
+namespace netgen
+{
+  static inline Point<2> P2( Point<3> p )
+  {
+    return {p[0], p[1]};
+  }
+
+  static inline Point<3> P3( Point<2> p )
+  {
+    return {p[0], p[1], 0};
+  }
+
+  class DelaunayTrig
+  {
+    PointIndex pnums[3];
+    Point<2> c;
+
+  public:
+    double r;
+    double rad2;
+    DelaunayTrig () = default;
+    DelaunayTrig (int p1, int p2, int p3)
+    {
+      pnums[0] = p1;
+      pnums[1] = p2;
+      pnums[2] = p3;
+    }
+
+    PointIndex & operator[] (int j) { return pnums[j]; }
+    const PointIndex & operator[] (int j) const { return pnums[j]; }
+
+    void CalcCenter (FlatArray<Point<2>, PointIndex> points);
+
+    Point<2> Center() const { return c; }
+    double Radius2() const { return rad2; }
+    Box<2> BoundingBox() const { return Box<2> (c-Vec<2>(r,r), c+Vec<2>(r,r)); }
+
+    mutable PointIndex visited_pi = -1;
+  };
+
+  class DelaunayMesh
+  {
+    ngcore::ClosedHashTable<INT<2>, INT<2>> edge_to_trig;
+    Array<DelaunayTrig> trigs;
+    unique_ptr<DelaunayTree<2>> tree;
+    Array<Point<2>, PointIndex> & points;
+
+    Array<int> closeels;
+    Array<int> intersecting;
+    Array<INT<2>> edges;
+
+    int GetNeighbour( int eli, int edge );
+
+    void SetNeighbour( int eli, int edge );
+
+    void UnsetNeighbours( int eli );
+
+    void AppendTrig( int pi0, int pi1, int pi2 );
+
+    public:
+    DelaunayMesh( Array<Point<2>, PointIndex> & points_, Box<2> box  );
+
+    void CalcIntersecting( PointIndex pi_new );
+    void CalcWeights( PointIndex pi_new, std::map<PointIndex, double> & weights );
+    void AddPoint( PointIndex pi_new, std::map<PointIndex, double> * weights = nullptr );
+    Array<DelaunayTrig> & GetElements() { return trigs; }
+
+  };
+
+} // namespace netgen