diff --git a/tests/test_PolynomialReconstruction_degree_1.cpp b/tests/test_PolynomialReconstruction_degree_1.cpp
index 927a5f4c3f965e53c1d940f012dc1bbddf2fbbe1..4e8b8ca57d4266acb26481e91c68c83d6139c951 100644
--- a/tests/test_PolynomialReconstruction_degree_1.cpp
+++ b/tests/test_PolynomialReconstruction_degree_1.cpp
@@ -12,7 +12,10 @@
 #include <analysis/GaussQuadratureDescriptor.hpp>
 #include <analysis/QuadratureFormula.hpp>
 #include <analysis/QuadratureManager.hpp>
+#include <geometry/CubeTransformation.hpp>
 #include <geometry/LineTransformation.hpp>
+#include <geometry/PrismTransformation.hpp>
+#include <geometry/PyramidTransformation.hpp>
 #include <geometry/SquareTransformation.hpp>
 #include <geometry/TetrahedronTransformation.hpp>
 #include <geometry/TriangleTransformation.hpp>
@@ -104,16 +107,36 @@ max_reconstruction_error(const MeshType& mesh,
         }
         break;
       }
-      // case CellType::Pyramid: {
-      // SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
-      //                                             xr[cell_nodes[3]]};
-      // auto qf = QuadratureManager::instance().getSquareFormula(GaussLobattoQuadratureDescriptor{dpk_f.degree() + 1});
-      // for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
-      //   auto x    = T(qf.point(i_quadrarture));
-      //   max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
-      // }
-      // break;
-      // }
+      case CellType::Pyramid: {
+        PyramidTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
+                                xr[cell_nodes[4]]};
+        auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+        for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+          auto x    = T(qf.point(i_quadrarture));
+          max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+        }
+        break;
+      }
+      case CellType::Prism: {
+        PrismTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+                              xr[cell_nodes[3]], xr[cell_nodes[4]], xr[cell_nodes[5]]};
+        auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+        for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+          auto x    = T(qf.point(i_quadrarture));
+          max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+        }
+        break;
+      }
+      case CellType::Hexahedron: {
+        CubeTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
+                             xr[cell_nodes[4]], xr[cell_nodes[5]], xr[cell_nodes[6]], xr[cell_nodes[7]]};
+        auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+        for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+          auto x    = T(qf.point(i_quadrarture));
+          max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+        }
+        break;
+      }
       default: {
         throw UnexpectedError("unexpected cell type");
       }
@@ -194,16 +217,45 @@ max_reconstruction_error(
         }
         break;
       }
-      // case CellType::Pyramid: {
-      // SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
-      //                                             xr[cell_nodes[3]]};
-      // auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
-      // for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
-      //   auto x    = T(qf.point(i_quadrarture));
-      //   max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
-      // }
-      // break;
-      // }
+      case CellType::Pyramid: {
+        PyramidTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
+                                xr[cell_nodes[4]]};
+        auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+        for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+          auto x = T(qf.point(i_quadrarture));
+          for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+            max_error =
+              std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+          }
+        }
+        break;
+      }
+      case CellType::Prism: {
+        PrismTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+                              xr[cell_nodes[3]], xr[cell_nodes[4]], xr[cell_nodes[5]]};
+        auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+        for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+          auto x = T(qf.point(i_quadrarture));
+          for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+            max_error =
+              std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+          }
+        }
+        break;
+      }
+      case CellType::Hexahedron: {
+        CubeTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
+                             xr[cell_nodes[4]], xr[cell_nodes[5]], xr[cell_nodes[6]], xr[cell_nodes[7]]};
+        auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+        for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+          auto x = T(qf.point(i_quadrarture));
+          for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+            max_error =
+              std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+          }
+        }
+        break;
+      }
       default: {
         throw UnexpectedError("unexpected cell type");
       }
@@ -366,12 +418,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
                 auto& mesh  = *p_mesh;
 
                 std::array<std::function<R3(const R1&)>, 2> vector_exact =
