diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
index 6cd2d798b6eaa2844d4eae0d33f8405d344593b7..5d9c23d6fee898b87808dcd9361018cf275130ce 100644
--- a/src/mesh/ItemToItemMatrix.hpp
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -20,8 +20,8 @@ class ItemToItemMatrix
    private:
     using IndexType = typename ConnectivityMatrix::IndexType;
 
-    const IndexType* const m_values;
     const size_t m_size;
+    const IndexType* const m_values;
 
    public:
     PUGS_INLINE
@@ -40,8 +40,9 @@ class ItemToItemMatrix
 
     PUGS_INLINE
     UnsafeSubItemArray(const ConnectivityMatrix& connectivity_matrix, SourceItemId source_item_id)
-      : m_values{&(connectivity_matrix.values()[connectivity_matrix.rowsMap()[source_item_id]])},
-        m_size(connectivity_matrix.rowsMap()[source_item_id + 1] - connectivity_matrix.rowsMap()[source_item_id])
+      : m_size(connectivity_matrix.rowsMap()[source_item_id + 1] - connectivity_matrix.rowsMap()[source_item_id]),
+        m_values{(m_size == 0) ? nullptr
+                               : (&(connectivity_matrix.values()[connectivity_matrix.rowsMap()[source_item_id]]))}
     {}
 
     PUGS_INLINE
