diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index d7ba065480be4aac1949b813df03fb34d94d3f56..e3739abe7afd3a625b6510d9be792bdb3b3a76dd 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -280,6 +280,7 @@ add_executable (mpi_unit_tests
test_ParallelChecker_read.cpp
test_Partitioner.cpp
test_PolynomialReconstruction_degree_1.cpp
+ test_PolynomialReconstruction_degree_2.cpp
test_PolynomialReconstructionDescriptor.cpp
test_RandomEngine.cpp
test_StencilBuilder_cell2cell.cpp
diff --git a/tests/test_PolynomialReconstruction_degree_2.cpp b/tests/test_PolynomialReconstruction_degree_2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..350a57426f4c70629e39edd3c2edb3754e350cec
--- /dev/null
+++ b/tests/test_PolynomialReconstruction_degree_2.cpp
@@ -0,0 +1,713 @@
+#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_2", "[scheme]")
+{
+ constexpr size_t degree = 2;
+
+ SECTION("without symmetries")
+ {
+ std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ {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] - 1.4 * x[0] * 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] - x[0] * x[0], //
+ +1.4 - 0.6 * x[0] + 2 * x[0] * x[0], //
+ -0.2 + 3.1 * x[0] + 1.4 * x[0] * 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] - 2.3 * x[0] * x[0], //
+ -1.7 + 2.1 * x[0] + 1.2 * x[0] * x[0], //
+ +1.4 - 0.6 * x[0] - 2.0 * x[0] * x[0], //
+ //
+ +2.4 - 2.3 * x[0] + 1.1 * x[0] * x[0], //
+ -0.2 + 3.1 * x[0] - 0.7 * x[0] * x[0], //
+ -3.2 - 3.6 * x[0] + 0.1 * x[0] * x[0], //
+ //
+ -4.1 + 3.1 * x[0] - 0.2 * x[0] * x[0], //
+ +0.8 + 2.9 * x[0] + 4.1 * x[0] * x[0], //
+ -1.6 + 2.3 * x[0] - 1.7 * x[0] * 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] + 1.2 * x[0] * x[0]; },
+ [](const R1& x) -> double { return -1.7 + 2.1 * x[0] + 2.1 * x[0] * x[0]; },
+ [](const R1& x) -> double { return +1.4 - 0.6 * x[0] - 1.3 * x[0] * 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] + 0.8 * x[0] * x[0], //
+ -1.7 + 2.1 * x[0] - 0.7 * x[0] * x[0], //
+ +1.4 - 0.6 * x[0] + 1.9 * x[0] * 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] + 2.3 * x[0] * x[0], //
+ +1.4 - 0.6 * x[0] - 1.5 * x[0] * x[0], //
+ -0.2 + 3.1 * x[0] + 2.9 * x[0] * x[0]};
+ };
+
+ auto R3x3_exact = [](const R1& x) -> R3x3 {
+ return R3x3{
+ +2.3 + 1.7 * x[0] - 0.3 * x[0] * x[0],
+ -1.7 + 2.1 * x[0] + 1.4 * x[0] * x[0],
+ +1.4 - 0.6 * x[0] - 2.3 * x[0] * x[0],
+ //
+ +2.4 - 2.3 * x[0] + 1.8 * x[0] * x[0],
+ -0.2 + 3.1 * x[0] - 1.7 * x[0] * x[0],
+ -3.2 - 3.6 * x[0] + 0.7 * x[0] * x[0],
+ //
+ -4.1 + 3.1 * x[0] - 1.9 * x[0] * x[0],
+ +0.8 + 2.9 * x[0] + 2.2 * x[0] * x[0],
+ -1.6 + 2.3 * x[0] - 1.3 * x[0] * x[0],
+ };
+ };
+
+ std::array<std::function<double(const R1&)>, 3> vector_exact =
+ {[](const R1& x) -> double { return +2.3 + 1.7 * x[0] + 2.2 * x[0] * x[0]; },
+ [](const R1& x) -> double { return -1.7 + 2.1 * x[0] - 1.9 * x[0] * x[0]; },
+ [](const R1& x) -> double { return +1.4 - 0.6 * x[0] + 3.1 * x[0] * 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] //
+ + 1.2 * x[0] * x[0] + 1.3 * x[0] * x[1] - 3.2 * x[1] * 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-12));
+ }
+ }
+ }
+
+ 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] //
+ - 2.1 * x[0] * x[0] - 2.3 * x[0] * x[1] - 3.2 * x[1] * x[1], //
+ +1.4 - 0.6 * x[0] + 1.3 * x[1] //
+ + 2.3 * x[0] * x[0] - 1.3 * x[0] * x[1] + 1.2 * x[1] * x[1], //
+ -0.2 + 3.1 * x[0] - 1.1 * x[1] //
+ - 2.1 * x[0] * x[0] + 1.3 * x[0] * x[1] - 1.1 * x[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-12));
+ }
+ }
+ }
+
+ 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] //
+ - 2.1 * x[0] * x[0] + 1.3 * x[0] * x[1] + 1.2 * x[1] * x[1], //
+ -1.7 + 2.1 * x[0] - 2.2 * x[1] //
+ - 1.2 * x[0] * x[0] + 2.1 * x[0] * x[1] - 1.3 * x[1] * x[1], //
+ +1.4 - 0.6 * x[0] - 2.1 * x[1] //
+ - 1.1 * x[0] * x[0] - 2.3 * x[0] * x[1] + 2.1 * x[1] * x[1],
+ +2.4 - 2.3 * x[0] + 1.3 * x[1] //
+ + 2.7 * x[0] * x[0] + 2.1 * x[0] * x[1] - 2.7 * x[1] * 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-12));
+ }
+ }
+ }
+
+ 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] + 1.2 * x[0] * x[0] - 2.1 * x[0] * x[1] + 3.1 * x[1] * x[1];
+ },
+ [](const R2& x) -> double {
+ return -1.7 + 2.1 * x[0] - 2.2 * x[1] - 0.7 * x[0] * x[0] + 2.2 * x[0] * x[1] - 1.6 * x[1] * x[1];
+ },
+ [](const R2& x) -> double {
+ return +1.4 - 0.6 * x[0] - 2.1 * x[1] + 2.3 * x[0] * x[0] + 2.3 * x[0] * x[1] - 2.9 * x[1] * x[1];
+ },
+ [](const R2& x) -> double {
+ return +2.4 - 2.3 * x[0] + 1.3 * x[1] - 2.7 * x[0] * x[0] - 1.2 * x[0] * x[1] - 0.7 * x[1] * 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-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 R_exact = [](const R3& x) {
+ return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2] //
+ + 1.7 * x[0] * x[0] + 1.4 * x[1] * x[1] + 1.7 * x[2] * x[2] //
+ - 2.3 * x[0] * x[1] + 1.6 * x[0] * x[2] - 1.9 * x[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.7 * x[0] * x[0] - 2.4 * x[1] * x[1] - 2.3 * x[2] * x[2] //
+ - 2.1 * x[0] * x[1] + 2.6 * x[0] * x[2] + 1.6 * x[1] * x[2], //
+ +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2] //
+ + 3.1 * x[0] * x[0] - 1.1 * x[1] * x[1] + 1.7 * x[2] * x[2] //
+ - 2.3 * x[0] * x[1] - 2.6 * x[0] * x[2] - 1.9 * x[1] * x[2], //
+ -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2] //
+ - 1.5 * x[0] * x[0] + 1.4 * x[1] * x[1] - 1.2 * x[2] * x[2] //
+ - 1.7 * x[0] * x[1] - 1.3 * x[0] * x[2] + 2.1 * x[1] * x[2]};
+ };
+
+ 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 * x[0] * x[0] + 1.4 * x[1] * x[1] + 1.7 * x[2] * x[2] //
+ - 1.3 * x[0] * x[1] + 1.6 * x[0] * x[2] - 1.9 * x[1] * x[2], //
+ -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2] //
+ + 3.7 * x[0] * x[0] + 1.3 * x[1] * x[1] + 1.6 * x[2] * x[2] //
+ - 2.1 * x[0] * x[1] - 1.5 * x[0] * x[2] - 1.7 * x[1] * x[2], //
+ //
+ +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2] //
+ - 2.1 * x[0] * x[0] + 1.7 * x[1] * x[1] + 1.8 * x[2] * x[2] //
+ - 1.4 * x[0] * x[1] + 1.3 * x[0] * x[2] - 2.9 * x[1] * x[2], //
+ -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2] //
+ + 1.6 * x[0] * x[0] + 2.1 * x[1] * x[1] - 2.1 * x[2] * x[2] //
+ - 1.1 * x[0] * x[1] - 1.3 * x[0] * x[2] + 1.6 * x[1] * 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;
+#warning continue here
+ 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("with symmetries")
+ {
+ // SECTION("1D")
+ // {
+ // std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ // {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
+ // std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+ // std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
+ // PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
+ // std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+ // std::make_shared<NamedBoundaryDescriptor>("XMIN")}}};
+
+ // using R1 = TinyVector<1>;
+
+ // 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::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::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));
+ // }
+ // }
+ // }
+ // }
+ }
+}