Add tests for situation that enlightened normal calculation issue

tests/test_MeshFlatNodeBoundary.cpp

@@ -1133,6 +1133,211 @@ TEST_CASE("MeshFlatNodeBoundary", "[mesh]")
     }
   }
 
+  SECTION("rotated diamond")
+  {
+    SECTION("2D")
+    {
+      static constexpr size_t Dimension = 2;
+
+      using ConnectivityType = Connectivity<Dimension>;
+      using MeshType         = Mesh<ConnectivityType>;
+
+      using R2 = TinyVector<2>;
+
+      auto T = [](const R2& x) -> R2 { return R2{x[0] + 0.1 * x[1], x[1] + 0.1 * x[0]}; };
+
+      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] = T(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 = 1. / std::sqrt(1.01) * R2{-1, 0.1};
+              break;
+            }
+            case 1: {
+              normal = 1. / std::sqrt(1.01) * R2{1, -0.1};
+              break;
+            }
+            case 2: {
+              normal = 1. / std::sqrt(1.01) * R2{0.1, -1};
+              break;
+            }
+            case 3: {
+              normal = 1. / std::sqrt(1.01) * R2{-0.1, 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 (const 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 = 1. / std::sqrt(1.01) * R2{-1, 0.1};
+            } else if (name == "XMAX") {
+              normal = 1. / std::sqrt(1.01) * R2{1, -0.1};
+            } else if (name == "YMIN") {
+              normal = 1. / std::sqrt(1.01) * R2{0.1, -1};
+            } else if (name == "YMAX") {
+              normal = 1. / std::sqrt(1.01) * R2{-0.1, 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>;
+
+      auto T = [](const R3& x) -> R3 {
+        return R3{x[0] + 0.1 * x[1] + 0.2 * x[2], x[1] + 0.1 * x[0] + 0.1 * x[2], x[2] + 0.1 * x[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] = T(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 = R3{-0.977717523265611, 0.0977717523265611, 0.185766329420466};
+              break;
+            }
+            case 1: {
+              normal = R3{0.977717523265611, -0.0977717523265612, -0.185766329420466};
+              break;
+            }
+            case 2: {
+              normal = R3{0.0911512175788074, -0.992535480302569, 0.0810233045144955};
+              break;
+            }
+            case 3: {
+              normal = R3{-0.0911512175788074, 0.992535480302569, -0.0810233045144955};
+              break;
+            }
+            case 4: {
+              normal = R3{0.100493631166705, -0.0100493631166705, -0.994886948550377};
+              break;
+            }
+            case 5: {
+              normal = R3{-0.100493631166705, 0.0100493631166705, 0.994886948550377};
+              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 (const 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{-0.977717523265611, 0.0977717523265611, 0.185766329420466};
+            } else if (name == "XMAX") {
+              normal = R3{0.977717523265611, -0.0977717523265612, -0.185766329420466};
+            } else if (name == "YMIN") {
+              normal = R3{0.0911512175788074, -0.992535480302569, 0.0810233045144955};
+            } else if (name == "YMAX") {
+              normal = R3{-0.0911512175788074, 0.992535480302569, -0.0810233045144955};
+            } else if (name == "ZMIN") {
+              normal = R3{0.100493631166705, -0.0100493631166705, -0.994886948550377};
+            } else if (name == "ZMAX") {
+              normal = R3{-0.100493631166705, 0.0100493631166705, 0.994886948550377};
+            } else {
+              FAIL("unexpected name: " + name);
+            }
+
+            REQUIRE(l2Norm(node_boundary.outgoingNormal() - normal) == Catch::Approx(0).margin(1E-13));
+          }
+        }
+      }
+    }
+  }
+
   SECTION("curved mesh")
   {
     SECTION("2D")