From 7a86aae0d158af439743f1d8e7c69960feff0f43 Mon Sep 17 00:00:00 2001
From: "Lackner, Christopher" <clackner@cerbsim.com>
Date: Thu, 10 Nov 2022 14:35:58 +0100
Subject: [PATCH] allow different materials in boundarylayer depending on bc
 name

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

diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp
index 863110c3..a9e31dc5 100644
--- a/libsrc/meshing/boundarylayer.cpp
+++ b/libsrc/meshing/boundarylayer.cpp
@@ -7,6 +7,7 @@
 #include "../geom2d/csg2d.hpp"
 
 #include <set>
+#include <regex>
 
 namespace netgen
 {
@@ -562,8 +563,22 @@ namespace netgen
       if(seg.edgenr > max_edge_nr)
         max_edge_nr = seg.edgenr;
 
-    new_mat_nr = mesh.GetNDomains() +1;
-    mesh.SetMaterial(new_mat_nr, params.new_mat);
+    int ndom = mesh.GetNDomains();
+    ndom_old = ndom;
+
+    new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1);
+    new_mat_nrs = -1;
+    for(auto [bcname, matname] : params.new_mat)
+      {
+        mesh.SetMaterial(++ndom, matname);
+        regex pattern(bcname);
+        for(auto i : Range(1, mesh.GetNFD()+1))
+          {
+            auto& fd = mesh.GetFaceDescriptor(i);
+            if(regex_match(fd.GetBCName(), pattern))
+              new_mat_nrs[i] = ndom;
+          }
+      }
 
     domains = params.domains;
     if(!params.outside)
@@ -606,8 +621,8 @@ namespace netgen
                 int new_si = mesh.GetNFD()+1;
                 surfacefacs[i] = isIn ? 1. : -1.;
                 // -1 surf nr is so that curving does not do anything
-                FaceDescriptor new_fd(-1, isIn ? new_mat_nr : fd.DomainIn(),
-                                      isIn ? fd.DomainOut() : new_mat_nr, -1);
+                FaceDescriptor new_fd(-1, isIn ? new_mat_nrs[i] : fd.DomainIn(),
+                                      isIn ? fd.DomainOut() : new_mat_nrs[i], -1);
                 new_fd.SetBCProperty(new_si);
                 mesh.AddFaceDescriptor(new_fd);
                 si_map[i] = new_si;
@@ -1050,7 +1065,7 @@ namespace netgen
                 if(surfacefacs[sel.GetIndex()] > 0) Swap(points[0], points[2]);
                 for(auto i : Range(points))
                   el[sel.PNums().Size() + i] = points[i];
-                el.SetIndex(new_mat_nr);
+                el.SetIndex(new_mat_nrs[sel.GetIndex()]);
                 mesh.AddVolumeElement(el);
               }
             Element2d newel = sel;
@@ -1238,10 +1253,10 @@ namespace netgen
     for(auto i : Range(1, nfd_old+1))
       if(si_map[i] != -1)
         {
-          if(mesh.GetFaceDescriptor(mesh.GetNFD()).DomainIn() == new_mat_nr)
-            mesh.GetFaceDescriptor(i).SetDomainOut(new_mat_nr);
+          if(auto dom = mesh.GetFaceDescriptor(si_map[i]).DomainIn(); dom > ndom_old)
+            mesh.GetFaceDescriptor(i).SetDomainOut(dom);
           else
-            mesh.GetFaceDescriptor(i).SetDomainIn(new_mat_nr);
+            mesh.GetFaceDescriptor(i).SetDomainIn(mesh.GetFaceDescriptor(si_map[i]).DomainOut());
         }
   }
 
diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp
index bd51718f..8e83cc7a 100644
--- a/libsrc/meshing/boundarylayer.hpp
+++ b/libsrc/meshing/boundarylayer.hpp
@@ -14,7 +14,7 @@ public:
   // parameters by Philippose ..
   Array<int> surfid;
   Array<double> heights;
-  string new_mat;
+  map<string, string> new_mat;
   BitArray domains;
   bool outside = false; // set the boundary layer on the outside
   bool grow_edges = false;
@@ -42,7 +42,8 @@ class BoundaryLayerTool
 
     BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved;
     Array<SegmentIndex> moved_segs;
-    int max_edge_nr, new_mat_nr, nfd_old;
+    int max_edge_nr, nfd_old, ndom_old;
+    Array<int> new_mat_nrs;
     int np, nseg, nse, ne;
     double height;
 
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index dd7b7910..50665ed3 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -1253,7 +1253,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array<int>{})
     .def ("BoundaryLayer", [](Mesh & self, variant<string, int> boundary,
                               variant<double, py::list> thickness,
-                              string material,
+                              variant<string, map<string, string>> material,
                               variant<string, int> domain, bool outside,
                               optional<string> project_boundaries,
                               bool grow_edges, bool limit_growth_vectors)
@@ -1298,7 +1298,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
              for(int i = 1; i<=self.GetNFD(); i++)
                if(boundaries.Test(i))
                  blp.surfid.Append(i);
-             blp.new_mat = material;
+             if(string* mat = get_if<string>(&material); mat)
+                 blp.new_mat = { { ".*", *mat } };
+             else
+                 blp.new_mat = *get_if<map<string, string>>(&material);
 
              if(project_boundaries.has_value())
                {