Fix surface bounds calculation

src/mesh/MeshNodeBoundary.cpp

@@ -85,23 +85,23 @@ MeshNodeBoundary<3>::_getBounds(const Mesh<Connectivity<3>>& mesh) const
   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]))) {
+        ((x[1] == ymax[1]) and (x[2] == ymax[2]) 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]))) {
+    if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[0] < zmin[0])) or
+        ((x[2] == zmin[2]) and (x[0] == zmin[0]) and (x[1] < zmin[1]))) {
       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]))) {
+    if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[0] > zmax[0])) or
+        ((x[2] == zmax[2]) and (x[0] == zmax[0]) and (x[1] > zmax[1]))) {
       zmax = x;
     }
   };
@@ -115,31 +115,34 @@ MeshNodeBoundary<3>::_getBounds(const Mesh<Connectivity<3>>& mesh) const
   R3& zmax = bounds[5];
 
   xmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  ymin = xmin;
-  zmin = xmin;
+  ymin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  zmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
 
-  xmax = -xmin;
-  ymax = xmax;
-  zmax = xmax;
+  xmax =
+    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  ymax =
+    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  zmax =
+    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
 
   const NodeValue<const R3>& xr = mesh.xr();
 
   for (size_t r = 0; r < m_node_list.size(); ++r) {
     const R3& x = xr[m_node_list[r]];
     update_xmin(x, xmin);
-    update_xmax(x, xmax);
     update_ymin(x, ymin);
-    update_ymax(x, ymax);
     update_zmin(x, zmin);
+    update_xmax(x, xmax);
+    update_ymax(x, ymax);
     update_zmax(x, zmax);
   }
 
-  if (parallel::size() > 0) {
+  if (parallel::size() > 1) {
     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> xmax_array = parallel::allGather(xmax);
+    Array<const R3> ymax_array = parallel::allGather(ymax);
     Array<const R3> zmax_array = parallel::allGather(zmax);
 
     for (size_t i = 0; i < xmin_array.size(); ++i) {
@@ -170,19 +173,8 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>
                                               const RefFaceList& ref_face_list)
   : m_boundary_name(ref_face_list.refId().tagName())
 {
-  const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
-
   const Array<const FaceId>& face_list = ref_face_list.list();
-
-  bool is_bad = false;
-  parallel_for(face_list.size(), [=, &is_bad](int l) {
-    const auto& face_cells = face_to_cell_matrix[face_list[l]];
-    if (face_cells.size() > 1) {
-      is_bad = true;
-    }
-  });
-
-  if (parallel::allReduceOr(is_bad)) {
+  if (not ref_face_list.isBoundary()) {
     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";
@@ -217,6 +209,11 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>
       face_list.size(), PUGS_LAMBDA(int r) { node_list[r] = static_cast<FaceId::base_type>(face_list[r]); });
     m_node_list = node_list;
   }
+
+  // This is quite dirty but it allows a non negligible performance
+  // improvement
+  const_cast<Connectivity<Dimension>&>(mesh.connectivity())
+    .addRefItemList(RefItemList<ItemType::node>(ref_face_list.refId(), m_node_list, ref_face_list.isBoundary()));
 }
 
 template <size_t Dimension>
@@ -225,44 +222,7 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>
   : m_boundary_name(ref_edge_list.refId().tagName())
 {
   const Array<const EdgeId>& edge_list = ref_edge_list.list();
-  const auto& edge_is_owned            = mesh.connectivity().edgeIsOwned();
-
-  bool is_bad = false;
-
-  if constexpr ((Dimension == 1) or (Dimension == 2)) {
-    const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();
-    parallel_for(edge_list.size(), [=, &is_bad](int l) {
-      const EdgeId edge_id = edge_list[l];
-      if (edge_is_owned[edge_id]) {
-        if (edge_to_cell_matrix[edge_id].size() != 1) {
-          is_bad = true;
-        }
-      }
-    });
-  } else {
-    static_assert(Dimension == 3);
-    const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix();
-    const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
-
-    parallel_for(edge_list.size(), [=, &is_bad](int l) {
-      const EdgeId edge_id = edge_list[l];
-      if (edge_is_owned[edge_id]) {
-        const auto& edge_faces             = edge_to_face_matrix[edge_id];
-        bool is_connected_to_boundary_face = false;
-        for (size_t i_edge_face = 0; i_edge_face < edge_faces.size(); ++i_edge_face) {
-          const FaceId edge_face_id = edge_faces[i_edge_face];
-          if (face_to_cell_matrix[edge_face_id].size() == 1) {
-            is_connected_to_boundary_face = true;
-          }
-        }
-        if (not is_connected_to_boundary_face) {
-          is_bad = true;
-        }
-      }
-    });
-  }
-
-  if (parallel::allReduceOr(is_bad)) {
+  if (not ref_edge_list.isBoundary()) {
     std::ostringstream ost;
     ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
         << ": inner edges cannot be used to define mesh boundaries";
@@ -296,12 +256,24 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>
       edge_list.size(), PUGS_LAMBDA(int r) { node_list[r] = static_cast<EdgeId::base_type>(edge_list[r]); });
     m_node_list = node_list;
   }
+
+  // This is quite dirty but it allows a non negligible performance
+  // improvement
+  const_cast<Connectivity<Dimension>&>(mesh.connectivity())
+    .addRefItemList(RefItemList<ItemType::node>(ref_edge_list.refId(), m_node_list, ref_edge_list.isBoundary()));
 }
 
 template <size_t Dimension>
 MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>>&, const RefNodeList& ref_node_list)
   : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName())
-{}
+{
+  if (not ref_node_list.isBoundary()) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": inner nodes cannot be used to define mesh boundaries";
+    throw NormalError(ost.str());
+  }
+}
 
 template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefFaceList&);
 template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);