From f69b27d8dfc78639e60e4f05c0ad72f33efadf84 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 26 Jul 2024 15:15:06 +0200
Subject: [PATCH] Add tests for PolynomialReconstruction using element
quadrature
---
tests/test_PolynomialReconstruction.cpp | 1542 ++++++++++++-----------
1 file changed, 777 insertions(+), 765 deletions(-)
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index 3efabce26..cdfb532f2 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -24,948 +24,960 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
SECTION("degree 1")
{
- PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::cell_center, 1};
+ std::vector<PolynomialReconstructionDescriptor> descriptor_list = {
+ PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1},
+ };
- SECTION("1D")
- {
- using R1 = TinyVector<1>;
-
- SECTION("R data")
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
{
- 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_affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<double> fh{p_mesh};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
-
- auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
-
- auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
+ SECTION("1D")
+ {
+ using R1 = TinyVector<1>;
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
-
- {
- double max_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const double reconstructed_slope =
- (dpk_fh[cell_id](R1{0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
-
- max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
- }
- REQUIRE(parallel::allReduceMax(max_slope_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())
+ SECTION("R data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
- auto& mesh = *p_mesh;
+ 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_affine = [](const R1& x) -> R3 {
- return R3{+2.3 + 1.7 * x[0], //
- +1.4 - 0.6 * x[0], //
- -0.2 + 3.1 * x[0]};
- };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+ auto R_affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<R3> uh{p_mesh};
+ DiscreteFunctionP0<double> fh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
- auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
- auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3>>();
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- {
- double max_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R3 reconstructed_slope =
- (1 / 0.2) * (dpk_uh[cell_id](R1{0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R1{0.1}));
+ {
+ double max_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const double reconstructed_slope =
+ (dpk_fh[cell_id](R1{0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
- max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+ max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
+ }
+ REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
+ }
}
- REQUIRE(parallel::allReduceMax(max_slope_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())
+ SECTION("R^3 data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
- auto& mesh = *p_mesh;
+ using R3 = TinyVector<3>;
- auto R3x3_affine = [](const R1& x) -> R3x3 {
- return R3x3{
- +2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0], //
- +2.4 - 2.3 * x[0], -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0],
- -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
- };
- };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+ for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
- DiscreteFunctionP0<R3x3> Ah{p_mesh};
+ auto R3_affine = [](const R1& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0], //
+ +1.4 - 0.6 * x[0], //
+ -0.2 + 3.1 * x[0]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); });
+ DiscreteFunctionP0<R3> uh{p_mesh};
- auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
- auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>();
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3>>();
- {
- double max_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R3x3 reconstructed_slope =
- (1 / 0.2) * (dpk_Ah[cell_id](R1{0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R1{0.1}));
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- R3x3 slops = R3x3{+1.7, +2.1, -0.6, //
- -2.3, +3.1, -3.6, //
- +3.1, +2.9, +2.3};
+ {
+ double max_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3 reconstructed_slope =
+ (1 / 0.2) * (dpk_uh[cell_id](R1{0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R1{0.1}));
- max_slope_error = std::max(max_slope_error, //
- frobeniusNorm(reconstructed_slope - slops));
+ max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+ }
+ REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
+ }
}
- REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-13));
}
}
- }
- }
-
- SECTION("R vector data")
- {
- using R3 = TinyVector<3>;
- for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
- SECTION(named_mesh.name())
+ SECTION("R^3x3 data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
- auto& mesh = *p_mesh;
-
- auto vector_affine = [](const R1& x) -> R3 {
- return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
- };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- auto vector = vector_affine(xj[cell_id]);
- for (size_t i = 0; i < vector.dimension(); ++i) {
- Vh[cell_id][i] = vector[i];
+ 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_affine = [](const R1& x) -> R3x3 {
+ return R3x3{
+ +2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0], //
+ +2.4 - 2.3 * x[0], -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0],
+ -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
+ };
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0<R3x3> Ah{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
- });
- auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+ {
+ double max_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3x3 reconstructed_slope =
+ (1 / 0.2) * (dpk_Ah[cell_id](R1{0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R1{0.1}));
- auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
+ R3x3 slops = R3x3{+1.7, +2.1, -0.6, //
+ -2.3, +3.1, -3.6, //
+ +3.1, +2.9, +2.3};
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- auto vector = vector_affine(xj[cell_id]);
- for (size_t i = 0; i < vector.dimension(); ++i) {
- max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+ max_slope_error = std::max(max_slope_error, //
+ frobeniusNorm(reconstructed_slope - slops));
+ }
+ REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-13));
}
}
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
+ }
+
+ SECTION("R 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;
+
+ auto vector_affine = [](const R1& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+ for (size_t i = 0; i < vector.dimension(); ++i) {
+ Vh[cell_id][i] = vector[i];
+ }
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+ for (size_t i = 0; i < vector.dimension(); ++i) {
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- {
- double max_slope_error = 0;
- const TinyVector<3> slope{+1.7, +2.1, -0.6};
+ {
+ double max_slope_error = 0;
+ const TinyVector<3> slope{+1.7, +2.1, -0.6};
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- for (size_t i = 0; i < slope.dimension(); ++i) {
- const double reconstructed_slope =
- (1 / 0.2) * (dpk_Vh(cell_id, i)(R1{0.1} + xj[cell_id]) - dpk_Vh(cell_id, i)(xj[cell_id] - R1{0.1}));
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < slope.dimension(); ++i) {
+ const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R1{0.1} + xj[cell_id]) -
+ dpk_Vh(cell_id, i)(xj[cell_id] - R1{0.1}));
- max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i]));
+ max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
}
}
- REQUIRE(parallel::allReduceMax(max_slope_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())
+ SECTION("list of various types")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
- auto& mesh = *p_mesh;
-
- auto R_affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
-
- auto R3_affine = [](const R1& x) -> R3 {
- return R3{+2.3 + 1.7 * x[0], //
- +1.4 - 0.6 * x[0], //
- -0.2 + 3.1 * x[0]};
- };
-
- auto R3x3_affine = [](const R1& x) -> R3x3 {
- return R3x3{
- +2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0], //
- +2.4 - 2.3 * x[0], -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0],
- -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
- };
- };
-
- auto vector_affine = [](const R1& x) -> R3 {
- return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
- };
-
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<double> fh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
-
- DiscreteFunctionP0<R3> uh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
-
- DiscreteFunctionP0<R3x3> Ah{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); });
-
- DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- auto vector = vector_affine(xj[cell_id]);
- for (size_t i = 0; i < vector.dimension(); ++i) {
- Vh[cell_id][i] = vector[i];
- }
- });
-
- 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_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
+ 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_affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
+
+ auto R3_affine = [](const R1& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0], //
+ +1.4 - 0.6 * x[0], //
+ -0.2 + 3.1 * x[0]};
+ };
+
+ auto R3x3_affine = [](const R1& x) -> R3x3 {
+ return R3x3{
+ +2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0], //
+ +2.4 - 2.3 * x[0], -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0],
+ -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0],
+ };
+ };
+
+ auto vector_affine = [](const R1& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
+ };
+
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0<double> fh{p_mesh};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
+
+ DiscreteFunctionP0<R3> uh{p_mesh};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+
+ DiscreteFunctionP0<R3x3> Ah{p_mesh};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); });
+
+ DiscreteFunctionP0Vector<double> Vh{p_mesh, 3};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+ for (size_t i = 0; i < vector.dimension(); ++i) {
+ Vh[cell_id][i] = vector[i];
+ }
+ });
+
+ 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_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- {
- double max_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const double reconstructed_slope =
- (dpk_fh[cell_id](R1{0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
+ {
+ double max_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const double reconstructed_slope =
+ (dpk_fh[cell_id](R1{0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
- max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
- }
- REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
- }
+ max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
+ }
+ REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
+ }
- auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>();
+ auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>();
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- {
- double max_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R3 reconstructed_slope =
- (1 / 0.2) * (dpk_uh[cell_id](R1{0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R1{0.1}));
+ {
+ double max_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3 reconstructed_slope =
+ (1 / 0.2) * (dpk_uh[cell_id](R1{0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R1{0.1}));
- max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
- }
- REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
- }
+ max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+ }
+ REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
+ }
- auto dpk_Ah = reconstructions[2]->get<DiscreteFunctionDPk<1, const R3x3>>();
+ auto dpk_Ah = reconstructions[2]->get<DiscreteFunctionDPk<1, const R3x3>>();
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- {
- double max_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R3x3 reconstructed_slope =
- (1 / 0.2) * (dpk_Ah[cell_id](R1{0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R1{0.1}));
+ {
+ double max_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3x3 reconstructed_slope =
+ (1 / 0.2) * (dpk_Ah[cell_id](R1{0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R1{0.1}));
- R3x3 slops = R3x3{+1.7, +2.1, -0.6, //
- -2.3, +3.1, -3.6, //
- +3.1, +2.9, +2.3};
+ R3x3 slops = R3x3{+1.7, +2.1, -0.6, //
+ -2.3, +3.1, -3.6, //
+ +3.1, +2.9, +2.3};
- max_slope_error = std::max(max_slope_error, //
- frobeniusNorm(reconstructed_slope - slops));
- }
- REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-13));
- }
-
- auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>();
+ max_slope_error = std::max(max_slope_error, //
+ frobeniusNorm(reconstructed_slope - slops));
+ }
+ REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- auto vector = vector_affine(xj[cell_id]);
- for (size_t i = 0; i < vector.dimension(); ++i) {
- max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+ auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+ for (size_t i = 0; i < vector.dimension(); ++i) {
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
- {
- double max_slope_error = 0;
- const TinyVector<3> slope{+1.7, +2.1, -0.6};
+ {
+ double max_slope_error = 0;
+ const TinyVector<3> slope{+1.7, +2.1, -0.6};
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- for (size_t i = 0; i < slope.dimension(); ++i) {
- const double reconstructed_slope =
- (1 / 0.2) * (dpk_Vh(cell_id, i)(R1{0.1} + xj[cell_id]) - dpk_Vh(cell_id, i)(xj[cell_id] - R1{0.1}));
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < slope.dimension(); ++i) {
+ const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R1{0.1} + xj[cell_id]) -
+ dpk_Vh(cell_id, i)(xj[cell_id] - R1{0.1}));
- max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i]));
+ max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
}
}
- REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
}
}
}
- }
- }
- SECTION("2D")
- {
- using R2 = TinyVector<2>;
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
- SECTION("R data")
- {
- for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
- SECTION(named_mesh.name())
+ SECTION("R data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
- auto& mesh = *p_mesh;
+ 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_affine = [](const R2& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1]; };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+ auto R_affine = [](const R2& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1]; };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<double> fh{p_mesh};
+ DiscreteFunctionP0<double> fh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
- auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
- auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- {
- double max_x_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const double reconstructed_slope =
- (dpk_fh[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0.1, 0})) / 0.2;
+ {
+ double max_x_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const double reconstructed_slope =
+ (dpk_fh[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0.1, 0})) / 0.2;
- max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
- }
- REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
- }
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
+ }
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
- {
- double max_y_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const double reconstructed_slope =
- (dpk_fh[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0, 0.1})) / 0.2;
+ {
+ double max_y_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const double reconstructed_slope =
+ (dpk_fh[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0, 0.1})) / 0.2;
- max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-14));
+ }
}
- REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-14));
}
}
- }
- }
-
- SECTION("R^3 data")
- {
- using R3 = TinyVector<3>;
- for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
- SECTION(named_mesh.name())
+ SECTION("R^3 data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
- auto& mesh = *p_mesh;
+ using R3 = TinyVector<3>;
- auto R3_affine = [](const R2& x) -> R3 {
- return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1], //
- +1.4 - 0.6 * x[0] + 1.3 * x[1], //
- -0.2 + 3.1 * x[0] - 1.1 * x[1]};
- };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+ for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+ auto& mesh = *p_mesh;
- DiscreteFunctionP0<R3> uh{p_mesh};
+ auto R3_affine = [](const R2& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1], //
+ +1.4 - 0.6 * x[0] + 1.3 * x[1], //
+ -0.2 + 3.1 * x[0] - 1.1 * x[1]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+ DiscreteFunctionP0<R3> uh{p_mesh};
- auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
- auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R3>>();
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R3>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- {
- double max_x_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R3 reconstructed_slope =
- (1 / 0.2) * (dpk_uh[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R2{0.1, 0}));
+ {
+ double max_x_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R2{0.1, 0} + xj[cell_id]) -
+ dpk_uh[cell_id](xj[cell_id] - R2{0.1, 0}));
- max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
- }
- REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
- }
+ max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+ }
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
- {
- double max_y_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R3 reconstructed_slope =
- (1 / 0.2) * (dpk_uh[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_uh[cell_id](xj[cell_id] - R2{0, 0.1}));
+ {
+ double max_y_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R2{0, 0.1} + xj[cell_id]) -
+ dpk_uh[cell_id](xj[cell_id] - R2{0, 0.1}));
- max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
+ max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
}
- REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
}
}
- }
- }
-
- SECTION("R^2x2 data")
- {
- using R2x2 = TinyMatrix<2, 2>;
- for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
- SECTION(named_mesh.name())
+ SECTION("R^2x2 data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
- auto& mesh = *p_mesh;
+ using R2x2 = TinyMatrix<2, 2>;
- auto R2x2_affine = [](const R2& x) -> R2x2 {
- return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1], //
- +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
- };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+ for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
+ auto& mesh = *p_mesh;
- DiscreteFunctionP0<R2x2> Ah{p_mesh};
+ auto R2x2_affine = [](const R2& x) -> R2x2 {
+ return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1], //
+ +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
+ DiscreteFunctionP0<R2x2> Ah{p_mesh};
- auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
- auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>();
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
-
- {
- double max_x_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R2x2 reconstructed_slope =
- (1 / 0.2) * (dpk_Ah[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R2{0.1, 0}));
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>();
- max_x_slope_error =
- std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, +2.1, //
- -0.6, -2.3}));
- }
- REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
- }
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- {
- double max_y_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R2x2 reconstructed_slope =
- (1 / 0.2) * (dpk_Ah[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_Ah[cell_id](xj[cell_id] - R2{0, 0.1}));
+ {
+ double max_x_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R2{0.1, 0} + xj[cell_id]) -
+ dpk_Ah[cell_id](xj[cell_id] - R2{0.1, 0}));
+
+ max_x_slope_error =
+ std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, +2.1, //
+ -0.6, -2.3}));
+ }
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
- max_y_slope_error =
- std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, //
- -2.1, +1.3}));
+ {
+ double max_y_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R2{0, 0.1} + xj[cell_id]) -
+ dpk_Ah[cell_id](xj[cell_id] - R2{0, 0.1}));
+
+ max_y_slope_error =
+ std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, //
+ -2.1, +1.3}));
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
}
- REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
}
}
- }
- }
- SECTION("vector data")
- {
- using R4 = TinyVector<4>;
-
- for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
- SECTION(named_mesh.name())
+ SECTION("vector data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
- auto& mesh = *p_mesh;
-
- auto vector_affine = [](const R2& x) -> R4 {
- return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1], //
- +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
- };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- auto vector = vector_affine(xj[cell_id]);
- for (size_t i = 0; i < vector.dimension(); ++i) {
- Vh[cell_id][i] = vector[i];
+ using R4 = TinyVector<4>;
+
+ 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 vector_affine = [](const R2& x) -> R4 {
+ return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1], //
+ +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+ for (size_t i = 0; i < vector.dimension(); ++i) {
+ Vh[cell_id][i] = vector[i];
+ }
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+ for (size_t i = 0; i < vector.dimension(); ++i) {
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
- });
-
- auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
- auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>();
-
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- auto vector = vector_affine(xj[cell_id]);
- for (size_t i = 0; i < vector.dimension(); ++i) {
- max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+ {
+ double max_x_slope_error = 0;
+ const R4 slope{+1.7, +2.1, -0.6, -2.3};
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < slope.dimension(); ++i) {
+ const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0.1, 0} + xj[cell_id]) -
+ dpk_Vh(cell_id, i)(xj[cell_id] - R2{0.1, 0}));
+
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
}
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
- {
- double max_x_slope_error = 0;
- const R4 slope{+1.7, +2.1, -0.6, -2.3};
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- for (size_t i = 0; i < slope.dimension(); ++i) {
- const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0.1, 0} + xj[cell_id]) -
- dpk_Vh(cell_id, i)(xj[cell_id] - R2{0.1, 0}));
-
- max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
- }
- }
- REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
- }
-
- {
- double max_y_slope_error = 0;
- const R4 slope{+1.2, -2.2, -2.1, +1.3};
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- for (size_t i = 0; i < slope.dimension(); ++i) {
- const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0, 0.1} + xj[cell_id]) -
- dpk_Vh(cell_id, i)(xj[cell_id] - R2{0, 0.1}));
-
- max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+ {
+ double max_y_slope_error = 0;
+ const R4 slope{+1.2, -2.2, -2.1, +1.3};
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < slope.dimension(); ++i) {
+ const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R2{0, 0.1} + xj[cell_id]) -
+ dpk_Vh(cell_id, i)(xj[cell_id] - R2{0, 0.1}));
+
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
}
}
- REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
}
}
}
- }
- }
- SECTION("3D")
- {
- using R3 = TinyVector<3>;
+ SECTION("3D")
+ {
+ using R3 = TinyVector<3>;
- SECTION("R data")
- {
- for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
- SECTION(named_mesh.name())
+ SECTION("R data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
- auto& mesh = *p_mesh;
+ 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_affine = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+ auto R_affine = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<double> fh{p_mesh};
+ DiscreteFunctionP0<double> fh{p_mesh};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
- auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
- auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- {
- double max_x_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const double reconstructed_slope =
- (dpk_fh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0.1, 0, 0})) / 0.2;
+ {
+ double max_x_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const double reconstructed_slope =
+ (dpk_fh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0.1, 0, 0})) /
+ 0.2;
- max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
- }
- REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
- }
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
+ }
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
- {
- double max_y_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const double reconstructed_slope =
- (dpk_fh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0.1, 0})) / 0.2;
+ {
+ double max_y_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const double reconstructed_slope =
+ (dpk_fh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0.1, 0})) /
+ 0.2;
- max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
- }
- REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
- }
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
- {
- double max_z_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const double reconstructed_slope =
- (dpk_fh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0, 0.1})) / 0.2;
+ {
+ double max_z_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const double reconstructed_slope =
+ (dpk_fh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0, 0.1})) /
+ 0.2;
- max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - 2.1));
+ max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - 2.1));
+ }
+ REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
+ }
}
- REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
}
}
- }
- }
- SECTION("R^3 data")
- {
- for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
- SECTION(named_mesh.name())
+ SECTION("R^3 data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
- auto& mesh = *p_mesh;
-
- auto R3_affine = [](const R3& x) -> R3 {
- return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2], //
- +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2], //
- -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]};
- };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<R3> uh{p_mesh};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
-
- auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
-
- auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
-
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
+ 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_affine = [](const R3& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2], //
+ +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2], //
+ -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0<R3> uh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
+
+ auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- {
- double max_x_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
- dpk_uh[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+ {
+ double max_x_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
+ dpk_uh[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
- max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
- }
- REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-12));
- }
+ max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+ }
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-12));
+ }
- {
- double max_y_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
- dpk_uh[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+ {
+ double max_y_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
+ dpk_uh[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
- max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
- }
- REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
- }
+ max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
+ }
- {
- double max_z_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
- dpk_uh[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+ {
+ double max_z_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
+ dpk_uh[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
- max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R3{1.8, -3.7, 1.9}));
+ max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R3{1.8, -3.7, 1.9}));
+ }
+ REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
+ }
}
- REQUIRE(parallel::allReduceMax(max_z_slope_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())
+ SECTION("R^2x2 data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
- auto& mesh = *p_mesh;
-
- auto R2x2_affine = [](const R3& x) -> R2x2 {
- return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
- //
- +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
- };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0<R2x2> Ah{p_mesh};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
-
- descriptor.setRowWeighting(false);
- auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
-
- auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
-
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- max_mean_error =
- std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
-
- {
- double max_x_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
- dpk_Ah[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
-
- max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, 2.1, //
- -2.3, +3.1}));
- }
- REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
- }
-
- {
- double max_y_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
- dpk_Ah[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+ 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_affine = [](const R3& x) -> R2x2 {
+ return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
+ //
+ +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0<R2x2> Ah{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
+
+ descriptor.setRowWeighting(false);
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
+
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
- max_y_slope_error =
- std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, //
- 1.3, +0.8}));
- }
- REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
- }
+ {
+ double max_x_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
+ dpk_Ah[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+
+ max_x_slope_error =
+ std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, 2.1, //
+ -2.3, +3.1}));
+ }
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
- {
- double max_z_slope_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
- dpk_Ah[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+ {
+ double max_y_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
+ dpk_Ah[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+
+ max_y_slope_error =
+ std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, //
+ 1.3, +0.8}));
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
+ }
- max_z_slope_error =
- std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R2x2{-1.3, -2.4, //
- +1.4, -1.8}));
+ {
+ double max_z_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
+ dpk_Ah[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+
+ max_z_slope_error =
+ std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R2x2{-1.3, -2.4, //
+ +1.4, -1.8}));
+ }
+ REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
+ }
}
- REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
}
}
- }
- }
-
- SECTION("vector data")
- {
- using R4 = TinyVector<4>;
- for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
- SECTION(named_mesh.name())
+ SECTION("vector data")
{
- auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
- auto& mesh = *p_mesh;
-
- auto vector_affine = [](const R3& x) -> R4 {
- return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
- //
- +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
- };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- auto vector = vector_affine(xj[cell_id]);
- for (size_t i = 0; i < vector.dimension(); ++i) {
- Vh[cell_id][i] = vector[i];
- }
- });
-
- descriptor.setPreconditioning(false);
- auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
-
- auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
-
- {
- double max_mean_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- auto vector = vector_affine(xj[cell_id]);
- for (size_t i = 0; i < vector.dimension(); ++i) {
- max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+ using R4 = TinyVector<4>;
+
+ for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
+ auto& mesh = *p_mesh;
+
+ auto vector_affine = [](const R3& x) -> R4 {
+ return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
+ //
+ +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+ for (size_t i = 0; i < vector.dimension(); ++i) {
+ Vh[cell_id][i] = vector[i];
+ }
+ });
+
+ descriptor.setPreconditioning(false);
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+ for (size_t i = 0; i < vector.dimension(); ++i) {
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
-
- {
- double max_x_slope_error = 0;
- const R4 slope{+1.7, 2.1, -2.3, +3.1};
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- for (size_t i = 0; i < slope.dimension(); ++i) {
- const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0.1, 0, 0} + xj[cell_id]) -
- dpk_Vh(cell_id, i)(xj[cell_id] - R3{0.1, 0, 0}));
- max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+ {
+ double max_x_slope_error = 0;
+ const R4 slope{+1.7, 2.1, -2.3, +3.1};
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < slope.dimension(); ++i) {
+ const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0.1, 0, 0} + xj[cell_id]) -
+ dpk_Vh(cell_id, i)(xj[cell_id] - R3{0.1, 0, 0}));
+
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
}
- }
- REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
- }
-
- {
- double max_y_slope_error = 0;
- const R4 slope{+1.2, -2.2, 1.3, +0.8};
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- for (size_t i = 0; i < slope.dimension(); ++i) {
- const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0.1, 0} + xj[cell_id]) -
- dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0.1, 0}));
- max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+ {
+ double max_y_slope_error = 0;
+ const R4 slope{+1.2, -2.2, 1.3, +0.8};
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < slope.dimension(); ++i) {
+ const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0.1, 0} + xj[cell_id]) -
+ dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0.1, 0}));
+
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
}
- }
- REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
- }
-
- {
- double max_z_slope_error = 0;
- const R4 slope{-1.3, -2.4, +1.4, -1.8};
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- for (size_t i = 0; i < slope.dimension(); ++i) {
- const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0, 0.1} + xj[cell_id]) -
- dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0, 0.1}));
- max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - slope[i]));
+ {
+ double max_z_slope_error = 0;
+ const R4 slope{-1.3, -2.4, +1.4, -1.8};
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < slope.dimension(); ++i) {
+ const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0, 0.1} + xj[cell_id]) -
+ dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0, 0.1}));
+
+ max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - slope[i]));
+ }
+ }
+ REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-13));
}
}
- REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-13));
}
}
}
- }
- }
- SECTION("errors")
- {
- auto p_mesh1 = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
- DiscreteFunctionP0<double> f1{p_mesh1};
+ SECTION("errors")
+ {
+ auto p_mesh1 = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+ DiscreteFunctionP0<double> f1{p_mesh1};
- auto p_mesh2 = MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
- DiscreteFunctionP0<double> f2{p_mesh2};
+ auto p_mesh2 = MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
+ DiscreteFunctionP0<double> f2{p_mesh2};
- REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(f1, f2),
- "error: cannot reconstruct functions living of different meshes simultaneously");
+ REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(f1, f2),
+ "error: cannot reconstruct functions living of different meshes simultaneously");
+ }
+ }
}
}
}
--
GitLab