diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index b01d11a792bcd250857758c52d1c4885b78d896f..2201e5f0f796047bdc3e1ab303ac6a614c816ba0 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -1243,6 +1243,7 @@ PolynomialReconstruction::_build(
 
         } else if ((m_descriptor.integrationMethodType() == IntegrationMethodType::element) or
                    (m_descriptor.integrationMethodType() == IntegrationMethodType::boundary)) {
+#warning implement boundary based reconstruction for 1d and 3d
           if ((m_descriptor.integrationMethodType() == IntegrationMethodType::boundary) and
               (MeshType::Dimension == 2)) {
             if constexpr (MeshType::Dimension == 2) {
diff --git a/tests/test_PolynomialReconstruction_degree_3.cpp b/tests/test_PolynomialReconstruction_degree_3.cpp
index d8badda38a09584a34696b7311120f7b859e4b19..e20fa89a3bdcba29f61a771deacba273c536f28b 100644
--- a/tests/test_PolynomialReconstruction_degree_3.cpp
+++ b/tests/test_PolynomialReconstruction_degree_3.cpp
@@ -34,10 +34,6 @@ TEST_CASE("PolynomialReconstruction_degree_3", "[scheme]")
        PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree}};
 
     for (auto descriptor : descriptor_list) {
-#warning remove
-      descriptor.setPreconditioning(false);
-      descriptor.setRowWeighting(false);
-
       SECTION(name(descriptor.integrationMethodType()))
       {
         SECTION("1D")
@@ -437,7 +433,7 @@ TEST_CASE("PolynomialReconstruction_degree_3", "[scheme]")
         {
           using R3 = TinyVector<3>;
 
-          auto p0 = [](const R3& X) -> double {
+          auto p = [](const R3& X, const std::array<double, 20>& a) -> double {
             const double x   = X[0];
             const double y   = X[1];
             const double z   = X[2];
@@ -458,730 +454,313 @@ TEST_CASE("PolynomialReconstruction_degree_3", "[scheme]")
             const double y3  = y2 * y;
             const double z3  = z2 * z;
 
-            return 2.3 + 1.7 * x - 1.3 * y + 2.1 * z     //
-                   + 1.7 * x2 + 1.4 * y2 + 1.7 * z2      //
-                   - 2.3 * xy + 1.6 * xz - 1.9 * yz      //
-                   + 1.2 * x3 - 2.1 * y3 - 1.1 * z3      //
-                   - 1.7 * x2y - 1.3 * x2z + 0.9 * xy2   //
-                   - 0.7 * y2z + 1.5 * xz2 - 0.7 * yz2   //
-                   + 2.8 * xyz                           //
+            return a[0] + a[1] * x + a[2] * y + a[3] * z       //
+                   + a[4] * x2 + a[5] * y2 + a[6] * z2         //
+                   + a[7] * xy + a[8] * xz + a[9] * yz         //
+                   + a[10] * x3 + a[11] * y3 + a[12] * z3      //
+                   + a[13] * x2y + a[14] * x2z + a[15] * xy2   //
+                   + a[16] * y2z + a[17] * xz2 + a[18] * yz2   //
+                   + a[19] * xyz                               //
               ;
           };
 
-          // auto p1 = [](const R3& x) {
-          //   return +2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2]                   //
-          //          + 1.7 * x[0] * x[0] - 2.4 * x[1] * x[1] - 2.3 * x[2] * x[2]   //
-          //          - 2.1 * x[0] * x[1] + 2.6 * x[0] * x[2] + 1.6 * x[1] * x[2];
-          // };
+          constexpr std::array<double, 20> a0 = {+2.3, +1.7, -1.3, +2.1, +1.7, +1.4, +1.7, -2.3, +1.6, -1.9,
+                                                 +1.2, -2.1, -1.1, -1.7, -1.3, +0.9, -0.7, +1.5, -0.7, +2.8};
 
-          // auto p2 = [](const R3& x) {
-          //   return +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2]                   //
-          //          + 3.1 * x[0] * x[0] - 1.1 * x[1] * x[1] + 1.7 * x[2] * x[2]   //
-          //          - 2.3 * x[0] * x[1] - 2.6 * x[0] * x[2] - 1.9 * x[1] * x[2];
-          // };
+          auto p0 = [&p, &a0](const R3& X) -> double { return p(X, a0); };
 
-          // auto p3 = [](const R3& x) {
-          //   return -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]                   //
-          //          - 1.5 * x[0] * x[0] + 1.4 * x[1] * x[1] - 1.2 * x[2] * x[2]   //
-          //          - 1.7 * x[0] * x[1] - 1.3 * x[0] * x[2] + 2.1 * x[1] * x[2];
-          // };
+          constexpr std::array<double, 20> a1 = {-1.3, +2.2, +0.1, -2.5, +0.2, -2.3, -1.4, +0.9, +0.2, -0.3,
+                                                 +2.4, -1.2, +1.7, -2.2, +0.6, +1.9, +1.0, -0.8, +2.4, +2.4};
 
-          SECTION("R data")
-          {
-            for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
-              SECTION(named_mesh.name())
-              {
-                auto p_initial_mesh = CartesianMeshBuilder{TinyVector<3>{-0.5, -0.5, -0.5},
-                                                           TinyVector<3>{3.5, 3.5, 3.5}, TinyVector<3, size_t>{4, 4, 4}}
-                                        .mesh()
-                                        ->get<Mesh<3>>();   // named_mesh.mesh()->get<Mesh<3>>();
-                auto& initial_mesh = *p_initial_mesh;
-
-                constexpr double theta = 0;
-                TinyMatrix<3> T0{std::cos(theta), -std::sin(theta), 0,
-                                 //
-                                 std::sin(theta), std::cos(theta), 0,
-                                 //
-                                 0, 0, 1};
-                constexpr double phi = 0;
-                TinyMatrix<3> T1{std::cos(phi), 0, -std::sin(phi),
-                                 //
-                                 0, 1, 0,
-                                 //
-                                 std::sin(phi), 0, std::cos(phi)};
-
-                TinyMatrix T = T0 * T1;
-
-                auto xr = initial_mesh.xr();
-
-                NodeValue<R3> new_xr{initial_mesh.connectivity()};
-                parallel_for(
-                  initial_mesh.numberOfNodes(),
-                  PUGS_LAMBDA(const NodeId node_id) { new_xr[node_id] = T * xr[node_id]; });
-
-                std::shared_ptr p_mesh = std::make_shared<const Mesh<3>>(initial_mesh.shared_connectivity(), new_xr);
-                const Mesh<3>& mesh    = *p_mesh;
-                // inverse rotation
-                TinyMatrix<3> inv_T{std::cos(theta), std::sin(theta), 0,
-                                    //
-                                    -std::sin(theta), std::cos(theta), 0,
-                                    //
-                                    0, 0, 1};
+          auto p1 = [&p, &a1](const R3& X) -> double { return p(X, a1); };
 
-                auto R_exact = p0;
+          constexpr std::array<double, 20> a2 = {+1.9, -1.2, -0.4, -1.2, -0.8, +1.4, +0.5, -1.6, +1.1, -0.7,
+                                                 +0.6, +2.3, -1.8, -1.9, -0.3, -2.4, -1.7, +0.2, -2.4, +1.9};
 
-                DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree + 3, std::function(R_exact));
+          auto p2 = [&p, &a2](const R3& X) -> double { return p(X, a2); };
 
-                auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
-                auto dpk_fh          = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
+          constexpr std::array<double, 20> a3 = {+0.8, +0.5, +1.3, -2.3, +0.9, -0.4, -2.0, +1.8, +0.5, +0.7,
+                                                 +1.0, -0.4, +1.1, +1.8, -0.4, +1.1, -0.0, +1.4, +1.9, -2.2};
 
-                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-12));
-              }
-            }
+          auto p3 = [&p, &a3](const R3& X) -> double { return p(X, a3); };
 
-#warning NOT FINISHED
-          }
+          auto p_mesh = CartesianMeshBuilder{TinyVector<3>{-0.5, -0.5, -0.5}, TinyVector<3>{3.5, 3.5, 3.5},
+                                             TinyVector<3, size_t>{4, 4, 4}}
+                          .mesh()
+                          ->get<Mesh<3>>();
+          const auto& mesh = *p_mesh;
+
+          SECTION("R data")
+          {
+            auto R_exact = p0;
 
-          // 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;
+            DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree + 3, std::function(R_exact));
 
-          //       auto R3_exact = [&](const R3& x) -> R3 { return R3{p1(x), p2(x), p3(x)}; };
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+            auto dpk_fh          = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
 
-          //       DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
+            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-12));
+          }
 
