diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp
index 90918353..4611a95f 100644
--- a/libsrc/meshing/boundarylayer.cpp
+++ b/libsrc/meshing/boundarylayer.cpp
@@ -3,6 +3,7 @@
 
 #include <regex>
 #include <set>
+#include <variant>
 
 #include "debugging.hpp"
 #include "global.hpp"
@@ -1095,6 +1096,8 @@ void BoundaryLayerTool ::SetDomInOutSides()
       if (done.Test(index))
         continue;
       done.SetBit(index);
+      if (index < nfd_old && moved_surfaces.Test(index))
+        continue;
       auto& fd = mesh.GetFaceDescriptor(index);
       if (fd.DomainIn() != -1)
         continue;
@@ -1110,7 +1113,7 @@ void BoundaryLayerTool ::SetDomInOutSides()
       if (e2)
         dom[1] = mesh.VolumeElement(e2).GetIndex();
 
-      auto& fd_old = mesh.GetFaceDescriptor(inv_si_map[index]);
+      const auto& fd_old = mesh.GetFaceDescriptor(inv_si_map[index]);
       int dom_old[2] = {fd_old.DomainIn(), fd_old.DomainOut()};
 
       for (auto i : Range(2))
@@ -1125,7 +1128,9 @@ void BoundaryLayerTool ::SetDomInOutSides()
             }
 
           // Check if the old domain adjacent to this face gets a new boundary layer domain, if so, use that number
-          int dom_new = domains.Test(dom_old[i]) ? new_mat_nrs[dom_old[i]] : dom_old[i];
+          int dom_new = dom_old[i];
+          if (domains.Test(dom_old[i]) && new_mat_nrs[dom_old[i]] > 0)
+            dom_new = new_mat_nrs[dom_old[i]];
 
           // This case is tested by test_boundarylayer.py::test_pyramids[False] -> look at the generated mesh to understand the text below :)
           // Special case check here: when growing "outside" the new face could have the same domain on both sides (before adding blayer elements).
@@ -1249,7 +1254,10 @@ void BoundaryLayerTool ::ProcessParameters()
       if (string* mat = get_if<string>(&*params.new_material); mat)
         par_new_mat = {{".*", *mat}};
       else
-        par_new_mat = *get_if<map<string, string>>(&*params.new_material);
+        {
+          par_new_mat = *get_if<map<string, string>>(&*params.new_material);
+          have_material_map = true;
+        }
     }
 
   if (params.project_boundaries.has_value())
@@ -1416,6 +1424,12 @@ void BoundaryLayerTool ::Perform()
   mesh.SetNextMajorTimeStamp();
   mesh.UpdateTopology();
   SetDomInOutSides();
+
+  if (have_material_map)
+    {
+      AddFacesBetweenDomains(mesh);
+      mesh.SplitFacesByAdjacentDomains();
+    }
 }
 
 void GenerateBoundaryLayer (Mesh& mesh, const BoundaryLayerParameters& blp)
diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp
index cdf82aa5..d504e2d2 100644
--- a/libsrc/meshing/boundarylayer.hpp
+++ b/libsrc/meshing/boundarylayer.hpp
@@ -67,6 +67,7 @@ public:
   Array<int> par_surfid;
   bool insert_only_volume_elements;
   map<string, string> par_new_mat;
+  bool have_material_map = false;
   Array<size_t> par_project_boundaries;
 
   bool have_single_segments;
diff --git a/libsrc/meshing/boundarylayer_interpolate.cpp b/libsrc/meshing/boundarylayer_interpolate.cpp
index 9933e826..fb819cca 100644
--- a/libsrc/meshing/boundarylayer_interpolate.cpp
+++ b/libsrc/meshing/boundarylayer_interpolate.cpp
@@ -151,6 +151,11 @@ void BoundaryLayerTool ::InterpolateGrowthVectors()
               no_angles = false;
             }
         }
