diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 2e6e1fabbf51ddb49d54ea8fdc1d9b072ff1075c..0f3982c4f3b6816996906550821308757e251a18 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -4,17 +4,6 @@
 #include <algebra/ShrinkMatrixView.hpp>
 #include <algebra/ShrinkVectorView.hpp>
 #include <algebra/SmallMatrix.hpp>
-#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
-#include <analysis/GaussQuadratureDescriptor.hpp>
-#include <analysis/QuadratureFormula.hpp>
-#include <analysis/QuadratureManager.hpp>
-#include <geometry/CubeTransformation.hpp>
-#include <geometry/LineTransformation.hpp>
-#include <geometry/PrismTransformation.hpp>
-#include <geometry/PyramidTransformation.hpp>
-#include <geometry/SquareTransformation.hpp>
-#include <geometry/TetrahedronTransformation.hpp>
-#include <geometry/TriangleTransformation.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshFlatFaceBoundary.hpp>
@@ -24,11 +13,12 @@
 #include <scheme/DiscreteFunctionDPkVariant.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
-
 #include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
 #include <scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp>
 #include <scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp>
+#include <scheme/reconstruction_utils/MutableDiscreteFunctionDPkVariant.hpp>
 
+#warning put in a file in geometry
 template <size_t Dimension>
 PUGS_INLINE auto
 symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimension>& u)
@@ -53,127 +43,6 @@ symmetrize_coordinates(const TinyVector<Dimension>& origin,
   return u - 2 * dot(u - origin, normal) * normal;
 }
 
