From f076698050244fd3174da29d537d2c4a98b803f1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 14 Apr 2022 16:47:57 +0200
Subject: [PATCH] Add tests for MeshFlatNodeBoundary

---
 src/mesh/MeshFlatNodeBoundary.cpp   |   12 +-
 tests/CMakeLists.txt                |    1 +
 tests/test_MeshFlatNodeBoundary.cpp | 1471 +++++++++++++++++++++++++++
 3 files changed, 1478 insertions(+), 6 deletions(-)
 create mode 100644 tests/test_MeshFlatNodeBoundary.cpp

diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index d8cea01af..f74ba5278 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 a907e37e4..85c99ead0 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 000000000..f0327319b
--- /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!");
+          }
+        }
+      }
+    }
+  }
+}
-- 
GitLab