diff --git a/src/scheme/DiscreteFunctionDPkVariant.hpp b/src/scheme/DiscreteFunctionDPkVariant.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0d683c1783c9aad299a81dcc4d7ef3817470d54a
--- /dev/null
+++ b/src/scheme/DiscreteFunctionDPkVariant.hpp
@@ -0,0 +1,108 @@
+#ifndef DISCRETE_FUNCTION_DPK_VARIANT_HPP
+#define DISCRETE_FUNCTION_DPK_VARIANT_HPP
+
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVector.hpp>
+
+#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>>,
+
+                               DiscreteFunctionDPkVector<1, const double>,
+                               DiscreteFunctionDPkVector<2, const double>,
+                               DiscreteFunctionDPkVector<3, const double>>;
+
+ 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");
+  }
+
+  template <size_t Dimension, typename DataType>
+  DiscreteFunctionDPkVariant(const DiscreteFunctionDPkVector<Dimension, DataType>& discrete_function_dpk)
+    : m_discrete_function_dpk{DiscreteFunctionDPkVector<Dimension, const DataType>{discrete_function_dpk}}
+  {
+    static_assert(std::is_same_v<std::remove_const_t<DataType>, double>,
+                  "DiscreteFunctionDPkVector 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;
+};
+
+#endif   // DISCRETE_FUNCTION_DPK_VARIANT_HPP
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 2e400c4e035282b3297bf5ee896c6c1e87c50397..0e10c1bf8c32f402e19eed70b48f434234d42632 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -3,7 +3,7 @@
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
 
-class MutableDiscreteFunctionDPkVariant
+class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
 {
  public:
   using Variant = std::variant<DiscreteFunctionDPk<1, double>,
@@ -102,9 +102,14 @@ class MutableDiscreteFunctionDPkVariant
 template <MeshConcept MeshType>
 [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
 PolynomialReconstruction::_build(
+  const size_t degree,
   const std::shared_ptr<const MeshType>& p_mesh,
   const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
 {
+  if (degree != 1) {
+    throw NotImplementedError("only degree 1 is available");
+  }
+
   const MeshType& mesh = *p_mesh;
 
   using Rd = TinyVector<MeshType::Dimension>;
@@ -139,7 +144,6 @@ PolynomialReconstruction::_build(
     return n;
   }();
 
-  const size_t degree = 1;
   const size_t basis_dimension =
     DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(degree);
 
@@ -392,6 +396,7 @@ PolynomialReconstruction::_build(
 
 std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
 PolynomialReconstruction::build(
+  const size_t degree,
   const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
 {
   if (not hasSameMesh(discrete_function_variant_list)) {
@@ -400,6 +405,6 @@ PolynomialReconstruction::build(
 
   auto mesh_v = getCommonMesh(discrete_function_variant_list);
 
-  return std::visit([&](auto&& p_mesh) { return this->_build(p_mesh, discrete_function_variant_list); },
+  return std::visit([&](auto&& p_mesh) { return this->_build(degree, p_mesh, discrete_function_variant_list); },
                     mesh_v->variant());
 }
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index 430ce0234368222420a3b4eb698a42142f1c98a9..454e77c9d748399adcdb2fe2bbfef3e7174b4a3f 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -1,16 +1,14 @@
 #ifndef POLYNOMIAL_RECONSTRUCTION_HPP
 #define POLYNOMIAL_RECONSTRUCTION_HPP
 
-#include <mesh/MeshTraits.hpp>
-#include <mesh/StencilManager.hpp>
-#include <scheme/DiscreteFunctionDPk.hpp>
-#include <scheme/DiscreteFunctionDPkVector.hpp>
-
 #include <algebra/ShrinkMatrixView.hpp>
 #include <algebra/ShrinkVectorView.hpp>
 #include <algebra/SmallMatrix.hpp>
 #include <algebra/SmallVector.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
 
+#warning MOVE TO .cpp WHEN FINISHED
 #include <algebra/Givens.hpp>
 #include <analysis/GaussLegendreQuadratureDescriptor.hpp>
 #include <analysis/QuadratureFormula.hpp>
@@ -18,116 +16,19 @@
 #include <geometry/LineTransformation.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <mesh/StencilManager.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>>,
-
-                               DiscreteFunctionDPkVector<1, const double>,
-                               DiscreteFunctionDPkVector<2, const double>,
-                               DiscreteFunctionDPkVector<3, const double>>;
-
- 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");
-  }
-
-  template <size_t Dimension, typename DataType>
-  DiscreteFunctionDPkVariant(const DiscreteFunctionDPkVector<Dimension, DataType>& discrete_function_dpk)
-    : m_discrete_function_dpk{DiscreteFunctionDPkVector<Dimension, const DataType>{discrete_function_dpk}}
-  {
-    static_assert(std::is_same_v<std::remove_const_t<DataType>, double>,
-                  "DiscreteFunctionDPkVector 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:
+  class MutableDiscreteFunctionDPkVariant;
+
   template <MeshConcept MeshType>
   [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build(
+    const size_t degree,
     const std::shared_ptr<const MeshType>& p_mesh,
     const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
 
@@ -332,11 +233,12 @@ class PolynomialReconstruction
   }
 
   [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> build(
+    const size_t degree,
     const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
 
   template <typename... DiscreteFunctionT>
   [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
-  build(DiscreteFunctionT... input) const
+  build(const size_t degree, DiscreteFunctionT... input) const
   {
     std::vector<std::shared_ptr<const DiscreteFunctionVariant>> variant_vector;
     auto convert = [&variant_vector](auto&& df) {
@@ -356,7 +258,7 @@ class PolynomialReconstruction
     };
 
     (convert(std::forward<DiscreteFunctionT>(input)), ...);
-    return this->build(variant_vector);
+    return this->build(degree, variant_vector);
   }
 
   PolynomialReconstruction() {}
diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
index 08527f96d0dd35838b3bbeebf5cb87db96c42b94..87b9bcc16193ed00a3d505f859c12c7ee14da214 100644
--- a/src/scheme/test_reconstruction.cpp
+++ b/src/scheme/test_reconstruction.cpp
@@ -55,9 +55,8 @@ test_reconstruction_one(const std::shared_ptr<const MeshVariant>& mesh_v,
 void
 test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list)
 {
-  std::cout << "** variable list one by one (30 times)\n";
-
   {
+    std::cout << "** variable list one by one (50 times)\n";
     Timer t;
     for (auto discrete_function_v : discrete_function_variant_list) {
       std::visit(
@@ -68,9 +67,11 @@ test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVari
               using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
               PolynomialReconstruction reconstruction;
               if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
-                for (size_t i = 0; i < 30; ++i) {
+                for (size_t i = 0; i < 50; ++i) {
                   auto res = reconstruction.build(*p_mesh, discrete_function);
                 }
+              } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+                std::cout << "Omitting discrete function p0 vector!\n";
               }
             },
             mesh_v->variant());
@@ -81,13 +82,12 @@ test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVari
     std::cout << "t = " << t << '\n';
   }
 
-  std::cout << "** variable list at once (30 times)\n";
-
   {
+    std::cout << "** variable list at once (50 times)\n";
     PolynomialReconstruction reconstruction;
     Timer t;
-    for (size_t i = 0; i < 30; ++i) {
-      auto res = reconstruction.build(discrete_function_variant_list);
+    for (size_t i = 0; i < 50; ++i) {
+      auto res = reconstruction.build(1, discrete_function_variant_list);
     }
     t.pause();
     std::cout << "t = " << t << '\n';
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index 029fa30e8aa88279e90896a5918dd024fb2e4474..7ce59a9109f1e13d040e46a69e3b2c2b76596790 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -21,930 +21,945 @@
 
 TEST_CASE("PolynomialReconstruction", "[scheme]")
 {
-  SECTION("1D")
+  SECTION("degree 1")
   {
-    using R1 = TinyVector<1>;
-
-    SECTION("R data")
+    SECTION("1D")
     {
-      for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
-          auto& mesh  = *p_mesh;
+      using R1 = TinyVector<1>;
+
+      SECTION("R data")
+      {
+        for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+            auto& mesh  = *p_mesh;
 
-          auto R_affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
-          auto xj       = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto R_affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
+            auto xj       = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0<double> fh{p_mesh};
+            DiscreteFunctionP0<double> fh{p_mesh};
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
 
-          auto reconstructions = PolynomialReconstruction{}.build(fh);
+            auto reconstructions = PolynomialReconstruction{}.build(1, fh);
 
-          auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
+            auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error = std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const double reconstructed_slope =
-                (dpk_fh[cell_id](R1{0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
+            {
+              double max_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const double reconstructed_slope =
+                  (dpk_fh[cell_id](R1{0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
 
-              max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
+                max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
+              }
+              REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
           }
         }
       }
-    }
 
-    SECTION("R^3 data")
-    {
-      using R3 = TinyVector<3>;
+      SECTION("R^3 data")
+      {
+        using R3 = TinyVector<3>;
 
-      for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
-          auto& mesh  = *p_mesh;
+        for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+            auto& mesh  = *p_mesh;
 
-          auto R3_affine = [](const R1& x) -> R3 {
-            return R3{+2.3 + 1.7 * x[0],   //
-                      +1.4 - 0.6 * x[0],   //
-                      -0.2 + 3.1 * x[0]};
-          };
-          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto R3_affine = [](const R1& x) -> R3 {
+              return R3{+2.3 + 1.7 * x[0],   //
+                        +1.4 - 0.6 * x[0],   //
+                        -0.2 + 3.1 * x[0]};
+            };
+            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0<R3> uh{p_mesh};
+            DiscreteFunctionP0<R3> uh{p_mesh};
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
 
-          auto reconstructions = PolynomialReconstruction{}.build(uh);
+            auto reconstructions = PolynomialReconstruction{}.build(1, uh);
 
-          auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3>>();
+            auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R3 reconstructed_slope =
-                (1 / 0.2) * (dpk_uh[cell_id](R1{0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R1{0.1}));
+            {
+              double max_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R3 reconstructed_slope =
+                  (1 / 0.2) * (dpk_uh[cell_id](R1{0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R1{0.1}));
 
-              max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+                max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+              }
+              REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
           }
         }
       }
-    }
 
-    SECTION("R^3x3 data")
-    {
-      using R3x3 = TinyMatrix<3, 3>;
-
-      for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
-          auto& mesh  = *p_mesh;
-
-          auto R3x3_affine = [](const R1& x) -> R3x3 {
-            return R3x3{
-              +2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0],   //
-              +2.4 - 2.3 * x[0], -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0],
-              -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
+      SECTION("R^3x3 data")
+      {
+        using R3x3 = TinyMatrix<3, 3>;
+
+        for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+            auto& mesh  = *p_mesh;
+
+            auto R3x3_affine = [](const R1& x) -> R3x3 {
+              return R3x3{
+                +2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0],   //
+                +2.4 - 2.3 * x[0], -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0],
+                -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
+              };
             };
-          };
-          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0<R3x3> Ah{p_mesh};
+            DiscreteFunctionP0<R3x3> Ah{p_mesh};
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); });
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); });
 
-          auto reconstructions = PolynomialReconstruction{}.build(Ah);
+            auto reconstructions = PolynomialReconstruction{}.build(1, Ah);
 
-          auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>();
+            auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error =
-                std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id])));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          {
-            double max_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R3x3 reconstructed_slope =
-                (1 / 0.2) * (dpk_Ah[cell_id](R1{0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R1{0.1}));
+            {
+              double max_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R3x3 reconstructed_slope =
+                  (1 / 0.2) * (dpk_Ah[cell_id](R1{0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R1{0.1}));
 
-              R3x3 slops = R3x3{+1.7, +2.1, -0.6,   //
-                                -2.3, +3.1, -3.6,   //
-                                +3.1, +2.9, +2.3};
+                R3x3 slops = R3x3{+1.7, +2.1, -0.6,   //
+                                  -2.3, +3.1, -3.6,   //
+                                  +3.1, +2.9, +2.3};
 
-              max_slope_error = std::max(max_slope_error,   //
-                                         frobeniusNorm(reconstructed_slope - slops));
+                max_slope_error = std::max(max_slope_error,   //
+                                           frobeniusNorm(reconstructed_slope - slops));
+              }
+              REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
           }
         }
       }
-    }
 
-    SECTION("R vector data")
-    {
-      using R3 = TinyVector<3>;
+      SECTION("R vector data")
+      {
+        using R3 = TinyVector<3>;
 
-      for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
-          auto& mesh  = *p_mesh;
+        for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+            auto& mesh  = *p_mesh;
 
-          auto vector_affine = [](const R1& x) -> R3 {
-            return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
-          };
-          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto vector_affine = [](const R1& x) -> R3 {
+              return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
+            };
+            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
+            DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-              auto vector = vector_affine(xj[cell_id]);
-              for (size_t i = 0; i < vector.dimension(); ++i) {
-                Vh[cell_id][i] = vector[i];
-              }
-            });
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                auto vector = vector_affine(xj[cell_id]);
+                for (size_t i = 0; i < vector.dimension(); ++i) {
+                  Vh[cell_id][i] = vector[i];
+                }
+              });
 
-          auto reconstructions = PolynomialReconstruction{}.build(Vh);
+            auto reconstructions = PolynomialReconstruction{}.build(1, Vh);
 
-          auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
+            auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              auto vector = vector_affine(xj[cell_id]);
-              for (size_t i = 0; i < vector.dimension(); ++i) {
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                auto vector = vector_affine(xj[cell_id]);
+                for (size_t i = 0; i < vector.dimension(); ++i) {
+                  max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+                }
               }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          {
-            double max_slope_error = 0;
-            const TinyVector<3> slope{+1.7, +2.1, -0.6};
+            {
+              double max_slope_error = 0;
+              const TinyVector<3> slope{+1.7, +2.1, -0.6};
 
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              for (size_t i = 0; i < slope.dimension(); ++i) {
-                const double reconstructed_slope =
-                  (1 / 0.2) * (dpk_Vh(cell_id, i)(R1{0.1} + xj[cell_id]) - dpk_Vh(cell_id, i)(xj[cell_id] - R1{0.1}));
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < slope.dimension(); ++i) {
+                  const double reconstructed_slope =
+                    (1 / 0.2) * (dpk_Vh(cell_id, i)(R1{0.1} + xj[cell_id]) - dpk_Vh(cell_id, i)(xj[cell_id] - R1{0.1}));
 
-                max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i]));
+                  max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i]));
+                }
+                REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
               }
-              REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
             }
           }
         }
       }
