cleanup, fixes

libsrc/meshing/boundarylayer.cpp

@@ -842,7 +842,7 @@ void BoundaryLayerTool ::InsertNewElements(
     auto &seg = segments[sei];
     if (is_boundary_moved.Test(seg.si)) {
       // cout << "moved setg " << seg << endl;
-       for (auto &p : seg.PNums())
+      for (auto &p : seg.PNums())
         if (hasMoved(p))
           p = newPoint(p);
     }
@@ -1155,18 +1155,16 @@ void BoundaryLayerTool ::Perform() {
   topo.SetBuildVertex2Element(true);
   mesh.UpdateTopology();
 
-  // cout << "GW 19911 " << growthvectors[19911] << endl;
   InterpolateGrowthVectors();
-  // cout << "GW 19911 " << growthvectors[19911] << endl;
   InterpolateSurfaceGrowthVectors();
-  // cout << "GW 19911 " << growthvectors[19911] << endl;
 
   AddSurfaceElements();
 
-
   if (params.limit_growth_vectors)
     LimitGrowthVectorLengths();
 
+  FixSurfaceElements();
+
   for (auto [pi, data] : growth_vector_map) {
     auto [gw, height] = data;
     mesh[pi] += height * (*gw);

libsrc/meshing/boundarylayer.hpp

@@ -89,6 +89,7 @@ class BoundaryLayerTool
     void InterpolateSurfaceGrowthVectors();
     void InterpolateGrowthVectors();
     void LimitGrowthVectorLengths();
+    void FixSurfaceElements();
 
     void InsertNewElements(FlatArray<Array<pair<SegmentIndex, int>>, SegmentIndex> segmap, const BitArray & in_surface_direction);
     void SetDomInOut();

libsrc/meshing/boundarylayer_interpolate.cpp

@@ -2,6 +2,65 @@
 
 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() {
   int new_max_edge_nr = max_edge_nr;
   for (const auto &seg : segments)
@@ -157,7 +216,7 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
   auto np_old = this->np;
   auto np = mesh.GetNP();
 
-  non_bl_growth_vectors.clear();
+  // non_bl_growth_vectors.clear();
   // for(const auto & sel : new_sels) {
   //   for(auto pi : sel.PNums()) {
   //     if(mesh[pi].Type() == INNERPOINT)
@@ -195,130 +254,34 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
 
   // cout << __FILE__ << ":" << __LINE__ << endl;
   std::set<PointIndex> points_set;
-  for (const auto &sel: mesh.SurfaceElements())
+  for (const auto &sel : mesh.SurfaceElements()) {
     for (auto pi : sel.PNums())
-      if(mesh[pi].Type() == SURFACEPOINT && hasMoved(pi))
+      if (mesh[pi].Type() == SURFACEPOINT && hasMoved(pi))
         points_set.insert(pi);
-
+  }
-  // Array<bool> has_moved_points(max_edge_nr + 1);
-  // has_moved_points = false;
-  // std::set<PointIndex> moved_edge_points;
-
-  // for (auto seg : segments) {
-  //   if (hasMoved(seg[0]) != hasMoved(seg[1]))
-  //     has_moved_points[seg.edgenr] = true;
-  // }
-  // cout << __FILE__ << ":" << __LINE__ << endl;
-
-  // for (auto seg : segments)
-  //   if (has_moved_points[seg.edgenr])
-  //     for (auto pi : seg.PNums())
-  //       if (mesh[pi].Type() == EDGEPOINT)
-  //         points_set.insert(pi);
 
   Array<PointIndex> points;
   for (auto pi : points_set)
     points.Append(pi);
   QuickSort(points);
 
-  // cout << __FILE__ << ":" << __LINE__ << endl;
-  // cout << "points to interpolate " << endl << points << endl;
-
-  // cout << __FILE__ << ":" << __LINE__ << endl;
-  auto p2sel = mesh.CreatePoint2SurfaceElementTable();
-  // cout << __FILE__ << ":" << __LINE__ << endl;
-  // auto p2sel = ngcore::CreateSortedTable<SurfaceElementIndex, PointIndex>(
-  //       new_sels.Range(),
-  //       [&](auto &table, SurfaceElementIndex ei) {
-  //         for (PointIndex pi : new_sels[ei].PNums())
-  //           table.Add(pi, ei);
-  //       },
-  //       mesh.GetNP());
   // smooth tangential part of growth vectors from edges to surface elements
   Array<Vec<3>, PointIndex> corrections(mesh.GetNP());
   corrections = 0.0;
   RegionTimer rtsmooth(tsmooth);
-  struct Neighbor {
+  auto neighbors = BuildNeighbors(points, mesh);
-    PointIndex pi;
-    SurfaceElementIndex sei;
-    double weight;
-  };
-  Array<ArrayMem<Neighbor, 20>> neighbors(points.Size());
-  // cout << __FILE__ << ":" << __LINE__ << endl;
-
-  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});
-        // if((mesh[pi1]-mesh[pi]).Length() < 1e-10) {
-        //   cout << "close points " << pi << "\t" << pi1 << endl;
-        // }
-        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;
-
-    // cout << "angles " << angles << endl;
-    // cout << "inv_dists " << inv_dists << endl;
-
-    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;
-
-  //   if(pi == 19911) {
-  //     cout << "pi " << pi << endl;
-  //     for(auto & nb : p_neighbors) {
-  //       cout << "neighbor " << nb.pi << "\t" << nb.weight << endl;
-  //     }
-  //     cout << "inv_dist " << inv_dists << endl;
-  //     cout << "angles " << angles << endl;
-  //   }
-  }
-  // cout << __FILE__ << ":" << __LINE__ << endl;
 
   Array<Vec<3>, SurfaceElementIndex> surf_normals(mesh.GetNSE());
   for (auto sei : mesh.SurfaceElements().Range())
     surf_normals[sei] = getNormal(mesh[sei]);
-  // cout << __FILE__ << ":" << __LINE__ << endl;
 
-  BitArray interpolate_tangent(mesh.GetNP()+1);
+  BitArray interpolate_tangent(mesh.GetNP() + 1);
-  // cout << __FILE__ << ":" << __LINE__ << endl;
   interpolate_tangent = false;
-  for(auto pi : points) {
+  for (auto pi : points) {
     for (auto sei : p2sel[pi])
-      if(is_boundary_moved[mesh[sei].GetIndex()])
+      if (is_boundary_moved[mesh[sei].GetIndex()])
         interpolate_tangent.SetBit(pi);
   }
-  // cout << __FILE__ << ":" << __LINE__ << endl;
 
   constexpr int N_STEPS = 64;
   for ([[maybe_unused]] auto i : Range(N_STEPS)) {
@@ -336,18 +299,14 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
       bool do_average_tangent = true;
       for (const auto &s : p_neighbors) {
         auto gw_other = getGW(s.pi) + corrections[s.pi];
-        // if(pi == 19911) cout << "neighbor " << s.pi << "\t" << s.weight << '\t' << gw_other << "\tdo avg: " << do_average_tangent << endl;
-        // if(pi == 19911) cout << "\tneighbor gw" << gw_other << endl;
         if (do_average_tangent) {
           auto n = surf_normals[s.sei];
           gw_other = gw_other - (gw_other * n) * n;
         }
-        // if(pi == 19911) cout << "\tneighbor gw" << gw_other << endl;
         auto v = gw_other;
         auto len = v.Length2();
         sum_len += len;
         max_len = max(max_len, len);
-        // if(pi == 19911) cout << "\tneighbor v" << v << endl;
         g_vectors.Append(v);
       }
 
@@ -365,13 +324,103 @@ void BoundaryLayerTool ::InterpolateSurfaceGrowthVectors() {
         auto v = g_vectors[i];
         double weight = lambda * p_neighbors[i].weight +
                         (1.0 - lambda) * v.Length2() / sum_len;
-        // if(pi == 19911) cout << "pi " << pi << "\tneighbor " << p_neighbors[i].pi << "\tweight " << weight << endl;
+        // if(pi == 19911) cout << "pi " << pi << "\tneighbor " <<
+        // p_neighbors[i].pi << "\tweight " << weight << endl;
         correction += weight * v;
       }
 
       if (!do_average_tangent)
         correction -= getGW(pi);
-      // if(pi == 19911) cout << "pi " << pi << "\tcorrection " << correction << endl;
+      // if(pi == 19911) cout << "pi " << pi << "\tcorrection " << correction <<
+      // endl;
+    }
+  }
+
+  for (auto pi : points)
+    addGW(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 (mesh[pi].Type() == 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;
+      }
     }
   }
 

