diff --git a/cmake/PugsDoc.cmake b/cmake/PugsDoc.cmake
index 2be770b3a3163c6ea490413b0f73a1f4f828c82f..7718046e21b5de20bc5f4a6aeb5c451ca9757890 100644
--- a/cmake/PugsDoc.cmake
+++ b/cmake/PugsDoc.cmake
@@ -20,13 +20,10 @@ add_custom_target(doc DEPENDS userdoc)
 
 if (EMACS AND GNUPLOT_FOUND AND GMSH)
 
-  add_custom_command(
-    OUTPUT "${PUGS_BINARY_DIR}/doc"
+  add_custom_target(pugsdoc-dir
     COMMAND ${CMAKE_COMMAND} -E make_directory "${PUGS_BINARY_DIR}/doc"
   )
 
-  add_custom_target(pugsdoc-dir DEPENDS "${PUGS_BINARY_DIR}/doc")
-
   set(ORG_GENERATOR_FILES
     "${PUGS_SOURCE_DIR}/doc/lisp/build-doc-config.el"
     "${PUGS_SOURCE_DIR}/doc/lisp/share/pugs.el"
@@ -102,16 +99,10 @@ if (EMACS AND GNUPLOT_FOUND AND GMSH)
       FILE_PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE
       @ONLY)
 
-    set_source_files_properties(
-      ${PUGS_BINARY_DIR}/build-pdf.sh2
-      PROPERTIES
-      GENERATED TRUE
-      HEADER_FILE_ONLY TRUE
-    )
-
     add_custom_target(userdoc-pdf DEPENDS pugsdoc-dir "${PUGS_BINARY_DIR}/doc/userdoc.pdf" "${PUGS_BINARY_DIR}/doc/build-userdoc-pdf.sh")
 
     add_dependencies(userdoc userdoc-pdf)
+    add_dependencies(userdoc-html userdoc-pdf)
 
   else()
     if (NOT LATEX_PDFLATEX_FOUND)
diff --git a/doc/userdoc.org b/doc/userdoc.org
index eaac7da46edb913aff2df4d527d7c9ac20456cfe..e9ba680359c371e480f525378931928127c1177d 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -21,6 +21,7 @@
 #+LATEX_HEADER: \usepackage{siunitx}
 #+LATEX_HEADER: \usepackage[table]{xcolor}
 #+LATEX_HEADER: \usepackage{colortbl}
+#+LATEX_HEADER: \hypersetup{linktoc = all, colorlinks = true, linkcolor = blue}
 #+LATEX_COMPILER: pdflatex --shell-escape
 
 #+LATEX_HEADER_EXTRA: \usepackage{amsmath}
@@ -4740,11 +4741,14 @@ Running this code produces the gnuplot file displayed in Figure
 
 ***** ~gnuplot~ writers <<gnuplot-writers>>
 
-There are two ~gnuplot~ writers. One is dedicated to output of dimension
-1 results (~gnuplot_1d_writer~) and the other allows post processing in
-dimension 1 and 2 (~gnuplot_writer~).
+There are three ~gnuplot~ writers. One is dedicated to output of
+dimension 1 results (~gnuplot_1d_writer~), the second allows post
+processing in dimension 1 and 2 (~gnuplot_writer~) and the third can be
+used to write arbitrary data in the ~gnuplot~ format, whichever the
+space dimension is (~gnuplot_raw_writer~). The last one is useful to
+create scatter plots for instance.
 
-Both of these writers can be used for single output or time series
+All of these writers can be used for single output or time series
 outputs. In the case of single output, the filename is completed by
 adding the extension ~.gnu~, in the case of time series, the filename is
 extending by adding ~.abcd.gnu~, where ~abcd~ is the number in the output
@@ -5006,6 +5010,58 @@ The gnuplot result is displayed on Figure [[fig:writer-gp-2d-cos-sin]].
 This is the time series function in the case of the ~gnuplot_writer~. It
 behaves similarly to [[gnuplot-1d-series]].
 
+****** ~gnuplot_raw_writer~ functions <<gnuplot-raw-writer-function>>
+
+These writers are used to store raw data to the ~gnuplot~ format. Doing
+so, it is the user responsibility to store coordinates in the file if
+needed. Therefore it can be used in any space dimension.
+
+******* ~gnuplot_raw_writer: string -> writer~
+
+Here is an illustrating example in dimension 3, see the result on
+Figure [[fig:writer-raw-sin]].
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0,0,0], [2,2,2], (10,10,10));
+
+  let pi:R, pi = acos(-1);
+  let sin_pi_x:R^3 -> R, x -> sin(pi*dot(x,x));
+
+  let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
+
+  let r:Vh, r = sqrt(dot(cell_center(m),cell_center(m)));
+
+  write(gnuplot_raw_writer("raw_sin"), (name_output(r,"r"), name_output(fh, "f")));
+#+END_SRC
+
+#+NAME: writer-raw-sin-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/raw-sin.png")
+  reset
+  set grid
+  set border
+  unset key
+  set xtics
+  set ytics
+  set square
+  set terminal png truecolor enhanced size 628,400
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/raw_sin.gnu)' lw 2 w lp
+#+END_SRC
+
+#+CAPTION: Example of produced gnuplot results from the ~gnuplot_raw_writer~.
+#+NAME: fig:writer-raw-sin
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: writer-raw-sin-img
+
+******* ~gnuplot_raw_writer: string*R -> writer~
+
+This is the time series function in the case of the ~gnuplot_raw_writer~. It
+behaves similarly to [[gnuplot-1d-series]].
+
 ***** ~vtk~ writers
 
 For more complex post processing (including 3d), ~pugs~ can generate ~vtk~
diff --git a/src/analysis/PyramidGaussQuadrature.cpp b/src/analysis/PyramidGaussQuadrature.cpp
index 679f3ded5b85764bca43ed062f379d8171cd5c65..e85a74ac1af9ec6893d15f5489cd568891e0ca92 100644
--- a/src/analysis/PyramidGaussQuadrature.cpp
+++ b/src/analysis/PyramidGaussQuadrature.cpp
@@ -24,6 +24,11 @@ PyramidGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
 
       const double w = (4. / 3) * unit_weight;
 
+      // gcc's bound checking are messed up due to the use of
+      // std::array and the following dynamic/general switch
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Warray-bounds"
+
       switch (id) {
       case 1: {
         Assert(value_list.size() == 1);
@@ -96,6 +101,8 @@ PyramidGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
       }
         // LCOV_EXCL_STOP
       }
+
+#pragma GCC diagnostic pop
     }
   };
 
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 13ce019b451e62b22332c9c4eeff9e277a57bf32..19beaad938376248256734f196f29e181def9451 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -51,7 +51,6 @@
 #include <scheme/RusanovEulerianCompositeSolver_v2_o2.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <scheme/VariableBCDescriptor.hpp>
-#include <scheme/test_reconstruction.hpp>
 #include <utils/Socket.hpp>
 
 #include <memory>