-    }
 
-    SECTION("list of various types")
-    {
-      using R3x3 = TinyMatrix<3>;
-      using R3   = TinyVector<3>;
-
-      for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
-          auto& mesh  = *p_mesh;
-
-          auto R_affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
-
-          auto R3_affine = [](const R1& x) -> R3 {
-            return R3{+2.3 + 1.7 * x[0],   //
-                      +1.4 - 0.6 * x[0],   //
-                      -0.2 + 3.1 * x[0]};
-          };
-
-          auto R3x3_affine = [](const R1& x) -> R3x3 {
-            return R3x3{
-              +2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0],   //
-              +2.4 - 2.3 * x[0], -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0],
-              -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
+      SECTION("list of various types")
+      {
+        using R3x3 = TinyMatrix<3>;
+        using R3   = TinyVector<3>;
+
+        for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+            auto& mesh  = *p_mesh;
+
+            auto R_affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
+
+            auto R3_affine = [](const R1& x) -> R3 {
+              return R3{+2.3 + 1.7 * x[0],   //
+                        +1.4 - 0.6 * x[0],   //
+                        -0.2 + 3.1 * x[0]};
             };
-          };
 
-          auto vector_affine = [](const R1& x) -> R3 {
-            return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
-          };
+            auto R3x3_affine = [](const R1& x) -> R3x3 {
+              return R3x3{
+                +2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0],   //
+                +2.4 - 2.3 * x[0], -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0],
+                -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
+              };
+            };
 
