From 9601f70c17a78cefc35c7b9490ad8956ba8e2ec5 Mon Sep 17 00:00:00 2001
From: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date: Thu, 27 Feb 2025 17:27:54 +0100
Subject: [PATCH] Keep free segments when generating boundary layers, also
 revert seg.si back to seg.edgenr+1

---
 libsrc/meshing/boundarylayer.cpp | 40 ++++++++++++++++++--------------
 libsrc/meshing/boundarylayer.hpp |  2 +-
 2 files changed, 24 insertions(+), 18 deletions(-)

diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp
index 4611a95f..be27b290 100644
--- a/libsrc/meshing/boundarylayer.cpp
+++ b/libsrc/meshing/boundarylayer.cpp
@@ -220,9 +220,8 @@ bool HaveSingleSegments (const Mesh& mesh)
 
 // duplicates segments (and sets seg.si accordingly) to have a unified data
 // structure for all geometry types
-Array<Segment> BuildSegments (Mesh& mesh)
+void BuildSegments (Mesh& mesh, bool have_single_segments, Array<Segment> & segments, Array<Segment> & free_segments)
 {
-  Array<Segment> segments;
   // auto& topo = mesh.GetTopology();
 
   NgArray<SurfaceElementIndex> surf_els;
@@ -230,6 +229,16 @@ Array<Segment> BuildSegments (Mesh& mesh)
   for (auto segi : Range(mesh.LineSegments()))
     {
       auto seg = mesh[segi];
+      if (seg.domin == seg.domout && seg.domin > 0)
+        {
+          free_segments.Append(seg);
+          continue;
+        }
+      if(!have_single_segments)
+      {
+        segments.Append(seg);
+        continue;
+      }
       mesh.GetTopology().GetSegmentSurfaceElements(segi + 1, surf_els);
       for (auto seli : surf_els)
         {
@@ -250,7 +259,6 @@ Array<Segment> BuildSegments (Mesh& mesh)
           segments.Append(seg);
         }
     }
-  return segments;
 }
 
 void MergeAndAddSegments (Mesh& mesh, FlatArray<Segment> segments, FlatArray<Segment> new_segments)
@@ -259,11 +267,11 @@ void MergeAndAddSegments (Mesh& mesh, FlatArray<Segment> segments, FlatArray<Seg
 
   mesh.LineSegments().SetSize0();
 
-  auto addSegment = [&] (const auto& seg) {
-    PointIndices<2> i2(seg[0], seg[1]);
-    i2.Sort();
+  auto addSegment = [&] (auto seg) {
+    SortedPointIndices<2> i2(seg[0], seg[1]);
     if (!already_added.Used(i2))
       {
+        seg.si = seg.edgenr + 1;
         mesh.AddSegment(seg);
         already_added.Set(i2, true);
       }
@@ -292,10 +300,7 @@ BoundaryLayerTool::BoundaryLayerTool(Mesh& mesh_,
   old_segments = mesh.LineSegments();
   have_single_segments = HaveSingleSegments(mesh);
 
-  if (have_single_segments)
-    segments = BuildSegments(mesh);
-  else
-    segments = mesh.LineSegments();
+  BuildSegments(mesh, have_single_segments, segments, free_segments);
 
   np = mesh.GetNP();
   first_new_pi = IndexBASE<PointIndex>() + np;
@@ -1157,8 +1162,11 @@ void BoundaryLayerTool ::AddSegments()
 
   if (params.disable_curving)
     {
-      for (auto& seg : old_segments)
-        if (mapto[seg[0]].Size() || mapto[seg[1]].Size())
+      auto is_mapped = [&] (PointIndex pi) {
+        return pi >= mapto.Range().Next() || mapto[pi].Size() > 0;
+      };
+      for (auto& seg : segments)
+        if (is_mapped(seg[0]) || is_mapped(seg[1]))
           {
             seg.epgeominfo[0].edgenr = -1;
             seg.epgeominfo[1].edgenr = -1;
@@ -1186,6 +1194,9 @@ void BoundaryLayerTool ::AddSegments()
       for (auto& seg : new_segs)
         mesh.AddSegment(seg);
     }
+
+  for (auto& seg : free_segments)
+    mesh.AddSegment(seg);
 }
 
 void BoundaryLayerTool ::AddSurfaceElements()
@@ -1399,11 +1410,6 @@ void BoundaryLayerTool ::Perform()
       mesh[pi] += height * (*gw);
     }
 
-  if (insert_only_volume_elements)
-    {
-      mesh.LineSegments() = old_segments;
-    }
-
   auto& identifications = mesh.GetIdentifications();
   NgArray<INDEX_2> pairs;
   for (auto nr : Range(0, identifications.GetMaxNr() + 1))
diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp
index d504e2d2..5a315cc2 100644
--- a/libsrc/meshing/boundarylayer.hpp
+++ b/libsrc/meshing/boundarylayer.hpp
@@ -71,7 +71,7 @@ public:
   Array<size_t> par_project_boundaries;
 
   bool have_single_segments;
-  Array<Segment> old_segments, segments, new_segments, new_segments_on_moved_bnd;
+  Array<Segment> old_segments, free_segments, segments, new_segments, new_segments_on_moved_bnd;
   Array<Element2d, SurfaceElementIndex> new_sels, new_sels_on_moved_bnd;
   Array<Array<PointIndex>, PointIndex> mapto;
   Array<PointIndex, PointIndex> mapfrom;