Fix stencil builder with symmetry and add better tests

src/mesh/StencilBuilder.cpp

-#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");
     }
   }
 }

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();
+        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));
       }
     }
   }