diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt index 091fb1c6c0de5e46839721ef9a3f91215fc446ab..afd3cb28a354332c35b0ca2757499c1e8da47d09 100644 --- a/src/mesh/CMakeLists.txt +++ b/src/mesh/CMakeLists.txt @@ -22,6 +22,7 @@ add_library( MeshBuilderBase.cpp MeshCellZone.cpp MeshDataManager.cpp + MeshEdgeBoundary.cpp MeshFaceBoundary.cpp MeshFlatFaceBoundary.cpp MeshFlatNodeBoundary.cpp diff --git a/src/mesh/MeshEdgeBoundary.cpp b/src/mesh/MeshEdgeBoundary.cpp new file mode 100644 index 0000000000000000000000000000000000000000..2542b0a3d6eaa1d47927e6c83147b4d06888cf06 --- /dev/null +++ b/src/mesh/MeshEdgeBoundary.cpp @@ -0,0 +1,112 @@ +#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&); diff --git a/src/mesh/MeshEdgeBoundary.hpp b/src/mesh/MeshEdgeBoundary.hpp new file mode 100644 index 0000000000000000000000000000000000000000..9ff40cc6e66070c6b1d3f6f2a5dcf3296cb9b791 --- /dev/null +++ b/src/mesh/MeshEdgeBoundary.hpp @@ -0,0 +1,54 @@ +#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 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6862833bf7d755f93ba2cd041519a6f967a91673..dfd6616c3d5f0a56b363b8e246c3e771566d5343 100644 --- a/tests/CMakeLists.txt +++ b/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 diff --git a/tests/test_MeshEdgeBoundary.cpp b/tests/test_MeshEdgeBoundary.cpp new file mode 100644 index 0000000000000000000000000000000000000000..32c1d1c6e33c2c109b5c18f67172e3f15136f0f9 --- /dev/null +++ b/tests/test_MeshEdgeBoundary.cpp @@ -0,0 +1,335 @@ +#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"); + } + } + } +}