From 1d505a118cdcb7be7bf431caf39b0ec264ac24b9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 14 Apr 2025 13:44:16 +0200
Subject: [PATCH] Add tests for PolynomialReconstruction of degree 2
---
...test_PolynomialReconstruction_degree_2.cpp | 758 +++++++++++++-----
1 file changed, 562 insertions(+), 196 deletions(-)
diff --git a/tests/test_PolynomialReconstruction_degree_2.cpp b/tests/test_PolynomialReconstruction_degree_2.cpp
index 435084af..a1c5c058 100644
--- a/tests/test_PolynomialReconstruction_degree_2.cpp
+++ b/tests/test_PolynomialReconstruction_degree_2.cpp
@@ -477,201 +477,567 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
SECTION("with symmetries")
{
-#warning continue here
- // 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));
- // }
- // }
- // }
- // }
+ 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));
+ }
+ }
+ }
+ }
+ }
+ }
}
}
--
GitLab