add project to boundary in boundarylayer and correctly treat inverse boundaries

libsrc/meshing/boundarylayer.cpp

@@ -215,6 +215,11 @@ namespace netgen
             if (blp.surfid.Contains(sel.GetIndex()))
               {
                 auto n2 = GetSurfaceNormal(mesh,sel);
+                auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
+                auto domin = fd.DomainIn();
+                auto domout = fd.DomainOut();
+                if(blp.domains.Test(domout) && !blp.domains.Test(domin))
+                  n2 *= -1;
                 if(!blp.outside)
                   n2 *= -1;
                 for(auto pi : sel.PNums())
@@ -238,24 +243,37 @@ namespace netgen
               }
 
           // project growthvector on surface for inner angles
-          // for(const auto& sel : mesh.SurfaceElements())
+          for(const auto& sel : mesh.SurfaceElements())
-          //   if(!blp.surfid.Contains(sel.GetIndex()))
+            // if(!blp.surfid.Contains(sel.GetIndex()))
-          //     {
+            //   {
-          //       auto n = GetSurfaceNormal(mesh, sel);
+            if(blp.project_boundaries.Contains(sel.GetIndex()))
-          //       for(auto pi : sel.PNums())
+              {
-          //         {
+                auto n = GetSurfaceNormal(mesh, sel);
-          //           if(growthvectors[pi].Length2() == 0.)
+                for(auto i : Range(sel.PNums()))
-          //             continue;
+                  {
-          //           auto& g = growthvectors[pi];
+                    auto pi = sel.PNums()[i];
-          //           auto ng = n * g;
+                    if(growthvectors[pi].Length2() == 0.)
-          //           auto gg = g * g;
+                      continue;
-          //           auto nn = n * n;
+                    auto next = sel.PNums()[(i+1)%sel.GetNV()];
-          //           if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue;
+                    auto prev = sel.PNums()[i == 0 ? sel.GetNV()-1 : i-1];
-          //           auto a = -ng*ng/(ng*ng-nn * gg);
+                    auto v1 = (mesh[next] - mesh[pi]).Normalize();
-          //           auto b = ng*gg/(ng*ng-nn*gg);
+                    auto v2 = (mesh[prev] - mesh[pi]).Normalize();
-          //           g += a*g + b*n;
+                    auto v3 = growthvectors[pi];
-          //         }
+                    v3.Normalize();
-          //     }
+                    // only project if sel goes in boundarylayer direction
+                    if((v1 * v3 < 1e-12) || (v2 * v3 < 1e-12))
+                      continue;
+
+                    auto& g = growthvectors[pi];
+                    auto ng = n * g;
+                    auto gg = g * g;
+                    auto nn = n * n;
+                    if(fabs(ng*ng-nn*gg) < 1e-12 || fabs(ng) < 1e-12) continue;
+                    auto a = -ng*ng/(ng*ng-nn * gg);
+                    auto b = ng*gg/(ng*ng-nn*gg);
+                    g += a*g + b*n;
+                  }
+              }
 
           if (!blp.grow_edges)
             {

libsrc/meshing/boundarylayer.hpp

@@ -18,6 +18,7 @@ public:
   BitArray domains;
   bool outside = false; // set the boundary layer on the outside
   bool grow_edges = false;
+  Array<size_t> project_boundaries;
 };
 
 DLL_HEADER void GenerateBoundaryLayer (Mesh & mesh,

libsrc/meshing/python_mesh.cpp

@@ -833,6 +833,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def("FaceDescriptor", static_cast<FaceDescriptor&(Mesh::*)(int)> (&Mesh::GetFaceDescriptor),
          py::return_value_policy::reference)
     .def("GetNFaceDescriptors", &Mesh::GetNFD)
+    .def("GetNDomains", &Mesh::GetNDomains)
 
     .def("GetVolumeNeighboursOfSurfaceElement", [](Mesh & self, size_t sel)
                                                 {
@@ -1012,7 +1013,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                               variant<double, py::list> thickness,
                               variant<string, py::list> material,
                               variant<string, int> domain, bool outside,
-                              bool grow_edges)
+                              bool grow_edges,
+                              optional<string> project_boundaries)
            {
              BoundaryLayerParameters blp;
              if(int* bc = get_if<int>(&boundary); bc)
@@ -1034,7 +1036,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                          if(dom_pattern)
                            {
                              regex pattern(*dom_pattern);
-                             if(regex_match(self.GetMaterial(fd.DomainIn()), pattern) || (fd.DomainOut() > 0 ? regex_match(self.GetMaterial(fd.DomainOut()), pattern) : false))
+                             if((fd.DomainIn() > 0 && regex_match(self.GetMaterial(fd.DomainIn()), pattern)) || (fd.DomainOut() > 0 && regex_match(self.GetMaterial(fd.DomainOut()), pattern)))
                                blp.surfid.Append(i);
                            }
                          else
@@ -1043,6 +1045,14 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                      }
                }
 
