diff --git a/src/algebra/ShrinkMatrixView.hpp b/src/algebra/ShrinkMatrixView.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a93efc07a2dd67103af7ba3bdb133f9fab3f59ec
--- /dev/null
+++ b/src/algebra/ShrinkMatrixView.hpp
@@ -0,0 +1,52 @@
+#ifndef SHRINK_MATRIX_VIEW_HPP
+#define SHRINK_MATRIX_VIEW_HPP
+
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <cstddef>
+
+template <typename MatrixType>
+class ShrinkMatrixView
+{
+ public:
+  using index_type = typename MatrixType::index_type;
+  using data_type  = typename MatrixType::data_type;
+
+ private:
+  MatrixType& m_matrix;
+  const size_t m_nb_rows;
+
+ public:
+  PUGS_INLINE size_t
+  numberOfRows() const noexcept
+  {
+    return m_nb_rows;
+  }
+
+  PUGS_INLINE size_t
+  numberOfColumns() const noexcept
+  {
+    return m_matrix.numberOfColumns();
+  }
+
+  PUGS_INLINE
+  data_type&
+  operator()(index_type i, index_type j) const noexcept(NO_ASSERT)
+  {
+    Assert(i < m_nb_rows and j < m_matrix.numberOfColumns(), "cannot access element: invalid indices");
+    return m_matrix(i, j);
+  }
+
+  ShrinkMatrixView(MatrixType& matrix, size_t nb_rows) noexcept(NO_ASSERT) : m_matrix{matrix}, m_nb_rows{nb_rows}
+  {
+    Assert(m_nb_rows <= matrix.numberOfRows(), "shrink number of rows must be smaller than original matrix's");
+  }
+
+  ShrinkMatrixView(const ShrinkMatrixView&) = delete;
+  ShrinkMatrixView(ShrinkMatrixView&&)      = delete;
+
+  ~ShrinkMatrixView() noexcept = default;
+};
+
+#endif   // SHRINK_MATRIX_VIEW_HPP
diff --git a/src/algebra/ShrinkVectorView.hpp b/src/algebra/ShrinkVectorView.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dc5faa800a1c0828b0569c33efee95add38f061a
--- /dev/null
+++ b/src/algebra/ShrinkVectorView.hpp
@@ -0,0 +1,46 @@
+#ifndef SHRINK_VECTOR_VIEW_HPP
+#define SHRINK_VECTOR_VIEW_HPP
+
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <cstddef>
+
+template <typename VectorType>
+class ShrinkVectorView
+{
+ public:
+  using index_type = typename VectorType::index_type;
+  using data_type  = typename VectorType::data_type;
+
+ private:
+  VectorType& m_vector;
+  const size_t m_dimension;
+
+ public:
+  PUGS_INLINE size_t
+  dimension() const noexcept
+  {
+    return m_dimension;
+  }
+
+  PUGS_INLINE
+  data_type&
+  operator[](index_type i) const noexcept(NO_ASSERT)
+  {
+    Assert(i < m_dimension, "cannot access element: invalid indices");
+    return m_vector[i];
+  }
+
+  ShrinkVectorView(VectorType& vector, size_t dimension) noexcept(NO_ASSERT) : m_vector{vector}, m_dimension{dimension}
+  {
+    Assert(m_dimension <= vector.dimension(), "shrink number of rows must be smaller than original vector's");
+  }
+
+  ShrinkVectorView(const ShrinkVectorView&) = delete;
+  ShrinkVectorView(ShrinkVectorView&&)      = delete;
+
+  ~ShrinkVectorView() noexcept = default;
+};
+
+#endif   // SHRINK_VECTOR_VIEW_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 7b404fcf670c5217721efba890728ce97e3e91a4..f96e96f13fc172b577f7f4e61cf77ae61518091a 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -43,6 +43,7 @@
 #include <scheme/OutflowBoundaryConditionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <scheme/VariableBCDescriptor.hpp>
+#include <scheme/test_reconstruction.hpp>
 #include <utils/Socket.hpp>
 
 #include <memory>
@@ -670,6 +671,28 @@ SchemeModule::SchemeModule()
 
                                                    ));
 
+#warning REMOVE AFTER TESTS FINISHED
+  this->_addBuiltinFunction("test_reconstruction",
+                            std::function(
+
+                              [](const std::shared_ptr<const MeshVariant> mesh_v,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function_v) -> void {
+                                test_reconstruction(mesh_v, discrete_function_v);
+                              }
+
+                              ));
+
+#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);
+                              }
+
+                              ));
+
   MathFunctionRegisterForVh{*this};
 }
 
diff --git a/src/main.cpp b/src/main.cpp
index 7d0112c80715b3cd0b69ea91b874527134b6a2ed..a9bec44f0e7c2912ec4e8a4cd8823b03928b5bbf 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -4,6 +4,7 @@
 #include <mesh/DualConnectivityManager.hpp>
 #include <mesh/DualMeshManager.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <mesh/StencilManager.hpp>
 #include <mesh/SynchronizerManager.hpp>
 #include <utils/ExecutionStatManager.hpp>
 #include <utils/GlobalVariableManager.hpp>
@@ -24,6 +25,7 @@ main(int argc, char* argv[])
   MeshDataManager::create();
   DualConnectivityManager::create();
   DualMeshManager::create();
+  StencilManager::create();
 
   GlobalVariableManager::create();
 
@@ -32,6 +34,7 @@ main(int argc, char* argv[])
 
   GlobalVariableManager::destroy();
 
+  StencilManager::destroy();
   DualMeshManager::destroy();
   DualConnectivityManager::destroy();
   MeshDataManager::destroy();
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index d03051cae50a4f2c39f95fbd9f9fb8e8cb06cdf9..b4d70359a5a4a7f8a4d2e352ab702ecc26761dbc 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -44,6 +44,7 @@ add_library(
   MeshUtils.cpp
   MeshVariant.cpp
   StencilBuilder.cpp
+  StencilManager.cpp
   SynchronizerManager.cpp
 )
 
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index bd1e0d85b6f94d1d2d00326857853f4da75eacfe..b28237ded4485de15e2b8c17233e95906752ab89 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -2,6 +2,7 @@
 
 #include <mesh/ConnectivityDescriptor.hpp>
 #include <mesh/ItemValueUtils.hpp>
+#include <mesh/StencilManager.hpp>
 #include <utils/GlobalVariableManager.hpp>
 #include <utils/Messenger.hpp>
 
@@ -310,6 +311,12 @@ Connectivity<Dimension>::_write(std::ostream& os) const
   return os;
 }
 
+template <size_t Dim>
+Connectivity<Dim>::~Connectivity()
+{
+  StencilManager::instance().deleteConnectivity(this->m_id);
+}
+
 template std::ostream& Connectivity<1>::_write(std::ostream&) const;
 template std::ostream& Connectivity<2>::_write(std::ostream&) const;
 template std::ostream& Connectivity<3>::_write(std::ostream&) const;
