From 9c1db9a6f32ea4b28842e2ab3b54ea824abde2af Mon Sep 17 00:00:00 2001
From: Christopher Lackner <clackner@cerbsim.com>
Date: Thu, 30 Mar 2023 23:21:58 +0200
Subject: [PATCH] correctly set domin and domout at sides

---
 libsrc/meshing/boundarylayer.cpp | 37 ++++++++++++++++++++++++++------
 libsrc/meshing/boundarylayer.hpp |  2 +-
 2 files changed, 32 insertions(+), 7 deletions(-)

diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp
index 55837d7d..ce837a08 100644
--- a/libsrc/meshing/boundarylayer.cpp
+++ b/libsrc/meshing/boundarylayer.cpp
@@ -613,9 +613,6 @@ namespace netgen
     surfacefacs.SetSize(nfd_old+1);
     surfacefacs = 0.0;
     // create new FaceDescriptors
-    domain_map.SetSize(mesh.GetNDomains()+1);
-    for(auto i : Range(mesh.GetNDomains()+1))
-      domain_map[i] = i;
     for(auto i : Range(1, nfd_old+1))
       {
         const auto& fd = mesh.GetFaceDescriptor(i);
@@ -629,7 +626,6 @@ namespace netgen
                 // -1 surf nr is so that curving does not do anything
                 FaceDescriptor new_fd(-1, isIn ? new_mat_nrs[i] : fd.DomainIn(),
                                       isIn ? fd.DomainOut() : new_mat_nrs[i], -1);
-                domain_map[isIn ? fd.DomainIn() : fd.DomainOut()] = new_mat_nrs[i];
                 new_fd.SetBCProperty(new_si);
                 mesh.AddFaceDescriptor(new_fd);
                 si_map[i] = new_si;
@@ -665,8 +661,9 @@ namespace netgen
             auto isIn = domains.Test(fd.DomainIn());
             auto isOut = domains.Test(fd.DomainOut());
             int si = params.sides_keep_surfaceindex ? facei : -1;
-            FaceDescriptor new_fd(si, domain_map[fd.DomainIn()],
-                                  domain_map[fd.DomainOut()], si);
+            // domin and domout can only be set later
+            FaceDescriptor new_fd(si, -1,
+                                  -1, si);
             new_fd.SetBCProperty(new_si);
             mesh.AddFaceDescriptor(new_fd);
             si_map[facei] = new_si;
@@ -1371,6 +1368,33 @@ namespace netgen
         }
   }
 
+  void BoundaryLayerTool :: SetDomInOutSides()
+  {
+    BitArray done(mesh.GetNFD()+1);
+    done.Clear();
+    for(auto sei : Range(mesh.SurfaceElements()))
+      {
+        auto& sel = mesh[sei];
+        auto index = sel.GetIndex();
+        if(done.Test(index))
+          continue;
+        done.SetBit(index);
+        auto& fd = mesh.GetFaceDescriptor(index);
+        if(fd.DomainIn() != -1)
+          continue;
+        int e1, e2;
+        mesh.GetTopology().GetSurface2VolumeElement(sei+1, e1, e2);
+        if(e1 == 0)
+          fd.SetDomainIn(0);
+        else
+          fd.SetDomainIn(mesh.VolumeElement(e1).GetIndex());
+        if(e2 == 0)
+          fd.SetDomainOut(0);
+        else
+          fd.SetDomainOut(mesh.VolumeElement(e2).GetIndex());
+      }
+  }
+
   void BoundaryLayerTool :: AddSegments()
   {
     if(have_single_segments)
@@ -1451,6 +1475,7 @@ namespace netgen
       mesh.GetTopology().ClearEdges();
       mesh.SetNextMajorTimeStamp();
       mesh.UpdateTopology();
+      SetDomInOutSides();
       MeshingParameters mp;
       mp.optimize3d ="m";
       mp.optsteps3d = 4;
diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp
index dbbbe90e..3b83dc05 100644
--- a/libsrc/meshing/boundarylayer.hpp
+++ b/libsrc/meshing/boundarylayer.hpp
@@ -45,7 +45,6 @@ class BoundaryLayerTool
     Array<SegmentIndex> moved_segs;
     int max_edge_nr, nfd_old, ndom_old;
     Array<int> new_mat_nrs;
-    Array<int> domain_map;
     BitArray moved_surfaces;
     int np, nseg, nse, ne;
     double height;
@@ -70,6 +69,7 @@ class BoundaryLayerTool
 
     void InsertNewElements(FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap, const BitArray & in_surface_direction);
     void SetDomInOut();
+    void SetDomInOutSides();
     void AddSegments();
     void FixVolumeElements();