diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 77fed256505a4862382e9e05691433a40ec3fa7e..386d741d6f1b73028e50e2d71d87a22cb4183764 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -1384,6 +1384,22 @@ PolynomialReconstruction::_build(
}
} else {
Givens::solveCollectionInPlace(A, X, B);
+#warning REMOVE
+ if (cell_j_id == 0) {
+ for (size_t l = 0; l < B.numberOfRows(); ++l) {
+ for (size_t i = 0; i < B.numberOfColumns(); ++i) {
+ std::clog << "B(" << l << ", " << i << ") =" << B(l, i) << '\n';
+ ;
+ }
+ }
+
+ for (size_t l = 0; l < X.numberOfRows(); ++l) {
+ for (size_t i = 0; i < X.numberOfColumns(); ++i) {
+ std::clog << "X(" << l << ", " << i << ") =" << X(l, i) << '\n';
+ ;
+ }
+ }
+ }
}
column_begin = 0;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index e3739abe7afd3a625b6510d9be792bdb3b3a76dd..fd8915b83320770bc36e8a83ad002db44b71bfd5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -281,6 +281,7 @@ add_executable (mpi_unit_tests
test_Partitioner.cpp
test_PolynomialReconstruction_degree_1.cpp
test_PolynomialReconstruction_degree_2.cpp
+ test_PolynomialReconstruction_degree_3.cpp
test_PolynomialReconstructionDescriptor.cpp
test_RandomEngine.cpp
test_StencilBuilder_cell2cell.cpp
diff --git a/tests/test_PolynomialReconstruction_degree_3.cpp b/tests/test_PolynomialReconstruction_degree_3.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d04ac14b7282e33e9b35ed2dc7fee18ed9aafb54
--- /dev/null
+++ b/tests/test_PolynomialReconstruction_degree_3.cpp
@@ -0,0 +1,1230 @@
+#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 <mesh/CartesianMeshBuilder.hpp>
+
+#include <DiscreteFunctionDPkForTests.hpp>
+#include <MeshDataBaseForTests.hpp>
+
+#include <NbGhostLayersTester.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PolynomialReconstruction_degree_3", "[scheme]")
+{
+ constexpr size_t degree = 3;
+
+ constexpr size_t nb_ghost_layers = 3;
+ NbGhostLayersTester t{nb_ghost_layers};
+
+ SECTION("without symmetries")
+ {
+ std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree}};
+
+ for (auto descriptor : descriptor_list) {
+#warning remove
+ descriptor.setPreconditioning(false);
+ descriptor.setRowWeighting(false);
+
+ SECTION(name(descriptor.integrationMethodType()))
+ {
+ SECTION("1D")
+ {
+ using R1 = TinyVector<1>;
+
+ auto p0 = [](const R1& X) {
+ const double x = X[0];
+ return +2.3 + (1.4 + (1.7 - 2.3 * x) * x) * x;
+ };
+ auto p1 = [](const R1& X) {
+ const double x = X[0];
+ return -1.2 - (2.3 - (1.1 + 2.1 * x) * x) * x;
+ };
+ auto p2 = [](const R1& X) {
+ const double x = X[0];
+ return +2.1 + (2.1 + (3.1 - 1.7 * x) * x) * x;
+ };
+ auto p3 = [](const R1& X) {
+ const double x = X[0];
+ return -1.7 + (1.4 + (1.6 - 3.1 * x) * x) * x;
+ };
+ auto p4 = [](const R1& X) {
+ const double x = X[0];
+ return +1.1 - (1.3 - (2.1 + 1.5 * x) * x) * x;
+ };
+ auto p5 = [](const R1& X) {
+ const double x = X[0];
+ return +1.9 - (2.1 + (1.6 - 2.7 * x) * x) * x;
+ };
+ auto p6 = [](const R1& X) {
+ const double x = X[0];
+ return -0.7 + (1.4 + (2.1 + 1.1 * x) * x) * x;
+ };
+ auto p7 = [](const R1& X) {
+ const double x = X[0];
+ return -1.4 - (1.2 + (1.5 - 2.1 * x) * x) * x;
+ };
+ auto p8 = [](const R1& X) {
+ const double x = X[0];
+ return -2.1 + (1.1 - (1.7 + 1.2 * x) * x) * x;
+ };
+ auto p9 = [](const R1& X) {
+ const double x = X[0];
+ return +1.8 - (3.1 + (2.1 - 2.4 * x) * x) * x;
+ };
+
+ 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 = p0;
+
+ 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-12));
+ }
+ }
+ }
+
+ 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{p2(x), p4(x), p1(x)}; };
+
+ 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-12));
+ }
+ }
+ }
+
+ 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{p1(x), p2(x), p3(x), //
+ p4(x), p5(x), p6(x), //
+ p7(x), p8(x), p9(x)};
+ };
+
+ 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-12));
+ }
+ }
+ }
+
+ 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 = {p1, p7, p9};
+
+ 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-12));
+ }
+ }
+ }
+
+ 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&)>, 3> vector_exact //
+ = {[&](const R1& x) -> R3 {
+ return R3{p1(x), p2(x), p3(x)};
+ },
+ [&](const R1& x) -> R3 {
+ return R3{p5(x), p7(x), p0(x)};
+ },
+ [&](const R1& x) -> R3 {
+ return R3{p9(x), p8(x), p4(x)};
+ }};
+
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const 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));
+ }
+ }
+ }
+
+ 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 = p0;
+
+ auto R3_exact = [&](const R1& x) -> R3 { return R3{p9(x), p4(x), p7(x)}; };
+
+ auto R3x3_exact = [&](const R1& x) -> R3x3 {
+ return R3x3{p2(x), p1(x), p0(x), //
+ p3(x), p2(x), p4(x), //
+ p6(x), p5(x), p9(x)};
+ };
+
+ std::array<std::function<double(const R1&)>, 3> vector_exact = {p1, p8, p7};
+
+ 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-12));
+ }
+
+ {
+ 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-12));
+ }
+
+ {
+ 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-12));
+ }
+
+ {
+ 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-12));
+ }
+ }
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
+ auto p0 = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double y2 = y * y;
+ const double x3 = x2 * x;
+ const double x2y = x2 * y;
+ const double xy2 = x * y2;
+ const double y3 = y * y2;
+
+ return +2.3 //
+ + 1.7 * x - 1.3 * y //
+ + 1.2 * x2 + 1.3 * xy - 3.2 * y2 //
+ - 1.3 * x3 + 2.1 * x2y - 1.6 * xy2 + 2.1 * y3;
+ };
+
+ auto p1 = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double y2 = y * y;
+ const double x3 = x2 * x;
+ const double x2y = x2 * y;
+ const double xy2 = x * y2;
+ const double y3 = y * y2;
+
+ return +1.4 //
+ + 1.6 * x - 2.1 * y //
+ + 1.3 * x2 + 2.6 * xy - 1.4 * y2 //
+ - 1.2 * x3 - 1.7 * x2y + 2.1 * xy2 - 2.2 * y3;
+ };
+
+ auto p2 = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double y2 = y * y;
+ const double x3 = x2 * x;
+ const double x2y = x2 * y;
+ const double xy2 = x * y2;
+ const double y3 = y * y2;
+
+ return -1.2 //
+ + 2.3 * x - 1.6 * y //
+ - 1.2 * x2 + 2.4 * xy - 1.9 * y2 //
+ + 0.9 * x3 + 1.5 * x2y - 2.3 * xy2 - 1.6 * y3;
+ };
+
+ auto p3 = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double y2 = y * y;
+ const double x3 = x2 * x;
+ const double x2y = x2 * y;
+ const double xy2 = x * y2;
+ const double y3 = y * y2;
+
+ return +2.4 //
+ + 2.5 * x + 1.4 * y //
+ - 2.7 * x2 + 1.9 * xy - 2.2 * y2 //
+ - 1.3 * x3 + 2.3 * x2y - 1.4 * xy2 + 2.2 * y3;
+ };
+
+ 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 = p0;
+
+ DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree + 1, 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{p1(x), p2(x), p3(x)}; };
+
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree + 1, 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{p0(x), p1(x), //
+ p2(x), p3(x)};
+ };
+
+ DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree + 1, 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 = {p0, p1, p2, p3};
+
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree + 1, 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>;
+
+ auto p0 = [](const R3& X) -> double {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ const double xy = x * y;
+ const double xz = x * z;
+ const double yz = y * z;
+ const double x2 = x * x;
+ const double y2 = y * y;
+ const double z2 = z * z;
+ const double xyz = x * y * z;
+ const double x2y = x2 * y;
+ const double x2z = x2 * z;
+ const double xy2 = x * y2;
+ const double y2z = y2 * z;
+ const double xz2 = x * z2;
+ const double yz2 = y * z2;
+ const double x3 = x2 * x;
+ const double y3 = y2 * y;
+ const double z3 = z2 * z;
+
+ return (y - 0.5) * (z - 0.5) * (z - 0.5);
+
+ // return 2.3 + 1.7 * x - 1.3 * y + 2.1 * z //
+ // + 1.7 * x2 + 1.4 * y2 + 1.7 * z2 //
+ // - 2.3 * xy + 1.6 * xz - 1.9 * yz //
+ // + 1.2 * x3 - 2.1 * y3 - 1.1 * z3 //
+ // - 1.7 * x2y - 1.3 * x2z + 0.9 * xy2 - 0.7 * y2z + 1.5 * xz2 //- 0.7 * yz2
+ // + 2.8 * xyz //
+ // ;
+ };
+
+ // auto p1 = [](const R3& x) {
+ // return +2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2] //
+ // + 1.7 * x[0] * x[0] - 2.4 * x[1] * x[1] - 2.3 * x[2] * x[2] //
+ // - 2.1 * x[0] * x[1] + 2.6 * x[0] * x[2] + 1.6 * x[1] * x[2];
+ // };
+
+ // auto p2 = [](const R3& x) {
+ // return +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2] //
+ // + 3.1 * x[0] * x[0] - 1.1 * x[1] * x[1] + 1.7 * x[2] * x[2] //
+ // - 2.3 * x[0] * x[1] - 2.6 * x[0] * x[2] - 1.9 * x[1] * x[2];
+ // };
+
+ // auto p3 = [](const R3& x) {
+ // return -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2] //
+ // - 1.5 * x[0] * x[0] + 1.4 * x[1] * x[1] - 1.2 * x[2] * x[2] //
+ // - 1.7 * x[0] * x[1] - 1.3 * x[0] * x[2] + 2.1 * x[1] * x[2];
+ // };
+
+ SECTION("R data")
+ {
+ for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_initial_mesh =
+ CartesianMeshBuilder{TinyVector<3>{0, 0, 0}, TinyVector<3>{4, 4, 4}, TinyVector<3, size_t>{4, 4, 4}}
+ .mesh()
+ ->get<Mesh<3>>(); // named_mesh.mesh()->get<Mesh<3>>();
+ auto& initial_mesh = *p_initial_mesh;
+
+ constexpr double theta = 0;
+ TinyMatrix<3> T0{std::cos(theta), -std::sin(theta), 0,
+ //
+ std::sin(theta), std::cos(theta), 0,
+ //
+ 0, 0, 1};
+ constexpr double phi = 0;
+ TinyMatrix<3> T1{std::cos(phi), 0, -std::sin(phi),
+ //
+ 0, 1, 0,
+ //
+ std::sin(phi), 0, std::cos(phi)};
+
+ TinyMatrix T = T0 * T1;
+
+ auto xr = initial_mesh.xr();
+
+ NodeValue<R3> new_xr{initial_mesh.connectivity()};
+ parallel_for(
+ initial_mesh.numberOfNodes(),
+ PUGS_LAMBDA(const NodeId node_id) { new_xr[node_id] = T * xr[node_id]; });
+
+ std::shared_ptr p_mesh = std::make_shared<const Mesh<3>>(initial_mesh.shared_connectivity(), new_xr);
+ const Mesh<3>& mesh = *p_mesh;
+ // inverse rotation
+ TinyMatrix<3> inv_T{std::cos(theta), std::sin(theta), 0,
+ //
+ -std::sin(theta), std::cos(theta), 0,
+ //
+ 0, 0, 1};
+
+ auto R_exact = p0;
+
+ DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree + 3, std::function(R_exact));
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
+
+ {
+ // PolynomialReconstructionDescriptor descriptor2{IntegrationMethodType::element, 2};
+ // DiscreteFunctionP0 f2h = test_only::exact_projection(mesh, 2, std::function(R_exact));
+ // auto reconstructions2 = PolynomialReconstruction{descriptor2}.build(fh);
+ // auto dpk_fh2 = reconstructions2[0]->get<DiscreteFunctionDPk<3, const double>>();
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_fh2,
+ // std::function(R_exact)); REQUIRE(parallel::allReduceMax(max_error) ==
+ // Catch::Approx(0).margin(1E-12));
+
+ auto cell_type = mesh.connectivity().cellType();
+ auto cell_number = mesh.connectivity().cellNumber();
+
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ std::clog << cell_id << "(" << cell_number[cell_id] << ")"
+ << ":---- " << name(cell_type[cell_id]) << "\n2| ";
+ // {
+ // auto coeffs = dpk_fh2.coefficients(cell_id);
+ // for (size_t i = 0; i < coeffs.size(); ++i) {
+ // if (std::abs(coeffs[i]) > 1e-12) {
+ // std::clog << ' ' << i << ':' << coeffs[i];
+ // }
+ // }
+ // }
+ std::clog << '\n';
+ std::clog << "3| ";
+ {
+ auto coeffs = dpk_fh.coefficients(cell_id);
+ for (size_t i = 0; i < coeffs.size(); ++i) {
+ if (std::abs(coeffs[i]) > 1e-12) {
+ std::clog << ' ' << i << ':' << coeffs[i];
+ }
+ }
+ }
+ std::clog << '\n';
+ break;
+ }
+ }
+
+ 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{p1(x), p2(x), p3(x)}; };
+
+ // 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{p0(x), p1(x), //
+ // p2(x), p3(x)};
+ // };
+
+ // 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 = {p0, p1, p2, p3};
+
+ // 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>;
+
+ // auto p0 = [](const R1& x) { return +1.7 * (x[0] + 1) * (x[0] + 1) - 1.1; };
+ // auto p1 = [](const R1& x) { return -1.2 * (x[0] + 1) * (x[0] + 1) + 1.3; };
+ // auto p2 = [](const R1& x) { return +1.4 * (x[0] + 1) * (x[0] + 1) - 0.6; };
+
+ // for (auto descriptor : descriptor_list) {
+ // SECTION(name(descriptor.integrationMethodType()))
+ // {
+ // SECTION("R data")
+ // {
+ // auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ // auto& mesh = *p_mesh;
+
+ // auto R_exact = p0;
+
+ // 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-12));
+ // }
+
+ // SECTION("R1x1 data")
+ // {
+ // using R1x1 = TinyMatrix<1>;
+
+ // auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ // auto& mesh = *p_mesh;
+
+ // auto R1x1_exact = [&](const R1& x) { return R1x1{p0(x)}; };
+
+ // DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R1x1_exact));
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+
+ // auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1x1>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R1x1_exact));
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // SECTION("R vector data")
+ // {
+ // auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+ // auto& mesh = *p_mesh;
+
+ // std::array<std::function<double(const R1&)>, 3> vector_exact //
+ // = {p0, p1, p2};
+
+ // 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-12));
+ // }
+
+ // SECTION("R1x1 vector data")
+ // {
+ // using R1x1 = TinyMatrix<1>;
+
+ // auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+ // auto& mesh = *p_mesh;
+
+ // std::array<std::function<R1x1(const R1&)>, 3> vector_exact //
+ // = {[&](const R1& x) { return R1x1{p2(x)}; }, //
+ // [&](const R1& x) { return R1x1{p0(x)}; }, //
+ // [&](const R1& x) { return R1x1{p1(x)}; }};
+
+ // DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ // auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R1x1>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // }
+ // }
+ // }
+
+ // SECTION("2D")
+ // {
+ // std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ // {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
+ // std::vector<std::shared_ptr<
+ // const IBoundaryDescriptor>>{std::
+ // make_shared<NamedBoundaryDescriptor>(
+ // "XMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "YMAX")}},
+ // PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
+ // std::vector<std::shared_ptr<
+ // const IBoundaryDescriptor>>{std::
+ // make_shared<NamedBoundaryDescriptor>(
+ // "XMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "YMAX")}}};
+
+ // using R2 = TinyVector<2>;
+
+ // auto p_initial_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+ // auto& initial_mesh = *p_initial_mesh;
+
+ // constexpr double theta = 1;
+ // TinyMatrix<2> T{std::cos(theta), -std::sin(theta), //
+ // std::sin(theta), std::cos(theta)};
+
+ // auto xr = initial_mesh.xr();
+
+ // NodeValue<R2> new_xr{initial_mesh.connectivity()};
+ // parallel_for(
+ // initial_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { new_xr[node_id] = T * xr[node_id]; });
+
+ // std::shared_ptr p_mesh = std::make_shared<const Mesh<2>>(initial_mesh.shared_connectivity(), new_xr);
+ // const Mesh<2>& mesh = *p_mesh;
+ // // inverse rotation
+ // TinyMatrix<2> inv_T{std::cos(theta), std::sin(theta), //
+ // -std::sin(theta), std::cos(theta)};
+
+ // auto p0 = [&inv_T](const R2& X) {
+ // const R2 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // return +1.7 * x * x + 2 * y * y - 1.1;
+ // };
+
+ // auto p1 = [&inv_T](const R2& X) {
+ // const R2 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // return -1.3 * x * x - 0.2 * y * y + 0.7;
+ // };
+
+ // auto p2 = [&inv_T](const R2& X) {
+ // const R2 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // return +2.6 * x * x - 1.4 * y * y - 1.9;
+ // };
+
+ // auto p3 = [&inv_T](const R2& X) {
+ // const R2 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // return -0.6 * x * x + 1.4 * y * y + 2.3;
+ // };
+
+ // auto q0 = [&inv_T](const R2& X) {
+ // const R2 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // return 3 * x * y;
+ // };
+
+ // auto q1 = [&inv_T](const R2& X) {
+ // const R2 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // return -1.3 * x * y;
+ // };
+
+ // for (auto descriptor : descriptor_list) {
+ // SECTION(name(descriptor.integrationMethodType()))
+ // {
+ // SECTION("R data")
+ // {
+ // auto R_exact = p0;
+
+ // DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ // auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // SECTION("R2x2 data")
+ // {
+ // using R2x2 = TinyMatrix<2>;
+ // auto R2x2_exact = [&](const R2& X) -> R2x2 {
+ // return T * TinyMatrix<2>{p0(X), q0(X), q1(X), p1(X)} * inv_T;
+ // };
+
+ // DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ // auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R2x2_exact));
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // SECTION("vector of R")
+ // {
+ // std::array<std::function<double(const R2&)>, 4> vector_exact //
+ // = {p0, p1, p2, p3};
+
+ // DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ // auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // SECTION("vector of R2x2")
+ // {
+ // using R2x2 = TinyMatrix<2>;
+
+ // std::array<std::function<R2x2(const R2&)>, 2> vector_R2x2_exact //
+ // = {[&](const R2& X) {
+ // return T * R2x2{p0(X), q0(X), q1(X), p1(X)} * inv_T;
+ // },
+ // [&](const R2& X) {
+ // return T * R2x2{p0(X), q1(X), 0, p1(X)} * inv_T;
+ // }};
+
+ // DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_R2x2_exact);
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ // auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2x2>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_R2x2_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // SECTION("list")
+ // {
+ // using R2x2 = TinyMatrix<2>;
+
+ // auto R_exact = p0;
+
+ // DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
+
+ // std::array<std::function<double(const R2&)>, 4> vector_exact //
+ // = {p0, p1, p2, p3};
+
+ // DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+ // std::array<std::function<R2x2(const R2&)>, 2> vector_R2x2_exact //
+ // = {[&](const R2& X) {
+ // return T * R2x2{p0(X), q0(X), q1(X), p1(X)} * inv_T;
+ // },
+ // [&](const R2& X) {
+ // return T * R2x2{p2(X), q1(X), 0, p3(X)} * inv_T;
+ // }};
+
+ // DiscreteFunctionP0Vector Wh = test_only::exact_projection(mesh, degree, vector_R2x2_exact);
+
+ // auto R2x2_exact = [&](const R2& X) -> R2x2 {
+ // return T * TinyMatrix<2>{p0(X), q1(X), q0(X), p3(X)} * inv_T;
+ // };
+
+ // DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(uh, Vh, Wh, Ah);
+
+ // auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
+ // auto dpk_Vh = reconstructions[1]->get<DiscreteFunctionDPkVector<2, const double>>();
+ // auto dpk_Wh = reconstructions[2]->get<DiscreteFunctionDPkVector<2, const R2x2>>();
+ // auto dpk_Ah = reconstructions[3]->get<DiscreteFunctionDPk<2, const R2x2>>();
+
+ // {
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // {
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // {
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Wh, vector_R2x2_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // {
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R2x2_exact));
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // }
+ // }
+ // }
+ // }
+
+ // SECTION("3D")
+ // {
+ // using R3 = TinyVector<3>;
+
+ // auto p_initial_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+ // auto& initial_mesh = *p_initial_mesh;
+
+ // constexpr double theta = 1;
+ // TinyMatrix<3> T{std::cos(theta), -std::sin(theta), 0,
+ // //
+ // std::sin(theta), std::cos(theta), 0,
+ // //
+ // 0, 0, 1};
+
+ // auto xr = initial_mesh.xr();
+
+ // NodeValue<R3> new_xr{initial_mesh.connectivity()};
+ // parallel_for(
+ // initial_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { new_xr[node_id] = T * xr[node_id]; });
+
+ // std::shared_ptr p_mesh = std::make_shared<const Mesh<3>>(initial_mesh.shared_connectivity(), new_xr);
+ // const Mesh<3>& mesh = *p_mesh;
+ // // inverse rotation
+ // TinyMatrix<3> inv_T{std::cos(theta), std::sin(theta), 0,
+ // //
+ // -std::sin(theta), std::cos(theta), 0,
+ // //
+ // 0, 0, 1};
+
+ // auto p0 = [&inv_T](const R3& X) {
+ // const R3 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // const double z = Y[2] - 1;
+ // return +1.7 * x * x + 2 * y * y + 1.3 * z * z - 1.1;
+ // };
+
+ // auto p1 = [&inv_T](const R3& X) {
+ // const R3 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // const double z = Y[2] - 1;
+ // return +2.1 * x * x - 1.4 * y * y - 3.1 * z * z + 2.2;
+ // };
+
+ // auto p2 = [&inv_T](const R3& X) {
+ // const R3 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // const double z = Y[2] - 1;
+ // return -1.1 * x * x - 1.2 * y * y + 1.3 * z * z - 1.7;
+ // };
+
+ // auto p3 = [&inv_T](const R3& X) {
+ // const R3 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // const double z = Y[2] - 1;
+ // return 1.9 * x * x + 2.1 * y * y - 3.1 * z * z + 1.6;
+ // };
+
+ // auto p4 = [&inv_T](const R3& X) {
+ // const R3 Y = inv_T * X;
+ // const double x = Y[0] - 2;
+ // const double y = Y[1] - 1;
+ // const double z = Y[2] - 1;
+ // return -2.4 * x * x + 3.3 * y * y - 1.7 * z * z + 2.1;
+ // };
+
+ // SECTION("3 symmetries")
+ // {
+ // std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ // {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
+ // std::vector<std::shared_ptr<
+ // const IBoundaryDescriptor>>{std::
+ // make_shared<NamedBoundaryDescriptor>(
+ // "XMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "YMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "ZMAX")}},
+ // PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
+ // std::vector<std::shared_ptr<
+ // const IBoundaryDescriptor>>{std::
+ // make_shared<NamedBoundaryDescriptor>(
+ // "XMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "YMAX"),
+ // std::make_shared<NamedBoundaryDescriptor>(
+ // "ZMAX")}}};
+
+ // for (auto descriptor : descriptor_list) {
+ // SECTION("R data")
+ // {
+ // auto R_exact = p0;
+
+ // DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ // auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // SECTION("vector of R")
+ // {
+ // std::array<std::function<double(const R3&)>, 4> vector_exact //
+ // = {p0, p1, p2, p3};
+
+ // DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ // auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // }
+ // }
+
+ // SECTION("1 symmetry")
+ // {
+ // // Matrix and their transformations are kept simple to
+ // // derive exact solutions
+ // std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ // {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
+ // std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+ // std::make_shared<NamedBoundaryDescriptor>("XMAX")}},
+ // PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
+ // std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+ // std::make_shared<NamedBoundaryDescriptor>("XMAX")}}};
+ // for (auto descriptor : descriptor_list) {
+ // SECTION(name(descriptor.integrationMethodType()))
+ // {
+ // SECTION("R3x3 data")
+ // {
+ // // Matrix and their transformations are kept simple to
+ // // derive exact solutions
+
+ // using R3x3 = TinyMatrix<3>;
+ // auto R2x2_exact = [&](const R3& X) -> R3x3 {
+ // return T * TinyMatrix<3>{p0(X), 0, 0, //
+ // 0, p1(X), p2(X), //
+ // 0, p3(X), p4(X)} *
+ // inv_T;
+ // };
+
+ // DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+ // auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3x3>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R2x2_exact));
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // SECTION("vector of R3x3")
+ // {
+ // using R3x3 = TinyMatrix<3>;
+
+ // std::array<std::function<R3x3(const R3&)>, 2> vector_R3x3_exact //
+ // = {[&](const R3& X) {
+ // return T * R3x3{p0(X), 0, 0, //
+ // 0, p1(X), p2(X), //
+ // 0, p3(X), p4(X)} *
+ // inv_T;
+ // },
+ // [&](const R3& X) {
+ // return T * R3x3{p0(X), 0, 0, //
+ // 0, p2(X), p4(X), //
+ // 0, p3(X), p1(X)} *
+ // inv_T;
+ // }};
+
+ // DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_R3x3_exact);
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ // auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3x3>>();
+
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_R3x3_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+
+ // SECTION("list")
+ // {
+ // using R3x3 = TinyMatrix<3>;
+
+ // auto R_exact = p0;
+
+ // DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
+
+ // std::array<std::function<double(const R3&)>, 4> vector_exact //
+ // = {p0, p1, p2, p3};
+
+ // DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
+
+ // std::array<std::function<R3x3(const R3&)>, 2> vector_R3x3_exact //
+ // = {[&](const R3& X) {
+ // return T * R3x3{p1(X), 0, 0, //
+ // 0, p3(X), p0(X), //
+ // 0, p2(X), p4(X)} *
+ // inv_T;
+ // },
+ // [&](const R3& X) {
+ // return T * R3x3{p2(X), 0, 0, //
+ // 0, p0(X), p1(X), //
+ // 0, p3(X), p4(X)} *
+ // inv_T;
+ // }};
+
+ // DiscreteFunctionP0Vector Wh = test_only::exact_projection(mesh, degree, vector_R3x3_exact);
+
+ // auto R3x3_exact = [&](const R3& X) -> R3x3 {
+ // return T * R3x3{p0(X), 0, 0, //
+ // 0, p1(X), p2(X), //
+ // 0, p3(X), p4(X)} *
+ // inv_T;
+ // };
+
+ // DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R3x3_exact));
+
+ // auto reconstructions = PolynomialReconstruction{descriptor}.build(uh, Vh, Wh, Ah);
+
+ // auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
+ // auto dpk_Vh = reconstructions[1]->get<DiscreteFunctionDPkVector<3, const double>>();
+ // auto dpk_Wh = reconstructions[2]->get<DiscreteFunctionDPkVector<3, const R3x3>>();
+ // auto dpk_Ah = reconstructions[3]->get<DiscreteFunctionDPk<3, const R3x3>>();
+
+ // {
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // {
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // {
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Wh, vector_R3x3_exact);
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // {
+ // double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R3x3_exact));
+ // REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
+ // }
+ // }
+ // }
+ // }
+ // }
+ // }
+ // }
+}