Merge branch 'boundarylayer_fixes' into 'master'

Boundary Layers - Automatic thickness limiation and fixes

See merge request ngsolve/netgen!690

.clang-format

@@ -0,0 +1,64 @@
+Language: Cpp
+BasedOnStyle: LLVM
+AlignAfterOpenBracket: Align
+AlignConsecutiveAssignments: false
+AlignConsecutiveDeclarations: false
+AlignEscapedNewlines: Left
+AlignOperands: true
+AlignTrailingComments: true
+AllowAllParametersOfDeclarationOnNextLine: false
+AllowShortBlocksOnASingleLine: false
+AllowShortCaseLabelsOnASingleLine: true
+AllowShortFunctionsOnASingleLine: InlineOnly
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: false
+AlwaysBreakAfterDefinitionReturnType: None
+AlwaysBreakTemplateDeclarations: Yes
+BinPackArguments: false
+BinPackParameters: false
+BreakBeforeBinaryOperators: NonAssignment
+BreakBeforeBraces: Custom
+BraceWrapping:
+  AfterClass: true
+  AfterControlStatement: true
+  AfterEnum: true
+  AfterFunction: true
+  AfterNamespace: true
+  AfterStruct: true
+  AfterUnion: true
+  BeforeCatch: true
+  BeforeElse: true
+  IndentBraces: true
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+BreakInheritanceList: AfterColon
+ColumnLimit: 0
+ConstructorInitializerAllOnOneLineOrOnePerLine: true
+ConstructorInitializerIndentWidth: 2
+ContinuationIndentWidth: 2
+Cpp11BracedListStyle: true
+EmptyLineBeforeAccessModifier: Never
+IndentCaseLabels: false
+IndentPPDirectives: None
+IndentWidth: 2
+KeepEmptyLinesAtTheStartOfBlocks: false
+MaxEmptyLinesToKeep: 1
+NamespaceIndentation: None
+PointerAlignment: Left
+ReflowComments: true
+SortIncludes: false
+SpaceAfterCStyleCast: false
+SpaceAfterLogicalNot: false
+SpaceBeforeParens: Custom
+SpaceBeforeParensOptions:
+  AfterControlStatements: true
+  AfterFunctionDefinitionName: true
+  AfterFunctionDeclarationName: true
+SpacesInAngles: false
+SpacesInContainerLiterals: false
+SpacesInCStyleCastParentheses: false
+SpacesInSquareBrackets: false
+Standard: Latest
+TabWidth: 2
+UseTab: Never
+

libsrc/csg/genmesh.cpp

@@ -826,6 +826,9 @@ namespace netgen
       {
 	multithread.task = "Volume meshing";
 
+	for (int i = 0; i < geom.GetNTopLevelObjects(); i++)
+	  mesh->SetMaterial (i+1, geom.GetTopLevelObject(i)->GetMaterial().c_str());
+
 	MESHING3_RESULT res =
 	  MeshVolume (mparam, *mesh);
 
@@ -838,10 +841,6 @@ namespace netgen
 
 	MeshQuality3d (*mesh);
       
-	for (int i = 0; i < geom.GetNTopLevelObjects(); i++)
-	  mesh->SetMaterial (i+1, geom.GetTopLevelObject(i)->GetMaterial().c_str());
-      
-
 #ifdef STAT_STREAM
 	(*statout) << GetTime() << " & ";
 #endif      

libsrc/meshing/CMakeLists.txt

@@ -12,7 +12,7 @@ target_sources(nglib PRIVATE
         parallelmesh.cpp  paralleltop.cpp  basegeom.cpp
         python_mesh.cpp surfacegeom.cpp
         debugging.cpp fieldlines.cpp visual_interface.cpp
-        boundarylayer2d.cpp
+        boundarylayer2d.cpp boundarylayer_interpolate.cpp
 )
 
 target_link_libraries( nglib PRIVATE $<BUILD_INTERFACE:netgen_metis> $<BUILD_INTERFACE:netgen_python> )

libsrc/meshing/basegeom.hpp

@@ -271,7 +271,7 @@ namespace netgen
 
   virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, EdgePointGeomInfo* gi = nullptr) const
   {
-    if(gi && gi->edgenr < edges.Size())
+    if(gi && gi->edgenr < edges.Size() && gi->edgenr >= 0)
       edges[gi->edgenr]->ProjectPoint(p, gi);
   }
 

libsrc/meshing/boundarylayer.hpp

@@ -1,85 +1,110 @@
 #ifndef NETGEN_BOUNDARYLAYER_HPP
 #define NETGEN_BOUNDARYLAYER_HPP
 
+#include <core/array.hpp>
+#include <mystdlib.h>
+#include <meshing.hpp>
+
 namespace netgen
 {
 
 ///
-DLL_HEADER extern void InsertVirtualBoundaryLayer (Mesh & mesh);
+DLL_HEADER extern void InsertVirtualBoundaryLayer (Mesh& mesh);
 
-/// Create a typical prismatic boundary layer on the given 
+/// Create a typical prismatic boundary layer on the given
 /// surfaces
 
-DLL_HEADER void GenerateBoundaryLayer (Mesh & mesh,
+struct SpecialBoundaryPoint
-                                       const BoundaryLayerParameters & blp);
+{
+  struct GrowthGroup
+  {
+    Array<int> faces;
+    Vec<3> growth_vector;
+    Array<PointIndex> new_points;
 
-DLL_HEADER int /* new_domain_number */ GenerateBoundaryLayer2 (Mesh & mesh, int domain, const Array<double> & thicknesses, bool should_make_new_domain=true, const Array<int> & boundaries=Array<int>{});
+    GrowthGroup(FlatArray<int> faces_, FlatArray<Vec<3>> normals);
+    GrowthGroup(const GrowthGroup&) = default;
+    GrowthGroup() = default;
+  };
+  Array<GrowthGroup> growth_groups;
+  Vec<3> separating_direction;
+
+  SpecialBoundaryPoint(const std::map<int, Vec<3>>& normals);
+  SpecialBoundaryPoint() = default;
+};
+
+DLL_HEADER void GenerateBoundaryLayer (Mesh& mesh,
+                                       const BoundaryLayerParameters& blp);
+
+DLL_HEADER int /* new_domain_number */ GenerateBoundaryLayer2 (Mesh& mesh, int domain, const Array<double>& thicknesses, bool should_make_new_domain = true, const Array<int>& boundaries = Array<int>{});
 
 class BoundaryLayerTool
 {
-  public:
+public:
-    BoundaryLayerTool(Mesh & mesh_, const BoundaryLayerParameters & params_);
+  BoundaryLayerTool(Mesh& mesh_, const BoundaryLayerParameters& params_);
-    void ProcessParameters();
+  void ProcessParameters ();
-    void Perform();
+  void Perform ();
 
-  protected:
+  Mesh& mesh;
-    Mesh & mesh;
+  MeshTopology& topo;
-    MeshTopology & topo;
+  BoundaryLayerParameters params;
-    BoundaryLayerParameters params;
+  Array<Vec<3>, PointIndex> growthvectors;
-    Array<Vec<3>, PointIndex> growthvectors;
+  std::map<PointIndex, Vec<3>> non_bl_growth_vectors;
-    Table<SurfaceElementIndex, PointIndex> p2sel;
+  Table<SurfaceElementIndex, PointIndex> p2sel;
 
-    BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved;
+  BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved;
-    Array<SegmentIndex> moved_segs;
+  Array<SegmentIndex> moved_segs;
-    int max_edge_nr, nfd_old, ndom_old;
+  int max_edge_nr, nfd_old, ndom_old;
-    Array<int> new_mat_nrs;
+  Array<int> new_mat_nrs;
-    BitArray moved_surfaces;
+  BitArray moved_surfaces;
-    int np, nseg, nse, ne;
+  int np, nseg, nse, ne;
-    double height;
+  double total_height;
+  Array<POINTTYPE, PointIndex> point_types;
 
-    // These parameters are derived from given BoundaryLayerParameters and the Mesh
+  // These parameters are derived from given BoundaryLayerParameters and the Mesh
-    Array<double> par_heights;
+  Array<double> par_heights;
-    Array<int> par_surfid;
+  Array<int> par_surfid;
-    map<string, string> par_new_mat;
+  bool insert_only_volume_elements;
-    Array<size_t> par_project_boundaries;
+  map<string, string> par_new_mat;
+  Array<size_t> par_project_boundaries;
 
-    bool have_single_segments;
+  bool have_single_segments;
-    Array<Segment> segments, new_segments;
+  Array<Segment> old_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;
 
-    Array<double> surfacefacs;
+  Array<double> surfacefacs;
-    Array<int> si_map;
+  Array<int> si_map;
-    Array<double, PointIndex> limits;
 
-    // major steps called in Perform()
+  std::map<PointIndex, SpecialBoundaryPoint> special_boundary_points;
-    void CreateNewFaceDescriptors();
+  std::map<PointIndex, std::tuple<Vec<3>*, double>> growth_vector_map;
-    void CreateFaceDescriptorsSides();
-    void CalculateGrowthVectors();
-    Array<Array<pair<SegmentIndex, int>>, SegmentIndex> BuildSegMap();
 
-    BitArray ProjectGrowthVectorsOnSurface();
+  // major steps called in Perform()
-    void InterpolateSurfaceGrowthVectors();
+  void CreateNewFaceDescriptors ();
-    void InterpolateGrowthVectors();
+  void CreateFaceDescriptorsSides ();
-    void LimitGrowthVectorLengths();
+  void CalculateGrowthVectors ();
+  Array<Array<pair<SegmentIndex, int>>, SegmentIndex> BuildSegMap ();
 
-    void InsertNewElements(FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap, const BitArray & in_surface_direction);
+  BitArray ProjectGrowthVectorsOnSurface ();
-    void SetDomInOut();
+  void InterpolateSurfaceGrowthVectors ();
-    void SetDomInOutSides();
+  void InterpolateGrowthVectors ();
-    void AddSegments();
+  void LimitGrowthVectorLengths ();
-    void FixVolumeElements();
+  void FixSurfaceElements ();
 
-    // utility functions
+  void InsertNewElements (FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap, const BitArray& in_surface_direction);
-    array<Point<3>, 2> GetMappedSeg( PointIndex pi );
+  void SetDomInOut ();
-    ArrayMem<Point<3>, 4> GetFace( SurfaceElementIndex sei );
+  void SetDomInOutSides ();
-    ArrayMem<Point<3>, 4> GetMappedFace( SurfaceElementIndex sei );
+  void AddSegments ();
-    ArrayMem<Point<3>, 4> GetMappedFace( SurfaceElementIndex sei, int face );
+  void AddSurfaceElements ();
 
-    Vec<3> getNormal(const Element2d & el)
+  Vec<3> getNormal (const Element2d& el)
-    {
+  {
-        auto v0 = mesh[el[0]];
+    auto v0 = mesh[el[0]];
-        return Cross(mesh[el[1]]-v0, mesh[el[2]]-v0).Normalize();
+    return Cross(mesh[el[1]] - v0, mesh[el[2]] - v0).Normalize();
-    }
+  }
 
-    Vec<3> getEdgeTangent(PointIndex pi, int edgenr);
+  Vec<3> getEdgeTangent (PointIndex pi, int edgenr, FlatArray<Segment*> segs);
 };
 
 } // namespace netgen