-          //       auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+          SECTION("R^3 data")
+          {
+            auto R3_exact = [&](const R3& x) -> R3 { return R3{p1(x), p2(x), p3(x)}; };
 
-          //       auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
+            DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
 
-          //       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));
-          //     }
-          //   }
-          // }
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
 
-          // SECTION("R^2x2 data")
-          // {
-          //   using R2x2 = TinyMatrix<2, 2>;
+            auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
 
-          //   for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
-          //     SECTION(named_mesh.name())
-          //     {
-          //       auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
-          //       auto& mesh  = *p_mesh;
+            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));
+          }
 
-          //       auto R2x2_exact = [&](const R3& x) -> R2x2 {
-          //         return R2x2{p0(x), p1(x),   //
-          //                     p2(x), p3(x)};
-          //       };
+          SECTION("R^2x2 data")
+          {
+            using R2x2 = TinyMatrix<2, 2>;
 
-          //       DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
+            auto R2x2_exact = [&](const R3& x) -> R2x2 {
+              return R2x2{p0(x), p1(x),   //
+                          p2(x), p3(x)};
+            };
 
-          //       descriptor.setRowWeighting(false);
-          //       auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+            DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
 
-          //       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));
-          //     }
-          //   }
-          // }
+            descriptor.setRowWeighting(false);
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
 
-          // SECTION("vector 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 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));
+          }
 
