Add origin calculation and storage for flat boundaries

src/mesh/MeshFlatEdgeBoundary.cpp

@@ -11,7 +11,7 @@ getMeshFlatEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundar
   MeshEdgeBoundary mesh_edge_boundary          = getMeshEdgeBoundary(mesh, boundary_descriptor);
   MeshFlatNodeBoundary mesh_flat_node_boundary = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
 
-  return MeshFlatEdgeBoundary<MeshType>{mesh, mesh_edge_boundary.refEdgeList(),
+  return MeshFlatEdgeBoundary<MeshType>{mesh, mesh_edge_boundary.refEdgeList(), mesh_flat_node_boundary.origin(),
                                         mesh_flat_node_boundary.outgoingNormal()};
 }
 

src/mesh/MeshFlatEdgeBoundary.hpp

@@ -12,9 +12,16 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=co
   using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
+  const Rd m_origin;
   const Rd m_outgoing_normal;
 
  public:
+  const Rd&
+  origin() const
+  {
+    return m_origin;
+  }
+
   const Rd&
   outgoingNormal() const
   {
@@ -29,8 +36,11 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=co
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
  private:
-  MeshFlatEdgeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list, const Rd& outgoing_normal)
-    : MeshEdgeBoundary(mesh, ref_edge_list), m_outgoing_normal(outgoing_normal)
+  MeshFlatEdgeBoundary(const MeshType& mesh,
+                       const RefEdgeList& ref_edge_list,
+                       const Rd& origin,
+                       const Rd& outgoing_normal)
+    : MeshEdgeBoundary(mesh, ref_edge_list), m_origin(origin), m_outgoing_normal(outgoing_normal)
   {}
 
  public:

src/mesh/MeshFlatFaceBoundary.cpp

@@ -11,7 +11,7 @@ getMeshFlatFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundar
   MeshFaceBoundary mesh_face_boundary          = getMeshFaceBoundary(mesh, boundary_descriptor);
   MeshFlatNodeBoundary mesh_flat_node_boundary = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
 
-  return MeshFlatFaceBoundary<MeshType>{mesh, mesh_face_boundary.refFaceList(),
+  return MeshFlatFaceBoundary<MeshType>{mesh, mesh_face_boundary.refFaceList(), mesh_flat_node_boundary.origin(),
                                         mesh_flat_node_boundary.outgoingNormal()};
 }
 

src/mesh/MeshFlatFaceBoundary.hpp

@@ -12,9 +12,16 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=co
   using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
+  const Rd m_origin;
   const Rd m_outgoing_normal;
 
  public:
+  const Rd&
+  origin() const
+  {
+    return m_origin;
+  }
+
   const Rd&
   outgoingNormal() const
   {
@@ -29,8 +36,11 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=co
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
  private:
-  MeshFlatFaceBoundary(const MeshType& mesh, const RefFaceList& ref_face_list, const Rd& outgoing_normal)
-    : MeshFaceBoundary(mesh, ref_face_list), m_outgoing_normal(outgoing_normal)
+  MeshFlatFaceBoundary(const MeshType& mesh,
+                       const RefFaceList& ref_face_list,
+                       const Rd& origin,
+                       const Rd& outgoing_normal)
+    : MeshFaceBoundary(mesh, ref_face_list), m_origin(origin), m_outgoing_normal(outgoing_normal)
   {}
 
  public:

src/mesh/MeshFlatNodeBoundary.cpp

@@ -34,14 +34,13 @@ MeshFlatNodeBoundary<MeshType>::_checkBoundaryIsFlat(const TinyVector<MeshType::
 
 template <>
 TinyVector<1, double>
-MeshFlatNodeBoundary<Mesh<1>>::_getNormal(const Mesh<1>& mesh)
+MeshFlatNodeBoundary<Mesh<1>>::_getOrigin(const Mesh<1>& mesh)
 {
-  using R = TinyVector<1, double>;
+  auto node_is_owned = mesh.connectivity().nodeIsOwned();
+  auto node_list     = m_ref_node_list.list();
 
   const size_t number_of_bc_nodes = [&]() {
     size_t nb_bc_nodes = 0;
-    auto node_is_owned = mesh.connectivity().nodeIsOwned();
-    auto node_list     = m_ref_node_list.list();
     for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
       nb_bc_nodes += (node_is_owned[node_list[i_node]]);
     }
@@ -55,7 +54,43 @@ MeshFlatNodeBoundary<Mesh<1>>::_getNormal(const Mesh<1>& mesh)
     throw NormalError(ost.str());
   }
 
-  return R{1};
+  auto xr = mesh.xr();
+
+  Array<Rd> origin;
+  for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+    const NodeId node_id = node_list[i_node];
+    if (node_is_owned[node_list[i_node]]) {
+      origin    = Array<Rd>(1);
+      origin[0] = xr[node_id];
+    }
+  }
+  origin = parallel::allGatherVariable(origin);
+  Assert(origin.size() == 1);
+  return origin[0];
+}
+
+template <>
+TinyVector<1, double>
+MeshFlatNodeBoundary<Mesh<1>>::_getNormal(const Mesh<1>&)
+{
+  using R1 = TinyVector<1, double>;
+
+  // The verification of the unicity of the boundary node is performed by _getOrigin()
+  return R1{1};
+}
+
+template <>
+TinyVector<2, double>
+MeshFlatNodeBoundary<Mesh<2>>::_getOrigin(const Mesh<2>& mesh)
+{
+  using R2 = TinyVector<2, double>;
+
+  std::array<R2, 2> bounds = getBounds(mesh, m_ref_node_list);
+
+  const R2& xmin = bounds[0];
+  const R2& xmax = bounds[1];
+
+  return 0.5 * (xmin + xmax);
 }
 
 template <>