-          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto vector_affine = [](const R1& x) -> R3 {
+              return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
+            };
 
-          DiscreteFunctionP0<double> fh{p_mesh};
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
+            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0<R3> uh{p_mesh};
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+            DiscreteFunctionP0<double> fh{p_mesh};
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
 
-          DiscreteFunctionP0<R3x3> Ah{p_mesh};
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); });
+            DiscreteFunctionP0<R3> uh{p_mesh};
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
 
-          DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-              auto vector = vector_affine(xj[cell_id]);
-              for (size_t i = 0; i < vector.dimension(); ++i) {
-                Vh[cell_id][i] = vector[i];
-              }
-            });
+            DiscreteFunctionP0<R3x3> Ah{p_mesh};
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); });
 
-          auto reconstructions = PolynomialReconstruction{}.build(std::make_shared<DiscreteFunctionVariant>(fh), uh,
-                                                                  std::make_shared<DiscreteFunctionP0<R3x3>>(Ah),
-                                                                  DiscreteFunctionVariant(Vh));
+            DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                auto vector = vector_affine(xj[cell_id]);
+                for (size_t i = 0; i < vector.dimension(); ++i) {
+                  Vh[cell_id][i] = vector[i];
+                }
+              });
 
-          auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
+            auto reconstructions = PolynomialReconstruction{}.build(1, std::make_shared<DiscreteFunctionVariant>(fh),
+                                                                    uh, std::make_shared<DiscreteFunctionP0<R3x3>>(Ah),
+                                                                    DiscreteFunctionVariant(Vh));
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error = std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+            auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
+
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const double reconstructed_slope =
-                (dpk_fh[cell_id](R1{0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
+            {
+              double max_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const double reconstructed_slope =
+                  (dpk_fh[cell_id](R1{0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
 
-              max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
+                max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
+              }
+              REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>();
+            auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R3 reconstructed_slope =
-                (1 / 0.2) * (dpk_uh[cell_id](R1{0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R1{0.1}));
+            {
+              double max_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R3 reconstructed_slope =
+                  (1 / 0.2) * (dpk_uh[cell_id](R1{0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R1{0.1}));
 
-              max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+                max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+              }
+              REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          auto dpk_Ah = reconstructions[2]->get<DiscreteFunctionDPk<1, const R3x3>>();
+            auto dpk_Ah = reconstructions[2]->get<DiscreteFunctionDPk<1, const R3x3>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error =
-                std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id])));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          {
-            double max_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R3x3 reconstructed_slope =
-                (1 / 0.2) * (dpk_Ah[cell_id](R1{0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R1{0.1}));
+            {
+              double max_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R3x3 reconstructed_slope =
+                  (1 / 0.2) * (dpk_Ah[cell_id](R1{0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R1{0.1}));
 
-              R3x3 slops = R3x3{+1.7, +2.1, -0.6,   //
-                                -2.3, +3.1, -3.6,   //
-                                +3.1, +2.9, +2.3};
+                R3x3 slops = R3x3{+1.7, +2.1, -0.6,   //
+                                  -2.3, +3.1, -3.6,   //
+                                  +3.1, +2.9, +2.3};
 
-              max_slope_error = std::max(max_slope_error,   //
-                                         frobeniusNorm(reconstructed_slope - slops));
+                max_slope_error = std::max(max_slope_error,   //
+                                           frobeniusNorm(reconstructed_slope - slops));
+              }
+              REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>();
+            auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              auto vector = vector_affine(xj[cell_id]);
-              for (size_t i = 0; i < vector.dimension(); ++i) {
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                auto vector = vector_affine(xj[cell_id]);
+                for (size_t i = 0; i < vector.dimension(); ++i) {
+                  max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+                }
               }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          {
-            double max_slope_error = 0;
-            const TinyVector<3> slope{+1.7, +2.1, -0.6};
+            {
+              double max_slope_error = 0;
+              const TinyVector<3> slope{+1.7, +2.1, -0.6};
 
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              for (size_t i = 0; i < slope.dimension(); ++i) {
-                const double reconstructed_slope =
-                  (1 / 0.2) * (dpk_Vh(cell_id, i)(R1{0.1} + xj[cell_id]) - dpk_Vh(cell_id, i)(xj[cell_id] - R1{0.1}));
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < slope.dimension(); ++i) {
+                  const double reconstructed_slope =
+                    (1 / 0.2) * (dpk_Vh(cell_id, i)(R1{0.1} + xj[cell_id]) - dpk_Vh(cell_id, i)(xj[cell_id] - R1{0.1}));
 
-                max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i]));
+                  max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i]));
+                }
+                REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
               }
-              REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
             }
           }
         }
       }
     }