-          //       std::array<std::function<double(const R3&)>, 4> vector_exact = {p0, p1, p2, p3};
+          SECTION("vector data")
+          {
+            std::array<std::function<double(const R3&)>, 4> vector_exact = {p0, p1, p2, p3};
 
-          //       DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+            DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
 
-          //       descriptor.setPreconditioning(false);
-          //       auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+            descriptor.setPreconditioning(false);
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
 
-          //       auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
+            auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
 
-          //       double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
-          //       REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-          //     }
-          //   }
-          // }
+            double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+          }
         }
       }
     }
   }
 
-  // SECTION("with symmetries")
-  // {
-  //   SECTION("1D")
-  //   {
-  //     std::vector<PolynomialReconstructionDescriptor> descriptor_list =
-  //       {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
-  //                                           std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
-  //                                             std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
-  //        PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
-  //                                           std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
-  //                                             std::make_shared<NamedBoundaryDescriptor>("XMIN")}}};
-  //     using R1 = TinyVector<1>;
-
-  //     auto p0 = [](const R1& x) { return +1.7 * (x[0] + 1) * (x[0] + 1) - 1.1; };
-  //     auto p1 = [](const R1& x) { return -1.2 * (x[0] + 1) * (x[0] + 1) + 1.3; };
-  //     auto p2 = [](const R1& x) { return +1.4 * (x[0] + 1) * (x[0] + 1) - 0.6; };
-
-  //     for (auto descriptor : descriptor_list) {
-  //       SECTION(name(descriptor.integrationMethodType()))
-  //       {
-  //         SECTION("R data")
-  //         {
-  //           auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
-
-  //           auto& mesh = *p_mesh;
+  SECTION("with symmetries")
+  {
+    SECTION("1D")
+    {
+      std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+        {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
+                                            std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+                                              std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
+         PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
+                                            std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+                                              std::make_shared<NamedBoundaryDescriptor>("XMIN")}}};
+      using R1 = TinyVector<1>;
+
+      auto p0 = [](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1) * (x[0] + 1) * (x[0] + 1)}; };
+      auto p1 = [](const R1& x) -> R1 { return R1{-1.2 * (x[0] + 1) * (x[0] + 1) * (x[0] + 1)}; };
+      auto p2 = [](const R1& x) -> R1 { return R1{+1.4 * (x[0] + 1) * (x[0] + 1) * (x[0] + 1)}; };
+
+      for (auto descriptor : descriptor_list) {
+        SECTION(name(descriptor.integrationMethodType()))
+        {
+          SECTION("R1 data")
+          {
+            auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
 
-  //           auto R_exact = p0;
+            auto& mesh = *p_mesh;
 
-  //           DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));
+            auto R1_exact = p0;
 
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+            DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R1_exact));
 
-  //           auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
 
-  //           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-12));
-  //         }
+            auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1>>();
 
-  //         SECTION("R1x1 data")
-  //         {
-  //           using R1x1 = TinyMatrix<1>;
+            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));
+          }
 