libsrc/meshing/boundarylayer2d.cpp

@@ -1,5 +1,5 @@
 #include <mystdlib.h>
-#include "meshing.hpp"
+#include "boundarylayer.hpp"
 #include "meshing2.hpp"
 #include "../geom2d/csg2d.hpp"
 

libsrc/meshing/boundarylayer_interpolate.cpp

@@ -0,0 +1,472 @@
+#include "boundarylayer.hpp"
+
+namespace netgen
+{
+
+namespace detail
+{
+struct Neighbor
+{
+  PointIndex pi;
+  SurfaceElementIndex sei;
+  double weight;
+};
+} // namespace detail
+
+Array<ArrayMem<detail::Neighbor, 20>>
+BuildNeighbors (FlatArray<PointIndex> points, const Mesh& mesh)
+{
+  auto p2sel = mesh.CreatePoint2SurfaceElementTable();
+
+  Array<ArrayMem<detail::Neighbor, 20>> neighbors(points.Size());
+
+  ArrayMem<double, 20> angles;
+  ArrayMem<double, 20> inv_dists;
+  for (auto i : points.Range())
+    {
+      auto& p_neighbors = neighbors[i];
+      auto pi = points[i];
+      angles.SetSize(0);
+      inv_dists.SetSize(0);
+      for (auto sei : p2sel[pi])
+        {
+          const auto& sel = mesh[sei];
+          for (auto pi1 : sel.PNums())
+            {
+              if (pi1 == pi)
+                continue;
+              auto pi2 = pi1;
+              for (auto pi_ : sel.PNums())
+                {
+                  if (pi_ != pi && pi_ != pi1)
+                    {
+                      pi2 = pi_;
+                      break;
+                    }
+                }
+              p_neighbors.Append({pi1, sei, 0.0});
+              inv_dists.Append(1.0 / (mesh[pi1] - mesh[pi]).Length());
+              auto dot = (mesh[pi1] - mesh[pi]).Normalize() * (mesh[pi2] - mesh[pi]).Normalize();
+              angles.Append(acos(dot));
+            }
+        }
+      double sum_inv_dist = 0.0;
+      for (auto inv_dist : inv_dists)
+        sum_inv_dist += inv_dist;
+      double sum_angle = 0.0;
+      for (auto angle : angles)
+        sum_angle += angle;
+
+      double sum_weight = 0.0;
+      for (auto i : Range(inv_dists))
+        {
+          p_neighbors[i].weight =
+            inv_dists[i] * angles[i] / sum_inv_dist / sum_angle;
+          sum_weight += p_neighbors[i].weight;
+        }
+      for (auto i : Range(inv_dists))
+        p_neighbors[i].weight /= sum_weight;
+    }
+  return neighbors;
+}
+
+void BoundaryLayerTool ::InterpolateGrowthVectors()
+{
+  point_types.SetSize(mesh.GetNP());
+  for (auto p : mesh.Points().Range())
+    point_types[p] = mesh[p].Type();
+
+  int new_max_edge_nr = max_edge_nr;
+  for (const auto& seg : segments)
+    if (seg.edgenr > new_max_edge_nr)
+      new_max_edge_nr = seg.edgenr;
+  for (const auto& seg : new_segments)
+    if (seg.edgenr > new_max_edge_nr)
+      new_max_edge_nr = seg.edgenr;
+
+  auto getGW = [&] (PointIndex pi) -> Vec<3> {
+    if (growth_vector_map.count(pi) == 0)
+      growth_vector_map[pi] = {&growthvectors[pi], total_height};
+    auto [gw, height] = growth_vector_map[pi];
+    return height * (*gw);
+  };
+  auto addGW = [&] (PointIndex pi, Vec<3> vec) {
+    if (growth_vector_map.count(pi) == 0)
+      growth_vector_map[pi] = {&growthvectors[pi], total_height};
+    auto [gw, height] = growth_vector_map[pi];
+    *gw += 1.0 / height * vec;
+  };
+
+  // interpolate tangential component of growth vector along edge
+  if (max_edge_nr >= new_max_edge_nr)
+    return;
+
+  auto edgenr2seg = ngcore::CreateSortedTable<Segment*, int>(
+    Range(segments.Size() + new_segments.Size()),
+    [&] (auto& table, size_t segi) {
+      auto& seg = segi < segments.Size()
+                    ? segments[segi]
+                    : new_segments[segi - segments.Size()];
+      table.Add(seg.edgenr, &seg);
+    },
+    new_max_edge_nr + 1);
+  auto point2seg = ngcore::CreateSortedTable<Segment*, PointIndex>(
+    Range(segments.Size() + new_segments.Size()),
+    [&] (auto& table, size_t segi) {
+      auto& seg = segi < segments.Size()
+                    ? segments[segi]
+                    : new_segments[segi - segments.Size()];
+      table.Add(seg[0], &seg);
+      table.Add(seg[1], &seg);
+    },
+    mesh.GetNP());
+
+  for (auto edgenr : Range(1, new_max_edge_nr + 1))
+    {
+      // "inner" edges between two flat faces are not treated as edges for interpolation
+      bool no_angles = true;
+      ArrayMem<SurfaceElementIndex, 4> faces;
+
+      for (auto* p_seg : edgenr2seg[edgenr])
+        {
+          auto& seg = *p_seg;
+          faces.SetSize(0);
+          if (seg[0] <= p2sel.Size())
+            {
+              for (auto sei : p2sel[seg[0]])
+                if (moved_surfaces.Test(mesh[sei].GetIndex()) && p2sel[seg[1]].Contains(sei))
+                  faces.Append(sei);
+            }
+
+          if (faces.Size() == 2)
+            {
+              auto n0 = getNormal(mesh[faces[0]]);
+              auto n1 = getNormal(mesh[faces[1]]);
+              if (n0 * n1 < 0.999)
+                no_angles = false;
+            }
+          else
+            {
+              no_angles = false;
+            }
+        }
+      if (no_angles)
+        {
+          for (auto* p_seg : edgenr2seg[edgenr])
+            for (auto pi : p_seg->PNums())
+              if (pi <= np && point_types[pi] == EDGEPOINT)
+                point_types[pi] = SURFACEPOINT;
+          continue;
+        }
+    }
+
+  for (auto edgenr : Range(max_edge_nr + 1, new_max_edge_nr + 1))
+    {
+      double edge_len = 0.;
+      bool any_grows = false;
+
+      auto is_end_point = [&] (PointIndex pi) {
+        auto segs = point2seg[pi];
+        if (segs.Size() == 1)
+          return true;
+        auto first_edgenr = (*segs[0]).edgenr;
+        for (auto* p_seg : segs)
+          if (p_seg->edgenr != first_edgenr)
+            return true;
+        return false;
+      };
+
+      Array<PointIndex> points;
+      for (auto* p_seg : edgenr2seg[edgenr])
+        {
+          auto& seg = *p_seg;
+
+          if (getGW(seg[0]).Length2() != 0 || getGW(seg[1]).Length2() != 0)
+            any_grows = true;
+
+          if (points.Size() == 0)
+            for (auto i : Range(2))
+              if (is_end_point(seg[i]))
+                {
+                  points.Append(seg[i]);
+                  points.Append(seg[1 - i]);
+                  edge_len += (mesh[seg[1]] - mesh[seg[0]]).Length();
+                  break;
+                }
+        }
+
+      if (!any_grows)
+        {
+          PrintMessage(1, "BLayer: skip interpolating growth vectors at edge ", edgenr + 1);
+          continue;
+        }
+
+      if (!points.Size())
+        {
+          if (debugparam.debugoutput)
+            cerr << "Could not find startpoint for edge " << edgenr << endl;
+          continue;
+        }
+
+      std::set<PointIndex> points_set;
+      points_set.insert(points[0]);
+      points_set.insert(points[1]);
+
+      bool point_found = true;
+      while (point_found)
+        {
+          if (is_end_point(points.Last()))
+            break;
+          point_found = false;
+          for (auto* p_seg : point2seg[points.Last()])
+            {
+              const auto& seg = *p_seg;
+              if (seg.edgenr != edgenr)
+                continue;
+              auto plast = points.Last();
+              if (plast != seg[0] && plast != seg[1])
+                continue;
+              auto pnew = plast == seg[0] ? seg[1] : seg[0];
+              if (pnew == points[0] && points.Size() > 1)
+                {
+                }
+              if (points_set.count(pnew) > 0 && (pnew != points[0] || points.Size() == 2))
+                continue;
+              edge_len += (mesh[points.Last()] - mesh[pnew]).Length();
+              points.Append(pnew);
+              points_set.insert(pnew);
+              point_found = true;
+              break;
+            }
+        }
+      if (!point_found)
+        {
+          if (debugparam.debugoutput)
+            {
+              cerr << "Could not find connected list of line segments for edge "
+                   << edgenr << endl;
+              cerr << "current points: " << endl
+                   << points << endl;
+            }
+          continue;
+        }
+
+      if (getGW(points[0]).Length2() == 0 && getGW(points.Last()).Length2() == 0)
+        continue;
+
+      // tangential part of growth vectors
+      auto t1 = (mesh[points[1]] - mesh[points[0]]).Normalize();
+      auto gt1 = getGW(points[0]) * t1 * t1;
+      auto t2 =
+        (mesh[points.Last()] - mesh[points[points.Size() - 2]]).Normalize();
+      auto gt2 = getGW(points.Last()) * t2 * t2;
+
+      double len = 0.;
+      for (auto i : IntRange(1, points.Size() - 1))
+        {
+          auto pi = points[i];
+          len += (mesh[pi] - mesh[points[i - 1]]).Length();
+          auto t = getEdgeTangent(pi, edgenr, point2seg[pi]);
+          auto lam = len / edge_len;
+          auto interpol = (1 - lam) * (gt1 * t) * t + lam * (gt2 * t) * t;
+          addGW(pi, interpol);
+        }
+    }
+}
+
+void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors()
+{
+  static Timer tall("InterpolateSurfaceGrowthVectors");
+  RegionTimer rtall(tall);
+  static Timer tsmooth("InterpolateSurfaceGrowthVectors-Smoothing");
+  auto np_old = this->np;
+  auto np = mesh.GetNP();
+
+  auto hasMoved = [&] (PointIndex pi) {
+    return (pi - PointIndex::BASE >= np_old) || mapto[pi].Size() > 0 || special_boundary_points.count(pi);
+  };
+
+  std::set<PointIndex> points_set;
+  for (const auto& sel : mesh.SurfaceElements())
+    {
+      for (auto pi : sel.PNums())
+        if (point_types[pi] == SURFACEPOINT && hasMoved(pi))
+          points_set.insert(pi);
+    }
+
+  Array<PointIndex> points;
+  for (auto pi : points_set)
+    points.Append(pi);
+  QuickSort(points);
+
+  // smooth tangential part of growth vectors from edges to surface elements
+  Array<Vec<3>, PointIndex> corrections(mesh.GetNP());
+  corrections = 0.0;
+  RegionTimer rtsmooth(tsmooth);
+  auto neighbors = BuildNeighbors(points, mesh);
+
+  Array<Vec<3>, SurfaceElementIndex> surf_normals(mesh.GetNSE());
+  for (auto sei : mesh.SurfaceElements().Range())
+    surf_normals[sei] = getNormal(mesh[sei]);
+
+  BitArray interpolate_tangent(mesh.GetNP() + 1);
+  interpolate_tangent = false;
+  for (auto pi : points)
+    {
+      for (auto sei : p2sel[pi])
+        if (is_boundary_moved[mesh[sei].GetIndex()])
+          interpolate_tangent.SetBit(pi);
+    }
+
+  constexpr int N_STEPS = 64;
+  for ([[maybe_unused]] auto i : Range(N_STEPS))
+    {
+      for (auto i : points.Range())
+        {
+          auto pi = points[i];
+          auto& p_neighbors = neighbors[i];
+
+          ArrayMem<Vec<3>, 20> g_vectors;
+          double max_len = 0.0;
+          double sum_len = 0.0;
+
+          // average only tangent component on new bl points, average whole growth
+          // vector otherwise
+          bool do_average_tangent = true;
+          for (const auto& s : p_neighbors)
+            {
+              auto gw_other = growthvectors[s.pi] + corrections[s.pi];
+              if (do_average_tangent)
+                {
+                  auto n = surf_normals[s.sei];
+                  gw_other = gw_other - (gw_other * n) * n;
+                }
+              auto v = gw_other;
+              auto len = v.Length2();
+              sum_len += len;
+              max_len = max(max_len, len);
+              g_vectors.Append(v);
+            }
+
+          if (max_len == 0.0)
+            continue;
+
+          double lambda = 0;
+          if (i > N_STEPS / 4.)
+            lambda = 2.0 * (i - N_STEPS / 4.) / (N_STEPS / 2.);
+          lambda = min(1.0, lambda);
+
+          auto& correction = corrections[pi];
+          correction = 0.0;
+          for (const auto i : p_neighbors.Range())
+            {
+              auto v = g_vectors[i];
+              double weight = lambda * p_neighbors[i].weight + (1.0 - lambda) * v.Length2() / sum_len;
+              correction += weight * v;
+            }
+
+          if (!do_average_tangent)
+            correction -= growthvectors[pi];
+        }
+    }
+
+  for (auto pi : points)
+    growthvectors[pi] += corrections[pi];
+}
+
+void BoundaryLayerTool ::FixSurfaceElements()
+{
+  static Timer tall("FixSurfaceElements");
+  RegionTimer rtall(tall);
+  auto np_old = this->np;
+  auto np = mesh.GetNP();
+
+  non_bl_growth_vectors.clear();
+
+  auto getGW = [&] (PointIndex pi) -> Vec<3> {
+    // return growthvectors[pi];
+    if (growth_vector_map.count(pi) == 0)
+      {
+        non_bl_growth_vectors[pi] = .0;
+        growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0};
+      }
+    auto [gw, height] = growth_vector_map[pi];
+    return height * (*gw);
+  };
+
+  auto addGW = [&] (PointIndex pi, Vec<3> vec) {
+    if (growth_vector_map.count(pi) == 0)
+      {
+        non_bl_growth_vectors[pi] = .0;
+        growth_vector_map[pi] = {&non_bl_growth_vectors[pi], 1.0};
+      }
+    auto [gw, height] = growth_vector_map[pi];
+    *gw += 1.0 / height * vec;
+  };
+
+  std::set<PointIndex> points_set;
+  // only smooth over old surface elements
+  for (SurfaceElementIndex sei : Range(nse))
+    {
+      const auto& sel = mesh[sei];
+      if (sel.GetNP() == 3 && is_boundary_moved[sel.GetIndex()])
+        for (auto pi : sel.PNums())
+          if (point_types[pi] == SURFACEPOINT)
+            points_set.insert(pi);
+    }
+
+  Array<PointIndex> points;
+  for (auto pi : points_set)
+    points.Append(pi);
+  QuickSort(points);
+
+  Array<Vec<3>, PointIndex> corrections(mesh.GetNP());
+  corrections = 0.0;
+
+  auto neighbors = BuildNeighbors(points, mesh);
+
+  constexpr int N_STEPS = 32;
+  for ([[maybe_unused]] auto i : Range(N_STEPS))
+    {
+      for (auto i : points.Range())
+        {
+          auto pi = points[i];
+          auto& p_neighbors = neighbors[i];
+
+          ArrayMem<Vec<3>, 20> g_vectors;
+          double max_len = 0.0;
+          double sum_len = 0.0;
+
+          for (const auto& s : p_neighbors)
+            {
+              auto v = getGW(s.pi) + corrections[s.pi];
+              auto len = v.Length2();
+              sum_len += len;
+              max_len = max(max_len, len);
+              g_vectors.Append(v);
+            }
+
+          if (max_len == 0.0)
+            continue;
+
+          double lambda = 0;
+          if (i > N_STEPS / 4.)
+            lambda = 2.0 * (i - N_STEPS / 4.) / (N_STEPS / 2.);
+          lambda = min(1.0, lambda);
+
+          auto& correction = corrections[pi];
+          correction = 0.0;
+          for (const auto i : p_neighbors.Range())
+            {
+              auto v = g_vectors[i];
+              double weight = lambda * p_neighbors[i].weight + (1.0 - lambda) * v.Length2() / sum_len;
+              correction += weight * v;
+            }
+        }
+    }
+
+  for (auto pi : points)
+    addGW(pi, corrections[pi]);
+}
+
+} // namespace netgen