-  }
-
-  SECTION("2D")
-  {
-    using R2 = TinyVector<2>;
 
-    SECTION("R data")
+    SECTION("2D")
     {
-      for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
-          auto& mesh  = *p_mesh;
+      using R2 = TinyVector<2>;
 
-          auto R_affine = [](const R2& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1]; };
-          auto xj       = MeshDataManager::instance().getMeshData(mesh).xj();
+      SECTION("R data")
+      {
+        for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+            auto& mesh  = *p_mesh;
 
-          DiscreteFunctionP0<double> fh{p_mesh};
+            auto R_affine = [](const R2& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1]; };
+            auto xj       = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
+            DiscreteFunctionP0<double> fh{p_mesh};
 
-          auto reconstructions = PolynomialReconstruction{}.build(fh);
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
 
-          auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
+            auto reconstructions = PolynomialReconstruction{}.build(1, fh);
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error = std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+            auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
+
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_x_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const double reconstructed_slope =
-                (dpk_fh[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0.1, 0})) / 0.2;
+            {
+              double max_x_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const double reconstructed_slope =
+                  (dpk_fh[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0.1, 0})) / 0.2;
 
-              max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
+                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
+              }
+              REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_y_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const double reconstructed_slope =
-                (dpk_fh[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0, 0.1})) / 0.2;
+            {
+              double max_y_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const double reconstructed_slope =
+                  (dpk_fh[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0, 0.1})) / 0.2;
 
