Add MeshEdgeBoundary

src/mesh/CMakeLists.txt

@@ -22,6 +22,7 @@ add_library(
   MeshBuilderBase.cpp
   MeshCellZone.cpp
   MeshDataManager.cpp
+  MeshEdgeBoundary.cpp
   MeshFaceBoundary.cpp
   MeshFlatFaceBoundary.cpp
   MeshFlatNodeBoundary.cpp

src/mesh/MeshEdgeBoundary.cpp

+#include <mesh/MeshEdgeBoundary.hpp>
+
+#include <Kokkos_Vector.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Messenger.hpp>
+
+template <size_t Dimension>
+MeshEdgeBoundary<Dimension>::MeshEdgeBoundary(const Mesh<Connectivity<Dimension>>&, const RefEdgeList& ref_edge_list)
+  : m_edge_list(ref_edge_list.list()), m_boundary_name(ref_edge_list.refId().tagName())
+{}
+
+template MeshEdgeBoundary<1>::MeshEdgeBoundary(const Mesh<Connectivity<1>>&, const RefEdgeList&);
+template MeshEdgeBoundary<2>::MeshEdgeBoundary(const Mesh<Connectivity<2>>&, const RefEdgeList&);
+template MeshEdgeBoundary<3>::MeshEdgeBoundary(const Mesh<Connectivity<3>>&, const RefEdgeList&);
+
+template <size_t Dimension>
+MeshEdgeBoundary<Dimension>::MeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+                                              const RefFaceList& ref_face_list)
+  : m_boundary_name(ref_face_list.refId().tagName())
+{
+  const Array<const FaceId>& face_list = ref_face_list.list();
+  static_assert(Dimension > 1, "conversion from to edge from face is valid in dimension > 1");
+
+  if constexpr (Dimension > 2) {
+    Kokkos::vector<unsigned int> edge_ids;
+    // not enough but should reduce significantly the number of resizing
+    edge_ids.reserve(Dimension * face_list.size());
+    const auto& face_to_edge_matrix = mesh.connectivity().faceToEdgeMatrix();
+
+    for (size_t l = 0; l < face_list.size(); ++l) {
+      const FaceId face_number = face_list[l];
+      const auto& face_edges   = face_to_edge_matrix[face_number];
+
+      for (size_t e = 0; e < face_edges.size(); ++e) {
+        edge_ids.push_back(face_edges[e]);
+      }
+    }
+    std::sort(edge_ids.begin(), edge_ids.end());
+    auto last = std::unique(edge_ids.begin(), edge_ids.end());
+    edge_ids.resize(std::distance(edge_ids.begin(), last));
+
+    Array<EdgeId> edge_list(edge_ids.size());
+    parallel_for(
+      edge_ids.size(), PUGS_LAMBDA(int r) { edge_list[r] = edge_ids[r]; });
+    m_edge_list = edge_list;
+  } else if constexpr (Dimension == 2) {
+    Array<EdgeId> edge_list(face_list.size());
+    parallel_for(
+      face_list.size(), PUGS_LAMBDA(int r) { edge_list[r] = static_cast<FaceId::base_type>(face_list[r]); });
+    m_edge_list = edge_list;
+  }
+
+  // This is quite dirty but it allows a non negligible performance
+  // improvement
+  const_cast<Connectivity<Dimension>&>(mesh.connectivity())
+    .addRefItemList(RefItemList<ItemType::edge>(ref_face_list.refId(), m_edge_list, ref_face_list.isBoundary()));
+}
+
+template MeshEdgeBoundary<2>::MeshEdgeBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);
+template MeshEdgeBoundary<3>::MeshEdgeBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&);
+
+template <size_t Dimension>
+MeshEdgeBoundary<Dimension>
+getMeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+{
+  for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
+       ++i_ref_edge_list) {
+    const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
+    const RefId& ref          = ref_edge_list.refId();
+
+    if (ref == boundary_descriptor) {
+      auto edge_list = ref_edge_list.list();
+      if (not ref_edge_list.isBoundary()) {
+        std::ostringstream ost;
+        ost << "invalid boundary " << rang::fgB::yellow << boundary_descriptor << rang::style::reset
+            << ": inner edges cannot be used to define mesh boundaries";
+        throw NormalError(ost.str());
+      }
+
+      return MeshEdgeBoundary<Dimension>{mesh, ref_edge_list};
+    }
+  }
+  if constexpr (Dimension > 1) {
+    for (size_t i_ref_face_list = 0;
+         i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
+      const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
+      const RefId& ref          = ref_face_list.refId();
+
+      if (ref == boundary_descriptor) {
+        auto face_list = ref_face_list.list();
+        if (not ref_face_list.isBoundary()) {
+          std::ostringstream ost;
+          ost << "invalid boundary " << rang::fgB::yellow << boundary_descriptor << rang::style::reset
+              << ": inner edges cannot be used to define mesh boundaries";
+          throw NormalError(ost.str());
+        }
+
+        return MeshEdgeBoundary<Dimension>{mesh, ref_face_list};
+      }
+    }
+  }
+
+  std::ostringstream ost;
+  ost << "cannot find edge list with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
+
+  throw NormalError(ost.str());
+}
+
+template MeshEdgeBoundary<1> getMeshEdgeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshEdgeBoundary<2> getMeshEdgeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshEdgeBoundary<3> getMeshEdgeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);