@@ -330,6 +337,10 @@ template Connectivity<1>::Connectivity();
 template Connectivity<2>::Connectivity();
 template Connectivity<3>::Connectivity();
 
+template Connectivity<1>::~Connectivity();
+template Connectivity<2>::~Connectivity();
+template Connectivity<3>::~Connectivity();
+
 template void Connectivity<1>::_buildFrom(const ConnectivityDescriptor&);
 template void Connectivity<2>::_buildFrom(const ConnectivityDescriptor&);
 template void Connectivity<3>::_buildFrom(const ConnectivityDescriptor&);
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 24d9ef656f0cdd0746db0ee2f27881ce5244c8b0..c04bdbdf9d215ec9324d5545dcf2a079d4590c7f 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -744,7 +744,7 @@ class Connectivity final : public IConnectivity
   void _buildFrom(const ConnectivityDescriptor& descriptor);
 
  public:
-  ~Connectivity() = default;
+  ~Connectivity();
 };
 
 template <size_t Dimension>
diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp
index a43dbef60cd76a280b8b3ab679c2d0c4e5b5406f..44313040017f24ee7d79a8dec35697b3ebae1bc7 100644
--- a/src/mesh/StencilBuilder.cpp
+++ b/src/mesh/StencilBuilder.cpp
@@ -104,7 +104,7 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar
 
 template <typename ConnectivityType>
 Stencil
-StencilBuilder::build(const ConnectivityType& connectivity) const
+StencilBuilder::_build(const ConnectivityType& connectivity) const
 {
   Array<const uint32_t> row_map        = this->_getRowMap(connectivity);
   Array<const uint32_t> column_indices = this->_getColumnIndices(connectivity, row_map);
@@ -112,6 +112,21 @@ StencilBuilder::build(const ConnectivityType& connectivity) const
   return ConnectivityMatrix{row_map, column_indices};
 }
 
-template Stencil StencilBuilder::build(const Connectivity<1>& connectivity) const;
-template Stencil StencilBuilder::build(const Connectivity<2>& connectivity) const;
-template Stencil StencilBuilder::build(const Connectivity<3>& connectivity) const;
+Stencil
+StencilBuilder::build(const IConnectivity& connectivity) const
+{
+  switch (connectivity.dimension()) {
+  case 1: {
+    return StencilBuilder::_build(dynamic_cast<const Connectivity<1>&>(connectivity));
+  }
+  case 2: {
+    return StencilBuilder::_build(dynamic_cast<const Connectivity<2>&>(connectivity));
+  }
+  case 3: {
+    return StencilBuilder::_build(dynamic_cast<const Connectivity<3>&>(connectivity));
+  }
+  default: {
+    throw UnexpectedError("invalid connectivity dimension");
+  }
+  }
+}
diff --git a/src/mesh/StencilBuilder.hpp b/src/mesh/StencilBuilder.hpp
index 3a42350795f0caf8e749bb6235874ea3fdb65eed..1743a2d719fd920d121b604e61732d9f6450197d 100644
--- a/src/mesh/StencilBuilder.hpp
+++ b/src/mesh/StencilBuilder.hpp
@@ -3,6 +3,7 @@
 
 #include <mesh/Stencil.hpp>
 