-              max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
+                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
+              }
+              REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
           }
         }
       }
-    }
 
-    SECTION("R^3 data")
-    {
-      using R3 = TinyVector<3>;
+      SECTION("R^3 data")
+      {
+        using R3 = TinyVector<3>;
 
-      for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
-          auto& mesh  = *p_mesh;
+        for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+            auto& mesh  = *p_mesh;
 
-          auto R3_affine = [](const R2& x) -> R3 {
-            return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1],   //
-                      +1.4 - 0.6 * x[0] + 1.3 * x[1],   //
-                      -0.2 + 3.1 * x[0] - 1.1 * x[1]};
-          };
-          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto R3_affine = [](const R2& x) -> R3 {
+              return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1],   //
+                        +1.4 - 0.6 * x[0] + 1.3 * x[1],   //
+                        -0.2 + 3.1 * x[0] - 1.1 * x[1]};
+            };
+            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0<R3> uh{p_mesh};
+            DiscreteFunctionP0<R3> uh{p_mesh};
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
 
-          auto reconstructions = PolynomialReconstruction{}.build(uh);
+            auto reconstructions = PolynomialReconstruction{}.build(1, uh);
 
-          auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R3>>();
+            auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R3>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_x_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R3 reconstructed_slope =
-                (1 / 0.2) * (dpk_uh[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R2{0.1, 0}));
+            {
+              double max_x_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R3 reconstructed_slope =
+                  (1 / 0.2) * (dpk_uh[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R2{0.1, 0}));
 
-              max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+                max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+              }
+              REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_y_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R3 reconstructed_slope =
-                (1 / 0.2) * (dpk_uh[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R2{0, 0.1}));
+            {
+              double max_y_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R3 reconstructed_slope =
+                  (1 / 0.2) * (dpk_uh[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R2{0, 0.1}));
 
-              max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
+                max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
+              }
+              REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
           }
         }
       }
-    }
 
-    SECTION("R^2x2 data")
-    {
-      using R2x2 = TinyMatrix<2, 2>;
+      SECTION("R^2x2 data")
+      {
+        using R2x2 = TinyMatrix<2, 2>;
 
-      for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
-          auto& mesh  = *p_mesh;
+        for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+            auto& mesh  = *p_mesh;
 
-          auto R2x2_affine = [](const R2& x) -> R2x2 {
-            return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1],   //
-                        +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
-          };
-          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto R2x2_affine = [](const R2& x) -> R2x2 {
+              return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1],   //
+                          +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
+            };
+            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0<R2x2> Ah{p_mesh};
+            DiscreteFunctionP0<R2x2> Ah{p_mesh};
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
 