@@ -1006,6 +1005,27 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("cell_center",
+                            std::function(
+
+                              [](const std::shared_ptr<const MeshVariant>& mesh_v)
+                                -> std::shared_ptr<const DiscreteFunctionVariant> {
+                                return std::visit(
+                                  [&](auto&& mesh) {
+                                    using MeshType = mesh_type_t<decltype(mesh)>;
+                                    if constexpr (is_polygonal_mesh_v<MeshType>) {
+                                      return std::make_shared<DiscreteFunctionVariant>(
+                                        DiscreteFunctionP0(mesh_v,
+                                                           MeshDataManager::instance().getMeshData(*mesh).xj()));
+                                    } else {
+                                      throw NormalError("unexpected mesh type");
+                                    }
+                                  },
+                                  mesh_v->variant());
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("hyperelastic_dt", std::function(
 
                                                  [](const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
@@ -1025,29 +1045,6 @@ SchemeModule::SchemeModule()
 
                                                    ));
 
-#warning REMOVE AFTER TESTS FINISHED
-  this->_addBuiltinFunction("test_reconstruction",
-                            std::function(
-
-                              [](const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>&
-                                   discrete_function_variant_list) -> void {
-                                test_reconstruction(discrete_function_variant_list);
-                              }
-
-                              ));
-
-#warning REMOVE AFTER TESTS FINISHED
-  this->_addBuiltinFunction("test_reconstruction",
-                            std::function(
-
-                              [](const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>&
-                                   discrete_function_variant_list,
-                                 uint64_t degree) -> void {
-                                test_reconstruction(discrete_function_variant_list, degree);
-                              }
-
-                              ));
-
   MathFunctionRegisterForVh{*this};
 }
 
diff --git a/src/language/modules/WriterModule.cpp b/src/language/modules/WriterModule.cpp
index fc9d5703ed793f23c7f9768f4f2137585c3be318..4db29e203492c098a4807665ee50eb68b3be8287 100644
--- a/src/language/modules/WriterModule.cpp
+++ b/src/language/modules/WriterModule.cpp
@@ -9,6 +9,7 @@
 #include <mesh/Mesh.hpp>
 #include <output/GnuplotWriter.hpp>
 #include <output/GnuplotWriter1D.hpp>
+#include <output/GnuplotWriterRaw.hpp>
 #include <output/IWriter.hpp>
 #include <output/NamedDiscreteFunction.hpp>
 #include <output/NamedItemArrayVariant.hpp>
@@ -73,6 +74,23 @@ WriterModule::WriterModule()
 
                               ));
 
+  this->_addBuiltinFunction("gnuplot_raw_writer", std::function(
+
+                                                    [](const std::string& filename) -> std::shared_ptr<const IWriter> {
+                                                      return std::make_shared<GnuplotWriterRaw>(filename);
+                                                    }
+
+                                                    ));
+
+  this->_addBuiltinFunction("gnuplot_raw_writer",
+                            std::function(
+
+                              [](const std::string& filename, const double& period) -> std::shared_ptr<const IWriter> {
+                                return std::make_shared<GnuplotWriterRaw>(filename, period);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("name_output", std::function(
 
                                              [](std::shared_ptr<const DiscreteFunctionVariant> discrete_function,
diff --git a/src/main.cpp b/src/main.cpp
index a9bec44f0e7c2912ec4e8a4cd8823b03928b5bbf..917addbe58d9fd90d6a5316315dc888f0a6bb97c 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -17,6 +17,8 @@ main(int argc, char* argv[])
   ExecutionStatManager::create();
   ParallelChecker::create();
 
+  GlobalVariableManager::create();
+
   std::string filename = initialize(argc, argv);
 
   SynchronizerManager::create();
@@ -27,13 +29,9 @@ main(int argc, char* argv[])
   DualMeshManager::create();
   StencilManager::create();
 
-  GlobalVariableManager::create();
-
   parser(filename);
   ExecutionStatManager::printInfo();
 
-  GlobalVariableManager::destroy();
-
   StencilManager::destroy();
   DualMeshManager::destroy();
   DualConnectivityManager::destroy();
@@ -44,6 +42,8 @@ main(int argc, char* argv[])
 
   finalize();
 
+  GlobalVariableManager::destroy();
+
   ParallelChecker::destroy();
 
   int return_code = ExecutionStatManager::getInstance().exitCode();
diff --git a/src/mesh/ConnectivityDispatcher.cpp b/src/mesh/ConnectivityDispatcher.cpp
index 65326b56bd49c1ce155dbf07fe32da8622c2e66c..e9491ad34c928b5c84e2f090057a46cebfaffba8 100644
--- a/src/mesh/ConnectivityDispatcher.cpp
+++ b/src/mesh/ConnectivityDispatcher.cpp
@@ -2,6 +2,7 @@
 
 #include <mesh/ItemOfItemType.hpp>
 #include <utils/CRSGraph.hpp>
+#include <utils/GlobalVariableManager.hpp>
 #include <utils/Partitioner.hpp>
 
 #include <iostream>
@@ -28,17 +29,17 @@ ConnectivityDispatcher<Dimension>::_buildNewOwner()
     using ItemId = ItemIdT<item_type>;
     ItemValue<int, item_type> item_new_owner(m_connectivity);
     parallel_for(
-      item_new_owner.numberOfItems(), PUGS_LAMBDA(const ItemId& l) {
-        const auto& item_to_cell = item_to_cell_matrix[l];
-        CellId Jmin              = item_to_cell[0];
-
-        for (size_t j = 1; j < item_to_cell.size(); ++j) {
-          const CellId J = item_to_cell[j];
-          if (cell_number[J] < cell_number[Jmin]) {
-            Jmin = J;
+      item_new_owner.numberOfItems(), PUGS_LAMBDA(const ItemId& item_id) {
+        const auto& item_to_cell  = item_to_cell_matrix[item_id];
+        CellId min_number_cell_id = item_to_cell[0];
+
+        for (size_t i_cell = 1; i_cell < item_to_cell.size(); ++i_cell) {
+          const CellId cell_id = item_to_cell[i_cell];
+          if (cell_number[cell_id] < cell_number[min_number_cell_id]) {
+            min_number_cell_id = cell_id;
           }
         }
-        item_new_owner[l] = cell_new_owner[Jmin];
+        item_new_owner[item_id] = cell_new_owner[min_number_cell_id];
       });
 
     synchronize(item_new_owner);
@@ -72,22 +73,65 @@ ConnectivityDispatcher<Dimension>::_buildItemListToSend()
 
     std::vector<std::vector<CellId>> cell_vector_to_send_by_proc(parallel::size());
     Array<bool> send_to_rank(parallel::size());
-    for (CellId j = 0; j < m_connectivity.numberOfCells(); ++j) {
+
+    NodeValue<bool> node_tag{m_connectivity};
+    node_tag.fill(false);
+    std::vector<NodeId> layer_node_id_list;
+    std::vector<NodeId> tagged_node_id_list;
+
+    const size_t nb_layers = GlobalVariableManager::instance().getNumberOfGhostLayers();
+
+    for (CellId cell_id = 0; cell_id < m_connectivity.numberOfCells(); ++cell_id) {
+      layer_node_id_list.clear();
       send_to_rank.fill(false);
-      const auto& cell_to_node = cell_to_node_matrix[j];
-
-      for (size_t R = 0; R < cell_to_node.size(); ++R) {
-        const NodeId& r          = cell_to_node[R];
-        const auto& node_to_cell = node_to_cell_matrix[r];
-        for (size_t K = 0; K < node_to_cell.size(); ++K) {
-          const CellId& k                 = node_to_cell[K];
-          send_to_rank[cell_new_owner[k]] = true;
+      const auto& cell_to_node = cell_to_node_matrix[cell_id];
+      for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+        const NodeId node_id = cell_to_node[i_node];
+        layer_node_id_list.push_back(node_id);
+        node_tag[node_id] = true;
+      }
+
+      for (size_t i_layer = 0; i_layer < nb_layers; ++i_layer) {
+        for (const auto node_id : layer_node_id_list) {
+          tagged_node_id_list.push_back(node_id);
+          const auto& node_to_cell = node_to_cell_matrix[node_id];
+          for (size_t i_node_cell = 0; i_node_cell < node_to_cell.size(); ++i_node_cell) {
+            const CellId& node_cell_id                 = node_to_cell[i_node_cell];
+            send_to_rank[cell_new_owner[node_cell_id]] = true;
+          }
+        }
+
+        if (i_layer + 1 < nb_layers) {
+          std::vector<NodeId> old_layer_node_id_list;
+          std::swap(layer_node_id_list, old_layer_node_id_list);
+
+          for (const auto node_id : old_layer_node_id_list) {
+            const auto& node_to_cell = node_to_cell_matrix[node_id];
+            for (size_t i_node_cell = 0; i_node_cell < node_to_cell.size(); ++i_node_cell) {
+              const CellId& node_cell_id    = node_to_cell[i_node_cell];
+              const auto& node_cell_to_node = cell_to_node_matrix[node_cell_id];
+              for (size_t i_node = 0; i_node < node_cell_to_node.size(); ++i_node) {
+                const NodeId cell_node_id = node_cell_to_node[i_node];
+                if (not node_tag[cell_node_id]) {
+                  layer_node_id_list.push_back(cell_node_id);
+                  node_tag[cell_node_id] = true;
+                }
+              }
+            }
+          }
+        } else {
+          layer_node_id_list.clear();
         }
       }
 
-      for (size_t k = 0; k < send_to_rank.size(); ++k) {
-        if (send_to_rank[k]) {
-          cell_vector_to_send_by_proc[k].push_back(j);
+      for (auto node_id : tagged_node_id_list) {
+        node_tag[node_id] = false;
+      }
+      tagged_node_id_list.clear();
+
+      for (size_t i_rank = 0; i_rank < send_to_rank.size(); ++i_rank) {
+        if (send_to_rank[i_rank]) {
+          cell_vector_to_send_by_proc[i_rank].push_back(cell_id);
         }
       }
     }
@@ -110,11 +154,11 @@ ConnectivityDispatcher<Dimension>::_buildItemListToSend()
       Array<bool> tag(m_connectivity.template numberOf<item_type>());
       tag.fill(false);
       std::vector<ItemId> item_id_vector;
-      for (size_t j = 0; j < cell_list_to_send_by_proc[i_rank].size(); ++j) {
-        const CellId& cell_id          = cell_list_to_send_by_proc[i_rank][j];
+      for (size_t i_cell = 0; i_cell < cell_list_to_send_by_proc[i_rank].size(); ++i_cell) {
+        const CellId& cell_id          = cell_list_to_send_by_proc[i_rank][i_cell];
         const auto& cell_sub_item_list = cell_to_sub_item_matrix[cell_id];
-        for (size_t r = 0; r < cell_sub_item_list.size(); ++r) {
-          const ItemId& item_id = cell_sub_item_list[r];
+        for (size_t i_item = 0; i_item < cell_sub_item_list.size(); ++i_item) {
+          const ItemId& item_id = cell_sub_item_list[i_item];
           if (not tag[item_id]) {
             item_id_vector.push_back(item_id);
             tag[item_id] = true;
@@ -155,9 +199,9 @@ ConnectivityDispatcher<Dimension>::_gatherFrom(const ItemValue<DataType, item_ty
   gathered_array = Array<std::remove_const_t<DataType>>(this->_dispatchedInfo<item_type>().m_number_to_id_map.size());
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
     Assert(recv_id_correspondance_by_proc[i_rank].size() == recv_item_data_by_proc[i_rank].size());
-    for (size_t r = 0; r < recv_id_correspondance_by_proc[i_rank].size(); ++r) {
-      const auto& item_id     = recv_id_correspondance_by_proc[i_rank][r];
-      gathered_array[item_id] = recv_item_data_by_proc[i_rank][r];
+    for (size_t i_item = 0; i_item < recv_id_correspondance_by_proc[i_rank].size(); ++i_item) {
+      const auto& item_id     = recv_id_correspondance_by_proc[i_rank][i_item];
+      gathered_array[item_id] = recv_item_data_by_proc[i_rank][i_item];
     }
   }
 }
@@ -180,11 +224,11 @@ ConnectivityDispatcher<Dimension>::_gatherFrom(
 
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
     std::vector<MutableDataType> data_by_item_vector;
-    for (size_t j = 0; j < item_list_to_send_by_proc[i_rank].size(); ++j) {
-      const ItemId& item_id = item_list_to_send_by_proc[i_rank][j];
+    for (size_t i_item = 0; i_item < item_list_to_send_by_proc[i_rank].size(); ++i_item) {
+      const ItemId& item_id = item_list_to_send_by_proc[i_rank][i_item];
       const auto& item_data = data_to_gather.itemArray(item_id);
-      for (size_t l = 0; l < item_data.size(); ++l) {
-        data_by_item_vector.push_back(item_data[l]);
+      for (size_t i_data = 0; i_data < item_data.size(); ++i_data) {
+        data_by_item_vector.push_back(item_data[i_data]);
       }
     }
     data_to_send_by_proc[i_rank] = convert_to_array(data_by_item_vector);
@@ -211,13 +255,13 @@ ConnectivityDispatcher<Dimension>::_gatherFrom(
 
   gathered_array = Array<std::remove_const_t<DataType>>(recv_array_size);
   {
-    size_t l = 0;
+    size_t i_value = 0;
     for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-      for (size_t j = 0; j < recv_data_to_gather_by_proc[i_rank].size(); ++j) {
-        gathered_array[l++] = recv_data_to_gather_by_proc[i_rank][j];
+      for (size_t i_rank_value = 0; i_rank_value < recv_data_to_gather_by_proc[i_rank].size(); ++i_rank_value) {
+        gathered_array[i_value++] = recv_data_to_gather_by_proc[i_rank][i_rank_value];
       }
     }
-    Assert(gathered_array.size() == l);
+    Assert(gathered_array.size() == i_value);
   }
 }
 
@@ -229,8 +273,8 @@ ConnectivityDispatcher<Dimension>::_buildCellNumberIdMap()
   auto& cell_number_id_map            = this->_dispatchedInfo<ItemType::cell>().m_number_to_id_map;
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
     CellId cell_id = 0;
-    for (size_t i = 0; i < recv_cell_number_by_proc[i_rank].size(); ++i) {
-      const int cell_number     = recv_cell_number_by_proc[i_rank][i];
+    for (size_t i_rank_cell = 0; i_rank_cell < recv_cell_number_by_proc[i_rank].size(); ++i_rank_cell) {
+      const int cell_number     = recv_cell_number_by_proc[i_rank][i_rank_cell];
       auto [iterator, inserted] = cell_number_id_map.insert(std::make_pair(cell_number, cell_id));
       if (inserted)
         ++cell_id;
@@ -252,8 +296,9 @@ ConnectivityDispatcher<Dimension>::_buildSubItemNumberToIdMap()
   auto& sub_item_number_id_map = this->_dispatchedInfo<ItemOfItemT::sub_item_type>().m_number_to_id_map;
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
     int sub_item_id = 0;
-    for (size_t i = 0; i < cell_sub_item_number_to_recv_by_proc[i_rank].size(); ++i) {
-      int sub_item_number       = cell_sub_item_number_to_recv_by_proc[i_rank][i];
+    for (size_t i_rank_sub_item = 0; i_rank_sub_item < cell_sub_item_number_to_recv_by_proc[i_rank].size();
+         ++i_rank_sub_item) {
+      int sub_item_number       = cell_sub_item_number_to_recv_by_proc[i_rank][i_rank_sub_item];
       auto [iterator, inserted] = sub_item_number_id_map.insert(std::make_pair(sub_item_number, sub_item_id));
       if (inserted)
         sub_item_id++;
@@ -273,8 +318,9 @@ ConnectivityDispatcher<Dimension>::_buildNumberOfSubItemPerItemToRecvByProc()
 
   using ItemId = ItemIdT<SubItemOfItemT::item_type>;
   parallel_for(
-    number_of_sub_item_per_item.numberOfItems(),
-    PUGS_LAMBDA(const ItemId& j) { number_of_sub_item_per_item[j] = item_to_sub_item_matrix[j].size(); });
+    number_of_sub_item_per_item.numberOfItems(), PUGS_LAMBDA(const ItemId& item_id) {
+      number_of_sub_item_per_item[item_id] = item_to_sub_item_matrix[item_id].size();
+    });
 
   this->_dispatchedInfo<SubItemOfItemT>().m_number_of_sub_item_per_item_to_recv_by_proc =
     this->exchange(number_of_sub_item_per_item);
@@ -299,11 +345,11 @@ ConnectivityDispatcher<Dimension>::_buildSubItemNumbersToRecvByProc()
       const auto& item_list_to_send_by_proc = this->_dispatchedInfo<SubItemOfItemT::item_type>().m_list_to_send_by_proc;
 
       std::vector<int> sub_item_numbers_by_item_vector;
-      for (size_t j = 0; j < item_list_to_send_by_proc[i_rank].size(); ++j) {
-        const ItemId& item_id     = item_list_to_send_by_proc[i_rank][j];
+      for (size_t i_rank_item = 0; i_rank_item < item_list_to_send_by_proc[i_rank].size(); ++i_rank_item) {
+        const ItemId& item_id     = item_list_to_send_by_proc[i_rank][i_rank_item];
         const auto& sub_item_list = item_to_sub_item_matrix[item_id];
-        for (size_t r = 0; r < sub_item_list.size(); ++r) {
-          const SubItemId& sub_item_id = sub_item_list[r];
+        for (size_t i_sub_item = 0; i_sub_item < sub_item_list.size(); ++i_sub_item) {
+          const SubItemId& sub_item_id = sub_item_list[i_sub_item];
           sub_item_numbers_by_item_vector.push_back(sub_item_number[sub_item_id]);
         }
       }
@@ -351,11 +397,14 @@ ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
   std::vector<std::vector<unsigned int>> item_to_subitem_legacy;
   size_t number_of_node_by_cell = 0;
   for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-    int l = 0;
-    for (size_t i = 0; i < item_list_to_recv_size_by_proc[i_rank]; ++i) {
+    int i_sub_item = 0;
+    for (size_t i_rank_item = 0; i_rank_item < item_list_to_recv_size_by_proc[i_rank]; ++i_rank_item) {
       std::vector<unsigned int> sub_item_vector;
-      for (int k = 0; k < number_of_sub_item_per_item_to_recv_by_proc[i_rank][i]; ++k) {
-        const auto& searched_sub_item_id = sub_item_number_id_map.find(recv_item_of_item_numbers_by_proc[i_rank][l++]);
+      for (int i_rank_sub_item_per_item = 0;
+           i_rank_sub_item_per_item < number_of_sub_item_per_item_to_recv_by_proc[i_rank][i_rank_item];
+           ++i_rank_sub_item_per_item) {
+        const auto& searched_sub_item_id =
+          sub_item_number_id_map.find(recv_item_of_item_numbers_by_proc[i_rank][i_sub_item++]);
         Assert(searched_sub_item_id != sub_item_number_id_map.end());
         sub_item_vector.push_back(searched_sub_item_id->second);
       }
@@ -371,14 +420,16 @@ ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
   item_to_subitem_list.fill(10000000);
 
   item_to_subitem_row_map[0] = 0;
-  for (size_t i = 0; i < item_to_subitem_legacy.size(); ++i) {
-    item_to_subitem_row_map[i + 1] = item_to_subitem_row_map[i] + item_to_subitem_legacy[i].size();
+  for (size_t i_rank = 0; i_rank < item_to_subitem_legacy.size(); ++i_rank) {
+    item_to_subitem_row_map[i_rank + 1] = item_to_subitem_row_map[i_rank] + item_to_subitem_legacy[i_rank].size();
   }
-  size_t l = 0;
-  for (size_t i = 0; i < item_to_subitem_legacy.size(); ++i) {
-    const auto& subitem_list = item_to_subitem_legacy[i];
-    for (size_t j = 0; j < subitem_list.size(); ++j, ++l) {
-      item_to_subitem_list[l] = subitem_list[j];
+  {
+    size_t i_sub_item = 0;
+    for (size_t i_rank = 0; i_rank < item_to_subitem_legacy.size(); ++i_rank) {
+      const auto& subitem_list = item_to_subitem_legacy[i_rank];
+      for (size_t i_rank_sub_item = 0; i_rank_sub_item < subitem_list.size(); ++i_rank_sub_item, ++i_sub_item) {
+        item_to_subitem_list[i_sub_item] = subitem_list[i_rank_sub_item];
+      }
     }
   }
 
@@ -401,7 +452,8 @@ ConnectivityDispatcher<Dimension>::_buildRecvItemIdCorrespondanceByProc()
     Array<int> send_item_number(item_list_to_send_by_proc[i_rank].size());
     const Array<const ItemId> send_item_id = item_list_to_send_by_proc[i_rank];
     parallel_for(
-      send_item_number.size(), PUGS_LAMBDA(size_t j) { send_item_number[j] = item_number[send_item_id[j]]; });
+      send_item_number.size(),
+      PUGS_LAMBDA(size_t i_item) { send_item_number[i_item] = item_number[send_item_id[i_item]]; });
     send_item_number_by_proc[i_rank] = send_item_number;
   }
 
@@ -415,11 +467,11 @@ ConnectivityDispatcher<Dimension>::_buildRecvItemIdCorrespondanceByProc()
   const auto& item_number_to_id_map = this->_dispatchedInfo<item_type>().m_number_to_id_map;
   for (size_t i_rank = 0; i_rank < item_list_to_recv_size_by_proc.size(); ++i_rank) {
     Array<ItemId> item_id_correspondance(item_list_to_recv_size_by_proc[i_rank]);
-    for (size_t l = 0; l < item_list_to_recv_size_by_proc[i_rank]; ++l) {
-      const int& recv_item_number  = recv_item_number_by_proc[i_rank][l];
+    for (size_t i_rank_item = 0; i_rank_item < item_list_to_recv_size_by_proc[i_rank]; ++i_rank_item) {
+      const int& recv_item_number  = recv_item_number_by_proc[i_rank][i_rank_item];
       const auto& searched_item_id = item_number_to_id_map.find(recv_item_number);
       Assert(searched_item_id != item_number_to_id_map.end());
-      item_id_correspondance[l] = searched_item_id->second;
+      item_id_correspondance[i_rank_item] = searched_item_id->second;
     }
     recv_item_id_correspondance_by_proc[i_rank] = item_id_correspondance;
   }
@@ -572,9 +624,9 @@ ConnectivityDispatcher<Dimension>::_buildItemReferenceList()
           this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
         std::vector<block_type> item_refs(m_new_descriptor.template itemNumberVector<item_type>().size());
         for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-          for (size_t r = 0; r < recv_item_refs_by_proc[i_rank].size(); ++r) {
-            const ItemId& item_id = recv_item_id_correspondance_by_proc[i_rank][r];
-            item_refs[item_id]    = recv_item_refs_by_proc[i_rank][r];
+          for (size_t i_item = 0; i_item < recv_item_refs_by_proc[i_rank].size(); ++i_item) {
+            const ItemId& item_id = recv_item_id_correspondance_by_proc[i_rank][i_item];
+            item_refs[item_id]    = recv_item_refs_by_proc[i_rank][i_item];
           }
         }
 
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 900c1b4d2e5dad1ff338c7c7714fcbe5f498cfb7..0000000000000000000000000000000000000000
--- a/src/mesh/Stencil.hpp
+++ /dev/null
@@ -1,44 +0,0 @@
-#ifndef STENCIL_HPP
-#define STENCIL_HPP
-
-#include <mesh/ConnectivityMatrix.hpp>
-#include <mesh/IBoundaryDescriptor.hpp>
-#include <mesh/ItemId.hpp>
-
-class Stencil
-{
- public:
-  using BoundaryDescriptorStencilList =
-    std::vector<std::pair<std::shared_ptr<const IBoundaryDescriptor>, ConnectivityMatrix>>;
-
- 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..9571c346bae6e4e345a16b09188dce15ab73b456
--- /dev/null
+++ b/src/mesh/StencilArray.hpp
@@ -0,0 +1,105 @@
+#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>;
+using CellToFaceStencilArray = StencilArray<ItemType::cell, ItemType::face>;
+using NodeToCellStencilArray = StencilArray<ItemType::node, ItemType::cell>;
+
+#endif   // STENCIL_ARRAY_HPP
diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp
index ddff37e0cb6a693999701b83eb5418ef35644af8..31db762a1002e3c633afeaabd4f06d6a98dff634 100644
--- a/src/mesh/StencilBuilder.cpp
+++ b/src/mesh/StencilBuilder.cpp
@@ -2,10 +2,189 @@
 
 #include <mesh/Connectivity.hpp>
 #include <mesh/ItemArray.hpp>
+#include <utils/GlobalVariableManager.hpp>
 #include <utils/Messenger.hpp>
 
 #include <set>
 
+template <ItemType item_type>
+class StencilBuilder::Layer
+{
+  using ItemId = ItemIdT<item_type>;
+  std::vector<ItemId> m_item_id_vector;
+  std::vector<int> m_item_number_vector;
+
+ public:
+  size_t
+  size() const
+  {
+    Assert(m_item_id_vector.size() == m_item_number_vector.size());
+    return m_item_id_vector.size();
+  }
+
+  bool
+  hasItemNumber(const int item_number) const
+  {
+    ssize_t begin = 0;
+    ssize_t end   = m_item_number_vector.size();
+
+    while (begin < end) {
+      const ssize_t mid      = (begin + end) / 2;
+      const auto& mid_number = m_item_number_vector[mid];
+      if (mid_number == item_number) {
+        return true;   // We found the value
+      } else if (mid_number < item_number) {
+        if (begin == mid) {
+          break;
+        }
+        begin = mid - 1;
+      } else {
+        if (end == mid) {
+          break;
+        }
+        end = mid + 1;
+      }
+    }
+    return false;
+  }
+
+  const std::vector<ItemId>&
+  itemIdList() const
+  {
+    return m_item_id_vector;
+  }
+
+  void
+  add(const ItemId item_id, const int item_number)
+  {
+    ssize_t begin = 0;
+    ssize_t end   = m_item_number_vector.size();
+
+    while (begin < end) {
+      const ssize_t mid      = (begin + end) / 2;
+      const auto& mid_number = m_item_number_vector[mid];
+
+      if (mid_number == item_number) {
+        return;   // We found the value
+      } else if (mid_number < item_number) {
+        if (begin == mid) {
+          break;
+        }
+        begin = mid;
+      } else {
+        if (end == mid) {
+          break;
+        }
+        end = mid;
+      }
+    }
+
+    m_item_id_vector.push_back(item_id);
+    m_item_number_vector.push_back(item_number);
+
+    const auto& begin_number = m_item_number_vector[begin];
+
+    if (begin_number > item_number) {
+      for (ssize_t i = m_item_number_vector.size() - 2; i >= begin; --i) {
+        std::swap(m_item_number_vector[i], m_item_number_vector[i + 1]);
+        std::swap(m_item_id_vector[i], m_item_id_vector[i + 1]);
+      }
+    } else if (begin_number < item_number) {
+      for (ssize_t i = m_item_number_vector.size() - 2; i > begin; --i) {
+        std::swap(m_item_number_vector[i], m_item_number_vector[i + 1]);
+        std::swap(m_item_id_vector[i], m_item_id_vector[i + 1]);
+      }
+    }
+  }
+
+  Layer() = default;
+
+  Layer(const Layer&) = default;
+  Layer(Layer&&)      = default;
+
+  ~Layer() = default;
+};
+
+template <ItemType connecting_item_type, typename ConnectivityType>
+auto
+StencilBuilder::_buildSymmetryConnectingItemList(const ConnectivityType& connectivity,
+                                                 const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
+{
+  static_assert(connecting_item_type != ItemType::cell, "cells cannot be used to define symmetry boundaries");
+
+  using ConnectingItemId = ItemIdT<connecting_item_type>;
+  ItemArray<bool, connecting_item_type> symmetry_connecting_item_list{connectivity,
+                                                                      symmetry_boundary_descriptor_list.size()};
+  symmetry_connecting_item_list.fill(false);
+
+  if constexpr (ConnectivityType::Dimension > 1) {
+    auto face_to_connecting_item_matrix =
+      connectivity.template getItemToItemMatrix<ItemType::face, connecting_item_type>();
+    size_t i_symmetry_boundary = 0;
+    for (auto p_boundary_descriptor : symmetry_boundary_descriptor_list) {
+      const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor;
+
+      bool found = false;
+      for (size_t i_ref_face_list = 0; i_ref_face_list < connectivity.template numberOfRefItemList<ItemType::face>();
+           ++i_ref_face_list) {
+        const auto& ref_face_list = connectivity.template refItemList<ItemType::face>(i_ref_face_list);
+        if (ref_face_list.refId() == boundary_descriptor) {
+          found = true;
+          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 connecting_item_list = face_to_connecting_item_matrix[face_id];
+            for (size_t i_connecting_item = 0; i_connecting_item < connecting_item_list.size(); ++i_connecting_item) {
+              const ConnectingItemId connecting_item_id = connecting_item_list[i_connecting_item];
+
+              symmetry_connecting_item_list[connecting_item_id][i_symmetry_boundary] = true;
+            }
+          }
+          break;
+        }
+      }
+      ++i_symmetry_boundary;
+      if (not found) {
+        std::ostringstream error_msg;
+        error_msg << "cannot find boundary '" << rang::fgB::yellow << boundary_descriptor << rang::fg::reset << '\'';
+        throw NormalError(error_msg.str());
+      }
+    }
+
+    return symmetry_connecting_item_list;
+  } else {
+    size_t i_symmetry_boundary = 0;
+    for (auto p_boundary_descriptor : symmetry_boundary_descriptor_list) {
+      const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor;
+
+      bool found = false;
+      for (size_t i_ref_connecting_item_list = 0;
+           i_ref_connecting_item_list < connectivity.template numberOfRefItemList<connecting_item_type>();
+           ++i_ref_connecting_item_list) {
+        const auto& ref_connecting_item_list =
+          connectivity.template refItemList<connecting_item_type>(i_ref_connecting_item_list);
+        if (ref_connecting_item_list.refId() == boundary_descriptor) {
+          found                     = true;
+          auto connecting_item_list = ref_connecting_item_list.list();
+          for (size_t i_connecting_item = 0; i_connecting_item < connecting_item_list.size(); ++i_connecting_item) {
+            const ConnectingItemId connecting_item_id = connecting_item_list[i_connecting_item];
+
+            symmetry_connecting_item_list[connecting_item_id][i_symmetry_boundary] = true;
+          }
+          break;
+        }
+      }
+      ++i_symmetry_boundary;
+      if (not found) {
+        std::ostringstream error_msg;
+        error_msg << "cannot find boundary '" << rang::fgB::yellow << boundary_descriptor << rang::fg::reset << '\'';
+        throw NormalError(error_msg.str());
+      }
+    }
+
+    return symmetry_connecting_item_list;
+  }
+}
+
 template <typename ConnectivityType>
 Array<const uint32_t>
 StencilBuilder::_getRowMap(const ConnectivityType& connectivity) const
@@ -79,8 +258,8 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar
               }
             }
             if (not found) {
-              int node_cell_number = cell_number[node_cell_id];
-              size_t i_index       = row_map[cell_id];
+              const auto node_cell_number = cell_number[node_cell_id];
+              size_t i_index              = row_map[cell_id];
               // search for position for index
               while ((i_index < max_index[cell_id])) {
                 if (node_cell_number > cell_number[CellId(column_indices[i_index])]) {
@@ -91,7 +270,7 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar
               }
 
               for (size_t i_destination = max_index[cell_id]; i_destination > i_index; --i_destination) {
-                size_t i_source = i_destination - 1;
+                const size_t i_source = i_destination - 1;
 
                 column_indices[i_destination] = column_indices[i_source];
               }
@@ -107,67 +286,118 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar
   return column_indices;
 }
 
-template <typename ConnectivityType>
-Stencil
-StencilBuilder::_build(const ConnectivityType& connectivity,
-                       size_t number_of_layers,
-                       const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
+template <ItemType item_type, ItemType connecting_item_type, typename ConnectivityType>
+StencilArray<item_type, item_type>
+StencilBuilder::_build_for_same_source_and_target(const ConnectivityType& connectivity,
+                                                  const size_t& number_of_layers,
+                                                  const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
 {
-#warning TEMPORARY
+  if (number_of_layers == 0) {
+    throw NormalError("number of layers must be greater than 0 to build stencils");
+  }
+  if (number_of_layers > 2) {
+    throw NotImplementedError("number of layers too large");
+  }
+  auto item_to_connecting_item_matrix = connectivity.template getItemToItemMatrix<item_type, connecting_item_type>();
+  auto connecting_item_to_item_matrix = connectivity.template getItemToItemMatrix<connecting_item_type, item_type>();
+
+  auto item_is_owned          = connectivity.template isOwned<item_type>();
+  auto item_number            = connectivity.template number<item_type>();
+  auto connecting_item_number = connectivity.template number<connecting_item_type>();
+
+  using ItemId           = ItemIdT<item_type>;
+  using ConnectingItemId = ItemIdT<connecting_item_type>;
+
   if (symmetry_boundary_descriptor_list.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);
+      Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1};
+      row_map[0] = 0;
 
-      return {ConnectivityMatrix{row_map, column_indices}, {}};
-    } else {
-      if (number_of_layers > 2) {
-        throw NotImplementedError("number of layers too large");
-      }
-      if (parallel::size() > 1) {
-        throw NotImplementedError("stencils with more than 1 layers are not defined in parallel");
+      std::vector<ItemId> column_indices_vector;
+
+      for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) {
+        // First layer is a special case
+
+        Layer<item_type> item_layer;
+        Layer<connecting_item_type> connecting_layer;
+
+        if (item_is_owned[item_id]) {
+          for (size_t i_connecting_item_1 = 0; i_connecting_item_1 < item_to_connecting_item_matrix[item_id].size();
+               ++i_connecting_item_1) {
+            const ConnectingItemId layer_1_connecting_item_id =
+              item_to_connecting_item_matrix[item_id][i_connecting_item_1];
+            connecting_layer.add(layer_1_connecting_item_id, connecting_item_number[layer_1_connecting_item_id]);
+          }
+
+          for (auto connecting_item_id : connecting_layer.itemIdList()) {
+            for (size_t i_item_1 = 0; i_item_1 < connecting_item_to_item_matrix[connecting_item_id].size();
+                 ++i_item_1) {
+              const ItemId layer_1_item_id = connecting_item_to_item_matrix[connecting_item_id][i_item_1];
+              if (layer_1_item_id != item_id) {
+                item_layer.add(layer_1_item_id, item_number[layer_1_item_id]);
+              }
+            }
+          }
+        }
+
+        for (auto layer_item_id : item_layer.itemIdList()) {
+          column_indices_vector.push_back(layer_item_id);
+        }
+        row_map[item_id + 1] = row_map[item_id] + item_layer.itemIdList().size();
       }
 
-      auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
-      auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+      if (row_map[row_map.size() - 1] != column_indices_vector.size()) {
+        throw UnexpectedError("incorrect stencil size");
+      }
+      Array<uint32_t> column_indices(row_map[row_map.size() - 1]);
+      column_indices.fill(std::numeric_limits<uint32_t>::max());
 
-      auto cell_is_owned = connectivity.cellIsOwned();
-      auto cell_number   = connectivity.cellNumber();
+      for (size_t i = 0; i < column_indices.size(); ++i) {
+        column_indices[i] = column_indices_vector[i];
+      }
 
-      Array<uint32_t> row_map{connectivity.numberOfCells() + 1};
+      return {ConnectivityMatrix{row_map, column_indices}, {}};
+    } else {
+      Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1};
       row_map[0] = 0;
 
-      std::vector<CellId> column_indices_vector;
+      std::vector<ItemId> column_indices_vector;
 
-      for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
-        if (cell_is_owned[cell_id]) {
-          std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
-            [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
+      for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) {
+        if (item_is_owned[item_id]) {
+          std::set<ItemId, std::function<bool(ItemId, ItemId)>> item_set(
+            [=](ItemId item_0, ItemId item_1) { return item_number[item_0] < item_number[item_1]; });
 
-          for (size_t i_node_1 = 0; i_node_1 < cell_to_node_matrix[cell_id].size(); ++i_node_1) {
-            const NodeId layer_1_node_id = cell_to_node_matrix[cell_id][i_node_1];
+          for (size_t i_connecting_item_1 = 0; i_connecting_item_1 < item_to_connecting_item_matrix[item_id].size();
+               ++i_connecting_item_1) {
+            const ConnectingItemId layer_1_connecting_item_id =
+              item_to_connecting_item_matrix[item_id][i_connecting_item_1];
 
-            for (size_t i_cell_1 = 0; i_cell_1 < node_to_cell_matrix[layer_1_node_id].size(); ++i_cell_1) {
-              CellId cell_1_id = node_to_cell_matrix[layer_1_node_id][i_cell_1];
+            for (size_t i_item_1 = 0; i_item_1 < connecting_item_to_item_matrix[layer_1_connecting_item_id].size();
+                 ++i_item_1) {
+              ItemId item_1_id = connecting_item_to_item_matrix[layer_1_connecting_item_id][i_item_1];
 
-              for (size_t i_node_2 = 0; i_node_2 < cell_to_node_matrix[cell_1_id].size(); ++i_node_2) {
-                const NodeId layer_2_node_id = cell_to_node_matrix[cell_1_id][i_node_2];
+              for (size_t i_connecting_item_2 = 0;
+                   i_connecting_item_2 < item_to_connecting_item_matrix[item_1_id].size(); ++i_connecting_item_2) {
+                const ConnectingItemId layer_2_connecting_item_id =
+                  item_to_connecting_item_matrix[item_1_id][i_connecting_item_2];
 
-                for (size_t i_cell_2 = 0; i_cell_2 < node_to_cell_matrix[layer_2_node_id].size(); ++i_cell_2) {
-                  CellId cell_2_id = node_to_cell_matrix[layer_2_node_id][i_cell_2];
+                for (size_t i_item_2 = 0; i_item_2 < connecting_item_to_item_matrix[layer_2_connecting_item_id].size();
+                     ++i_item_2) {
+                  ItemId item_2_id = connecting_item_to_item_matrix[layer_2_connecting_item_id][i_item_2];
 
-                  if (cell_2_id != cell_id) {
-                    cell_set.insert(cell_2_id);
+                  if (item_2_id != item_id) {
+                    item_set.insert(item_2_id);
                   }
                 }
               }
             }
           }
 
-          for (auto stencil_cell_id : cell_set) {
-            column_indices_vector.push_back(stencil_cell_id);
+          for (auto stencil_item_id : item_set) {
+            column_indices_vector.push_back(stencil_item_id);
           }
-          row_map[cell_id + 1] = row_map[cell_id] + cell_set.size();
+          row_map[item_id + 1] = row_map[item_id] + item_set.size();
         }
       }
 
@@ -186,173 +416,222 @@ StencilBuilder::_build(const ConnectivityType& connectivity,
       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_descriptor_list.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_symmetry_boundary = 0;
-        for (auto p_boundary_descriptor : symmetry_boundary_descriptor_list) {
-          const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor;
-
-          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() == boundary_descriptor) {
-              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_symmetry_boundary] = true;
-                }
-              }
-              break;
+    ItemArray<bool, connecting_item_type> symmetry_item_list =
+      this->_buildSymmetryConnectingItemList<connecting_item_type>(connectivity, symmetry_boundary_descriptor_list);
+
+    Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1};
+    row_map[0] = 0;
+    std::vector<Array<uint32_t>> symmetry_row_map_list(symmetry_boundary_descriptor_list.size());
+    for (auto&& symmetry_row_map : symmetry_row_map_list) {
+      symmetry_row_map    = Array<uint32_t>{connectivity.template numberOf<item_type>() + 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_descriptor_list.size());
+
+    for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) {
+      std::set<ItemId> item_set;
+      std::vector<std::set<ItemId>> by_boundary_symmetry_item(symmetry_boundary_descriptor_list.size());
+
+      if (item_is_owned[item_id]) {
+        auto item_to_connecting_item_list = item_to_connecting_item_matrix[item_id];
+        for (size_t i_connecting_item_of_item = 0; i_connecting_item_of_item < item_to_connecting_item_list.size();
+             ++i_connecting_item_of_item) {
+          const ConnectingItemId connecting_item_id_of_item = item_to_connecting_item_list[i_connecting_item_of_item];
+          auto connecting_item_to_item_list = connecting_item_to_item_matrix[connecting_item_id_of_item];
+          for (size_t i_item_of_connecting_item = 0; i_item_of_connecting_item < connecting_item_to_item_list.size();
+               ++i_item_of_connecting_item) {
+            const ItemId item_id_of_connecting_item = connecting_item_to_item_list[i_item_of_connecting_item];
+            if (item_id != item_id_of_connecting_item) {
+              item_set.insert(item_id_of_connecting_item);
             }
           }
-          ++i_symmetry_boundary;
-          if (not found) {
-            std::ostringstream error_msg;
-            error_msg << "cannot find boundary '" << rang::fgB::yellow << boundary_descriptor << rang::fg::reset
-                      << '\'';
-            throw NormalError(error_msg.str());
-          }
         }
-      }
-
-      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_descriptor_list.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_descriptor_list.size());
-
-      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_descriptor_list.size());
-
-        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);
-              }
-            }
+        {
+          std::vector<ItemId> item_vector;
+          for (auto&& set_item_id : item_set) {
+            item_vector.push_back(set_item_id);
           }
+          std::sort(item_vector.begin(), item_vector.end(),
+                    [&item_number](const ItemId& item0_id, const ItemId& item1_id) {
+                      return item_number[item0_id] < item_number[item1_id];
+                    });
 
-          {
-            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);
-            }
+          for (auto&& vector_item_id : item_vector) {
+            column_indices_vector.push_back(vector_item_id);
           }
+        }
 
-          for (size_t i = 0; i < symmetry_boundary_descriptor_list.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);
-                }
+        for (size_t i = 0; i < symmetry_boundary_descriptor_list.size(); ++i) {
+          std::set<ItemId> symmetry_item_set;
+          for (size_t i_connecting_item_of_item = 0; i_connecting_item_of_item < item_to_connecting_item_list.size();
+               ++i_connecting_item_of_item) {
+            const ConnectingItemId connecting_item_id_of_item = item_to_connecting_item_list[i_connecting_item_of_item];
+            if (symmetry_item_list[connecting_item_id_of_item][i]) {
+              auto item_of_connecting_item_list = connecting_item_to_item_matrix[connecting_item_id_of_item];
+              for (size_t i_item_of_connecting_item = 0;
+                   i_item_of_connecting_item < item_of_connecting_item_list.size(); ++i_item_of_connecting_item) {
+                const ItemId item_id_of_connecting_item = item_of_connecting_item_list[i_item_of_connecting_item];
+                symmetry_item_set.insert(item_id_of_connecting_item);
               }
             }
-            by_boundary_symmetry_cell[i] = symmetry_cell_set;
+          }
+          by_boundary_symmetry_item[i] = symmetry_item_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];
-                      });
+          std::vector<ItemId> item_vector;
+          for (auto&& set_item_id : symmetry_item_set) {
+            item_vector.push_back(set_item_id);
+          }
+          std::sort(item_vector.begin(), item_vector.end(),
+                    [&item_number](const ItemId& item0_id, const ItemId& item1_id) {
+                      return item_number[item0_id] < item_number[item1_id];
+                    });
 
-            for (auto&& vector_cell_id : cell_vector) {
-              symmetry_column_indices_vector[i].push_back(vector_cell_id);
-            }
+          for (auto&& vector_item_id : item_vector) {
+            symmetry_column_indices_vector[i].push_back(vector_item_id);
           }
         }
-        row_map[cell_id + 1] = row_map[cell_id] + cell_set.size();
+      }
+      row_map[item_id + 1] = row_map[item_id] + item_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();
-        }
+      for (size_t i = 0; i < symmetry_row_map_list.size(); ++i) {
+        symmetry_row_map_list[i][item_id + 1] = symmetry_row_map_list[i][item_id] + by_boundary_symmetry_item[i].size();
       }
-      ConnectivityMatrix primal_stencil{row_map, convert_to_array(column_indices_vector)};
-
-      Stencil::BoundaryDescriptorStencilList symmetry_boundary_stencil_list;
-      {
-        size_t i = 0;
-        for (auto&& p_boundary_descriptor : symmetry_boundary_descriptor_list) {
-          symmetry_boundary_stencil_list.emplace_back(
-            std::make_pair(p_boundary_descriptor,
-                           ConnectivityMatrix{symmetry_row_map_list[i],
-                                              convert_to_array(symmetry_column_indices_vector[i])}));
-          ++i;
-        }
+    }
+    ConnectivityMatrix primal_stencil{row_map, convert_to_array(column_indices_vector)};
+
+    typename StencilArray<item_type, item_type>::BoundaryDescriptorStencilArrayList symmetry_boundary_stencil_list;
+    {
+      size_t i = 0;
+      for (auto&& p_boundary_descriptor : symmetry_boundary_descriptor_list) {
+        symmetry_boundary_stencil_list.emplace_back(
+          typename StencilArray<item_type, item_type>::
+            BoundaryDescriptorStencilArray{p_boundary_descriptor,
+                                           ConnectivityMatrix{symmetry_row_map_list[i],
+                                                              convert_to_array(symmetry_column_indices_vector[i])}});
+        ++i;
       }
+    }
+
+    return {{primal_stencil}, {symmetry_boundary_stencil_list}};
+  }
+}
 
-      return {{primal_stencil}, {symmetry_boundary_stencil_list}};
+template <ItemType source_item_type,
+          ItemType connecting_item_type,
+          ItemType target_item_type,
+          typename ConnectivityType>
+StencilArray<source_item_type, target_item_type>
+StencilBuilder::_build_for_different_source_and_target(
+  const ConnectivityType& connectivity,
+  const size_t& number_of_layers,
+  const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
+{
+  static_assert(source_item_type != target_item_type);
+  throw NotImplementedError("different source target");
+}
 
+template <ItemType source_item_type,
+          ItemType connecting_item_type,
+          ItemType target_item_type,
+          typename ConnectivityType>
+StencilArray<source_item_type, target_item_type>
+StencilBuilder::_build(const ConnectivityType& connectivity,
+                       const size_t& number_of_layers,
+                       const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
+{
+  if constexpr (connecting_item_type != target_item_type) {
+    if constexpr (source_item_type == target_item_type) {
+      return this
+        ->_build_for_same_source_and_target<source_item_type, connecting_item_type>(connectivity, number_of_layers,
+                                                                                    symmetry_boundary_descriptor_list);
     } else {
-      throw NotImplementedError("Only implemented in 2D/3D");
+      return this->_build_for_different_source_and_target<source_item_type, connecting_item_type,
+                                                          target_item_type>(connectivity, number_of_layers,
+                                                                            symmetry_boundary_descriptor_list);
     }
+  } else {
+    std::ostringstream error_msg;
+    error_msg << "cannot build stencil of " << rang::fgB::yellow << itemName(target_item_type) << rang::fg::reset
+              << " using " << rang::fgB::yellow << itemName(connecting_item_type) << rang::fg::reset
+              << " for connectivity";
+    throw UnexpectedError(error_msg.str());
+  }
+}
+
+template <ItemType source_item_type, ItemType target_item_type, typename ConnectivityType>
+StencilArray<source_item_type, target_item_type>
+StencilBuilder::_build(const ConnectivityType& connectivity,
+                       const StencilDescriptor& stencil_descriptor,
+                       const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
+{
+  switch (stencil_descriptor.connectionType()) {
+  case StencilDescriptor::ConnectionType::by_nodes: {
+    return this->_build<source_item_type, ItemType::node, target_item_type>(connectivity,
+                                                                            stencil_descriptor.numberOfLayers(),
+                                                                            symmetry_boundary_descriptor_list);
+  }
+  case StencilDescriptor::ConnectionType::by_edges: {
+    return this->_build<source_item_type, ItemType::edge, target_item_type>(connectivity,
+                                                                            stencil_descriptor.numberOfLayers(),
+                                                                            symmetry_boundary_descriptor_list);
+  }
+  case StencilDescriptor::ConnectionType::by_faces: {
+    return this->_build<source_item_type, ItemType::face, target_item_type>(connectivity,
+                                                                            stencil_descriptor.numberOfLayers(),
+                                                                            symmetry_boundary_descriptor_list);
+  }
+  case StencilDescriptor::ConnectionType::by_cells: {
+    return this->_build<source_item_type, ItemType::cell, target_item_type>(connectivity,
+                                                                            stencil_descriptor.numberOfLayers(),
+                                                                            symmetry_boundary_descriptor_list);
+  }
+    // LCOV_EXCL_START
+  default: {
+    throw UnexpectedError("invalid connection type");
+  }
+    // LCOV_EXCL_STOP
   }
 }
 
-Stencil
-StencilBuilder::build(const IConnectivity& connectivity,
-                      size_t number_of_layers,
-                      const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
+CellToCellStencilArray
+StencilBuilder::buildC2C(const IConnectivity& connectivity,
+                         const StencilDescriptor& stencil_descriptor,
+                         const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
 {
+  if ((parallel::size() > 1) and
+      (stencil_descriptor.numberOfLayers() > GlobalVariableManager::instance().getNumberOfGhostLayers())) {
+    std::ostringstream error_msg;
+    error_msg << "Stencil builder requires" << rang::fgB::yellow << stencil_descriptor.numberOfLayers()
+              << rang::fg::reset << " layers while parallel number of ghost layer is "
+              << GlobalVariableManager::instance().getNumberOfGhostLayers() << ".\n";
+    error_msg << "Increase the number of ghost layers (using the '--number-of-ghost-layers' option).";
+    throw NormalError(error_msg.str());
+  }
+
   switch (connectivity.dimension()) {
   case 1: {
-    return StencilBuilder::_build(dynamic_cast<const Connectivity<1>&>(connectivity), number_of_layers,
-                                  symmetry_boundary_descriptor_list);
+    return this->_build<ItemType::cell, ItemType::cell>(dynamic_cast<const Connectivity<1>&>(connectivity),
+                                                        stencil_descriptor, symmetry_boundary_descriptor_list);
   }
   case 2: {
-    return StencilBuilder::_build(dynamic_cast<const Connectivity<2>&>(connectivity), number_of_layers,
-                                  symmetry_boundary_descriptor_list);
+    return this->_build<ItemType::cell, ItemType::cell>(dynamic_cast<const Connectivity<2>&>(connectivity),
+                                                        stencil_descriptor, symmetry_boundary_descriptor_list);
   }
   case 3: {
-    return StencilBuilder::_build(dynamic_cast<const Connectivity<3>&>(connectivity), number_of_layers,
-                                  symmetry_boundary_descriptor_list);
+    return this->_build<ItemType::cell, ItemType::cell>(dynamic_cast<const Connectivity<3>&>(connectivity),
+                                                        stencil_descriptor, symmetry_boundary_descriptor_list);
   }
   default: {
     throw UnexpectedError("invalid connectivity dimension");
   }
   }
 }
+
+NodeToCellStencilArray
+StencilBuilder::buildN2C(const IConnectivity&, const StencilDescriptor&, const BoundaryDescriptorList&) const
+{
+  throw NotImplementedError("node to cell stencil");
+}
diff --git a/src/mesh/StencilBuilder.hpp b/src/mesh/StencilBuilder.hpp
index 7f11c0686a49adebeb4c7ffd5a0c945c6b6fa07d..4f7ead00b618db999dd5386258ac058b664c65a7 100644
--- a/src/mesh/StencilBuilder.hpp
+++ b/src/mesh/StencilBuilder.hpp
@@ -2,18 +2,23 @@
 #define STENCIL_BUILDER_HPP
 
 #include <mesh/IBoundaryDescriptor.hpp>
-#include <mesh/Stencil.hpp>
+#include <mesh/StencilArray.hpp>
+#include <mesh/StencilDescriptor.hpp>
 
 #include <string>
 #include <vector>
 
 class IConnectivity;
+
 class StencilBuilder
 {
  public:
   using BoundaryDescriptorList = std::vector<std::shared_ptr<const IBoundaryDescriptor>>;
 
  private:
+  template <ItemType item_type>
+  class Layer;
+
   template <typename ConnectivityType>
   Array<const uint32_t> _getRowMap(const ConnectivityType& connectivity) const;
 
@@ -21,15 +26,68 @@ class StencilBuilder
   Array<const uint32_t> _getColumnIndices(const ConnectivityType& connectivity,
                                           const Array<const uint32_t>& row_map) const;
 
+  template <ItemType connecting_item, typename ConnectivityType>
+  auto _buildSymmetryConnectingItemList(const ConnectivityType& connectivity,
+                                        const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
+
+  template <ItemType item_type, ItemType connecting_item_type, typename ConnectivityType>
+  StencilArray<item_type, item_type> _build_for_same_source_and_target(
+    const ConnectivityType& connectivity,
+    const size_t& number_of_layers,
+    const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
+
+  template <ItemType source_item_type,
+            ItemType connecting_item_type,
+            ItemType target_item_type,
+            typename ConnectivityType>
+  StencilArray<source_item_type, target_item_type> _build_for_different_source_and_target(
+    const ConnectivityType& connectivity,
+    const size_t& number_of_layers,
+    const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
+
+  template <ItemType source_item_type,
+            ItemType connecting_item_type,
+            ItemType target_item_type,
+            typename ConnectivityType>
+  StencilArray<source_item_type, target_item_type> _build(
+    const ConnectivityType& connectivity,
+    const size_t& number_of_layers,
+    const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
+
+  template <ItemType source_item_type, ItemType target_item_type, typename ConnectivityType>
+  StencilArray<source_item_type, target_item_type> _build(
+    const ConnectivityType& connectivity,
+    const StencilDescriptor& stencil_descriptor,
+    const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
+
   template <typename ConnectivityType>
-  Stencil _build(const ConnectivityType& connectivity,
-                 size_t number_of_layers,
-                 const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
+  NodeToCellStencilArray _buildN2C(const ConnectivityType& connectivity,
+                                   size_t number_of_layers,
+                                   const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
+
+  CellToCellStencilArray buildC2C(const IConnectivity& connectivity,
+                                  const StencilDescriptor& stencil_descriptor,
+                                  const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const;
+
+  NodeToCellStencilArray buildN2C(const IConnectivity& connectivity,
+                                  const StencilDescriptor& stencil_descriptor,
+                                  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;
+  template <ItemType source_item_type, ItemType target_item_type>
+  StencilArray<source_item_type, target_item_type>
+  build(const IConnectivity& connectivity,
+        const StencilDescriptor& stencil_descriptor,
+        const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const
+  {
+    if constexpr ((source_item_type == ItemType::cell) and (target_item_type == ItemType::cell)) {
+      return buildC2C(connectivity, stencil_descriptor, symmetry_boundary_descriptor_list);
+    } else if constexpr ((source_item_type == ItemType::node) and (target_item_type == ItemType::cell)) {
+      return buildN2C(connectivity, stencil_descriptor, symmetry_boundary_descriptor_list);
+    } else {
+      static_assert(is_false_item_type_v<source_item_type>, "invalid stencil type");
+    }
+  }
 
  public:
   StencilBuilder()                      = default;
diff --git a/src/mesh/StencilDescriptor.hpp b/src/mesh/StencilDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d099f00864ea7f2485664fc4ed6928a69bf95bfc
--- /dev/null
+++ b/src/mesh/StencilDescriptor.hpp
@@ -0,0 +1,62 @@
+#ifndef STENCIL_DESCRIPTOR_HPP
+#define STENCIL_DESCRIPTOR_HPP
+
+#include <utils/PugsMacros.hpp>
+
+#include <cstddef>
+
+class StencilDescriptor
+{
+ public:
+  enum class ConnectionType : int
+  {
+    by_nodes = 1,
+    by_edges = 2,
+    by_faces = 3,
+    by_cells = 4
+  };
+
+ private:
+  size_t m_number_of_layers;
+  ConnectionType m_connection_type;
+
+ public:
+  PUGS_INLINE
+  const size_t&
+  numberOfLayers() const
+  {
+    return m_number_of_layers;
+  };
+
+  PUGS_INLINE
+  const ConnectionType&
+  connectionType() const
+  {
+    return m_connection_type;
+  };
+
+  PUGS_INLINE
+  friend bool
+  operator==(const StencilDescriptor& sd0, const StencilDescriptor& sd1)
+  {
+    return (sd0.m_number_of_layers == sd1.m_number_of_layers) and (sd0.m_connection_type == sd1.m_connection_type);
+  }
+
+  PUGS_INLINE
+  friend bool
+  operator!=(const StencilDescriptor& sd0, const StencilDescriptor& sd1)
+  {
+    return not(sd0 == sd1);
+  }
+
+  StencilDescriptor(const size_t number_of_layers, const ConnectionType type)
+    : m_number_of_layers{number_of_layers}, m_connection_type{type}
+  {}
+
+  StencilDescriptor(const StencilDescriptor&) = default;
+  StencilDescriptor(StencilDescriptor&&)      = default;
+
+  ~StencilDescriptor() = default;
+};
+
+#endif   // STENCIL_DESCRIPTOR_HPP
diff --git a/src/mesh/StencilManager.cpp b/src/mesh/StencilManager.cpp
index cc99577f23a5d9b85e344d39a1dc0e2260dcb22a..842f8858f9d1a571e28d3f0397f70e97682da0af 100644
--- a/src/mesh/StencilManager.cpp
+++ b/src/mesh/StencilManager.cpp
@@ -18,44 +18,82 @@ StencilManager::destroy()
   Assert(m_instance != nullptr, "StencilManager was not created");
 
   // LCOV_EXCL_START
-  if (m_instance->m_stored_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) {
-      error << " - connectivity of id " << rang::fgB::magenta << key.connectivity_id << rang::style::reset << '\n';
+  auto check_still_registered = [](auto& stored_stencil_map) {
+    if (stored_stencil_map.size() > 0) {
+      std::stringstream error;
+      error << ": some connectivities are still registered\n";
+      for (const auto& [key, stencil_p] : stored_stencil_map) {
+        error << " - connectivity of id " << rang::fgB::magenta << key.connectivity_id << rang::style::reset << '\n';
+      }
+      throw UnexpectedError(error.str());
     }
-    throw UnexpectedError(error.str());
-    // LCOV_EXCL_STOP
-  }
+  };
+
+  check_still_registered(m_instance->m_stored_cell_to_cell_stencil_map);
+  check_still_registered(m_instance->m_stored_node_to_cell_stencil_map);
+  // LCOV_EXCL_STOP
+
   delete m_instance;
   m_instance = nullptr;
 }
 
-const Stencil&
-StencilManager::getStencil(const IConnectivity& connectivity,
-                           size_t degree,
-                           const BoundaryDescriptorList& symmetry_boundary_descriptor_list)
+template <ItemType source_item_type, ItemType target_item_type>
+const StencilArray<source_item_type, target_item_type>&
+StencilManager::_getStencilArray(
+  const IConnectivity& connectivity,
+  const StencilDescriptor& stencil_descriptor,
+  const BoundaryDescriptorList& symmetry_boundary_descriptor_list,
+  StoredStencilTMap<source_item_type, target_item_type>& stored_source_to_target_stencil_map)
 {
-  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 stored_source_to_target_stencil_map.contains(
+        Key{connectivity.id(), stencil_descriptor, symmetry_boundary_descriptor_list})) {
+    stored_source_to_target_stencil_map[Key{connectivity.id(), stencil_descriptor, symmetry_boundary_descriptor_list}] =
+      std::make_shared<StencilArray<source_item_type, target_item_type>>(
+        StencilBuilder{}.template build<source_item_type, target_item_type>(connectivity, stencil_descriptor,
+                                                                            symmetry_boundary_descriptor_list));
   }
 
-  return *m_stored_stencil_map.at(Key{connectivity.id(), degree, symmetry_boundary_descriptor_list});
+  return *stored_source_to_target_stencil_map.at(
+    Key{connectivity.id(), stencil_descriptor, symmetry_boundary_descriptor_list});
+}
+
+const CellToCellStencilArray&
+StencilManager::getCellToCellStencilArray(const IConnectivity& connectivity,
+                                          const StencilDescriptor& stencil_descriptor,
+                                          const BoundaryDescriptorList& symmetry_boundary_descriptor_list)
+{
+  return this->_getStencilArray<ItemType::cell, ItemType::cell>(connectivity, stencil_descriptor,
+                                                                symmetry_boundary_descriptor_list,
+                                                                m_stored_cell_to_cell_stencil_map);
+}
+
+const NodeToCellStencilArray&
+StencilManager::getNodeToCellStencilArray(const IConnectivity& connectivity,
+                                          const StencilDescriptor& stencil_descriptor,
+                                          const BoundaryDescriptorList& symmetry_boundary_descriptor_list)
+{
+  return this->_getStencilArray<ItemType::node, ItemType::cell>(connectivity, stencil_descriptor,
+                                                                symmetry_boundary_descriptor_list,
+                                                                m_stored_node_to_cell_stencil_map);
 }
 
 void
 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) {
-      if (connectivity_id == key.connectivity_id) {
-        m_stored_stencil_map.erase(key);
-        has_removed = true;
-        break;
+  auto delete_connectivity_stencil = [&connectivity_id](auto& stored_stencil_map) {
+    bool has_removed = false;
+    do {
+      has_removed = false;
+      for (const auto& [key, p_stencil] : stored_stencil_map) {
+        if (connectivity_id == key.connectivity_id) {
+          stored_stencil_map.erase(key);
+          has_removed = true;
+          break;
+        }
       }
-    }
-  } while (has_removed);
+    } while (has_removed);
+  };
+
+  delete_connectivity_stencil(m_stored_cell_to_cell_stencil_map);
+  delete_connectivity_stencil(m_stored_node_to_cell_stencil_map);
 }
diff --git a/src/mesh/StencilManager.hpp b/src/mesh/StencilManager.hpp
index 814d8e476d71ada4ec66322e1f3ed80035f808e2..8f07e282848af23bca8604e5bbc68754c63251a7 100644
--- a/src/mesh/StencilManager.hpp
+++ b/src/mesh/StencilManager.hpp
@@ -3,7 +3,8 @@
 
 #include <mesh/IBoundaryDescriptor.hpp>
 #include <mesh/IConnectivity.hpp>
-#include <mesh/Stencil.hpp>
+#include <mesh/StencilArray.hpp>
+#include <mesh/StencilDescriptor.hpp>
 
 #include <memory>
 #include <set>
@@ -24,13 +25,13 @@ class StencilManager
   struct Key
   {
     size_t connectivity_id;
-    size_t degree;
+    StencilDescriptor stencil_descriptor;
     BoundaryDescriptorList symmetry_boundary_descriptor_list;
 
     PUGS_INLINE bool
     operator==(const Key& k) const
     {
-      if ((connectivity_id != k.connectivity_id) or (degree != k.degree) or
+      if ((connectivity_id != k.connectivity_id) or (stencil_descriptor != k.stencil_descriptor) or
           (symmetry_boundary_descriptor_list.size() != k.symmetry_boundary_descriptor_list.size())) {
         return false;
       }
@@ -50,18 +51,33 @@ class StencilManager
       return boundary_descriptor_set == k_boundary_descriptor_set;
     }
   };
+
   struct HashKey
   {
     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);
+      return ((std::hash<decltype(Key::connectivity_id)>()(key.connectivity_id)) ^
+              (std::hash<std::decay_t<decltype(Key::stencil_descriptor.numberOfLayers())>>()(
+                 key.stencil_descriptor.numberOfLayers())
+               << 1));
     }
   };
 
-  std::unordered_map<Key, std::shared_ptr<const Stencil>, HashKey> m_stored_stencil_map;
+  template <ItemType source_item_type, ItemType target_item_type>
+  using StoredStencilTMap =
+    std::unordered_map<Key, std::shared_ptr<const StencilArray<source_item_type, target_item_type>>, HashKey>;
+
+  StoredStencilTMap<ItemType::cell, ItemType::cell> m_stored_cell_to_cell_stencil_map;
+  StoredStencilTMap<ItemType::node, ItemType::cell> m_stored_node_to_cell_stencil_map;
+
+  template <ItemType source_item_type, ItemType target_item_type>
+  const StencilArray<source_item_type, target_item_type>& _getStencilArray(
+    const IConnectivity& connectivity,
+    const StencilDescriptor& stencil_descriptor,
+    const BoundaryDescriptorList& symmetry_boundary_descriptor_list,
+    StoredStencilTMap<source_item_type, target_item_type>& stored_source_to_target_stencil_map);
 
  public:
   static void create();
@@ -77,9 +93,15 @@ 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,
+    const StencilDescriptor& stencil_descriptor,
+    const BoundaryDescriptorList& symmetry_boundary_descriptor_list = {});
+
+  const NodeToCellStencilArray& getNodeToCellStencilArray(
+    const IConnectivity& i_connectivity,
+    const StencilDescriptor& stencil_descriptor,
+    const BoundaryDescriptorList& symmetry_boundary_descriptor_list = {});
 
   StencilManager(const StencilManager&) = delete;
   StencilManager(StencilManager&&)      = delete;
diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt
index 051e96480431aaa733059e2cf9c9518b29216812..fc975ba9b8234d8f70bd850a39bef58360e64108 100644
--- a/src/output/CMakeLists.txt
+++ b/src/output/CMakeLists.txt
@@ -4,6 +4,8 @@ add_library(
   PugsOutput
   GnuplotWriter.cpp
   GnuplotWriter1D.cpp
+  GnuplotWriterBase.cpp
+  GnuplotWriterRaw.cpp
   VTKWriter.cpp
   WriterBase.cpp)
 
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index 991ec74537fffaaa2bef855f5c4af88f7acfa03f..83e2d8effaae57605d0bcf6f09f0a521c6524ff7 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -7,94 +7,14 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
+#include <utils/Demangle.hpp>
 #include <utils/Filesystem.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PugsTraits.hpp>
 #include <utils/RevisionInfo.hpp>
 #include <utils/Stringify.hpp>
 
-#include <utils/Demangle.hpp>
-
 #include <fstream>
-#include <iomanip>
-
-std::string
-GnuplotWriter::_getDateAndVersionComment() const
-{
-  std::ostringstream os;
-
-  std::time_t now = std::time(nullptr);
-  os << "#  Generated by pugs: " << std::ctime(&now);
-  os << "#  version: " << RevisionInfo::version() << '\n';
-  os << "#  tag:  " << RevisionInfo::gitTag() << '\n';
-  os << "#  HEAD: " << RevisionInfo::gitHead() << '\n';
-  os << "#  hash: " << RevisionInfo::gitHash() << " (" << ((RevisionInfo::gitIsClean()) ? "clean" : "dirty") << ")\n";
-  os << '\n';
-
-  return os.str();
-}
-
-std::string
-GnuplotWriter::_getFilename() const
-{
-  std::ostringstream sout;
-  sout << m_base_filename;
-  if (m_period_manager.has_value()) {
-    sout << '.' << std::setfill('0') << std::setw(4) << m_period_manager->nbSavedTimes();
-  }
-  sout << ".gnu";
-  return sout.str();
-}
-
-template <size_t Dimension>
-void
-GnuplotWriter::_writePreamble(const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const
-{
-  fout << "# list of data\n";
-  fout << "# 1:x";
-  if constexpr (Dimension > 1) {
-    fout << " 2:y";
-  }
-  uint64_t i = Dimension + 1;
-  for (const auto& i_named_item_data : output_named_item_data_set) {
-    const std::string name        = i_named_item_data.first;
-    const auto& item_data_variant = i_named_item_data.second;
-    std::visit(
-      [&](auto&& item_data) {
-        using ItemDataType = std::decay_t<decltype(item_data)>;
-        using DataType     = std::decay_t<typename ItemDataType::data_type>;
-        if constexpr (is_item_value_v<ItemDataType>) {
-          if constexpr (std::is_arithmetic_v<DataType>) {
-            fout << ' ' << i++ << ':' << name;
-          } else if constexpr (is_tiny_vector_v<DataType>) {
-            for (size_t j = 0; j < DataType{}.dimension(); ++j) {
-              fout << ' ' << i++ << ':' << name << '[' << j << ']';
-            }
-          } else if constexpr (is_tiny_matrix_v<DataType>) {
-            for (size_t j = 0; j < DataType{}.numberOfRows(); ++j) {
-              for (size_t k = 0; k < DataType{}.numberOfColumns(); ++k) {
-                fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
-              }
-            }
-          } else {
-            throw UnexpectedError("invalid data type");
-          }
-        } else if constexpr (is_item_array_v<ItemDataType>) {
-          if constexpr (std::is_arithmetic_v<DataType>) {
-            for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
-              fout << ' ' << i++ << ':' << name << '[' << j << ']';
-            }
-          } else {
-            throw UnexpectedError("invalid data type");
-          }
-        } else {
-          throw UnexpectedError("invalid ItemData type");
-        }
-      },
-      item_data_variant);
-  }
-  fout << "\n\n";
-}
 
 template <typename DataType>
 void
@@ -263,7 +183,7 @@ GnuplotWriter::_write(const MeshType& mesh,
     if (time.has_value()) {
       fout << "# time = " << *time << "\n\n";
     }
-    this->_writePreamble<MeshType::Dimension>(output_named_item_data_set, fout);
+    this->_writePreamble(MeshType::Dimension, output_named_item_data_set, true /*store coordinates*/, fout);
   }
 
   for (const auto& [name, item_data_variant] : output_named_item_data_set) {
diff --git a/src/output/GnuplotWriter.hpp b/src/output/GnuplotWriter.hpp
index 3069b79da716d2640e8446822e23fcab37b98316..0c14204d930bc947832a484c617a9aa5965b04e6 100644
--- a/src/output/GnuplotWriter.hpp
+++ b/src/output/GnuplotWriter.hpp
@@ -1,7 +1,7 @@
 #ifndef GNUPLOT_WRITER_HPP
 #define GNUPLOT_WRITER_HPP
 
-#include <output/WriterBase.hpp>
+#include <output/GnuplotWriterBase.hpp>
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
@@ -13,13 +13,9 @@ class MeshVariant;
 #include <optional>
 #include <string>
 
-class GnuplotWriter final : public WriterBase
+class GnuplotWriter final : public GnuplotWriterBase
 {
  private:
-  std::string _getDateAndVersionComment() const;
-
-  std::string _getFilename() const;
-
   template <typename DataType>
   void _writeCellData(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const;
 
@@ -32,9 +28,6 @@ class GnuplotWriter final : public WriterBase
   template <typename DataType>
   void _writeNodeData(const NodeArray<DataType>& node_value, NodeId cell_id, std::ostream& fout) const;
 
-  template <size_t Dimension>
-  void _writePreamble(const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const;
-
   template <typename ItemDataT>
   void _writeData(const ItemDataT& item_data,
                   [[maybe_unused]] CellId cell_id,
@@ -61,9 +54,11 @@ class GnuplotWriter final : public WriterBase
   void _writeMesh(const MeshVariant& mesh_v) const final;
 
  public:
-  GnuplotWriter(const std::string& base_filename) : WriterBase(base_filename) {}
+  GnuplotWriter(const std::string& base_filename) : GnuplotWriterBase(base_filename) {}
 
-  GnuplotWriter(const std::string& base_filename, const double time_period) : WriterBase(base_filename, time_period) {}
+  GnuplotWriter(const std::string& base_filename, const double time_period)
+    : GnuplotWriterBase(base_filename, time_period)
+  {}
 
   ~GnuplotWriter() = default;
 };
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index 84e8aa080c6802eb546a4d0183fb446201e31ac8..799a8d51d7dceafeea7f3555cbd7d73088d116c0 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -7,118 +7,13 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
+#include <utils/Demangle.hpp>
 #include <utils/Filesystem.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PugsTraits.hpp>
-#include <utils/RevisionInfo.hpp>
 #include <utils/Stringify.hpp>
 
-#include <utils/Demangle.hpp>
-
 #include <fstream>
-#include <iomanip>
-
-std::string
-GnuplotWriter1D::_getDateAndVersionComment() const
-{
-  std::ostringstream os;
-
-  std::time_t now = std::time(nullptr);
-  os << "#  Generated by pugs: " << std::ctime(&now);
-  os << "#  version: " << RevisionInfo::version() << '\n';
-  os << "#  tag:  " << RevisionInfo::gitTag() << '\n';
-  os << "#  HEAD: " << RevisionInfo::gitHead() << '\n';
-  os << "#  hash: " << RevisionInfo::gitHash() << " (" << ((RevisionInfo::gitIsClean()) ? "clean" : "dirty") << ")\n";
-  os << '\n';
-
-  return os.str();
-}
-
-std::string
-GnuplotWriter1D::_getFilename() const
-{
-  std::ostringstream sout;
-  sout << m_base_filename;
-  if (m_period_manager.has_value()) {
-    sout << '.' << std::setfill('0') << std::setw(4) << m_period_manager->nbSavedTimes();
-  }
-  sout << ".gnu";
-  return sout.str();
-}
-
-template <typename ItemDataType>
-bool
-GnuplotWriter1D::_is_cell_data(const ItemDataType&) const
-{
-  return ItemDataType::item_t == ItemType::cell;
-}
-
-template <typename ItemDataType>
-bool
-GnuplotWriter1D::_is_face_data(const ItemDataType&) const
-{
-  return ItemDataType::item_t == ItemType::face;
-}
-
-template <typename ItemDataType>
-bool
-GnuplotWriter1D::_is_edge_data(const ItemDataType&) const
-{
-  return ItemDataType::item_t == ItemType::edge;
-}
-
-template <typename ItemDataType>
-bool
-GnuplotWriter1D::_is_node_data(const ItemDataType&) const
-{
-  return ItemDataType::item_t == ItemType::node;
-}
-
-void
-GnuplotWriter1D::_writePreamble(const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const
-{
-  fout << "# list of data\n";
-  fout << "# 1:x";
-  uint64_t i = 2;
-  for (const auto& i_named_item_data : output_named_item_data_set) {
-    const std::string name        = i_named_item_data.first;
-    const auto& item_data_variant = i_named_item_data.second;
-    std::visit(
-      [&](auto&& item_data) {
-        using ItemDataType = std::decay_t<decltype(item_data)>;
-        using DataType     = std::decay_t<typename ItemDataType::data_type>;
-        if constexpr (is_item_value_v<ItemDataType>) {
-          if constexpr (std::is_arithmetic_v<DataType>) {
-            fout << ' ' << i++ << ':' << name;
-          } else if constexpr (is_tiny_vector_v<DataType>) {
-            for (size_t j = 0; j < DataType{}.dimension(); ++j) {
-              fout << ' ' << i++ << ':' << name << '[' << j << ']';
-            }
-          } else if constexpr (is_tiny_matrix_v<DataType>) {
-            for (size_t j = 0; j < DataType{}.numberOfRows(); ++j) {
-              for (size_t k = 0; k < DataType{}.numberOfColumns(); ++k) {
-                fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
-              }
-            }
-          } else {
-            throw UnexpectedError("invalid data type");
-          }
-        } else if constexpr (is_item_array_v<ItemDataType>) {
-          if constexpr (std::is_arithmetic_v<DataType>) {
-            for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
-              fout << ' ' << i++ << ':' << name << '[' << j << ']';
-            }
-          } else {
-            throw UnexpectedError("invalid data type");
-          }
-        } else {
-          throw UnexpectedError("invalid ItemData type");
-        }
-      },
-      item_data_variant);
-  }
-  fout << "\n\n";
-}
 
 template <typename DataType, ItemType item_type>
 size_t
@@ -339,7 +234,7 @@ GnuplotWriter1D::_write(const MeshType& mesh,
       fout << "# time = " << *time << "\n\n";
     }
 
-    _writePreamble(output_named_item_data_set, fout);
+    _writePreamble(MeshType::Dimension, output_named_item_data_set, true /*store coordinates*/, fout);
   }
 
   if (has_cell_data) {
@@ -355,7 +250,7 @@ GnuplotWriter1D::_writeMesh(const MeshVariant&) const
   std::ostringstream errorMsg;
   errorMsg << "gnuplot_1d_writer does not write meshes\n"
            << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue
-           << "gnuplot_writer" << rang::fg::reset << "  instead";
+           << "gnuplot_writer" << rang::fg::reset << " instead";
   throw NormalError(errorMsg.str());
 }
 
diff --git a/src/output/GnuplotWriter1D.hpp b/src/output/GnuplotWriter1D.hpp
index 66661918252856d8b6a069773c931c705a9b9330..32e4501873b175fb9af10123935d4666cd331404 100644
--- a/src/output/GnuplotWriter1D.hpp
+++ b/src/output/GnuplotWriter1D.hpp
@@ -1,7 +1,7 @@
 #ifndef GNUPLOT_WRITER_1D_HPP
 #define GNUPLOT_WRITER_1D_HPP
 
-#include <output/WriterBase.hpp>
+#include <output/GnuplotWriterBase.hpp>
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
@@ -13,25 +13,9 @@ class MeshVariant;
 #include <optional>
 #include <string>
 
-class GnuplotWriter1D final : public WriterBase
+class GnuplotWriter1D final : public GnuplotWriterBase
 {
  private:
-  std::string _getDateAndVersionComment() const;
-
-  std::string _getFilename() const;
-
-  template <typename ItemDataType>
-  bool _is_cell_data(const ItemDataType&) const;
-
-  template <typename ItemDataType>
-  bool _is_face_data(const ItemDataType&) const;
-
-  template <typename ItemDataType>
-  bool _is_edge_data(const ItemDataType&) const;
-
-  template <typename ItemDataType>
-  bool _is_node_data(const ItemDataType&) const;
-
   template <typename DataType, ItemType item_type>
   size_t _itemDataNbRow(const ItemValue<DataType, item_type>&) const;
 
@@ -43,8 +27,6 @@ class GnuplotWriter1D final : public WriterBase
                        const OutputNamedItemDataSet& output_named_item_data_set,
                        std::ostream& fout) const;
 
-  void _writePreamble(const OutputNamedItemDataSet& output_named_item_value_set, std::ostream& fout) const;
-
   template <MeshConcept MeshType>
   void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_value_set,
@@ -60,9 +42,10 @@ class GnuplotWriter1D final : public WriterBase
   void _writeMesh(const MeshVariant& mesh_v) const final;
 
  public:
-  GnuplotWriter1D(const std::string& base_filename) : WriterBase(base_filename) {}
+  GnuplotWriter1D(const std::string& base_filename) : GnuplotWriterBase(base_filename) {}
 
-  GnuplotWriter1D(const std::string& base_filename, const double time_period) : WriterBase(base_filename, time_period)
+  GnuplotWriter1D(const std::string& base_filename, const double time_period)
+    : GnuplotWriterBase(base_filename, time_period)
   {}
 
   ~GnuplotWriter1D() = default;
diff --git a/src/output/GnuplotWriterBase.cpp b/src/output/GnuplotWriterBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..110c03f77406c3e68e0ffb35b8528c5639948a94
--- /dev/null
+++ b/src/output/GnuplotWriterBase.cpp
@@ -0,0 +1,94 @@
+#include <output/GnuplotWriterBase.hpp>
+
+#include <output/OutputNamedItemValueSet.hpp>
+#include <utils/RevisionInfo.hpp>
+
+#include <ctime>
+#include <iomanip>
+#include <sstream>
+
+std::string
+GnuplotWriterBase::_getDateAndVersionComment() const
+{
+  std::ostringstream os;
+
+  std::time_t now = std::time(nullptr);
+  os << "#  Generated by pugs: " << std::ctime(&now);
+  os << "#  version: " << RevisionInfo::version() << '\n';
+  os << "#  tag:  " << RevisionInfo::gitTag() << '\n';
+  os << "#  HEAD: " << RevisionInfo::gitHead() << '\n';
+  os << "#  hash: " << RevisionInfo::gitHash() << " (" << ((RevisionInfo::gitIsClean()) ? "clean" : "dirty") << ")\n";
+  os << '\n';
+
+  return os.str();
+}
+
+std::string
+GnuplotWriterBase::_getFilename() const
+{
+  std::ostringstream sout;
+  sout << m_base_filename;
+  if (m_period_manager.has_value()) {
+    sout << '.' << std::setfill('0') << std::setw(4) << m_period_manager->nbSavedTimes();
+  }
+  sout << ".gnu";
+  return sout.str();
+}
+
+void
+GnuplotWriterBase::_writePreamble(const size_t& dimension,
+                                  const OutputNamedItemDataSet& output_named_item_data_set,
+                                  const bool& store_coordinates,
+                                  std::ostream& fout) const
+{
+  fout << "# list of data\n";
+  uint64_t i = 0;
+
+  fout << "#";
+  if (store_coordinates) {
+    fout << " 1:x";
+    if (dimension > 1) {
+      fout << " 2:y";
+    }
+    i = dimension;
+  }
+
+  for (const auto& i_named_item_data : output_named_item_data_set) {
+    const std::string name        = i_named_item_data.first;
+    const auto& item_data_variant = i_named_item_data.second;
+    std::visit(
+      [&](auto&& item_data) {
+        using ItemDataType = std::decay_t<decltype(item_data)>;
+        using DataType     = std::decay_t<typename ItemDataType::data_type>;
+        if constexpr (is_item_value_v<ItemDataType>) {
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            fout << ' ' << i++ << ':' << name;
+          } else if constexpr (is_tiny_vector_v<DataType>) {
+            for (size_t j = 0; j < DataType{}.dimension(); ++j) {
+              fout << ' ' << i++ << ':' << name << '[' << j << ']';
+            }
+          } else if constexpr (is_tiny_matrix_v<DataType>) {
+            for (size_t j = 0; j < DataType{}.numberOfRows(); ++j) {
+              for (size_t k = 0; k < DataType{}.numberOfColumns(); ++k) {
+                fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
+              }
+            }
+          } else {
+            throw UnexpectedError("invalid data type");
+          }
+        } else if constexpr (is_item_array_v<ItemDataType>) {
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
+              fout << ' ' << i++ << ':' << name << '[' << j << ']';
+            }
+          } else {
+            throw UnexpectedError("invalid data type");
+          }
+        } else {
+          throw UnexpectedError("invalid ItemData type");
+        }
+      },
+      item_data_variant);
+  }
+  fout << "\n\n";
+}
diff --git a/src/output/GnuplotWriterBase.hpp b/src/output/GnuplotWriterBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a99390facb8d2bd67cc323d676d305f26401dd6e
--- /dev/null
+++ b/src/output/GnuplotWriterBase.hpp
@@ -0,0 +1,58 @@
+#ifndef GNUPLOT_WRITER_BASE_HPP
+#define GNUPLOT_WRITER_BASE_HPP
+
+#include <output/WriterBase.hpp>
+
+#include <utils/PugsMacros.hpp>
+
+class GnuplotWriterBase : public WriterBase
+{
+ protected:
+  std::string _getDateAndVersionComment() const;
+
+  std::string _getFilename() const;
+
+  void _writePreamble(const size_t& dimension,
+                      const OutputNamedItemDataSet& output_named_item_data_set,
+                      const bool& store_coordinates,
+                      std::ostream& fout) const;
+
+  template <typename ItemDataType>
+  PUGS_INLINE bool
+  _is_cell_data(const ItemDataType&) const
+  {
+    return ItemDataType::item_t == ItemType::cell;
+  }
+
+  template <typename ItemDataType>
+  PUGS_INLINE bool
+  _is_face_data(const ItemDataType&) const
+  {
+    return ItemDataType::item_t == ItemType::face;
+  }
+
+  template <typename ItemDataType>
+  PUGS_INLINE bool
+  _is_edge_data(const ItemDataType&) const
+  {
+    return ItemDataType::item_t == ItemType::edge;
+  }
+
+  template <typename ItemDataType>
+  PUGS_INLINE bool
+  _is_node_data(const ItemDataType&) const
+  {
+    return ItemDataType::item_t == ItemType::node;
+  }
+
+ public:
+  GnuplotWriterBase(const std::string& base_filename, const double& time_period)
+    : WriterBase(base_filename, time_period)
+  {}
+
+  GnuplotWriterBase(const std::string& base_filename) : WriterBase(base_filename) {}
+
+  virtual ~GnuplotWriterBase() = default;
+};
+
+#endif   // GNUPLOT_WRITER_BASE_HPP
diff --git a/src/output/GnuplotWriterRaw.cpp b/src/output/GnuplotWriterRaw.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ac93ba0029f3d73583c98d62c85beb7901f3569f
--- /dev/null
+++ b/src/output/GnuplotWriterRaw.cpp
@@ -0,0 +1,269 @@
+#include <output/GnuplotWriterRaw.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <utils/Filesystem.hpp>
+#include <utils/Messenger.hpp>
+#include <utils/PugsTraits.hpp>
+#include <utils/RevisionInfo.hpp>
+#include <utils/Stringify.hpp>
+
+#include <utils/Demangle.hpp>
+
+#include <fstream>
+#include <iomanip>
+
+template <typename DataType, ItemType item_type>
+size_t
+GnuplotWriterRaw::_itemDataNbRow(const ItemValue<DataType, item_type>&) const
+{
+  if constexpr (std::is_arithmetic_v<DataType>) {
+    return 1;
+  } else if constexpr (is_tiny_vector_v<std::decay_t<DataType>>) {
+    return DataType::Dimension;
+  } else if constexpr (is_tiny_matrix_v<std::decay_t<DataType>>) {
+    return DataType{}.dimension();
+  } else {
+    throw UnexpectedError("invalid data type for cell value output: " + demangle<DataType>());
+  }
+}
+
+template <typename DataType, ItemType item_type>
+size_t
+GnuplotWriterRaw::_itemDataNbRow(const ItemArray<DataType, item_type>& item_array) const
+{
+  return item_array.sizeOfArrays();
+}
+
+template <MeshConcept MeshType, ItemType item_type>
+void
+GnuplotWriterRaw::_writeItemDatas(const MeshType& mesh,
+                                  const OutputNamedItemDataSet& output_named_item_data_set,
+                                  std::ostream& fout) const
+{
+  using ItemId = ItemIdT<item_type>;
+
+  const size_t number_of_columns = [&] {
+    size_t nb_columns = 0;
+    for (const auto& [name, item_data] : output_named_item_data_set) {
+      std::visit([&](auto&& value) { nb_columns += _itemDataNbRow(value); }, item_data);
+    }
+    return nb_columns;
+  }();
+
+  auto is_owned = mesh.connectivity().template isOwned<item_type>();
+
+  const size_t& number_of_owned_lines = [&]() {
+    if (parallel::size() > 1) {
+      size_t number_of_owned_items = 0;
+      for (ItemId item_id = 0; item_id < mesh.template numberOf<item_type>(); ++item_id) {
+        if (is_owned[item_id]) {
+          ++number_of_owned_items;
+        }
+      }
+
+      return number_of_owned_items;
+    } else {
+      return mesh.template numberOf<item_type>();
+    }
+  }();
+
+  Array<double> values{number_of_columns * number_of_owned_lines};
+
+  size_t column_number = 0;
+  for (const auto& [name, output_item_data] : output_named_item_data_set) {
+    std::visit(
+      [&](auto&& item_data) {
+        using ItemDataT = std::decay_t<decltype(item_data)>;
+        if constexpr (ItemDataT::item_t == item_type) {
+          if constexpr (is_item_value_v<ItemDataT>) {
+            using DataT  = std::decay_t<typename ItemDataT::data_type>;
+            size_t index = 0;
+            for (ItemId item_id = 0; item_id < item_data.numberOfItems(); ++item_id) {
+              if (is_owned[item_id]) {
+                if constexpr (std::is_arithmetic_v<DataT>) {
+                  values[number_of_columns * index + column_number] = item_data[item_id];
+                } else if constexpr (is_tiny_vector_v<DataT>) {
+                  const size_t k = number_of_columns * index + column_number;
+                  for (size_t j = 0; j < DataT::Dimension; ++j) {
+                    values[k + j] = item_data[item_id][j];
+                  }
+                } else if constexpr (is_tiny_matrix_v<DataT>) {
+                  size_t k = number_of_columns * index + column_number;
+                  for (size_t i = 0; i < DataT{}.numberOfRows(); ++i) {
+                    for (size_t j = 0; j < DataT{}.numberOfColumns(); ++j) {
+                      values[k++] = item_data[item_id](i, j);
+                    }
+                  }
+                }
+                ++index;
+              }
+            }
+          } else {
+            using DataT  = std::decay_t<typename ItemDataT::data_type>;
+            size_t index = 0;
+            for (ItemId item_id = 0; item_id < item_data.numberOfItems(); ++item_id) {
+              if (is_owned[item_id]) {
+                if constexpr (std::is_arithmetic_v<DataT>) {
+                  const size_t k = number_of_columns * index + column_number;
+                  for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
+                    values[k + j] = item_data[item_id][j];
+                  }
+                }
+                ++index;
+              }
+            }
+          }
+        }
+        column_number += _itemDataNbRow(item_data);
+      },
+      output_item_data);
+  }
+
+  if (parallel::size() > 1) {
+    values = parallel::gatherVariable(values, 0);
+  }
+
+  if (parallel::rank() == 0) {
+    Assert(values.size() % number_of_columns == 0);
+
+    std::vector<size_t> line_numbers(values.size() / number_of_columns);
+    for (size_t i = 0; i < line_numbers.size(); ++i) {
+      line_numbers[i] = i;
+    }
+
+    std::sort(line_numbers.begin(), line_numbers.end(),
+              [&](size_t i, size_t j) { return values[i * number_of_columns] < values[j * number_of_columns]; });
+
+    for (auto i_line : line_numbers) {
+      fout << values[i_line * number_of_columns];
+      for (size_t j = 1; j < number_of_columns; ++j) {
+        fout << ' ' << values[i_line * number_of_columns + j];
+      }
+      fout << '\n';
+    }
+  }
+}
+
+template <MeshConcept MeshType>
+void
+GnuplotWriterRaw::_write(const MeshType& mesh,
+                         const OutputNamedItemDataSet& output_named_item_data_set,
+                         std::optional<double> time) const
+{
+  bool has_cell_data = false;
+  for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+    has_cell_data |= std::visit([&](auto&& item_data) { return this->_is_cell_data(item_data); }, item_data_variant);
+  }
+
+  for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+    std::visit(
+      [&, name = name](auto&& item_data) {
+        if (this->_is_face_data(item_data)) {
+          std::ostringstream error_msg;
+          error_msg << "gnuplot_raw_writer does not support face data, cannot save variable \"" << rang::fgB::yellow
+                    << name << rang::fg::reset << '"';
+          throw NormalError(error_msg.str());
+        }
+      },
+      item_data_variant);
+  }
+
+  for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+    std::visit(
+      [&, name = name](auto&& item_data) {
+        if (this->_is_edge_data(item_data)) {
+          std::ostringstream error_msg;
+          error_msg << "gnuplot_1d_writer does not support edge data, cannot save variable \"" << rang::fgB::yellow
+                    << name << rang::fg::reset << '"';
+          throw NormalError(error_msg.str());
+        }
+      },
+      item_data_variant);
+  }
+
+  bool has_node_data = false;
+  for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+    has_node_data |= std::visit([&](auto&& item_data) { return this->_is_node_data(item_data); }, item_data_variant);
+  }
+
+  if (has_cell_data and has_node_data) {
+    throw NormalError("cannot store both node and cell data in the same gnuplot file");
+  }
+
+  std::ofstream fout;
+
+  if (parallel::rank() == 0) {
+    fout.open(_getFilename());
+    if (not fout) {
+      std::ostringstream error_msg;
+      error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilename() << rang::fg::reset << '"';
+      throw NormalError(error_msg.str());
+    }
+
+    fout.precision(15);
+    fout.setf(std::ios_base::scientific);
+    fout << _getDateAndVersionComment();
+
+    if (time.has_value()) {
+      fout << "# time = " << *time << "\n\n";
+    }
+
+    _writePreamble(MeshType::Dimension, output_named_item_data_set, false /*do not store coordinates*/, fout);
+  }
+
+  if (has_cell_data) {
+    this->_writeItemDatas<MeshType, ItemType::cell>(mesh, output_named_item_data_set, fout);
+  } else {   // has_node_value
+    this->_writeItemDatas<MeshType, ItemType::node>(mesh, output_named_item_data_set, fout);
+  }
+}
+
+void
+GnuplotWriterRaw::_writeMesh(const MeshVariant&) const
+{
+  std::ostringstream errorMsg;
+  errorMsg << "gnuplot_raw_writer does not write meshes\n"
+           << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue
+           << "gnuplot_writer" << rang::fg::reset << " instead";
+  throw NormalError(errorMsg.str());
+}
+
+void
+GnuplotWriterRaw::_writeAtTime(const MeshVariant& mesh_v,
+                               const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+                               double time) const
+{
+  OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
+
+  std::visit(
+    [&](auto&& p_mesh) {
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
+      if constexpr (MeshType::Dimension == 1) {
+        this->_write(*p_mesh, output_named_item_data_set, time);
+      } else if constexpr (MeshType::Dimension == 2) {
+        std::ostringstream errorMsg;
+        errorMsg << "gnuplot_raw_writer is not available in dimension " << stringify(MeshType::Dimension) << '\n'
+                 << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue
+                 << "gnuplot_writer" << rang::fg::reset << " in dimension 2";
+        throw NormalError(errorMsg.str());
+      } else {
+        throw NormalError("gnuplot format is not available in dimension " + stringify(MeshType::Dimension));
+      }
+    },
+    mesh_v.variant());
+}
+
+void
+GnuplotWriterRaw::_write(const MeshVariant& mesh_v,
+                         const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
+{
+  OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
+
+  std::visit([&](auto&& p_mesh) { this->_write(*p_mesh, output_named_item_data_set, {}); }, mesh_v.variant());
+}
diff --git a/src/output/GnuplotWriterRaw.hpp b/src/output/GnuplotWriterRaw.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..481bb2a4bfc548c89ea684ec7a841df67bede57b
--- /dev/null
+++ b/src/output/GnuplotWriterRaw.hpp
@@ -0,0 +1,54 @@
+#ifndef GNUPLOT_WRITER_RAW_HPP
+#define GNUPLOT_WRITER_RAW_HPP
+
+#include <output/GnuplotWriterBase.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <output/OutputNamedItemValueSet.hpp>
+
+class MeshVariant;
+
+#include <optional>
+#include <string>
+
+class GnuplotWriterRaw final : public GnuplotWriterBase
+{
+ private:
+  template <typename DataType, ItemType item_type>
+  size_t _itemDataNbRow(const ItemValue<DataType, item_type>&) const;
+
+  template <typename DataType, ItemType item_type>
+  size_t _itemDataNbRow(const ItemArray<DataType, item_type>&) const;
+
+  template <MeshConcept MeshType, ItemType item_type>
+  void _writeItemDatas(const MeshType& mesh,
+                       const OutputNamedItemDataSet& output_named_item_data_set,
+                       std::ostream& fout) const;
+
+  template <MeshConcept MeshType>
+  void _write(const MeshType& mesh,
+              const OutputNamedItemDataSet& output_named_item_value_set,
+              std::optional<double> time) const;
+
+  void _writeAtTime(const MeshVariant& mesh_v,
+                    const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+                    double time) const final;
+
+  void _write(const MeshVariant& mesh_v,
+              const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
+
+  void _writeMesh(const MeshVariant& mesh_v) const final;
+
+ public:
+  GnuplotWriterRaw(const std::string& base_filename) : GnuplotWriterBase(base_filename) {}
+
+  GnuplotWriterRaw(const std::string& base_filename, const double time_period)
+    : GnuplotWriterBase(base_filename, time_period)
+  {}
+
+  ~GnuplotWriterRaw() = default;
+};
+
+#endif   // GNUPLOT_WRITER_RAW_HPP
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index e4db63bd7ec57502c08a1e3f212386ebe8215859..5c8e894c0e2d804813a727f18fd5de44ecdc9cd8 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -17,6 +17,7 @@
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/GlobalVariableManager.hpp>
 #include <utils/Socket.hpp>
 
 #include <variant>
@@ -371,6 +372,10 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
                  const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
+    if ((parallel::size() > 1) and (GlobalVariableManager::instance().getNumberOfGhostLayers() == 0)) {
+      throw NormalError("Acoustic solver requires 1 layer of ghost cells in parallel");
+    }
+
     std::shared_ptr mesh_v = getCommonMesh({rho_v, c_v, u_v, p_v});
     if (not mesh_v) {
       throw NormalError("discrete functions are not defined on the same mesh");
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 14739b777f524ce12dffdfcfeed01b06cdaa4ad7..13f1117862a81296211c5cead802b49268132a70 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -21,6 +21,7 @@ add_library(
   DiscreteFunctionVectorIntegrator.cpp
   DiscreteFunctionVectorInterpoler.cpp
   FluxingAdvectionSolver.cpp
+  HyperelasticSolver.cpp
   PolynomialReconstruction.cpp
   RusanovEulerianCompositeSolverTools.cpp
   RusanovEulerianCompositeSolver.cpp
@@ -29,8 +30,6 @@ add_library(
   RusanovEulerianCompositeSolver_o2.cpp
   RusanovEulerianCompositeSolver_v2_o2.cpp
   RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
-
-  test_reconstruction.cpp
 )
 
 target_link_libraries(
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index a2ebb2b256a3957169cf52457fd87a3d3236e8f3..d57cb9add43530ff8bbbaee3ff9e4786a61bb57c 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -23,6 +23,7 @@
 #include <scheme/InflowBoundaryConditionDescriptor.hpp>
 #include <scheme/OutflowBoundaryConditionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/GlobalVariableManager.hpp>
 
 #include <variant>
 #include <vector>
@@ -201,6 +202,10 @@ class FluxingAdvectionSolver
       throw NormalError("old and new meshes must share the same connectivity");
     }
 
+    if ((parallel::size() > 1) and (GlobalVariableManager::instance().getNumberOfGhostLayers() == 0)) {
+      throw NormalError("Fluxing advection requires 1 layer of ghost cells in parallel");
+    }
+
     this->_computeGeometricalData();
   }
 
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index b0cacf9d9768a73e5b1127d91569c0832254a523..2af32e157475b0fe99b9f2cc9cc6c445e05681b3 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -19,6 +19,7 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshFlatFaceBoundary.hpp>
 #include <mesh/NamedBoundaryDescriptor.hpp>
+#include <mesh/StencilDescriptor.hpp>
 #include <mesh/StencilManager.hpp>
 #include <scheme/DiscreteFunctionDPkVariant.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
@@ -769,8 +770,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.stencilDescriptor(),
+                                                         m_descriptor.symmetryBoundaryDescriptorList());
 
   auto xr = mesh.xr();
   auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
@@ -782,10 +784,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].second;
+      auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i].stencilArray();
       stencil_size += ghost_stencil[cell_id].size();
     }
 
@@ -814,6 +816,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 +866,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 +900,9 @@ PolynomialReconstruction::_build(
                   }
                 }
 
-                for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryBoundaryStencilList().size(); ++i_symmetry) {
-                  auto& ghost_stencil  = stencil.symmetryBoundaryStencilList()[i_symmetry].second;
+                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 +957,9 @@ PolynomialReconstruction::_build(
                   }
                 }
 
-                for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryBoundaryStencilList().size(); ++i_symmetry) {
-                  auto& ghost_stencil  = stencil.symmetryBoundaryStencilList()[i_symmetry].second;
+                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 +993,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].second;
+          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 +1039,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].second;
+              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 +1085,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].second;
+            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 +1127,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].second;
+          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/src/scheme/PolynomialReconstructionDescriptor.hpp b/src/scheme/PolynomialReconstructionDescriptor.hpp
index 8a51a85714a0f58749e56cc00e8c31664b2d36f1..44fdd2082e0efdf550bc1490a6df18d5658ff9bd 100644
--- a/src/scheme/PolynomialReconstructionDescriptor.hpp
+++ b/src/scheme/PolynomialReconstructionDescriptor.hpp
@@ -2,6 +2,7 @@
 #define POLYNOMIAL_RECONSTRUCTION_DESCRIPTOR_HPP
 
 #include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/StencilDescriptor.hpp>
 #include <scheme/IntegrationMethodType.hpp>
 #include <utils/PugsMacros.hpp>
 
@@ -17,6 +18,7 @@ class PolynomialReconstructionDescriptor
  private:
   IntegrationMethodType m_integration_method;
   size_t m_degree;
+  StencilDescriptor m_stencil_descriptor;
 
   BoundaryDescriptorList m_symmetry_boundary_descriptor_list;
 
@@ -37,6 +39,13 @@ class PolynomialReconstructionDescriptor
     return m_degree;
   }
 
+  PUGS_INLINE
+  const StencilDescriptor&
+  stencilDescriptor() const
+  {
+    return m_stencil_descriptor;
+  }
+
   PUGS_INLINE
   const BoundaryDescriptorList&
   symmetryBoundaryDescriptorList() const
@@ -58,13 +67,6 @@ class PolynomialReconstructionDescriptor
     return m_row_weighting;
   }
 
-  PUGS_INLINE
-  void
-  setDegree(const size_t degree)
-  {
-    m_degree = degree;
-  }
-
   PUGS_INLINE
   void
   setPreconditioning(const bool preconditioning)
@@ -79,11 +81,34 @@ class PolynomialReconstructionDescriptor
     m_row_weighting = row_weighting;
   }
 
+  PolynomialReconstructionDescriptor(const IntegrationMethodType integration_method, const size_t degree)
+    : m_integration_method{integration_method},
+      m_degree{degree},
+      m_stencil_descriptor(degree, StencilDescriptor::ConnectionType::by_nodes)
+  {}
+
+  PolynomialReconstructionDescriptor(const IntegrationMethodType integration_method,
+                                     const size_t degree,
+                                     const BoundaryDescriptorList& symmetry_boundary_descriptor_list)
+    : m_integration_method{integration_method},
+      m_degree{degree},
+      m_stencil_descriptor(degree, StencilDescriptor::ConnectionType::by_nodes),
+      m_symmetry_boundary_descriptor_list(symmetry_boundary_descriptor_list)
+  {}
+
+  PolynomialReconstructionDescriptor(const IntegrationMethodType integration_method,
+                                     const size_t degree,
+                                     const StencilDescriptor& stencil_descriptor)
+    : m_integration_method{integration_method}, m_degree{degree}, m_stencil_descriptor{stencil_descriptor}
+  {}
+
   PolynomialReconstructionDescriptor(const IntegrationMethodType integration_method,
                                      const size_t degree,
-                                     const BoundaryDescriptorList& symmetry_boundary_descriptor_list = {})
+                                     const StencilDescriptor& stencil_descriptor,
+                                     const BoundaryDescriptorList& symmetry_boundary_descriptor_list)
     : m_integration_method{integration_method},
       m_degree{degree},
+      m_stencil_descriptor{stencil_descriptor},
       m_symmetry_boundary_descriptor_list(symmetry_boundary_descriptor_list)
   {}
 
diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
deleted file mode 100644
index 629d1337ebca06a338dd7bf5bdd8417d2316250d..0000000000000000000000000000000000000000
--- a/src/scheme/test_reconstruction.cpp
+++ /dev/null
@@ -1,35 +0,0 @@
-#include <scheme/test_reconstruction.hpp>
-
-#include <mesh/NamedBoundaryDescriptor.hpp>
-#include <scheme/IntegrationMethodType.hpp>
-#include <scheme/PolynomialReconstruction.hpp>
-#include <utils/Timer.hpp>
-
-void
-test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
-                    uint64_t degree)
-{
-  std::vector descriptor_list = {
-    PolynomialReconstructionDescriptor{IntegrationMethodType::boundary,
-                                       degree,
-                                       {std::make_shared<NamedBoundaryDescriptor>("XMIN"),
-                                        std::make_shared<NamedBoundaryDescriptor>("XMAX"),
-                                        std::make_shared<NamedBoundaryDescriptor>("YMIN"),
-                                        std::make_shared<NamedBoundaryDescriptor>("YMAX")}}};
-
-  for (auto&& descriptor : descriptor_list) {
-    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";
-
-    PolynomialReconstruction reconstruction{descriptor};
-    Timer t;
-    for (size_t i = 0; i < nb_loops; ++i) {
-      auto res = reconstruction.build(discrete_function_variant_list);
-    }
-    t.pause();
-    std::cout << "t = " << t << '\n';
-  }
-
-  std::cout << "finished!\n";
-}
diff --git a/src/scheme/test_reconstruction.hpp b/src/scheme/test_reconstruction.hpp
deleted file mode 100644
index d2fc5bd031989eecdf4d837f566c137b3f874888..0000000000000000000000000000000000000000
--- a/src/scheme/test_reconstruction.hpp
+++ /dev/null
@@ -1,14 +0,0 @@
-#ifndef TEST_RECONSTRUCTION_HPP
-#define TEST_RECONSTRUCTION_HPP
-
-#warning REMOVE AFTER TESTS FINISHED
-#include <mesh/MeshVariant.hpp>
-#include <scheme/DiscreteFunctionVariant.hpp>
-
-#include <vector>
-
-void test_reconstruction(
-  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
-  uint64_t degree = 1);
-
-#endif   // TEST_RECONSTRUCTION_HPP
diff --git a/src/utils/GlobalVariableManager.hpp b/src/utils/GlobalVariableManager.hpp
index f720252fdd7c04b1ddd8e234d890bf5b1956c899..8b7e0a339fd91c3569901ba0d179b0713b6ea391 100644
--- a/src/utils/GlobalVariableManager.hpp
+++ b/src/utils/GlobalVariableManager.hpp
@@ -1,15 +1,23 @@
 #ifndef GLOBAL_VARIABLE_MANAGER_HPP
 #define GLOBAL_VARIABLE_MANAGER_HPP
 
+#include <utils/Exceptions.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 
+#include <optional>
+
 class GlobalVariableManager
 {
  private:
+  // Give some special access for testing
+  friend class NbGhostLayersTester;
+
   size_t m_connectivity_id = 0;
   size_t m_mesh_id         = 0;
 
+  std::optional<size_t> m_number_of_ghost_layers;
+
   static GlobalVariableManager* m_instance;
 
   explicit GlobalVariableManager()                    = default;
@@ -32,6 +40,24 @@ class GlobalVariableManager
     return m_mesh_id++;
   }
 
+  PUGS_INLINE
+  void
+  setNumberOfGhostLayers(const size_t number)
+  {
+    if (m_number_of_ghost_layers.has_value()) {
+      throw UnexpectedError("changing number of ghost layers is forbidden");
+    }
+    m_number_of_ghost_layers = number;
+  }
+
+  PUGS_INLINE
+  size_t
+  getNumberOfGhostLayers()
+  {
+    Assert(m_number_of_ghost_layers.has_value());
+    return m_number_of_ghost_layers.value();
+  }
+
   PUGS_INLINE
   static GlobalVariableManager&
   instance()
diff --git a/src/utils/Messenger.cpp b/src/utils/Messenger.cpp
index 14891477a2b804c4f0c8698fdfa51a8ea2470c8b..c66102669709be2f4d5aa19677b212611391ea16 100644
--- a/src/utils/Messenger.cpp
+++ b/src/utils/Messenger.cpp
@@ -28,7 +28,7 @@ Messenger::destroy()
   }
 }
 
-Messenger::Messenger([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[], bool parallel_output)
+Messenger::Messenger([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[], [[maybe_unused]] bool parallel_output)
 {
 #ifdef PUGS_HAS_MPI
   MPI_Init(&argc, &argv);
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index f89f543f92bc3a175e4f055b3286e4a45d985396..e19e6a29a694d4f9971683a99050c0e41f9ea190 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -129,6 +129,10 @@ initialize(int& argc, char* argv[])
     bool pause_on_error = false;
     app.add_flag("-p,--pause-on-error", pause_on_error, "Pause for debugging on unexpected error [default: false]");
 
+    int nb_ghost_layers = 1;
+    app.add_option("--number-of-ghost-layers", nb_ghost_layers, "Number of ghost layers of cells [default: 1]")
+      ->check(CLI::Range(0, std::numeric_limits<decltype(nb_ghost_layers)>::max()));
+
     bool reproducible_sums = true;
     app.add_flag("--reproducible-sums,!--no-reproducible-sums", reproducible_sums,
                  "Special treatment of array sums to ensure reproducibility [default: true]");
@@ -165,6 +169,7 @@ initialize(int& argc, char* argv[])
       CommunicatorManager::setSplitColor(mpi_split_color);
     }
 
+    GlobalVariableManager::instance().setNumberOfGhostLayers(nb_ghost_layers);
     ExecutionStatManager::getInstance().setPrint(print_exec_stat);
     BacktraceManager::setShow(show_backtrace);
     ConsoleManager::setShowPreamble(show_preamble);
diff --git a/src/utils/ReproducibleSumMPI.hpp b/src/utils/ReproducibleSumMPI.hpp
index d4b80eb5d32384ba14cf49ade607cbf66fcf5632..fcb9bdc1e39db434581937287b24d5b0c8bbaae5 100644
--- a/src/utils/ReproducibleSumMPI.hpp
+++ b/src/utils/ReproducibleSumMPI.hpp
@@ -43,6 +43,9 @@ allReduceReproducibleSum(const ReproducibleScalarSum<ArrayT, BoolArrayT>& s)
   BinType sum_bin = zero;
   MPI_Allreduce(&local_bin, &sum_bin, 1, mpi_bin_type, mpi_bin_sum_op, parallel::Messenger::getInstance().comm());
 
+  MPI_Type_free(&mpi_bin_type);
+  MPI_Op_free(&mpi_bin_sum_op);
+
   return sum_bin;
 }
 
@@ -68,6 +71,9 @@ allReduceReproducibleSum(const ReproducibleTinyVectorSum<ArrayT, BoolArrayT>& s)
   BinType sum_bin = zero;
   MPI_Allreduce(&local_bin, &sum_bin, 1, mpi_bin_type, mpi_bin_sum_op, parallel::Messenger::getInstance().comm());
 
+  MPI_Op_free(&mpi_bin_sum_op);
+  MPI_Type_free(&mpi_bin_type);
+
   return sum_bin;
 }
 
@@ -93,6 +99,9 @@ allReduceReproducibleSum(const ReproducibleTinyMatrixSum<ArrayT, BoolArrayT>& s)
   BinType sum_bin = zero;
   MPI_Allreduce(&local_bin, &sum_bin, 1, mpi_bin_type, mpi_bin_sum_op, parallel::Messenger::getInstance().comm());
 
+  MPI_Op_free(&mpi_bin_sum_op);
+  MPI_Type_free(&mpi_bin_type);
+
   return sum_bin;
 }
 
diff --git a/src/utils/ReproducibleSumUtils.hpp b/src/utils/ReproducibleSumUtils.hpp
index f2d52aa20122452e96130c0c29b6d12ccc32a4aa..2ce3d763cfdf8f8af2d2b1dcb975640dad6290b6 100644
--- a/src/utils/ReproducibleSumUtils.hpp
+++ b/src/utils/ReproducibleSumUtils.hpp
@@ -168,7 +168,7 @@ class ReproducibleScalarSum
   PUGS_INLINE static void
   _update(const DataType& m, Bin& bin) noexcept(NO_ASSERT)
   {
-    Assert(m > 0);
+    Assert(m >= 0);
 
     using namespace reproducible_sum_utils;
     if (m >= std::pow(DataType{2}, W - 1.) * ulp(bin.S[0])) {
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index f7e3a8d4f6b75fff1feda88bcdcaebbf81958755..19a1aaae48a3aaa092129062746841606706e23c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -162,6 +162,7 @@ add_executable (unit_tests
 add_executable (mpi_unit_tests
   mpi_test_main.cpp
   test_Connectivity.cpp
+  test_ConnectivityDispatcher.cpp
   test_DiscreteFunctionDPk.cpp
   test_DiscreteFunctionDPkVector.cpp
   test_DiscreteFunctionIntegrator.cpp
diff --git a/tests/mpi_test_main.cpp b/tests/mpi_test_main.cpp
index 657cccc8ce4a6129de59f1df111414fdd5236799..7947364942fff1b8213a63b1972f910f12dfe63d 100644
--- a/tests/mpi_test_main.cpp
+++ b/tests/mpi_test_main.cpp
@@ -103,6 +103,7 @@ main(int argc, char* argv[])
       StencilManager::create();
 
       GlobalVariableManager::create();
+      GlobalVariableManager::instance().setNumberOfGhostLayers(1);
 
       MeshDataBaseForTests::create();
 
diff --git a/tests/test_ConnectivityDispatcher.cpp b/tests/test_ConnectivityDispatcher.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..abe07cc21e21a1b34ae3f20dfe13ff7f52ac83b9
--- /dev/null
+++ b/tests/test_ConnectivityDispatcher.cpp
@@ -0,0 +1,214 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <mesh/CartesianMeshBuilder.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/GmshReader.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <utils/Messenger.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+
+#include <filesystem>
+
+// clazy:excludeall=non-pod-global-static
+
+class NbGhostLayersTester
+{
+ private:
+  const size_t m_original_number_of_ghost_layers;
+
+ public:
+  NbGhostLayersTester(const size_t number_of_ghost_layers)
+    : m_original_number_of_ghost_layers{GlobalVariableManager::instance().getNumberOfGhostLayers()}
+  {
+    GlobalVariableManager::instance().m_number_of_ghost_layers = number_of_ghost_layers;
+  }
+
+  ~NbGhostLayersTester()
+  {
+    GlobalVariableManager::instance().m_number_of_ghost_layers = m_original_number_of_ghost_layers;
+  }
+};
+
+TEST_CASE("ConnectivityDispatcher", "[mesh]")
+{
+  auto check_number_of_ghost_layers = [](const auto& connectivity, const size_t number_of_layers) {
+    // We assume that the specified number of layers can be built
+    // (there are enough non owned layer of cells in the connectivity)
+    const auto cell_is_owned = connectivity.cellIsOwned();
+
+    CellValue<size_t> cell_layer{connectivity};
+    cell_layer.fill(number_of_layers + 1);
+
+    NodeValue<size_t> node_layer{connectivity};
+    node_layer.fill(number_of_layers + 1);
+
+    auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+    auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+
+    parallel_for(
+      connectivity.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        if (cell_is_owned[cell_id]) {
+          cell_layer[cell_id] = 0;
+        }
+      });
+
+    for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) {
+      parallel_for(
+        connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+          auto node_cell_list = node_to_cell_matrix[node_id];
+          size_t min_layer    = cell_layer[node_cell_list[0]];
+          for (size_t i_cell = 1; i_cell < node_cell_list.size(); ++i_cell) {
+            min_layer = std::min(min_layer, cell_layer[node_cell_list[i_cell]]);
+          }
+          if (min_layer < number_of_layers + 1) {
+            node_layer[node_id] = min_layer;
+          }
+        });
+
+      parallel_for(
+        connectivity.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+          auto cell_node_list = cell_to_node_matrix[cell_id];
+          size_t min_layer    = node_layer[cell_node_list[0]];
+          size_t max_layer    = min_layer;
+          for (size_t i_node = 1; i_node < cell_node_list.size(); ++i_node) {
+            min_layer = std::min(min_layer, node_layer[cell_node_list[i_node]]);
+            max_layer = std::max(max_layer, node_layer[cell_node_list[i_node]]);
+          }
+          if ((min_layer != max_layer) or
+              ((min_layer < number_of_layers + 1) and (cell_layer[cell_id] == number_of_layers + 1))) {
+            cell_layer[cell_id] = min_layer + 1;
+          }
+        });
+    }
+
+    auto is_boundary_face    = connectivity.isBoundaryFace();
+    auto face_to_cell_matrix = connectivity.faceToCellMatrix();
+
+    bool has_required_number_of_ghost_layers = true;
+    for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+      auto face_cell_list = face_to_cell_matrix[face_id];
+      if ((face_cell_list.size() == 1) and (not is_boundary_face[face_id])) {
+        const CellId face_cell_id = face_cell_list[0];
+        has_required_number_of_ghost_layers &= (cell_layer[face_cell_id] == number_of_layers);
+      }
+    }
+
+    REQUIRE(parallel::allReduceAnd(has_required_number_of_ghost_layers));
+    bool first_ghost_layer_is_1 = true;
+    for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+      auto face_cell_list = face_to_cell_matrix[face_id];
+      if (face_cell_list.size() == 2) {
+        const CellId face_cell0_id = face_cell_list[0];
+        const CellId face_cell1_id = face_cell_list[1];
+        if (cell_is_owned[face_cell0_id] xor cell_is_owned[face_cell1_id]) {
+          for (size_t i_cell = 0; i_cell < face_cell_list.size(); ++i_cell) {
+            const CellId face_cell_id = face_cell_list[i_cell];
+            if (not cell_is_owned[face_cell_id]) {
+              first_ghost_layer_is_1 &= (cell_layer[face_cell_id] == 1);
+            }
+          }
+        }
+      }
+    }
+
+    REQUIRE(parallel::allReduceAnd(first_ghost_layer_is_1));
+  };
+
+  SECTION("1 layer meshes")
+  {
+    for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+      SECTION(named_mesh.name())
+      {
+        const std::shared_ptr p_mesh = named_mesh.mesh()->get<const Mesh<1>>();
+        check_number_of_ghost_layers(p_mesh->connectivity(), 1);
+      }
+    }
+    for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+      SECTION(named_mesh.name())
+      {
+        const std::shared_ptr p_mesh = named_mesh.mesh()->get<const Mesh<2>>();
+        check_number_of_ghost_layers(p_mesh->connectivity(), 1);
+      }
+    }
+    for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+      SECTION(named_mesh.name())
+      {
+        const std::shared_ptr p_mesh = named_mesh.mesh()->get<const Mesh<3>>();
+        check_number_of_ghost_layers(p_mesh->connectivity(), 1);
+      }
+    }
+  }
+
+  for (size_t nb_ghost_layers = 2; nb_ghost_layers < 5; ++nb_ghost_layers) {
+    std::stringstream os;
+    os << nb_ghost_layers << " layer meshes";
+
+    SECTION(os.str())
+    {
+      REQUIRE(GlobalVariableManager::instance().getNumberOfGhostLayers() == 1);
+
+      NbGhostLayersTester nb_ghost_layers_tester(nb_ghost_layers);
+
+      REQUIRE(GlobalVariableManager::instance().getNumberOfGhostLayers() == nb_ghost_layers);
+
+      SECTION("Cartesian 1D mesh")
+      {
+        auto cartesian_1d_mesh =
+          CartesianMeshBuilder{TinyVector<1>{-1}, TinyVector<1>{3}, TinyVector<1, size_t>{23}}.mesh();
+        const std::shared_ptr p_mesh = cartesian_1d_mesh->get<const Mesh<1>>();
+        check_number_of_ghost_layers(p_mesh->connectivity(), nb_ghost_layers);
+      }
+
+      SECTION("Cartesian 2D mesh")
+      {
+        auto cartesian_2d_mesh =
+          CartesianMeshBuilder{TinyVector<2>{0, -1}, TinyVector<2>{3, 2}, TinyVector<2, size_t>{6, 7}}.mesh();
+        const std::shared_ptr p_mesh = cartesian_2d_mesh->get<const Mesh<2>>();
+        check_number_of_ghost_layers(p_mesh->connectivity(), nb_ghost_layers);
+      }
+
+      SECTION("Cartesian 3D mesh")
+      {
+        auto cartesian_3d_mesh =
+          CartesianMeshBuilder{TinyVector<3>{0, 1, 0}, TinyVector<3>{2, -1, 3}, TinyVector<3, size_t>{6, 7, 4}}.mesh();
+        const std::shared_ptr p_mesh = cartesian_3d_mesh->get<const Mesh<3>>();
+        check_number_of_ghost_layers(p_mesh->connectivity(), nb_ghost_layers);
+      }
+
+      SECTION("unordered 1d mesh")
+      {
+        const std::string filename = std::filesystem::path{PUGS_BINARY_DIR}.append("tests").append("unordered-1d.msh");
+
+        auto mesh_v = GmshReader{filename}.mesh();
+
+        const std::shared_ptr p_mesh = mesh_v->get<const Mesh<1>>();
+        check_number_of_ghost_layers(p_mesh->connectivity(), nb_ghost_layers);
+      }
+
+      SECTION("hybrid 2d mesh")
+      {
+        const std::string filename = std::filesystem::path{PUGS_BINARY_DIR}.append("tests").append("hybrid-2d.msh");
+
+        auto mesh_v = GmshReader{filename}.mesh();
+
+        const std::shared_ptr p_mesh = mesh_v->get<const Mesh<2>>();
+        check_number_of_ghost_layers(p_mesh->connectivity(), nb_ghost_layers);
+      }
+
+      SECTION("hybrid 3d mesh")
+      {
+        const std::string filename = std::filesystem::path{PUGS_BINARY_DIR}.append("tests").append("hybrid-3d.msh");
+
+        auto mesh_v = GmshReader{filename}.mesh();
+
+        const std::shared_ptr p_mesh = mesh_v->get<const Mesh<3>>();
+        check_number_of_ghost_layers(p_mesh->connectivity(), nb_ghost_layers);
+      }
+    }
+
+    REQUIRE(GlobalVariableManager::instance().getNumberOfGhostLayers() == 1);
+  }
+}
diff --git a/tests/test_PolynomialReconstructionDescriptor.cpp b/tests/test_PolynomialReconstructionDescriptor.cpp
index cadc6aaf4904d502ba77396366282fbe89ba1e27..32d3c4af12c48d6c7f69b2b250838bea7e3cf19e 100644
--- a/tests/test_PolynomialReconstructionDescriptor.cpp
+++ b/tests/test_PolynomialReconstructionDescriptor.cpp
@@ -2,18 +2,23 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_all.hpp>
 
+#include <mesh/NamedBoundaryDescriptor.hpp>
+#include <mesh/NumberedBoundaryDescriptor.hpp>
 #include <scheme/PolynomialReconstructionDescriptor.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
 TEST_CASE("PolynomialReconstructionDescriptor", "[scheme]")
 {
-#warning tests are not completed
   SECTION("degree 1")
   {
     PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::cell_center, 1};
 
     REQUIRE(descriptor.degree() == 1);
+    REQUIRE(descriptor.stencilDescriptor().numberOfLayers() == 1);
+    REQUIRE(descriptor.stencilDescriptor().connectionType() == StencilDescriptor::ConnectionType::by_nodes);
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList().size() == 0);
+
     REQUIRE(descriptor.preconditioning() == true);
     REQUIRE(descriptor.rowWeighting() == true);
   }
@@ -23,32 +28,82 @@ TEST_CASE("PolynomialReconstructionDescriptor", "[scheme]")
     PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::cell_center, 4};
 
     REQUIRE(descriptor.degree() == 4);
+    REQUIRE(descriptor.stencilDescriptor().numberOfLayers() == 4);
+    REQUIRE(descriptor.stencilDescriptor().connectionType() == StencilDescriptor::ConnectionType::by_nodes);
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList().size() == 0);
+
     REQUIRE(descriptor.preconditioning() == true);
     REQUIRE(descriptor.rowWeighting() == true);
   }
 
-  SECTION("utlities")
+  SECTION("degree and stencil")
   {
-    PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::cell_center, 3};
+    StencilDescriptor sd{2, StencilDescriptor::ConnectionType::by_faces};
+
+    PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::cell_center, 1, sd};
+
+    REQUIRE(descriptor.degree() == 1);
+    REQUIRE(descriptor.stencilDescriptor().numberOfLayers() == 2);
+    REQUIRE(descriptor.stencilDescriptor().connectionType() == StencilDescriptor::ConnectionType::by_faces);
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList().size() == 0);
 
-    REQUIRE(descriptor.degree() == 3);
     REQUIRE(descriptor.preconditioning() == true);
     REQUIRE(descriptor.rowWeighting() == true);
+  }
 
-    SECTION("set degree")
-    {
-      descriptor.setDegree(2);
+  SECTION("degree and symmetries")
+  {
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> bc_list;
+    bc_list.push_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+    bc_list.push_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
+    bc_list.push_back(std::make_shared<NumberedBoundaryDescriptor>(2));
 
-      REQUIRE(descriptor.degree() == 2);
-      REQUIRE(descriptor.preconditioning() == true);
-      REQUIRE(descriptor.rowWeighting() == true);
+    PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::cell_center, 2, bc_list};
 
-      descriptor.setDegree(4);
+    REQUIRE(descriptor.degree() == 2);
+    REQUIRE(descriptor.stencilDescriptor().numberOfLayers() == 2);
+    REQUIRE(descriptor.stencilDescriptor().connectionType() == StencilDescriptor::ConnectionType::by_nodes);
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList().size() == 3);
 
-      REQUIRE(descriptor.degree() == 4);
-      REQUIRE(descriptor.preconditioning() == true);
-      REQUIRE(descriptor.rowWeighting() == true);
-    }
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList()[0]->type() == IBoundaryDescriptor::Type::named);
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList()[1]->type() == IBoundaryDescriptor::Type::named);
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList()[2]->type() == IBoundaryDescriptor::Type::numbered);
+
+    REQUIRE(descriptor.preconditioning() == true);
+    REQUIRE(descriptor.rowWeighting() == true);
+  }
+
+  SECTION("degree, stencil and symmetries")
+  {
+    StencilDescriptor sd{3, StencilDescriptor::ConnectionType::by_edges};
+
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> bc_list;
+    bc_list.push_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+    bc_list.push_back(std::make_shared<NumberedBoundaryDescriptor>(2));
+    bc_list.push_back(std::make_shared<NamedBoundaryDescriptor>("YMIN"));
+
+    PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::cell_center, 1, sd, bc_list};
+
+    REQUIRE(descriptor.degree() == 1);
+    REQUIRE(descriptor.stencilDescriptor().numberOfLayers() == 3);
+    REQUIRE(descriptor.stencilDescriptor().connectionType() == StencilDescriptor::ConnectionType::by_edges);
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList().size() == 3);
+
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList()[0]->type() == IBoundaryDescriptor::Type::named);
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList()[1]->type() == IBoundaryDescriptor::Type::numbered);
+    REQUIRE(descriptor.symmetryBoundaryDescriptorList()[2]->type() == IBoundaryDescriptor::Type::named);
+
+    REQUIRE(descriptor.preconditioning() == true);
+    REQUIRE(descriptor.rowWeighting() == true);
+  }
+
+  SECTION("utlities")
+  {
+    PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::cell_center, 3};
+
+    REQUIRE(descriptor.degree() == 3);
+    REQUIRE(descriptor.preconditioning() == true);
+    REQUIRE(descriptor.rowWeighting() == true);
 
     SECTION("set preconditioning")
     {
diff --git a/tests/test_StencilBuilder.cpp b/tests/test_StencilBuilder.cpp
index 0ef8d97ba9ef690ca4420a49905881e5ccbc49b3..b99ba218aa1e3eb5d8ec884e974cfad2f4a157db 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,11 @@ TEST_CASE("StencilBuilder", "[mesh]")
 
         const Connectivity<1>& connectivity = mesh.connectivity();
 
-        Stencil stencil = StencilManager::instance().getStencil(connectivity);
+        auto stencil_array =
+          StencilManager::instance()
+            .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes});
 