libsrc/meshing/boundarylayer_limiter.hpp

@@ -0,0 +1,752 @@
+#include "boundarylayer.hpp"
+#include <core/array.hpp>
+
+namespace netgen
+{
+
+struct Intersection_
+{
+  bool is_intersecting = false;
+  double lam0 = -1, lam1 = -1;
+  Point<3> p;
+  double bary[3];
+  operator bool() const { return is_intersecting; }
+};
+
+struct GrowthVectorLimiter
+{
+  typedef std::array<Point<3>, 2> Seg;
+  typedef std::array<Point<3>, 3> Trig;
+
+  BoundaryLayerTool& tool;
+  const BoundaryLayerParameters& params;
+  Mesh& mesh;
+  double height;
+  Array<double, PointIndex> limits;
+  FlatArray<Vec<3>, PointIndex> growthvectors;
+  BitArray changed_domains;
+  unique_ptr<BoxTree<3>> tree;
+  Array<PointIndex, PointIndex> map_from;
+  Table<SurfaceElementIndex, PointIndex> p2sel;
+
+  GrowthVectorLimiter(BoundaryLayerTool& tool_)
+    : tool(tool_), params(tool_.params), mesh(tool_.mesh), height(tool_.total_height), growthvectors(tool_.growthvectors), map_from(mesh.Points().Size())
+  {
+    changed_domains = tool.domains;
+    if (!params.outside)
+      changed_domains.Invert();
+
+    map_from = tool.mapfrom;
+    p2sel = ngcore::CreateSortedTable<SurfaceElementIndex, PointIndex>(
+      tool.new_sels.Range(),
+      [&] (auto& table, SurfaceElementIndex ei) {
+        for (PointIndex pi : tool.new_sels[ei].PNums())
+          table.Add(pi, ei);
+      },
+      mesh.GetNP());
+  }
+
+  auto SurfaceElementsRange () { return Range(tool.nse + tool.new_sels.Size()); }
+
+  void WriteErrorMesh (string name)
+  {
+    if (!debugparam.write_mesh_on_error)
+      return;
+    Mesh out_mesh;
+    out_mesh = mesh;
+    for (auto [pi, data] : tool.growth_vector_map)
+      {
+        auto [gw, height] = data;
+        out_mesh[pi] += limits[pi] * height * (*gw);
+      }
+    out_mesh.Save(name);
+  }
+
+  const auto& Get (SurfaceElementIndex sei)
+  {
+    if (sei < tool.nse)
+      return mesh[sei];
+    return tool.new_sels[sei - tool.nse];
+  }
+
+  std::pair<double, double> GetMinMaxLimit (SurfaceElementIndex sei)
+  {
+    const auto& sel = Get(sei);
+    double min_limit = GetLimit(sel[0]);
+    double max_limit = min_limit;
+    for (auto i : IntRange(1, sel.GetNP()))
+      {
+        auto limit = GetLimit(sel[i]);
+        min_limit = min(min_limit, limit);
+        max_limit = max(max_limit, limit);
+      }
+    return {min_limit, max_limit};
+  }
+
+  double GetLimit (PointIndex pi)
+  {
+    if (pi <= tool.np)
+      return limits[pi];
+    return limits[map_from[pi]];
+  }
+
+  bool SetLimit (PointIndex pi, double new_limit)
+  {
+    double& limit = (pi <= tool.np) ? limits[pi] : limits[map_from[pi]];
+    if (limit <= new_limit)
+      return false;
+    limit = new_limit;
+    return true;
+  }
+
+  bool ScaleLimit (PointIndex pi, double factor)
+  {
+    double& limit = (pi <= tool.np) ? limits[pi] : limits[map_from[pi]];
+    return SetLimit(pi, limit * factor);
+  }
+
+  Vec<3> GetVector (PointIndex pi_to, double shift = 1., bool apply_limit = false)
+  {
+    auto [gw, height] = tool.growth_vector_map[pi_to];
+    if (apply_limit)
+      shift *= GetLimit(pi_to);
+    return shift * height * (*gw);
+  }
+
+  Point<3> GetPoint (PointIndex pi_to, double shift = 1., bool apply_limit = false)
+  {
+    if (pi_to <= tool.np || tool.growth_vector_map.count(pi_to) == 0)
+      return mesh[pi_to];
+
+    return mesh[pi_to] + GetVector(pi_to, shift, apply_limit);
+  }
+
+  Point<3> GetMappedPoint (PointIndex pi_from, double shift = 1., bool apply_limit = false)
+  {
+    auto pi_to = tool.mapto[pi_from].Last();
+    return GetPoint(pi_to, shift, apply_limit);
+  }
+
+  Seg GetMappedSeg (PointIndex pi_from, double shift = 1.)
+  {
+    return {mesh[pi_from], GetMappedPoint(pi_from, shift)};
+  }
+
+  Seg GetSeg (PointIndex pi_to, double shift = 1., bool apply_limit = false)
+  {
+    return {GetPoint(pi_to, 0), GetPoint(pi_to, shift, apply_limit)};
+  }
+
+  Trig GetTrig (SurfaceElementIndex sei, double shift = 0.0, bool apply_limit = false)
+  {
+    auto sel = Get(sei);
+    Trig trig;
+    for (auto i : Range(3))
+      trig[i] = GetPoint(sel[i], shift, apply_limit);
+    return trig;
+  }
+
+  Trig GetMappedTrig (SurfaceElementIndex sei, double shift = 0.0)
+  {
+    auto sel = Get(sei);
+    Trig trig;
+    for (auto i : Range(3))
+      trig[i] = GetMappedPoint(sel[i], shift);
+    return trig;
+  }
+
+  Trig GetSideTrig (SurfaceElementIndex sei, int index, double shift = 0.0, bool grow_first_vertex = true)
+  {
+    auto trig = GetMappedTrig(sei, 0.0);
+    auto sel = Get(sei);
+    auto index1 = (index + 1) % 3;
+    if (!grow_first_vertex)
+      index1 = (index + 2) % 3;
+    trig[index] = GetMappedPoint(sel[index1], shift, true);
+    return trig;
+  }
+
+  array<Trig, 4> GetSideTrigs (SurfaceElementIndex sei, int i0, double shift = 0.0)
+  {
+    auto trig = GetMappedTrig(sei, 0.0);
+    array<Trig, 4> trigs{trig, trig, trig, trig};
+
+    auto sel = Get(sei);
+    auto i1 = (i0 + 1) % 3;
+    auto i2 = (i0 + 2) % 3;
+    auto p1 = GetMappedPoint(sel[i1], shift, true);
+    auto p2 = GetMappedPoint(sel[i2], shift, true);
+
+    // create four trigs to span the quad from i1,i2 and their shifted points
+    // i1, i2, shifted i1
+    trigs[0][i0] = p1;
+
+    // i1, i2, shifted i2
+    trigs[1][i0] = p2;
+
+    // i1, shifted i1, shifted i2
+    trigs[2][i0] = p1;
+    trigs[2][i2] = p2;
+
+    // i2, shifted i1, shifted i2
+    trigs[2][i0] = p2;
+    trigs[2][i1] = p1;
+
+    return trigs;
+  }
+
+  static constexpr double INTERSECTION_SAFETY = .9;
+  bool LimitGrowthVector (PointIndex pi_to, SurfaceElementIndex sei, double trig_shift, double seg_shift, bool check_prism_sides = false)
+  {
+    auto pi_from = map_from[pi_to];
+    if (!pi_from.IsValid())
+      return false;
+
+    for (auto pi : Get(sei).PNums())
+      {
+        if (pi == pi_from)
+          return false;
+        if (map_from[pi] == pi_from)
+          return false;
+      }
+
+    if (check_prism_sides || trig_shift > .0)
+      {
+        auto [trig_min_limit, trig_max_limit] = GetMinMaxLimit(sei);
+        if (GetLimit(pi_to) < trig_min_limit)
+          return false;
+
+        auto getTrigs = [&] (double scaling = 1.0) -> ArrayMem<Trig, 3> {
+          ArrayMem<Trig, 12> trigs;
+          if (check_prism_sides)
+            for (auto i : Range(3))
+              for (auto trig : GetSideTrigs(sei, i, scaling * trig_shift))
+                trigs.Append(trig);
+          else
+            trigs.Append(GetTrig(sei, scaling * trig_shift, true));
+          return trigs;
+        };
+
+        if (!check_prism_sides)
+          {
+            // If the growth vectors of all points are pointing in the same direction,
+            // an intersection means, we also have an intersection with a prism side face
+            // this is an extra check and handled later
+            auto seg = GetSeg(pi_to, 1.0, false);
+            auto gw = seg[1] - seg[0];
+
+            bool have_same_growth_direction = true;
+            for (auto pi : Get(sei).PNums())
+              {
+                auto p_seg = GetSeg(pi, 1.0, false);
+                auto p_gw = p_seg[1] - p_seg[0];
+                have_same_growth_direction &= (gw * p_gw) > 0;
+              }
+            if (have_same_growth_direction)
+              return false;
+          }
+
+        double scaling = 1.0;
+        while (true)
+          {
+            bool have_intersection = false;
+            auto seg = GetSeg(pi_to, scaling * seg_shift, true);
+            for (auto trig : getTrigs(scaling))
+              have_intersection |= isIntersectingTrig(seg, trig);
+            if (!have_intersection)
+              break;
+            scaling *= 0.9;
+          }
+        if (scaling == 1.0)
+          return false;
+
+        double new_limit = scaling * max(GetLimit(pi_to), trig_max_limit);
+        SetLimit(pi_to, new_limit);
+        for (auto pi : Get(sei).PNums())
+          SetLimit(pi, new_limit);
+        return true;
+      }
+    else
+      {
+        auto seg = GetSeg(pi_to, seg_shift, false);
+        auto trig = GetTrig(sei, 0.0);
+        auto intersection = isIntersectingTrig(seg, trig);
+        // checking with original surface elements -> allow only half the distance
+        auto new_seg_limit = 0.40 * intersection.lam0 * seg_shift;
+        if (intersection && new_seg_limit < GetLimit(pi_from))
+          return SetLimit(pi_from, new_seg_limit);
+        return false;
+      }
+  }
+
+  void EqualizeLimits (double factor = .5)
+  {
+    static Timer t("GrowthVectorLimiter::EqualizeLimits");
+    PrintMessage(5, "GrowthVectorLimiter - equalize limits");
+    RegionTimer reg(t);
+    if (factor == 0.0)
+      return;
+    for (PointIndex pi : IntRange(tool.np, mesh.GetNP()))
+      {
+        auto pi_from = map_from[pi];
+        std::set<PointIndex> pis;
+        for (auto sei : p2sel[pi])
+          for (auto pi_ : tool.new_sels[sei].PNums())
+            pis.insert(pi_);
+        ArrayMem<double, 20> limits;
+        for (auto pi1 : pis)
+          {
+            auto limit = GetLimit(pi1);
+            if (limit > 0.0)
+              limits.Append(GetLimit(pi1));
+          }
+
+        if (limits.Size() == 0)
+          continue;
+
+        double average = 0.0;
+        for (auto l : limits)
+          average += l;
+        average /= limits.Size();
+
+        SetLimit(pi, factor * average + (1.0 - factor) * GetLimit(pi));
+      }
+  }
+
+  void LimitSelfIntersection (double safety = 1.4)
+  {
+    static Timer t("GrowthVectorLimiter::LimitSelfIntersection");
+    PrintMessage(5, "GrowthVectorLimiter - self intersection");
+    RegionTimer reg(t);
+    // check for self-intersection within new elements (prisms/hexes)
+    auto isIntersecting = [&] (SurfaceElementIndex sei, double shift) {
+      // checks if surface element is self intersecting when growing with factor
+      // shift
+
+      // ignore new surface elements, side trigs are only built
+      // from original surface elements
+      if (sei >= tool.nse)
+        return false;
+      const auto sel = Get(sei);
+      auto np = sel.GetNP();
+      for (auto i : Range(np))
+        {
+          if (sel[i] > tool.np)
+            return false;
+          if (tool.mapto[sel[i]].Size() == 0)
+            return false;
+        }
+      for (auto i : Range(np))
+        {
+          auto seg = GetMappedSeg(sel[i], shift * limits[sel[i]]);
+          for (auto fi : Range(np - 2))
+            {
+              for (auto side : {true, false})
+                {
+                  auto trig = GetSideTrig(sei, i + fi, 1.0, side);
+                  if (isIntersectingPlane(seg, trig))
+                    return true;
+                }
+            }
+        }
+      return false;
+    };
+
+    for (SurfaceElementIndex sei : mesh.SurfaceElements().Range())
+      {
+        auto sel = mesh[sei];
+        if (!tool.moved_surfaces[sel.GetIndex()])
+          continue;
+        if (sel.GetNP() == 4)
+          continue;
+
+        const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
+        auto np = sel.GetNP();
+
+        double shift = 1.0;
+        const double step_factor = 0.9;
+        while (isIntersecting(sei, shift * safety))
+          {
+            shift *= step_factor;
+            double max_limit = 0;
+            for (auto i : Range(np))
+              max_limit = max(max_limit, limits[sel[i]]);
+            for (auto i : Range(np))
+              if (max_limit == limits[sel[i]])
+                ScaleLimit(sel[i], step_factor);
+            // if (max_limit < 0.01) break;
+          }
+      }
+  }
+
+  // checks if a segment is intersecting a plane, spanned by three points, lam
+  // will be set s.t. p_intersect = seg[0] + lam * (seg[1]-seg[0])
+  Intersection_ isIntersectingPlane (const Seg& seg,
+                                     const Trig& trig)
+  {
+    auto t1 = trig[1] - trig[0];
+    auto t2 = trig[2] - trig[0];
+    auto n = Cross(t1, t2);
+    auto v0n = (seg[0] - trig[0]) * n;
+    auto v1n = (seg[1] - trig[0]) * n;
+
+    Intersection_ intersection;
+    intersection.lam0 = -v0n / (v1n - v0n);
+    intersection.p = seg[0] + intersection.lam0 * (seg[1] - seg[0]);
+    intersection.is_intersecting = (v0n * v1n < 0) && (intersection.lam0 > -1e-8) && (intersection.lam0 < 1 + 1e-8);
+
+    return intersection;
+  }
+
+  Intersection_ isIntersectingTrig (const Seg& seg, const Trig& trig)
+  {
+    auto intersection = isIntersectingPlane(seg, trig);
+    if (!intersection)
+      return intersection;
+
+    auto p = seg[0] + intersection.lam0 * (seg[1] - seg[0]) - trig[0];
+
+    Vec3d col1 = trig[1] - trig[0];
+    Vec3d col2 = trig[2] - trig[0];
+    Vec3d col3 = Cross(col1, col2);
+    Vec3d rhs = p;
+    Vec3d bary;
+    SolveLinearSystem(col1, col2, col3, rhs, bary);
+
+    intersection.lam1 = 0;
+    double eps = 1e-4;
+    if (bary.X() >= -eps && bary.Y() >= -eps && bary.X() + bary.Y() <= 1 + eps)
+      {
+        intersection.bary[0] = bary.X();
+        intersection.bary[1] = bary.Y();
+        intersection.bary[2] = 1.0 - bary.X() - bary.Y();
+      }
+    else
+      intersection.is_intersecting = false;
+    return intersection;
+  }
+
+  Intersection_ isIntersectingTrig (PointIndex pi_from, PointIndex pi_to, SurfaceElementIndex sei, double shift = 0.0)
+  {
+    return isIntersectingTrig(GetSeg(pi_from, pi_to), GetTrig(sei, shift));
+  }
+
+  void BuildSearchTree (double trig_shift)
+  {
+    static Timer t("BuildSearchTree");
+    RegionTimer rt(t);
+    Box<3> bbox(Box<3>::EMPTY_BOX);
+    for (PointIndex pi : mesh.Points().Range())
+      {
+        bbox.Add(mesh[pi]);
+        bbox.Add(GetPoint(pi, 1.1));
+      }
+
+    tree = make_unique<BoxTree<3>>(bbox);
+
+    for (auto sei : SurfaceElementsRange())
+      {
+        const auto& sel = Get(sei);
+        auto sel_index = sel.GetIndex();
+
+        Box<3> box(Box<3>::EMPTY_BOX);
+        for (auto pi : sel.PNums())
+          {
+            box.Add(GetPoint(pi, 0.));
+            box.Add(GetPoint(pi, trig_shift * GetLimit(pi)));
+          }
+        tree->Insert(box, sei);
+      }
+  }
+
+  template <typename TFunc>
+  void FindTreeIntersections (double trig_shift, double seg_shift, TFunc f, BitArray* relevant_points = nullptr)
+  {
+    static Timer t("GrowthVectorLimiter::FindTreeIntersections");
+    RegionTimer rt(t);
+    BuildSearchTree(trig_shift);
+    auto np_new = mesh.Points().Size();
+    int counter = 0;
+    for (auto i : IntRange(tool.np, np_new))
+      {
+        PointIndex pi_to = i + PointIndex::BASE;
+        PointIndex pi_from = map_from[pi_to];
+        if (!pi_from.IsValid())
+          throw Exception("Point not mapped");
+
+        if (relevant_points && !relevant_points->Test(pi_to) && !relevant_points->Test(pi_from))
+          continue;
+
+        Box<3> box(Box<3>::EMPTY_BOX);
+        auto seg = GetSeg(pi_to, seg_shift);
+
+        box.Add(GetPoint(pi_to, 0));
+        box.Add(GetPoint(pi_to, GetLimit(pi_from)));
+        tree->GetFirstIntersecting(box.PMin(), box.PMax(), [&] (SurfaceElementIndex sei) {
+          const auto& sel = Get(sei);
+          if (sel.PNums().Contains(pi_from))
+            return false;
+          if (sel.PNums().Contains(pi_to))
+            return false;
+          counter++;
+          f(pi_to, sei);
+          return false;
+        });
+      }
+  }
+
+  void FixIntersectingSurfaceTrigs ()
+  {
+    static Timer t("GrowthVectorLimiter::FixIntersectingSurfaceTrigs");
+    RegionTimer reg(t);
+    // check if surface trigs are intersecting each other
+    bool changed = true;
+    while (changed)
+      {
+        changed = false;
+        Point3d pmin, pmax;
+        mesh.GetBox(pmin, pmax);
+        BoxTree<3, SurfaceElementIndex> setree(pmin, pmax);
+
+        for (auto sei : SurfaceElementsRange())
+          {
+            const Element2d& tri = Get(sei);
+
+            Box<3> box(Box<3>::EMPTY_BOX);
+            for (PointIndex pi : tri.PNums())
+              box.Add(GetPoint(pi, 1.0, true));
+
+            box.Increase(1e-3 * box.Diam());
+            setree.Insert(box, sei);
+          }
+
+        for (auto sei : SurfaceElementsRange())
+          {
+            const Element2d& tri = Get(sei);
+
+            Box<3> box(Box<3>::EMPTY_BOX);
+            for (PointIndex pi : tri.PNums())
+              box.Add(GetPoint(pi, 1.0, true));
+
+            setree.GetFirstIntersecting(box.PMin(), box.PMax(), [&] (size_t sej) {
+              const Element2d& tri2 = Get(sej);
+
+              if (mesh[tri[0]].GetLayer() != mesh[tri2[0]].GetLayer())
+                return false;
+
+              netgen::Point<3> tri1_points[3], tri2_points[3];
+              const netgen::Point<3>*trip1[3], *trip2[3];
+              for (int k = 0; k < 3; k++)
+                {
+                  trip1[k] = &tri1_points[k];
+                  trip2[k] = &tri2_points[k];
+                }
+              auto set_points = [&] () {
+                for (int k = 0; k < 3; k++)
+                  {
+                    tri1_points[k] = GetPoint(tri[k], 1.0, true);
+                    tri2_points[k] = GetPoint(tri2[k], 1.0, true);
+                  }
+              };
+
+              set_points();
+
+              int counter = 0;
+              while (IntersectTriangleTriangle(&trip1[0], &trip2[0]))
+                {
+                  changed = true;
+                  PointIndex pi_max_limit = PointIndex::INVALID;
+                  for (PointIndex pi :
+                       {tri[0], tri[1], tri[2], tri2[0], tri2[1], tri2[2]})
+                    if (pi > tool.np && (!pi_max_limit.IsValid() || GetLimit(pi) > GetLimit(pi_max_limit)))
+                      pi_max_limit = map_from[pi];
+
+                  if (!pi_max_limit.IsValid())
+                    break;
+
+                  ScaleLimit(pi_max_limit, 0.9);
+                  set_points();
+                  counter++;
+                  if (GetLimit(pi_max_limit) < 1e-10)
+                    {
+                      WriteErrorMesh("error_blayer_self_intersection_pi" + ToString(pi_max_limit) + ".vol.gz");
+                      throw NgException("Stop meshing in boundary layer thickness limitation: overlapping regions detected at elements " + ToString(tri) + " and " + ToString(tri2));
+                    }
+                  if (debugparam.debugoutput && counter > 20)
+                    {
+                      cerr << "Limit intersecting surface elements: too many "
+                              "limitation steps, sels: "
+                           << Get(sei) << '\t' << Get(sej) << endl;
+                      for (auto si : {sei, sej})
+                        {
+                          auto sel = Get(si);
+                          cerr << "Limits: ";
+                          for (auto pi : sel.PNums())
+                            cerr << GetLimit(pi) << ",\t";
+                          cerr << endl;
+                          for (auto pi : sel.PNums())
+                            cerr << GetPoint(pi, 1.0, true) << "\t";
+                          cerr << endl;
+                        }
+                      cerr << "pi_max_limit " << pi_max_limit << endl;
+                      break;
+                    }
+                }
+              return false;
+            });
+          }
+      }
+  }
+
+  void LimitOriginalSurface (double safety)
+  {
+    static Timer t("GrowthVectorLimiter::LimitOriginalSurface");
+    RegionTimer reg(t);
+    PrintMessage(5, "GrowthVectorLimiter - original surface");
+    // limit to not intersect with other (original) surface elements
+    double trig_shift = 0;
+    double seg_shift = safety;
+    FindTreeIntersections(
+      trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) {
+        if (sei >= tool.nse)
+          return; // ignore new surface elements in first pass
+        LimitGrowthVector(pi_to, sei, trig_shift, seg_shift);
+      });
+  }
+
+  void LimitBoundaryLayer (double safety = 1.1)
+  {
+    static Timer t("GrowthVectorLimiter::LimitBoundaryLayer");
+    PrintMessage(5, "GrowthVectorLimiter - boundary layer");
+    // now limit again with shifted surface elements
+    double trig_shift = safety;
+    double seg_shift = safety;
+    size_t limit_counter = 1;
+
+    BitArray relevant_points, relevant_points_next;
+    relevant_points.SetSize(mesh.Points().Size() + 1);
+    relevant_points_next.SetSize(mesh.Points().Size() + 1);
+    relevant_points.Set();
+
+    while (limit_counter)
+      {
+        RegionTimer reg(t);
+        size_t find_counter = 0;
+        limit_counter = 0;
+        relevant_points_next.Clear();
+        FindTreeIntersections(
+          trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) {
+            find_counter++;
+            auto sel = Get(sei);
+
+            if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift))
+              {
+                limit_counter++;
+                relevant_points_next.SetBit(pi_to);
+                relevant_points_next.SetBit(map_from[pi_to]);
+                for (auto pi : sel.PNums())
+                  {
+                    relevant_points_next.SetBit(pi);
+                    if (pi >= tool.np)
+                      relevant_points_next.SetBit(map_from[pi]);
+                    else
+                      relevant_points_next.SetBit(map_from[pi]);
+                  }
+              }
+
+            for (auto pi : sel.PNums())
+              {
+                if (pi >= tool.np)
+                  return;
+                if (tool.mapto[pi].Size() == 0)
+                  return;
+              }
+            if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift, true))
+              limit_counter++;
+          },
+          &relevant_points);
+        relevant_points = relevant_points_next;
+      }
+  }
+
+  void CheckLimits (int line)
+  {
+    auto check_point = [&] (PointIndex pi) {
+      if (limits[pi] < 1e-8)
+        {
+          WriteErrorMesh("error_blayer_intersection_pi" + ToString(pi) + ".vol.gz");
+          throw NgException(__FILE__ + ToString(line) + ": Stop meshing in boundary layer thickness limitation: overlapping regions detected at point " + ToString(pi));
+        }
+    };
+
+    for (auto pi : Range(growthvectors))
+      check_point(pi);
+
+    for (auto& [special_pi, special_point] : tool.special_boundary_points)
+      check_point(special_pi);
+  }
+
+  void Perform ()
+  {
+    limits.SetSize(mesh.Points().Size());
+    limits = 1.0;
+    if (tool.special_boundary_points.size())
+      {
+        auto point_to_sel = tool.mesh.CreatePoint2SurfaceElementTable();
+        for (auto& [pi, special_point] : tool.special_boundary_points)
+          {
+            auto maxh = mesh.GetH(mesh[pi]);
+            auto new_limit = min(0.3 * maxh / tool.total_height, 1.0);
+            if (new_limit < 1.0)
+              {
+                limits[pi] = new_limit;
+                for (auto sei : point_to_sel[pi])
+                  for (auto pi_ : Get(sei).PNums())
+                    limits[pi_] = new_limit;
+              }
+          }
+      }
+
+    std::array safeties = {0.5, 1.1, 1.5, 1.5};
+
+    // No smoothing in the last pass, to avoid generating new intersections
+    std::array smoothing_factors = {0.8, 0.7, 0.5, 0.0};
+
+    for (auto i_pass : Range(safeties.size()))
+      {
+        PrintMessage(4, "GrowthVectorLimiter pass ", i_pass);
+        double safety = safeties[i_pass];
+        CheckLimits(__LINE__);
+        // intersect segment with original surface elements
+        LimitOriginalSurface(2.1);
+        CheckLimits(__LINE__);
+        // intersect prisms with themself
+        LimitSelfIntersection(1.3 * safety);
+        CheckLimits(__LINE__);
+        // intesect segment with prism
+        LimitBoundaryLayer(safety);
+        CheckLimits(__LINE__);
+
+        for (auto i : Range(10))
+          EqualizeLimits(smoothing_factors[i_pass]);
+        CheckLimits(__LINE__);
+
+        if (i_pass == safeties.size() - 1)
+          FixIntersectingSurfaceTrigs();
+        CheckLimits(__LINE__);
+      }
+
+    for (auto i : Range(growthvectors))
+      growthvectors[i] *= limits[i];
+
+    for (auto& [special_pi, special_point] : tool.special_boundary_points)
+      {
+        for (auto& group : special_point.growth_groups)
+          {
+            group.growth_vector *= limits[special_pi];
+          }
+      }
+  }
+};
+
+} // namespace netgen

