diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ede705d518895f33e40790818a5ea2ccb7099892..3b0e5e84185f8bc0c248b8ebca9c6f1713fd1126 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -176,6 +176,7 @@ add_executable (mpi_unit_tests
test_ItemValue.cpp
test_ItemValueUtils.cpp
test_MeshFaceBoundary.cpp
+ test_MeshFlatFaceBoundary.cpp
test_Messenger.cpp
test_OFStream.cpp
test_Partitioner.cpp
diff --git a/tests/test_MeshFlatFaceBoundary.cpp b/tests/test_MeshFlatFaceBoundary.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ffd81828ed6251e30778f25ae70e03ae1bd0b174
--- /dev/null
+++ b/tests/test_MeshFlatFaceBoundary.cpp
@@ -0,0 +1,1251 @@
+#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/MeshFlatFaceBoundary.hpp>
+#include <mesh/NamedBoundaryDescriptor.hpp>
+#include <mesh/NumberedBoundaryDescriptor.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("MeshFlatFaceBoundary", "[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("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& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ REQUIRE(is_same(face_boundary.faceList(), face_list));
+
+ R1 normal = zero;
+
+ if (name == "XMIN") {
+ normal = R1{-1};
+ } else if (name == "XMAX") {
+ normal = R1{1};
+ } else {
+ FAIL("unexpected name: " + name);
+ }
+
+ REQUIRE(l2Norm(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ REQUIRE(is_same(face_boundary.faceList(), face_list));
+
+ R1 normal = zero;
+
+ if (name == "XMIN") {
+ normal = R1{-1};
+ } else if (name == "XMAX") {
+ normal = R1{1};
+ } else {
+ FAIL("unexpected name: " + name);
+ }
+
+ REQUIRE(l2Norm(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+
+ 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 = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ R2 normal = zero;
+
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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>;
+
+ 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& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+
+ 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& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+ }
+
+ 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 = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ R2 normal = zero;
+
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_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& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ NumberedBoundaryDescriptor numbered_boundary_descriptor(3);
+ REQUIRE_THROWS_WITH(getMeshFlatFaceBoundary(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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ R2 normal = zero;
+
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ NamedBoundaryDescriptor named_boundary_descriptor("YMAX");
+ REQUIRE_THROWS_WITH(getMeshFlatFaceBoundary(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 face_list = get_face_list_from_tag(tag, connectivity);
+
+ NumberedBoundaryDescriptor numbered_boundary_descriptor(tag);
+ const auto& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ NumberedBoundaryDescriptor numbered_boundary_descriptor(3);
+ REQUIRE_THROWS_WITH(getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor),
+ "error: invalid boundary YMAX: boundary is not flat!");
+ }
+
+ {
+ NumberedBoundaryDescriptor numbered_boundary_descriptor(5);
+ REQUIRE_THROWS_WITH(getMeshFlatFaceBoundary(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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ NamedBoundaryDescriptor named_boundary_descriptor("YMAX");
+ REQUIRE_THROWS_WITH(getMeshFlatFaceBoundary(mesh, named_boundary_descriptor),
+ "error: invalid boundary YMAX: boundary is not flat!");
+ }
+
+ {
+ NamedBoundaryDescriptor named_boundary_descriptor("ZMAX");
+ REQUIRE_THROWS_WITH(getMeshFlatFaceBoundary(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& face_boundary = getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor);
+
+ auto face_list = get_face_list_from_tag(tag, connectivity);
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ NumberedBoundaryDescriptor numbered_boundary_descriptor(24);
+ REQUIRE_THROWS_WITH(getMeshFlatFaceBoundary(mesh, numbered_boundary_descriptor),
+ "error: invalid boundary ZMAX: boundary is not flat!");
+ }
+
+ {
+ NumberedBoundaryDescriptor numbered_boundary_descriptor(26);
+ REQUIRE_THROWS_WITH(getMeshFlatFaceBoundary(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& face_boundary = getMeshFlatFaceBoundary(mesh, named_boundary_descriptor);
+
+ auto face_list = get_face_list_from_name(name, connectivity);
+
+ REQUIRE(is_same(face_boundary.faceList(), face_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(face_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ NamedBoundaryDescriptor named_boundary_descriptor("YMAX");
+ REQUIRE_THROWS_WITH(getMeshFlatFaceBoundary(mesh, named_boundary_descriptor),
+ "error: invalid boundary YMAX: boundary is not flat!");
+ }
+
+ {
+ NamedBoundaryDescriptor named_boundary_descriptor("ZMAX");
+ REQUIRE_THROWS_WITH(getMeshFlatFaceBoundary(mesh, named_boundary_descriptor),
+ "error: invalid boundary ZMAX: boundary is not flat!");
+ }
+ }
+ }
+ }
+ }
+}