-                  {[](const R1& x) -> R3 {
-                     return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
-                   },
-                   [](const R1& x) -> R3 {
-                     return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]};
-                   }};
+                  {[](const R1& x) -> R3 { return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]}; },
+                   [](const R1& x) -> R3 { return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]}; }};
 
                 auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
@@ -619,62 +667,20 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
                 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();
+                auto R_exact = [](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]); });
+                  mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_exact(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));
-                }
+                double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact));
+                REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
               }
             }
           }
@@ -687,7 +693,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
                 auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
                 auto& mesh  = *p_mesh;
 
-                auto R3_affine = [](const R3& x) -> R3 {
+                auto R3_exact = [](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]};
@@ -697,53 +703,14 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
                 DiscreteFunctionP0<R3> uh{p_mesh};
 
                 parallel_for(
-                  mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+                  mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_exact(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));
-                }
+                double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact));
+                REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
               }
             }
           }
@@ -758,7 +725,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
                 auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
                 auto& mesh  = *p_mesh;
 
-                auto R2x2_affine = [](const R3& x) -> R2x2 {
+                auto R2x2_exact = [](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]};
@@ -768,88 +735,40 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
                 DiscreteFunctionP0<R2x2> Ah{p_mesh};
 
                 parallel_for(
-                  mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
+                  mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_exact(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));
-                }
+                auto dpk_Ah      = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
+                double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R2x2_exact));
+                REQUIRE(parallel::allReduceMax(max_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]};
-                };
+                std::array<std::function<double(const R3&)>, 4> vector_exact =
+                  {[](const R3& x) -> double { return +2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2]; },
+                   [](const R3& x) -> double { return -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2]; },
+                   [](const R3& x) -> double { return +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2]; },
+                   [](const R3& x) -> double { return -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];
+                    for (size_t i = 0; i < vector_exact.size(); ++i) {
+                      Vh[cell_id][i] = vector_exact[i](xj[cell_id]);
                     }
                   });
 
@@ -858,58 +777,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
 
                 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));
-                }
+                double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+                REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
               }
             }
           }
@@ -1070,37 +939,20 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
 
             auto& mesh = *p_mesh;
 
-            auto R1_affine = [](const R1& x) { return R1{1.7 * (x[0] + 1)}; };
-            auto xj        = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto R1_exact = [](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]); });
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R1_exact(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));
-            }
+            double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R1_exact));
+            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
           }
 
           SECTION("R1 vector data")
@@ -1111,9 +963,9 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
                 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)}; };
+                std::array<std::function<R1(const R1&)>, 2> vector_exact   //
+                  = {[](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1)}; },
+                     [](const R1& x) -> R1 { return R1{-0.3 * (x[0] + 1)}; }};
 
                 auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
@@ -1121,54 +973,16 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
 
                 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]);
+                    Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
+                    Vh[cell_id][1] = vector_exact[1](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));
-                }
+                double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+                REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
               }
             }
           }
@@ -1211,128 +1025,47 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
             auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
             auto& mesh  = *p_mesh;
 
-            auto R2_affine = [](const R2& x) -> R2 { return R2{2.3 * (x[0] - 2), -1.3 * (x[1] - 1)}; };
-            auto xj        = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto R2_exact = [](const R2& x) -> R2 { return R2{2.3 * (x[0] - 2), -1.3 * (x[1] - 1)}; };
+            auto xj       = MeshDataManager::instance().getMeshData(mesh).xj();
 
             DiscreteFunctionP0<R2> fh{p_mesh};
 
             parallel_for(
-              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R2_affine(xj[cell_id]); });
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R2_exact(xj[cell_id]); });
 
             auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
 
             auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2>>();
 
