diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 0bc65fc85c5e8f9a650a212c0517a2c0d21b8e6b..1a2b1e072734154516b01001cb9831bfaff8ea73 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -1085,10 +1085,32 @@ PolynomialReconstruction::_build(
 
                   for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
                        ++i_symmetry) {
-                    auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
-                    auto ghost_cell_list = ghost_stencil[cell_j_id];
-                    if (ghost_cell_list.size() > 0) {
-                      throw NotImplementedError("TinyVector symmetries for reconstruction of DiscreteFunctionP0Vector");
+                    if constexpr (DataType::Dimension == MeshType::Dimension) {
+                      auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+                      auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+                      const Rd& normal = symmetry_normal_list[i_symmetry];
+
+                      for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+                        const CellId cell_i_id = ghost_cell_list[i];
+
+                        for (size_t l = 0; l < qj_vector.size(); ++l) {
+                          const DataType& qj    = qj_vector[l];
+                          const DataType& qi    = discrete_function[cell_i_id][l];
+                          const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
+                          for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
+                               ++k, ++kB) {
+                            B(index, kB) = qi_qj[k];
+                          }
+                        }
+                      }
+                    } else {
+                      // LCOV_EXCL_START
+                      std::stringstream error_msg;
+                      error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+                                << " using a mesh of dimension " << MeshType::Dimension;
+                      throw UnexpectedError(error_msg.str());
+                      // LCOV_EXCL_STOP
                     }
                   }
                 } else if constexpr (is_tiny_matrix_v<DataType>) {
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 690030bb8aaafde19d6a7f5129a5d676f4d599c7..6136f26ac21e62e7d8a07449a8f54a8ca3725849 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -275,7 +275,7 @@ add_executable (mpi_unit_tests
   test_OFStream.cpp
   test_ParallelChecker_read.cpp
   test_Partitioner.cpp
-  test_PolynomialReconstruction.cpp
+  test_PolynomialReconstruction_degree_1.cpp
   test_PolynomialReconstructionDescriptor.cpp
   test_RandomEngine.cpp
   test_StencilBuilder_cell2cell.cpp
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction_degree_1.cpp
similarity index 52%
rename from tests/test_PolynomialReconstruction.cpp
rename to tests/test_PolynomialReconstruction_degree_1.cpp
index 3fd31d9527da6954cdbcd74bebc86871c01599cc..7c4f07266514eaf8b277af5f26751c2ec17bc8ae 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction_degree_1.cpp
@@ -9,8 +9,13 @@
 
 #include <algebra/SmallMatrix.hpp>
 #include <algebra/SmallVector.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureFormula.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <mesh/NamedBoundaryDescriptor.hpp>
 #include <scheme/DiscreteFunctionDPkVariant.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
@@ -20,14 +25,14 @@
 
 // clazy:excludeall=non-pod-global-static
 
-TEST_CASE("PolynomialReconstruction", "[scheme]")
+TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
 {
-  SECTION("degree 1")
+  SECTION("without symmetries")
   {
-    std::vector<PolynomialReconstructionDescriptor> descriptor_list = {
-      PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1},
-      PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1},
-    };
+    std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+      {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1},
+       PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1},
+       PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, 1}};
 
     for (auto descriptor : descriptor_list) {
       SECTION(name(descriptor.integrationMethodType()))
@@ -46,6 +51,9 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
 
                 auto R_affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
                 auto xj       = MeshDataManager::instance().getMeshData(mesh).xj();
+                auto xr       = mesh.xr();
+
+                auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
                 DiscreteFunctionP0<double> fh{p_mesh};
 
@@ -56,6 +64,21 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
 
                 auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
 
+                auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{2});
+                {
+                  double max_error = 0;
+                  for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                    auto cell_nodes = cell_to_node_matrix[cell_id];
+                    LineTransformation<1> T{xr[cell_nodes[0]], xr[cell_nodes[1]]};
+                    for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+                      max_error = std::max(max_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+                    }
+                  }
+
+                  REQUIRE(parallel::allReduceMax(max_error) == 0   // Catch::Approx(0).margin(0 * 1E-14)
+                  );
+                }
+
                 {
                   double max_mean_error = 0;
                   for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
@@ -1056,4 +1079,834 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
       }
     }
   }
