diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp index d8cea01af7c3c88a73af75cbf51c5a43f3e5ebd1..f74ba52787244cb2afa38d7c5e09cbe5397c1733 100644 --- a/src/mesh/MeshFlatNodeBoundary.cpp +++ b/src/mesh/MeshFlatNodeBoundary.cpp @@ -24,8 +24,8 @@ MeshFlatNodeBoundary<Dimension>::_checkBoundaryIsFlat(const TinyVector<Dimension if (parallel::allReduceOr(is_bad)) { std::ostringstream ost; - ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset - << ": boundary is not flat!"; + ost << "invalid boundary \"" << rang::fgB::yellow << this->m_boundary_name << rang::style::reset + << "\": boundary is not flat!"; throw NormalError(ost.str()); } } @@ -68,8 +68,8 @@ MeshFlatNodeBoundary<2>::_getNormal(const Mesh<Connectivity<2>>& mesh) if (xmin == xmax) { std::ostringstream ost; - ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset - << ": unable to compute normal"; + ost << "invalid boundary \"" << rang::fgB::yellow << this->m_boundary_name << rang::style::reset + << "\": unable to compute normal"; throw NormalError(ost.str()); } @@ -126,8 +126,8 @@ MeshFlatNodeBoundary<3>::_getNormal(const Mesh<Connectivity<3>>& mesh) if (normal_l2 == 0) { std::ostringstream ost; - ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset - << ": unable to compute normal"; + ost << "invalid boundary \"" << rang::fgB::yellow << this->m_boundary_name << rang::style::reset + << "\": unable to compute normal"; throw NormalError(ost.str()); } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index a907e37e43b49f318d45cdb5e02183341ea56015..85c99ead032dccef8ad4bb68269495ac3da88951 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -177,6 +177,7 @@ add_executable (mpi_unit_tests test_ItemValueUtils.cpp test_MeshFaceBoundary.cpp test_MeshFlatFaceBoundary.cpp + test_MeshFlatNodeBoundary.cpp test_MeshNodeBoundary.cpp test_Messenger.cpp test_OFStream.cpp diff --git a/tests/test_MeshFlatNodeBoundary.cpp b/tests/test_MeshFlatNodeBoundary.cpp new file mode 100644 index 0000000000000000000000000000000000000000..f0327319b9b256845d20472f0fa0e8a66b788eda --- /dev/null +++ b/tests/test_MeshFlatNodeBoundary.cpp @@ -0,0 +1,1471 @@ +#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/MeshFlatNodeBoundary.hpp> +#include <mesh/NamedBoundaryDescriptor.hpp> +#include <mesh/NumberedBoundaryDescriptor.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("MeshFlatNodeBoundary", "[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_node_list_from_tag = [](const size_t tag, const auto& connectivity) -> Array<const NodeId> { + for (size_t i = 0; i < connectivity.template numberOfRefItemList<ItemType::node>(); ++i) { + const auto& ref_node_list = connectivity.template refItemList<ItemType::node>(i); + const RefId ref_id = ref_node_list.refId(); + if (ref_id.tagNumber() == tag) { + return ref_node_list.list(); + } + } + return {}; + }; + + auto get_node_list_from_name = [](const std::string& name, const auto& connectivity) -> Array<const NodeId> { + for (size_t i = 0; i < connectivity.template numberOfRefItemList<ItemType::node>(); ++i) { + const auto& ref_node_list = connectivity.template refItemList<ItemType::node>(i); + const RefId ref_id = ref_node_list.refId(); + if (ref_id.tagName() == name) { + return ref_node_list.list(); + } + } + return {}; + }; + + SECTION("aligned axis") + { + SECTION("1D") + { + static constexpr size_t Dimension = 1; + + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + + using R1 = TinyVector<1>; + + 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& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R1 normal = zero; + + switch (tag) { + case 0: { + normal = R1{-1}; + break; + } + case 1: { + normal = R1{1}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R1 normal = zero; + + if (name == "XMIN") { + normal = R1{-1}; + } else if (name == "XMAX") { + normal = R1{1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + } + } + + 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& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R1 normal = zero; + + switch (tag) { + case 1: { + normal = R1{-1}; + break; + } + case 2: { + normal = R1{1}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R1 normal = zero; + + if (name == "XMIN") { + normal = R1{-1}; + } else if (name == "XMAX") { + normal = R1{1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + } + } + } + + 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& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R2 normal = zero; + + switch (tag) { + case 0: { + normal = R2{-1, 0}; + break; + } + case 1: { + normal = R2{1, 0}; + break; + } + case 2: { + normal = R2{0, -1}; + break; + } + case 3: { + normal = R2{0, 1}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(10)), + "error: invalid boundary \"XMINYMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(11)), + "error: invalid boundary \"XMAXYMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(12)), + "error: invalid boundary \"XMAXYMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(13)), + "error: invalid boundary \"XMINYMAX\": unable to compute normal"); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R2 normal = zero; + + if (name == "XMIN") { + normal = R2{-1, 0}; + } else if (name == "XMAX") { + normal = R2{1, 0}; + } else if (name == "YMIN") { + normal = R2{0, -1}; + } else if (name == "YMAX") { + normal = R2{0, 1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMIN")), + "error: invalid boundary \"XMINYMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAX")), + "error: invalid boundary \"XMINYMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMIN")), + "error: invalid boundary \"XMAXYMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAX")), + "error: invalid boundary \"XMAXYMAX\": unable to compute normal"); + } + } + } + + 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& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R2 normal = zero; + + switch (tag) { + case 1: { + normal = R2{-1, 0}; + break; + } + case 2: { + normal = R2{1, 0}; + break; + } + case 3: { + normal = R2{0, 1}; + break; + } + case 4: { + normal = R2{0, -1}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(8)), + "error: invalid boundary \"XMINYMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(9)), + "error: invalid boundary \"XMINYMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(10)), + "error: invalid boundary \"XMAXYMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(11)), + "error: invalid boundary \"XMAXYMAX\": unable to compute normal"); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + R2 normal = zero; + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + if (name == "XMIN") { + normal = R2{-1, 0}; + } else if (name == "XMAX") { + normal = R2{1, 0}; + } else if (name == "YMIN") { + normal = R2{0, -1}; + } else if (name == "YMAX") { + normal = R2{0, 1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMIN")), + "error: invalid boundary \"XMINYMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAX")), + "error: invalid boundary \"XMINYMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMIN")), + "error: invalid boundary \"XMAXYMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAX")), + "error: invalid boundary \"XMAXYMAX\": unable to compute normal"); + } + } + } + } + + 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 = {0, 1, 2, 3, 4, 5}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + switch (tag) { + case 0: { + normal = R3{-1, 0, 0}; + break; + } + case 1: { + normal = R3{1, 0, 0}; + break; + } + case 2: { + normal = R3{0, -1, 0}; + break; + } + case 3: { + normal = R3{0, 1, 0}; + break; + } + case 4: { + normal = R3{0, 0, -1}; + break; + } + case 5: { + normal = R3{0, 0, 1}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(20)), + "error: cannot find surface with name \"20\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(21)), + "error: cannot find surface with name \"21\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(22)), + "error: cannot find surface with name \"22\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(23)), + "error: cannot find surface with name \"23\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(24)), + "error: cannot find surface with name \"24\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(25)), + "error: cannot find surface with name \"25\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(26)), + "error: cannot find surface with name \"26\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(27)), + "error: cannot find surface with name \"27\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(28)), + "error: cannot find surface with name \"28\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(29)), + "error: cannot find surface with name \"29\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(30)), + "error: cannot find surface with name \"30\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(31)), + "error: cannot find surface with name \"31\""); + + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(10)), + "error: invalid boundary \"XMINYMINZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(11)), + "error: invalid boundary \"XMAXYMINZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(12)), + "error: invalid boundary \"XMAXYMAXZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(13)), + "error: invalid boundary \"XMINYMAXZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(14)), + "error: invalid boundary \"XMINYMINZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(15)), + "error: invalid boundary \"XMAXYMINZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(16)), + "error: invalid boundary \"XMAXYMAXZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(17)), + "error: invalid boundary \"XMINYMAXZMAX\": unable to compute normal"); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + if (name == "XMIN") { + normal = R3{-1, 0, 0}; + } else if (name == "XMAX") { + normal = R3{1, 0, 0}; + } else if (name == "YMIN") { + normal = R3{0, -1, 0}; + } else if (name == "YMAX") { + normal = R3{0, 1, 0}; + } else if (name == "ZMIN") { + normal = R3{0, 0, -1}; + } else if (name == "ZMAX") { + normal = R3{0, 0, 1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMIN")), + "error: cannot find surface with name \"XMINYMIN\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAX")), + "error: cannot find surface with name \"XMINYMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAX")), + "error: cannot find surface with name \"XMAXYMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAX")), + "error: cannot find surface with name \"XMINYMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINZMIN")), + "error: cannot find surface with name \"XMINZMIN\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINZMAX")), + "error: cannot find surface with name \"XMINZMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXZMAX")), + "error: cannot find surface with name \"XMAXZMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINZMAX")), + "error: cannot find surface with name \"XMINZMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("YMINZMIN")), + "error: cannot find surface with name \"YMINZMIN\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("YMINZMAX")), + "error: cannot find surface with name \"YMINZMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("YMAXZMAX")), + "error: cannot find surface with name \"YMAXZMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("YMINZMAX")), + "error: cannot find surface with name \"YMINZMAX\""); + + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMINZMIN")), + "error: invalid boundary \"XMINYMINZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMINZMIN")), + "error: invalid boundary \"XMAXYMINZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAXZMIN")), + "error: invalid boundary \"XMAXYMAXZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAXZMIN")), + "error: invalid boundary \"XMINYMAXZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMINZMAX")), + "error: invalid boundary \"XMINYMINZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMINZMAX")), + "error: invalid boundary \"XMAXYMINZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAXZMAX")), + "error: invalid boundary \"XMAXYMAXZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAXZMAX")), + "error: invalid boundary \"XMINYMAXZMAX\": unable to compute normal"); + } + } + } + + 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& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + switch (tag) { + case 22: { + normal = R3{-1, 0, 0}; + break; + } + case 23: { + normal = R3{1, 0, 0}; + break; + } + case 24: { + normal = R3{0, 0, 1}; + break; + } + case 25: { + normal = R3{0, 0, -1}; + break; + } + case 26: { + normal = R3{0, 1, 0}; + break; + } + case 27: { + normal = R3{0, -1, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(28)), + "error: cannot find surface with name \"28\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(29)), + "error: cannot find surface with name \"29\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(30)), + "error: cannot find surface with name \"30\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(31)), + "error: cannot find surface with name \"31\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(32)), + "error: cannot find surface with name \"32\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(33)), + "error: cannot find surface with name \"33\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(34)), + "error: cannot find surface with name \"34\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(35)), + "error: cannot find surface with name \"35\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(36)), + "error: cannot find surface with name \"36\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(37)), + "error: cannot find surface with name \"37\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(38)), + "error: cannot find surface with name \"38\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(39)), + "error: cannot find surface with name \"39\""); + + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(40)), + "error: invalid boundary \"XMINYMINZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(41)), + "error: invalid boundary \"XMAXYMINZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(42)), + "error: invalid boundary \"XMINYMAXZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(43)), + "error: invalid boundary \"XMINYMAXZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(44)), + "error: invalid boundary \"XMINYMINZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(45)), + "error: invalid boundary \"XMAXYMINZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(47)), + "error: invalid boundary \"XMAXYMAXZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NumberedBoundaryDescriptor(51)), + "error: invalid boundary \"XMAXYMAXZMIN\": unable to compute normal"); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + if (name == "XMIN") { + normal = R3{-1, 0, 0}; + } else if (name == "XMAX") { + normal = R3{1, 0, 0}; + } else if (name == "YMIN") { + normal = R3{0, -1, 0}; + } else if (name == "YMAX") { + normal = R3{0, 1, 0}; + } else if (name == "ZMIN") { + normal = R3{0, 0, -1}; + } else if (name == "ZMAX") { + normal = R3{0, 0, 1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMIN")), + "error: cannot find surface with name \"XMINYMIN\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAX")), + "error: cannot find surface with name \"XMINYMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAX")), + "error: cannot find surface with name \"XMAXYMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAX")), + "error: cannot find surface with name \"XMINYMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINZMIN")), + "error: cannot find surface with name \"XMINZMIN\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINZMAX")), + "error: cannot find surface with name \"XMINZMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXZMAX")), + "error: cannot find surface with name \"XMAXZMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINZMAX")), + "error: cannot find surface with name \"XMINZMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("YMINZMIN")), + "error: cannot find surface with name \"YMINZMIN\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("YMINZMAX")), + "error: cannot find surface with name \"YMINZMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("YMAXZMAX")), + "error: cannot find surface with name \"YMAXZMAX\""); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("YMINZMAX")), + "error: cannot find surface with name \"YMINZMAX\""); + + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMINZMIN")), + "error: invalid boundary \"XMINYMINZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMINZMIN")), + "error: invalid boundary \"XMAXYMINZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAXZMIN")), + "error: invalid boundary \"XMAXYMAXZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAXZMIN")), + "error: invalid boundary \"XMINYMAXZMIN\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMINZMAX")), + "error: invalid boundary \"XMINYMINZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMINZMAX")), + "error: invalid boundary \"XMAXYMINZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMAXYMAXZMAX")), + "error: invalid boundary \"XMAXYMAXZMAX\": unable to compute normal"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor("XMINYMAXZMAX")), + "error: invalid boundary \"XMINYMAXZMAX\": unable to compute normal"); + } + } + } + } + } + + 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& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R2 normal = zero; + + switch (tag) { + case 0: { + normal = R * R2{-1, 0}; + break; + } + case 1: { + normal = R * R2{1, 0}; + break; + } + case 2: { + normal = R * R2{0, -1}; + break; + } + case 3: { + normal = R * R2{0, 1}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == 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& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R2 normal = zero; + + if (name == "XMIN") { + normal = R * R2{-1, 0}; + } else if (name == "XMAX") { + normal = R * R2{1, 0}; + } else if (name == "YMIN") { + normal = R * R2{0, -1}; + } else if (name == "YMAX") { + normal = R * R2{0, 1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == 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& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R2 normal = zero; + + switch (tag) { + case 1: { + normal = R * R2{-1, 0}; + break; + } + case 2: { + normal = R * R2{1, 0}; + break; + } + case 3: { + normal = R * R2{0, 1}; + break; + } + case 4: { + normal = R * R2{0, -1}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == 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& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + R2 normal = zero; + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + if (name == "XMIN") { + normal = R * R2{-1, 0}; + } else if (name == "XMAX") { + normal = R * R2{1, 0}; + } else if (name == "YMIN") { + normal = R * R2{0, -1}; + } else if (name == "YMAX") { + normal = R * R2{0, 1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == 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 = {0, 1, 2, 3, 4, 5}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + switch (tag) { + case 0: { + normal = R * R3{-1, 0, 0}; + break; + } + case 1: { + normal = R * R3{1, 0, 0}; + break; + } + case 2: { + normal = R * R3{0, -1, 0}; + break; + } + case 3: { + normal = R * R3{0, 1, 0}; + break; + } + case 4: { + normal = R * R3{0, 0, -1}; + break; + } + case 5: { + normal = R * R3{0, 0, 1}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + if (name == "XMIN") { + normal = R * R3{-1, 0, 0}; + } else if (name == "XMAX") { + normal = R * R3{1, 0, 0}; + } else if (name == "YMIN") { + normal = R * R3{0, -1, 0}; + } else if (name == "YMAX") { + normal = R * R3{0, 1, 0}; + } else if (name == "ZMIN") { + normal = R * R3{0, 0, -1}; + } else if (name == "ZMAX") { + normal = R * R3{0, 0, 1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == 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 = {22, 23, 24, 25, 26, 27}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + switch (tag) { + case 22: { + normal = R * R3{-1, 0, 0}; + break; + } + case 23: { + normal = R * R3{1, 0, 0}; + break; + } + case 24: { + normal = R * R3{0, 0, 1}; + break; + } + case 25: { + normal = R * R3{0, 0, -1}; + break; + } + case 26: { + normal = R * R3{0, 1, 0}; + break; + } + case 27: { + normal = R * R3{0, -1, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + if (name == "XMIN") { + normal = R * R3{-1, 0, 0}; + } else if (name == "XMAX") { + normal = R * R3{1, 0, 0}; + } else if (name == "YMIN") { + normal = R * R3{0, -1, 0}; + } else if (name == "YMAX") { + normal = R * R3{0, 1, 0}; + } else if (name == "ZMIN") { + normal = R * R3{0, 0, -1}; + } else if (name == "ZMAX") { + normal = R * R3{0, 0, 1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == 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& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R2 normal = zero; + + switch (tag) { + case 1: { + normal = R2{-1, 0}; + break; + } + case 2: { + normal = R2{1, 0}; + break; + } + case 4: { + normal = R2{0, -1}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + NumberedBoundaryDescriptor numbered_boundary_descriptor(3); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor), + "error: invalid boundary \"YMAX\": boundary is not flat!"); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + R2 normal = zero; + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + if (name == "XMIN") { + normal = R2{-1, 0}; + } else if (name == "XMAX") { + normal = R2{1, 0}; + } else if (name == "YMIN") { + normal = R2{0, -1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + NamedBoundaryDescriptor named_boundary_descriptor("YMAX"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, named_boundary_descriptor), + "error: invalid boundary \"YMAX\": boundary is not flat!"); + } + } + } + } + + 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 = {0, 1, 2, 4}; + + for (auto tag : tag_set) { + auto node_list = get_node_list_from_tag(tag, connectivity); + + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + switch (tag) { + case 0: { + normal = R3{-1, 0, 0}; + break; + } + case 1: { + normal = R3{1, 0, 0}; + break; + } + case 2: { + normal = R3{0, -1, 0}; + break; + } + case 4: { + normal = R3{0, 0, -1}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + NumberedBoundaryDescriptor numbered_boundary_descriptor(3); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor), + "error: invalid boundary \"YMAX\": boundary is not flat!"); + } + + { + NumberedBoundaryDescriptor numbered_boundary_descriptor(5); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor), + "error: invalid boundary \"ZMAX\": boundary is not flat!"); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "ZMIN"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + if (name == "XMIN") { + normal = R3{-1, 0, 0}; + } else if (name == "XMAX") { + normal = R3{1, 0, 0}; + } else if (name == "YMIN") { + normal = R3{0, -1, 0}; + } else if (name == "ZMIN") { + normal = R3{0, 0, -1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + NamedBoundaryDescriptor named_boundary_descriptor("YMAX"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, named_boundary_descriptor), + "error: invalid boundary \"YMAX\": boundary is not flat!"); + } + + { + NamedBoundaryDescriptor named_boundary_descriptor("ZMAX"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, named_boundary_descriptor), + "error: invalid boundary \"ZMAX\": boundary is not flat!"); + } + } + } + + 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 = {22, 23, 25, 27}; + + for (auto tag : tag_set) { + NumberedBoundaryDescriptor numbered_boundary_descriptor(tag); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor); + + auto node_list = get_node_list_from_tag(tag, connectivity); + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + switch (tag) { + case 22: { + normal = R3{-1, 0, 0}; + break; + } + case 23: { + normal = R3{1, 0, 0}; + break; + } + case 25: { + normal = R3{0, 0, -1}; + break; + } + case 27: { + normal = R3{0, -1, 0}; + break; + } + default: { + FAIL("unexpected tag number"); + } + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + NumberedBoundaryDescriptor numbered_boundary_descriptor(24); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor), + "error: invalid boundary \"ZMAX\": boundary is not flat!"); + } + + { + NumberedBoundaryDescriptor numbered_boundary_descriptor(26); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, numbered_boundary_descriptor), + "error: invalid boundary \"YMAX\": boundary is not flat!"); + } + } + + { + const std::set<std::string> name_set = {"XMIN", "XMAX", "YMIN", "ZMIN"}; + + for (auto name : name_set) { + NamedBoundaryDescriptor named_boundary_descriptor(name); + const auto& node_boundary = getMeshFlatNodeBoundary(mesh, named_boundary_descriptor); + + auto node_list = get_node_list_from_name(name, connectivity); + + REQUIRE(is_same(node_boundary.nodeList(), node_list)); + + R3 normal = zero; + + if (name == "XMIN") { + normal = R3{-1, 0, 0}; + } else if (name == "XMAX") { + normal = R3{1, 0, 0}; + } else if (name == "YMIN") { + normal = R3{0, -1, 0}; + } else if (name == "ZMIN") { + normal = R3{0, 0, -1}; + } else { + FAIL("unexpected name: " + name); + } + + REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13)); + } + + { + NamedBoundaryDescriptor named_boundary_descriptor("YMAX"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, named_boundary_descriptor), + "error: invalid boundary \"YMAX\": boundary is not flat!"); + } + + { + NamedBoundaryDescriptor named_boundary_descriptor("ZMAX"); + REQUIRE_THROWS_WITH(getMeshFlatNodeBoundary(mesh, named_boundary_descriptor), + "error: invalid boundary \"ZMAX\": boundary is not flat!"); + } + } + } + } + } +}