src/mesh/MeshEdgeBoundary.hpp

+#ifndef MESH_EDGE_BOUNDARY_HPP
+#define MESH_EDGE_BOUNDARY_HPP
+
+#include <algebra/TinyVector.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/RefItemList.hpp>
+#include <utils/Array.hpp>
+
+template <size_t Dimension>
+class Connectivity;
+
+template <typename ConnectivityType>
+class Mesh;
+
+template <size_t Dimension>
+class [[nodiscard]] MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
+{
+ protected:
+  Array<const EdgeId> m_edge_list;
+  std::string m_boundary_name;
+
+  std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds(const Mesh<Connectivity<Dimension>>& mesh)
+    const;
+
+ public:
+  template <size_t MeshDimension>
+  friend MeshEdgeBoundary<MeshDimension> getMeshEdgeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
+                                                             const IBoundaryDescriptor& boundary_descriptor);
+
+  MeshEdgeBoundary& operator=(const MeshEdgeBoundary&) = default;
+  MeshEdgeBoundary& operator=(MeshEdgeBoundary&&) = default;
+
+  const Array<const EdgeId>& edgeList() const
+  {
+    return m_edge_list;
+  }
+
+ protected:
+  MeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefEdgeList& ref_edge_list);
+  MeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list);
+
+ public:
+  MeshEdgeBoundary(const MeshEdgeBoundary&) = default;   // LCOV_EXCL_LINE
+  MeshEdgeBoundary(MeshEdgeBoundary &&)     = default;   // LCOV_EXCL_LINE
+
+  MeshEdgeBoundary()          = default;
+  virtual ~MeshEdgeBoundary() = default;
+};
+
+template <size_t Dimension>
+MeshEdgeBoundary<Dimension> getMeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+                                                const IBoundaryDescriptor& boundary_descriptor);
+
+#endif   // MESH_EDGE_BOUNDARY_HPP

tests/CMakeLists.txt