+
+      if (no_angles && faces.Size() == 2 && have_material_map)
+        if (par_new_mat[mesh.GetBCName(mesh[faces[0]].GetIndex() - 1)] != par_new_mat[mesh.GetBCName(mesh[faces[1]].GetIndex() - 1)])
+          no_angles = false;
+
       if (no_angles)
         {
           for (auto* p_seg : edgenr2seg[edgenr])
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 49ba8518..719f7432 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -1,7 +1,9 @@
+#include <algorithm>
 #include <mystdlib.h>
 #include <atomic>
 #include <regex>
 #include <set>
+#include "core/array.hpp"
 #include "meshing.hpp"
 #include "../general/gzstream.h"
 
@@ -7292,15 +7294,23 @@ namespace netgen
         int eli0, eli1;
         GetTopology().GetSurface2VolumeElement(sei+1, eli0, eli1);
         // auto [ei0,ei1] = GetTopology().GetSurface2VolumeElement(sei); // the way to go
+        if(eli0 == 0)
+          continue;
         auto & sel = (*this)[sei];
         int face = sel.GetIndex();
         int domin = VolumeElement(eli0).GetIndex();
         int domout = eli1 ? VolumeElement(eli1).GetIndex() : 0;
         if(domin < domout)
           swap(domin, domout);
+
         auto key = std::make_tuple(face, domin, domout);
         if(face_doms_2_new_face.find(key) == face_doms_2_new_face.end())
           {
+            {
+              auto & fd = FaceDescriptors()[face-1];
+              if(domout == 0 && min(fd.DomainIn(), fd.DomainOut()) > 0)
+                continue;
+            }
             if(!first_visit[face-1]) {
               nfaces++;
               FaceDescriptor new_fd = FaceDescriptors()[face-1];
@@ -7817,4 +7827,108 @@ namespace netgen
     return nm_;
   }
 
+  void AddFacesBetweenDomains(Mesh & mesh)
+  {
+    static Timer timer("AddFacesBetweenDomains"); RegionTimer rt(timer);
+    auto & topo = mesh.GetTopology();
+    auto p2el = mesh.CreatePoint2ElementTable();
+
+    Array<size_t> els_per_domain(mesh.GetNDomains()+1);
+    els_per_domain = 0;
+
+    for(const auto & el : mesh.VolumeElements())
+      els_per_domain[el.GetIndex()]++;
+
+    std::map<tuple<int,int>, int> doms_2_new_face;
+
+    for(const auto & [facei, fd]: Enumerate(mesh.FaceDescriptors()))
+    {
+      auto dom0 = fd.DomainIn();
+      auto dom1 = fd.DomainOut();
+      if(dom0 > dom1)
+        swap(dom0, dom1);
+
+      doms_2_new_face[{dom0, dom1}] = facei+1;
+    }
+
+    for(auto dom : Range(1, 1+mesh.GetNDomains()))
+    {
+      if(els_per_domain[dom] == 0)
+        continue;
+
+      mesh.UpdateTopology();
+
+      mesh.FindOpenElements(dom);
+      for(const auto & openel : mesh.OpenElements())
+      {
+        std::set<ElementIndex> has_p1, has_p2, has_p3;
+        for (auto ei: topo.GetVertexElements(openel[0]))
+          has_p1.insert(ei);
+        for (auto ei: topo.GetVertexElements(openel[1]))
+          has_p2.insert(ei);
+        for (auto ei: topo.GetVertexElements(openel[2]))
+          has_p3.insert(ei);
+
+        std::set<ElementIndex> has_p12, has_all;
+        set_intersection(has_p1.begin(), has_p1.end(),
+                         has_p2.begin(), has_p2.end(),
+                         inserter(has_p12, has_p12.begin()));
+        set_intersection(has_p12.begin(), has_p12.end(),
+                         has_p3.begin(), has_p3.end(),
+                         inserter(has_all, has_all.begin()));
+
+        ArrayMem<ElementIndex, 5> els;
+        for(auto ei : has_all)
+          els.Append(ei);
+
+        if(els.Size() == 2 && mesh[els[0]].GetIndex() != mesh[els[1]].GetIndex())
+        {
+          auto dom0 = mesh[els[0]].GetIndex();
+          auto dom1 = mesh[els[1]].GetIndex();
+          ElementIndex ei0 = els[0];
+          if(dom0 > dom1)
+          {
+            Swap(dom0, dom1);
+            ei0 = els[1];
+          }
+
+          if(dom1 == dom)
+            continue;
+
+          if(doms_2_new_face.count({dom0, dom1}) == 0)
+          {
+            auto fd = FaceDescriptor(-1, dom0, dom1, -1);
+            auto new_si = mesh.GetNFD()+1;
+            fd.SetBCProperty(new_si);
+            auto new_face = mesh.AddFaceDescriptor(fd);
+            mesh.SetBCName(new_si - 1, "default");
+            doms_2_new_face[{dom0, dom1}] = new_face;
+          }
+          for(auto face : topo.GetFaces(ei0)) {
+            auto verts = topo.GetFaceVertices(face);
+            if(verts.Contains(openel[0]) && verts.Contains(openel[1]) && verts.Contains(openel[2])) {
+              Element2d sel(static_cast<int>(verts.Size()));
+              sel.SetIndex(doms_2_new_face[{dom0, dom1}]);
+
+              for(auto j : Range(verts.Size()))
+                sel[j] = verts[j];
+              auto normal = Cross(mesh[sel[1]]-mesh[sel[0]], mesh[sel[2]]-mesh[sel[0]]);
+              Vec<3> surf_center = Vec<3>(Center(mesh[sel[0]] , mesh[sel[1]] , mesh[sel[2]]));
+              Vec<3> center(0., 0., 0.);
+              for(auto pi : mesh[ei0].PNums())
+                center += Vec<3>(mesh[pi]);
+              center *= 1.0/mesh[ei0].GetNP();
+              if((normal * (center - surf_center)) < 0)
+                sel.Invert();
+              mesh.AddSurfaceElement(sel);
+              break;
+            }
+          }
+        }
+      }
+    }
+
+
+  }
+
 }
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index c83ee575..ec6ba270 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -1035,6 +1035,7 @@ namespace netgen
     return FlatArray<T_FACE>(GetNFaces ( (*mesh)[elnr].GetType()), &faces[elnr][0]);
   }
 
+  DLL_HEADER void AddFacesBetweenDomains(Mesh & mesh);
   
 }
 
diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp
index 2d1f24da..ce6a4699 100644
--- a/libsrc/meshing/meshfunc.cpp
+++ b/libsrc/meshing/meshfunc.cpp
@@ -105,7 +105,7 @@ namespace netgen
 
         for( auto dom : {dom_in, dom_out} )
         {
-            if(dom==0)
+            if(dom<=0)
               continue;
 
             auto & sels = ret[dom-1].mesh->SurfaceElements();