From bc194027a282606043e9c4b2dd26aedd424b8d4a Mon Sep 17 00:00:00 2001
From: Christopher Lackner <clackner@cerbsim.com>
Date: Tue, 4 Mar 2025 09:59:56 +0100
Subject: [PATCH] fix edges and new domains in BoundaryLayer2d with
 make_new_domain=True

---
 libsrc/meshing/boundarylayer2d.cpp | 26 +++++++++++---------------
 1 file changed, 11 insertions(+), 15 deletions(-)

diff --git a/libsrc/meshing/boundarylayer2d.cpp b/libsrc/meshing/boundarylayer2d.cpp
index 8a5ae979..bb5ab6e6 100644
--- a/libsrc/meshing/boundarylayer2d.cpp
+++ b/libsrc/meshing/boundarylayer2d.cpp
@@ -243,10 +243,6 @@ namespace netgen
 
      Array<SegmentIndex> segments;
 
-    // surface index map
-    Array<int> si_map(mesh.GetNFD()+2);
-    si_map = -1;
-
     // int fd_old = mesh.GetNFD();
 
     int max_edge_nr = -1;
@@ -254,8 +250,8 @@ namespace netgen
 
     for(const auto& seg : line_segments)
     {
-      if(seg.epgeominfo[0].edgenr > max_edge_nr)
-        max_edge_nr = seg.epgeominfo[0].edgenr;
+      if(seg.edgenr > max_edge_nr)
+        max_edge_nr = seg.edgenr;
       if(seg.domin > max_domain)
          max_domain = seg.domin;
       if(seg.domout > max_domain)
@@ -263,6 +259,7 @@ namespace netgen
     }
 
     int new_domain = max_domain+1;
+    int new_edge_nr = max_edge_nr+1;
 
     BitArray active_boundaries(max_edge_nr+1);
     BitArray active_segments(nseg);
@@ -288,7 +285,10 @@ namespace netgen
         // int new_fd_index =
         mesh.AddFaceDescriptor(new_fd);
         if(should_make_new_domain)
-           mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(domain-1));
+          {
+           mesh.SetMaterial(new_domain, "layer_" + mesh.GetMaterial(domain));
+           mesh.SetBCName(new_edge_nr - 1, "moved");
+          }
     }
 
     for(auto segi : Range(line_segments))
@@ -618,8 +618,6 @@ namespace netgen
            }
         }
 
-     map<pair<PointIndex, PointIndex>, int> seg2edge;
-
      // insert new elements ( and move old ones )
      for(auto si : moved_segs)
      {
@@ -629,8 +627,6 @@ namespace netgen
         auto & pm0 = mapto[seg[0]];
         auto & pm1 = mapto[seg[1]];
 
-        // auto newindex = si_map[domain];
-
         Segment s = seg;
         s.geominfo[0] = {};
         s.geominfo[1] = {};
@@ -638,10 +634,10 @@ namespace netgen
         s[1] = pm1.Last();
         s[2] = PointIndex::INVALID;
         auto pair = s[0] < s[1] ? make_pair(s[0], s[1]) : make_pair(s[1], s[0]);
-        if(seg2edge.find(pair) == seg2edge.end())
-           seg2edge[pair] = ++max_edge_nr;
-        s.edgenr = seg2edge[pair];
-        s.si = seg.si;
+        s.edgenr = new_edge_nr;
+        s.epgeominfo[0].edgenr = -1;
+        s.epgeominfo[1].edgenr = -1;
+        s.si = s.edgenr;
         mesh.AddSegment(s);
 
         for ( auto i : Range(thicknesses))