diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt index 7d7caaf31cdf3d791de1915660491939efcb32fd..091fb1c6c0de5e46839721ef9a3f91215fc446ab 100644 --- a/src/mesh/CMakeLists.txt +++ b/src/mesh/CMakeLists.txt @@ -26,6 +26,7 @@ add_library( MeshFlatFaceBoundary.cpp MeshFlatNodeBoundary.cpp MeshRelaxer.cpp + MeshLineFaceBoundary.cpp MeshLineNodeBoundary.cpp MeshNodeBoundary.cpp MeshRandomizer.cpp diff --git a/src/mesh/MeshLineFaceBoundary.cpp b/src/mesh/MeshLineFaceBoundary.cpp new file mode 100644 index 0000000000000000000000000000000000000000..546c18050296e4b224e3fadb1653a698dfcc548d --- /dev/null +++ b/src/mesh/MeshLineFaceBoundary.cpp @@ -0,0 +1,29 @@ +#include <mesh/MeshLineFaceBoundary.hpp> + +#include <mesh/Connectivity.hpp> +#include <mesh/Mesh.hpp> +#include <mesh/MeshLineNodeBoundary.hpp> +#include <utils/Messenger.hpp> + +template <size_t Dimension> +MeshLineFaceBoundary<Dimension> +getMeshLineFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor) +{ + 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) { + MeshLineNodeBoundary<Dimension> mesh_line_node_boundary = getMeshLineNodeBoundary(mesh, boundary_descriptor); + + return MeshLineFaceBoundary<Dimension>{mesh, ref_face_list, mesh_line_node_boundary.direction()}; + } + } + + std::ostringstream ost; + ost << "cannot find face list with name " << rang::fgB::red << boundary_descriptor << rang::style::reset; + + throw NormalError(ost.str()); +} + +template MeshLineFaceBoundary<2> getMeshLineFaceBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&); diff --git a/src/mesh/MeshLineFaceBoundary.hpp b/src/mesh/MeshLineFaceBoundary.hpp new file mode 100644 index 0000000000000000000000000000000000000000..5dba2d9bd1119dbbae97c8229581fc6615975b6a --- /dev/null +++ b/src/mesh/MeshLineFaceBoundary.hpp @@ -0,0 +1,49 @@ +#ifndef MESH_LINE_FACE_BOUNDARY_HPP +#define MESH_LINE_FACE_BOUNDARY_HPP + +#include <algebra/TinyMatrix.hpp> +#include <mesh/MeshFaceBoundary.hpp> + +template <size_t Dimension> +class [[nodiscard]] MeshLineFaceBoundary final + : public MeshFaceBoundary<Dimension> // clazy:exclude=copyable-polymorphic +{ + public: + static_assert(Dimension == 2, "MeshFaceBoundary makes only sense in dimension 2"); + + using Rd = TinyVector<Dimension, double>; + + private: + const Rd m_direction; + + public: + template <size_t MeshDimension> + friend MeshLineFaceBoundary<MeshDimension> getMeshLineFaceBoundary(const Mesh<Connectivity<MeshDimension>>& mesh, + const IBoundaryDescriptor& boundary_descriptor); + + PUGS_INLINE + const Rd& direction() const + { + return m_direction; + } + + MeshLineFaceBoundary& operator=(const MeshLineFaceBoundary&) = default; + MeshLineFaceBoundary& operator=(MeshLineFaceBoundary&&) = default; + + private: + MeshLineFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list, const Rd& direction) + : MeshFaceBoundary<Dimension>(mesh, ref_face_list), m_direction(direction) + {} + + public: + MeshLineFaceBoundary() = default; + MeshLineFaceBoundary(const MeshLineFaceBoundary&) = default; + MeshLineFaceBoundary(MeshLineFaceBoundary &&) = default; + ~MeshLineFaceBoundary() = default; +}; + +template <size_t Dimension> +MeshLineFaceBoundary<Dimension> getMeshLineFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, + const IBoundaryDescriptor& boundary_descriptor); + +#endif // MESH_LINE_FACE_BOUNDARY_HPP diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6db69f436ba9fcae7ae5957a4edddf14be12cf05..6862833bf7d755f93ba2cd041519a6f967a91673 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -178,6 +178,7 @@ add_executable (mpi_unit_tests test_MeshFaceBoundary.cpp test_MeshFlatFaceBoundary.cpp test_MeshFlatNodeBoundary.cpp + test_MeshLineFaceBoundary.cpp test_MeshLineNodeBoundary.cpp test_MeshNodeBoundary.cpp test_Messenger.cpp diff --git a/tests/test_MeshLineFaceBoundary.cpp b/tests/test_MeshLineFaceBoundary.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a41b9e80ac2308955cca4e71c0d6b0af2e23d60a --- /dev/null +++ b/tests/test_MeshLineFaceBoundary.cpp @@ -0,0 +1,536 @@ +#include <catch2/catch_test_macros.hpp> +#include <catch2/matchers/catch_matchers_all.hpp> + +#include <MeshDataBaseForTests.hpp> + +#include <algebra/TinyMatrix.hpp> +#include <mesh/Connectivity.hpp> +#include <mesh/Mesh.hpp> +#include <mesh/MeshLineFaceBoundary.hpp> +#include <mesh/NamedBoundaryDescriptor.hpp> +#include <mesh/NumberedBoundaryDescriptor.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("MeshLineFaceBoundary", "[mesh]") +{ + auto is_same = [](const auto& a, const auto& b) -> bool { + if (a.size() > 0 and b.size() > 0) { + return (a[0] == b[0]); + } else { + return (a.size() == b.size()); + } + }; + + auto get_face_list_from_tag = [](const size_t tag, const auto& connectivity) -> Array<const FaceId> { + for (size_t i = 0; i < connectivity.template numberOfRefItemList<ItemType::face>(); ++i) { + const auto& ref_face_list = connectivity.template refItemList<ItemType::face>(i); + const RefId ref_id = ref_face_list.refId(); + if (ref_id.tagNumber() == tag) { + return ref_face_list.list(); + } + } + return {}; + }; + + auto get_face_list_from_name = [](const std::string& name, const auto& connectivity) -> Array<const FaceId> { + for (size_t i = 0; i < connectivity.template numberOfRefItemList<ItemType::face>(); ++i) { + const auto& ref_face_list = connectivity.template refItemList<ItemType::face>(i); + const RefId ref_id = ref_face_list.refId(); + if (ref_id.tagName() == name) { + return ref_face_list.list(); + } + } + return {}; + }; + + SECTION("aligned axis") + { + SECTION("2D") + { + static constexpr size_t Dimension = 2; + + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + + using R2 = TinyVector<2>; + + 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& face_boundary = getMeshLineFaceBoundary(mesh, numbered_boundary_descriptor); + + auto face_list = get_face_list_from_tag(tag, connectivity); + REQUIRE(is_same(face_boundary.faceList(), face_list)); + + R2 direction = zero; + + switch (tag) { + case 0: { + direction = R2{0, 1}; + break; + } + case 1: { + direction = R2{0, 1}; + break; + } + case 2: { + direction = R2{1, 0}; + break; + } + case 3: { + direction = R2{1, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(face_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NumberedBoundaryDescriptor(10)), + "error: cannot find face list with name \"10\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NumberedBoundaryDescriptor(11)), + "error: cannot find face list with name \"11\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NumberedBoundaryDescriptor(12)), + "error: cannot find face list with name \"12\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NumberedBoundaryDescriptor(13)), + "error: cannot find face list with name \"13\""); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& face_boundary = getMeshLineFaceBoundary(mesh, named_boundary_descriptor); + + auto face_list = get_face_list_from_name(name, connectivity); + REQUIRE(is_same(face_boundary.faceList(), face_list)); + + R2 direction = zero; + + if (name == "XMIN") { + direction = R2{0, 1}; + } else if (name == "XMAX") { + direction = R2{0, 1}; + } else if (name == "YMIN") { + direction = R2{1, 0}; + } else if (name == "YMAX") { + direction = R2{1, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(face_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NamedBoundaryDescriptor("XMINYMIN")), + "error: cannot find face list with name \"XMINYMIN\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NamedBoundaryDescriptor("XMINYMAX")), + "error: cannot find face list with name \"XMINYMAX\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NamedBoundaryDescriptor("XMAXYMIN")), + "error: cannot find face list with name \"XMAXYMIN\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAX")), + "error: cannot find face list with name \"XMAXYMAX\""); + } + } + } + + 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& face_boundary = getMeshLineFaceBoundary(mesh, numbered_boundary_descriptor); + + auto face_list = get_face_list_from_tag(tag, connectivity); + REQUIRE(is_same(face_boundary.faceList(), face_list)); + + R2 direction = zero; + + switch (tag) { + case 1: { + direction = R2{0, 1}; + break; + } + case 2: { + direction = R2{0, 1}; + break; + } + case 3: { + direction = R2{1, 0}; + break; + } + case 4: { + direction = R2{1, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(face_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NumberedBoundaryDescriptor(8)), + "error: cannot find face list with name \"8\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NumberedBoundaryDescriptor(9)), + "error: cannot find face list with name \"9\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NumberedBoundaryDescriptor(10)), + "error: cannot find face list with name \"10\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NumberedBoundaryDescriptor(11)), + "error: cannot find face list with name \"11\""); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& face_boundary = getMeshLineFaceBoundary(mesh, named_boundary_descriptor); + + auto face_list = get_face_list_from_name(name, connectivity); + + R2 direction = zero; + + REQUIRE(is_same(face_boundary.faceList(), face_list)); + + if (name == "XMIN") { + direction = R2{0, 1}; + } else if (name == "XMAX") { + direction = R2{0, 1}; + } else if (name == "YMIN") { + direction = R2{1, 0}; + } else if (name == "YMAX") { + direction = R2{1, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(face_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NamedBoundaryDescriptor("XMINYMIN")), + "error: cannot find face list with name \"XMINYMIN\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NamedBoundaryDescriptor("XMINYMAX")), + "error: cannot find face list with name \"XMINYMAX\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NamedBoundaryDescriptor("XMAXYMIN")), + "error: cannot find face list with name \"XMAXYMIN\""); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAX")), + "error: cannot find face list with name \"XMAXYMAX\""); + } + } + } + } + } + + SECTION("rotated axis") + { + SECTION("2D") + { + static constexpr size_t Dimension = 2; + + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + + using R2 = TinyVector<2>; + + const double theta = 0.3; + const TinyMatrix<2> R{std::cos(theta), -std::sin(theta), // + std::sin(theta), std::cos(theta)}; + + SECTION("cartesian 2d") + { + std::shared_ptr p_mesh = MeshDataBaseForTests::get().cartesian2DMesh(); + + const ConnectivityType& connectivity = p_mesh->connectivity(); + + auto xr = p_mesh->xr(); + + NodeValue<R2> rotated_xr{connectivity}; + + parallel_for( + connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { rotated_xr[node_id] = R * xr[node_id]; }); + + MeshType mesh{p_mesh->shared_connectivity(), rotated_xr}; + + { + const std::set<size_t> tag_set = {0, 1, 2, 3}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& face_boundary = getMeshLineFaceBoundary(mesh, numbered_boundary_descriptor); + + auto face_list = get_face_list_from_tag(tag, connectivity); + REQUIRE(is_same(face_boundary.faceList(), face_list)); + + R2 direction = zero; + + switch (tag) { + case 0: { + direction = R * R2{0, -1}; + break; + } + case 1: { + direction = R * R2{0, -1}; + break; + } + case 2: { + direction = R * R2{1, 0}; + break; + } + case 3: { + direction = R * R2{1, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(face_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& face_boundary = getMeshLineFaceBoundary(mesh, named_boundary_descriptor); + + auto face_list = get_face_list_from_name(name, connectivity); + REQUIRE(is_same(face_boundary.faceList(), face_list)); + + R2 direction = zero; + + if (name == "XMIN") { + direction = R * R2{0, -1}; + } else if (name == "XMAX") { + direction = R * R2{0, -1}; + } else if (name == "YMIN") { + direction = R * R2{1, 0}; + } else if (name == "YMAX") { + direction = R * R2{1, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(face_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + } + } + + SECTION("hybrid 2d") + { + std::shared_ptr p_mesh = MeshDataBaseForTests::get().hybrid2DMesh(); + + const ConnectivityType& connectivity = p_mesh->connectivity(); + + auto xr = p_mesh->xr(); + + NodeValue<R2> rotated_xr{connectivity}; + parallel_for( + connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { rotated_xr[node_id] = R * xr[node_id]; }); + + MeshType mesh{p_mesh->shared_connectivity(), rotated_xr}; + + { + const std::set<size_t> tag_set = {1, 2, 3, 4}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& face_boundary = getMeshLineFaceBoundary(mesh, numbered_boundary_descriptor); + + auto face_list = get_face_list_from_tag(tag, connectivity); + REQUIRE(is_same(face_boundary.faceList(), face_list)); + + R2 direction = zero; + + switch (tag) { + case 1: { + direction = R * R2{0, -1}; + break; + } + case 2: { + direction = R * R2{0, -1}; + break; + } + case 3: { + direction = R * R2{1, 0}; + break; + } + case 4: { + direction = R * R2{1, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(face_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& face_boundary = getMeshLineFaceBoundary(mesh, named_boundary_descriptor); + + auto face_list = get_face_list_from_name(name, connectivity); + + R2 direction = zero; + + REQUIRE(is_same(face_boundary.faceList(), face_list)); + + if (name == "XMIN") { + direction = R * R2{0, -1}; + } else if (name == "XMAX") { + direction = R * R2{0, -1}; + } else if (name == "YMIN") { + direction = R * R2{1, 0}; + } else if (name == "YMAX") { + direction = R * R2{1, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(face_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + } + } + } + } + + SECTION("curved mesh") + { + SECTION("2D") + { + static constexpr size_t Dimension = 2; + + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + + using R2 = TinyVector<2>; + + auto curve = [](const R2& X) -> R2 { return R2{X[0], (1 + X[0] * X[0]) * X[1]}; }; + + SECTION("hybrid 2d") + { + std::shared_ptr p_mesh = MeshDataBaseForTests::get().hybrid2DMesh(); + + const ConnectivityType& connectivity = p_mesh->connectivity(); + + auto xr = p_mesh->xr(); + + NodeValue<TinyVector<2>> curved_xr{connectivity}; + parallel_for( + connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { curved_xr[node_id] = curve(xr[node_id]); }); + + MeshType mesh{p_mesh->shared_connectivity(), curved_xr}; + + { + const std::set<size_t> tag_set = {1, 2, 4}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& face_boundary = getMeshLineFaceBoundary(mesh, numbered_boundary_descriptor); + + auto face_list = get_face_list_from_tag(tag, connectivity); + REQUIRE(is_same(face_boundary.faceList(), face_list)); + + R2 direction = zero; + + switch (tag) { + case 1: { + direction = R2{0, 1}; + break; + } + case 2: { + direction = R2{0, 1}; + break; + } + case 4: { + direction = R2{1, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(face_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + NumberedBoundaryDescriptor numbered_boundary_descriptor(3); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, numbered_boundary_descriptor), + "error: invalid boundary \"YMAX\": boundary is not a line!"); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& face_boundary = getMeshLineFaceBoundary(mesh, named_boundary_descriptor); + + auto face_list = get_face_list_from_name(name, connectivity); + + R2 direction = zero; + + REQUIRE(is_same(face_boundary.faceList(), face_list)); + + if (name == "XMIN") { + direction = R2{0, 1}; + } else if (name == "XMAX") { + direction = R2{0, 1}; + } else if (name == "YMIN") { + direction = R2{1, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(face_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + NamedBoundaryDescriptor named_boundary_descriptor("YMAX"); + REQUIRE_THROWS_WITH(getMeshLineFaceBoundary(mesh, named_boundary_descriptor), + "error: invalid boundary \"YMAX\": boundary is not a line!"); + } + } + } + } + } +}