-        REQUIRE(is_valid(connectivity, stencil));
+        REQUIRE(is_valid(connectivity, stencil_array));
       }
 
       SECTION("unordered")
@@ -73,9 +75,11 @@ 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, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes});
 
-        REQUIRE(is_valid(connectivity, stencil));
+        REQUIRE(is_valid(connectivity, stencil_array));
       }
     }
 
@@ -86,7 +90,11 @@ 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,
+                                                StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes})));
       }
 
       SECTION("hybrid")
@@ -94,7 +102,11 @@ 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,
+                                                StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes})));
       }
     }
 
@@ -105,7 +117,11 @@ 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,
+                                                StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes})));
       }
 
       SECTION("hybrid")
@@ -113,24 +129,30 @@ 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,
+                                                StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes})));
       }
     }
   }
 
   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].second[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 +160,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 +168,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) {
-        auto p_boundary_descriptor = stencil.symmetryBoundaryStencilList()[i_symmetry_stencil].first;
-        const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor;
+        const IBoundaryDescriptor& boundary_descriptor =
+          stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].boundaryDescriptor();
 
-        auto boundary_stencil   = stencil.symmetryBoundaryStencilList()[i_symmetry_stencil].second;
+        auto boundary_stencil   = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray();
         auto boundary_node_list = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
 
         CellValue<bool> boundary_cell{mesh.connectivity()};
