From 827233cf2d4123de07bf5151af314a3b07ade9d2 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Mon, 12 Aug 2024 18:37:07 +0200
Subject: [PATCH] Fix stencil builder with symmetry and add better tests
---
src/mesh/StencilBuilder.cpp | 24 ++++---
tests/test_StencilBuilder.cpp | 123 +++++++++++++++++++++++++++++++++-
2 files changed, 135 insertions(+), 12 deletions(-)
diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp
index 92738a643..209007a2d 100644
--- a/src/mesh/StencilBuilder.cpp
+++ b/src/mesh/StencilBuilder.cpp
@@ -1,10 +1,9 @@
-#include <mesh/Connectivity.hpp>
#include <mesh/StencilBuilder.hpp>
+#include <mesh/Connectivity.hpp>
#include <mesh/ItemArray.hpp>
#include <utils/Messenger.hpp>
-#warning remove after rework
#include <set>
template <typename ConnectivityType>
@@ -221,9 +220,10 @@ StencilBuilder::_build(const ConnectivityType& connectivity,
}
++i_boundary_name;
if (not found) {
-#warning IMPROVE ERROR MESSAGE
- throw NormalError("cannot find boundary " + boundary_name);
- };
+ std::ostringstream error_msg;
+ error_msg << "cannot find boundary '" << rang::fgB::yellow << boundary_name << rang::fg::reset << '\'';
+ throw NormalError(error_msg.str());
+ }
}
}
@@ -243,6 +243,7 @@ StencilBuilder::_build(const ConnectivityType& connectivity,
for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
std::set<CellId> cell_set;
+ std::vector<std::set<CellId>> by_boundary_symmetry_cell(symmetry_boundary_names.size());
if (cell_is_owned[cell_id]) {
auto cell_node_list = cell_to_node_matrix[cell_id];
@@ -256,7 +257,6 @@ StencilBuilder::_build(const ConnectivityType& connectivity,
}
}
}
- row_map[cell_id + 1] = row_map[cell_id] + cell_set.size();
{
std::vector<CellId> cell_vector;
@@ -273,7 +273,6 @@ StencilBuilder::_build(const ConnectivityType& connectivity,
}
}
- std::vector<std::set<CellId>> by_boundary_symmetry_cell;
for (size_t i = 0; i < symmetry_boundary_names.size(); ++i) {
std::set<CellId> symmetry_cell_set;
for (size_t i_cell_node = 0; i_cell_node < cell_node_list.size(); ++i_cell_node) {
@@ -286,8 +285,7 @@ StencilBuilder::_build(const ConnectivityType& connectivity,
}
}
}
- symmetry_row_map_list[i][cell_id + 1] = symmetry_row_map_list[i][cell_id] + symmetry_cell_set.size();
- by_boundary_symmetry_cell.emplace_back(symmetry_cell_set);
+ by_boundary_symmetry_cell[i] = symmetry_cell_set;
std::vector<CellId> cell_vector;
for (auto&& set_cell_id : symmetry_cell_set) {
@@ -303,6 +301,12 @@ StencilBuilder::_build(const ConnectivityType& connectivity,
}
}
}
+ row_map[cell_id + 1] = row_map[cell_id] + cell_set.size();
+
+ for (size_t i = 0; i < symmetry_row_map_list.size(); ++i) {
+ symmetry_row_map_list[i][cell_id + 1] =
+ symmetry_row_map_list[i][cell_id] + by_boundary_symmetry_cell[i].size();
+ }
}
ConnectivityMatrix primal_stencil{row_map, convert_to_array(column_indices_vector)};
@@ -320,7 +324,7 @@ StencilBuilder::_build(const ConnectivityType& connectivity,
return {{primal_stencil}, {symmetry_name_stencil}};
} else {
- throw NotImplementedError("Only implemented in 2D");
+ throw NotImplementedError("Only implemented in 2D/3D");
}
}
}
diff --git a/tests/test_StencilBuilder.cpp b/tests/test_StencilBuilder.cpp
index 54dd97556..32ec498e6 100644
--- a/tests/test_StencilBuilder.cpp
+++ b/tests/test_StencilBuilder.cpp
@@ -7,7 +7,9 @@
#include <mesh/ItemValue.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshVariant.hpp>
+#include <mesh/NamedBoundaryDescriptor.hpp>
#include <mesh/StencilManager.hpp>
#include <utils/Messenger.hpp>
@@ -118,6 +120,89 @@ TEST_CASE("StencilBuilder", "[mesh]")
SECTION("Stencil using symmetries")
{
+ auto check_ghost_cells_have_empty_stencils = [](const auto& stencil, const auto& connecticity) {
+ bool is_empty = true;
+
+ auto cell_is_owned = connecticity.cellIsOwned();
+
+ for (CellId cell_id = 0; cell_id < connecticity.numberOfCells(); ++cell_id) {
+ if (not cell_is_owned[cell_id]) {
+ is_empty &= (stencil[cell_id].size() == 0);
+ for (size_t i_symmetry_stencil = 0; i_symmetry_stencil < stencil.symmetryNameStencil().size();
+ ++i_symmetry_stencil) {
+ is_empty &= (stencil.symmetryNameStencil()[i_symmetry_stencil].second[cell_id].size() == 0);
+ }
+ }
+ }
+
+ return is_empty;
+ };
+
+ auto symmetry_stencils_are_valid = [](const auto& stencil, const auto& mesh) {
+ bool is_valid = true;
+
+ auto node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ auto cell_is_owned = mesh.connectivity().cellIsOwned();
+ auto cell_number = mesh.connectivity().cellNumber();
+
+ for (size_t i_symmetry_stencil = 0; i_symmetry_stencil < stencil.symmetryNameStencil().size();
+ ++i_symmetry_stencil) {
+ auto boundary_name = stencil.symmetryNameStencil()[i_symmetry_stencil].first;
+ auto boundary_stencil = stencil.symmetryNameStencil()[i_symmetry_stencil].second;
+ auto boundary_node_list = getMeshFlatNodeBoundary(mesh, NamedBoundaryDescriptor(boundary_name));
+
+ CellValue<bool> boundary_cell{mesh.connectivity()};
+ boundary_cell.fill(false);
+ auto node_list = boundary_node_list.nodeList();
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+ auto node_cell_list = node_to_cell_matrix[node_id];
+ for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+ const CellId cell_id = node_cell_list[i_cell];
+ boundary_cell[cell_id] = true;
+ }
+ }
+
+ std::set<NodeId> symmetry_node;
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+ symmetry_node.insert(node_id);
+ }
+
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ if ((not boundary_cell[cell_id]) or (not cell_is_owned[cell_id])) {
+ is_valid &= (boundary_stencil[cell_id].size() == 0);
+ } else {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
+ [&](CellId cell0_id, CellId cell1_id) { return cell_number[cell0_id] < cell_number[cell1_id]; });
+ for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
+ const NodeId node_id = cell_node_list[i_node];
+ if (symmetry_node.contains(node_id)) {
+ auto node_cell_list = node_to_cell_matrix[node_id];
+ for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+ const CellId node_cell_id = node_cell_list[i_node_cell];
+ cell_set.insert(node_cell_id);
+ }
+ }
+ }
+
+ if (cell_set.size() == boundary_stencil[cell_id].size()) {
+ size_t i = 0;
+ for (auto&& id : cell_set) {
+ is_valid &= (id == boundary_stencil[cell_id][i++]);
+ }
+ } else {
+ is_valid = false;
+ }
+ }
+ }
+ }
+
+ return is_valid;
+ };
+
SECTION("2D")
{
SECTION("cartesian")
@@ -125,7 +210,23 @@ TEST_CASE("StencilBuilder", "[mesh]")
const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
const Connectivity<2>& connectivity = mesh.connectivity();
- StencilManager::instance().getStencil(connectivity, 1, {"XMIN", "XMAX", "YMIN", "YMAX"});
+ auto stencil = StencilManager::instance().getStencil(connectivity, 1, {"XMIN", "XMAX", "YMIN", "YMAX"});
+
+ REQUIRE(stencil.symmetryNameStencil().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
+ REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ auto stencil = StencilManager::instance().getStencil(connectivity, 1, {"XMIN", "XMAX", "YMIN", "YMAX"});
+
+ REQUIRE(stencil.symmetryNameStencil().size() == 4);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
+ REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
}
}
@@ -136,7 +237,25 @@ TEST_CASE("StencilBuilder", "[mesh]")
const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
const Connectivity<3>& connectivity = mesh.connectivity();
- StencilManager::instance().getStencil(connectivity, 1, {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"});
+ auto stencil =
+ StencilManager::instance().getStencil(connectivity, 1, {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"});
+
+ REQUIRE(stencil.symmetryNameStencil().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
+ REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ auto stencil =
+ StencilManager::instance().getStencil(connectivity, 1, {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"});
+
+ REQUIRE(stencil.symmetryNameStencil().size() == 6);
+ REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
+ REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
}
}
}
--
GitLab