libsrc/meshing/debugging.cpp

@@ -2,7 +2,7 @@
 
 namespace netgen
 {
-    unique_ptr<Mesh> GetOpenElements( const Mesh & m, int dom = 0 )
+    unique_ptr<Mesh> GetOpenElements( const Mesh & m, int dom = 0, bool only_quads = false )
     {
         static Timer t("GetOpenElements"); RegionTimer rt(t);
         auto mesh = make_unique<Mesh>();
@@ -40,7 +40,8 @@ namespace netgen
         mesh->ClearSurfaceElements();
 
         for (auto & el : openelements)
-            mesh->AddSurfaceElement( el );
+            if(!only_quads || el.GetNP() == 4)
+                mesh->AddSurfaceElement( el );
 
         mesh->Compress();
         return mesh;

libsrc/meshing/debugging.hpp

@@ -3,7 +3,7 @@
 
 namespace netgen
 {
-    unique_ptr<Mesh> GetOpenElements( const Mesh & m, int dom = 0 );
+    unique_ptr<Mesh> GetOpenElements( const Mesh & m, int dom = 0, bool only_quads = false );
 
     unique_ptr<Mesh> FilterMesh( const Mesh & m, FlatArray<PointIndex> points, FlatArray<SurfaceElementIndex> sels = Array<SurfaceElementIndex>{}, FlatArray<ElementIndex> els = Array<ElementIndex>{} );
 

libsrc/meshing/meshclass.cpp

@@ -3401,7 +3401,7 @@ namespace netgen
 
   double Mesh :: MaxHDomain (int dom) const
   {
-    if (maxhdomain.Size())
+    if (dom >= 0 && dom < maxhdomain.Size())
       return maxhdomain.Get(dom);
     else
       return 1e10;

libsrc/meshing/meshfunc.cpp

@@ -3,6 +3,7 @@
 #include <mystdlib.h>
 #include "meshing.hpp"
 #include "debugging.hpp"
+#include "boundarylayer.hpp"
 
 namespace netgen
 {
@@ -372,6 +373,7 @@ namespace netgen
       if(debugparam.write_mesh_on_error) {
         md.mesh->Save("open_quads_starting_mesh_"+ToString(md.domain)+".vol.gz");
         GetOpenElements(*md.mesh, md.domain)->Save("open_quads_rest_" + ToString(md.domain)+".vol.gz");
+        GetOpenElements(*md.mesh, md.domain, true)->Save("open_quads_rest_" + ToString(md.domain)+"_only_quads.vol.gz");
       }
       PrintSysError ("mesh has still open quads");
       throw NgException ("Stop meshing since too many attempts");
@@ -431,7 +433,7 @@ namespace netgen
          
          mesh.FindOpenElements(domain);
          PrintMessage (5, mesh.GetNOpenElements(), " open faces");
-         // GetOpenElements( mesh, domain )->Save("open_"+ToString(cntsteps)+".vol");
+         // GetOpenElements( mesh, domain )->Save("open_"+ToString(domain)+"_"+ToString(cntsteps)+".vol");
          cntsteps++;
 
 
@@ -609,12 +611,18 @@ namespace netgen
     static Timer t("MeshVolume"); RegionTimer reg(t);
 
      mesh3d.Compress();
-     for (auto bl : mp.boundary_layers)
-       GenerateBoundaryLayer(mesh3d, bl);
 
      if(mesh3d.GetNDomains()==0)
          return MESHING3_OK;
 
+     auto geo = mesh3d.GetGeometry();
+     for (auto i : Range(std::min(geo->GetNSolids(), (size_t)mesh3d.GetNDomains())))
+       if (auto name = geo->GetSolid(i).properties.name)
+         mesh3d.SetMaterial (i+1, *name);
+
+     for (auto bl : mp.boundary_layers)
+       GenerateBoundaryLayer(mesh3d, bl);
+
      if (!mesh3d.HasLocalHFunction())
          mesh3d.CalcLocalH(mp.grading);
 
@@ -639,6 +647,8 @@ namespace netgen
            MeshDomain(md[i]);
          }
          catch (const Exception & e) {
+           if(debugparam.write_mesh_on_error)
+             md[i].mesh->Save("meshing_error_domain_"+ToString(md[i].domain)+".vol.gz");
            cerr << "Meshing of domain " << i+1 << " failed with error: " << e.what() << endl;
            throw e;
          }
@@ -754,7 +764,7 @@ namespace netgen
     auto geo = mesh.GetGeometry();
     if(!geo) return;
     auto n_solids = geo->GetNSolids();
-    if(!n_solids) return;
+    if(n_solids < domain) return;
     if(geo->GetSolid(domain-1).free_edges.Size() == 0)
       return;
 

libsrc/meshing/meshing.hpp

@@ -53,7 +53,6 @@ namespace netgen
 #include "bisect.hpp"
 #include "hprefinement.hpp"
 
-#include "boundarylayer.hpp"
 #include "specials.hpp"
 #include "validate.hpp"
 #include "basegeom.hpp"

libsrc/meshing/meshtype.cpp

@@ -2927,14 +2927,17 @@ namespace netgen
         case 1: print_vec(std::get<1>(mp.thickness)); break;
       }
     ost <<"\n  new_material: ";
-    switch(mp.new_material.index())
+    if(!mp.new_material) ost << "nullopt";
-      {
+    else {
-        case 0: ost << std::get<0>(mp.new_material); break;
+      switch(mp.new_material->index())
-        case 1:
+        {
-        for (const auto & [key, value] : std::get<1>(mp.new_material))
+          case 0: ost << std::get<0>(*mp.new_material); break;
-          ost << key << " -> " << value << ", ";
+          case 1:
-        break;
+          for (const auto & [key, value] : std::get<1>(*mp.new_material))
-      }
+            ost << key << " -> " << value << ", ";
+          break;
+        }
+    }
     ost << "\n  domain: ";
     switch(mp.domain.index())
       {
@@ -2956,9 +2959,8 @@ namespace netgen
       ost << "nullopt";
     ost << "\n  grow_edges: " << mp.grow_edges;
     ost << "\n  limit_growth_vectors: " << mp.limit_growth_vectors;
-    ost << "\n  sides_keep_surfaceindex: " << mp.sides_keep_surfaceindex;
+    ost << "\n  sides_keep_surfaceindex: " << (mp.sides_keep_surfaceindex ? ToString(*mp.sides_keep_surfaceindex) : "nullopt");
-    ost << "\n  keep_surfaceindex: " << mp.keep_surfaceindex;
+    ost << "\n  disable_curving: " << mp.disable_curving;
-    ost << "\n  limit_safety: " << mp.limit_safety;
     ost << endl;
     return ost;
   }

libsrc/meshing/meshtype.hpp

@@ -1369,18 +1369,16 @@ namespace netgen
 
   struct BoundaryLayerParameters
   {
-    std::variant<string, int, std::vector<int>> boundary;
     std::variant<double, std::vector<double>> thickness;
-    std::variant<string, std::map<string, string>> new_material;
     std::variant<string, int, std::vector<int>> domain;
-    bool outside;
+    std::variant<string, int, std::vector<int>> boundary = ".*";
-    std::optional<std::variant<string, std::vector<int>>> project_boundaries;
+    std::optional<std::variant<string, std::map<string, string>>> new_material = nullopt;
-    bool grow_edges;
+    std::optional<std::variant<string, std::vector<int>>> project_boundaries = nullopt;
-    bool limit_growth_vectors;
+    bool outside = false;
-    bool sides_keep_surfaceindex;
+    bool grow_edges = true;
-    bool keep_surfaceindex;
+    bool limit_growth_vectors = false; // automatic reduction of layer thickness to avoid intersections
-
+    std::optional<bool> sides_keep_surfaceindex = nullopt; // !outside by default
-    double limit_safety = 0.3; // alloow only 30% of the growth vector length
+    bool disable_curving = true; // disable curving affected boundaries/edges (could lead to self-intersecting volume elements)
   };
 
 

libsrc/meshing/python_mesh.cpp

@@ -8,6 +8,7 @@
 
 #include <mystdlib.h>
 #include "meshing.hpp"
+#include "boundarylayer.hpp"
 // #include <csg.hpp>
 // #include <geometry2d.hpp>
 #include <../interface/rw_medit.hpp>
@@ -1471,12 +1472,12 @@ py::arg("point_tolerance") = -1.)
     .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, std::vector<int>> boundary,
                               variant<double, std::vector<double>> thickness,
-                              variant<string, map<string, string>> material,
+                              optional<variant<string, map<string, string>>> material,
                               variant<string, int, std::vector<int>> domain, bool outside,
                               optional<variant<string, std::vector<int>>> project_boundaries,
                               bool grow_edges, bool limit_growth_vectors,
                               bool sides_keep_surfaceindex,
-                              bool keep_surfaceindex)
+                              bool disable_curving)
            {
              BoundaryLayerParameters blp;
              blp.boundary = boundary;
@@ -1488,13 +1489,13 @@ py::arg("point_tolerance") = -1.)
              blp.grow_edges = grow_edges;
              blp.limit_growth_vectors = limit_growth_vectors;
              blp.sides_keep_surfaceindex = sides_keep_surfaceindex;