@@ -175,6 +175,7 @@ add_executable (mpi_unit_tests
   test_ItemArrayUtils.cpp
   test_ItemValue.cpp
   test_ItemValueUtils.cpp
+  test_MeshEdgeBoundary.cpp
   test_MeshFaceBoundary.cpp
   test_MeshFlatFaceBoundary.cpp
   test_MeshFlatNodeBoundary.cpp

tests/test_MeshEdgeBoundary.cpp

+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshEdgeBoundary.hpp>
+#include <mesh/NamedBoundaryDescriptor.hpp>
+#include <mesh/NumberedBoundaryDescriptor.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("MeshEdgeBoundary", "[mesh]")
+{
+  auto is_same = [](const auto& a, const auto& b) -> bool {
+    if (a.size() > 0 and b.size() > 0) {
+      return (a.size() == b.size()) and (&(a[0]) == &(b[0]));
+    } else {
+      return (a.size() == b.size());
+    }
+  };
+
+  auto get_edge_list_from_tag = [](const size_t tag, const auto& connectivity) -> Array<const EdgeId> {
+    for (size_t i = 0; i < connectivity.template numberOfRefItemList<ItemType::edge>(); ++i) {
+      const auto& ref_edge_list = connectivity.template refItemList<ItemType::edge>(i);
+      const RefId ref_id        = ref_edge_list.refId();
+      if (ref_id.tagNumber() == tag) {
+        return ref_edge_list.list();
+      }
+    }
+    return {};
+  };
+
+  auto get_edge_list_from_name = [](const std::string& name, const auto& connectivity) -> Array<const EdgeId> {
+    for (size_t i = 0; i < connectivity.template numberOfRefItemList<ItemType::edge>(); ++i) {
+      const auto& ref_edge_list = connectivity.template refItemList<ItemType::edge>(i);
+      const RefId ref_id        = ref_edge_list.refId();
+      if (ref_id.tagName() == name) {
+        return ref_edge_list.list();
+      }
+    }
+    return {};
+  };
+
+  SECTION("1D")
+  {
+    static constexpr size_t Dimension = 1;
+
+    using ConnectivityType = Connectivity<Dimension>;
+    using MeshType         = Mesh<ConnectivityType>;
+
+    SECTION("cartesian 1d")
+    {
+      std::shared_ptr p_mesh = MeshDataBaseForTests::get().cartesian1DMesh();
+      const MeshType& mesh   = *p_mesh;
+
+      const ConnectivityType& connectivity = mesh.connectivity();
+
+      {
+        const std::set<size_t> tag_set = {0, 1};
+
+        for (auto tag : tag_set) {
+          NumberedBoundaryDescriptor numbered_boundary_descriptor(tag);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, numbered_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_tag(tag, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+
+      {
+        const std::set<std::string> name_set = {"XMIN", "XMAX"};
+
+        for (auto name : name_set) {
+          NamedBoundaryDescriptor named_boundary_descriptor(name);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, named_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_name(name, connectivity);
+
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+    }
+
+    SECTION("unordered 1d")
+    {
+      std::shared_ptr p_mesh = MeshDataBaseForTests::get().unordered1DMesh();
+      const MeshType& mesh   = *p_mesh;
+
+      const ConnectivityType& connectivity = mesh.connectivity();
+
+      {
+        const std::set<size_t> tag_set = {1, 2};
+
+        for (auto tag : tag_set) {
+          NumberedBoundaryDescriptor numbered_boundary_descriptor(tag);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, numbered_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_tag(tag, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+
+      {
+        const std::set<std::string> name_set = {"XMIN", "XMAX"};
+
+        for (auto name : name_set) {
+          NamedBoundaryDescriptor named_boundary_descriptor(name);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, named_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_name(name, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+    }
+  }
+
+  SECTION("2D")
+  {
+    static constexpr size_t Dimension = 2;
+
+    using ConnectivityType = Connectivity<Dimension>;
+    using MeshType         = Mesh<ConnectivityType>;
+
+    SECTION("cartesian 2d")
+    {
+      std::shared_ptr p_mesh = MeshDataBaseForTests::get().cartesian2DMesh();
+      const MeshType& mesh   = *p_mesh;
+
+      const ConnectivityType& connectivity = mesh.connectivity();
+
+      {
+        const std::set<size_t> tag_set = {0, 1, 2, 3};
+
+        for (auto tag : tag_set) {
+          NumberedBoundaryDescriptor numbered_boundary_descriptor(tag);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, numbered_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_tag(tag, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+
+      {
+        const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX"};
+
+        for (auto name : name_set) {
+          NamedBoundaryDescriptor numbered_boundary_descriptor(name);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, numbered_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_name(name, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+    }
+
+    SECTION("hybrid 2d")
+    {
+      std::shared_ptr p_mesh = MeshDataBaseForTests::get().hybrid2DMesh();
+      const MeshType& mesh   = *p_mesh;
+
+      const ConnectivityType& connectivity = mesh.connectivity();
+
+      {
+        const std::set<size_t> tag_set = {1, 2, 3, 4};
+
+        for (auto tag : tag_set) {
+          NumberedBoundaryDescriptor numbered_boundary_descriptor(tag);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, numbered_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_tag(tag, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+
+      {
+        const std::set<std::string> name_set = {"XMIN", "YMIN", "XMAX", "YMIN"};
+
+        for (auto name : name_set) {
+          NamedBoundaryDescriptor numbered_boundary_descriptor(name);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, numbered_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_name(name, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+    }
+  }
+
+  SECTION("3D")
+  {
+    static constexpr size_t Dimension = 3;
+
+    using ConnectivityType = Connectivity<Dimension>;
+    using MeshType         = Mesh<ConnectivityType>;
+
+    SECTION("cartesian 3d")
+    {
+      std::shared_ptr p_mesh = MeshDataBaseForTests::get().cartesian3DMesh();
+      const MeshType& mesh   = *p_mesh;
+
+      const ConnectivityType& connectivity = mesh.connectivity();
+
+      {
+        const std::set<size_t> tag_set = {0, 1, 2, 3, 4, 5};
+
+        for (auto tag : tag_set) {
+          NumberedBoundaryDescriptor numbered_boundary_descriptor(tag);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, numbered_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_tag(tag, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+
+      {
+        const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
+
+        for (auto name : name_set) {
+          NamedBoundaryDescriptor numbered_boundary_descriptor(name);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, numbered_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_name(name, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+    }
+
+    SECTION("hybrid 3d")
+    {
+      std::shared_ptr p_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+      const MeshType& mesh   = *p_mesh;
+
+      const ConnectivityType& connectivity = mesh.connectivity();
+
+      {
+        const std::set<size_t> tag_set = {22, 23, 24, 25, 26, 27};
+
+        for (auto tag : tag_set) {
+          NumberedBoundaryDescriptor numbered_boundary_descriptor(tag);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, numbered_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_tag(tag, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+
+      {
+        const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
+
+        for (auto name : name_set) {
+          NamedBoundaryDescriptor numbered_boundary_descriptor(name);
+          const auto& edge_boundary = getMeshEdgeBoundary(mesh, numbered_boundary_descriptor);
+
+          auto edge_list = get_edge_list_from_name(name, connectivity);
+          REQUIRE(is_same(edge_boundary.edgeList(), edge_list));
+        }
+      }
+    }
+  }
+
+  SECTION("errors")
+  {
+    SECTION("cannot find boundary")
+    {
+      static constexpr size_t Dimension = 3;
+
+      using ConnectivityType = Connectivity<Dimension>;
+      using MeshType         = Mesh<ConnectivityType>;
+
+      std::shared_ptr p_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+      const MeshType& mesh   = *p_mesh;
+
+      NamedBoundaryDescriptor named_boundary_descriptor("invalid_boundary");
+
+      REQUIRE_THROWS_WITH(getMeshEdgeBoundary(mesh, named_boundary_descriptor),
+                          "error: cannot find edge list with name \"invalid_boundary\"");
+    }
+
+    SECTION("suredge is inside")
+    {
+      SECTION("1D")
+      {
+        static constexpr size_t Dimension = 1;
+
+        using ConnectivityType = Connectivity<Dimension>;
+        using MeshType         = Mesh<ConnectivityType>;
+
+        std::shared_ptr p_mesh = MeshDataBaseForTests::get().unordered1DMesh();
+        const MeshType& mesh   = *p_mesh;
+
+        NamedBoundaryDescriptor named_boundary_descriptor("INTERFACE");
+
+        REQUIRE_THROWS_WITH(getMeshEdgeBoundary(mesh, named_boundary_descriptor),
+                            "error: invalid boundary \"INTERFACE\": inner edges cannot be used to define mesh "
+                            "boundaries");
+      }
+
+      SECTION("2D")
+      {
+        static constexpr size_t Dimension = 2;
+
+        using ConnectivityType = Connectivity<Dimension>;
+        using MeshType         = Mesh<ConnectivityType>;
+
+        std::shared_ptr p_mesh = MeshDataBaseForTests::get().hybrid2DMesh();
+        const MeshType& mesh   = *p_mesh;
+
+        NamedBoundaryDescriptor named_boundary_descriptor("INTERFACE");
+
+        REQUIRE_THROWS_WITH(getMeshEdgeBoundary(mesh, named_boundary_descriptor),
+                            "error: invalid boundary \"INTERFACE\": inner edges cannot be used to define mesh "
+                            "boundaries");
+      }
+
+      SECTION("3D")
+      {
+        static constexpr size_t Dimension = 3;
+
+        using ConnectivityType = Connectivity<Dimension>;
+        using MeshType         = Mesh<ConnectivityType>;
+
+        std::shared_ptr p_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+        const MeshType& mesh   = *p_mesh;
+
+        NamedBoundaryDescriptor named_boundary_descriptor("INTERFACE1");
+
+        REQUIRE_THROWS_WITH(getMeshEdgeBoundary(mesh, named_boundary_descriptor),
+                            "error: invalid boundary \"INTERFACE1\": inner edges cannot be used to define mesh "
+                            "boundaries");
+      }
+    }
+  }
+}