-  //           auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+          SECTION("R1 vector data")
+          {
+            auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+            auto& mesh  = *p_mesh;
 
-  //           auto& mesh = *p_mesh;
+            std::array<std::function<R1(const R1&)>, 3> vector_exact   //
+              = {p0, p1, p2};
 
-  //           auto R1x1_exact = [&](const R1& x) { return R1x1{p0(x)}; };
+            DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
 
-  //           DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R1x1_exact));
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
 
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+            auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R1>>();
 
-  //           auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1x1>>();
+            double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+          }
+        }
+      }
+    }
 
-  //           double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R1x1_exact));
-  //           REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //         }
+    SECTION("2D")
+    {
+      std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+        {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
+                                            std::vector<std::shared_ptr<
+                                              const IBoundaryDescriptor>>{std::
+                                                                            make_shared<NamedBoundaryDescriptor>(
+                                                                              "XMAX"),
+                                                                          std::make_shared<NamedBoundaryDescriptor>(
+                                                                            "YMAX")}},
+         PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
+                                            std::vector<std::shared_ptr<
+                                              const IBoundaryDescriptor>>{std::
+                                                                            make_shared<NamedBoundaryDescriptor>(
+                                                                              "XMAX"),
+                                                                          std::make_shared<NamedBoundaryDescriptor>(
+                                                                            "YMAX")}}};
+
+      using R2 = TinyVector<2>;
+
+      auto p_mesh = CartesianMeshBuilder{TinyVector<2>{0, 0}, TinyVector<2>{2, 1}, TinyVector<2, size_t>{3, 3}}
+                      .mesh()
+                      ->get<Mesh<2>>();
+
+      auto& mesh = *p_mesh;
+
+      for (auto descriptor : descriptor_list) {
+        SECTION(name(descriptor.integrationMethodType()))
+        {
+          SECTION("R^2 data")
+          {
+            auto R2_exact = [](const R2& x) -> R2 {
+              return R2{+2.3 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2),   //
+                        -1.3 * (x[1] - 1) * (x[1] - 1) * (x[1] - 1)};
+            };
 
-  //         SECTION("R vector data")
-  //         {
-  //           auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
-  //           auto& mesh  = *p_mesh;
+            DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R2_exact));
 
-  //           std::array<std::function<double(const R1&)>, 3> vector_exact   //
-  //             = {p0, p1, p2};
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
 
-  //           DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+            auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2>>();
 
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+            double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R2_exact));
+            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+          }
 
-  //           auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
-
-  //           double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
-  //           REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //         }
-
-  //         SECTION("R1x1 vector data")
-  //         {
-  //           using R1x1 = TinyMatrix<1>;
-
-  //           auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
-  //           auto& mesh  = *p_mesh;
-
-  //           std::array<std::function<R1x1(const R1&)>, 3> vector_exact   //
-  //             = {[&](const R1& x) { return R1x1{p2(x)}; },               //
-  //                [&](const R1& x) { return R1x1{p0(x)}; },               //
-  //                [&](const R1& x) { return R1x1{p1(x)}; }};
-
-  //           DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
-
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
-
-  //           auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R1x1>>();
-
-  //           double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
-  //           REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //         }
-  //       }
-  //     }
-  //   }
-
-  //   SECTION("2D")
-  //   {
-  //     std::vector<PolynomialReconstructionDescriptor> descriptor_list =
-  //       {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
-  //                                           std::vector<std::shared_ptr<
-  //                                             const IBoundaryDescriptor>>{std::
-  //                                                                           make_shared<NamedBoundaryDescriptor>(
-  //                                                                             "XMAX"),
-  //                                                                         std::make_shared<NamedBoundaryDescriptor>(
-  //                                                                           "YMAX")}},
-  //        PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
-  //                                           std::vector<std::shared_ptr<
-  //                                             const IBoundaryDescriptor>>{std::
-  //                                                                           make_shared<NamedBoundaryDescriptor>(
-  //                                                                             "XMAX"),
-  //                                                                         std::make_shared<NamedBoundaryDescriptor>(
-  //                                                                           "YMAX")}}};
+          SECTION("vector of R2")
+          {
+            std::array<std::function<R2(const R2&)>, 2> vector_exact   //
+              = {[](const R2& x) -> R2 {
+                   return R2{+1.7 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2),   //
+                             -0.6 * (x[1] - 1) * (x[1] - 1) * (x[1] - 1)};
+                 },
+                 [](const R2& x) -> R2 {
+                   return R2{-2.3 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2),   //
+                             +1.1 * (x[1] - 1) * (x[1] - 1) * (x[1] - 1)};
+                 }};
 
