From bb25aa656a9c1ad82ca71cfcc007e677dda66308 Mon Sep 17 00:00:00 2001
From: Christopher Lackner <clackner@cerbsim.com>
Date: Mon, 16 Jan 2023 19:57:48 +0100
Subject: [PATCH] more stable occ parameter space projection

occ uv parameters are not very accurate sometimes, therefore new
strategy for finding same uv points in parameter space:

if uv coordinates are closer together than 0.9 * min(seg_length) and
global point numbers are the same -> same uv point
---
 libsrc/gprim/adtree.hpp     |  1 +
 libsrc/meshing/meshing2.cpp | 13 +++++++++----
 libsrc/meshing/meshing2.hpp |  3 ++-
 libsrc/occ/occgenmesh.cpp   | 31 ++++++++++++++++++++++---------
 4 files changed, 34 insertions(+), 14 deletions(-)

diff --git a/libsrc/gprim/adtree.hpp b/libsrc/gprim/adtree.hpp
index 5bbca0e6..594cf76d 100644
--- a/libsrc/gprim/adtree.hpp
+++ b/libsrc/gprim/adtree.hpp
@@ -817,6 +817,7 @@ public:
       : BoxTree(box.PMin(), box.PMax())
     { }
 
+  void SetTolerance(double _tol) { tol = _tol; }
   double GetTolerance() { return tol; }
 
   size_t GetNLeaves()
diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp
index 56209458..deadcd10 100644
--- a/libsrc/meshing/meshing2.cpp
+++ b/libsrc/meshing/meshing2.cpp
@@ -79,12 +79,17 @@ namespace netgen
   Meshing2 :: ~Meshing2 ()
   { ; } 
 
-  void Meshing2 :: AddPoint (const Point3d & p, PointIndex globind, 
-			     MultiPointGeomInfo * mgi,
-			     bool pointonsurface)
+  int Meshing2 :: AddPoint (const Point3d & p, PointIndex globind,
+                            MultiPointGeomInfo * mgi,
+                            bool pointonsurface)
   {
     //(*testout) << "add point " << globind << endl;
-    adfront.AddPoint (p, globind, mgi, pointonsurface);
+    return adfront.AddPoint (p, globind, mgi, pointonsurface);
+  }
+
+  PointIndex Meshing2 :: GetGlobalIndex(int pi) const
+  {
+    return adfront.GetGlobalIndex(pi);
   }
 
   void Meshing2 :: AddBoundaryElement (int i1, int i2,
diff --git a/libsrc/meshing/meshing2.hpp b/libsrc/meshing/meshing2.hpp
index a05cd3dd..2af29296 100644
--- a/libsrc/meshing/meshing2.hpp
+++ b/libsrc/meshing/meshing2.hpp
@@ -66,8 +66,9 @@ public:
 
 
   ///
-  DLL_HEADER void AddPoint (const Point3d & p, PointIndex globind, MultiPointGeomInfo * mgi = NULL,
+  DLL_HEADER int AddPoint (const Point3d & p, PointIndex globind, MultiPointGeomInfo * mgi = NULL,
 		 bool pointonsurface = true);
+  DLL_HEADER PointIndex GetGlobalIndex(int pi) const;
 
   ///
   DLL_HEADER void AddBoundaryElement (INDEX i1, INDEX i2,
diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp
index 2ed46b68..d2903765 100644
--- a/libsrc/occ/occgenmesh.cpp
+++ b/libsrc/occ/occgenmesh.cpp
@@ -350,6 +350,14 @@ namespace netgen
                 uv_box.Add( {seg.epgeominfo[i].u, seg.epgeominfo[i].v } );
 
         BoxTree<2> uv_tree(uv_box);
+        double tol = 1e99;
+        for(auto& seg : segments)
+          {
+            Point<2> p1 = { seg.epgeominfo[0].u, seg.epgeominfo[0].v };
+            Point<2> p2 = { seg.epgeominfo[1].u, seg.epgeominfo[1].v };
+            tol = min2(tol, Dist(p1, p2));
+          }
+        uv_tree.SetTolerance(0.9 * tol);
         Array<int> found_points;
 
         for(auto & seg : segments)
@@ -368,15 +376,21 @@ namespace netgen
                 Point<2> uv = {gi[j].u, gi[j].v};
                 uv_tree.GetIntersecting(uv, uv, found_points);
 
-                if(found_points.Size())
-                    locpnum[j] = found_points[0];
-                else
+                bool found = false;
+                for(auto& fp : found_points)
+                  {
+                    if(meshing.GetGlobalIndex(fp - 1) == seg[j])
+                      {
+                        locpnum[j] = fp;
+                        found = true;
+                      }
+                  }
+                if(!found)
                 {
                     PointIndex pi = seg[j];
-                    meshing.AddPoint (mesh.Point(pi), pi);
-                    glob2loc[pi] = ++cntpt;
+                    locpnum[j] = meshing.AddPoint (mesh.Point(pi), pi) + 1;
+                    glob2loc[pi] = locpnum[j];
                     gis.Append (gi[j]);
-                    locpnum[j] = gis.Size();
                     uv_tree.Insert(uv, locpnum[j]);
                 }
             }
@@ -391,11 +405,10 @@ namespace netgen
                 auto gi = occface.Project(mesh[pi]);
                 MultiPointGeomInfo mgi;
                 mgi.AddPointGeomInfo(gi);
-                meshing.AddPoint(mesh[pi], pi, &mgi);
+                glob2loc[pi] = meshing.AddPoint(mesh[pi], pi, &mgi) + 1;
                 gis.Append(gi);
                 Point<2> uv = { gi.u, gi.v };
-                uv_tree.Insert(uv, gis.Size());
-                glob2loc[pi] = ++cntpt;
+                uv_tree.Insert(uv, glob2loc[pi]);
               }
           }
       }