@@ -137,6 +172,46 @@ MeshFlatNodeBoundary<Mesh<3>>::_getFarestNode(const Mesh<3>& mesh, const Rd& x0,
   return farest_x;
 }
 
+template <>
+TinyVector<3, double>
+MeshFlatNodeBoundary<Mesh<3>>::_getOrigin(const Mesh<3>& mesh)
+{
+  using R3 = TinyVector<3, double>;
+
+  std::array<R3, 2> diagonal = [](const std::array<R3, 6>& bounds) {
+    size_t max_i      = 0;
+    size_t max_j      = 0;
+    double max_length = 0;
+
+    for (size_t i = 0; i < bounds.size(); ++i) {
+      for (size_t j = i + 1; j < bounds.size(); ++j) {
+        double length = l2Norm(bounds[i] - bounds[j]);
+        if (length > max_length) {
+          max_i      = i;
+          max_j      = j;
+          max_length = length;
+        }
+      }
+    }
+
+    return std::array<R3, 2>{bounds[max_i], bounds[max_j]};
+  }(getBounds(mesh, m_ref_node_list));
+
+  const R3& x0 = diagonal[0];
+  const R3& x1 = diagonal[1];
+
+  if (x0 == x1) {
+    std::ostringstream ost;
+    ost << "invalid boundary \"" << rang::fgB::yellow << m_ref_node_list.refId() << rang::style::reset
+        << "\": unable to compute normal";
+    throw NormalError(ost.str());
+  }
+
+  const R3 x2 = this->_getFarestNode(mesh, x0, x1);
+
+  return 1. / 3. * (x0 + x1 + x2);
+}
+
 template <>
 TinyVector<3, double>
 MeshFlatNodeBoundary<Mesh<3>>::_getNormal(const Mesh<3>& mesh)

src/mesh/MeshFlatNodeBoundary.hpp

@@ -13,6 +13,7 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
   using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
+  const Rd m_origin;
   const Rd m_outgoing_normal;
 
   Rd _getFarestNode(const MeshType& mesh, const Rd& x0, const Rd& x1);
@@ -24,9 +25,17 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
                             const double length,
                             const MeshType& mesh) const;
 
+  Rd _getOrigin(const MeshType& mesh);
+
   Rd _getOutgoingNormal(const MeshType& mesh);
 
  public:
+  const Rd&
+  origin() const
+  {
+    return m_origin;
+  }
+
   const Rd&
   outgoingNormal() const
   {
@@ -42,11 +51,11 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
 
  private:
   MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
-    : MeshNodeBoundary(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
+    : MeshNodeBoundary(mesh, ref_face_list), m_origin{_getOrigin(mesh)}, m_outgoing_normal(_getOutgoingNormal(mesh))
   {}
 
   MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
-    : MeshNodeBoundary(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
+    : MeshNodeBoundary(mesh, ref_node_list), m_origin{_getOrigin(mesh)}, m_outgoing_normal(_getOutgoingNormal(mesh))
   {}
 
  public: