From f3150c5bbeb310b1eb72a40a6ddb28898375f703 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 12 Jul 2021 23:08:50 +0200
Subject: [PATCH] Improve MeshNodeBoundary and MeshFlatNodeBoundary

- improve error messages
- fixes #8 (parallel/sequential reproducibility and straight boundaries)
---
 src/mesh/MeshNodeBoundary.hpp | 248 +++++++++++++++++++---------------
 1 file changed, 140 insertions(+), 108 deletions(-)

diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 5846781f5..3e23be746 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -16,6 +16,7 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  protected:
   Array<const NodeId> m_node_list;
+  std::string m_boundary_name;
 
  public:
   MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
@@ -29,6 +30,7 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 
   template <typename MeshType>
   MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
+    : m_boundary_name(ref_face_list.refId().tagName())
   {
     static_assert(Dimension == MeshType::Dimension);
     const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
@@ -38,7 +40,10 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
       face_list.size(), PUGS_LAMBDA(int l) {
         const auto& face_cells = face_to_cell_matrix[face_list[l]];
         if (face_cells.size() > 1) {
-          throw NormalError("internal faces cannot be used to define mesh boundaries");
+          std::ostringstream ost;
+          ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+              << ": inner faces cannot be used to define mesh boundaries";
+          throw NormalError(ost.str());
         }
       });
 
@@ -66,7 +71,9 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   }
 
   template <typename MeshType>
-  MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list) : m_node_list(ref_node_list.list())
+  MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
+    : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName())
+
   {
     static_assert(Dimension == MeshType::Dimension);
   }
@@ -92,10 +99,27 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclu
   PUGS_INLINE Rd _getNormal(const MeshType& mesh);
 
   template <typename MeshType>
-  PUGS_INLINE void _checkBoundaryIsFlat(TinyVector<2, double> normal,
-                                        TinyVector<2, double> xmin,
-                                        TinyVector<2, double> xmax,
-                                        const MeshType& mesh) const;
+  PUGS_INLINE void
+  _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
+                       const TinyVector<Dimension, double>& origin,
+                       const double length,
+                       const MeshType& mesh) const
+  {
+    static_assert(MeshType::Dimension == Dimension);
+
+    const NodeValue<const Rd>& xr = mesh.xr();
+
+    parallel_for(
+      this->m_node_list.size(), PUGS_LAMBDA(const size_t r) {
+        const Rd& x = xr[this->m_node_list[r]];
+        if (dot(x - origin, normal) > 1E-13 * length) {
+          std::ostringstream ost;
+          ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+              << ": boundary is not flat!";
+          throw NormalError(ost.str());
+        }
+      });
+  }
 
   template <typename MeshType>
   PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh);
@@ -130,31 +154,6 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclu
   virtual ~MeshFlatNodeBoundary()                   = default;
 };
 
-template <>
-template <typename MeshType>
-void
-MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
-                                              TinyVector<2, double> xmin,
-                                              TinyVector<2, double> xmax,
-                                              const MeshType& mesh) const
-{
-  static_assert(MeshType::Dimension == 2);
-  using R2 = TinyVector<2, double>;
-
-  const R2 origin     = 0.5 * (xmin + xmax);
-  const double length = l2Norm(xmax - xmin);
-
-  const NodeValue<const R2>& xr = mesh.xr();
-
-  parallel_for(
-    m_node_list.size(), PUGS_LAMBDA(size_t r) {
-      const R2& x = xr[m_node_list[r]];
-      if (dot(x - origin, normal) > 1E-13 * length) {
-        throw NormalError("this FlatBoundary is not flat!");
-      }
-    });
-}
-
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<1, double>
@@ -173,7 +172,10 @@ MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
   }();
 
   if (number_of_bc_nodes != 1) {
-    throw NormalError("Node boundaries in 1D require to have exactly 1 node");
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << m_boundary_name << rang::style::reset
+        << ": node boundaries in 1D require to have exactly 1 node";
+    throw NormalError(ost.str());
   }
 
   return R{1};
@@ -190,46 +192,50 @@ MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
   const NodeValue<const R2>& xr = mesh.xr();
 
   R2 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
-
   R2 xmax(-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max());
 