diff --git a/src/mesh/Stencil.hpp b/src/mesh/Stencil.hpp
deleted file mode 100644
index eb4864c7d064ab45b934aa70639697de18cc989e..0000000000000000000000000000000000000000
--- a/src/mesh/Stencil.hpp
+++ /dev/null
@@ -1,74 +0,0 @@
-#ifndef STENCIL_HPP
-#define STENCIL_HPP
-
-#include <mesh/ConnectivityMatrix.hpp>
-#include <mesh/IBoundaryDescriptor.hpp>
-#include <mesh/ItemId.hpp>
-
-class Stencil
-{
- public:
-  class BoundaryDescriptorStencil
-  {
-   private:
-    std::shared_ptr<const IBoundaryDescriptor> m_iboundary_descriptor;
-    ConnectivityMatrix m_stencil;
-
-   public:
-    const IBoundaryDescriptor&
-    boundaryDescriptor() const
-    {
-      return *m_iboundary_descriptor;
-    }
-
-    const ConnectivityMatrix&
-    stencil() const
-    {
-      return m_stencil;
-    }
-
-    BoundaryDescriptorStencil(const std::shared_ptr<const IBoundaryDescriptor>& iboundary_descriptor,
-                              const ConnectivityMatrix& stencil)
-      : m_iboundary_descriptor{iboundary_descriptor}, m_stencil{stencil}
-    {}
-
-    BoundaryDescriptorStencil(const BoundaryDescriptorStencil&) = default;
-    BoundaryDescriptorStencil(BoundaryDescriptorStencil&&)      = default;
-
-    BoundaryDescriptorStencil()  = default;
-    ~BoundaryDescriptorStencil() = default;
-  };
-
-  using BoundaryDescriptorStencilList = std::vector<BoundaryDescriptorStencil>;
-
- private:
-  ConnectivityMatrix m_stencil;
-  BoundaryDescriptorStencilList m_symmetry_boundary_stencil_list;
-
- public:
-#warning REWORK INTERFACE
-  PUGS_INLINE
-  const auto&
-  symmetryBoundaryStencilList() const
-  {
-    return m_symmetry_boundary_stencil_list;
-  }
-
-  PUGS_INLINE
-  auto
-  operator[](CellId cell_id) const
-  {
-    return m_stencil[cell_id];
-  }
-
-  Stencil(const ConnectivityMatrix& stencil, const BoundaryDescriptorStencilList& symmetry_name_stencil)
-    : m_stencil{stencil}, m_symmetry_boundary_stencil_list{symmetry_name_stencil}
-  {}
-
-  Stencil(const Stencil&) = default;
-  Stencil(Stencil&&)      = default;
-
-  ~Stencil() = default;
-};
-
-#endif   // STENCIL_HPP
diff --git a/src/mesh/StencilArray.hpp b/src/mesh/StencilArray.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d1b35def0fc92c19af436cd192eeab7fae846444
--- /dev/null
+++ b/src/mesh/StencilArray.hpp
@@ -0,0 +1,103 @@
+#ifndef STENCIL_ARRAY_HPP
+#define STENCIL_ARRAY_HPP
+
+#include <mesh/ConnectivityMatrix.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/ItemId.hpp>
+#include <mesh/ItemToItemMatrix.hpp>
+
+template <ItemType source_item_type, ItemType target_item_type>
+class StencilArray
+{
+ public:
+  using ItemToItemMatrixT = ItemToItemMatrix<source_item_type, target_item_type>;
+
+  class BoundaryDescriptorStencilArray
+  {
+   private:
+    std::shared_ptr<const IBoundaryDescriptor> m_iboundary_descriptor;
+    const ConnectivityMatrix m_connectivity_matrix;
+    const ItemToItemMatrixT m_stencil_array;
+
+   public:
+    const IBoundaryDescriptor&
+    boundaryDescriptor() const
+    {
+      return *m_iboundary_descriptor;
+    }
+
+    const auto&
+    stencilArray() const
+    {
+      return m_stencil_array;
+    }
+
+    BoundaryDescriptorStencilArray(const std::shared_ptr<const IBoundaryDescriptor>& iboundary_descriptor,
+                                   const ConnectivityMatrix& connectivity_matrix)
+      : m_iboundary_descriptor{iboundary_descriptor},
+        m_connectivity_matrix{connectivity_matrix},
+        m_stencil_array{m_connectivity_matrix}
+    {}
+
+    BoundaryDescriptorStencilArray(const BoundaryDescriptorStencilArray& bdsa)
+      : m_iboundary_descriptor{bdsa.m_iboundary_descriptor},
+        m_connectivity_matrix{bdsa.m_connectivity_matrix},
+        m_stencil_array{m_connectivity_matrix}
+    {}
+
+    BoundaryDescriptorStencilArray(BoundaryDescriptorStencilArray&& bdsa)
+      : m_iboundary_descriptor{std::move(bdsa.m_iboundary_descriptor)},
+        m_connectivity_matrix{std::move(bdsa.m_connectivity_matrix)},
+        m_stencil_array{m_connectivity_matrix}
+    {}
+
+    ~BoundaryDescriptorStencilArray() = default;
+  };
+
+  using BoundaryDescriptorStencilArrayList = std::vector<BoundaryDescriptorStencilArray>;
+
+ private:
+  const ConnectivityMatrix m_connectivity_matrix;
+  const ItemToItemMatrixT m_stencil_array;
+  BoundaryDescriptorStencilArrayList m_symmetry_boundary_stencil_array_list;
+
+ public:
+  PUGS_INLINE
+  const auto&
+  symmetryBoundaryStencilArrayList() const
+  {
+    return m_symmetry_boundary_stencil_array_list;
+  }
+
+  PUGS_INLINE
+  auto
+  operator[](CellId cell_id) const
+  {
+    return m_stencil_array[cell_id];
+  }
+
+  StencilArray(const ConnectivityMatrix& connectivity_matrix,
+               const BoundaryDescriptorStencilArrayList& symmetry_boundary_descriptor_stencil_array_list)
+    : m_connectivity_matrix{connectivity_matrix},
+      m_stencil_array{m_connectivity_matrix},
+      m_symmetry_boundary_stencil_array_list{symmetry_boundary_descriptor_stencil_array_list}
+  {}
+
+  StencilArray(const StencilArray& stencil_array)
+    : m_connectivity_matrix{stencil_array.m_connectivity_matrix},
+      m_stencil_array{m_connectivity_matrix},
+      m_symmetry_boundary_stencil_array_list{stencil_array.m_symmetry_boundary_stencil_array_list}
+  {}
+
+  StencilArray(StencilArray&& stencil_array)
+    : m_connectivity_matrix{std::move(stencil_array.m_connectivity_matrix)},
+      m_stencil_array{m_connectivity_matrix},
+      m_symmetry_boundary_stencil_array_list{std::move(stencil_array.m_symmetry_boundary_stencil_array_list)}
+  {}
+
+  ~StencilArray() = default;
+};
+
+using CellToCellStencilArray = StencilArray<ItemType::cell, ItemType::cell>;
+
+#endif   // STENCIL_ARRAY_HPP
diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp
index 30e922ec24b8303b9dd279ac702c154648a26b29..920dad9e248ef97c2f995cabbc72cff6a24aac7e 100644
--- a/src/mesh/StencilBuilder.cpp
+++ b/src/mesh/StencilBuilder.cpp
@@ -108,12 +108,11 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar
 }
 
 template <typename ConnectivityType>
