Add dirty implementation of stencil dealing with symmetric BC

This is a really quick and dirty implementation. Should not be used by now.

src/mesh/Stencil.hpp

@@ -8,8 +8,17 @@ class Stencil
 {
  private:
   ConnectivityMatrix m_stencil;
+#warning TEMPORARY IMPLEMENTATION
+  std::vector<std::pair<std::string, ConnectivityMatrix>> m_symmetry_name_stencil;
 
  public:
+  PUGS_INLINE
+  const auto&
+  symmetryNameStencil() const
+  {
+    return m_symmetry_name_stencil;
+  }
+
   PUGS_INLINE
   auto
   operator[](CellId cell_id) const
@@ -17,7 +26,10 @@ class Stencil
     return m_stencil[cell_id];
   }
 
-  Stencil(const ConnectivityMatrix& stencil) : m_stencil(stencil) {}
+  Stencil(const ConnectivityMatrix& stencil,
+          const std::vector<std::pair<std::string, ConnectivityMatrix>>& symmetry_name_stencil)
+    : m_stencil{stencil}, m_symmetry_name_stencil{symmetry_name_stencil}
+  {}
   Stencil(const Stencil&) = default;
   Stencil(Stencil&&)      = default;
   ~Stencil()              = default;

src/mesh/StencilBuilder.cpp

 #include <mesh/Connectivity.hpp>
 #include <mesh/StencilBuilder.hpp>
 
+#include <mesh/ItemArray.hpp>
 #include <utils/Messenger.hpp>
 
 #warning remove after rework
@@ -109,13 +110,17 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar
 
 template <typename ConnectivityType>
 Stencil