@@ -205,6 +227,44 @@ TEST_CASE("StencilBuilder", "[mesh]")
       return is_valid;
     };
 
+    SECTION("1D")
+    {
+      StencilManager::BoundaryDescriptorList list;
+      list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN"));
+      list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX"));
+
+      SECTION("cartesian")
+      {
+        const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
+
+        const Connectivity<1>& connectivity = mesh.connectivity();
+
+        auto stencil_array =
+          StencilManager::instance()
+            .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
+                                       list);
+
+        REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2);
+        REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity));
+        REQUIRE(symmetry_stencils_are_valid(stencil_array, mesh));
+      }
+
+      SECTION("hybrid")
+      {
+        const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+        const Connectivity<1>& connectivity = mesh.connectivity();
+        auto stencil =
+          StencilManager::instance()
+            .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
+                                       list);
+
+        REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 2);
+        REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity));
+        REQUIRE(symmetry_stencils_are_valid(stencil, mesh));
+      }
+    }
+
     SECTION("2D")
     {
       StencilManager::BoundaryDescriptorList list;
@@ -219,11 +279,14 @@ 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, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
+                                       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 +294,12 @@ 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, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
+                                       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 +320,12 @@ 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, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
+                                       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 +335,12 @@ 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, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
+                                       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));
       }
diff --git a/tests/test_main.cpp b/tests/test_main.cpp
index 3f81068df53dc2708110b1dd719cad2009cdbfb3..504541ee0951895bd73899eb8aa0b5d89f7d46f6 100644
--- a/tests/test_main.cpp
+++ b/tests/test_main.cpp
@@ -62,6 +62,7 @@ main(int argc, char* argv[])
       StencilManager::create();
 
       GlobalVariableManager::create();
+      GlobalVariableManager::instance().setNumberOfGhostLayers(1);
 
       MeshDataBaseForTests::create();