-Stencil
+CellToCellStencilArray
 StencilBuilder::_build(const ConnectivityType& connectivity,
                        size_t number_of_layers,
                        const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
 {
-#warning TEMPORARY
   if (symmetry_boundary_descriptor_list.size() == 0) {
     if (number_of_layers == 1) {
       Array<const uint32_t> row_map        = this->_getRowMap(connectivity);
@@ -313,15 +312,15 @@ StencilBuilder::_build(const ConnectivityType& connectivity,
       }
       ConnectivityMatrix primal_stencil{row_map, convert_to_array(column_indices_vector)};
 
-      Stencil::BoundaryDescriptorStencilList symmetry_boundary_stencil_list;
+      CellToCellStencilArray::BoundaryDescriptorStencilArrayList symmetry_boundary_stencil_list;
       {
         size_t i = 0;
         for (auto&& p_boundary_descriptor : symmetry_boundary_descriptor_list) {
           symmetry_boundary_stencil_list.emplace_back(
-            Stencil::BoundaryDescriptorStencil{p_boundary_descriptor,
-                                               ConnectivityMatrix{symmetry_row_map_list[i],
-                                                                  convert_to_array(
-                                                                    symmetry_column_indices_vector[i])}});
+            CellToCellStencilArray::
+              BoundaryDescriptorStencilArray{p_boundary_descriptor,
+                                             ConnectivityMatrix{symmetry_row_map_list[i],
+                                                                convert_to_array(symmetry_column_indices_vector[i])}});
           ++i;
         }
       }
@@ -334,7 +333,7 @@ StencilBuilder::_build(const ConnectivityType& connectivity,
   }
 }
 
-Stencil
+CellToCellStencilArray
 StencilBuilder::build(const IConnectivity& connectivity,
                       size_t number_of_layers,
                       const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
diff --git a/src/mesh/StencilBuilder.hpp b/src/mesh/StencilBuilder.hpp
index 7f11c0686a49adebeb4c7ffd5a0c945c6b6fa07d..aaf5076f15cca27b4d199b1b5772dcb573478b31 100644
--- a/src/mesh/StencilBuilder.hpp
+++ b/src/mesh/StencilBuilder.hpp
@@ -2,7 +2,7 @@
 #define STENCIL_BUILDER_HPP
 
 #include <mesh/IBoundaryDescriptor.hpp>
-#include <mesh/Stencil.hpp>
+#include <mesh/StencilArray.hpp>
 
 #include <string>
 #include <vector>
@@ -22,14 +22,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 BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
+  CellToCellStencilArray _build(const ConnectivityType& connectivity,
+                                size_t number_of_layers,
+                                const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
 
   friend class StencilManager;
-  Stencil build(const IConnectivity& connectivity,
-                size_t number_of_layers,
-                const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
+  CellToCellStencilArray build(const IConnectivity& connectivity,
+                               size_t number_of_layers,
+                               const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
 
  public:
   StencilBuilder()                      = default;
diff --git a/src/mesh/StencilManager.cpp b/src/mesh/StencilManager.cpp
index cc99577f23a5d9b85e344d39a1dc0e2260dcb22a..ed1f9ef76ae69937f8c188c9b848d08d3b1dde5c 100644
--- a/src/mesh/StencilManager.cpp
+++ b/src/mesh/StencilManager.cpp
@@ -18,10 +18,10 @@ StencilManager::destroy()
   Assert(m_instance != nullptr, "StencilManager was not created");
 
   // LCOV_EXCL_START
-  if (m_instance->m_stored_stencil_map.size() > 0) {
+  if (m_instance->m_stored_cell_to_cell_stencil_map.size() > 0) {
     std::stringstream error;
     error << ": some connectivities are still registered\n";
-    for (const auto& [key, stencil_p] : m_instance->m_stored_stencil_map) {
+    for (const auto& [key, stencil_p] : m_instance->m_stored_cell_to_cell_stencil_map) {
       error << " - connectivity of id " << rang::fgB::magenta << key.connectivity_id << rang::style::reset << '\n';
     }
     throw UnexpectedError(error.str());
@@ -31,17 +31,19 @@ StencilManager::destroy()
   m_instance = nullptr;
 }
 
-const Stencil&
-StencilManager::getStencil(const IConnectivity& connectivity,
-                           size_t degree,
-                           const BoundaryDescriptorList& symmetry_boundary_descriptor_list)
+const CellToCellStencilArray&
+StencilManager::getCellToCellStencilArray(const IConnectivity& connectivity,
+                                          size_t degree,
+                                          const BoundaryDescriptorList& symmetry_boundary_descriptor_list)
 {
-  if (not m_stored_stencil_map.contains(Key{connectivity.id(), degree, symmetry_boundary_descriptor_list})) {
-    m_stored_stencil_map[Key{connectivity.id(), degree, symmetry_boundary_descriptor_list}] =
-      std::make_shared<Stencil>(StencilBuilder{}.build(connectivity, degree, symmetry_boundary_descriptor_list));
+  if (not m_stored_cell_to_cell_stencil_map.contains(
+        Key{connectivity.id(), degree, symmetry_boundary_descriptor_list})) {
+    m_stored_cell_to_cell_stencil_map[Key{connectivity.id(), degree, symmetry_boundary_descriptor_list}] =
+      std::make_shared<CellToCellStencilArray>(
+        StencilBuilder{}.build(connectivity, degree, symmetry_boundary_descriptor_list));
   }
 
-  return *m_stored_stencil_map.at(Key{connectivity.id(), degree, symmetry_boundary_descriptor_list});
+  return *m_stored_cell_to_cell_stencil_map.at(Key{connectivity.id(), degree, symmetry_boundary_descriptor_list});
 }
 
 void
@@ -50,9 +52,9 @@ StencilManager::deleteConnectivity(const size_t connectivity_id)
   bool has_removed = false;
   do {
     has_removed = false;
-    for (const auto& [key, p_stencil] : m_stored_stencil_map) {
+    for (const auto& [key, p_stencil] : m_stored_cell_to_cell_stencil_map) {
       if (connectivity_id == key.connectivity_id) {
-        m_stored_stencil_map.erase(key);
+        m_stored_cell_to_cell_stencil_map.erase(key);
         has_removed = true;
         break;
       }
diff --git a/src/mesh/StencilManager.hpp b/src/mesh/StencilManager.hpp
index 814d8e476d71ada4ec66322e1f3ed80035f808e2..de0007631dc6a7b0d788cdd1f635275fa35b64e5 100644
--- a/src/mesh/StencilManager.hpp
+++ b/src/mesh/StencilManager.hpp
@@ -3,7 +3,7 @@
 
 #include <mesh/IBoundaryDescriptor.hpp>
 #include <mesh/IConnectivity.hpp>
-#include <mesh/Stencil.hpp>
+#include <mesh/StencilArray.hpp>
 
 #include <memory>
 #include <set>
@@ -50,6 +50,7 @@ class StencilManager
       return boundary_descriptor_set == k_boundary_descriptor_set;
     }
   };
+
   struct HashKey
   {
     size_t
@@ -61,7 +62,7 @@ class StencilManager
     }
   };
 
-  std::unordered_map<Key, std::shared_ptr<const Stencil>, HashKey> m_stored_stencil_map;
+  std::unordered_map<Key, std::shared_ptr<const CellToCellStencilArray>, HashKey> m_stored_cell_to_cell_stencil_map;
 
  public:
   static void create();
@@ -77,9 +78,10 @@ class StencilManager
 
   void deleteConnectivity(const size_t connectivity_id);
 
-  const Stencil& getStencil(const IConnectivity& i_connectivity,
-                            size_t degree                                                   = 1,
-                            const BoundaryDescriptorList& symmetry_boundary_descriptor_list = {});
+  const CellToCellStencilArray& getCellToCellStencilArray(
+    const IConnectivity& i_connectivity,
+    size_t degree                                                   = 1,
+    const BoundaryDescriptorList& symmetry_boundary_descriptor_list = {});
 
   StencilManager(const StencilManager&) = delete;
   StencilManager(StencilManager&&)      = delete;
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 7bcf3a5c68407bc0fb50d6cf46560de4178b635f..82f5081a772dd39cc814436b9adace162a670a9e 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -769,8 +769,9 @@ PolynomialReconstruction::_build(
   const size_t basis_dimension =
     DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(m_descriptor.degree());
 
-  const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity(), m_descriptor.degree(),
-                                                                 m_descriptor.symmetryBoundaryDescriptorList());
+  const auto& stencil_array =
+    StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), m_descriptor.degree(),
+                                                         m_descriptor.symmetryBoundaryDescriptorList());
 
   auto xr = mesh.xr();
   auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
@@ -782,10 +783,10 @@ PolynomialReconstruction::_build(
   auto cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
 
   auto full_stencil_size = [&](const CellId cell_id) {
-    auto stencil_cell_list = stencil[cell_id];
+    auto stencil_cell_list = stencil_array[cell_id];
     size_t stencil_size    = stencil_cell_list.size();
     for (size_t i = 0; i < m_descriptor.symmetryBoundaryDescriptorList().size(); ++i) {
-      auto& ghost_stencil = stencil.symmetryBoundaryStencilList()[i].stencil();
+      auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i].stencilArray();
       stencil_size += ghost_stencil[cell_id].size();
     }
 
@@ -814,6 +815,7 @@ PolynomialReconstruction::_build(
     }
     return normal_list;
   }();
+
   SmallArray<const Rd> symmetry_origin_list = [&] {
     SmallArray<Rd> origin_list(m_descriptor.symmetryBoundaryDescriptorList().size());
     size_t i_symmetry_boundary = 0;
@@ -863,7 +865,7 @@ PolynomialReconstruction::_build(
       if (cell_is_owned[cell_j_id]) {
         const int32_t t = tokens.acquire();
 
-        auto stencil_cell_list = stencil[cell_j_id];
+        auto stencil_cell_list = stencil_array[cell_j_id];
 
         ShrinkMatrixView B(B_pool[t], full_stencil_size(cell_j_id));
 
@@ -897,8 +899,9 @@ PolynomialReconstruction::_build(
                   }
                 }
 
-                for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryBoundaryStencilList().size(); ++i_symmetry) {
-                  auto& ghost_stencil  = stencil.symmetryBoundaryStencilList()[i_symmetry].stencil();
+                for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                     ++i_symmetry) {
+                  auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
                   auto ghost_cell_list = ghost_stencil[cell_j_id];
                   for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
                     const CellId cell_i_id = ghost_cell_list[i];
@@ -953,8 +956,9 @@ PolynomialReconstruction::_build(
                   }
                 }
 
-                for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryBoundaryStencilList().size(); ++i_symmetry) {
-                  auto& ghost_stencil  = stencil.symmetryBoundaryStencilList()[i_symmetry].stencil();
+                for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                     ++i_symmetry) {
+                  auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
                   auto ghost_cell_list = ghost_stencil[cell_j_id];
                   for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
                     const CellId cell_i_id = ghost_cell_list[i];
@@ -988,8 +992,9 @@ PolynomialReconstruction::_build(
               A(index, l) = Xi_Xj[l];
             }
           }
-          for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryBoundaryStencilList().size(); ++i_symmetry) {
-            auto& ghost_stencil  = stencil.symmetryBoundaryStencilList()[i_symmetry].stencil();
+          for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+               ++i_symmetry) {
+            auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
             auto ghost_cell_list = ghost_stencil[cell_j_id];
 
             const Rd& origin = symmetry_origin_list[i_symmetry];
@@ -1033,8 +1038,9 @@ PolynomialReconstruction::_build(
                   A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
                 }
               }
-              for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryBoundaryStencilList().size(); ++i_symmetry) {
-                auto& ghost_stencil  = stencil.symmetryBoundaryStencilList()[i_symmetry].stencil();
+              for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                   ++i_symmetry) {
+                auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
                 auto ghost_cell_list = ghost_stencil[cell_j_id];
 
                 const Rd& origin = symmetry_origin_list[i_symmetry];
@@ -1078,8 +1084,9 @@ PolynomialReconstruction::_build(
                 A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
               }
             }
-            for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryBoundaryStencilList().size(); ++i_symmetry) {
-              auto& ghost_stencil  = stencil.symmetryBoundaryStencilList()[i_symmetry].stencil();
+            for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                 ++i_symmetry) {
+              auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
               auto ghost_cell_list = ghost_stencil[cell_j_id];
 
               const Rd& origin = symmetry_origin_list[i_symmetry];
@@ -1119,8 +1126,9 @@ PolynomialReconstruction::_build(
               B(index, l) *= weight;
             }
           }
-          for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryBoundaryStencilList().size(); ++i_symmetry) {
-            auto& ghost_stencil  = stencil.symmetryBoundaryStencilList()[i_symmetry].stencil();
+          for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+               ++i_symmetry) {
+            auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
             auto ghost_cell_list = ghost_stencil[cell_j_id];
 
             const Rd& origin = symmetry_origin_list[i_symmetry];
diff --git a/tests/test_StencilBuilder.cpp b/tests/test_StencilBuilder.cpp
index 754a443a7d63d391807d8133110e04ec8b2b224f..b68f0883d2319fe61a4e42db7bd57bfab59d877f 100644
--- a/tests/test_StencilBuilder.cpp
+++ b/tests/test_StencilBuilder.cpp
@@ -19,7 +19,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
 {
   SECTION("inner stencil")
   {
-    auto is_valid = [](const auto& connectivity, const Stencil& stencil) {
+    auto is_valid = [](const auto& connectivity, const auto& stencil_array) {
       auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
       auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
 
@@ -42,7 +42,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
             }
           }
 
-          auto cell_stencil = stencil[cell_id];
+          auto cell_stencil = stencil_array[cell_id];
 
           auto i_set_cell = cell_set.begin();
           for (size_t index = 0; index < cell_stencil.size(); ++index, ++i_set_cell) {
@@ -63,9 +63,9 @@ TEST_CASE("StencilBuilder", "[mesh]")
 
         const Connectivity<1>& connectivity = mesh.connectivity();
 
-        Stencil stencil = StencilManager::instance().getStencil(connectivity);
+        auto stencil_array = StencilManager::instance().getCellToCellStencilArray(connectivity);
 
-        REQUIRE(is_valid(connectivity, stencil));
+        REQUIRE(is_valid(connectivity, stencil_array));
       }
 
       SECTION("unordered")
@@ -73,9 +73,9 @@ TEST_CASE("StencilBuilder", "[mesh]")
         const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
 
         const Connectivity<1>& connectivity = mesh.connectivity();
-        Stencil stencil                     = StencilManager::instance().getStencil(connectivity);
+        auto stencil_array                  = StencilManager::instance().getCellToCellStencilArray(connectivity);
 
-        REQUIRE(is_valid(connectivity, stencil));
+        REQUIRE(is_valid(connectivity, stencil_array));
       }
     }
 
@@ -86,7 +86,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
         const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
 
         const Connectivity<2>& connectivity = mesh.connectivity();
-        REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+        REQUIRE(is_valid(connectivity, StencilManager::instance().getCellToCellStencilArray(connectivity)));
       }
 
       SECTION("hybrid")
@@ -94,7 +94,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
         const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
 
         const Connectivity<2>& connectivity = mesh.connectivity();
-        REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+        REQUIRE(is_valid(connectivity, StencilManager::instance().getCellToCellStencilArray(connectivity)));
       }
     }
 
@@ -105,7 +105,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
         const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
 
         const Connectivity<3>& connectivity = mesh.connectivity();
-        REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+        REQUIRE(is_valid(connectivity, StencilManager::instance().getCellToCellStencilArray(connectivity)));
       }
 
       SECTION("hybrid")