-  for (size_t r = 0; r < m_node_list.size(); ++r) {
-    const R2& x = xr[m_node_list[r]];
+  auto update_xmin = [](const R2& x, R2& xmin) {
     if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
       xmin = x;
     }
+  };
+
+  auto update_xmax = [](const R2& x, R2& xmax) {
     if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
       xmax = x;
     }
+  };
+
+  for (size_t r = 0; r < m_node_list.size(); ++r) {
+    const R2& x = xr[m_node_list[r]];
+    update_xmin(x, xmin);
+    update_xmax(x, xmax);
   }
 
-  Array<R2> xmin_array = parallel::allGather(xmin);
-  Array<R2> xmax_array = parallel::allGather(xmax);
-  for (size_t i = 0; i < xmin_array.size(); ++i) {
-    const R2& x = xmin_array[i];
-    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
-      xmin = x;
+  if (parallel::size() > 1) {
+    Array<R2> xmin_array = parallel::allGather(xmin);
+    Array<R2> xmax_array = parallel::allGather(xmax);
+    for (size_t i = 0; i < xmin_array.size(); ++i) {
+      update_xmin(xmin_array[i], xmin);
     }
-  }
-  for (size_t i = 0; i < xmax_array.size(); ++i) {
-    const R2& x = xmax_array[i];
-    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
-      xmax = x;
+    for (size_t i = 0; i < xmax_array.size(); ++i) {
+      update_xmax(xmax_array[i], xmax);
     }
   }
 
   if (xmin == xmax) {
-    std::stringstream os;
-    os << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";
-    throw NormalError(os.str());
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": unable to compute normal";
+    throw NormalError(ost.str());
   }
 
   R2 dx = xmax - xmin;
   dx *= 1. / l2Norm(dx);
 
-  R2 normal(-dx[1], dx[0]);
+  const R2 normal(-dx[1], dx[0]);
 
-  this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
+  this->_checkBoundaryIsFlat(normal, 0.5 * (xmin + xmax), l2Norm(xmax - xmin), mesh);
 
   return normal;
 }
@@ -242,6 +248,54 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
   static_assert(MeshType::Dimension == 3);
   using R3 = TinyVector<3, double>;
 
+  auto update_xmin = [](const R3& x, R3& xmin) {
+    // XMIN: X.xmin X.ymax X.zmax
+    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] > xmin[1])) or
+        ((x[0] == xmin[0]) and (x[1] == xmin[1]) and (x[2] > xmin[2]))) {
+      xmin = x;
+    }
+  };
+
+  auto update_xmax = [](const R3& x, R3& xmax) {
+    // XMAX: X.xmax X.ymin X.zmin
+    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] < xmax[1])) or
+        ((x[0] == xmax[0]) and (x[1] == xmax[1]) and (x[2] < xmax[2]))) {
+      xmax = x;
+    }
+  };
+
+  auto update_ymin = [](const R3& x, R3& ymin) {
+    // YMIN: X.ymin X.zmax X.xmin
+    if ((x[1] < ymin[1]) or ((x[1] == ymin[1]) and (x[2] > ymin[2])) or
+        ((x[1] == ymin[1]) and (x[2] == ymin[2]) and (x[0] < ymin[0]))) {
+      ymin = x;
+    }
+  };
+
+  auto update_ymax = [](const R3& x, R3& ymax) {
+    // YMAX: X.ymax X.zmin X.xmax
+    if ((x[1] > ymax[1]) or ((x[1] == ymax[1]) and (x[2] < ymax[2])) or
+        ((x[1] == ymax[1]) and (x[2] == ymax[1]) and (x[0] > ymax[0]))) {
+      ymax = x;
+    }
+  };
+
+  auto update_zmin = [](const R3& x, R3& zmin) {
+    // ZMIN: X.zmin X.xmin X.ymin
+    if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[2] < zmin[2])) or
+        ((x[1] == zmin[1]) and (x[2] == zmin[1]) and (x[0] < zmin[0]))) {
+      zmin = x;
+    }
+  };
+
+  auto update_zmax = [](const R3& x, R3& zmax) {
+    // ZMAX: X.zmax X.xmax X.ymax
+    if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[2] > zmax[2])) or
+        ((x[1] == zmax[1]) and (x[2] == zmax[1]) and (x[0] > zmax[0]))) {
+      zmax = x;
+    }
+  };
+
   R3 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
   R3 ymin = xmin;
   R3 zmin = xmin;
