diff --git a/src/algebra/Givens.hpp b/src/algebra/Givens.hpp
index 648086e1ea5942dc00ca5b80051ec20c093a89e5..dc69a453764659cffa085810846338fc401647a9 100644
--- a/src/algebra/Givens.hpp
+++ b/src/algebra/Givens.hpp
@@ -40,10 +40,10 @@ class Givens
static void
_lift(const MatrixType& A, const RHSVectorType& b, VectorType& x)
{
- for (ssize_t i = x.size() - 1; i >= 0; --i) {
+ for (ssize_t i = x.dimension() - 1; i >= 0; --i) {
const double inv_Aii = 1 / A(i, i);
double sum = 0;
- for (size_t j = i + 1; j < x.size(); ++j) {
+ for (size_t j = i + 1; j < x.dimension(); ++j) {
sum += A(i, j) * x[j];
}
x[i] = (b[i] - sum) * inv_Aii;
@@ -71,12 +71,13 @@ class Givens
static void
solveInPlace(MatrixType& A, VectorType& x, RHSVectorType& b)
{
- Assert(x.size() == A.numberOfColumns(), "The number of columns of A must be the size of x");
- Assert(b.size() == A.numberOfRows(), "The number of rows of A must be the size of b");
+ Assert(x.dimension() == A.numberOfColumns(), "The number of columns of A must be the size of x");
+ Assert(b.dimension() == A.numberOfRows(), "The number of rows of A must be the size of b");
for (size_t j = 0; j < A.numberOfColumns(); ++j) {
for (size_t i = A.numberOfRows() - 1; i > j; --i) {
- double c(0), s(0);
+ double c;
+ double s;
_givens(A(i - 1, j), A(i, j), c, s);
for (size_t k = j; k < A.numberOfColumns(); ++k) {
double xi_1(0), xi(0);
@@ -98,8 +99,8 @@ class Givens
static void
solve(const MatrixType& A, VectorType& x, const RHSVectorType& b)
{
- Assert(x.size() == A.numberOfColumns(), "The number of columns of A must be the size of x");
- Assert(b.size() == A.numberOfRows(), "The number of rows of A must be the size of b");
+ Assert(x.dimension() == A.numberOfColumns(), "The number of columns of A must be the size of x");
+ Assert(b.dimension() == A.numberOfRows(), "The number of rows of A must be the size of b");
MatrixType rotateA = copy(A);
RHSVectorType rotateb = copy(b);
@@ -109,7 +110,7 @@ class Givens
template <typename MatrixType, typename UnknownMatrixType, typename RHSMatrixType>
static void
- solveCollectionInPlace(const MatrixType& A, UnknownMatrixType& X, const RHSMatrixType& B)
+ solveCollectionInPlace(MatrixType& A, UnknownMatrixType& X, RHSMatrixType& B)
{
Assert(X.numberOfRows() == A.numberOfColumns(), "The number of columns of A must be the number of rows of X");
Assert(B.numberOfRows() == A.numberOfRows(), "The number of rows of A must be the rows of B");
@@ -117,7 +118,8 @@ class Givens
for (size_t j = 0; j < A.numberOfColumns(); ++j) {
for (size_t i = A.numberOfRows() - 1; i > j; --i) {
- double c(0), s(0);
+ double c;
+ double s;
_givens(A(i - 1, j), A(i, j), c, s);
for (size_t k = j; k < A.numberOfColumns(); ++k) {
double xi_1(0), xi(0);
diff --git a/src/algebra/SmallVector.hpp b/src/algebra/SmallVector.hpp
index d4a0bf84268bdf1d164cf6216153117adafaab56..bc3686c75e227c738129c532f4ce4683b5953e49 100644
--- a/src/algebra/SmallVector.hpp
+++ b/src/algebra/SmallVector.hpp
@@ -166,6 +166,13 @@ class SmallVector // LCOV_EXCL_LINE
return m_values.size();
}
+ PUGS_INLINE
+ size_t
+ dimension() const noexcept
+ {
+ return m_values.size();
+ }
+
PUGS_INLINE SmallVector&
fill(const DataType& value) noexcept
{
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index dadbe73bb24ad84e491248b89016baf85f0dbc31..c093c34f1608e91fda3c2a8ba4011e21648b4381 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -10,6 +10,8 @@ add_library(
DiscreteFunctionVectorIntegrator.cpp
DiscreteFunctionVectorInterpoler.cpp
FluxingAdvectionSolver.cpp
+
+ test_reconstruction.cpp
)
target_link_libraries(
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index 7010d3d196b012b2f39a434c0c81111844f4b608..90b95c1db6eed6aeae60873969045a0f9846ebf8 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -20,55 +20,149 @@ class PolynomialReconstruction
{
public:
template <MeshConcept MeshType, typename DataType>
- PUGS_INLINE DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>>
+ [[nodiscard]] PUGS_INLINE DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>>
build(const MeshType& mesh, const DiscreteFunctionP0<DataType> p0_function)
{
+ using Rd = TinyVector<MeshType::Dimension>;
+
const size_t degree = 1;
const Stencil stencil = StencilBuilder{}.build(mesh.connectivity());
- auto xr = mesh.xr();
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
auto cell_is_owned = mesh.connectivity().cellIsOwned();
DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>> dpk{mesh.meshVariant(), degree};
- if constexpr ((is_polygonal_mesh_v<MeshType>)and(MeshType::Dimension == 1)) {
- auto qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(1));
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
- if (cell_is_owned[cell_j_id]) {
- auto stencil_cell_list = stencil[cell_j_id];
- SmallVector<double> b(stencil_cell_list.size());
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
+ if (cell_is_owned[cell_j_id]) {
+ auto stencil_cell_list = stencil[cell_j_id];
+ SmallVector<double> b(stencil_cell_list.size());
- const double p_j = p0_function[cell_j_id];
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- const CellId cell_i_id = stencil_cell_list[i];
+ const double pj = p0_function[cell_j_id];
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
- b[i] = p0_function[cell_i_id] - p_j;
- }
+ b[i] = p0_function[cell_i_id] - pj;
+ }
- SmallMatrix<double> A(stencil_cell_list.size(), degree);
+ SmallMatrix<double> A(stencil_cell_list.size(), MeshType::Dimension);
- using R1 = TinyVector<1>;
+ const Rd& Xj = xj[cell_j_id];
- const R1& Xj = xj[cell_j_id];
- auto X_Xj = [&](const R1& X) { return X - Xj; };
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const Rd Xi_Xj = xj[cell_i_id] - Xj;
+ for (size_t l = 0; l < MeshType::Dimension; ++l) {
+ A(i, l) = Xi_Xj[l];
+ }
+ }
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- const CellId cell_i_id = stencil_cell_list[i];
+ SmallVector<double> x(MeshType::Dimension);
+ Givens::solveInPlace(A, x, b);
- A(i, 0) = X_Xj(xj[cell_i_id])[0];
- }
+ auto dpk_j = dpk.coefficients(cell_j_id);
- SmallVector<double> x(1);
- Givens::solveInPlace(A, x, b);
+ dpk_j[0] = pj;
- dpk.coefficients(cell_j_id)[0] = p_j;
- dpk.coefficients(cell_j_id)[1] = x[0];
- }
- });
+ for (size_t l = 0; l < MeshType::Dimension; ++l) {
+ dpk_j[1 + l] = x[l];
+ }
+ }
+ });
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
+ if (cell_is_owned[cell_j_id]) {
+ auto stencil_cell_list = stencil[cell_j_id];
+ SmallMatrix<double> B(stencil_cell_list.size(), DataType::Dimension);
+
+ const DataType& pj = p0_function[cell_j_id];
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const DataType& pi_pj = p0_function[cell_i_id] - pj;
+ for (size_t k = 0; k < DataType::Dimension; ++k) {
+ B(i, k) = pi_pj[k];
+ }
+ }
+
+ SmallMatrix<double> A(stencil_cell_list.size(), MeshType::Dimension);
+
+ const Rd& Xj = xj[cell_j_id];
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const Rd Xi_Xj = xj[cell_i_id] - Xj;
+ for (size_t l = 0; l < MeshType::Dimension; ++l) {
+ A(i, l) = Xi_Xj[l];
+ }
+ }
+
+ SmallMatrix<double> X(MeshType::Dimension, DataType::Dimension);
+ Givens::solveCollectionInPlace(A, X, B);
+
+ auto dpk_j = dpk.coefficients(cell_j_id);
+
+ dpk_j[0] = pj;
+ 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, k);
+ }
+ }
+ }
+ });
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
+ if (cell_is_owned[cell_j_id]) {
+ auto stencil_cell_list = stencil[cell_j_id];
+ SmallMatrix<double> B(stencil_cell_list.size(), DataType::Dimension);
+
+ const DataType& pj = p0_function[cell_j_id];
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const DataType& pi_pj = p0_function[cell_i_id] - pj;
+ for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+ for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+ B(i, k * DataType::NumberOfColumns + l) = pi_pj(k, l);
+ }
+ }
+ }
+
+ SmallMatrix<double> A(stencil_cell_list.size(), MeshType::Dimension);
+
+ const Rd& Xj = xj[cell_j_id];
+
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const Rd Xi_Xj = xj[cell_i_id] - Xj;
+ for (size_t l = 0; l < MeshType::Dimension; ++l) {
+ A(i, l) = Xi_Xj[l];
+ }
+ }
+
+ SmallMatrix<double> X(MeshType::Dimension, DataType::Dimension);
+ Givens::solveCollectionInPlace(A, X, B);
+
+ auto dpk_j = dpk.coefficients(cell_j_id);
+ dpk_j[0] = pj;
+
+ 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, k * DataType::NumberOfColumns + l);
+ }
+ }
+ }
+ }
+ });
+ } else {
+ throw NotImplementedError("dealing with " + demangle<DataType>());
+ }
}
synchronize(dpk.cellArrays());
diff --git a/src/utils/PugsTraits.hpp b/src/utils/PugsTraits.hpp
index 2e1d409578c9295a1ad032466842c2f336ba272f..ab2244b69095f5bba89dc3d823590524fc7a49e5 100644
--- a/src/utils/PugsTraits.hpp
+++ b/src/utils/PugsTraits.hpp
@@ -106,6 +106,12 @@ inline constexpr bool is_tiny_vector_v = false;
template <size_t N, typename T>
inline constexpr bool is_tiny_vector_v<TinyVector<N, T>> = true;
+template <typename T>
+inline constexpr bool is_tiny_vector_v<const T> = is_tiny_vector_v<std::remove_cvref_t<T>>;
+
+template <typename T>
+inline constexpr bool is_tiny_vector_v<T&> = is_tiny_vector_v<std::remove_cvref_t<T>>;
+
// Traits is_tiny_matrix
template <typename T>
@@ -114,6 +120,12 @@ inline constexpr bool is_tiny_matrix_v = false;
template <size_t M, size_t N, typename T>
inline constexpr bool is_tiny_matrix_v<TinyMatrix<M, N, T>> = true;
+template <typename T>
+inline constexpr bool is_tiny_matrix_v<const T> = is_tiny_matrix_v<std::remove_cvref_t<T>>;
+
+template <typename T>
+inline constexpr bool is_tiny_matrix_v<T&> = is_tiny_matrix_v<std::remove_cvref_t<T>>;
+
// Trais is ItemValue
template <typename T>
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index 9a1df325788e1277b4e0996c0316bb0fe9081ce1..72036d09fd9995cbf8ef8a1635dfbb7fa0c17983 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -24,40 +24,529 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
using R1 = TinyVector<1>;
- for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
- SECTION(named_mesh.name())
- {
- auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
- auto& mesh = *p_mesh;
+ SECTION("R data")
+ {
+ for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ SECTION(named_mesh.name())
+ {
+ auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
+
+ auto 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]); });
+
+ 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, std::abs(dpk[cell_id](xj[cell_id]) - 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[cell_id](R1{0.1} + xj[cell_id]) - dpk[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));
+ }
+ }
+ }
+ }
+
+ SECTION("R^d 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 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};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
- auto affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
- auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+ PolynomialReconstruction reconstruction;
+ auto dpk = reconstruction.build(mesh, fh);
- DiscreteFunctionP0<double> fh{p_mesh};
+ {
+ 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])));
+ }
+ 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[cell_id](R1{0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R1{0.1}));
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+ 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));
+ }
+ }
+ }
+ }
- PolynomialReconstruction reconstruction;
- auto dpk = reconstruction.build(mesh, fh);
+ SECTION("R^mn data")
+ {
+ using R32 = TinyMatrix<3, 2>;
+ for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
+ SECTION(named_mesh.name())
{
- 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])));
+ auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
+
+ auto affine = [](const R1& x) -> R32 {
+ return R32{+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]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0<R32> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+ 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, frobeniusNorm(dpk[cell_id](xj[cell_id]) - 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 R32 reconstructed_slope =
+ (1 / 0.2) * (dpk[cell_id](R1{0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R1{0.1}));
+
+ max_slope_error = std::max(max_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.7, +2.1, //
+ -0.6, -2.3, //
+ +3.1, -3.6}));
+ }
+ REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13));
}
- REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
}
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
+ SECTION("R data")
+ {
+ for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
+ SECTION(named_mesh.name())
{
- 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;
+ 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();
+
+ DiscreteFunctionP0<double> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+ 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, std::abs(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ }
+ REQUIRE(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[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk[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(max_x_slope_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 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;
+
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
+ }
+ REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
+ }
+ }
+ }
+ }
+
+ SECTION("R^d data")
+ {
+ using R3 = TinyVector<3>;
+
+ 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 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};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+ 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[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ }
+ REQUIRE(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[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk[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(max_x_slope_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 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}));
+
+ 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-14));
+ }
+ }
+ }
+ }
+
+ SECTION("R^mn data")
+ {
+ using R32 = TinyMatrix<3, 2>;
+
+ 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 affine = [](const R2& x) -> R32 {
+ return R32{+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], //
+ -0.2 + 3.1 * x[0] + 0.8 * x[1], -3.2 - 3.6 * x[0] - 1.7 * x[1]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0<R32> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+ 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, frobeniusNorm(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ }
+ REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ double max_x_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R32 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}));
+
+ max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.7, 2.1, //
+ -0.6, -2.3, //
+ +3.1, -3.6}));
+ }
+ REQUIRE(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 R32 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}));
+
+ max_y_slope_error = std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.2, -2.2, //
+ -2.1, 1.3, //
+ +0.8, -1.7}));
+ }
+ REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
+ }
+ }
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ using R3 = TinyVector<3>;
+
+ SECTION("R data")
+ {
+ 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) { 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]); });
+
+ 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, std::abs(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ }
+ REQUIRE(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[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk[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(max_x_slope_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 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;
+
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
+ }
+ REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
+ }
+
+ {
+ 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;
+
+ max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - 2.1));
+ }
+ REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-14));
+ }
+ }
+ }
+ }
+
+ SECTION("R^d 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], //
+ +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]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0<R4> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+ 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[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ }
+ REQUIRE(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 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}));
+
+ max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R4{1.7, -0.6, 3.1, -1.5}));
+ }
+ REQUIRE(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 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}));
+
+ max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R4{-2.2, 1.3, -1.1, -0.1}));
+ }
+ REQUIRE(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 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}));
+
+ max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R4{1.8, -3.7, 1.9, -1.6}));
+ }
+ REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-13));
+ }
+ }
+ }
+ }
+
+ SECTION("R^mn data")
+ {
+ using R32 = TinyMatrix<3, 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 affine = [](const R3& x) -> R32 {
+ return R32{+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], //
+ +1.4 - 0.6 * x[0] - 2.1 * x[1] + 2.9 * 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], -3.2 - 3.6 * x[0] - 1.7 * x[1] + 0.7 * x[2]};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0<R32> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
+
+ 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, frobeniusNorm(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
+ }
+ REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ double max_x_slope_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R32 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}));
+
+ max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.7, 2.1, //
+ -0.6, -2.3, //
+ +3.1, -3.6}));
+ }
+ REQUIRE(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 R32 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}));
+
+ max_y_slope_error = std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.2, -2.2, //
+ -2.1, 1.3, //
+ +0.8, -1.7}));
+ }
+ REQUIRE(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 R32 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}));
- max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
+ max_z_slope_error = std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R32{-1.3, -2.4, //
+ +2.9, +1.4, //
+ -1.8, +0.7}));
+ }
+ REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12));
}
- REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
}
}
}