-          auto reconstructions = PolynomialReconstruction{}.build(Ah);
+            auto reconstructions = PolynomialReconstruction{}.build(1, Ah);
 
-          auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>();
+            auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error =
-                std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          {
-            double max_x_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R2x2 reconstructed_slope =
-                (1 / 0.2) * (dpk_Ah[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R2{0.1, 0}));
+            {
+              double max_x_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R2x2 reconstructed_slope =
+                  (1 / 0.2) * (dpk_Ah[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R2{0.1, 0}));
 
-              max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, +2.1,   //
-                                                                                                       -0.6, -2.3}));
+                max_x_slope_error =
+                  std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, +2.1,   //
+                                                                                       -0.6, -2.3}));
+              }
+              REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
             }
-            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
-          }
 
-          {
-            double max_y_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R2x2 reconstructed_slope =
-                (1 / 0.2) * (dpk_Ah[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R2{0, 0.1}));
+            {
+              double max_y_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R2x2 reconstructed_slope =
+                  (1 / 0.2) * (dpk_Ah[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R2{0, 0.1}));
 
-              max_y_slope_error = std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2,   //
-                                                                                                       -2.1, +1.3}));
+                max_y_slope_error =
+                  std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2,   //
+                                                                                       -2.1, +1.3}));
+              }
+              REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
             }
-            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
           }
         }
       }
-    }
 
-    SECTION("vector data")
-    {
-      using R4 = TinyVector<4>;
+      SECTION("vector data")
+      {
+        using R4 = TinyVector<4>;
 
-      for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
-          auto& mesh  = *p_mesh;
+        for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+            auto& mesh  = *p_mesh;
 
-          auto vector_affine = [](const R2& x) -> R4 {
-            return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1],   //
-                      +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
-          };
-          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto vector_affine = [](const R2& x) -> R4 {
+              return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1],   //
+                        +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
+            };
+            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
+            DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-              auto vector = vector_affine(xj[cell_id]);
-              for (size_t i = 0; i < vector.dimension(); ++i) {
-                Vh[cell_id][i] = vector[i];
-              }
-            });
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                auto vector = vector_affine(xj[cell_id]);
+                for (size_t i = 0; i < vector.dimension(); ++i) {
+                  Vh[cell_id][i] = vector[i];
+                }
+              });
 
-          auto reconstructions = PolynomialReconstruction{}.build(Vh);
+            auto reconstructions = PolynomialReconstruction{}.build(1, Vh);
 
-          auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>();
+            auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              auto vector = vector_affine(xj[cell_id]);
-              for (size_t i = 0; i < vector.dimension(); ++i) {
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                auto vector = vector_affine(xj[cell_id]);
+                for (size_t i = 0; i < vector.dimension(); ++i) {
+                  max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+                }
               }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          {
-            double max_x_slope_error = 0;
-            const R4 slope{+1.7, +2.1, -0.6, -2.3};
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              for (size_t i = 0; i < slope.dimension(); ++i) {
-                const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0.1, 0} + xj[cell_id]) -
-                                                                dpk_Vh(cell_id, i)(xj[cell_id] - R2{0.1, 0}));
+            {
+              double max_x_slope_error = 0;
+              const R4 slope{+1.7, +2.1, -0.6, -2.3};
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < slope.dimension(); ++i) {
+                  const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0.1, 0} + xj[cell_id]) -
+                                                                  dpk_Vh(cell_id, i)(xj[cell_id] - R2{0.1, 0}));
 
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+                  max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+                }
               }
+              REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
             }
-            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
-          }
 
-          {
-            double max_y_slope_error = 0;
-            const R4 slope{+1.2, -2.2, -2.1, +1.3};
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              for (size_t i = 0; i < slope.dimension(); ++i) {
-                const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0, 0.1} + xj[cell_id]) -
-                                                                dpk_Vh(cell_id, i)(xj[cell_id] - R2{0, 0.1}));
+            {
+              double max_y_slope_error = 0;
+              const R4 slope{+1.2, -2.2, -2.1, +1.3};
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < slope.dimension(); ++i) {
+                  const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0, 0.1} + xj[cell_id]) -
+                                                                  dpk_Vh(cell_id, i)(xj[cell_id] - R2{0, 0.1}));
 
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+                  max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+                }
               }
+              REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
             }
-            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
           }
         }
       }
     }