+
+  SECTION("with symmetries")
+  {
+    SECTION("errors")
+    {
+      SECTION("1D")
+      {
+        auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+        PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, 1,
+                                                      std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+                                                        std::make_shared<NamedBoundaryDescriptor>("XMIN")}};
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<2>>{p_mesh}),
+                            "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 1");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<3>>{p_mesh}),
+                            "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 1");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyVector<2>>{p_mesh, 1}),
+                            "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 1");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyVector<3>>{p_mesh, 1}),
+                            "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 1");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<2>>{p_mesh}),
+                            "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 1");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<3>>{p_mesh}),
+                            "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 1");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyMatrix<2>>{p_mesh, 1}),
+                            "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 1");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyMatrix<3>>{p_mesh, 1}),
+                            "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 1");
+      }
+
+      SECTION("2D")
+      {
+        auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+
+        PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, 1,
+                                                      std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+                                                        std::make_shared<NamedBoundaryDescriptor>("XMIN")}};
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<1>>{p_mesh}),
+                            "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 2");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<3>>{p_mesh}),
+                            "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 2");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyVector<1>>{p_mesh, 1}),
+                            "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 2");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyVector<3>>{p_mesh, 1}),
+                            "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 2");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<1>>{p_mesh}),
+                            "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 2");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<3>>{p_mesh}),
+                            "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 2");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyMatrix<1>>{p_mesh, 1}),
+                            "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 2");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyMatrix<3>>{p_mesh, 1}),
+                            "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 2");
+      }
+
+      SECTION("3D")
+      {
+        auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+
+        PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, 1,
+                                                      std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+                                                        std::make_shared<NamedBoundaryDescriptor>("XMIN")}};
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<1>>{p_mesh}),
+                            "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 3");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<2>>{p_mesh}),
+                            "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 3");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyVector<1>>{p_mesh, 1}),
+                            "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 3");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyVector<2>>{p_mesh, 1}),
+                            "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 3");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<1>>{p_mesh}),
+                            "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 3");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<2>>{p_mesh}),
+                            "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 3");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyMatrix<1>>{p_mesh, 1}),
+                            "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 3");
+
+        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
+                              DiscreteFunctionP0Vector<TinyMatrix<2>>{p_mesh, 1}),
+                            "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 3");
+      }
+    }
+
+    std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+      {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1,
+                                          std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+                                            std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
+       PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1,
+                                          std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+                                            std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
+       PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, 1,
+                                          std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+                                            std::make_shared<NamedBoundaryDescriptor>("XMIN")}}};
+
+    for (auto descriptor : descriptor_list) {
+      SECTION(name(descriptor.integrationMethodType()))
+      {
+        SECTION("1D")
+        {
+          using R1 = TinyVector<1>;
+
+          SECTION("R^1 data")
+          {
+            auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+            auto& mesh = *p_mesh;
+
+            auto R1_affine = [](const R1& x) { return R1{1.7 * (x[0] + 1)}; };
+            auto xj        = MeshDataManager::instance().getMeshData(mesh).xj();
+
+            DiscreteFunctionP0<R1> fh{p_mesh};
+
+            parallel_for(
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R1_affine(xj[cell_id]); });
+
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+
+            auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1>>();
+
+            {
+              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]) - R1_affine(xj[cell_id]))[0]));
+              }
+              REQUIRE(parallel::allReduceMax(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] / 0.2;
+
+                max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
+              }
+              REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
+            }
+          }
+
+          SECTION("R1 vector 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 vector_affine0 = [](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1)}; };
+
+                auto vector_affine1 = [](const R1& x) -> R1 { return R1{-0.3 * (x[0] + 1)}; };
+
+                auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+                DiscreteFunctionP0Vector<R1> Vh{p_mesh, 2};
+
+                parallel_for(
+                  mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+                    Vh[cell_id][0] = vector_affine0(xj[cell_id]);
+                    Vh[cell_id][1] = vector_affine1(xj[cell_id]);
+                  });
+
+                auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+                auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R1>>();
+
+                {
+                  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_Vh(cell_id, 0)(xj[cell_id]) - vector_affine0(xj[cell_id])));
+                    max_mean_error =
+                      std::max(max_mean_error, l2Norm(dpk_Vh(cell_id, 1)(xj[cell_id]) - vector_affine1(xj[cell_id])));
+                  }
+                  REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+                }
+
+                {
+                  double max_slope_error = 0;
+                  {
+                    const TinyVector<1> slope0{+1.7};
+
+                    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                      for (size_t i = 0; i < R1::Dimension; ++i) {
+                        const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R1{0.1} + xj[cell_id])[i] -
+                                                                        dpk_Vh(cell_id, 0)(xj[cell_id] - R1{0.1})[i]);
+
+                        max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope0[i]));
+                      }
+                    }
+                  }
+
+                  {
+                    const TinyVector<1> slope1{-0.3};
+
+                    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+                      for (size_t i = 0; i < R1::Dimension; ++i) {
+                        const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R1{0.1} + xj[cell_id])[i] -
+                                                                        dpk_Vh(cell_id, 1)(xj[cell_id] - R1{0.1})[i]);
+
+                        max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope1[i]));
+                      }
+                    }
+                  }
+                  REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
+                }
+              }
+            }
+          }
+        }
+
+        SECTION("2D")
+        {
+#warning not done
+          using R2 = TinyVector<2>;
+
+          // 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;
+
+          //       auto R_affine = [](const R2& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1]; };
+          //       auto xj       = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          //       DiscreteFunctionP0<double> fh{p_mesh};
+
+          //       parallel_for(
+          //         mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
+
+          //       auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+
+          //       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(parallel::allReduceMax(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;
+
+          //           max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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 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)));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-14));
+          //       }
+          //     }
+          //   }
+          // }
+
+          // 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;
+
+          //       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};
+
+          //       parallel_for(
+          //         mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+
+          //       auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+          //       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])));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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}));
+
+          //           max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7,
+          //           -0.6, 3.1}));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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](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}));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
+          //       }
+          //     }
+          //   }
+          // }
+
+          // 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;
+
+          //       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};
+
+          //       parallel_for(
+          //         mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]);
+          //         });
+
+          //       auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+          //       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])));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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 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}));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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 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}));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
+          //       }
+          //     }
+          //   }
+          // }
+
+          // 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;
+
+          //       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};
+
+          //       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{descriptor}.build(Vh);
+
+          //       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]));
+          //           }
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+          //       }
+
+          //       {
+          //         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]));
+          //           }
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+          //       }
+
+          //       {
+          //         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]));
+          //           }
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
+          //       }
+          //     }
+          //   }
+          // }
+        }
+
+        SECTION("3D")
+        {
+          using R3 = TinyVector<3>;
+#warning not done
+
+          // 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;
+
+          //       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();
+
+          //       DiscreteFunctionP0<double> fh{p_mesh};
+
+          //       parallel_for(
+          //         mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
+
+          //       auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+
+          //       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(parallel::allReduceMax(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;
+
+          //           max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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 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)));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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 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));
+          //         }//         REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
+          //       }
+          //     }
+          //   }
+          // }
+
+          // 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();
+
+          //       DiscreteFunctionP0<R3> uh{p_mesh};
+
+          //       parallel_for(
+          //         mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+
+          //       auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+          //       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])));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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}));
+
+          //           max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7,
+          //           -0.6, 3.1}));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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 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}));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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 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}));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
+          //       }
+          //     }
+          //   }
+          // }
+
+          // 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;
+
+          //       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};
+
+          //       parallel_for(
+          //         mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]);
+          //         });
+
+          //       descriptor.setRowWeighting(false);
+          //       auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+          //       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])));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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 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}));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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 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}));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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}));
+
+          //           max_z_slope_error =
+          //             std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R2x2{-1.3, -2.4,   //
+          //                                                                                  +1.4, -1.8}));
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
+          //       }
+          //     }
+          //   }
+          // }
+
+          // 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;
+
+          //       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};
+
+          //       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];
+          //           }
+          //         });
+
+          //       descriptor.setPreconditioning(false);
+          //       auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+          //       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]));
+          //           }
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+          //       }
+
+          //       {
+          //         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]));
+          //           }
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+          //       }
+
+          //       {
+          //         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]));
+          //           }
+          //         }
+          //         REQUIRE(parallel::allReduceMax(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}));
+
+          //             max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - slope[i]));
+          //           }
+          //         }
+          //         REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-13));
+          //       }
+          //     }
+          //   }
+          // }
+        }
+      }
+    }
+  }
 }