-StencilBuilder::_build(const ConnectivityType& connectivity, size_t number_of_layers) const
+StencilBuilder::_build(const ConnectivityType& connectivity,
+                       size_t number_of_layers,
+                       const std::set<std::string>& symmetry_boundary_names) const
 {
+#warning TEMPORARY
+  if (symmetry_boundary_names.size() == 0) {
     if (number_of_layers == 1) {
       Array<const uint32_t> row_map        = this->_getRowMap(connectivity);
       Array<const uint32_t> column_indices = this->_getColumnIndices(connectivity, row_map);
 
-    return ConnectivityMatrix{row_map, column_indices};
+      return {ConnectivityMatrix{row_map, column_indices}, {}};
     } else {
       if (number_of_layers > 2) {
         throw NotImplementedError("number of layers too large");
@@ -177,23 +182,166 @@ StencilBuilder::_build(const ConnectivityType& connectivity, size_t number_of_la
       for (size_t i = 0; i < column_indices.size(); ++i) {
         column_indices[i] = column_indices_vector[i];
       }
+      ConnectivityMatrix primal_stencil{row_map, column_indices};
 
-    return ConnectivityMatrix{row_map, column_indices};
+      return {primal_stencil, {}};
+    }
+  } else {
+    if constexpr (ConnectivityType::Dimension > 1) {
+      std::vector<Array<const FaceId>> boundary_node_list;
+
+      NodeArray<bool> symmetry_node_list(connectivity, symmetry_boundary_names.size());
+      symmetry_node_list.fill(0);
+
+      auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+      auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+      auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+
+      {
+        size_t i_boundary_name = 0;
+        for (auto boundary_name : symmetry_boundary_names) {
+          bool found = false;
+          for (size_t i_ref_node_list = 0;
+               i_ref_node_list < connectivity.template numberOfRefItemList<ItemType::face>(); ++i_ref_node_list) {
+            const auto& ref_face_list = connectivity.template refItemList<ItemType::face>(i_ref_node_list);
+            if (ref_face_list.refId().tagName() == boundary_name) {
+              found = true;
+              boundary_node_list.push_back(ref_face_list.list());
+              for (size_t i_face = 0; i_face < ref_face_list.list().size(); ++i_face) {
+                const FaceId face_id = ref_face_list.list()[i_face];
+                auto node_list       = face_to_node_matrix[face_id];
+                for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+                  const NodeId node_id = node_list[i_node];
+
+                  symmetry_node_list[node_id][i_boundary_name] = true;
+                }
+              }
+              break;
+            }
+          }
+          ++i_boundary_name;
+          if (not found) {
+#warning IMPROVE ERROR MESSAGE
+            throw NormalError("cannot find boundary " + boundary_name);
+          };
+        }
+      }
+
+      auto cell_is_owned = connectivity.cellIsOwned();
+      auto cell_number   = connectivity.cellNumber();
+
+      Array<uint32_t> row_map{connectivity.numberOfCells() + 1};
+      row_map[0] = 0;
+      std::vector<Array<uint32_t>> symmetry_row_map_list(symmetry_boundary_names.size());
+      for (auto&& symmetry_row_map : symmetry_row_map_list) {
+        symmetry_row_map    = Array<uint32_t>{connectivity.numberOfCells() + 1};
+        symmetry_row_map[0] = 0;
+      }
+
+      std::vector<uint32_t> column_indices_vector;
+      std::vector<std::vector<uint32_t>> symmetry_column_indices_vector(symmetry_boundary_names.size());
+
+      for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+        std::set<CellId> cell_set;
+
+        if (cell_is_owned[cell_id]) {
+          auto cell_node_list = cell_to_node_matrix[cell_id];
+          for (size_t i_cell_node = 0; i_cell_node < cell_node_list.size(); ++i_cell_node) {
+            const NodeId cell_node_id = cell_node_list[i_cell_node];
+            auto node_cell_list       = node_to_cell_matrix[cell_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];
+              if (cell_id != node_cell_id) {
+                cell_set.insert(node_cell_id);
+              }
+            }
+          }
+          row_map[cell_id + 1] = row_map[cell_id] + cell_set.size();
+
+          {
+            std::vector<CellId> cell_vector;
+            for (auto&& set_cell_id : cell_set) {
+              cell_vector.push_back(set_cell_id);
+            }
+            std::sort(cell_vector.begin(), cell_vector.end(),
+                      [&cell_number](const CellId& cell0_id, const CellId& cell1_id) {
+                        return cell_number[cell0_id] < cell_number[cell1_id];
+                      });
+
+            for (auto&& vector_cell_id : cell_vector) {
+              column_indices_vector.push_back(vector_cell_id);
+            }
+          }
+
+          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) {
+              const NodeId cell_node_id = cell_node_list[i_cell_node];
+              if (symmetry_node_list[cell_node_id][i]) {
+                auto node_cell_list = node_to_cell_matrix[cell_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];
+                  symmetry_cell_set.insert(node_cell_id);
+                }
+              }
+            }
+            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);
+
+            std::vector<CellId> cell_vector;
+            for (auto&& set_cell_id : symmetry_cell_set) {
+              cell_vector.push_back(set_cell_id);
+            }
+            std::sort(cell_vector.begin(), cell_vector.end(),
+                      [&cell_number](const CellId& cell0_id, const CellId& cell1_id) {
+                        return cell_number[cell0_id] < cell_number[cell1_id];
+                      });
+
+            for (auto&& vector_cell_id : cell_vector) {
+              symmetry_column_indices_vector[i].push_back(vector_cell_id);
+            }
+          }
+        }
+      }
+      ConnectivityMatrix primal_stencil{row_map, convert_to_array(column_indices_vector)};
+
+      std::vector<std::pair<std::string, ConnectivityMatrix>> symmetry_name_stencil;
+      {
+        size_t i = 0;
+        for (auto&& boundary_name : symmetry_boundary_names) {
+          symmetry_name_stencil.emplace_back(
+            std::make_pair(boundary_name, ConnectivityMatrix{symmetry_row_map_list[i],
+                                                             convert_to_array(symmetry_column_indices_vector[i])}));
+          ++i;
+        }
+      }
+
+      return {{primal_stencil}, {symmetry_name_stencil}};
+
+    } else {
+      throw NotImplementedError("Only implemented in 2D");
+    }
   }
 }
 
 Stencil