@@ -254,66 +308,39 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
 
   for (size_t r = 0; r < m_node_list.size(); ++r) {
     const R3& x = xr[m_node_list[r]];
-    if (x[0] < xmin[0]) {
-      xmin = x;
-    }
-    if (x[1] < ymin[1]) {
-      ymin = x;
-    }
-    if (x[2] < zmin[2]) {
-      zmin = x;
-    }
-    if (x[0] > xmax[0]) {
-      xmax = x;
-    }
-    if (x[1] > ymax[1]) {
-      ymax = x;
-    }
-    if (x[2] > zmax[2]) {
-      zmax = x;
-    }
+    update_xmin(x, xmin);
+    update_xmax(x, xmax);
+    update_ymin(x, ymin);
+    update_ymax(x, ymax);
+    update_zmin(x, zmin);
+    update_zmax(x, zmax);
   }
-  Array<R3> xmin_array = parallel::allGather(xmin);
-  Array<R3> xmax_array = parallel::allGather(xmax);
-  Array<R3> ymin_array = parallel::allGather(ymin);
-  Array<R3> ymax_array = parallel::allGather(ymax);
-  Array<R3> zmin_array = parallel::allGather(zmin);
-  Array<R3> zmax_array = parallel::allGather(zmax);
-
-  for (size_t i = 0; i < xmin_array.size(); ++i) {
-    const R3& x = xmin_array[i];
-    if (x[0] < xmin[0]) {
-      xmin = x;
+
+  if (parallel::size() > 0) {
+    Array<const R3> xmin_array = parallel::allGather(xmin);
+    Array<const R3> xmax_array = parallel::allGather(xmax);
+    Array<const R3> ymin_array = parallel::allGather(ymin);
+    Array<const R3> ymax_array = parallel::allGather(ymax);
+    Array<const R3> zmin_array = parallel::allGather(zmin);
+    Array<const R3> zmax_array = parallel::allGather(zmax);
+
+    for (size_t i = 0; i < xmin_array.size(); ++i) {
+      update_xmin(xmin_array[i], xmin);
     }
-  }
-  for (size_t i = 0; i < ymin_array.size(); ++i) {
-    const R3& x = ymin_array[i];
-    if (x[1] < ymin[1]) {
-      ymin = x;
+    for (size_t i = 0; i < ymin_array.size(); ++i) {
+      update_ymin(ymin_array[i], ymin);
     }
-  }
-  for (size_t i = 0; i < zmin_array.size(); ++i) {
-    const R3& x = zmin_array[i];
-    if (x[2] < zmin[2]) {
-      zmin = x;
+    for (size_t i = 0; i < zmin_array.size(); ++i) {
+      update_zmin(zmin_array[i], zmin);
     }
-  }
-  for (size_t i = 0; i < xmax_array.size(); ++i) {
-    const R3& x = xmax_array[i];
-    if (x[0] > xmax[0]) {
-      xmax = x;
+    for (size_t i = 0; i < xmax_array.size(); ++i) {
+      update_xmax(xmax_array[i], xmax);
     }
-  }
-  for (size_t i = 0; i < ymax_array.size(); ++i) {
-    const R3& x = ymax_array[i];
-    if (x[1] > ymax[1]) {
-      ymax = x;
+    for (size_t i = 0; i < ymax_array.size(); ++i) {
+      update_ymax(ymax_array[i], ymax);
     }
-  }
-  for (size_t i = 0; i < zmax_array.size(); ++i) {
-    const R3& x = zmax_array[i];
-    if (x[2] > zmax[2]) {
-      zmax = x;
+    for (size_t i = 0; i < zmax_array.size(); ++i) {
+      update_zmax(zmax_array[i], zmax);
     }
   }
 
@@ -344,12 +371,17 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
   }
 
   if (normal_l2 == 0) {
-    throw NormalError("cannot to compute normal!");
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": unable to compute normal";
+    throw NormalError(ost.str());
   }
 
-  normal *= 1. / sqrt(normal_l2);
+  const double length = sqrt(normal_l2);
+
+  normal *= 1. / length;
 
-  // this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
+  this->_checkBoundaryIsFlat(normal, 1. / 6. * (xmin + xmax + ymin + ymax + zmin + zmax), length, mesh);
 
   return normal;
 }
-- 
GitLab