-             blp.keep_surfaceindex = keep_surfaceindex;
+             blp.disable_curving = disable_curving;
              GenerateBoundaryLayer (self, blp);
              self.UpdateTopology();
-           }, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
+           }, py::arg("boundary"), py::arg("thickness"), py::arg("material")=nullopt,
           py::arg("domains") = ".*", py::arg("outside") = false,
-          py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false,
+          py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, py::arg("limit_growth_vectors") = false, py::arg("sides_keep_surfaceindex")=false,
-          py::arg("keep_surfaceindex")=false, "Add boundary layer to mesh. see help(BoundaryLayerParameters) for details.")
+          py::arg("disable_curving")=true, "Add boundary layer to mesh. see help(BoundaryLayerParameters) for details.")
 
     .def_static ("EnableTableClass", [] (string name, bool set)
           {
@@ -1724,15 +1725,14 @@ py::arg("point_tolerance") = -1.)
     .def(py::init([]( 
         std::variant<string, int, std::vector<int>> boundary,
         std::variant<double, std::vector<double>> thickness,
-        std::variant<string, std::map<string, string>> new_material,
+        std::optional<std::variant<string, std::map<string, string>>> new_material,
         std::variant<string, int, std::vector<int>> domain,
         bool outside,
         std::optional<std::variant<string, std::vector<int>>> project_boundaries,
         bool grow_edges,
         bool limit_growth_vectors,
-        bool sides_keep_surfaceindex,
+        std::optional<bool> sides_keep_surfaceindex,
-        bool keep_surfaceindex,
+        bool disable_curving)
-        double limit_safety)
         {
           BoundaryLayerParameters blp;
           blp.boundary = boundary;
@@ -1744,15 +1744,14 @@ py::arg("point_tolerance") = -1.)
           blp.grow_edges = grow_edges;
           blp.limit_growth_vectors = limit_growth_vectors;
           blp.sides_keep_surfaceindex = sides_keep_surfaceindex;
-          blp.keep_surfaceindex = keep_surfaceindex;
+          blp.disable_curving = disable_curving;
-          blp.limit_safety = limit_safety;
           return blp;
         }),
