From 70558d9576558ce0c303a4807da878f263e0f587 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 5 Jul 2024 19:09:07 +0200
Subject: [PATCH] Fix degree 1 reconstruction for discrete function p0 vectors
Also add tests
---
src/scheme/PolynomialReconstruction.cpp | 106 +++++----
tests/test_PolynomialReconstruction.cpp | 273 +++++++++++++++++++++++-
2 files changed, 339 insertions(+), 40 deletions(-)
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index c770f5a31..2cc86d14a 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -122,12 +122,16 @@ PolynomialReconstruction::_build(
} else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
return data_type::Dimension;
} else {
+ // LCOV_EXCL_START
throw UnexpectedError("unexpected data type " + demangle<data_type>());
+ // LCOV_EXCL_STOP
}
} else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
return discrete_function.size();
} else {
+ // LCOV_EXCL_START
throw UnexpectedError("unexpected discrete function type");
+ // LCOV_EXCL_STOP
}
},
discrete_function_variant_p->discreteFunction());
@@ -171,9 +175,11 @@ PolynomialReconstruction::_build(
} else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
mutable_discrete_function_dpk_variant_list.push_back(
- DiscreteFunctionDPkVector<MeshType::Dimension, DataType>(p_mesh, discrete_function.size(), degree));
+ DiscreteFunctionDPkVector<MeshType::Dimension, DataType>(p_mesh, degree, discrete_function.size()));
} else {
+ // LCOV_EXCL_START
throw UnexpectedError("unexpected discrete function type");
+ // LCOV_EXCL_STOP
}
},
discrete_function_variant->discreteFunction());
@@ -245,7 +251,9 @@ PolynomialReconstruction::_build(
}
column_begin += pj_vector.size();
} else {
+ // LCOV_EXCL_START
throw UnexpectedError("invalid discrete function type");
+ // LCOV_EXCL_STOP
}
},
discrete_function_variant->discreteFunction());
@@ -282,59 +290,79 @@ PolynomialReconstruction::_build(
if constexpr (std::is_same_v<DataType, P0DataType>) {
if constexpr (is_discrete_function_P0_v<P0FunctionT>) {
- auto dpk_j = dpk_function.coefficients(cell_j_id);
- dpk_j[0] = p0_function[cell_j_id];
+ if constexpr (is_discrete_function_dpk_scalar_v<DPkFunctionT>) {
+ auto dpk_j = dpk_function.coefficients(cell_j_id);
+ dpk_j[0] = p0_function[cell_j_id];
- if constexpr (std::is_arithmetic_v<DataType>) {
- for (size_t i = 0; i < MeshType::Dimension; ++i) {
- auto& dpk_j_ip1 = dpk_j[i + 1];
- dpk_j_ip1 = X(i, column_begin);
- }
- ++column_begin;
- } else if constexpr (is_tiny_vector_v<DataType>) {
- for (size_t i = 0; i < MeshType::Dimension; ++i) {
- auto& dpk_j_ip1 = dpk_j[i + 1];
- for (size_t k = 0; k < DataType::Dimension; ++k) {
- dpk_j_ip1[k] = X(i, column_begin + k);
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ for (size_t i = 0; i < MeshType::Dimension; ++i) {
+ auto& dpk_j_ip1 = dpk_j[i + 1];
+ dpk_j_ip1 = X(i, column_begin);
}
- }
- column_begin += DataType::Dimension;
- } else if constexpr (is_tiny_matrix_v<DataType>) {
- for (size_t i = 0; i < MeshType::Dimension; ++i) {
- auto& dpk_j_ip1 = dpk_j[i + 1];
- for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
+ ++column_begin;
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ for (size_t i = 0; i < MeshType::Dimension; ++i) {
+ auto& dpk_j_ip1 = dpk_j[i + 1];
+ for (size_t k = 0; k < DataType::Dimension; ++k) {
+ dpk_j_ip1[k] = X(i, column_begin + k);
+ }
+ }
+ column_begin += DataType::Dimension;
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ for (size_t i = 0; i < MeshType::Dimension; ++i) {
+ auto& dpk_j_ip1 = dpk_j[i + 1];
+ for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+ for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+ dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
+ }
}
}
+ column_begin += DataType::Dimension;
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected data type");
+ // LCOV_EXCL_STOP
}
- column_begin += DataType::Dimension;
} else {
- throw UnexpectedError("unexpected data type");
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected discrete dpk function type");
+ // LCOV_EXCL_STOP
}
} else if constexpr (is_discrete_function_P0_vector_v<P0FunctionT>) {
- auto dpk_j = dpk_function.coefficients(cell_j_id);
- auto cell_vector = p0_function[cell_j_id];
- const size_t size = MeshType::Dimension + 1;
-
- for (size_t l = 0; l < cell_vector.size(); ++l) {
- const size_t component_begin = l * size;
- dpk_j[component_begin] = cell_vector[l];
- if constexpr (std::is_arithmetic_v<DataType>) {
- for (size_t i = 0; i < MeshType::Dimension; ++i) {
- auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
- dpk_j_ip1 = X(i, column_begin);
+ if constexpr (is_discrete_function_dpk_vector_v<DPkFunctionT>) {
+ auto dpk_j = dpk_function.coefficients(cell_j_id);
+ auto cell_vector = p0_function[cell_j_id];
+ const size_t size = MeshType::Dimension + 1;
+
+ for (size_t l = 0; l < cell_vector.size(); ++l) {
+ const size_t component_begin = l * size;
+ dpk_j[component_begin] = cell_vector[l];
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ for (size_t i = 0; i < MeshType::Dimension; ++i) {
+ auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
+ dpk_j_ip1 = X(i, column_begin);
+ }
+ ++column_begin;
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected data type");
+ // LCOV_EXCL_STOP
}
- ++column_begin;
- } else {
- throw UnexpectedError("unexpected data type");
}
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected discrete dpk function type");
+ // LCOV_EXCL_STOP
}
} else {
+ // LCOV_EXCL_START
throw UnexpectedError("unexpected discrete function type");
+ // LCOV_EXCL_STOP
}
} else {
+ // LCOV_EXCL_START
throw UnexpectedError("incompatible data types");
+ // LCOV_EXCL_STOP
}
},
dpk_variant.mutableDiscreteFunctionDPk(), discrete_function_variant->discreteFunction());
@@ -364,7 +392,7 @@ PolynomialReconstruction::build(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
if (not hasSameMesh(discrete_function_variant_list)) {
- throw NormalError("cannot reconstruct functions living of different mesh simultaneously");
+ throw NormalError("cannot reconstruct functions living of different meshes simultaneously");
}
auto mesh_v = getCommonMesh(discrete_function_variant_list);
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index 123072266..029fa30e8 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -171,6 +171,64 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
}
}
+ 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{}.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(max_mean_error == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ 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(max_slope_error == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+ }
+
SECTION("list of various types")
{
using R3x3 = TinyMatrix<3>;
@@ -198,6 +256,10 @@ TEST_CASE("PolynomialReconstruction", "[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]};
+ };
+
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
DiscreteFunctionP0<double> fh{p_mesh};
@@ -212,8 +274,18 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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{}.build(std::make_shared<DiscreteFunctionVariant>(fh), uh,
- std::make_shared<DiscreteFunctionP0<R3x3>>(Ah));
+ std::make_shared<DiscreteFunctionP0<R3x3>>(Ah),
+ DiscreteFunctionVariant(Vh));
auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
@@ -283,6 +355,34 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
}
REQUIRE(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(max_mean_error == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ 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(max_slope_error == Catch::Approx(0).margin(1E-13));
+ }
+ }
}
}
}
@@ -464,6 +564,78 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
}
}
}
+
+ 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]};
+ };
+ 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{}.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(max_mean_error == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ 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(max_x_slope_error == Catch::Approx(0).margin(1E-12));
+ }
+
+ {
+ 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(max_y_slope_error == Catch::Approx(0).margin(1E-12));
+ }
+ }
+ }
+ }
}
SECTION("3D")
@@ -675,5 +847,104 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
}
}
}
+
+ SECTION("vector 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 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];
+ }
+ });
+
+ auto reconstructions = PolynomialReconstruction{}.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(max_mean_error == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ 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(max_x_slope_error == Catch::Approx(0).margin(1E-12));
+ }
+
+ {
+ 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(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]));
+ }
+ }
+ REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12));
+ }
+ }
+ }
+ }
+ }
+
+ 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};
+
+ REQUIRE_THROWS_WITH(PolynomialReconstruction{}.build(f1, f2),
+ "error: cannot reconstruct functions living of different meshes simultaneously");
}
}
--
GitLab