+class IConnectivity;
 class StencilBuilder
 {
  private:
@@ -13,10 +14,13 @@ class StencilBuilder
   Array<const uint32_t> _getColumnIndices(const ConnectivityType& connectivity,
                                           const Array<const uint32_t>& row_map) const;
 
- public:
   template <typename ConnectivityType>
-  Stencil build(const ConnectivityType& connectivity) const;
+  Stencil _build(const ConnectivityType& connectivity) const;
+
+  friend class StencilManager;
+  Stencil build(const IConnectivity& connectivity) const;
 
+ public:
   StencilBuilder()                      = default;
   StencilBuilder(const StencilBuilder&) = default;
   StencilBuilder(StencilBuilder&&)      = default;
diff --git a/src/mesh/StencilManager.cpp b/src/mesh/StencilManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..58190f002739ca0070ac14c35214598c4476d8d9
--- /dev/null
+++ b/src/mesh/StencilManager.cpp
@@ -0,0 +1,59 @@
+#include <mesh/StencilManager.hpp>
+
+#include <mesh/StencilBuilder.hpp>
+#include <utils/Exceptions.hpp>
+
+StencilManager* StencilManager::m_instance = nullptr;
+
+void
+StencilManager::create()
+{
+  Assert(m_instance == nullptr, "StencilManager is already created");
+  m_instance = new StencilManager;
+}
+
+void
+StencilManager::destroy()
+{
+  Assert(m_instance != nullptr, "StencilManager was not created");
+
+  // LCOV_EXCL_START
+  if (m_instance->m_connectivity_id_to_stencil_map.size() > 0) {
+    std::stringstream error;
+    error << ": some connectivities are still registered\n";
+    for (const auto& [id, stencil_p] : m_instance->m_connectivity_id_to_stencil_map) {
+      error << " - connectivity of id " << rang::fgB::magenta << id << rang::style::reset << '\n';
+    }
+    throw UnexpectedError(error.str());
+    // LCOV_EXCL_STOP
+  }
+  delete m_instance;
+  m_instance = nullptr;
+}
+
+const Stencil&
+StencilManager::getStencil(const IConnectivity& connectivity)
+{
+  if (not m_connectivity_id_to_stencil_map.contains(connectivity.id())) {
+    m_connectivity_id_to_stencil_map[connectivity.id()] =
+      std::make_shared<Stencil>(StencilBuilder{}.build(connectivity));
+  }
+
+  return *m_connectivity_id_to_stencil_map.at(connectivity.id());
+}
+
+void
+StencilManager::deleteConnectivity(const size_t connectivity_id)
+{
+  bool has_removed = false;
+  do {
+    has_removed = false;
+    for (const auto& [id, p_stencil] : m_connectivity_id_to_stencil_map) {
+      if (connectivity_id == id) {
+        m_connectivity_id_to_stencil_map.erase(connectivity_id);
+        has_removed = true;
+        break;
+      }
+    }
+  } while (has_removed);
+}
diff --git a/src/mesh/StencilManager.hpp b/src/mesh/StencilManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d1503e60f64eb44e18ed03c228d5c754afcbe248
--- /dev/null
+++ b/src/mesh/StencilManager.hpp
@@ -0,0 +1,40 @@
+#ifndef STENCIL_MANAGER_HPP
+#define STENCIL_MANAGER_HPP
+
+#include <mesh/IConnectivity.hpp>
+#include <mesh/Stencil.hpp>
+
+#include <memory>
+#include <unordered_map>
+
+class StencilManager
+{
+ private:
+  StencilManager()  = default;
+  ~StencilManager() = default;
+
+  static StencilManager* m_instance;
+
+  std::unordered_map<size_t, std::shared_ptr<const Stencil>> m_connectivity_id_to_stencil_map;
+
+ public:
+  static void create();
+  static void destroy();
+
+  PUGS_INLINE
+  static StencilManager&
+  instance()
+  {
+    Assert(m_instance != nullptr, "StencilManager was not created!");
+    return *m_instance;
+  }
+
+  void deleteConnectivity(const size_t connectivity_id);
+
+  const Stencil& getStencil(const IConnectivity& i_connectivity);
+
+  StencilManager(const StencilManager&) = delete;
+  StencilManager(StencilManager&&)      = delete;
+};
+
+#endif   // STENCIL_MANAGER_HPP
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index c093c34f1608e91fda3c2a8ba4011e21648b4381..a759fca6af69a2f2b3e1622adef687443f15f0d1 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -10,6 +10,7 @@ add_library(
   DiscreteFunctionVectorIntegrator.cpp
   DiscreteFunctionVectorInterpoler.cpp
   FluxingAdvectionSolver.cpp
+  PolynomialReconstruction.cpp
 
   test_reconstruction.cpp
 )
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fc9b88d0e76aba7b030712be95a0ad14e0855943
--- /dev/null
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -0,0 +1,558 @@
+#include <scheme/PolynomialReconstruction.hpp>
+
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+class MutableDiscreteFunctionDPkVariant
+{
+ public:
+  using Variant = std::variant<DiscreteFunctionDPk<1, double>,
+                               DiscreteFunctionDPk<1, TinyVector<1>>,
+                               DiscreteFunctionDPk<1, TinyVector<2>>,
+                               DiscreteFunctionDPk<1, TinyVector<3>>,
+                               DiscreteFunctionDPk<1, TinyMatrix<1>>,
+                               DiscreteFunctionDPk<1, TinyMatrix<2>>,
+                               DiscreteFunctionDPk<1, TinyMatrix<3>>,
+
+                               DiscreteFunctionDPk<2, double>,
+                               DiscreteFunctionDPk<2, TinyVector<1>>,
+                               DiscreteFunctionDPk<2, TinyVector<2>>,
+                               DiscreteFunctionDPk<2, TinyVector<3>>,
+                               DiscreteFunctionDPk<2, TinyMatrix<1>>,
+                               DiscreteFunctionDPk<2, TinyMatrix<2>>,
+                               DiscreteFunctionDPk<2, TinyMatrix<3>>,
+
+                               DiscreteFunctionDPk<3, double>,
+                               DiscreteFunctionDPk<3, TinyVector<1>>,
+                               DiscreteFunctionDPk<3, TinyVector<2>>,
+                               DiscreteFunctionDPk<3, TinyVector<3>>,
+                               DiscreteFunctionDPk<3, TinyMatrix<1>>,
+                               DiscreteFunctionDPk<3, TinyMatrix<2>>,
+                               DiscreteFunctionDPk<3, TinyMatrix<3>>>;
+
+#warning add DiscreteFunctionDPkVector to the variant
+
+ private:
+  Variant m_mutable_discrete_function_dpk;
+
+ public:
+  PUGS_INLINE
+  const Variant&
+  mutableDiscreteFunctionDPk() const
+  {
+    return m_mutable_discrete_function_dpk;
+  }
+
+  template <typename DiscreteFunctionDPkT>
+  PUGS_INLINE auto&&
+  get() const
+  {
+    static_assert(is_discrete_function_dpk_v<DiscreteFunctionDPkT>, "invalid template argument");
+
+#ifndef NDEBUG
+    if (not std::holds_alternative<DiscreteFunctionDPkT>(this->m_mutable_discrete_function_dpk)) {
+      std::ostringstream error_msg;
+      error_msg << "invalid discrete function type\n";
+      error_msg << "- required " << rang::fgB::red << demangle<DiscreteFunctionDPkT>() << rang::fg::reset << '\n';
+      error_msg << "- contains " << rang::fgB::yellow
+                << std::visit([](auto&& f) -> std::string { return demangle<decltype(f)>(); },
+                              this->m_mutable_discrete_function_dpk)
+                << rang::fg::reset;
+      throw NormalError(error_msg.str());
+    }
+#endif   // NDEBUG
+
+    return std::get<DiscreteFunctionDPkT>(this->mutableDiscreteFunctionDPk());
+  }
+
+  template <size_t Dimension, typename DataType>
+  MutableDiscreteFunctionDPkVariant(const DiscreteFunctionDPk<Dimension, DataType>& discrete_function_dpk)
+    : m_mutable_discrete_function_dpk{DiscreteFunctionDPk<Dimension, DataType>{discrete_function_dpk}}
+  {
+    static_assert(std::is_same_v<std::remove_const_t<DataType>, double> or                       //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<1, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<2, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<3, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<1, 1, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<2, 2, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<3, 3, double>>,
+                  "DiscreteFunctionDPk with this DataType is not allowed in variant");
+  }
+
+  MutableDiscreteFunctionDPkVariant& operator=(MutableDiscreteFunctionDPkVariant&&)      = default;
+  MutableDiscreteFunctionDPkVariant& operator=(const MutableDiscreteFunctionDPkVariant&) = default;
+
+  MutableDiscreteFunctionDPkVariant(const MutableDiscreteFunctionDPkVariant&) = default;
+  MutableDiscreteFunctionDPkVariant(MutableDiscreteFunctionDPkVariant&&)      = default;
+
+  MutableDiscreteFunctionDPkVariant()  = delete;
+  ~MutableDiscreteFunctionDPkVariant() = default;
+};
+
+template <MeshConcept MeshType>
+[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
+PolynomialReconstruction::_build(
+  const std::shared_ptr<const MeshType>& p_mesh,
+  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+  const MeshType& mesh = *p_mesh;
+
+  using Rd = TinyVector<MeshType::Dimension>;
+
+  const size_t number_of_columns = [&]() {
+    size_t n = 0;
+    for (auto discrete_function_variant_p : discrete_function_variant_list) {
+      n += std::visit(
+        [](auto&& discrete_function) -> size_t {
+          using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+          if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+            using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
+            if constexpr (std::is_same_v<data_type, double>) {
+              return 1;
+            } else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
+              return data_type::Dimension;
+            } else {
+              throw UnexpectedError("unexpected data type " + demangle<data_type>());
+            }
+          } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+            return discrete_function.size();
+          } else {
+            throw UnexpectedError("unexpected discrete function type");
+          }
+        },
+        discrete_function_variant_p->discreteFunction());
+    }
+    return n;
+  }();
+
+  const size_t degree    = 1;
+  const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity());
+
+  auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+  auto cell_is_owned = mesh.connectivity().cellIsOwned();
+
+  const size_t max_stencil_size = [&]() {
+    size_t max_size = 0;
+    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+      auto stencil_cell_list = stencil[cell_id];
+      if (cell_is_owned[cell_id] and stencil_cell_list.size() > max_size) {
+        max_size = stencil_cell_list.size();
+      }
+    }
+    return max_size;
+  }();
+
+  Kokkos::Experimental::UniqueToken<Kokkos::DefaultExecutionSpace::execution_space,
+                                    Kokkos::Experimental::UniqueTokenScope::Global>
+    tokens;
+
+  std::vector<MutableDiscreteFunctionDPkVariant> mutable_discrete_function_dpk_variant_list;
+  for (size_t i_discrete_function_variant = 0; i_discrete_function_variant < discrete_function_variant_list.size();
+       ++i_discrete_function_variant) {
+    auto discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
+
+    std::visit(
+      [&](auto&& discrete_function) {
+        using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+        if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+          using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
+          mutable_discrete_function_dpk_variant_list.push_back(
+            DiscreteFunctionDPk<MeshType::Dimension, DataType>(p_mesh, degree));
+        } else {
+          throw NotImplementedError("discrete function type");
+        }
+      },
+      discrete_function_variant->discreteFunction());
+  }
+
+  SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
+  SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency());
+  SmallArray<SmallMatrix<double>> X_pool(Kokkos::DefaultExecutionSpace::concurrency());
+
+  for (size_t i = 0; i < A_pool.size(); ++i) {
+    A_pool[i] = SmallMatrix<double>(max_stencil_size, MeshType::Dimension);
+    B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns);
+    X_pool[i] = SmallMatrix<double>(MeshType::Dimension, number_of_columns);
+  }
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
+      if (cell_is_owned[cell_j_id]) {
+        const int32_t t = tokens.acquire();
+
+        auto stencil_cell_list = stencil[cell_j_id];
+        ShrinkMatrixView B(B_pool[t], stencil_cell_list.size());
+
+        size_t column_begin = 0;
+        for (size_t i_discrete_function_variant = 0;
+             i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
+          auto discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
+
+          std::visit(
+            [&](auto&& discrete_function) {
+              using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+              if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+                using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+
+                const DataType& pj = discrete_function[cell_j_id];
+                for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+                  const CellId cell_i_id = stencil_cell_list[i];
+                  const DataType& pi_pj  = discrete_function[cell_i_id] - pj;
+                  if constexpr (std::is_arithmetic_v<DataType>) {
+                    B(i, column_begin) = pi_pj;
+                  } else if constexpr (is_tiny_vector_v<DataType>) {
+                    for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
+                      B(i, kB) = pi_pj[k];
+                    }
+                  } else if constexpr (is_tiny_matrix_v<DataType>) {
+                    for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+                      for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+                        B(i, column_begin + k * DataType::NumberOfColumns + l) = pi_pj(k, l);
+                      }
+                    }
+                  }
+                }
+
+                if constexpr (std::is_arithmetic_v<DataType>) {
+                  ++column_begin;
+                } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
+                  column_begin += DataType::Dimension;
+                }
+              }
+            },
+            discrete_function_variant->discreteFunction());
+        }
+
+        ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
+
+        const Rd& Xj = xj[cell_j_id];
+
+        for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+          const CellId cell_i_id = stencil_cell_list[i];
+          const Rd Xi_Xj         = xj[cell_i_id] - Xj;
+          for (size_t l = 0; l < MeshType::Dimension; ++l) {
+            A(i, l) = Xi_Xj[l];
+          }
+        }
+
+        SmallMatrix<double> X = X_pool[t];
+        Givens::solveCollectionInPlace(A, X, B);
+
+        column_begin = 0;
+        for (size_t i_discrete_function_variant = 0;
+             i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
+          auto discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
+
+          std::visit(
+            [&](auto&& discrete_function) {
+              using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+              if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+                using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+
+                const DataType& pj         = discrete_function[cell_j_id];
+                auto discrete_function_dpk = mutable_discrete_function_dpk_variant_list[i_discrete_function_variant]
+                                               .get<DiscreteFunctionDPk<MeshType::Dimension, DataType>>();
+                auto dpk_j = discrete_function_dpk.coefficients(cell_j_id);
+                dpk_j[0]   = pj;
+
+                if constexpr (std::is_arithmetic_v<DataType>) {
+                  for (size_t i = 0; i < MeshType::Dimension; ++i) {
+                    auto& dpk_j_ip1 = dpk_j[i + 1];
+                    dpk_j_ip1       = X(i, column_begin);
+                  }
+                  ++column_begin;
+                } else if constexpr (is_tiny_vector_v<DataType>) {
+                  for (size_t i = 0; i < MeshType::Dimension; ++i) {
+                    auto& dpk_j_ip1 = dpk_j[i + 1];
+                    for (size_t k = 0; k < DataType::Dimension; ++k) {
+                      dpk_j_ip1[k] = X(i, column_begin + k);
+                    }
+                  }
+                  column_begin += DataType::Dimension;
+                } else if constexpr (is_tiny_matrix_v<DataType>) {
+                  for (size_t i = 0; i < MeshType::Dimension; ++i) {
+                    auto& dpk_j_ip1 = dpk_j[i + 1];
+                    for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+                      for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+                        dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
+                      }
+                    }
+                  }
+                  column_begin += DataType::Dimension;
+                }
+
+              } else {
+                throw NotImplementedError("discrete function type");
+              }
+            },
+            discrete_function_variant->discreteFunction());
+        }
+
+        tokens.release(t);
+      }
+    });
+
+  std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> discrete_function_dpk_variant_list;
+
+  // for (auto discrete_function_dpk_variant_p : mutable_discrete_function_dpk_variant_list) {
+  //   const_discrete_function_dpk_variant_list.push_back(discrete_function_dpk_variant_p);
+  // }
+
+  return discrete_function_dpk_variant_list;
+}
+
+std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
+PolynomialReconstruction::build(
+  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+  if (not hasSameMesh(discrete_function_variant_list)) {
+    throw NormalError("cannot reconstruct functions living of different mesh simultaneously");
+  }
+
+  auto mesh_v = getCommonMesh(discrete_function_variant_list);
+
+  return std::visit([&](auto&& p_mesh) { return this->_build(p_mesh, discrete_function_variant_list); },
+                    mesh_v->variant());
+}
+
+/**************************************/
+
+template <MeshConcept MeshType>
+[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
+PolynomialReconstruction::_build2(
+  const std::shared_ptr<const MeshType>& p_mesh,
+  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+  const MeshType& mesh = *p_mesh;
+
+  using Rd = TinyVector<MeshType::Dimension>;
+
+  const size_t number_of_columns = [&]() {
+    size_t n = 0;
+    for (auto discrete_function_variant_p : discrete_function_variant_list) {
+      n += std::visit(
+        [](auto&& discrete_function) -> size_t {
+          using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+          if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+            using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
+            if constexpr (std::is_same_v<data_type, double>) {
+              return 1;
+            } else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
+              return data_type::Dimension;
+            } else {
+              throw UnexpectedError("unexpected data type " + demangle<data_type>());
+            }
+          } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+            return discrete_function.size();
+          } else {
+            throw UnexpectedError("unexpected discrete function type");
+          }
+        },
+        discrete_function_variant_p->discreteFunction());
+    }
+    return n;
+  }();
+
+  const size_t degree    = 1;
+  const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity());
+
+  Array<const double*> value_begining(discrete_function_variant_list.size());
+  Array<size_t> value_size(discrete_function_variant_list.size());
+
+  for (size_t i = 0; i < discrete_function_variant_list.size(); ++i) {
+    DiscreteFunctionVariant discrete_function_v = *(discrete_function_variant_list[i]);
+    std::visit(
+      [=](auto&& discrete_function) {
+        using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+        if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+          using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
+          if constexpr (std::is_same_v<data_type, double>) {
+            value_begining[i] = &(discrete_function[CellId{0}]);
+            value_size[i]     = 1;
+          } else if constexpr (is_tiny_vector_v<data_type>) {
+            value_begining[i] = &(discrete_function[CellId{0}][0]);
+            value_size[i]     = data_type::Dimension;
+          } else if constexpr (is_tiny_matrix_v<data_type>) {
+            value_begining[i] = &(discrete_function[CellId{0}](0, 0));
+            value_size[i]     = data_type::Dimension;
+          } else {
+            throw UnexpectedError("unexpected data type " + demangle<data_type>());
+          }
+        } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+          throw NotImplementedError("Discrete Function P0 Vector");
+        } else {
+          throw UnexpectedError("unexpected discrete function type");
+        }
+      },
+      discrete_function_v.discreteFunction());
+  }
+
+  auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+  auto cell_is_owned = mesh.connectivity().cellIsOwned();
+
+  const size_t max_stencil_size = [&]() {
+    size_t max_size = 0;
+    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+      auto stencil_cell_list = stencil[cell_id];
+      if (cell_is_owned[cell_id] and stencil_cell_list.size() > max_size) {
+        max_size = stencil_cell_list.size();
+      }
+    }
+    return max_size;
+  }();
+
+  Kokkos::Experimental::UniqueToken<Kokkos::DefaultExecutionSpace::execution_space,
+                                    Kokkos::Experimental::UniqueTokenScope::Global>
+    tokens;
+
+  std::vector<MutableDiscreteFunctionDPkVariant> mutable_discrete_function_dpk_variant_list;
+  for (size_t i_discrete_function_variant = 0; i_discrete_function_variant < discrete_function_variant_list.size();
+       ++i_discrete_function_variant) {
+    auto discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
+
+    std::visit(
+      [&](auto&& discrete_function) {
+        using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+        if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+          using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
+          mutable_discrete_function_dpk_variant_list.push_back(
+            DiscreteFunctionDPk<MeshType::Dimension, DataType>(p_mesh, degree));
+        } else {
+          throw NotImplementedError("discrete function type");
+        }
+      },
+      discrete_function_variant->discreteFunction());
+  }
+
+  SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
+  SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency());
+  SmallArray<SmallMatrix<double>> X_pool(Kokkos::DefaultExecutionSpace::concurrency());
+
+  for (size_t i = 0; i < A_pool.size(); ++i) {
+    A_pool[i] = SmallMatrix<double>(max_stencil_size, MeshType::Dimension);
+    B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns);
+    X_pool[i] = SmallMatrix<double>(MeshType::Dimension, number_of_columns);
+  }
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
+      if (cell_is_owned[cell_j_id]) {
+        const int32_t t = tokens.acquire();
+
+        auto stencil_cell_list = stencil[cell_j_id];
+        ShrinkMatrixView B(B_pool[t], stencil_cell_list.size());
+
+        size_t i_column = 0;
+
+        for (size_t i_value = 0; i_value < value_begining.size(); ++i_value) {
+          const size_t i_value_size = value_size[i_value];
+
+          const size_t cell_j_value_index = cell_j_id * i_value_size;
+
+          for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+            const CellId cell_i_id          = stencil_cell_list[i];
+            const size_t cell_i_value_index = cell_i_id * i_value_size;
+            size_t j_index                  = cell_j_value_index;
+            size_t i_index                  = cell_i_value_index;
+
+            for (size_t shift = 0; shift < i_value_size; ++shift, ++i_index, ++j_index) {
+              B(i, i_column++) = value_begining[i_value][i_index] - value_begining[i_value][j_index];
+            }
+          }
+        }
+
+        ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
+
+        const Rd& Xj = xj[cell_j_id];
+
+        for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+          const CellId cell_i_id = stencil_cell_list[i];
+          const Rd Xi_Xj         = xj[cell_i_id] - Xj;
+          for (size_t l = 0; l < MeshType::Dimension; ++l) {
+            A(i, l) = Xi_Xj[l];
+          }
+        }
+
+        //        SmallMatrix<double> X = X_pool[t];
+
+        TinyMatrix<MeshType::Dimension, 1> X;
+        Givens::solveCollectionInPlace(A, X, B);
+
+        size_t column_begin = 0;
+        for (size_t i_discrete_function_variant = 0;
+             i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
+          auto discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
+
+          std::visit(
+            [&](auto&& discrete_function) {
+              using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+              if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+                using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+
+                const DataType& pj         = discrete_function[cell_j_id];
+                auto discrete_function_dpk = mutable_discrete_function_dpk_variant_list[i_discrete_function_variant]
+                                               .get<DiscreteFunctionDPk<MeshType::Dimension, DataType>>();
+                auto dpk_j = discrete_function_dpk.coefficients(cell_j_id);
+                dpk_j[0]   = pj;
+
+                if constexpr (std::is_arithmetic_v<DataType>) {
+                  for (size_t i = 0; i < MeshType::Dimension; ++i) {
+                    auto& dpk_j_ip1 = dpk_j[i + 1];
+                    dpk_j_ip1       = X(i, column_begin);
+                  }
+                  ++column_begin;
+                } else if constexpr (is_tiny_vector_v<DataType>) {
+                  for (size_t i = 0; i < MeshType::Dimension; ++i) {
+                    auto& dpk_j_ip1 = dpk_j[i + 1];
+                    for (size_t k = 0; k < DataType::Dimension; ++k) {
+                      dpk_j_ip1[k] = X(i, column_begin + k);
+                    }
+                  }
+                  column_begin += DataType::Dimension;
+                } else if constexpr (is_tiny_matrix_v<DataType>) {
+                  for (size_t i = 0; i < MeshType::Dimension; ++i) {
+                    auto& dpk_j_ip1 = dpk_j[i + 1];
+                    for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+                      for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+                        dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
+                      }
+                    }
+                  }
+                  column_begin += DataType::Dimension;
+                }
+
+              } else {
+                throw NotImplementedError("discrete function type");
+              }
+            },
+            discrete_function_variant->discreteFunction());
+        }
+
+        tokens.release(t);
+      }
+    });
+
+  std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> discrete_function_dpk_variant_list;
+
+  // for (auto discrete_function_dpk_variant_p : mutable_discrete_function_dpk_variant_list) {
+  //   const_discrete_function_dpk_variant_list.push_back(discrete_function_dpk_variant_p);
+  // }
+
+  return discrete_function_dpk_variant_list;
+}
+
+std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
+PolynomialReconstruction::build2(
+  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+  if (not hasSameMesh(discrete_function_variant_list)) {
+    throw NormalError("cannot reconstruct functions living of different mesh simultaneously");
+  }
+
+  auto mesh_v = getCommonMesh(discrete_function_variant_list);
+
+  return std::visit([&](auto&& p_mesh) { return this->_build(p_mesh, discrete_function_variant_list); },
+                    mesh_v->variant());
+}
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index 90b95c1db6eed6aeae60873969045a0f9846ebf8..0612e605b55f9a580e8e3bd13e6f46eba18c7dde 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -2,9 +2,11 @@
 #define POLYNOMIAL_RECONSTRUCTION_HPP
 
 #include <mesh/MeshTraits.hpp>
-#include <mesh/StencilBuilder.hpp>
+#include <mesh/StencilManager.hpp>
 #include <scheme/DiscreteFunctionDPk.hpp>
 
+#include <algebra/ShrinkMatrixView.hpp>
+#include <algebra/ShrinkVectorView.hpp>
 #include <algebra/SmallMatrix.hpp>
 #include <algebra/SmallVector.hpp>
 
@@ -16,31 +18,161 @@
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
 
+class DiscreteFunctionDPkVariant;
+class DiscreteFunctionVariant;
+
+#include <memory>
+#include <variant>
+
+class DiscreteFunctionDPkVariant
+{
+ public:
+  using Variant = std::variant<DiscreteFunctionDPk<1, const double>,
+                               DiscreteFunctionDPk<1, const TinyVector<1>>,
+                               DiscreteFunctionDPk<1, const TinyVector<2>>,
+                               DiscreteFunctionDPk<1, const TinyVector<3>>,
+                               DiscreteFunctionDPk<1, const TinyMatrix<1>>,
+                               DiscreteFunctionDPk<1, const TinyMatrix<2>>,
+                               DiscreteFunctionDPk<1, const TinyMatrix<3>>,
+
+                               DiscreteFunctionDPk<2, const double>,
+                               DiscreteFunctionDPk<2, const TinyVector<1>>,
+                               DiscreteFunctionDPk<2, const TinyVector<2>>,
+                               DiscreteFunctionDPk<2, const TinyVector<3>>,
+                               DiscreteFunctionDPk<2, const TinyMatrix<1>>,
+                               DiscreteFunctionDPk<2, const TinyMatrix<2>>,
+                               DiscreteFunctionDPk<2, const TinyMatrix<3>>,
+
+                               DiscreteFunctionDPk<3, const double>,
+                               DiscreteFunctionDPk<3, const TinyVector<1>>,
+                               DiscreteFunctionDPk<3, const TinyVector<2>>,
+                               DiscreteFunctionDPk<3, const TinyVector<3>>,
+                               DiscreteFunctionDPk<3, const TinyMatrix<1>>,
+                               DiscreteFunctionDPk<3, const TinyMatrix<2>>,
+                               DiscreteFunctionDPk<3, const TinyMatrix<3>>>;
+
+#warning add DiscreteFunctionDPkVector to the variant
+
+ private:
+  Variant m_discrete_function_dpk;
+
+ public:
+  PUGS_INLINE
+  const Variant&
+  discreteFunctionDPk() const
+  {
+    return m_discrete_function_dpk;
+  }
+
+  template <typename DiscreteFunctionDPkT>
+  PUGS_INLINE auto
+  get() const
+  {
+    static_assert(is_discrete_function_dpk_v<DiscreteFunctionDPkT>, "invalid template argument");
+    using DataType = typename DiscreteFunctionDPkT::data_type;
+    static_assert(std::is_const_v<DataType>, "data type of extracted discrete function must be const");
+
+#ifndef NDEBUG
+    if (not std::holds_alternative<DiscreteFunctionDPkT>(this->m_discrete_function_dpk)) {
+      std::ostringstream error_msg;
+      error_msg << "invalid discrete function type\n";
+      error_msg << "- required " << rang::fgB::red << demangle<DiscreteFunctionDPkT>() << rang::fg::reset << '\n';
+      error_msg << "- contains " << rang::fgB::yellow
+                << std::visit([](auto&& f) -> std::string { return demangle<decltype(f)>(); },
+                              this->m_discrete_function_dpk)
+                << rang::fg::reset;
+      throw NormalError(error_msg.str());
+    }
+#endif   // NDEBUG
+
+    return std::get<DiscreteFunctionDPkT>(this->discreteFunctionDPk());
+  }
+
+  template <size_t Dimension, typename DataType>
+  DiscreteFunctionDPkVariant(const DiscreteFunctionDPk<Dimension, DataType>& discrete_function_dpk)
+    : m_discrete_function_dpk{DiscreteFunctionDPk<Dimension, const DataType>{discrete_function_dpk}}
+  {
+    static_assert(std::is_same_v<std::remove_const_t<DataType>, double> or                       //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<1, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<2, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyVector<3, double>> or      //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<1, 1, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<2, 2, double>> or   //
+                    std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<3, 3, double>>,
+                  "DiscreteFunctionDPk with this DataType is not allowed in variant");
+  }
+
+  DiscreteFunctionDPkVariant& operator=(DiscreteFunctionDPkVariant&&)      = default;
+  DiscreteFunctionDPkVariant& operator=(const DiscreteFunctionDPkVariant&) = default;
+
+  DiscreteFunctionDPkVariant(const DiscreteFunctionDPkVariant&) = default;
+  DiscreteFunctionDPkVariant(DiscreteFunctionDPkVariant&&)      = default;
+
+  DiscreteFunctionDPkVariant()  = delete;
+  ~DiscreteFunctionDPkVariant() = default;
+};
+
 class PolynomialReconstruction
 {
+ private:
+  template <MeshConcept MeshType>
+  [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build(
+    const std::shared_ptr<const MeshType>& p_mesh,
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+
+  template <MeshConcept MeshType>
+  [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build2(
+    const std::shared_ptr<const MeshType>& p_mesh,
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+
  public:
   template <MeshConcept MeshType, typename DataType>
   [[nodiscard]] PUGS_INLINE DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>>
-  build(const MeshType& mesh, const DiscreteFunctionP0<DataType> p0_function)
+  build(const MeshType& mesh, const DiscreteFunctionP0<DataType> p0_function) const
   {
     using Rd = TinyVector<MeshType::Dimension>;
 
-    const size_t degree   = 1;
-    const Stencil stencil = StencilBuilder{}.build(mesh.connectivity());
+    const size_t degree    = 1;
+    const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity());
 
     auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
     auto cell_is_owned = mesh.connectivity().cellIsOwned();
 
+    const size_t max_stencil_size = [&]() {
+      size_t max_size = 0;
+      for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+        auto stencil_cell_list = stencil[cell_id];
+        if (cell_is_owned[cell_id] and stencil_cell_list.size() > max_size) {
+          max_size = stencil_cell_list.size();
+        }
+      }
+      return max_size;
+    }();
+
+    Kokkos::Experimental::UniqueToken<Kokkos::DefaultExecutionSpace::execution_space,
+                                      Kokkos::Experimental::UniqueTokenScope::Global>
+      tokens;
+
     DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>> dpk{mesh.meshVariant(), degree};
 
     if constexpr (is_polygonal_mesh_v<MeshType>) {
       if constexpr (std::is_arithmetic_v<DataType>) {
+        SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
+        SmallArray<SmallVector<double>> b_pool(Kokkos::DefaultExecutionSpace::concurrency());
+        for (size_t i = 0; i < A_pool.size(); ++i) {
+          A_pool[i] = SmallMatrix<double>(max_stencil_size, MeshType::Dimension);
+          b_pool[i] = SmallVector<double>(max_stencil_size);
+        }
+
         parallel_for(
           mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
             if (cell_is_owned[cell_j_id]) {
+              const int32_t t = tokens.acquire();
+
               auto stencil_cell_list = stencil[cell_j_id];
-              SmallVector<double> b(stencil_cell_list.size());
+
+              ShrinkVectorView b(b_pool[t], stencil_cell_list.size());
 
               const double pj = p0_function[cell_j_id];
               for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
@@ -49,7 +181,7 @@ class PolynomialReconstruction
                 b[i] = p0_function[cell_i_id] - pj;
               }
 
-              SmallMatrix<double> A(stencil_cell_list.size(), MeshType::Dimension);
+              ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
 
               const Rd& Xj = xj[cell_j_id];
 
@@ -61,7 +193,7 @@ class PolynomialReconstruction
                 }
               }
 
-              SmallVector<double> x(MeshType::Dimension);
+              TinyVector<MeshType::Dimension> x;
               Givens::solveInPlace(A, x, b);
 
               auto dpk_j = dpk.coefficients(cell_j_id);
@@ -71,14 +203,25 @@ class PolynomialReconstruction
               for (size_t l = 0; l < MeshType::Dimension; ++l) {
                 dpk_j[1 + l] = x[l];
               }
+
+              tokens.release(t);
             }
           });
       } else if constexpr (is_tiny_vector_v<DataType>) {
+        SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
+        SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency());
+        for (size_t i = 0; i < A_pool.size(); ++i) {
+          A_pool[i] = SmallMatrix<double>(max_stencil_size, MeshType::Dimension);
+          B_pool[i] = SmallMatrix<double>(max_stencil_size, DataType::Dimension);
+        }
+
         parallel_for(
           mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
             if (cell_is_owned[cell_j_id]) {
+              const int32_t t = tokens.acquire();
+
               auto stencil_cell_list = stencil[cell_j_id];
-              SmallMatrix<double> B(stencil_cell_list.size(), DataType::Dimension);
+              ShrinkMatrixView B(B_pool[t], stencil_cell_list.size());
 
               const DataType& pj = p0_function[cell_j_id];
               for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
@@ -89,7 +232,7 @@ class PolynomialReconstruction
                 }
               }
 
-              SmallMatrix<double> A(stencil_cell_list.size(), MeshType::Dimension);
+              ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
 
               const Rd& Xj = xj[cell_j_id];
               for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
@@ -100,7 +243,7 @@ class PolynomialReconstruction
                 }
               }
 
