diff --git a/tests/test_PolynomialReconstruction_degree_1.cpp b/tests/test_PolynomialReconstruction_degree_1.cpp
index 6c333913e85fa811abe34af119176d3d7209f156..927a5f4c3f965e53c1d940f012dc1bbddf2fbbe1 100644
--- a/tests/test_PolynomialReconstruction_degree_1.cpp
+++ b/tests/test_PolynomialReconstruction_degree_1.cpp
@@ -13,6 +13,9 @@
#include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/LineTransformation.hpp>
+#include <geometry/SquareTransformation.hpp>
+#include <geometry/TetrahedronTransformation.hpp>
+#include <geometry/TriangleTransformation.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/NamedBoundaryDescriptor.hpp>
@@ -23,8 +26,195 @@
#include <MeshDataBaseForTests.hpp>
+#include <type_traits>
+
// clazy:excludeall=non-pod-global-static
+namespace test_only
+{
+
+template <typename DataType>
+PUGS_INLINE double
+get_max_error(const DataType& x, const DataType& y)
+{
+ if constexpr (is_tiny_matrix_v<DataType>) {
+ return frobeniusNorm(x - y);
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ return l2Norm(x - y);
+ } else {
+ static_assert(std::is_arithmetic_v<DataType>, "expecting arithmetic type");
+ return std::abs(x - y);
+ }
+}
+
+template <MeshConcept MeshType, typename DataType>
+double
+max_reconstruction_error(const MeshType& mesh,
+ DiscreteFunctionDPk<MeshType::Dimension, const DataType> dpk_f,
+ std::function<DataType(const TinyVector<MeshType::Dimension>&)> exact)
+{
+ auto xr = mesh.xr();
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ double max_error = 0;
+ auto cell_type = mesh.connectivity().cellType();
+
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto cell_nodes = cell_to_node_matrix[cell_id];
+ if constexpr (MeshType::Dimension == 1) {
+ Assert(cell_type[cell_id] == CellType::Line);
+ LineTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]]};
+ auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ } else if constexpr (MeshType::Dimension == 2) {
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ TriangleTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]]};
+ auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+ xr[cell_nodes[3]]};
+ auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("unexpected cell type");
+ }
+ }
+ } else if constexpr (MeshType::Dimension == 3) {
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ TetrahedronTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]]};
+ auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ }
+ break;
+ }
+ // case CellType::Pyramid: {
+ // SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+ // xr[cell_nodes[3]]};
+ // auto qf = QuadratureManager::instance().getSquareFormula(GaussLobattoQuadratureDescriptor{dpk_f.degree() + 1});
+ // for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ // auto x = T(qf.point(i_quadrarture));
+ // max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ // }
+ // break;
+ // }
+ default: {
+ throw UnexpectedError("unexpected cell type");
+ }
+ }
+ }
+ }
+ return max_error;
+}
+
+template <MeshConcept MeshType, typename DataType, size_t NbComponents>
+double
+max_reconstruction_error(
+ const MeshType& mesh,
+ DiscreteFunctionDPkVector<MeshType::Dimension, const DataType> dpk_v,
+ const std::array<std::function<DataType(const TinyVector<MeshType::Dimension>&)>, NbComponents>& vector_exact)
+{
+ auto xr = mesh.xr();
+ auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ double max_error = 0;
+ auto cell_type = mesh.connectivity().cellType();
+
+ REQUIRE(NbComponents == dpk_v.numberOfComponents());
+
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto cell_nodes = cell_to_node_matrix[cell_id];
+ if constexpr (MeshType::Dimension == 1) {
+ Assert(cell_type[cell_id] == CellType::Line);
+ LineTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]]};
+ auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error = std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ } else if constexpr (MeshType::Dimension == 2) {
+ switch (cell_type[cell_id]) {
+ case CellType::Triangle: {
+ TriangleTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]]};
+ auto qf = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error =
+ std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ break;
+ }
+ case CellType::Quadrangle: {
+ SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+ xr[cell_nodes[3]]};
+ auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error =
+ std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ break;
+ }
+ default: {
+ throw UnexpectedError("unexpected cell type");
+ }
+ }
+ } else if constexpr (MeshType::Dimension == 3) {
+ switch (cell_type[cell_id]) {
+ case CellType::Tetrahedron: {
+ TetrahedronTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]]};
+ auto qf = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
+ for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ auto x = T(qf.point(i_quadrarture));
+ for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
+ max_error =
+ std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
+ }
+ }
+ break;
+ }
+ // case CellType::Pyramid: {
+ // SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
+ // xr[cell_nodes[3]]};
+ // auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
+ // for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
+ // auto x = T(qf.point(i_quadrarture));
+ // max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
+ // }
+ // break;
+ // }
+ default: {
+ throw UnexpectedError("unexpected cell type");
+ }
+ }
+ }
+ }
+ return max_error;
+}
+
+} // namespace test_only
+
TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
{
SECTION("without symmetries")
@@ -49,55 +239,20 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
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();
- auto xr = mesh.xr();
-
- auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ auto R_exact = [](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]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_exact(xj[cell_id]); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
- auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{2});
- {
- double max_error = 0;
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- auto cell_nodes = cell_to_node_matrix[cell_id];
- LineTransformation<1> T{xr[cell_nodes[0]], xr[cell_nodes[1]]};
- for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
- max_error = std::max(max_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
- }
- }
-
- REQUIRE(parallel::allReduceMax(max_error) == 0 // Catch::Approx(0).margin(0 * 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 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));
- }
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
}
}
@@ -112,7 +267,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
auto& mesh = *p_mesh;
- auto R3_affine = [](const R1& x) -> R3 {
+ auto R3_exact = [](const R1& x) -> R3 {
return R3{+2.3 + 1.7 * x[0], //
+1.4 - 0.6 * x[0], //
-0.2 + 3.1 * x[0]};
@@ -122,31 +277,14 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
DiscreteFunctionP0<R3> uh{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_exact(xj[cell_id]); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
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_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}));
-
- 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));
- }
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
}
}
@@ -161,7 +299,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
auto& mesh = *p_mesh;
- auto R3x3_affine = [](const R1& x) -> R3x3 {
+ auto R3x3_exact = [](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],
@@ -173,94 +311,46 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
DiscreteFunctionP0<R3x3> Ah{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_exact(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));
- }
-
- {
- 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(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-13));
- }
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R3x3_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
}
}
SECTION("R vector data")
{
- using R3 = TinyVector<3>;
-
for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
SECTION(named_mesh.name())
{
auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
auto& mesh = *p_mesh;
+ std::array<std::function<double(const R1&)>, 3> vector_exact =
+ {[](const R1& x) -> double { return +2.3 + 1.7 * x[0]; },
+ [](const R1& x) -> double { return -1.7 + 2.1 * x[0]; },
+ [](const R1& x) -> double { return +1.4 - 0.6 * 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();
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];
+ for (size_t i = 0; i < vector_exact.size(); ++i) {
+ Vh[cell_id][i] = vector_exact[i](xj[cell_id]);
}
});
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+ auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
- 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};
-
- 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]));
- }
- }
- REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
- }
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
}
}
@@ -275,13 +365,13 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
auto& mesh = *p_mesh;
- auto vector_affine0 = [](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 vector_affine1 = [](const R1& x) -> R3 {
- return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]};
- };
+ std::array<std::function<R3(const R1&)>, 2> vector_exact =
+ {[](const R1& x) -> R3 {
+ return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]};
+ },
+ [](const R1& x) -> R3 {
+ return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]};
+ }};
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
@@ -289,54 +379,16 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- Vh[cell_id][0] = vector_affine0(xj[cell_id]);
- Vh[cell_id][1] = vector_affine1(xj[cell_id]);
+ Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
+ Vh[cell_id][1] = vector_exact[1](xj[cell_id]);
});
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<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_Vh(cell_id, 0)(xj[cell_id]) - vector_affine0(xj[cell_id])));
- max_mean_error =
- std::max(max_mean_error, l2Norm(dpk_Vh(cell_id, 1)(xj[cell_id]) - vector_affine1(xj[cell_id])));
- }
- REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- }
-
- {
- double max_slope_error = 0;
- {
- const TinyVector<3> slope0{+1.7, +2.1, -0.6};
-
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- for (size_t i = 0; i < R3::Dimension; ++i) {
- const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R1{0.1} + xj[cell_id])[i] -
- dpk_Vh(cell_id, 0)(xj[cell_id] - R1{0.1})[i]);
-
- max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope0[i]));
- }
- }
- }
-
- {
- const TinyVector<3> slope1{+0.7, +1.2, -0.3};
-
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- for (size_t i = 0; i < R3::Dimension; ++i) {
- const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R1{0.1} + xj[cell_id])[i] -
- dpk_Vh(cell_id, 1)(xj[cell_id] - R1{0.1})[i]);
-
- max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope1[i]));
- }
- }
- }
- REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
- }
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
}
}
@@ -352,15 +404,15 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
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 R_exact = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
- auto R3_affine = [](const R1& x) -> R3 {
+ auto R3_exact = [](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 {
+ auto R3x3_exact = [](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],
@@ -368,30 +420,30 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
};
};
- 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]};
- };
+ std::array<std::function<double(const R1&)>, 3> vector_exact =
+ {[](const R1& x) -> double { return +2.3 + 1.7 * x[0]; },
+ [](const R1& x) -> double { return -1.7 + 2.1 * x[0]; },
+ [](const R1& x) -> double { return +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]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_exact(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]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_exact(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]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_exact(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];
+ for (size_t i = 0; i < vector_exact.size(); ++i) {
+ Vh[cell_id][i] = vector_exact[i](xj[cell_id]);
}
});
@@ -400,103 +452,28 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
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;
-
- 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>>();
-
- {
- 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}));
-
- 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_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
- auto dpk_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));
+ auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
{
- 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(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-13));
- }
-
- 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));
+ auto dpk_Ah = reconstructions[2]->get<DiscreteFunctionDPk<1, const R3x3>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R3x3_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
{
- 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}));
-
- 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));
+ auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
}
}
@@ -515,48 +492,20 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
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_exact = [](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] = R_affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_exact(xj[cell_id]); });
auto reconstructions = PolynomialReconstruction{descriptor}.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_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;
-
- 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;
-
- 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));
- }
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
}
}
@@ -571,7 +520,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
auto& mesh = *p_mesh;
- auto R3_affine = [](const R2& x) -> R3 {
+ auto R3_exact = [](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]};
@@ -581,42 +530,14 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
DiscreteFunctionP0<R3> uh{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_exact(xj[cell_id]); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
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}));
-
- 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}));
-
- 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));
- }
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
}
}
@@ -631,7 +552,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<2>>();
auto& mesh = *p_mesh;
- auto R2x2_affine = [](const R2& x) -> R2x2 {
+ auto R2x2_exact = [](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]};
};
@@ -640,118 +561,47 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
DiscreteFunctionP0<R2x2> Ah{p_mesh};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_exact(xj[cell_id]); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
- 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_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}));
-
- 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_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));
- }
+ auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R2x2_exact));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
}
}
SECTION("vector data")
{
- 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]};
- };
+ std::array<std::function<double(const R2&)>, 4> vector_exact =
+ {[](const R2& x) -> double { return +2.3 + 1.7 * x[0] + 1.2 * x[1]; },
+ [](const R2& x) -> double { return -1.7 + 2.1 * x[0] - 2.2 * x[1]; },
+ [](const R2& x) -> double { return +1.4 - 0.6 * x[0] - 2.1 * x[1]; },
+ [](const R2& x) -> double { return +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];
+ for (size_t i = 0; i < vector_exact.size(); ++i) {
+ Vh[cell_id][i] = vector_exact[i](xj[cell_id]);
}
});
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));
- }
-
- {
- 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]));
- }
- }
- REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
- }
+ auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>();
+ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
}
}
}