#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>

#include <utils/PugsAssert.hpp>

#include <mesh/Mesh.hpp>
#include <mesh/NamedBoundaryDescriptor.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/PolynomialReconstruction.hpp>

#include <DiscreteFunctionDPkForTests.hpp>
#include <MeshDataBaseForTests.hpp>

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

TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
{
  constexpr size_t degree = 1;

  SECTION("without symmetries")
  {
    std::vector<PolynomialReconstructionDescriptor> descriptor_list =
      {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree},
       PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree},
       PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree}};

    for (auto descriptor : descriptor_list) {
      SECTION(name(descriptor.integrationMethodType()))
      {
        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 R_exact = [](const R1& x) { return 2.3 + 1.7 * x[0]; };

                DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));

                auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);

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

                double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact));
                REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
              }
            }
          }

          SECTION("R^3 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 R3_exact = [](const R1& x) -> R3 {
                  return R3{+2.3 + 1.7 * x[0],   //
                            +1.4 - 0.6 * x[0],   //
                            -0.2 + 3.1 * x[0]};
                };

                DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));

                auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);

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

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

          SECTION("R^3x3 data")
          {
            using R3x3 = TinyMatrix<3, 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 R3x3_exact = [](const R1& x) -> R3x3 {
                  return R3x3{
                    +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],
                    -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
                  };
                };

                DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R3x3_exact));

                auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);

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

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

          SECTION("R vector data")
          {
            for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
              SECTION(named_mesh.name())
              {
                auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
                auto& mesh  = *p_mesh;
                std::array<std::function<double(const R1&)>, 3> vector_exact =
                  {[](const R1& x) -> double { return +2.3 + 1.7 * x[0]; },
                   [](const R1& x) -> double { return -1.7 + 2.1 * x[0]; },
                   [](const R1& x) -> double { return +1.4 - 0.6 * x[0]; }};

                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 double>>();

                double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
                REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
              }
            }
          }

          SECTION("R3 vector 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;

                std::array<std::function<R3(const R1&)>, 2> vector_exact =
                  {[](const R1& x) -> R3 {
                     return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
                   },
                   [](const R1& x) -> R3 {
                     return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]};
                   }};

                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 R3>>();

                double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
                REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
              }
            }
          }

          SECTION("list of various types")
          {
            using R3x3 = TinyMatrix<3>;
            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 R_exact = [](const R1& x) { return 2.3 + 1.7 * x[0]; };

                auto R3_exact = [](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 R3x3_exact = [](const R1& x) -> R3x3 {
                  return R3x3{
                    +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],
                    -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
                  };
                };

                std::array<std::function<double(const R1&)>, 3> vector_exact =
                  {[](const R1& x) -> double { return +2.3 + 1.7 * x[0]; },
                   [](const R1& x) -> double { return -1.7 + 2.1 * x[0]; },
                   [](const R1& x) -> double { return +1.4 - 0.6 * x[0]; }};

                DiscreteFunctionP0 fh       = test_only::exact_projection(mesh, degree, std::function(R_exact));
                DiscreteFunctionP0 uh       = test_only::exact_projection(mesh, degree, std::function(R3_exact));
                DiscreteFunctionP0 Ah       = test_only::exact_projection(mesh, degree, std::function(R3x3_exact));
                DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);

                auto reconstructions =
                  PolynomialReconstruction{descriptor}.build(std::make_shared<DiscreteFunctionVariant>(fh), uh,
                                                             std::make_shared<DiscreteFunctionP0<R3x3>>(Ah),
                                                             DiscreteFunctionVariant(Vh));

                {
                  auto dpk_fh      = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
                  double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact));
                  REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
                }

                {
                  auto dpk_uh      = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>();
                  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-14));
                }

                {
                  auto dpk_Ah      = reconstructions[2]->get<DiscreteFunctionDPk<1, const R3x3>>();
                  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-14));
                }

                {
                  auto dpk_Vh      = reconstructions[3]->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-14));
                }
              }
            }
          }
        }

        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 R_exact = [](const R2& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1]; };

                DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));

                auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);

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

                double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact));
                REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
              }
            }
          }

          SECTION("R^3 data")
          {
            using R3 = TinyVector<3>;

            for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
              SECTION(named_mesh.name())
              {
                auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
                auto& mesh  = *p_mesh;

                auto R3_exact = [](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]};
                };

                DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));

                auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);

                auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R3>>();

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

          SECTION("R^2x2 data")
          {
            using R2x2 = TinyMatrix<2, 2>;

            for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
              SECTION(named_mesh.name())
              {
                auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
                auto& mesh  = *p_mesh;

                auto R2x2_exact = [](const R2& x) -> R2x2 {
                  return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1],   //
                              +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
                };

                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<2, 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-14));
              }
            }
          }

          SECTION("vector 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;

                std::array<std::function<double(const R2&)>, 4> vector_exact =
                  {[](const R2& x) -> double { return +2.3 + 1.7 * x[0] + 1.2 * x[1]; },
                   [](const R2& x) -> double { return -1.7 + 2.1 * x[0] - 2.2 * x[1]; },
                   [](const R2& x) -> double { return +1.4 - 0.6 * x[0] - 2.1 * x[1]; },
                   [](const R2& x) -> double { return +2.4 - 2.3 * x[0] + 1.3 * x[1]; }};

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

        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 R_exact = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; };

                DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));

                auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);

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

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

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

                auto R3_exact = [](const R3& x) -> R3 {
                  return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2],   //
                            +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2],   //
                            -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]};
                };

                DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));

                auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);

                auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();

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

          SECTION("R^2x2 data")
          {
            using R2x2 = TinyMatrix<2, 2>;

            for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
              SECTION(named_mesh.name())
              {
                auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
                auto& mesh  = *p_mesh;

                auto R2x2_exact = [](const R3& x) -> R2x2 {
                  return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
                              //
                              +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
                };

                DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));

                descriptor.setRowWeighting(false);
                auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);

                auto dpk_Ah      = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
                double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R2x2_exact));
                REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
              }
            }
          }

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

                std::array<std::function<double(const R3&)>, 4> vector_exact =
                  {[](const R3& x) -> double { return +2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2]; },
                   [](const R3& x) -> double { return -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2]; },
                   [](const R3& x) -> double { return +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2]; },
                   [](const R3& x) -> double { return -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]; }};

                DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);

                descriptor.setPreconditioning(false);
                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("errors")
        {
          auto p_mesh1 = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
          DiscreteFunctionP0<double> f1{p_mesh1};

          auto p_mesh2 = MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
          DiscreteFunctionP0<double> f2{p_mesh2};

          REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(f1, f2),
                              "error: cannot reconstruct functions living of different meshes simultaneously");
        }
      }
    }
  }

  SECTION("with symmetries")
  {
    SECTION("errors")
    {
      SECTION("1D")
      {
        auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();

        PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree,
                                                      std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
                                                        std::make_shared<NamedBoundaryDescriptor>("XMIN")}};

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<2>>{p_mesh}),
                            "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 1");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<3>>{p_mesh}),
                            "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 1");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyVector<2>>{p_mesh, 1}),
                            "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 1");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyVector<3>>{p_mesh, 1}),
                            "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 1");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<2>>{p_mesh}),
                            "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 1");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<3>>{p_mesh}),
                            "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 1");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyMatrix<2>>{p_mesh, 1}),
                            "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 1");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyMatrix<3>>{p_mesh, 1}),
                            "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 1");
      }

      SECTION("2D")
      {
        auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();

        PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree,
                                                      std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
                                                        std::make_shared<NamedBoundaryDescriptor>("XMIN")}};

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<1>>{p_mesh}),
                            "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 2");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<3>>{p_mesh}),
                            "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 2");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyVector<1>>{p_mesh, 1}),
                            "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 2");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyVector<3>>{p_mesh, 1}),
                            "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 2");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<1>>{p_mesh}),
                            "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 2");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<3>>{p_mesh}),
                            "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 2");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyMatrix<1>>{p_mesh, 1}),
                            "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 2");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyMatrix<3>>{p_mesh, 1}),
                            "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 2");
      }

      SECTION("3D")
      {
        auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();

        PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree,
                                                      std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
                                                        std::make_shared<NamedBoundaryDescriptor>("XMIN")}};

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<1>>{p_mesh}),
                            "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 3");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<2>>{p_mesh}),
                            "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 3");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyVector<1>>{p_mesh, 1}),
                            "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 3");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyVector<2>>{p_mesh, 1}),
                            "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 3");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<1>>{p_mesh}),
                            "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 3");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<2>>{p_mesh}),
                            "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 3");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyMatrix<1>>{p_mesh, 1}),
                            "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 3");

        REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(
                              DiscreteFunctionP0Vector<TinyMatrix<2>>{p_mesh, 1}),
                            "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 3");
      }
    }

    SECTION("1D")
    {
      std::vector<PolynomialReconstructionDescriptor> descriptor_list =
        {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree,
                                            std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
                                              std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
         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>;

      for (auto descriptor : descriptor_list) {
        SECTION(name(descriptor.integrationMethodType()))
        {
          SECTION("R^1 data")
          {
            auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();

            auto& mesh = *p_mesh;

            auto R1_exact = [](const R1& x) { return R1{1.7 * (x[0] + 1)}; };

            DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R1_exact));

            auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);

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

            double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R1_exact));
            REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
          }

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

                std::array<std::function<R1(const R1&)>, 2> vector_exact   //
                  = {[](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1)}; },
                     [](const R1& x) -> R1 { return R1{-0.3 * (x[0] + 1)}; }};

                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 R1>>();

                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::cell_center, degree,
                                            std::vector<std::shared_ptr<
                                              const IBoundaryDescriptor>>{std::
                                                                            make_shared<NamedBoundaryDescriptor>(
                                                                              "XMAX"),
                                                                          std::make_shared<NamedBoundaryDescriptor>(
                                                                            "YMAX")}},
         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>;

      for (auto descriptor : descriptor_list) {
        SECTION(name(descriptor.integrationMethodType()))
        {
          SECTION("R^2 data")
          {
            auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
            auto& mesh  = *p_mesh;

            auto R2_exact = [](const R2& x) -> R2 { return R2{2.3 * (x[0] - 2), -1.3 * (x[1] - 1)}; };

            DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R2_exact));

            auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);

            auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2>>();

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

          SECTION("vector of R2")
          {
            auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
            auto& mesh  = *p_mesh;

            std::array<std::function<R2(const R2&)>, 2> vector_exact   //
              = {[](const R2& x) -> R2 {
                   return R2{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1)};
                 },
                 [](const R2& x) -> R2 {
                   return R2{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1)};
                 }};

            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 R2>>();

            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("3D")
    {
      std::vector<PolynomialReconstructionDescriptor> descriptor_list =
        {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 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::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>;

      for (auto descriptor : descriptor_list) {
        SECTION(name(descriptor.integrationMethodType()))
        {
          SECTION("R^3 data")
          {
            auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
            auto& mesh  = *p_mesh;

            auto R3_exact = [](const R3& x) -> R3 { return R3{2.3 * (x[0] - 2), -1.3 * (x[1] - 1), 1.4 * (x[2] - 1)}; };

            DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));

            auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);

            auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();

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

          SECTION("vector of R3")
          {
            auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
            auto& mesh  = *p_mesh;

            std::array<std::function<R3(const R3&)>, 2> vector_exact   //
              = {[](const R3& x) -> R3 {
                   return R3{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1), +1.2 * (x[2] - 1)};
                 },
                 [](const R3& x) -> R3 {
                   return R3{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1), -0.3 * (x[2] - 1)};
                 }};

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