-              SmallMatrix<double> X(MeshType::Dimension, DataType::Dimension);
+              TinyMatrix<MeshType::Dimension, DataType::Dimension> X;
               Givens::solveCollectionInPlace(A, X, B);
 
               auto dpk_j = dpk.coefficients(cell_j_id);
@@ -112,14 +255,25 @@ class PolynomialReconstruction
                   dpk_j_ip1[k] = X(i, k);
                 }
               }
+
+              tokens.release(t);
             }
           });
       } else if constexpr (is_tiny_matrix_v<DataType>) {
+        SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
+        SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency());
+        for (size_t i = 0; i < A_pool.size(); ++i) {
+          A_pool[i] = SmallMatrix<double>(max_stencil_size, MeshType::Dimension);
+          B_pool[i] = SmallMatrix<double>(max_stencil_size, DataType::Dimension);
+        }
+
         parallel_for(
           mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
             if (cell_is_owned[cell_j_id]) {
+              const int32_t t = tokens.acquire();
+
               auto stencil_cell_list = stencil[cell_j_id];
-              SmallMatrix<double> B(stencil_cell_list.size(), DataType::Dimension);
+              ShrinkMatrixView B(B_pool[t], stencil_cell_list.size());
 
               const DataType& pj = p0_function[cell_j_id];
               for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
@@ -132,7 +286,7 @@ class PolynomialReconstruction
                 }
               }
 