-  }
-
-  SECTION("3D")
-  {
-    using R3 = TinyVector<3>;
 
-    SECTION("R data")
+    SECTION("3D")
     {
-      for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
-          auto& mesh  = *p_mesh;
+      using R3 = TinyVector<3>;
 
-          auto R_affine = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; };
-          auto xj       = MeshDataManager::instance().getMeshData(mesh).xj();
+      SECTION("R data")
+      {
+        for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+            auto& mesh  = *p_mesh;
 
-          DiscreteFunctionP0<double> fh{p_mesh};
+            auto R_affine = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; };
+            auto xj       = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
+            DiscreteFunctionP0<double> fh{p_mesh};
 
-          auto reconstructions = PolynomialReconstruction{}.build(fh);
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
 
-          auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
+            auto reconstructions = PolynomialReconstruction{}.build(1, fh);
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error = std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+            auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
+
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_x_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const double reconstructed_slope =
-                (dpk_fh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0.1, 0, 0})) / 0.2;
+            {
+              double max_x_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const double reconstructed_slope =
+                  (dpk_fh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0.1, 0, 0})) / 0.2;
 
-              max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
+                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
+              }
+              REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_y_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const double reconstructed_slope =
-                (dpk_fh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0.1, 0})) / 0.2;
+            {
+              double max_y_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const double reconstructed_slope =
+                  (dpk_fh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0.1, 0})) / 0.2;
 
-              max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
+                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
+              }
+              REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_z_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const double reconstructed_slope =
-                (dpk_fh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0, 0.1})) / 0.2;
+            {
+              double max_z_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const double reconstructed_slope =
+                  (dpk_fh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0, 0.1})) / 0.2;
 
-              max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - 2.1));
+                max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - 2.1));
+              }
+              REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-14));
           }
         }
       }
-    }
 
-    SECTION("R^3 data")
-    {
-      for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
-          auto& mesh  = *p_mesh;
+      SECTION("R^3 data")
+      {
+        for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+            auto& mesh  = *p_mesh;
 
-          auto R3_affine = [](const R3& x) -> R3 {
-            return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2],   //
-                      +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2],   //
-                      -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]};
-          };
-          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto R3_affine = [](const R3& x) -> R3 {
+              return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2],   //
+                        +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2],   //
+                        -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]};
+            };
+            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0<R3> uh{p_mesh};
+            DiscreteFunctionP0<R3> uh{p_mesh};
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
 
-          auto reconstructions = PolynomialReconstruction{}.build(uh);
+            auto reconstructions = PolynomialReconstruction{}.build(1, uh);
 
-          auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
+            auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
-          }
 
-          {
-            double max_x_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
-                                                          dpk_uh[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+            {
+              double max_x_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
+                                                            dpk_uh[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
 
-              max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+                max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+              }
+              REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          {
-            double max_y_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
-                                                          dpk_uh[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+            {
+              double max_y_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
+                                                            dpk_uh[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
 
-              max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
+                max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
+              }
+              REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          {
-            double max_z_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
-                                                          dpk_uh[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+            {
+              double max_z_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
+                                                            dpk_uh[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
 
-              max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R3{1.8, -3.7, 1.9}));
+                max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R3{1.8, -3.7, 1.9}));
+              }
+              REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-13));
           }
         }
       }
-    }
 
-    SECTION("R^2x2 data")
-    {
-      using R2x2 = TinyMatrix<2, 2>;
+      SECTION("R^2x2 data")
+      {
+        using R2x2 = TinyMatrix<2, 2>;
 
-      for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
-          auto& mesh  = *p_mesh;
+        for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+            auto& mesh  = *p_mesh;
 
-          auto R2x2_affine = [](const R3& x) -> R2x2 {
-            return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
-                        //
-                        +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
-          };
-          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto R2x2_affine = [](const R3& x) -> R2x2 {
+              return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
+                          //
+                          +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
+            };
+            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0<R2x2> Ah{p_mesh};
+            DiscreteFunctionP0<R2x2> Ah{p_mesh};
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
 
-          auto reconstructions = PolynomialReconstruction{}.build(Ah);
+            auto reconstructions = PolynomialReconstruction{}.build(1, Ah);
 
-          auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
+            auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              max_mean_error =
-                std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                max_mean_error =
+                  std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
+              }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          {
-            double max_x_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
-                                                            dpk_Ah[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+            {
+              double max_x_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
+                                                              dpk_Ah[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
 
-              max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, 2.1,   //
-                                                                                                       -2.3, +3.1}));
+                max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, 2.1,   //
+                                                                                                         -2.3, +3.1}));
+              }
+              REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
             }
-            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
-          }
 
-          {
-            double max_y_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
-                                                            dpk_Ah[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+            {
+              double max_y_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
+                                                              dpk_Ah[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
 
-              max_y_slope_error = std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2,   //
-                                                                                                       1.3, +0.8}));
+                max_y_slope_error =
+                  std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2,   //
+                                                                                       1.3, +0.8}));
+              }
+              REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
             }
-            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
-          }
 
-          {
-            double max_z_slope_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
-                                                            dpk_Ah[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+            {
+              double max_z_slope_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
+                                                              dpk_Ah[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
 
-              max_z_slope_error = std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R2x2{-1.3, -2.4,   //
-                                                                                                       +1.4, -1.8}));
+                max_z_slope_error =
+                  std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R2x2{-1.3, -2.4,   //
+                                                                                       +1.4, -1.8}));
+              }
+              REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12));
             }