+             if(project_boundaries.has_value())
+               {
+                 regex pattern(*project_boundaries);
+                 for(int i = 1; i<=self.GetNFD(); i++)
+                   if(regex_match(self.GetFaceDescriptor(i).GetBCName(), pattern))
+                     blp.project_boundaries.Append(i);
+               }
+
              if(double* pthickness = get_if<double>(&thickness); pthickness)
                {
                  blp.heights.Append(*pthickness);
@@ -1102,7 +1112,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
              self.UpdateTopology();
            }, py::arg("boundary"), py::arg("thickness"), py::arg("material"),
           py::arg("domains") = ".*", py::arg("outside") = false,
-          py::arg("grow_edges") = false,
+          py::arg("grow_edges") = false, py::arg("project_boundaries")=nullopt,
           R"delimiter(
 Add boundary layer to mesh.
 
@@ -1127,6 +1137,11 @@ outside : bool = False
 grow_edges : bool = False
   Grow boundary layer over edges.
 
+project_boundaries : Optional[str] = None
+  Project boundarylayer to these boundaries if they meet them. Set
+  to boundaries that meet boundarylayer at a non-orthogonal edge and
+  layer-ending should be projected to that boundary.
+
 )delimiter")
 
     .def ("EnableTable", [] (Mesh & self, string name, bool set)

tests/pytest/test_boundarylayer.py

@@ -54,6 +54,29 @@ def test_boundarylayer2(outside, version, capfd):
     mesh = geo.GenerateMesh()
     should_ne = mesh.ne + 2 * GetNSurfaceElements(mesh, ["default"], "part")
     layersize = 0.05
-    mesh.BoundaryLayer("default", [0.5 * layersize, layersize], "layer", domains="part", outside=outside, grow_edges=True)
+    mesh.BoundaryLayer("default", [0.5 * layersize, layersize], "part", domains="part", outside=outside, grow_edges=True)
     assert mesh.ne == should_ne
     assert not "elements are not matching" in capfd.readouterr().out
+    import netgen.gui
+    ngs = pytest.importorskip("ngsolve")
+    ngs.Draw(ngs.Mesh(mesh))
+    mesh = ngs.Mesh(mesh)
+    assert ngs.Integrate(1, mesh.Materials("part")) == pytest.approx(0.5*2.05*2.05 if outside else 0.4*2*2)
+    assert ngs.Integrate(1, mesh) == pytest.approx(3**3)
+
+
+@pytest.mark.parametrize("outside", [True, False])
+def test_wrong_orientation(outside):
+    geo = CSGeometry()
+    brick = OrthoBrick((-1,0,0),(1,1,1)) - Plane((0,0,0), (1,0,0))
+    geo.Add(brick.mat("air"))
+
+    mesh = geo.GenerateMesh()
+
+    mesh.BoundaryLayer(".*", 0.1, "air", domains="air", outside=outside,
+                       grow_edges=True)
+    ngs = pytest.importorskip("ngsolve")
+    mesh = ngs.Mesh(mesh)
+    assert ngs.Integrate(1, mesh) == pytest.approx(1.2**3 if outside else 1)
+
+