From 46461fd166d3eefc566d98965d62bc24333fbe34 Mon Sep 17 00:00:00 2001
From: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date: Thu, 12 Dec 2024 20:11:18 +0100
Subject: [PATCH] Boundary layers - interpolate accross "inner" edges

Inner edges meaning that it's adjacent to two moved faces which have the
same normal vector (there is no kink at this edge)
---
 libsrc/meshing/boundarylayer.hpp             |  1 +
 libsrc/meshing/boundarylayer_interpolate.cpp | 50 ++++++++++++++++++--
 2 files changed, 47 insertions(+), 4 deletions(-)

diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp
index 1c2b0484..68311e97 100644
--- a/libsrc/meshing/boundarylayer.hpp
+++ b/libsrc/meshing/boundarylayer.hpp
@@ -60,6 +60,7 @@ public:
   BitArray moved_surfaces;
   int np, nseg, nse, ne;
   double total_height;
+  Array<POINTTYPE, PointIndex> point_types;
 
   // These parameters are derived from given BoundaryLayerParameters and the Mesh
   Array<double> par_heights;
diff --git a/libsrc/meshing/boundarylayer_interpolate.cpp b/libsrc/meshing/boundarylayer_interpolate.cpp
index b424c985..5537c3d1 100644
--- a/libsrc/meshing/boundarylayer_interpolate.cpp
+++ b/libsrc/meshing/boundarylayer_interpolate.cpp
@@ -72,6 +72,10 @@ BuildNeighbors (FlatArray<PointIndex> points, const Mesh& mesh)
 
 void BoundaryLayerTool ::InterpolateGrowthVectors()
 {
+  point_types.SetSize(mesh.GetNP());
+  for (auto p : mesh.Points().Range())
+    point_types[p] = mesh[p].Type();
+
   int new_max_edge_nr = max_edge_nr;
   for (const auto& seg : segments)
     if (seg.edgenr > new_max_edge_nr)
@@ -117,9 +121,49 @@ void BoundaryLayerTool ::InterpolateGrowthVectors()
     },
     mesh.GetNP());
 
+  for (auto edgenr : Range(1, new_max_edge_nr + 1))
+    {
+      // "inner" edges between two flat faces are not treated as edges for interpolation
+      bool no_angles = true;
+      ArrayMem<SurfaceElementIndex, 4> faces;
+
+      for (auto* p_seg : edgenr2seg[edgenr])
+        {
+          auto& seg = *p_seg;
+          faces.SetSize(0);
+          if (seg[0] <= p2sel.Size())
+            {
+              for (auto sei : p2sel[seg[0]])
+                if (moved_surfaces.Test(mesh[sei].GetIndex()) && p2sel[seg[1]].Contains(sei))
+                  faces.Append(sei);
+            }
+
+          if (faces.Size() == 2)
+            {
+              auto n0 = getNormal(mesh[faces[0]]);
+              auto n1 = getNormal(mesh[faces[1]]);
+              if (n0 * n1 < 0.999)
+                no_angles = false;
+            }
+          else
+            {
+              no_angles = false;
+            }
+        }
+      if (no_angles)
+        {
+          for (auto* p_seg : edgenr2seg[edgenr])
+            for (auto pi : p_seg->PNums())
+              if (pi <= np && point_types[pi] == EDGEPOINT)
+                point_types[pi] = SURFACEPOINT;
+          continue;
+        }
+    }
+
   for (auto edgenr : Range(max_edge_nr + 1, new_max_edge_nr + 1))
     {
       double edge_len = 0.;
+      bool any_grows = false;
 
       auto is_end_point = [&] (PointIndex pi) {
         auto segs = point2seg[pi];
@@ -132,8 +176,6 @@ void BoundaryLayerTool ::InterpolateGrowthVectors()
         return false;
       };
 
-      bool any_grows = false;
-
       Array<PointIndex> points;
       for (auto* p_seg : edgenr2seg[edgenr])
         {
@@ -280,7 +322,7 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors()
   for (const auto& sel : mesh.SurfaceElements())
     {
       for (auto pi : sel.PNums())
-        if (mesh[pi].Type() == SURFACEPOINT && hasMoved(pi))
+        if (point_types[pi] == SURFACEPOINT && hasMoved(pi))
           points_set.insert(pi);
     }
 
@@ -406,7 +448,7 @@ void BoundaryLayerTool ::FixSurfaceElements()
       const auto& sel = mesh[sei];
       if (sel.GetNP() == 3 && is_boundary_moved[sel.GetIndex()])
         for (auto pi : sel.PNums())
-          if (mesh[pi].Type() == SURFACEPOINT)
+          if (point_types[pi] == SURFACEPOINT)
             points_set.insert(pi);
     }