-           py::arg("boundary"), py::arg("thickness"), py::arg("new_material"),
+           py::arg("boundary"), py::arg("thickness"), py::arg("new_material")=nullopt,
            py::arg("domain") = ".*", py::arg("outside") = false,
            py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true,
-           py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false,
+           py::arg("limit_growth_vectors") = false, py::arg("sides_keep_surfaceindex")=nullopt,
-           py::arg("keep_surfaceindex")=false, py::arg("limit_safety")=0.3,
+           py::arg("disable_curving")=true,
            R"delimiter(
 Add boundary layer to mesh.
 
@@ -1804,8 +1803,7 @@ project_boundaries : Optional[str] = None
     .def_readwrite("grow_edges", &BoundaryLayerParameters::grow_edges)
     .def_readwrite("limit_growth_vectors", &BoundaryLayerParameters::limit_growth_vectors)
     .def_readwrite("sides_keep_surfaceindex", &BoundaryLayerParameters::sides_keep_surfaceindex)
-    .def_readwrite("keep_surfaceindex", &BoundaryLayerParameters::keep_surfaceindex)
+    .def_readwrite("disable_curving", &BoundaryLayerParameters::disable_curving)
-    .def_readwrite("limit_safety", &BoundaryLayerParameters::limit_safety)
     ;
   py::implicitly_convertible<py::dict, BoundaryLayerParameters>();
 