-              SmallMatrix<double> A(stencil_cell_list.size(), MeshType::Dimension);
+              ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
 
               const Rd& Xj = xj[cell_j_id];
 
@@ -144,7 +298,7 @@ class PolynomialReconstruction
                 }
               }
 
-              SmallMatrix<double> X(MeshType::Dimension, DataType::Dimension);
+              TinyMatrix<MeshType::Dimension, DataType::Dimension> X;
               Givens::solveCollectionInPlace(A, X, B);
 
               auto dpk_j = dpk.coefficients(cell_j_id);
@@ -158,6 +312,8 @@ class PolynomialReconstruction
                   }
                 }
               }
+
+              tokens.release(t);
             }
           });
       } else {
@@ -169,6 +325,12 @@ class PolynomialReconstruction
     return dpk;
   }
 
+  [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> build(
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+
+  [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> build2(
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+
   PolynomialReconstruction() {}
 };
 
diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9ec19396e28f00119412582c9e42cd2d4ed25bcb
--- /dev/null
+++ b/src/scheme/test_reconstruction.cpp
@@ -0,0 +1,64 @@
+#include <scheme/test_reconstruction.hpp>
+
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/Timer.hpp>
+
+void
+test_reconstruction(const std::shared_ptr<const MeshVariant>& mesh_v,
+                    const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function_v)
+{
+  std::cout << "** one variable at a time (30 times)\n";
+
+  std::visit(
+    [&](auto&& p_mesh) {
+      std::visit(
+        [&](auto&& discrete_function) {
+          using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+          PolynomialReconstruction reconstruction;
+          if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+            Timer t;
+            for (size_t i = 0; i < 30; ++i) {
+              auto res = reconstruction.build(*p_mesh, discrete_function);
+            }
+            t.pause();
+            std::cout << "t = " << t << '\n';
+          }
+        },
+        discrete_function_v->discreteFunction());
+    },
+    mesh_v->variant());
+
+  std::cout << "finished!\n";
+}
+
+void
+test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list)
+{
+  std::cout << "** variable list at once (30 times)\n";
+
+  {
+    PolynomialReconstruction reconstruction;
+    Timer t;
+    for (size_t i = 0; i < 30; ++i) {
+      auto res = reconstruction.build(discrete_function_variant_list);
+    }
+    t.pause();
+    std::cout << "t = " << t << '\n';
+  }
+
+  std::cout << "finished!\n";
+
+  std::cout << "** variable list at once (30 times)\n";
+
+  {
+    PolynomialReconstruction reconstruction;
+    Timer t;
+    for (size_t i = 0; i < 30; ++i) {
+      auto res = reconstruction.build2(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
new file mode 100644
index 0000000000000000000000000000000000000000..f74bf078161e73b6d7267eda0e4e540f7f5e18ec
--- /dev/null
+++ b/src/scheme/test_reconstruction.hpp
@@ -0,0 +1,16 @@
+#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::shared_ptr<const MeshVariant>& mesh_v,
+                         const std::shared_ptr<const DiscreteFunctionVariant>& discrete_function_variant_list);
+
+void test_reconstruction(
+  const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list);
+
+#endif   // TEST_RECONSTRUCTION_HPP
diff --git a/src/utils/PugsTraits.hpp b/src/utils/PugsTraits.hpp
index ab2244b69095f5bba89dc3d823590524fc7a49e5..cb8f154e3b344ae34bb2b769b5a0419e37d44e4a 100644
--- a/src/utils/PugsTraits.hpp
+++ b/src/utils/PugsTraits.hpp
@@ -26,6 +26,12 @@ class DiscreteFunctionP0;
 template <typename DataType>
 class DiscreteFunctionP0Vector;
 
+template <size_t Dimension, typename DataType, typename BasisView>
+class DiscreteFunctionDPk;
+
+template <size_t Dimension, typename DataType, typename BasisView>
+class DiscreteFunctionDPkVector;
+
 // Traits is_trivially_castable
 
 template <typename T>
@@ -163,6 +169,29 @@ constexpr inline bool is_discrete_function_P0_vector_v<DiscreteFunctionP0Vector<
 template <typename T>
 constexpr inline bool is_discrete_function_v = is_discrete_function_P0_v<T> or is_discrete_function_P0_vector_v<T>;
 
+// Trais is DiscreteFunctionDPk
+
+template <typename T>
+constexpr inline bool is_discrete_function_dpk_scalar_v = false;
+
+template <size_t Dimension, typename DataType, typename BasisView>
+constexpr inline bool is_discrete_function_dpk_scalar_v<DiscreteFunctionDPk<Dimension, DataType, BasisView>> = true;
+
+// Trais is DiscreteFunctionDPkVector
+
+template <typename T>
+constexpr inline bool is_discrete_function_dpk_vector_v = false;
+
+template <size_t Dimension, typename DataType, typename BasisView>
+constexpr inline bool is_discrete_function_dpk_vector_v<DiscreteFunctionDPkVector<Dimension, DataType, BasisView>> =
+  true;
+
+// Trais is DiscreteFunction
+
+template <typename T>
+constexpr inline bool is_discrete_function_dpk_v =
+  is_discrete_function_dpk_scalar_v<T> or is_discrete_function_dpk_vector_v<T>;
+
 // helper to check if a type is part of a variant
 
 template <typename T, typename V>
diff --git a/tests/mpi_test_main.cpp b/tests/mpi_test_main.cpp
index ef3fdfcc00569ea3042de6724171fbd24e2ed680..657cccc8ce4a6129de59f1df111414fdd5236799 100644
--- a/tests/mpi_test_main.cpp
+++ b/tests/mpi_test_main.cpp
@@ -7,6 +7,7 @@
 #include <mesh/DualConnectivityManager.hpp>
 #include <mesh/DualMeshManager.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <mesh/StencilManager.hpp>
 #include <mesh/SynchronizerManager.hpp>
 #include <utils/GlobalVariableManager.hpp>
 #include <utils/Messenger.hpp>
@@ -99,6 +100,8 @@ main(int argc, char* argv[])
       MeshDataManager::create();
       DualConnectivityManager::create();
       DualMeshManager::create();
+      StencilManager::create();
+
       GlobalVariableManager::create();
 
       MeshDataBaseForTests::create();
@@ -111,6 +114,7 @@ main(int argc, char* argv[])
 
       MeshDataBaseForTests::destroy();
 
+      StencilManager::destroy();
       GlobalVariableManager::destroy();
       DualMeshManager::destroy();
       DualConnectivityManager::destroy();
diff --git a/tests/test_StencilBuilder.cpp b/tests/test_StencilBuilder.cpp
index af6031ee499253c283683e7e15d7e779327def17..0cb54f2623541fc97c1ce5058313c4c2f3b2ab5d 100644
--- a/tests/test_StencilBuilder.cpp
+++ b/tests/test_StencilBuilder.cpp
@@ -8,7 +8,7 @@
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshVariant.hpp>
-#include <mesh/StencilBuilder.hpp>
+#include <mesh/StencilManager.hpp>
 #include <utils/Messenger.hpp>
 
 // clazy:excludeall=non-pod-global-static
@@ -59,7 +59,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
 
       const Connectivity<1>& connectivity = mesh.connectivity();
 
-      Stencil stencil = StencilBuilder{}.build(connectivity);
+      Stencil stencil = StencilManager::instance().getStencil(connectivity);
 
       REQUIRE(is_valid(connectivity, stencil));
     }
@@ -69,7 +69,9 @@ TEST_CASE("StencilBuilder", "[mesh]")
       const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
 
       const Connectivity<1>& connectivity = mesh.connectivity();
-      REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
+      Stencil stencil                     = StencilManager::instance().getStencil(connectivity);
+
+      REQUIRE(is_valid(connectivity, stencil));
     }
   }
 
@@ -80,7 +82,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
       const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
 
       const Connectivity<2>& connectivity = mesh.connectivity();
-      REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
+      REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
     }
 
     SECTION("hybrid")
