From a2eecb6dd6d10c6f2d66dcaba6676df3d0bd7fb0 Mon Sep 17 00:00:00 2001 From: Stephane Del Pino <stephane.delpino44@gmail.com> Date: Thu, 7 Sep 2023 09:08:01 +0200 Subject: [PATCH] Add tests for situation that enlightened normal calculation issue --- tests/test_MeshFlatNodeBoundary.cpp | 205 ++++++++++++++++++++++++++++ 1 file changed, 205 insertions(+) diff --git a/tests/test_MeshFlatNodeBoundary.cpp b/tests/test_MeshFlatNodeBoundary.cpp index 9b9169f73..760f55642 100644 --- a/tests/test_MeshFlatNodeBoundary.cpp +++ b/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") -- GitLab