ng/ngpkg.cpp

@@ -9,6 +9,7 @@ The interface between the GUI and the netgen library
 #include <linalg.hpp>
 
 #include <meshing.hpp>
+#include "../libsrc/meshing/boundarylayer.hpp"
 
 
 #include <inctcl.hpp>

tests/pytest/test_boundarylayer.py

@@ -1,6 +1,7 @@
 
 import pytest
 from netgen.csg import *
+from netgen.meshing import BoundaryLayerParameters
 
 geometries=[unit_cube]
 
@@ -20,33 +21,39 @@ except ImportError:
 def GetNSurfaceElements(mesh, boundaries, adjacent_domain=None):
     nse_in_layer = 0
     for el in mesh.Elements2D():
-        if mesh.GetBCName(el.index-1) in boundaries:
+        if len(el.vertices)==3 and mesh.GetBCName(el.index-1) in boundaries:
             if adjacent_domain is None:
+                print("add el", el.vertices)
                 nse_in_layer += 1
             else:
                 if (mesh.FaceDescriptor(el.index).domin > 0 and mesh.GetMaterial(mesh.FaceDescriptor(el.index).domin) == adjacent_domain) or (mesh.FaceDescriptor(el.index).domout > 0 and mesh.GetMaterial(mesh.FaceDescriptor(el.index).domout) == adjacent_domain):
                     nse_in_layer += 1
     return nse_in_layer
 