-StencilBuilder::build(const IConnectivity& connectivity, size_t number_of_layers) const
+StencilBuilder::build(const IConnectivity& connectivity,
+                      size_t number_of_layers,
+                      const std::set<std::string>& symmetry_boundary_names) const
 {
   switch (connectivity.dimension()) {
   case 1: {
-    return StencilBuilder::_build(dynamic_cast<const Connectivity<1>&>(connectivity), number_of_layers);
+    return StencilBuilder::_build(dynamic_cast<const Connectivity<1>&>(connectivity), number_of_layers,
+                                  symmetry_boundary_names);
   }
   case 2: {
-    return StencilBuilder::_build(dynamic_cast<const Connectivity<2>&>(connectivity), number_of_layers);
+    return StencilBuilder::_build(dynamic_cast<const Connectivity<2>&>(connectivity), number_of_layers,
+                                  symmetry_boundary_names);
   }
   case 3: {
-    return StencilBuilder::_build(dynamic_cast<const Connectivity<3>&>(connectivity), number_of_layers);
+    return StencilBuilder::_build(dynamic_cast<const Connectivity<3>&>(connectivity), number_of_layers,
+                                  symmetry_boundary_names);
   }
   default: {
     throw UnexpectedError("invalid connectivity dimension");

src/mesh/StencilBuilder.hpp

@@ -3,6 +3,9 @@
 
 #include <mesh/Stencil.hpp>
 
+#include <set>
+#include <string>
+
 class IConnectivity;
 class StencilBuilder
 {
@@ -15,10 +18,14 @@ class StencilBuilder
                                           const Array<const uint32_t>& row_map) const;
 
   template <typename ConnectivityType>
-  Stencil _build(const ConnectivityType& connectivity, size_t number_of_layers) const;
+  Stencil _build(const ConnectivityType& connectivity,
+                 size_t number_of_layers,
+                 const std::set<std::string>& symmetry_boundary_names) const;
 
   friend class StencilManager;
-  Stencil build(const IConnectivity& connectivity, size_t number_of_layers) const;
+  Stencil build(const IConnectivity& connectivity,
+                size_t number_of_layers,
+                const std::set<std::string>& symmetry_boundary_names) const;
 
  public:
   StencilBuilder()                      = default;

src/mesh/StencilManager.cpp

@@ -32,14 +32,16 @@ StencilManager::destroy()
 }
 
 const Stencil&
-StencilManager::getStencil(const IConnectivity& connectivity, size_t degree)
+StencilManager::getStencil(const IConnectivity& connectivity,
+                           size_t degree,
+                           const std::set<std::string>& symmetry_boundary_names)
 {
-  if (not m_stored_stencil_map.contains(Key{connectivity.id(), degree})) {
-    m_stored_stencil_map[Key{connectivity.id(), degree}] =
-      std::make_shared<Stencil>(StencilBuilder{}.build(connectivity, degree));
+  if (not m_stored_stencil_map.contains(Key{connectivity.id(), degree, symmetry_boundary_names})) {
+    m_stored_stencil_map[Key{connectivity.id(), degree, symmetry_boundary_names}] =
+      std::make_shared<Stencil>(StencilBuilder{}.build(connectivity, degree, symmetry_boundary_names));
   }
 
-  return *m_stored_stencil_map.at(Key{connectivity.id(), degree});
+  return *m_stored_stencil_map.at(Key{connectivity.id(), degree, symmetry_boundary_names});
 }
 
 void

src/mesh/StencilManager.hpp

@@ -5,6 +5,7 @@
 #include <mesh/Stencil.hpp>
 
 #include <memory>
+#include <set>
 #include <unordered_map>
 
 class StencilManager
@@ -19,11 +20,26 @@ class StencilManager
   {
     size_t connectivity_id;
     size_t degree;
+    std::set<std::string> symmetry_boundary_names;
 
     PUGS_INLINE bool
     operator==(const Key& k) const
     {
-      return (connectivity_id == k.connectivity_id) and (degree == k.degree);
+      if ((connectivity_id != k.connectivity_id) or (degree != k.degree) or
+          (symmetry_boundary_names.size() != k.symmetry_boundary_names.size())) {
+        return false;
+      }
+
+      auto i_name   = symmetry_boundary_names.begin();
+      auto i_k_name = k.symmetry_boundary_names.begin();
+
+      for (; i_name != symmetry_boundary_names.end(); ++i_name, ++i_k_name) {
+        if (*i_name != *i_k_name) {
+          return false;
+        }
+      }
+
+      return true;
     }
   };
   struct HashKey
@@ -31,6 +47,7 @@ class StencilManager
     size_t
     operator()(const Key& key) const
     {
+      // We do not use the symmetry boundary set by now
       return (std::hash<decltype(Key::connectivity_id)>()(key.connectivity_id)) ^
              (std::hash<decltype(Key::degree)>()(key.degree) << 1);
     }
@@ -52,7 +69,9 @@ class StencilManager
 
   void deleteConnectivity(const size_t connectivity_id);
 
-  const Stencil& getStencil(const IConnectivity& i_connectivity, size_t degree = 1);
+  const Stencil& getStencil(const IConnectivity& i_connectivity,
+                            size_t degree                                        = 1,
+                            const std::set<std::string>& symmetry_boundary_names = {});
 
   StencilManager(const StencilManager&) = delete;
   StencilManager(StencilManager&&)      = delete;

src/scheme/PolynomialReconstructionDescriptor.hpp

@@ -5,6 +5,7 @@
 #include <utils/PugsMacros.hpp>
 
 #include <cstddef>
+#include <set>
 
 class PolynomialReconstructionDescriptor
 {
@@ -12,6 +13,9 @@ class PolynomialReconstructionDescriptor
  private:
   IntegrationMethodType m_integration_method;
   size_t m_degree;
+
+  std::set<std::string> m_symmetry_boundary_names;
+
   bool m_preconditioning = true;
   bool m_row_weighting   = true;
 
@@ -29,6 +33,13 @@ class PolynomialReconstructionDescriptor
     return m_degree;
   }
 
+  PUGS_INLINE
+  const std::set<std::string>&
+  symmetryBoundaryNames() const
+  {
+    return m_symmetry_boundary_names;
+  }
+
   PUGS_INLINE
   bool
   preconditioning() const
@@ -64,8 +75,11 @@ class PolynomialReconstructionDescriptor
     m_row_weighting = row_weighting;
   }
 
-  PolynomialReconstructionDescriptor(const IntegrationMethodType integration_method, const size_t degree)
-    : m_integration_method{integration_method}, m_degree{degree}
+#warning TEMPORARY
+  PolynomialReconstructionDescriptor(const IntegrationMethodType integration_method,
+                                     const size_t degree,
+                                     const std::set<std::string>& symmetry_boundary_names = {})
+    : m_integration_method{integration_method}, m_degree{degree}, m_symmetry_boundary_names(symmetry_boundary_names)
   {}
 
   PolynomialReconstructionDescriptor() = delete;

src/scheme/test_reconstruction.cpp

@@ -8,16 +8,11 @@ void
 test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
                     uint64_t degree)
 {
-  std::vector descriptor_list = {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree},
-                                 PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree},
-                                 PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree}};
-
-  [[maybe_unused]] auto x =
-    PolynomialReconstruction{PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree}}.build(
-      discrete_function_variant_list);
+  std::vector descriptor_list = {
+    PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree, {"XMIN", "XMAX", "YMIN", "YMAX"}}};
 
   for (auto&& descriptor : descriptor_list) {
-    const size_t nb_loops = 50;
+    const size_t nb_loops = 1;
     std::cout << "** variable list at once (" << nb_loops << " times) using " << rang::fgB::yellow
               << name(descriptor.integrationMethodType()) << rang::fg::reset << "\n";
 

tests/test_StencilBuilder.cpp

@@ -14,6 +14,8 @@
 // clazy:excludeall=non-pod-global-static
 
 TEST_CASE("StencilBuilder", "[mesh]")
+{
+  SECTION("inner stencil")
   {
     auto is_valid = [](const auto& connectivity, const Stencil& stencil) {
       auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
@@ -113,3 +115,29 @@ TEST_CASE("StencilBuilder", "[mesh]")
       }
     }
   }
+
+  SECTION("Stencil using symmetries")
+  {
+    SECTION("2D")
+    {
+      SECTION("cartesian")
+      {
+        const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
+
+        const Connectivity<2>& connectivity = mesh.connectivity();
+        StencilManager::instance().getStencil(connectivity, 1, {"XMIN", "XMAX", "YMIN", "YMAX"});
+      }
+    }
+
+    SECTION("3D")
+    {
+      SECTION("cartesian")
+      {
+        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"});
+      }
+    }
+  }
+}