Skip to content
Snippets Groups Projects
Select Git revision
  • 346c39d4a5f043e6330e33301d6173e23671b3a7
  • develop default protected
  • feature/advection
  • feature/composite-scheme-other-fluxes
  • origin/stage/bouguettaia
  • save_clemence
  • feature/local-dt-fsi
  • feature/variational-hydro
  • feature/gmsh-reader
  • feature/reconstruction
  • feature/kinetic-schemes
  • feature/composite-scheme-sources
  • feature/serraille
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

test_AffectationProcessor.cpp

Blame
    • test_PolynomialReconstruction.cpp 23.02 KiB 
    #include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>

#include <Kokkos_Core.hpp>

#include <utils/PugsAssert.hpp>
#include <utils/Types.hpp>

#include <algebra/SmallMatrix.hpp>
#include <algebra/SmallVector.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshDataManager.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/PolynomialReconstruction.hpp>

#include <MeshDataBaseForTests.hpp>

template <typename... DiscreteFunctionT>
std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
build_list(DiscreteFunctionT... input)
{
  std::vector<std::shared_ptr<const DiscreteFunctionVariant>> variant_vector;
  auto convert = [&variant_vector](auto&& df) {
    using DF_T = std::decay_t<decltype(df)>;
    if constexpr (is_discrete_function_v<DF_T> or std::is_same_v<DiscreteFunctionVariant, DF_T>) {
      variant_vector.push_back(std::make_shared<DiscreteFunctionVariant>(df));
    } else if constexpr (is_shared_ptr_v<DF_T>) {
      using DF_Value_T = std::decay_t<typename DF_T::element_type>;
      if constexpr (is_discrete_function_v<DF_Value_T> or std::is_same_v<DiscreteFunctionVariant, DF_Value_T>) {
        variant_vector.push_back(std::make_shared<DiscreteFunctionVariant>(*df));
      } else {
        static_assert(is_false_v<DF_T>, "unexpected type");
      }
    } else {
      static_assert(is_false_v<DF_T>, "unexpected type");
    }
  };

  (convert(input), ...);
  return variant_vector;
}

// clazy:excludeall=non-pod-global-static

TEST_CASE("PolynomialReconstruction", "[scheme]")
{
  SECTION("1D")
  {
    using R1 = TinyVector<1>;

    SECTION("R data")
    {
      for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
        SECTION(named_mesh.name())
        {
          auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
          auto& mesh  = *p_mesh;

          auto affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
          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]); });

          DiscreteFunctionVariant fh_v(fh);
          std::shared_ptr<const DiscreteFunctionVariant> p_fh_v = std::make_shared<DiscreteFunctionVariant>(fh);
    
          auto reconstructions = PolynomialReconstruction{}.build(build_list(p_fh_v, fh_v, fh));

          auto dpk = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();

          {
            double max_mean_error = 0;
            for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
              max_mean_error = std::max(max_mean_error, std::abs(dpk[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 double reconstructed_slope =
                (dpk[cell_id](R1{0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R1{0.1})) / 0.2;

              max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
            }
            REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
          }
        }
      }
    }

    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));
          }
        }
      }
    }
  }
}