-  //     using R2 = TinyVector<2>;
+            DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
 
-  //     auto p_initial_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
-  //     auto& initial_mesh  = *p_initial_mesh;
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
 
-  //     constexpr double theta = 1;
-  //     TinyMatrix<2> T{std::cos(theta), -std::sin(theta),   //
-  //                     std::sin(theta), std::cos(theta)};
+            auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2>>();
 
-  //     auto xr = initial_mesh.xr();
+            double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+          }
+        }
+      }
+    }
 
-  //     NodeValue<R2> new_xr{initial_mesh.connectivity()};
-  //     parallel_for(
-  //       initial_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { new_xr[node_id] = T * xr[node_id]; });
+    SECTION("3D")
+    {
+      std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+        {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
+                                            std::vector<std::shared_ptr<
+                                              const IBoundaryDescriptor>>{std::
+                                                                            make_shared<NamedBoundaryDescriptor>(
+                                                                              "XMAX"),
+                                                                          std::make_shared<NamedBoundaryDescriptor>(
+                                                                            "YMAX"),
+                                                                          std::make_shared<NamedBoundaryDescriptor>(
+                                                                            "ZMAX")}},
+         PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
+                                            std::vector<std::shared_ptr<
+                                              const IBoundaryDescriptor>>{std::
+                                                                            make_shared<NamedBoundaryDescriptor>(
+                                                                              "XMAX"),
+                                                                          std::make_shared<NamedBoundaryDescriptor>(
+                                                                            "YMAX"),
+                                                                          std::make_shared<NamedBoundaryDescriptor>(
+                                                                            "ZMAX")}}};
+
+      using R3 = TinyVector<3>;
+
+      auto p_mesh = CartesianMeshBuilder{TinyVector<3>{0, 0, 0}, TinyVector<3>{2, 1, 1}, TinyVector<3, size_t>{3, 3, 3}}
+                      .mesh()
+                      ->get<Mesh<3>>();
+
+      auto& mesh = *p_mesh;
+
+      for (auto descriptor : descriptor_list) {
+        SECTION(name(descriptor.integrationMethodType()))
+        {
+          SECTION("R^3 data")
+          {
+            auto R3_exact = [](const R3& x) -> R3 {
+              return R3{+2.3 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2),   //
+                        -1.3 * (x[1] - 1) * (x[1] - 1) * (x[1] - 1),   //
+                        +1.4 * (x[2] - 1) * (x[2] - 1) * (x[2] - 1)};
+            };
 
-  //     std::shared_ptr p_mesh = std::make_shared<const Mesh<2>>(initial_mesh.shared_connectivity(), new_xr);
-  //     const Mesh<2>& mesh    = *p_mesh;
-  //     // inverse rotation
-  //     TinyMatrix<2> inv_T{std::cos(theta), std::sin(theta),   //
-  //                         -std::sin(theta), std::cos(theta)};
+            DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
 
-  //     auto p0 = [&inv_T](const R2& X) {
-  //       const R2 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       return +1.7 * x * x + 2 * y * y - 1.1;
-  //     };
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
 
-  //     auto p1 = [&inv_T](const R2& X) {
-  //       const R2 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       return -1.3 * x * x - 0.2 * y * y + 0.7;
-  //     };
-
-  //     auto p2 = [&inv_T](const R2& X) {
-  //       const R2 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       return +2.6 * x * x - 1.4 * y * y - 1.9;
-  //     };
+            auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
 
-  //     auto p3 = [&inv_T](const R2& X) {
-  //       const R2 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       return -0.6 * x * x + 1.4 * y * y + 2.3;
-  //     };
+            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));
+          }
 
