Add affine reconstruction for DiscreteFunctionP0 in 1, 2 and 3d

The possible data types are double, TinyVector and TinyMatrix

src/algebra/Givens.hpp

@@ -40,10 +40,10 @@ class Givens
   static void
   _lift(const MatrixType& A, const RHSVectorType& b, VectorType& x)
   {
-    for (ssize_t i = x.size() - 1; i >= 0; --i) {
+    for (ssize_t i = x.dimension() - 1; i >= 0; --i) {
       const double inv_Aii = 1 / A(i, i);
       double sum           = 0;
-      for (size_t j = i + 1; j < x.size(); ++j) {
+      for (size_t j = i + 1; j < x.dimension(); ++j) {
         sum += A(i, j) * x[j];
       }
       x[i] = (b[i] - sum) * inv_Aii;
@@ -71,12 +71,13 @@ class Givens
   static void
   solveInPlace(MatrixType& A, VectorType& x, RHSVectorType& b)
   {
-    Assert(x.size() == A.numberOfColumns(), "The number of columns of A must be the size of x");
-    Assert(b.size() == A.numberOfRows(), "The number of rows of A must be the size of b");
+    Assert(x.dimension() == A.numberOfColumns(), "The number of columns of A must be the size of x");
+    Assert(b.dimension() == A.numberOfRows(), "The number of rows of A must be the size of b");
 
     for (size_t j = 0; j < A.numberOfColumns(); ++j) {
       for (size_t i = A.numberOfRows() - 1; i > j; --i) {
-        double c(0), s(0);
+        double c;
+        double s;
         _givens(A(i - 1, j), A(i, j), c, s);
         for (size_t k = j; k < A.numberOfColumns(); ++k) {
           double xi_1(0), xi(0);
@@ -98,8 +99,8 @@ class Givens
   static void
   solve(const MatrixType& A, VectorType& x, const RHSVectorType& b)
   {
-    Assert(x.size() == A.numberOfColumns(), "The number of columns of A must be the size of x");
-    Assert(b.size() == A.numberOfRows(), "The number of rows of A must be the size of b");
+    Assert(x.dimension() == A.numberOfColumns(), "The number of columns of A must be the size of x");
+    Assert(b.dimension() == A.numberOfRows(), "The number of rows of A must be the size of b");
 
     MatrixType rotateA    = copy(A);
     RHSVectorType rotateb = copy(b);
@@ -109,7 +110,7 @@ class Givens
 
   template <typename MatrixType, typename UnknownMatrixType, typename RHSMatrixType>
   static void
-  solveCollectionInPlace(const MatrixType& A, UnknownMatrixType& X, const RHSMatrixType& B)
+  solveCollectionInPlace(MatrixType& A, UnknownMatrixType& X, RHSMatrixType& B)
   {
     Assert(X.numberOfRows() == A.numberOfColumns(), "The number of columns of A must be the number of rows of X");
     Assert(B.numberOfRows() == A.numberOfRows(), "The number of rows of A must be the rows of B");
@@ -117,7 +118,8 @@ class Givens
 
     for (size_t j = 0; j < A.numberOfColumns(); ++j) {
       for (size_t i = A.numberOfRows() - 1; i > j; --i) {
-        double c(0), s(0);
+        double c;
+        double s;
         _givens(A(i - 1, j), A(i, j), c, s);
         for (size_t k = j; k < A.numberOfColumns(); ++k) {
           double xi_1(0), xi(0);

src/algebra/SmallVector.hpp

@@ -166,6 +166,13 @@ class SmallVector   // LCOV_EXCL_LINE
     return m_values.size();
   }
 
+  PUGS_INLINE
+  size_t
+  dimension() const noexcept
+  {
+    return m_values.size();
+  }
+
   PUGS_INLINE SmallVector&
   fill(const DataType& value) noexcept
   {

src/scheme/CMakeLists.txt

@@ -10,6 +10,8 @@ add_library(
   DiscreteFunctionVectorIntegrator.cpp
   DiscreteFunctionVectorInterpoler.cpp
   FluxingAdvectionSolver.cpp
+
+  test_reconstruction.cpp
 )
 
 target_link_libraries(

src/scheme/PolynomialReconstruction.hpp

@@ -20,55 +20,149 @@ class PolynomialReconstruction
 {
  public:
   template <MeshConcept MeshType, typename DataType>
-  PUGS_INLINE DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>>
+  [[nodiscard]] PUGS_INLINE DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>>
   build(const MeshType& mesh, const DiscreteFunctionP0<DataType> p0_function)
   {
+    using Rd = TinyVector<MeshType::Dimension>;
+
     const size_t degree   = 1;
     const Stencil stencil = StencilBuilder{}.build(mesh.connectivity());
 
-    auto xr = mesh.xr();
     auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
 
     auto cell_is_owned = mesh.connectivity().cellIsOwned();
 
     DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>> dpk{mesh.meshVariant(), degree};
 
-    if constexpr ((is_polygonal_mesh_v<MeshType>)and(MeshType::Dimension == 1)) {
-      auto qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(1));
-
+    if constexpr (is_polygonal_mesh_v<MeshType>) {
+      if constexpr (std::is_arithmetic_v<DataType>) {
         parallel_for(
           mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
             if (cell_is_owned[cell_j_id]) {
               auto stencil_cell_list = stencil[cell_j_id];
               SmallVector<double> b(stencil_cell_list.size());
 
-            const double p_j = p0_function[cell_j_id];
+              const double pj = p0_function[cell_j_id];
               for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
                 const CellId cell_i_id = stencil_cell_list[i];
 
-              b[i] = p0_function[cell_i_id] - p_j;
+                b[i] = p0_function[cell_i_id] - pj;
+              }
+
+              SmallMatrix<double> A(stencil_cell_list.size(), MeshType::Dimension);
+
+              const Rd& Xj = xj[cell_j_id];
+
+              for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+                const CellId cell_i_id = stencil_cell_list[i];
+                const Rd Xi_Xj         = xj[cell_i_id] - Xj;
+                for (size_t l = 0; l < MeshType::Dimension; ++l) {
+                  A(i, l) = Xi_Xj[l];
+                }
               }
 
-            SmallMatrix<double> A(stencil_cell_list.size(), degree);
+              SmallVector<double> x(MeshType::Dimension);
+              Givens::solveInPlace(A, x, b);
 
-            using R1 = TinyVector<1>;
+              auto dpk_j = dpk.coefficients(cell_j_id);
 
-            const R1& Xj = xj[cell_j_id];
-            auto X_Xj    = [&](const R1& X) { return X - Xj; };
+              dpk_j[0] = pj;
 
+              for (size_t l = 0; l < MeshType::Dimension; ++l) {
+                dpk_j[1 + l] = x[l];
+              }
+            }
+          });
+      } else if constexpr (is_tiny_vector_v<DataType>) {
+        parallel_for(
+          mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
+            if (cell_is_owned[cell_j_id]) {
+              auto stencil_cell_list = stencil[cell_j_id];
+              SmallMatrix<double> B(stencil_cell_list.size(), DataType::Dimension);
+
+              const DataType& pj = p0_function[cell_j_id];
               for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
                 const CellId cell_i_id = stencil_cell_list[i];
+                const DataType& pi_pj  = p0_function[cell_i_id] - pj;
+                for (size_t k = 0; k < DataType::Dimension; ++k) {
+                  B(i, k) = pi_pj[k];
+                }
+              }
+
+              SmallMatrix<double> A(stencil_cell_list.size(), MeshType::Dimension);
 
-              A(i, 0) = X_Xj(xj[cell_i_id])[0];
+              const Rd& Xj = xj[cell_j_id];
+              for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+                const CellId cell_i_id = stencil_cell_list[i];
+                const Rd Xi_Xj         = xj[cell_i_id] - Xj;
+                for (size_t l = 0; l < MeshType::Dimension; ++l) {
+                  A(i, l) = Xi_Xj[l];
+                }
               }
 
-            SmallVector<double> x(1);
-            Givens::solveInPlace(A, x, b);
+              SmallMatrix<double> X(MeshType::Dimension, DataType::Dimension);
+              Givens::solveCollectionInPlace(A, X, B);
+
+              auto dpk_j = dpk.coefficients(cell_j_id);
+
+              dpk_j[0] = pj;
+              for (size_t i = 0; i < MeshType::Dimension; ++i) {
+                auto& dpk_j_ip1 = dpk_j[i + 1];
+                for (size_t k = 0; k < DataType::Dimension; ++k) {
+                  dpk_j_ip1[k] = X(i, k);
+                }
+              }
+            }
+          });
+      } else if constexpr (is_tiny_matrix_v<DataType>) {
+        parallel_for(
+          mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
+            if (cell_is_owned[cell_j_id]) {
+              auto stencil_cell_list = stencil[cell_j_id];
+              SmallMatrix<double> B(stencil_cell_list.size(), DataType::Dimension);
 
-            dpk.coefficients(cell_j_id)[0] = p_j;
-            dpk.coefficients(cell_j_id)[1] = x[0];
+              const DataType& pj = p0_function[cell_j_id];
+              for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+                const CellId cell_i_id = stencil_cell_list[i];
+                const DataType& pi_pj  = p0_function[cell_i_id] - pj;
+                for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+                  for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+                    B(i, k * DataType::NumberOfColumns + l) = pi_pj(k, l);
+                  }
+                }
+              }
+
+              SmallMatrix<double> A(stencil_cell_list.size(), MeshType::Dimension);
+
+              const Rd& Xj = xj[cell_j_id];
+
+              for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+                const CellId cell_i_id = stencil_cell_list[i];
+                const Rd Xi_Xj         = xj[cell_i_id] - Xj;
+                for (size_t l = 0; l < MeshType::Dimension; ++l) {
+                  A(i, l) = Xi_Xj[l];
+                }
+              }
+
+              SmallMatrix<double> X(MeshType::Dimension, DataType::Dimension);
+              Givens::solveCollectionInPlace(A, X, B);
+
+              auto dpk_j = dpk.coefficients(cell_j_id);
+              dpk_j[0]   = pj;
+
+              for (size_t i = 0; i < MeshType::Dimension; ++i) {
+                auto& dpk_j_ip1 = dpk_j[i + 1];
+                for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+                  for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+                    dpk_j_ip1(k, l) = X(i, k * DataType::NumberOfColumns + l);
+                  }
+                }
+              }
             }
           });
+      } else {
+        throw NotImplementedError("dealing with " + demangle<DataType>());
+      }
     }
 
     synchronize(dpk.cellArrays());

src/utils/PugsTraits.hpp

@@ -106,6 +106,12 @@ inline constexpr bool is_tiny_vector_v = false;
 template <size_t N, typename T>
 inline constexpr bool is_tiny_vector_v<TinyVector<N, T>> = true;
 
+template <typename T>
+inline constexpr bool is_tiny_vector_v<const T> = is_tiny_vector_v<std::remove_cvref_t<T>>;
+
+template <typename T>
+inline constexpr bool is_tiny_vector_v<T&> = is_tiny_vector_v<std::remove_cvref_t<T>>;
+
 // Traits is_tiny_matrix
 
 template <typename T>
@@ -114,6 +120,12 @@ inline constexpr bool is_tiny_matrix_v = false;
 template <size_t M, size_t N, typename T>
 inline constexpr bool is_tiny_matrix_v<TinyMatrix<M, N, T>> = true;
 
+template <typename T>
+inline constexpr bool is_tiny_matrix_v<const T> = is_tiny_matrix_v<std::remove_cvref_t<T>>;
+
+template <typename T>
+inline constexpr bool is_tiny_matrix_v<T&> = is_tiny_matrix_v<std::remove_cvref_t<T>>;
+
 // Trais is ItemValue
 
 template <typename T>

tests/test_PolynomialReconstruction.cpp

@@ -24,6 +24,8 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
   {
     using R1 = TinyVector<1>;
 
+    SECTION("R data")
+    {
       for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
         SECTION(named_mesh.name())
         {
@@ -62,4 +64,491 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
         }
       }
     }
+
+    SECTION("R^d data")
+    {
+      using R3 = TinyVector<3>;
+
+      for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+        SECTION(named_mesh.name())
+        {
+          auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+          auto& mesh  = *p_mesh;
+
+          auto affine = [](const R1& x) -> R3 {
+            return R3{+2.3 + 1.7 * x[0],   //
+                      +1.4 - 0.6 * x[0],   //
+                      -0.2 + 3.1 * x[0]};
+          };
+          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0<R3> fh{p_mesh};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+          PolynomialReconstruction reconstruction;
+          auto dpk = reconstruction.build(mesh, fh);
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              max_mean_error = std::max(max_mean_error, l2Norm(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
+          }
+
+          {
+            double max_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R3 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R1{0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R1{0.1}));
+
+              max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+            }
+            REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
+          }
+        }
+      }
+    }
+
+    SECTION("R^mn data")
+    {
+      using R32 = TinyMatrix<3, 2>;
+
+      for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+        SECTION(named_mesh.name())
+        {
+          auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+          auto& mesh  = *p_mesh;
+
+          auto affine = [](const R1& x) -> R32 {
+            return R32{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0],   //
+                       +1.4 - 0.6 * x[0], +2.4 - 2.3 * x[0],   //
+                       -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0]};
+          };
+          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0<R32> fh{p_mesh};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+          PolynomialReconstruction reconstruction;
+          auto dpk = reconstruction.build(mesh, fh);
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              max_mean_error = std::max(max_mean_error, frobeniusNorm(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
+          }
+
+          {
+            double max_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R32 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R1{0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R1{0.1}));
+
+              max_slope_error = std::max(max_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.7, +2.1,   //
+                                                                                                  -0.6, -2.3,   //
+                                                                                                  +3.1, -3.6}));
+            }
+            REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("2D")
+  {
+    using R2 = TinyVector<2>;
+
+    SECTION("R data")
+    {
+      for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+        SECTION(named_mesh.name())
+        {
+          auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+          auto& mesh  = *p_mesh;
+
+          auto affine = [](const R2& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1]; };
+          auto xj     = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0<double> fh{p_mesh};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+          PolynomialReconstruction reconstruction;
+          auto dpk = reconstruction.build(mesh, fh);
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              max_mean_error = std::max(max_mean_error, std::abs(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
+          }
+
+          {
+            double max_x_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const double reconstructed_slope =
+                (dpk[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0.1, 0})) / 0.2;
+
+              max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
+            }
+            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14));
+          }
+
+          {
+            double max_y_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const double reconstructed_slope =
+                (dpk[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0, 0.1})) / 0.2;
+
+              max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
+            }
+            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
+          }
+        }
+      }
+    }
+
+    SECTION("R^d data")
+    {
+      using R3 = TinyVector<3>;
+
+      for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+        SECTION(named_mesh.name())
+        {
+          auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+          auto& mesh  = *p_mesh;
+
+          auto affine = [](const R2& x) -> R3 {
+            return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1],   //
+                      +1.4 - 0.6 * x[0] + 1.3 * x[1],   //
+                      -0.2 + 3.1 * x[0] - 1.1 * x[1]};
+          };
+          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0<R3> fh{p_mesh};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+          PolynomialReconstruction reconstruction;
+          auto dpk = reconstruction.build(mesh, fh);
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              max_mean_error = std::max(max_mean_error, l2Norm(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
+          }
+
+          {
+            double max_x_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R3 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0.1, 0}));
+
+              max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+            }
+            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14));
+          }
+
+          {
+            double max_y_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R3 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0, 0.1}));
+
+              max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
+            }
+            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
+          }
+        }
+      }
+    }
+
+    SECTION("R^mn data")
+    {
+      using R32 = TinyMatrix<3, 2>;
+
+      for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+        SECTION(named_mesh.name())
+        {
+          auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+          auto& mesh  = *p_mesh;
+
+          auto affine = [](const R2& x) -> R32 {
+            return R32{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1],   //
+                       +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1],   //
+                       -0.2 + 3.1 * x[0] + 0.8 * x[1], -3.2 - 3.6 * x[0] - 1.7 * x[1]};
+          };
+          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0<R32> fh{p_mesh};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+          PolynomialReconstruction reconstruction;
+          auto dpk = reconstruction.build(mesh, fh);
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              max_mean_error = std::max(max_mean_error, frobeniusNorm(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
+          }
+
+          {
+            double max_x_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R32 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0.1, 0}));
+
+              max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.7, 2.1,    //
+                                                                                                      -0.6, -2.3,   //
+                                                                                                      +3.1, -3.6}));
+            }
+            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
+          }
+
+          {
+            double max_y_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R32 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0, 0.1}));
+
+              max_y_slope_error = std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.2, -2.2,   //
+                                                                                                      -2.1, 1.3,    //
+                                                                                                      +0.8, -1.7}));
+            }
+            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("3D")
+  {
+    using R3 = TinyVector<3>;
+
+    SECTION("R data")
+    {
+      for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+        SECTION(named_mesh.name())
+        {
+          auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+          auto& mesh  = *p_mesh;
+
+          auto affine = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; };
+          auto xj     = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0<double> fh{p_mesh};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+          PolynomialReconstruction reconstruction;
+          auto dpk = reconstruction.build(mesh, fh);
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              max_mean_error = std::max(max_mean_error, std::abs(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
+          }
+
+          {
+            double max_x_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const double reconstructed_slope =
+                (dpk[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk[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(max_x_slope_error == Catch::Approx(0).margin(1E-14));
+          }
+
+          {
+            double max_y_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const double reconstructed_slope =
+                (dpk[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk[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(max_y_slope_error == Catch::Approx(0).margin(1E-14));
+          }
+
+          {
+            double max_z_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const double reconstructed_slope =
+                (dpk[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk[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(max_z_slope_error == Catch::Approx(0).margin(1E-14));
+          }
+        }
+      }
+    }
+
+    SECTION("R^d 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 affine = [](const R3& x) -> R4 {
+            return R4{+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],   //
+                      -1.2 - 1.5 * x[0] - 0.1 * x[1] - 1.6 * x[2]};
+          };
+          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0<R4> fh{p_mesh};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+          PolynomialReconstruction reconstruction;
+          auto dpk = reconstruction.build(mesh, fh);
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              max_mean_error = std::max(max_mean_error, l2Norm(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
+          }
+
+          {
+            double max_x_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R4 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+
+              max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R4{1.7, -0.6, 3.1, -1.5}));
+            }
+            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-13));
+          }
+
+          {
+            double max_y_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R4 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+
+              max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R4{-2.2, 1.3, -1.1, -0.1}));
+            }
+            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-13));
+          }
+
+          {
+            double max_z_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R4 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+
+              max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R4{1.8, -3.7, 1.9, -1.6}));
+            }
+            REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-13));
+          }
+        }
+      }
+    }
+
+    SECTION("R^mn data")
+    {
+      using R32 = TinyMatrix<3, 2>;
+
+      for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+        SECTION(named_mesh.name())
+        {
+          auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+          auto& mesh  = *p_mesh;
+
+          auto affine = [](const R3& x) -> R32 {
+            return R32{+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],   //
+                       +1.4 - 0.6 * x[0] - 2.1 * x[1] + 2.9 * 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], -3.2 - 3.6 * x[0] - 1.7 * x[1] + 0.7 * x[2]};
+          };
+          auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+          DiscreteFunctionP0<R32> fh{p_mesh};
+
+          parallel_for(
+            mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+          PolynomialReconstruction reconstruction;
+          auto dpk = reconstruction.build(mesh, fh);
+
+          {
+            double max_mean_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              max_mean_error = std::max(max_mean_error, frobeniusNorm(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+            }
+            REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
+          }
+
+          {
+            double max_x_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R32 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+
+              max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.7, 2.1,    //
+                                                                                                      -0.6, -2.3,   //
+                                                                                                      +3.1, -3.6}));
+            }
+            REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
+          }
+
+          {
+            double max_y_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R32 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+
+              max_y_slope_error = std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.2, -2.2,   //
+                                                                                                      -2.1, 1.3,    //
+                                                                                                      +0.8, -1.7}));
+            }
+            REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
+          }
+
+          {
+            double max_z_slope_error = 0;
+            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+              const R32 reconstructed_slope =
+                (1 / 0.2) * (dpk[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+
+              max_z_slope_error = std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R32{-1.3, -2.4,   //
+                                                                                                      +2.9, +1.4,   //
+                                                                                                      -1.8, +0.7}));
+            }
+            REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12));
+          }
+        }
+      }
+    }
+  }
 }