-            {
-              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_fh[cell_id](xj[cell_id]) - R2_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 R2 reconstructed_slope =
-                  1. / 0.2 * (dpk_fh[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0.1, 0}));
-
-                max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R2{2.3, 0}));
-              }
-              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 R2 reconstructed_slope =
-                  1 / 0.2 * (dpk_fh[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0, 0.1}));
-
-                max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - (R2{0, -1.3})));
-              }
-              REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
-            }
+            double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R2_exact));
+            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
           }
 
           SECTION("vector of R2")
           {
-            using R4 = TinyVector<4>;
-
             auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
             auto& mesh  = *p_mesh;
 
-            auto vector_affine = [](const R2& x) -> R4 {
-              return R4{+1.7 * (x[0] - 2),   //
-                        -0.6 * (x[1] - 1),   //
-                        -2.3 * (x[0] - 2),   //
-                        +1.1 * (x[1] - 1)};
-            };
+            std::array<std::function<R2(const R2&)>, 2> vector_exact   //
+              = {[](const R2& x) -> R2 { return R2{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1)}; },
+                 [](const R2& x) -> R2 { return R2{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1)}; }};
+
             auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
             DiscreteFunctionP0Vector<R2> Vh{p_mesh, 2};
 
             parallel_for(
               mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-                auto vector       = vector_affine(xj[cell_id]);
-                Vh[cell_id][0][0] = vector[0];
-                Vh[cell_id][0][1] = vector[1];
-                Vh[cell_id][1][0] = vector[2];
-                Vh[cell_id][1][1] = vector[3];
+                Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
+                Vh[cell_id][1] = vector_exact[1](xj[cell_id]);
               });
 
             auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
 
             auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2>>();
 
-            {
-              double max_mean_error = 0;
-              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-                auto vector = vector_affine(xj[cell_id]);
-
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[0] - vector[0]));
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[1] - vector[1]));
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[0] - vector[2]));
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[1] - vector[3]));
-              }
-              REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
-            }
-
-            {
-              double max_x_slope_error = 0;
-              const R4 slope{+1.7, 0, -2.3, 0};
-              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-                const R2 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R2{0.1, 0} + xj[cell_id]) -
-                                                             dpk_Vh(cell_id, 0)(xj[cell_id] - R2{0.1, 0}));
-
-                const R2 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R2{0.1, 0} + xj[cell_id]) -
-                                                             dpk_Vh(cell_id, 1)(xj[cell_id] - R2{0.1, 0}));
-
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[0] - slope[2]));
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[1] - slope[3]));
-              }
-              REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
-            }
-
-            {
-              double max_y_slope_error = 0;
-              const R4 slope{0, -0.6, 0, 1.1};
-              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-                const R2 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R2{0, 0.1} + xj[cell_id]) -
-                                                             dpk_Vh(cell_id, 0)(xj[cell_id] - R2{0, 0.1}));
-
-                const R2 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R2{0, 0.1} + xj[cell_id]) -
-                                                             dpk_Vh(cell_id, 1)(xj[cell_id] - R2{0, 0.1}));
-
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[0] - slope[2]));
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[1] - slope[3]));
-              }
-              REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
-            }
+            double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
           }
         }
       }
@@ -1379,174 +1112,47 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
             auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
             auto& mesh  = *p_mesh;
 
-            auto R3_affine = [](const R3& x) -> R3 {
-              return R3{2.3 * (x[0] - 2), -1.3 * (x[1] - 1), 1.4 * (x[2] - 1)};
-            };
-            auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+            auto R3_exact = [](const R3& x) -> R3 { return R3{2.3 * (x[0] - 2), -1.3 * (x[1] - 1), 1.4 * (x[2] - 1)}; };
+            auto xj       = MeshDataManager::instance().getMeshData(mesh).xj();
 
             DiscreteFunctionP0<R3> fh{p_mesh};
 
             parallel_for(
-              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R3_affine(xj[cell_id]); });
+              mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R3_exact(xj[cell_id]); });
 
             auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
 
             auto dpk_fh = 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_fh[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_fh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
-
-                max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{2.3, 0, 0}));
-              }
-              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_fh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
-
-                max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - (R3{0, -1.3, 0})));
-              }
-              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_fh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
-
-                max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - (R3{0, 0, 1.4})));
-              }
-              REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
-            }
+            double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R3_exact));
+            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
           }
 