-  //     auto q0 = [&inv_T](const R2& X) {
-  //       const R2 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       return 3 * x * y;
-  //     };
-
-  //     auto q1 = [&inv_T](const R2& X) {
-  //       const R2 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       return -1.3 * x * y;
-  //     };
-
-  //     for (auto descriptor : descriptor_list) {
-  //       SECTION(name(descriptor.integrationMethodType()))
-  //       {
-  //         SECTION("R data")
-  //         {
-  //           auto R_exact = p0;
-
-  //           DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
-
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
-
-  //           auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
-
-  //           double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
-  //           REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //         }
-
-  //         SECTION("R2x2 data")
-  //         {
-  //           using R2x2      = TinyMatrix<2>;
-  //           auto R2x2_exact = [&](const R2& X) -> R2x2 {
-  //             return T * TinyMatrix<2>{p0(X), q0(X), q1(X), p1(X)} * inv_T;
-  //           };
-
-  //           DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
-
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
-
-  //           auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>();
-
-  //           double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R2x2_exact));
-  //           REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //         }
-
-  //         SECTION("vector of R")
-  //         {
-  //           std::array<std::function<double(const R2&)>, 4> vector_exact   //
-  //             = {p0, p1, p2, p3};
-
-  //           DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
-
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
-
-  //           auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>();
-
-  //           double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
-  //           REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //         }
-
-  //         SECTION("vector of R2x2")
-  //         {
-  //           using R2x2 = TinyMatrix<2>;
-
-  //           std::array<std::function<R2x2(const R2&)>, 2> vector_R2x2_exact   //
-  //             = {[&](const R2& X) {
-  //                  return T * R2x2{p0(X), q0(X), q1(X), p1(X)} * inv_T;
-  //                },
-  //                [&](const R2& X) {
-  //                  return T * R2x2{p0(X), q1(X), 0, p1(X)} * inv_T;
-  //                }};
-
-  //           DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_R2x2_exact);
-
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
-
-  //           auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2x2>>();
-
-  //           double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_R2x2_exact);
-  //           REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //         }
-
-  //         SECTION("list")
-  //         {
-  //           using R2x2 = TinyMatrix<2>;
-
-  //           auto R_exact = p0;
-
-  //           DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
-
-  //           std::array<std::function<double(const R2&)>, 4> vector_exact   //
-  //             = {p0, p1, p2, p3};
-
-  //           DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
-
-  //           std::array<std::function<R2x2(const R2&)>, 2> vector_R2x2_exact   //
-  //             = {[&](const R2& X) {
-  //                  return T * R2x2{p0(X), q0(X), q1(X), p1(X)} * inv_T;
-  //                },
-  //                [&](const R2& X) {
-  //                  return T * R2x2{p2(X), q1(X), 0, p3(X)} * inv_T;
-  //                }};
-
-  //           DiscreteFunctionP0Vector Wh = test_only::exact_projection(mesh, degree, vector_R2x2_exact);
-
-  //           auto R2x2_exact = [&](const R2& X) -> R2x2 {
-  //             return T * TinyMatrix<2>{p0(X), q1(X), q0(X), p3(X)} * inv_T;
-  //           };
-
-  //           DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
-
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(uh, Vh, Wh, Ah);
-
-  //           auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
-  //           auto dpk_Vh = reconstructions[1]->get<DiscreteFunctionDPkVector<2, const double>>();
-  //           auto dpk_Wh = reconstructions[2]->get<DiscreteFunctionDPkVector<2, const R2x2>>();
-  //           auto dpk_Ah = reconstructions[3]->get<DiscreteFunctionDPk<2, const R2x2>>();
-
-  //           {
-  //             double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
-  //             REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //           }
-  //           {
-  //             double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
-  //             REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //           }
-  //           {
-  //             double max_error = test_only::max_reconstruction_error(mesh, dpk_Wh, vector_R2x2_exact);
-  //             REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //           }
-  //           {
-  //             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("3D")
-  //   {
-  //     using R3 = TinyVector<3>;
-
-  //     auto p_initial_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
-  //     auto& initial_mesh  = *p_initial_mesh;
-
-  //     constexpr double theta = 1;
-  //     TinyMatrix<3> T{std::cos(theta), -std::sin(theta), 0,
-  //                     //
-  //                     std::sin(theta), std::cos(theta), 0,
-  //                     //
-  //                     0, 0, 1};
-
-  //     auto xr = initial_mesh.xr();
-
-  //     NodeValue<R3> new_xr{initial_mesh.connectivity()};
-  //     parallel_for(
-  //       initial_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { new_xr[node_id] = T * xr[node_id]; });
-
-  //     std::shared_ptr p_mesh = std::make_shared<const Mesh<3>>(initial_mesh.shared_connectivity(), new_xr);
-  //     const Mesh<3>& mesh    = *p_mesh;
-  //     // inverse rotation
-  //     TinyMatrix<3> inv_T{std::cos(theta), std::sin(theta), 0,
-  //                         //
-  //                         -std::sin(theta), std::cos(theta), 0,
-  //                         //
-  //                         0, 0, 1};
-
-  //     auto p0 = [&inv_T](const R3& X) {
-  //       const R3 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       const double z = Y[2] - 1;
-  //       return +1.7 * x * x + 2 * y * y + 1.3 * z * z - 1.1;
-  //     };
-
-  //     auto p1 = [&inv_T](const R3& X) {
-  //       const R3 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       const double z = Y[2] - 1;
-  //       return +2.1 * x * x - 1.4 * y * y - 3.1 * z * z + 2.2;
-  //     };
-
-  //     auto p2 = [&inv_T](const R3& X) {
-  //       const R3 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       const double z = Y[2] - 1;
-  //       return -1.1 * x * x - 1.2 * y * y + 1.3 * z * z - 1.7;
-  //     };
-
-  //     auto p3 = [&inv_T](const R3& X) {
-  //       const R3 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       const double z = Y[2] - 1;
-  //       return 1.9 * x * x + 2.1 * y * y - 3.1 * z * z + 1.6;
-  //     };
-
-  //     auto p4 = [&inv_T](const R3& X) {
-  //       const R3 Y     = inv_T * X;
-  //       const double x = Y[0] - 2;
-  //       const double y = Y[1] - 1;
-  //       const double z = Y[2] - 1;
-  //       return -2.4 * x * x + 3.3 * y * y - 1.7 * z * z + 2.1;
-  //     };
-
-  //     SECTION("3 symmetries")
-  //     {
-  //       std::vector<PolynomialReconstructionDescriptor> descriptor_list =
-  //         {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
-  //                                             std::vector<std::shared_ptr<
-  //                                               const IBoundaryDescriptor>>{std::
-  //                                                                             make_shared<NamedBoundaryDescriptor>(
-  //                                                                               "XMAX"),
-  //                                                                           std::make_shared<NamedBoundaryDescriptor>(
-  //                                                                             "YMAX"),
-  //                                                                           std::make_shared<NamedBoundaryDescriptor>(
-  //                                                                             "ZMAX")}},
-  //          PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
-  //                                             std::vector<std::shared_ptr<
-  //                                               const IBoundaryDescriptor>>{std::
-  //                                                                             make_shared<NamedBoundaryDescriptor>(
-  //                                                                               "XMAX"),
-  //                                                                           std::make_shared<NamedBoundaryDescriptor>(
-  //                                                                             "YMAX"),
-  //                                                                           std::make_shared<NamedBoundaryDescriptor>(
-  //                                                                             "ZMAX")}}};
-
-  //       for (auto descriptor : descriptor_list) {
-  //         SECTION("R data")
-  //         {
-  //           auto R_exact = p0;
-
-  //           DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
-
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
-
-  //           auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
-
-  //           double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
-  //           REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //         }
-
-  //         SECTION("vector of R")
-  //         {
-  //           std::array<std::function<double(const R3&)>, 4> vector_exact   //
-  //             = {p0, p1, p2, p3};
-
-  //           DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
-
-  //           auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
-
-  //           auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
-
-  //           double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
-  //           REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //         }
-  //       }
-  //     }
-
-  //     SECTION("1 symmetry")
-  //     {
-  //       // Matrix and their transformations are kept simple to
-  //       // derive exact solutions
-  //       std::vector<PolynomialReconstructionDescriptor> descriptor_list =
-  //         {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
-  //                                             std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
-  //                                               std::make_shared<NamedBoundaryDescriptor>("XMAX")}},
-  //          PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
-  //                                             std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
-  //                                               std::make_shared<NamedBoundaryDescriptor>("XMAX")}}};
-  //       for (auto descriptor : descriptor_list) {
-  //         SECTION(name(descriptor.integrationMethodType()))
-  //         {
-  //           SECTION("R3x3 data")
-  //           {
-  //             // Matrix and their transformations are kept simple to
-  //             // derive exact solutions
-
-  //             using R3x3      = TinyMatrix<3>;
-  //             auto R2x2_exact = [&](const R3& X) -> R3x3 {
-  //               return T * TinyMatrix<3>{p0(X), 0,     0,       //
-  //                                        0,     p1(X), p2(X),   //
-  //                                        0,     p3(X), p4(X)} *
-  //                      inv_T;
-  //             };
-
-  //             DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
-
-  //             auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
-
-  //             auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3x3>>();
-
-  //             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 of R3x3")
-  //           {
-  //             using R3x3 = TinyMatrix<3>;
-
-  //             std::array<std::function<R3x3(const R3&)>, 2> vector_R3x3_exact   //
-  //               = {[&](const R3& X) {
-  //                    return T * R3x3{p0(X), 0,     0,       //
-  //                                    0,     p1(X), p2(X),   //
-  //                                    0,     p3(X), p4(X)} *
-  //                           inv_T;
-  //                  },
-  //                  [&](const R3& X) {
-  //                    return T * R3x3{p0(X), 0,     0,       //
-  //                                    0,     p2(X), p4(X),   //
-  //                                    0,     p3(X), p1(X)} *
-  //                           inv_T;
-  //                  }};
-
-  //             DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_R3x3_exact);
-
-  //             auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
-
-  //             auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3x3>>();
-
-  //             double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_R3x3_exact);
-  //             REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //           }
-
-  //           SECTION("list")
-  //           {
-  //             using R3x3 = TinyMatrix<3>;
-
-  //             auto R_exact = p0;
-
-  //             DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
-
-  //             std::array<std::function<double(const R3&)>, 4> vector_exact   //
-  //               = {p0, p1, p2, p3};
-
-  //             DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
-
-  //             std::array<std::function<R3x3(const R3&)>, 2> vector_R3x3_exact   //
-  //               = {[&](const R3& X) {
-  //                    return T * R3x3{p1(X), 0,     0,       //
-  //                                    0,     p3(X), p0(X),   //
-  //                                    0,     p2(X), p4(X)} *
-  //                           inv_T;
-  //                  },
-  //                  [&](const R3& X) {
-  //                    return T * R3x3{p2(X), 0,     0,       //
-  //                                    0,     p0(X), p1(X),   //
-  //                                    0,     p3(X), p4(X)} *
-  //                           inv_T;
-  //                  }};
-
-  //             DiscreteFunctionP0Vector Wh = test_only::exact_projection(mesh, degree, vector_R3x3_exact);
-
-  //             auto R3x3_exact = [&](const R3& X) -> R3x3 {
-  //               return T * R3x3{p0(X), 0,     0,       //
-  //                               0,     p1(X), p2(X),   //
-  //                               0,     p3(X), p4(X)} *
-  //                      inv_T;
-  //             };
-
-  //             DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R3x3_exact));
-
-  //             auto reconstructions = PolynomialReconstruction{descriptor}.build(uh, Vh, Wh, Ah);
-
-  //             auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
-  //             auto dpk_Vh = reconstructions[1]->get<DiscreteFunctionDPkVector<3, const double>>();
-  //             auto dpk_Wh = reconstructions[2]->get<DiscreteFunctionDPkVector<3, const R3x3>>();
-  //             auto dpk_Ah = reconstructions[3]->get<DiscreteFunctionDPk<3, const R3x3>>();
-
-  //             {
-  //               double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
-  //               REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //             }
-  //             {
-  //               double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
-  //               REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //             }
-  //             {
-  //               double max_error = test_only::max_reconstruction_error(mesh, dpk_Wh, vector_R3x3_exact);
-  //               REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //             }
-  //             {
-  //               double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R3x3_exact));
-  //               REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
-  //             }
-  //           }
-  //         }
-  //       }
-  //     }
-  //   }
-  // }
+          SECTION("vector of R3")
+          {
+            std::array<std::function<R3(const R3&)>, 2> vector_exact   //
+              = {[](const R3& x) -> R3 {
+                   return R3{+1.7 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2),   //
+                             -0.6 * (x[1] - 1) * (x[1] - 1) * (x[1] - 1),   //
+                             +1.2 * (x[2] - 1) * (x[2] - 1) * (x[2] - 1)};
+                 },
+                 [](const R3& x) -> R3 {
+                   return R3{-2.3 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2),   //
+                             +1.1 * (x[1] - 1) * (x[1] - 1) * (x[1] - 1),   //
+                             -0.3 * (x[2] - 1) * (x[2] - 1) * (x[2] - 1)};
+                 }};
+
+            DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+            auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+            auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3>>();
+
+            double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+          }
+        }
+      }
+    }
+  }
 }