-class PolynomialReconstruction::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>>,
-
-                               DiscreteFunctionDPkVector<1, double>,
-                               DiscreteFunctionDPkVector<1, TinyVector<1>>,
-                               DiscreteFunctionDPkVector<1, TinyVector<2>>,
-                               DiscreteFunctionDPkVector<1, TinyVector<3>>,
-                               DiscreteFunctionDPkVector<1, TinyMatrix<1>>,
-                               DiscreteFunctionDPkVector<1, TinyMatrix<2>>,
-                               DiscreteFunctionDPkVector<1, TinyMatrix<3>>,
-
-                               DiscreteFunctionDPkVector<2, double>,
-                               DiscreteFunctionDPkVector<2, TinyVector<1>>,
-                               DiscreteFunctionDPkVector<2, TinyVector<2>>,
-                               DiscreteFunctionDPkVector<2, TinyVector<3>>,
-                               DiscreteFunctionDPkVector<2, TinyMatrix<1>>,
-                               DiscreteFunctionDPkVector<2, TinyMatrix<2>>,
-                               DiscreteFunctionDPkVector<2, TinyMatrix<3>>,
-
-                               DiscreteFunctionDPkVector<3, double>,
-                               DiscreteFunctionDPkVector<3, TinyVector<1>>,
-                               DiscreteFunctionDPkVector<3, TinyVector<2>>,
-                               DiscreteFunctionDPkVector<3, TinyVector<3>>,
-                               DiscreteFunctionDPkVector<3, TinyMatrix<1>>,
-                               DiscreteFunctionDPkVector<3, TinyMatrix<2>>,
-                               DiscreteFunctionDPkVector<3, TinyMatrix<3>>>;
-
- 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{discrete_function_dpk}
-  {
-    static_assert(std::is_same_v<DataType, double> or                       //
-                    std::is_same_v<DataType, TinyVector<1, double>> or      //
-                    std::is_same_v<DataType, TinyVector<2, double>> or      //
-                    std::is_same_v<DataType, TinyVector<3, double>> or      //
-                    std::is_same_v<DataType, TinyMatrix<1, 1, double>> or   //
-                    std::is_same_v<DataType, TinyMatrix<2, 2, double>> or   //
-                    std::is_same_v<DataType, TinyMatrix<3, 3, double>>,
-                  "DiscreteFunctionDPk with this DataType is not allowed in variant");
-  }
-
-  template <size_t Dimension, typename DataType>
-  MutableDiscreteFunctionDPkVariant(const DiscreteFunctionDPkVector<Dimension, DataType>& discrete_function_dpk)
-    : m_mutable_discrete_function_dpk{discrete_function_dpk}
-  {
-    static_assert(std::is_same_v<DataType, double> or                       //
-                    std::is_same_v<DataType, TinyVector<1, double>> or      //
-                    std::is_same_v<DataType, TinyVector<2, double>> or      //
-                    std::is_same_v<DataType, TinyVector<3, double>> or      //
-                    std::is_same_v<DataType, TinyMatrix<1, 1, double>> or   //
-                    std::is_same_v<DataType, TinyMatrix<2, 2, double>> or   //
-                    std::is_same_v<DataType, TinyMatrix<3, 3, double>>,
-                  "DiscreteFunctionDPkVector 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;
-};
-
 size_t
 PolynomialReconstruction::_getNumberOfColumns(
   const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
diff --git a/src/scheme/reconstruction_utils/MutableDiscreteFunctionDPkVariant.hpp b/src/scheme/reconstruction_utils/MutableDiscreteFunctionDPkVariant.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..761f3c1129b00609b1be9a23b5d256db4b48c084
--- /dev/null
+++ b/src/scheme/reconstruction_utils/MutableDiscreteFunctionDPkVariant.hpp
@@ -0,0 +1,130 @@
+#ifndef MUTABLE_DISCRETE_FUNCTION_D_PK_VARIANT_HPP
+#define MUTABLE_DISCRETE_FUNCTION_D_PK_VARIANT_HPP
+
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+
+#include <variant>
+
+class PolynomialReconstruction::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>>,
+
+                               DiscreteFunctionDPkVector<1, double>,
+                               DiscreteFunctionDPkVector<1, TinyVector<1>>,
+                               DiscreteFunctionDPkVector<1, TinyVector<2>>,
+                               DiscreteFunctionDPkVector<1, TinyVector<3>>,
+                               DiscreteFunctionDPkVector<1, TinyMatrix<1>>,
+                               DiscreteFunctionDPkVector<1, TinyMatrix<2>>,
+                               DiscreteFunctionDPkVector<1, TinyMatrix<3>>,
+
+                               DiscreteFunctionDPkVector<2, double>,
+                               DiscreteFunctionDPkVector<2, TinyVector<1>>,
+                               DiscreteFunctionDPkVector<2, TinyVector<2>>,
+                               DiscreteFunctionDPkVector<2, TinyVector<3>>,
+                               DiscreteFunctionDPkVector<2, TinyMatrix<1>>,
+                               DiscreteFunctionDPkVector<2, TinyMatrix<2>>,
+                               DiscreteFunctionDPkVector<2, TinyMatrix<3>>,
+
+                               DiscreteFunctionDPkVector<3, double>,
+                               DiscreteFunctionDPkVector<3, TinyVector<1>>,
+                               DiscreteFunctionDPkVector<3, TinyVector<2>>,
+                               DiscreteFunctionDPkVector<3, TinyVector<3>>,
+                               DiscreteFunctionDPkVector<3, TinyMatrix<1>>,
+                               DiscreteFunctionDPkVector<3, TinyMatrix<2>>,
+                               DiscreteFunctionDPkVector<3, TinyMatrix<3>>>;
+
+ 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{discrete_function_dpk}
+  {
+    static_assert(std::is_same_v<DataType, double> or                       //
+                    std::is_same_v<DataType, TinyVector<1, double>> or      //
+                    std::is_same_v<DataType, TinyVector<2, double>> or      //
+                    std::is_same_v<DataType, TinyVector<3, double>> or      //
+                    std::is_same_v<DataType, TinyMatrix<1, 1, double>> or   //
+                    std::is_same_v<DataType, TinyMatrix<2, 2, double>> or   //
+                    std::is_same_v<DataType, TinyMatrix<3, 3, double>>,
+                  "DiscreteFunctionDPk with this DataType is not allowed in variant");
+  }
+
+  template <size_t Dimension, typename DataType>
+  MutableDiscreteFunctionDPkVariant(const DiscreteFunctionDPkVector<Dimension, DataType>& discrete_function_dpk)
+    : m_mutable_discrete_function_dpk{discrete_function_dpk}
+  {
+    static_assert(std::is_same_v<DataType, double> or                       //
+                    std::is_same_v<DataType, TinyVector<1, double>> or      //
+                    std::is_same_v<DataType, TinyVector<2, double>> or      //
+                    std::is_same_v<DataType, TinyVector<3, double>> or      //
+                    std::is_same_v<DataType, TinyMatrix<1, 1, double>> or   //
+                    std::is_same_v<DataType, TinyMatrix<2, 2, double>> or   //
+                    std::is_same_v<DataType, TinyMatrix<3, 3, double>>,
+                  "DiscreteFunctionDPkVector 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;
+};
+
+#endif   // MUTABLE_DISCRETE_FUNCTION_D_PK_VARIANT_HPP