From 803eb73d2d958aeaace749af89e38057eb0bada6 Mon Sep 17 00:00:00 2001
From: Christopher Lackner <clackner@cerbsim.com>
Date: Tue, 13 Feb 2024 09:35:44 +0100
Subject: [PATCH] allow internal edges on boundarylayer (for sphere, ellipsoid)

---
 libsrc/meshing/boundarylayer.cpp | 17 +++++++++++------
 libsrc/meshing/boundarylayer.hpp |  1 +
 libsrc/meshing/python_mesh.cpp   |  5 ++++-
 3 files changed, 16 insertions(+), 7 deletions(-)

diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp
index e54464d5..e5a927e9 100644
--- a/libsrc/meshing/boundarylayer.cpp
+++ b/libsrc/meshing/boundarylayer.cpp
@@ -679,7 +679,8 @@ namespace netgen
             // result in pushed through elements, since we do not (yet)
             // curvature through layers.
             // Therefore we disable curving for these surfaces.
-            mesh.GetFaceDescriptor(i).SetSurfNr(-1);
+            if(!params.keep_surfaceindex)
+              mesh.GetFaceDescriptor(i).SetSurfNr(-1);
           }
       }
 
@@ -984,7 +985,7 @@ namespace netgen
         Array<PointIndex> points;
         // find first vertex on edge
         double edge_len = 0.;
-        auto is_end_point = [&] (PointIndex pi)
+        auto is_end_point = [&] (PointIndex pi, const Segment& testseg)
         {
           // if(mesh[pi].Type() == FIXEDPOINT)
           //   return true;
@@ -992,13 +993,14 @@ namespace netgen
             auto segs = topo.GetVertexSegments(pi);
             auto first_edgenr = mesh[segs[0]].edgenr;
             for(auto segi : segs)
-                if(mesh[segi].edgenr != first_edgenr)
-                    return true;
+              if(auto& seg = mesh[segi]; seg.edgenr != first_edgenr || (testseg[0] == seg[1] && testseg[1] == seg[0]))
+                return true;
             return false;
         };
 
         bool any_grows = false;
         
+        const Segment* last_seg = nullptr;
         for(const auto& seg : segments)
           {
             if(seg.edgenr-1 == edgenr)
@@ -1006,8 +1008,9 @@ namespace netgen
                 if(growthvectors[seg[0]].Length2() != 0 ||
                    growthvectors[seg[1]].Length2() != 0)
                   any_grows = true;
-                if(points.Size() == 0 && is_end_point(seg[0]))
+                if(points.Size() == 0 && is_end_point(seg[0], seg))
                   {
+                    last_seg = &seg;
                     points.Append(seg[0]);
                     points.Append(seg[1]);
                     edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length();
@@ -1034,6 +1037,7 @@ namespace netgen
                     edge_len += (mesh[points.Last()] - mesh[seg[1]]).Length();
                     points.Append(seg[1]);
                     point_found = true;
+                    last_seg = &seg;
                     break;
                   }
                 else if(seg[1] == points.Last() &&
@@ -1042,10 +1046,11 @@ namespace netgen
                     edge_len += (mesh[points.Last()] - mesh[seg[0]]).Length();
                     points.Append(seg[0]);
                     point_found = true;
+                    last_seg = &seg;
                     break;
                   }
               }
-            if(is_end_point(points.Last()))
+            if(last_seg && is_end_point(points.Last(), *last_seg))
               break;
             if(!point_found)
               {
diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp
index 87790a2d..c879e7bb 100644
--- a/libsrc/meshing/boundarylayer.hpp
+++ b/libsrc/meshing/boundarylayer.hpp
@@ -23,6 +23,7 @@ public:
   bool limit_growth_vectors = true;
   double limit_safety = 0.3; // alloow only 30% of the growth vector length
   bool sides_keep_surfaceindex = false;
+  bool keep_surfaceindex = false;
   Array<size_t> project_boundaries;
 };
 
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 886122e6..910a5e13 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -1463,7 +1463,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                               variant<string, int> domain, bool outside,
                               optional<string> project_boundaries,
                               bool grow_edges, bool limit_growth_vectors,
-                              bool sides_keep_surfaceindex)
+                              bool sides_keep_surfaceindex,
+                              bool keep_surfaceindex)
            {
              BoundaryLayerParameters blp;
              BitArray boundaries(self.GetNFD()+1);
@@ -1549,12 +1550,14 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
              blp.grow_edges = grow_edges;
              blp.limit_growth_vectors = limit_growth_vectors;
              blp.sides_keep_surfaceindex = sides_keep_surfaceindex;
+             blp.keep_surfaceindex = keep_surfaceindex;
 
              GenerateBoundaryLayer (self, blp);
              self.UpdateTopology();
            }, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
           py::arg("domains") = ".*", py::arg("outside") = false,
           py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false,
+          py::arg("keep_surfaceindex")=false,
           R"delimiter(
 Add boundary layer to mesh.