diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt index 113e2527dbe4c3896979c3bd7f68a00369e62029..66ce42ae6f9e490a0b3e6e406934478b0f026f4a 100644 --- a/src/mesh/CMakeLists.txt +++ b/src/mesh/CMakeLists.txt @@ -28,6 +28,7 @@ add_library( MeshFlatFaceBoundary.cpp MeshFlatNodeBoundary.cpp MeshRelaxer.cpp + MeshLineEdgeBoundary.cpp MeshLineFaceBoundary.cpp MeshLineNodeBoundary.cpp MeshNodeBoundary.cpp diff --git a/src/mesh/MeshLineEdgeBoundary.cpp b/src/mesh/MeshLineEdgeBoundary.cpp new file mode 100644 index 0000000000000000000000000000000000000000..16e2e3e9764d1465cd767c130fbf27b0e7589afd --- /dev/null +++ b/src/mesh/MeshLineEdgeBoundary.cpp @@ -0,0 +1,19 @@ +#include <mesh/MeshLineEdgeBoundary.hpp> + +#include <mesh/Connectivity.hpp> +#include <mesh/Mesh.hpp> +#include <mesh/MeshLineNodeBoundary.hpp> +#include <utils/Messenger.hpp> + +template <size_t Dimension> +MeshLineEdgeBoundary<Dimension> +getMeshLineEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor) +{ + MeshEdgeBoundary<Dimension> mesh_edge_boundary = getMeshEdgeBoundary(mesh, boundary_descriptor); + MeshLineNodeBoundary<Dimension> mesh_line_node_boundary = getMeshLineNodeBoundary(mesh, boundary_descriptor); + + return MeshLineEdgeBoundary<Dimension>{mesh, mesh_edge_boundary.refEdgeList(), mesh_line_node_boundary.direction()}; +} + +template MeshLineEdgeBoundary<2> getMeshLineEdgeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&); +template MeshLineEdgeBoundary<3> getMeshLineEdgeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&); diff --git a/src/mesh/MeshLineEdgeBoundary.hpp b/src/mesh/MeshLineEdgeBoundary.hpp new file mode 100644 index 0000000000000000000000000000000000000000..160444648524f266941a8d2903124bf91b77ee8b --- /dev/null +++ b/src/mesh/MeshLineEdgeBoundary.hpp @@ -0,0 +1,49 @@ +#ifndef MESH_LINE_EDGE_BOUNDARY_HPP +#define MESH_LINE_EDGE_BOUNDARY_HPP + +#include <algebra/TinyMatrix.hpp> +#include <mesh/MeshEdgeBoundary.hpp> + +template <size_t Dimension> +class [[nodiscard]] MeshLineEdgeBoundary final + : public MeshEdgeBoundary<Dimension> // clazy:exclude=copyable-polymorphic +{ + public: + static_assert(Dimension > 1, "MeshLineEdgeBoundary makes only sense in dimension 1"); + + using Rd = TinyVector<Dimension, double>; + + private: + const Rd m_direction; + + public: + template <size_t MeshDimension> + friend MeshLineEdgeBoundary<MeshDimension> getMeshLineEdgeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh, + const IBoundaryDescriptor& boundary_descriptor); + + PUGS_INLINE + const Rd& direction() const + { + return m_direction; + } + + MeshLineEdgeBoundary& operator=(const MeshLineEdgeBoundary&) = default; + MeshLineEdgeBoundary& operator=(MeshLineEdgeBoundary&&) = default; + + private: + MeshLineEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefEdgeList& ref_edge_list, const Rd& direction) + : MeshEdgeBoundary<Dimension>(mesh, ref_edge_list), m_direction(direction) + {} + + public: + MeshLineEdgeBoundary() = default; + MeshLineEdgeBoundary(const MeshLineEdgeBoundary&) = default; + MeshLineEdgeBoundary(MeshLineEdgeBoundary &&) = default; + ~MeshLineEdgeBoundary() = default; +}; + +template <size_t Dimension> +MeshLineEdgeBoundary<Dimension> getMeshLineEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, + const IBoundaryDescriptor& boundary_descriptor); + +#endif // MESH_LINE_EDGE_BOUNDARY_HPP diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 0b2849df2e0da1394b80800c150ded28b4cf1ada..196938400cab3735d0633b1d1ba41efc00358f4b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -180,6 +180,7 @@ add_executable (mpi_unit_tests test_MeshFlatEdgeBoundary.cpp test_MeshFlatFaceBoundary.cpp test_MeshFlatNodeBoundary.cpp + test_MeshLineEdgeBoundary.cpp test_MeshLineFaceBoundary.cpp test_MeshLineNodeBoundary.cpp test_MeshNodeBoundary.cpp diff --git a/tests/test_MeshLineEdgeBoundary.cpp b/tests/test_MeshLineEdgeBoundary.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b4b21d0dc7e692216f80d0119f97d423518eb884 --- /dev/null +++ b/tests/test_MeshLineEdgeBoundary.cpp @@ -0,0 +1,1258 @@ +#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/MeshLineEdgeBoundary.hpp> +#include <mesh/NamedBoundaryDescriptor.hpp> +#include <mesh/NumberedBoundaryDescriptor.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("MeshLineEdgeBoundary", "[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_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("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& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + auto edge_list = get_edge_list_from_tag(tag, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_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(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(10)), + "error: cannot find edge list with name \"10\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(11)), + "error: cannot find edge list with name \"11\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(12)), + "error: cannot find edge list with name \"12\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(13)), + "error: cannot find edge 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& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_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(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMIN")), + "error: cannot find edge list with name \"XMINYMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAX")), + "error: cannot find edge list with name \"XMINYMAX\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMIN")), + "error: cannot find edge list with name \"XMAXYMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAX")), + "error: cannot find edge 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& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + auto edge_list = get_edge_list_from_tag(tag, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_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(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(8)), + "error: cannot find edge list with name \"8\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(9)), + "error: cannot find edge list with name \"9\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(10)), + "error: cannot find edge list with name \"10\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(11)), + "error: cannot find edge 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& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + + R2 direction = zero; + + REQUIRE(is_same(edge_boundary.edgeList(), edge_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(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMIN")), + "error: cannot find edge list with name \"XMINYMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAX")), + "error: cannot find edge list with name \"XMINYMAX\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMIN")), + "error: cannot find edge list with name \"XMAXYMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAX")), + "error: cannot find edge list with name \"XMAXYMAX\""); + } + } + } + } + + SECTION("3D") + { + static constexpr size_t Dimension = 3; + + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + + using R3 = TinyVector<3>; + + 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 = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + auto edge_list = get_edge_list_from_tag(tag, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + switch (tag) { + case 20: + case 21: + case 22: + case 23: { + direction = R3{0, 0, -1}; + break; + } + case 24: + case 25: + case 26: + case 27: { + direction = R3{0, -1, 0}; + break; + } + case 28: + case 29: + case 30: + case 31: { + direction = R3{1, 0, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(0)), + "error: invalid boundary \"XMIN(0)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(1)), + "error: invalid boundary \"XMAX(1)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(2)), + "error: invalid boundary \"YMIN(2)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(3)), + "error: invalid boundary \"YMAX(3)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(4)), + "error: invalid boundary \"ZMIN(4)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(5)), + "error: invalid boundary \"ZMAX(5)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(10)), + "error: cannot find edge list with name \"10\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(11)), + "error: cannot find edge list with name \"11\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(12)), + "error: cannot find edge list with name \"12\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(13)), + "error: cannot find edge list with name \"13\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(14)), + "error: cannot find edge list with name \"14\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(15)), + "error: cannot find edge list with name \"15\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(16)), + "error: cannot find edge list with name \"16\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(17)), + "error: cannot find edge list with name \"17\""); + } + + { + const std::set<std::string> name_set = {"XMINYMIN", "XMINYMAX", "XMAXYMIN", "XMAXYMAX", + "XMINZMIN", "XMINZMAX", "XMAXZMAX", "XMINZMAX", + "YMINZMIN", "YMINZMAX", "YMAXZMAX", "YMAXZMIN"}; + + for (const auto& name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + if ((name == "XMINYMIN") or (name == "XMINYMAX") or (name == "XMAXYMIN") or (name == "XMAXYMAX")) { + direction = R3{0, 0, -1}; + } else if ((name == "XMINZMIN") or (name == "XMINZMAX") or (name == "XMAXZMIN") or (name == "XMAXZMAX")) { + direction = R3{0, -1, 0}; + } else if ((name == "YMINZMIN") or (name == "YMINZMAX") or (name == "YMAXZMIN") or (name == "YMAXZMAX")) { + direction = R3{1, 0, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMIN")), + "error: invalid boundary \"XMIN(0)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAX")), + "error: invalid boundary \"XMAX(1)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("YMIN")), + "error: invalid boundary \"YMIN(2)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("YMAX")), + "error: invalid boundary \"YMAX(3)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("ZMIN")), + "error: invalid boundary \"ZMIN(4)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("ZMAX")), + "error: invalid boundary \"ZMAX(5)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMINZMIN")), + "error: cannot find edge list with name \"XMINYMINZMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMINZMIN")), + "error: cannot find edge list with name \"XMAXYMINZMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAXZMIN")), + "error: cannot find edge list with name \"XMAXYMAXZMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAXZMIN")), + "error: cannot find edge list with name \"XMINYMAXZMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMINZMAX")), + "error: cannot find edge list with name \"XMINYMINZMAX\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMINZMAX")), + "error: cannot find edge list with name \"XMAXYMINZMAX\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAXZMAX")), + "error: cannot find edge list with name \"XMAXYMAXZMAX\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAXZMAX")), + "error: cannot find edge list with name \"XMINYMAXZMAX\""); + } + } + } + + 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 = {28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + auto edge_list = get_edge_list_from_tag(tag, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + switch (tag) { + case 30: + case 31: + case 34: + case 35: { + direction = R3{0, 0, -1}; + break; + } + case 28: + case 29: + case 32: + case 33: { + direction = R3{0, -1, 0}; + break; + } + case 36: + case 37: + case 38: + case 39: { + direction = R3{1, 0, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(22)), + "error: invalid boundary \"XMIN(22)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(23)), + "error: invalid boundary \"XMAX(23)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(24)), + "error: invalid boundary \"ZMAX(24)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(25)), + "error: invalid boundary \"ZMIN(25)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(26)), + "error: invalid boundary \"YMAX(26)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(27)), + "error: invalid boundary \"YMIN(27)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(40)), + "error: cannot find edge list with name \"40\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(41)), + "error: cannot find edge list with name \"41\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(42)), + "error: cannot find edge list with name \"42\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(43)), + "error: cannot find edge list with name \"43\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(44)), + "error: cannot find edge list with name \"44\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(45)), + "error: cannot find edge list with name \"45\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(47)), + "error: cannot find edge list with name \"47\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(51)), + "error: cannot find edge list with name \"51\""); + } + } + + { + const std::set<std::string> name_set = {"XMINYMIN", "XMINYMAX", "XMAXYMIN", "XMAXYMAX", + "XMINZMIN", "XMINZMAX", "XMAXZMAX", "XMINZMAX", + "YMINZMIN", "YMINZMAX", "YMAXZMAX", "YMAXZMIN"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + if ((name == "XMINYMIN") or (name == "XMINYMAX") or (name == "XMAXYMIN") or (name == "XMAXYMAX")) { + direction = R3{0, 0, -1}; + } else if ((name == "XMINZMIN") or (name == "XMINZMAX") or (name == "XMAXZMIN") or (name == "XMAXZMAX")) { + direction = R3{0, -1, 0}; + } else if ((name == "YMINZMIN") or (name == "YMINZMAX") or (name == "YMAXZMIN") or (name == "YMAXZMAX")) { + direction = R3{1, 0, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMIN")), + "error: invalid boundary \"XMIN(22)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAX")), + "error: invalid boundary \"XMAX(23)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("YMIN")), + "error: invalid boundary \"YMIN(27)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("YMAX")), + "error: invalid boundary \"YMAX(26)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("ZMIN")), + "error: invalid boundary \"ZMIN(25)\": boundary is not a line!"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("ZMAX")), + "error: invalid boundary \"ZMAX(24)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMINZMIN")), + "error: cannot find edge list with name \"XMINYMINZMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMINZMIN")), + "error: cannot find edge list with name \"XMAXYMINZMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAXZMIN")), + "error: cannot find edge list with name \"XMAXYMAXZMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAXZMIN")), + "error: cannot find edge list with name \"XMINYMAXZMIN\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMINZMAX")), + "error: cannot find edge list with name \"XMINYMINZMAX\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMINZMAX")), + "error: cannot find edge list with name \"XMAXYMINZMAX\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAXZMAX")), + "error: cannot find edge list with name \"XMAXYMAXZMAX\""); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAXZMAX")), + "error: cannot find edge list with name \"XMINYMAXZMAX\""); + } + } + } + } + } + + 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& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + auto edge_list = get_edge_list_from_tag(tag, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_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(edge_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& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_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(edge_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& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + auto edge_list = get_edge_list_from_tag(tag, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_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(edge_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& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + + R2 direction = zero; + + REQUIRE(is_same(edge_boundary.edgeList(), edge_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(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + } + } + } + + SECTION("3D") + { + static constexpr size_t Dimension = 3; + + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + + using R3 = TinyVector<3>; + + const double theta = 0.3; + const double phi = 0.4; + const TinyMatrix<3> R = + TinyMatrix<3>{std::cos(theta), -std::sin(theta), 0, std::sin(theta), std::cos(theta), 0, 0, 0, 1} * + TinyMatrix<3>{0, std::cos(phi), -std::sin(phi), 0, std::sin(phi), std::cos(phi), 1, 0, 0}; + + SECTION("cartesian 3d") + { + std::shared_ptr p_mesh = MeshDataBaseForTests::get().cartesian3DMesh(); + + const ConnectivityType& connectivity = p_mesh->connectivity(); + + auto xr = p_mesh->xr(); + + NodeValue<R3> 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 = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + auto edge_list = get_edge_list_from_tag(tag, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + switch (tag) { + case 20: + case 21: + case 22: + case 23: { + direction = R * R3{0, 0, -1}; + break; + } + case 24: + case 25: + case 26: + case 27: { + direction = R * R3{0, 1, 0}; + break; + } + case 28: + case 29: + case 30: + case 31: { + direction = R * R3{-1, 0, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + } + + { + const std::set<std::string> name_set = {"XMINYMIN", "XMINYMAX", "XMAXYMIN", "XMAXYMAX", + "XMINZMIN", "XMINZMAX", "XMAXZMAX", "XMINZMAX", + "YMINZMIN", "YMINZMAX", "YMAXZMAX", "YMAXZMIN"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + if ((name == "XMINYMIN") or (name == "XMINYMAX") or (name == "XMAXYMIN") or (name == "XMAXYMAX")) { + direction = R * R3{0, 0, -1}; + } else if ((name == "XMINZMIN") or (name == "XMINZMAX") or (name == "XMAXZMIN") or (name == "XMAXZMAX")) { + direction = R * R3{0, 1, 0}; + } else if ((name == "YMINZMIN") or (name == "YMINZMAX") or (name == "YMAXZMIN") or (name == "YMAXZMAX")) { + direction = R * R3{-1, 0, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + } + } + + SECTION("hybrid 3d") + { + std::shared_ptr p_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + + const ConnectivityType& connectivity = p_mesh->connectivity(); + + auto xr = p_mesh->xr(); + + NodeValue<R3> 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 = {28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + auto edge_list = get_edge_list_from_tag(tag, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + switch (tag) { + case 30: + case 31: + case 34: + case 35: { + direction = R * R3{0, 0, -1}; + break; + } + case 28: + case 29: + case 32: + case 33: { + direction = R * R3{0, 1, 0}; + break; + } + case 36: + case 37: + case 38: + case 39: { + direction = R * R3{-1, 0, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + } + + { + const std::set<std::string> name_set = {"XMINYMIN", "XMINYMAX", "XMAXYMIN", "XMAXYMAX", + "XMINZMIN", "XMINZMAX", "XMAXZMAX", "XMINZMAX", + "YMINZMIN", "YMINZMAX", "YMAXZMAX", "YMAXZMIN"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + if ((name == "XMINYMIN") or (name == "XMINYMAX") or (name == "XMAXYMIN") or (name == "XMAXYMAX")) { + direction = R * R3{0, 0, -1}; + } else if ((name == "XMINZMIN") or (name == "XMINZMAX") or (name == "XMAXZMIN") or (name == "XMAXZMAX")) { + direction = R * R3{0, 1, 0}; + } else if ((name == "YMINZMIN") or (name == "YMINZMAX") or (name == "YMAXZMIN") or (name == "YMAXZMAX")) { + direction = R * R3{-1, 0, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(edge_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& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + auto edge_list = get_edge_list_from_tag(tag, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_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(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + NumberedBoundaryDescriptor numbered_boundary_descriptor(3); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor), + "error: invalid boundary \"YMAX(3)\": 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& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + + R2 direction = zero; + + REQUIRE(is_same(edge_boundary.edgeList(), edge_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(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + NamedBoundaryDescriptor named_boundary_descriptor("YMAX"); + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, named_boundary_descriptor), + "error: invalid boundary \"YMAX(3)\": boundary is not a line!"); + } + } + } + } + + SECTION("3D") + { + static constexpr size_t Dimension = 3; + + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + + using R3 = TinyVector<3>; + + SECTION("cartesian 3d") + { + std::shared_ptr p_mesh = MeshDataBaseForTests::get().cartesian3DMesh(); + + const ConnectivityType& connectivity = p_mesh->connectivity(); + + auto xr = p_mesh->xr(); + + auto curve = [](const R3& X) -> R3 { + return R3{X[0], (1 + X[0] * X[0]) * (X[1] + 1), (1 - 0.2 * X[0] * X[0]) * X[2]}; + }; + + NodeValue<R3> 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 = {20, 21, 22, 23, 24, 25, 26, 27, 28}; + + for (auto tag : tag_set) { + auto edge_list = get_edge_list_from_tag(tag, connectivity); + + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + switch (tag) { + case 20: + case 21: + case 22: + case 23: { + direction = R3{0, 0, -1}; + break; + } + case 24: + case 25: + case 26: + case 27: { + direction = R3{0, -1, 0}; + break; + } + case 28: { + direction = R3{1, 0, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(29)), + "error: invalid boundary \"YMINZMAX(29)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(30)), + "error: invalid boundary \"YMAXZMAX(30)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(31)), + "error: invalid boundary \"YMAXZMIN(31)\": boundary is not a line!"); + } + } + + { + const std::set<std::string> name_set = {"XMINYMIN", "XMINYMAX", "XMAXYMIN", "XMAXYMAX", "XMINZMIN", + "XMINZMAX", "XMAXZMAX", "XMINZMAX", "YMINZMIN"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + if ((name == "XMINYMIN") or (name == "XMINYMAX") or (name == "XMAXYMIN") or (name == "XMAXYMAX")) { + direction = R3{0, 0, -1}; + } else if ((name == "XMINZMIN") or (name == "XMINZMAX") or (name == "XMAXZMIN") or (name == "XMAXZMAX")) { + direction = R3{0, -1, 0}; + } else if ((name == "YMINZMIN")) { + direction = R3{1, 0, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("YMAXZMAX")), + "error: invalid boundary \"YMAXZMAX(30)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("YMINZMAX")), + "error: invalid boundary \"YMINZMAX(29)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("YMAXZMIN")), + "error: invalid boundary \"YMAXZMIN(31)\": boundary is not a line!"); + } + } + } + + SECTION("hybrid 3d") + { + std::shared_ptr p_mesh = MeshDataBaseForTests::get().hybrid3DMesh(); + + const ConnectivityType& connectivity = p_mesh->connectivity(); + + auto xr = p_mesh->xr(); + + auto curve = [](const R3& X) -> R3 { + return R3{X[0], (1 + X[0] * X[0]) * X[1], (1 - 0.2 * X[0] * X[0]) * X[2]}; + }; + + NodeValue<R3> 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 = {28, 29, 30, 31, 32, 33, 34, 35, 36}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, numbered_boundary_descriptor); + + auto edge_list = get_edge_list_from_tag(tag, connectivity); + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + switch (tag) { + case 30: + case 31: + case 34: + case 35: { + direction = R3{0, 0, -1}; + break; + } + case 28: + case 29: + case 32: + case 33: { + direction = R3{0, -1, 0}; + break; + } + case 36: + // case 37: + // case 38: + case 39: { + direction = R3{1, 0, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(37)), + "error: invalid boundary \"YMINZMAX(37)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(38)), + "error: invalid boundary \"YMAXZMIN(38)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NumberedBoundaryDescriptor(39)), + "error: invalid boundary \"YMAXZMAX(39)\": boundary is not a line!"); + } + } + + { + const std::set<std::string> name_set = {"XMINYMIN", "XMINYMAX", "XMAXYMIN", "XMAXYMAX", "XMINZMIN", + "XMINZMAX", "XMAXZMAX", "XMINZMAX", "YMINZMIN"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& edge_boundary = getMeshLineEdgeBoundary(mesh, named_boundary_descriptor); + + auto edge_list = get_edge_list_from_name(name, connectivity); + + REQUIRE(is_same(edge_boundary.edgeList(), edge_list)); + + R3 direction = zero; + + if ((name == "XMINYMIN") or (name == "XMINYMAX") or (name == "XMAXYMIN") or (name == "XMAXYMAX")) { + direction = R3{0, 0, -1}; + } else if ((name == "XMINZMIN") or (name == "XMINZMAX") or (name == "XMAXZMIN") or (name == "XMAXZMAX")) { + direction = R3{0, -1, 0}; + } else if ((name == "YMINZMIN")) { + direction = R3{1, 0, 0}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(edge_boundary.direction() - direction) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("YMAXZMAX")), + "error: invalid boundary \"YMAXZMAX(39)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("YMINZMAX")), + "error: invalid boundary \"YMINZMAX(37)\": boundary is not a line!"); + + REQUIRE_THROWS_WITH(getMeshLineEdgeBoundary(mesh, NamedBoundaryDescriptor("YMAXZMIN")), + "error: invalid boundary \"YMAXZMIN(38)\": boundary is not a line!"); + } + } + } + } + } +}