@@ -113,24 +113,26 @@ TEST_CASE("StencilBuilder", "[mesh]")
         const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
 
         const Connectivity<3>& connectivity = mesh.connectivity();
-        REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+        REQUIRE(is_valid(connectivity, StencilManager::instance().getCellToCellStencilArray(connectivity)));
       }
     }
   }
 
   SECTION("Stencil using symmetries")
   {
-    auto check_ghost_cells_have_empty_stencils = [](const auto& stencil, const auto& connecticity) {
+    auto check_ghost_cells_have_empty_stencils = [](const auto& stencil_array, 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.symmetryBoundaryStencilList().size();
-               ++i_symmetry_stencil) {
-            is_empty &= (stencil.symmetryBoundaryStencilList()[i_symmetry_stencil].stencil()[cell_id].size() == 0);
+          is_empty &= (stencil_array[cell_id].size() == 0);
+          for (size_t i_symmetry_stencil = 0;
+               i_symmetry_stencil < stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry_stencil) {
+            is_empty &=
+              (stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray()[cell_id].size() ==
+               0);
           }
         }
       }
@@ -138,7 +140,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
       return is_empty;
     };
 
-    auto symmetry_stencils_are_valid = [](const auto& stencil, const auto& mesh) {
+    auto symmetry_stencils_are_valid = [](const auto& stencil_array, const auto& mesh) {
       bool is_valid = true;
 
       auto node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
@@ -146,12 +148,12 @@ TEST_CASE("StencilBuilder", "[mesh]")
       auto cell_is_owned       = mesh.connectivity().cellIsOwned();
       auto cell_number         = mesh.connectivity().cellNumber();
 
-      for (size_t i_symmetry_stencil = 0; i_symmetry_stencil < stencil.symmetryBoundaryStencilList().size();
+      for (size_t i_symmetry_stencil = 0; i_symmetry_stencil < stencil_array.symmetryBoundaryStencilArrayList().size();
            ++i_symmetry_stencil) {
         const IBoundaryDescriptor& boundary_descriptor =
-          stencil.symmetryBoundaryStencilList()[i_symmetry_stencil].boundaryDescriptor();
+          stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].boundaryDescriptor();
 
-        auto boundary_stencil   = stencil.symmetryBoundaryStencilList()[i_symmetry_stencil].stencil();
+        auto boundary_stencil   = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray();
         auto boundary_node_list = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
 
         CellValue<bool> boundary_cell{mesh.connectivity()};
@@ -219,11 +221,11 @@ TEST_CASE("StencilBuilder", "[mesh]")
 
         const Connectivity<2>& connectivity = mesh.connectivity();
 
-        auto stencil = StencilManager::instance().getStencil(connectivity, 1, list);
+        auto stencil_array = StencilManager::instance().getCellToCellStencilArray(connectivity, 1, list);
 
-        REQUIRE(stencil.symmetryBoundaryStencilList().size() == 4);
-        REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
-        REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
+        REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4);
+        REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+        REQUIRE(symmetry_stencils_are_valid(stencil_array, mesh));
       }
 
       SECTION("hybrid")
@@ -231,9 +233,9 @@ TEST_CASE("StencilBuilder", "[mesh]")
         const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
 
         const Connectivity<2>& connectivity = mesh.connectivity();
-        auto stencil                        = StencilManager::instance().getStencil(connectivity, 1, list);
+        auto stencil = StencilManager::instance().getCellToCellStencilArray(connectivity, 1, list);
 
-        REQUIRE(stencil.symmetryBoundaryStencilList().size() == 4);
+        REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 4);
         REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
         REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
       }
@@ -254,9 +256,9 @@ 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, list);
+        auto stencil = StencilManager::instance().getCellToCellStencilArray(connectivity, 1, list);
 
-        REQUIRE(stencil.symmetryBoundaryStencilList().size() == 6);
+        REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 6);
         REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
         REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
       }
@@ -266,9 +268,9 @@ TEST_CASE("StencilBuilder", "[mesh]")
         const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
 
         const Connectivity<3>& connectivity = mesh.connectivity();
-        auto stencil                        = StencilManager::instance().getStencil(connectivity, 1, list);
+        auto stencil = StencilManager::instance().getCellToCellStencilArray(connectivity, 1, list);
 
-        REQUIRE(stencil.symmetryBoundaryStencilList().size() == 6);
+        REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 6);
         REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
         REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
       }