-          SECTION("vector of R2")
+          SECTION("vector of R3")
           {
-            using R6 = TinyVector<6>;
-
             auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
             auto& mesh  = *p_mesh;
 
-            auto vector_affine = [](const R3& x) -> R6 {
-              return R6{+1.7 * (x[0] - 2),   //
-                        -0.6 * (x[1] - 1),   //
-                        +1.2 * (x[2] - 1),   //
-                        -2.3 * (x[0] - 2),   //
-                        +1.1 * (x[1] - 1),   //
-                        -0.3 * (x[2] - 1)};
-            };
+            std::array<std::function<R3(const R3&)>, 2> vector_exact   //
+              = {[](const R3& x) -> R3 { return R3{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1), +1.2 * (x[2] - 1)}; },
+                 [](const R3& x) -> R3 { return R3{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1), -0.3 * (x[2] - 1)}; }};
+
             auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
             DiscreteFunctionP0Vector<R3> Vh{p_mesh, 2};
 
             parallel_for(
               mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-                auto vector       = vector_affine(xj[cell_id]);
-                Vh[cell_id][0][0] = vector[0];
-                Vh[cell_id][0][1] = vector[1];
-                Vh[cell_id][0][2] = vector[2];
-                Vh[cell_id][1][0] = vector[3];
-                Vh[cell_id][1][1] = vector[4];
-                Vh[cell_id][1][2] = vector[5];
+                Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
+                Vh[cell_id][1] = vector_exact[1](xj[cell_id]);
               });
 
             auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
 
             auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3>>();
 
-            {
-              double max_mean_error = 0;
-              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-                auto vector = vector_affine(xj[cell_id]);
-
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[0] - vector[0]));
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[1] - vector[1]));
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[2] - vector[2]));
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[0] - vector[3]));
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[1] - vector[4]));
-                max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[2] - vector[5]));
-              }
-              REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
-            }
-
-            {
-              double max_x_slope_error = 0;
-              const R6 slope{+1.7, 0, 0, -2.3, 0, 0};
-              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-                const R3 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R3{0.1, 0, 0} + xj[cell_id]) -
-                                                             dpk_Vh(cell_id, 0)(xj[cell_id] - R3{0.1, 0, 0}));
-
-                const R3 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R3{0.1, 0, 0} + xj[cell_id]) -
-                                                             dpk_Vh(cell_id, 1)(xj[cell_id] - R3{0.1, 0, 0}));
-
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[2] - slope[2]));
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[0] - slope[3]));
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[1] - slope[4]));
-                max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[2] - slope[5]));
-              }
-              REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
-            }
-
-            {
-              double max_y_slope_error = 0;
-              const R6 slope{0, -0.6, 0, 0, 1.1, 0};
-              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-                const R3 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R3{0, 0.1, 0} + xj[cell_id]) -
-                                                             dpk_Vh(cell_id, 0)(xj[cell_id] - R3{0, 0.1, 0}));
-
-                const R3 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R3{0, 0.1, 0} + xj[cell_id]) -
-                                                             dpk_Vh(cell_id, 1)(xj[cell_id] - R3{0, 0.1, 0}));
-
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[2] - slope[2]));
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[0] - slope[3]));
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[1] - slope[4]));
-                max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[2] - slope[5]));
-              }
-              REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
-            }
-
-            {
-              double max_z_slope_error = 0;
-              const R6 slope{0, 0, 1.2, 0, 0, -0.3};
-              for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-                const R3 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R3{0, 0, 0.1} + xj[cell_id]) -
-                                                             dpk_Vh(cell_id, 0)(xj[cell_id] - R3{0, 0, 0.1}));
-
-                const R3 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R3{0, 0, 0.1} + xj[cell_id]) -
-                                                             dpk_Vh(cell_id, 1)(xj[cell_id] - R3{0, 0, 0.1}));
-
-                max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
-                max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
-                max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope0[2] - slope[2]));
-                max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope1[0] - slope[3]));
-                max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope1[1] - slope[4]));
-                max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope1[2] - slope[5]));
-              }
-              REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-13));
-            }
+            double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
           }
         }
       }