From 0db0579a0692eac1f48fc5ed94571c38ee90a4ec Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 2 Jul 2024 07:47:38 +0200
Subject: [PATCH] Use new polynomial reconstruction interface in tests
---
tests/test_PolynomialReconstruction.cpp | 369 ++++++++++++++++--------
1 file changed, 241 insertions(+), 128 deletions(-)
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index 3c3d3e0dc..123072266 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -33,33 +33,22 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
auto& mesh = *p_mesh;
- auto affine = [](const R1& x) { return 2.3 + 1.7 * 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<double> fh{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
- DiscreteFunctionP0<double> gh{p_mesh};
+ auto reconstructions = PolynomialReconstruction{}.build(fh);
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { gh[cell_id] = 0; });
-
- DiscreteFunctionVariant gh_v(gh);
- std::shared_ptr<const DiscreteFunctionVariant> p_gh_v = std::make_shared<DiscreteFunctionVariant>(fh);
-
- std::shared_ptr<DiscreteFunctionP0<const double>> g_fh =
- std::make_shared<DiscreteFunctionP0<const double>>(fh);
-
- auto reconstructions = PolynomialReconstruction{}.build(g_fh, p_gh_v, gh_v, fh);
-
- auto dpk = reconstructions[3]->get<DiscreteFunctionDPk<1, const double>>();
+ 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[cell_id](xj[cell_id]) - affine(xj[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(max_mean_error == Catch::Approx(0).margin(1E-14));
}
@@ -68,7 +57,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
double max_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
- (dpk[cell_id](R1{0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
+ (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));
}
@@ -78,7 +67,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
}
}
- SECTION("R^d data")
+ SECTION("R^3 data")
{
using R3 = TinyVector<3>;
@@ -88,25 +77,26 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
auto& mesh = *p_mesh;
- auto affine = [](const R1& x) -> R3 {
+ 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();
- DiscreteFunctionP0<R3> fh{p_mesh};
+ DiscreteFunctionP0<R3> uh{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+
+ auto reconstructions = PolynomialReconstruction{}.build(uh);
- PolynomialReconstruction reconstruction;
- auto dpk = reconstruction.build(mesh, fh);
+ auto dpk_uh = reconstructions[0]->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[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
}
REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
}
@@ -115,7 +105,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
double max_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R3 reconstructed_slope =
- (1 / 0.2) * (dpk[cell_id](R1{0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R1{0.1}));
+ (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}));
}
@@ -125,9 +115,66 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
}
}
- SECTION("R^mn data")
+ 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_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{}.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(max_mean_error == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ 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};
+
+ max_slope_error = std::max(max_slope_error, //
+ frobeniusNorm(reconstructed_slope - slops));
+ }
+ REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+
+ SECTION("list of various types")
{
- using R33 = TinyMatrix<3, 3>;
+ using R3x3 = TinyMatrix<3>;
+ using R3 = TinyVector<3>;
for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
SECTION(named_mesh.name())
@@ -135,27 +182,88 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
auto& mesh = *p_mesh;
- auto affine = [](const R1& x) -> R33 {
- return R33{
+ 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 xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<R33> 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]); });
+ DiscreteFunctionP0<R3> uh{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+ 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]); });
+
+ auto reconstructions = PolynomialReconstruction{}.build(std::make_shared<DiscreteFunctionVariant>(fh), uh,
+ std::make_shared<DiscreteFunctionP0<R3x3>>(Ah));
+
+ 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(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(max_slope_error == Catch::Approx(0).margin(1E-14));
+ }
+
+ auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>();
- PolynomialReconstruction reconstruction;
- auto dpk = reconstruction.build(mesh, fh);
+ {
+ 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(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}));
+
+ max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
+ }
+ REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
+ }
+
+ 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[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ max_mean_error =
+ std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id])));
}
REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
}
@@ -163,12 +271,12 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R33 reconstructed_slope =
- (1 / 0.2) * (dpk[cell_id](R1{0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R1{0.1}));
+ 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}));
- R33 slops = R33{+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));
@@ -192,21 +300,22 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
auto& mesh = *p_mesh;
- auto 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};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
- PolynomialReconstruction reconstruction;
- auto dpk = reconstruction.build(mesh, fh);
+ auto reconstructions = PolynomialReconstruction{}.build(fh);
+
+ 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[cell_id](xj[cell_id]) - affine(xj[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(max_mean_error == Catch::Approx(0).margin(1E-14));
}
@@ -215,7 +324,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
- (dpk[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0.1, 0})) / 0.2;
+ (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));
}
@@ -226,7 +335,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
- (dpk[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0, 0.1})) / 0.2;
+ (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)));
}
@@ -236,7 +345,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
}
}
- SECTION("R^d data")
+ SECTION("R^3 data")
{
using R3 = TinyVector<3>;
@@ -246,25 +355,26 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
auto& mesh = *p_mesh;
- auto affine = [](const R2& x) -> R3 {
+ 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();
- DiscreteFunctionP0<R3> fh{p_mesh};
+ DiscreteFunctionP0<R3> uh{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+
+ auto reconstructions = PolynomialReconstruction{}.build(uh);
- PolynomialReconstruction reconstruction;
- auto dpk = reconstruction.build(mesh, fh);
+ 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[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
}
REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
}
@@ -273,7 +383,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0.1, 0}));
+ (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}));
}
@@ -284,7 +394,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0, 0.1}));
+ (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}));
}
@@ -294,9 +404,9 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
}
}
- SECTION("R^mn data")
+ SECTION("R^2x2 data")
{
- using R22 = TinyMatrix<2, 2>;
+ using R2x2 = TinyMatrix<2, 2>;
for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
SECTION(named_mesh.name())
@@ -304,24 +414,26 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
auto& mesh = *p_mesh;
- auto affine = [](const R2& x) -> R22 {
- return R22{+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 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();
- DiscreteFunctionP0<R22> fh{p_mesh};
+ DiscreteFunctionP0<R2x2> Ah{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
+
+ auto reconstructions = PolynomialReconstruction{}.build(Ah);
- PolynomialReconstruction reconstruction;
- auto dpk = reconstruction.build(mesh, fh);
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<2, 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[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ max_mean_error =
+ std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
}
REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
}
@@ -329,11 +441,11 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R22 reconstructed_slope =
- (1 / 0.2) * (dpk[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0.1, 0}));
+ 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 - R22{+1.7, +2.1, //
- -0.6, -2.3}));
+ max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, +2.1, //
+ -0.6, -2.3}));
}
REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
}
@@ -341,11 +453,11 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R22 reconstructed_slope =
- (1 / 0.2) * (dpk[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R2{0, 0.1}));
+ 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 - R22{+1.2, -2.2, //
- -2.1, +1.3}));
+ max_y_slope_error = std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, //
+ -2.1, +1.3}));
}
REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
}
@@ -366,21 +478,22 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
auto& mesh = *p_mesh;
- auto 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};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
- PolynomialReconstruction reconstruction;
- auto dpk = reconstruction.build(mesh, fh);
+ auto reconstructions = PolynomialReconstruction{}.build(fh);
+
+ 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[cell_id](xj[cell_id]) - affine(xj[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(max_mean_error == Catch::Approx(0).margin(1E-14));
}
@@ -389,7 +502,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
- (dpk[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0.1, 0, 0})) / 0.2;
+ (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));
}
@@ -400,7 +513,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
- (dpk[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0, 0.1, 0})) / 0.2;
+ (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)));
}
@@ -411,7 +524,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
double max_z_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
- (dpk[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0, 0, 0.1})) / 0.2;
+ (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));
}
@@ -421,36 +534,34 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
}
}
- SECTION("R^d data")
+ SECTION("R^3 data")
{
- 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 affine = [](const R3& x) -> R4 {
- return R4{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2], //
+ 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], //
- -1.2 - 1.5 * x[0] - 0.1 * x[1] - 1.6 * x[2]};
+ -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]};
};
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<R4> fh{p_mesh};
+ DiscreteFunctionP0<R3> uh{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+
+ auto reconstructions = PolynomialReconstruction{}.build(uh);
- PolynomialReconstruction reconstruction;
- auto dpk = reconstruction.build(mesh, fh);
+ 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[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
}
REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
}
@@ -458,10 +569,10 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R4 reconstructed_slope =
- (1 / 0.2) * (dpk[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+ 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 - R4{1.7, -0.6, 3.1, -1.5}));
+ max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
}
REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-13));
}
@@ -469,10 +580,10 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R4 reconstructed_slope =
- (1 / 0.2) * (dpk[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+ 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 - R4{-2.2, 1.3, -1.1, -0.1}));
+ max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
}
REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-13));
}
@@ -480,10 +591,10 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_z_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R4 reconstructed_slope =
- (1 / 0.2) * (dpk[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+ 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 - R4{1.8, -3.7, 1.9, -1.6}));
+ max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R3{1.8, -3.7, 1.9}));
}
REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-13));
}
@@ -491,9 +602,9 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
}
}
- SECTION("R^mn data")
+ SECTION("R^2x2 data")
{
- using R22 = TinyMatrix<2, 2>;
+ using R2x2 = TinyMatrix<2, 2>;
for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
SECTION(named_mesh.name())
@@ -501,25 +612,27 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
auto& mesh = *p_mesh;
- auto affine = [](const R3& x) -> R22 {
- return R22{+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 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<R22> fh{p_mesh};
+ DiscreteFunctionP0<R2x2> Ah{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
+
+ auto reconstructions = PolynomialReconstruction{}.build(Ah);
- PolynomialReconstruction reconstruction;
- auto dpk = reconstruction.build(mesh, fh);
+ 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[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ max_mean_error =
+ std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
}
REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
}
@@ -527,11 +640,11 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R22 reconstructed_slope =
- (1 / 0.2) * (dpk[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+ 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 - R22{+1.7, 2.1, //
- -2.3, +3.1}));
+ max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, 2.1, //
+ -2.3, +3.1}));
}
REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
}
@@ -539,11 +652,11 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R22 reconstructed_slope =
- (1 / 0.2) * (dpk[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+ 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 - R22{+1.2, -2.2, //
- 1.3, +0.8}));
+ max_y_slope_error = std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, //
+ 1.3, +0.8}));
}
REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
}
@@ -551,11 +664,11 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_z_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R22 reconstructed_slope =
- (1 / 0.2) * (dpk[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+ 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 - R22{-1.3, -2.4, //
- +1.4, -1.8}));
+ max_z_slope_error = std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R2x2{-1.3, -2.4, //
+ +1.4, -1.8}));
}
REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12));
}
--
GitLab