-            REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12));
           }
         }
       }
-    }
 
-    SECTION("vector data")
-    {
-      using R4 = TinyVector<4>;
+      SECTION("vector data")
+      {
+        using R4 = TinyVector<4>;
 
-      for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
-        SECTION(named_mesh.name())
-        {
-          auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
-          auto& mesh  = *p_mesh;
+        for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+          SECTION(named_mesh.name())
+          {
+            auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+            auto& mesh  = *p_mesh;
 
-          auto vector_affine = [](const R3& x) -> R4 {
-            return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
-                      //
-                      +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
-          };
-          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto vector_affine = [](const R3& x) -> R4 {
+              return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
+                        //
+                        +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
+            };
+            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
-          DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
+            DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
 
-          parallel_for(
-            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-              auto vector = vector_affine(xj[cell_id]);
-              for (size_t i = 0; i < vector.dimension(); ++i) {
-                Vh[cell_id][i] = vector[i];
-              }
-            });
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                auto vector = vector_affine(xj[cell_id]);
+                for (size_t i = 0; i < vector.dimension(); ++i) {
+                  Vh[cell_id][i] = vector[i];
+                }
+              });
 
-          auto reconstructions = PolynomialReconstruction{}.build(Vh);
+            auto reconstructions = PolynomialReconstruction{}.build(1, Vh);
 
-          auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
+            auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
 
-          {
-            double max_mean_error = 0;
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              auto vector = vector_affine(xj[cell_id]);
-              for (size_t i = 0; i < vector.dimension(); ++i) {
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+            {
+              double max_mean_error = 0;
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                auto vector = vector_affine(xj[cell_id]);
+                for (size_t i = 0; i < vector.dimension(); ++i) {
+                  max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+                }
               }
+              REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
             }
-            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
-          }
 
-          {
-            double max_x_slope_error = 0;
-            const R4 slope{+1.7, 2.1, -2.3, +3.1};
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              for (size_t i = 0; i < slope.dimension(); ++i) {
-                const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0.1, 0, 0} + xj[cell_id]) -
-                                                                dpk_Vh(cell_id, i)(xj[cell_id] - R3{0.1, 0, 0}));
+            {
+              double max_x_slope_error = 0;
+              const R4 slope{+1.7, 2.1, -2.3, +3.1};
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < slope.dimension(); ++i) {
+                  const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0.1, 0, 0} + xj[cell_id]) -
+                                                                  dpk_Vh(cell_id, i)(xj[cell_id] - R3{0.1, 0, 0}));
 
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+                  max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+                }
               }
+              REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
             }
-            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
-          }
 
-          {
-            double max_y_slope_error = 0;
-            const R4 slope{+1.2, -2.2, 1.3, +0.8};
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              for (size_t i = 0; i < slope.dimension(); ++i) {
-                const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0.1, 0} + xj[cell_id]) -
-                                                                dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0.1, 0}));
+            {
+              double max_y_slope_error = 0;
+              const R4 slope{+1.2, -2.2, 1.3, +0.8};
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < slope.dimension(); ++i) {
+                  const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0.1, 0} + xj[cell_id]) -
+                                                                  dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0.1, 0}));
 
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+                  max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+                }
               }
+              REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
             }
-            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
-          }
 
-          {
-            double max_z_slope_error = 0;
-            const R4 slope{-1.3, -2.4, +1.4, -1.8};
-            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-              for (size_t i = 0; i < slope.dimension(); ++i) {
-                const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0, 0.1} + xj[cell_id]) -
-                                                                dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0, 0.1}));
+            {
+              double max_z_slope_error = 0;
+              const R4 slope{-1.3, -2.4, +1.4, -1.8};
+              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                for (size_t i = 0; i < slope.dimension(); ++i) {
+                  const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0, 0.1} + xj[cell_id]) -
+                                                                  dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0, 0.1}));
 
-                max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - slope[i]));
+                  max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - slope[i]));
+                }
               }
+              REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12));
             }
-            REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12));
           }
         }
       }
     }
-  }
 
-  SECTION("errors")
-  {
-    auto p_mesh1 = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
-    DiscreteFunctionP0<double> f1{p_mesh1};
+    SECTION("errors")
+    {
+      auto p_mesh1 = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+      DiscreteFunctionP0<double> f1{p_mesh1};
 
-    auto p_mesh2 = MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
-    DiscreteFunctionP0<double> f2{p_mesh2};
+      auto p_mesh2 = MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
+      DiscreteFunctionP0<double> f2{p_mesh2};
 
-    REQUIRE_THROWS_WITH(PolynomialReconstruction{}.build(f1, f2),
-                        "error: cannot reconstruct functions living of different meshes simultaneously");
+      REQUIRE_THROWS_WITH(PolynomialReconstruction{}.build(1, f1, f2),
+                          "error: cannot reconstruct functions living of different meshes simultaneously");
+    }
   }
 }