libsrc/meshing/boundarylayer_limiter.hpp

@@ -1,5 +1,5 @@
-#include <core/array.hpp>
 #include "boundarylayer.hpp"
+#include <core/array.hpp>
 
 namespace netgen {
 
@@ -214,7 +214,8 @@ struct GrowthVectorLimiter {
       return;
     for (PointIndex pi : IntRange(tool.np, mesh.GetNP())) {
       auto pi_from = map_from[pi];
-      // if(pi_from == 21110) cout << "equalize " << pi << "\tfactor " << factor << endl;
+      // if(pi_from == 21110) cout << "equalize " << pi << "\tfactor " << factor
+      // << endl;
       std::set<PointIndex> pis;
       for (auto sei : p2sel[pi])
         for (auto pi_ : tool.new_sels[sei].PNums())
@@ -224,7 +225,8 @@ struct GrowthVectorLimiter {
         auto limit = GetLimit(pi1);
         if (limit > 0.0)
           limits.Append(GetLimit(pi1));
-        // if(pi_from == 21110) cout << "\tneighbor point " << map_from[pi1] << " -> " << pi1  << " with limit " << limit << endl;
+        // if(pi_from == 21110) cout << "\tneighbor point " << map_from[pi1] <<
+        // " -> " << pi1  << " with limit " << limit << endl;
       }
       // if(pi_from == 21110) cout << "\town limit " << GetLimit(pi) << endl;
       if (limits.Size() == 0)
@@ -351,7 +353,8 @@ struct GrowthVectorLimiter {
   }
 
   void BuildSearchTree(double trig_shift) {
-    static Timer t("BuildSearchTree"); RegionTimer rt(t);
+    static Timer t("BuildSearchTree");
+    RegionTimer rt(t);
     Box<3> bbox(Box<3>::EMPTY_BOX);
     for (PointIndex pi : mesh.Points().Range()) {
       bbox.Add(mesh[pi]);
@@ -470,7 +473,6 @@ struct GrowthVectorLimiter {
             set_points();
             counter++;
             if (counter > 20) {
-              break;
               cerr << "Limit intersecting surface elements: too many "
                       "limitation steps, sels: "
                    << Get(sei) << '\t' << Get(sej) << endl;