From d36b3d8b4e314c74630200c6a03bdf8231f9c966 Mon Sep 17 00:00:00 2001
From: Christopher Lackner <clackner@cerbsim.com>
Date: Wed, 30 Mar 2022 12:47:07 +0200
Subject: [PATCH] [occ] inner point of surface -> surface mesh

---
 libsrc/meshing/basegeom.cpp | 10 ++++++++++
 libsrc/meshing/basegeom.hpp |  2 ++
 libsrc/occ/occ_face.hpp     |  1 +
 libsrc/occ/occgenmesh.cpp   | 37 +++++++++++++++++++++++++++++++++++--
 libsrc/occ/occgeom.cpp      |  8 ++++++++
 libsrc/occ/occgeom.hpp      |  2 ++
 6 files changed, 58 insertions(+), 2 deletions(-)

diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp
index 4df8e44a..a340ccb6 100644
--- a/libsrc/meshing/basegeom.cpp
+++ b/libsrc/meshing/basegeom.cpp
@@ -704,6 +704,16 @@ namespace netgen
               }
           }
       }
+    for(const auto& vert : GetFaceVertices(face))
+      {
+        PointIndex pi = vert->nr + 1;
+        if(glob2loc[pi] == 0)
+          {
+            meshing.AddPoint(mesh[pi], pi);
+            cntp++;
+            glob2loc[pi] = cntp;
+          }
+      }
     for(auto & seg : segments)
       {
         PointGeomInfo gi0, gi1;
diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp
index e0bf18b6..36c923c4 100644
--- a/libsrc/meshing/basegeom.hpp
+++ b/libsrc/meshing/basegeom.hpp
@@ -203,6 +203,8 @@ namespace netgen
     const GeometryEdge & GetEdge(int i) const { return *edges[i]; }
     const GeometryVertex & GetVertex(int i) const { return *vertices[i]; }
 
+    virtual Array<GeometryVertex*> GetFaceVertices(const GeometryFace& face) const { return Array<GeometryVertex*>{}; }
+
     void Clear();
 
     virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
diff --git a/libsrc/occ/occ_face.hpp b/libsrc/occ/occ_face.hpp
index 63cc7715..ed749461 100644
--- a/libsrc/occ/occ_face.hpp
+++ b/libsrc/occ/occ_face.hpp
@@ -25,6 +25,7 @@ namespace netgen
         public:
         OCCFace(TopoDS_Shape dshape);
 
+        const TopoDS_Face Shape() const { return face; }
         T_Shape TShape() { return tface; }
 
         size_t GetHash() const override;
diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp
index d36b5d88..e0164177 100644
--- a/libsrc/occ/occgenmesh.cpp
+++ b/libsrc/occ/occgenmesh.cpp
@@ -4,6 +4,7 @@
 #include <meshing.hpp>
 
 #include "occgeom.hpp"
+#include "occ_face.hpp"
 #include "occmeshsurf.hpp"
 
 #include <BRepAdaptor_Curve.hxx>
@@ -250,6 +251,7 @@ namespace netgen
 
     FaceDescriptor & fd = mesh.GetFaceDescriptor(k);
     auto face = TopoDS::Face(geom.fmap(k));
+    const auto& occface = dynamic_cast<const OCCFace&>(geom.GetFace(k-1));
     auto fshape = face.TShape();
 
     int oldnf = mesh.GetNSE();
@@ -298,6 +300,20 @@ namespace netgen
                     glob2loc[pi] = cntp;
                   }
               }
+        for(const auto& vert : geom.GetFaceVertices(geom.GetFace(k-1)))
+          {
+            PointIndex pi = vert->nr + 1;
+            if(glob2loc[pi] == 0)
+              {
+                auto gi = occface.Project(mesh[pi]);
+                MultiPointGeomInfo mgi;
+                mgi.AddPointGeomInfo(gi);
+                meshing.AddPoint(mesh[pi], pi, &mgi);
+                cntp++;
+                glob2loc[pi] = cntp;
+              }
+          }
+
 
         /*
         for (int i = 1; i <= mesh.GetNSeg(); i++)
@@ -323,9 +339,11 @@ namespace netgen
     else
       {
         static Timer t("MeshSurface: Find edges and points - Parameter"); RegionTimer r(t);
-        
+
         Array<PointGeomInfo> gis(2*segments.Size());
         gis.SetSize (0);
+        glob2loc = 0;
+        int cntpt = 0;
 
         Box<2> uv_box(Box<2>::EMPTY_BOX);
         for(auto & seg : segments)
@@ -357,7 +375,7 @@ namespace netgen
                 {
                     PointIndex pi = seg[j];
                     meshing.AddPoint (mesh.Point(pi), pi);
-
+                    glob2loc[pi] = ++cntpt;
                     gis.Append (gi[j]);
                     locpnum[j] = gis.Size();
                     uv_tree.Insert(uv, locpnum[j]);
@@ -366,6 +384,21 @@ namespace netgen
 
             meshing.AddBoundaryElement (locpnum[0], locpnum[1], gi[0], gi[1]);
         }
+        for(const auto& vert : geom.GetFaceVertices(geom.GetFace(k-1)))
+          {
+            PointIndex pi = vert->nr + 1;
+            if(glob2loc[pi] == 0)
+              {
+                auto gi = occface.Project(mesh[pi]);
+                MultiPointGeomInfo mgi;
+                mgi.AddPointGeomInfo(gi);
+                meshing.AddPoint(mesh[pi], pi, &mgi);
+                gis.Append(gi);
+                Point<2> uv = { gi.u, gi.v };
+                uv_tree.Insert(uv, gis.Size());
+                glob2loc[pi] = ++cntpt;
+              }
+          }
       }
 
 
diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp
index a163baf3..a4800c02 100644
--- a/libsrc/occ/occgeom.cpp
+++ b/libsrc/occ/occgeom.cpp
@@ -1316,6 +1316,14 @@ namespace netgen
    }
 
 
+  Array<GeometryVertex*> OCCGeometry :: GetFaceVertices(const GeometryFace& face) const
+  {
+    Array<GeometryVertex*> verts;
+    const auto& occface = dynamic_cast<const OCCFace&>(face);
+    for(auto& vert : GetVertices(occface.Shape()))
+      verts.Append(vertices[vertex_map.at(vert.TShape())].get());
+    return move(verts);
+  }
 
 
    void OCCGeometry :: BuildVisualizationMesh (double deflection)
diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp
index 45859026..03c8d7d7 100644
--- a/libsrc/occ/occgeom.hpp
+++ b/libsrc/occ/occgeom.hpp
@@ -247,6 +247,8 @@ namespace netgen
 
     void MakeSolid();
 
+    Array<GeometryVertex*> GetFaceVertices(const GeometryFace& face) const override;
+
     void HealGeometry();
     void GlueGeometry();