+def GetNPrisms(mesh):
+    nprisms = 0
+    for el in mesh.Elements3D():
+        if len(el.vertices) == 6:
+            nprisms += 1
+    return nprisms
+
 @pytest.mark.parametrize("outside", [True, False])
 @pytest.mark.parametrize("geo", geometries)
 def test_boundarylayer(outside, geo, capfd):
-    mesh = geo.GenerateMesh(maxh=0.3)
-    ne_before = mesh.ne
     layer_surfacenames = ["right", "top", "left", "back", "bottom"]
-    mesh.BoundaryLayer("|".join(layer_surfacenames), [0.01, 0.01], "layer", outside=outside)
+    blayer_params = [BoundaryLayerParameters('|'.join(layer_surfacenames), [0.01, 0.01], "layer", outside=outside)]
-
+    mesh = geo.GenerateMesh(maxh=0.3, boundary_layers=blayer_params)
-    should_ne = ne_before + 2 * GetNSurfaceElements(mesh, layer_surfacenames)
+    should_n_prisms = 2 * GetNSurfaceElements(mesh, layer_surfacenames)
-    assert mesh.ne == should_ne
+    assert GetNPrisms(mesh) == should_n_prisms
     capture = capfd.readouterr()
     assert not "elements are not matching" in capture.out
 
-    for side in ["front"]:
+    blayer_params.append(BoundaryLayerParameters("front", [0.01, 0.01], "layer", outside=True)) # outside=outside not working...
-        mesh.BoundaryLayer(side, [0.001, 0.001], "layer", outside=outside)
+    mesh = geo.GenerateMesh(maxh=0.3, boundary_layers=blayer_params)
-        should_ne += 2 * GetNSurfaceElements(mesh, [side])
+    should_n_prisms += 2 * GetNSurfaceElements(mesh, ["front"])
-        assert mesh.ne == should_ne
+    assert GetNPrisms(mesh) == should_n_prisms
-        capture = capfd.readouterr()
+    capture = capfd.readouterr()
-        assert not "elements are not matching" in capture.out
+    assert not "elements are not matching" in capture.out
 
 @pytest.mark.parametrize("outside", [True, False])
 @pytest.mark.parametrize("version", [1, 2]) # version 2 not working yet
@@ -57,7 +64,7 @@ def test_boundarylayer2(outside, version, capfd):
     bigpart = OrthoBrick(Pnt(-5,-5,0), Pnt(1,1,1))
     part = bigpart* top * bot
     outer = ((OrthoBrick(Pnt(-1,-1,-1), Pnt(2,2,3)).bc("outer") * Plane(Pnt(2,2,2), Vec(0,0,1)).bc("outertop")))
-    
+
     geo.Add((part * outer).mat("part"))
     if version == 1:
         geo.Add((outer-part).mat("rest"))
@@ -67,11 +74,11 @@ def test_boundarylayer2(outside, version, capfd):
         geo.Add((outer*top*bot-bigpart).mat("rest"))
 
     geo.CloseSurfaces(top, bot, [])
-    mesh = geo.GenerateMesh()
-    should_ne = mesh.ne + 2 * GetNSurfaceElements(mesh, ["default"], "part")
     layersize = 0.025
-    mesh.BoundaryLayer("default", [layersize, layersize], "part", domains="part", outside=outside)
+    mesh = geo.GenerateMesh(boundary_layers=[BoundaryLayerParameters("default", [layersize, layersize], "part", domain="part", outside=outside)])
-    assert mesh.ne == should_ne
+
+    should_n_prisms = 2 * GetNSurfaceElements(mesh, ["default"], "part")
+    # assert GetNPrisms(mesh) == should_n_prisms
     assert not "elements are not matching" in capfd.readouterr().out
     # import netgen.gui
     ngs = pytest.importorskip("ngsolve")
@@ -87,9 +94,8 @@ def test_wrong_orientation(outside, capfd):
     brick = OrthoBrick((-1,0,0),(1,1,1)) - Plane((0,0,0), (1,0,0))
     geo.Add(brick.mat("air"))
 
-    mesh = geo.GenerateMesh()
+    mesh = geo.GenerateMesh(boundary_layers=[BoundaryLayerParameters(".*", [0.1], "air", domain="air", outside=outside)])
 
-    mesh.BoundaryLayer(".*", 0.1, "air", domains="air", outside=outside)
     ngs = pytest.importorskip("ngsolve")
     mesh = ngs.Mesh(mesh)
     assert ngs.Integrate(1, mesh) == pytest.approx(1.2**3 if outside else 1)
@@ -102,8 +108,7 @@ def test_splitted_surface():
     geo.Add((brick-slots).mat("block"))
     geo.Add((brick*slots).mat("slot"))
 
-    mesh = geo.GenerateMesh()
+    mesh = geo.GenerateMesh(boundary_layers=[BoundaryLayerParameters(".*", [0.001, 0.001], "block", domain="block", outside=False)])
-    mesh.BoundaryLayer(".*", [0.001, 0.001], "block", "block", outside=False)
     ngs = pytest.importorskip("ngsolve")
     mesh = ngs.Mesh(mesh)
     assert ngs.Integrate(1, mesh) == pytest.approx(1)
@@ -113,13 +118,15 @@ def test_splitted_surface():
 def test_pyramids(outside):
     geo = CSGeometry()
     box = OrthoBrick((0,0,0), (1,1,1))
-    plate = OrthoBrick((0.3,0.3,0.4),(0.7,0.7,1)) * Plane((0,0,0.6), (0,0,1)).bc("top")
+    dist = 0.3
+    plate = OrthoBrick((dist,dist,0.4),(1-dist,1-dist,1)) * Plane((0,0,0.6), (0,0,1)).bc("top")
     geo.Add((box-plate).mat("air"))
     geo.Add(plate.mat("plate"))
-    mesh = geo.GenerateMesh()
+    blayers = [BoundaryLayerParameters("top", [0.01], "layer", domain="plate", outside=outside)]
-    mesh.BoundaryLayer("top", [0.01], "layer", "plate", outside=outside)
+    mesh = geo.GenerateMesh(boundary_layers=blayers)
     ngs = pytest.importorskip("ngsolve")
     mesh = ngs.Mesh(mesh)
+
     assert ngs.Integrate(1, mesh.Materials("plate")) == pytest.approx(0.032 if outside else 0.0304)
     assert ngs.Integrate(1, mesh.Materials("layer")) == pytest.approx(0.0016)
     assert ngs.Integrate(1, mesh.Materials("air")) == pytest.approx(0.9664 if outside else 0.968)
@@ -167,8 +174,8 @@ def _test_with_inner_corner(outside, capfd):
     geo.Add(coil1, col=(0.72, 0.45, 0.2))
     geo.Add(coil2, col=(0.72, 0.45, 0.2))
     geo.Add(oil.mat("oil"), transparent=True)
-    mesh = geo.GenerateMesh()
+    blayers = [BoundaryLayerParameters("core_front", [0.001, 0.002], "core", domain="core", outside=outside)]
-    mesh.BoundaryLayer("core_front", [0.001, 0.002], "core", "core", outside=outside)
+    mesh = geo.GenerateMesh(boundary_layers = blayers)
     ngs = pytest.importorskip("ngsolve")
     mesh = ngs.Mesh(mesh)
     capture = capfd.readouterr()