@@ -88,7 +90,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
       const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
 
       const Connectivity<2>& connectivity = mesh.connectivity();
-      REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
+      REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
     }
   }
 
@@ -99,7 +101,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
       const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
 
       const Connectivity<3>& connectivity = mesh.connectivity();
-      REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
+      REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
     }
 
     SECTION("hybrid")
@@ -107,7 +109,7 @@ TEST_CASE("StencilBuilder", "[mesh]")
       const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
 
       const Connectivity<3>& connectivity = mesh.connectivity();
-      REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
+      REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
     }
   }
 }
diff --git a/tests/test_main.cpp b/tests/test_main.cpp
index d8641d318f243eb624915882a92fc15d57a71346..3f81068df53dc2708110b1dd719cad2009cdbfb3 100644
--- a/tests/test_main.cpp
+++ b/tests/test_main.cpp
@@ -7,6 +7,7 @@
 #include <mesh/DualConnectivityManager.hpp>
 #include <mesh/DualMeshManager.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <mesh/StencilManager.hpp>
 #include <mesh/SynchronizerManager.hpp>
 #include <utils/GlobalVariableManager.hpp>
 #include <utils/Messenger.hpp>
@@ -58,6 +59,8 @@ main(int argc, char* argv[])
       MeshDataManager::create();
       DualConnectivityManager::create();
       DualMeshManager::create();
+      StencilManager::create();
+
       GlobalVariableManager::create();
 
       MeshDataBaseForTests::create();
@@ -71,6 +74,7 @@ main(int argc, char* argv[])
       MeshDataBaseForTests::destroy();
 
       GlobalVariableManager::destroy();
+      StencilManager::destroy();
       DualMeshManager::destroy();
       DualConnectivityManager::destroy();
       MeshDataManager::destroy();