diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index b55e84fee171e21b762e4934054f6cb2b04b82e7..4c88d7476b34e92e40971a6962e649358ab5a5b6 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -323,6 +323,31 @@ class PolynomialReconstruction
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> build(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+ template <typename... DiscreteFunctionT>
+ [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
+ build(DiscreteFunctionT... input) const
+ {
+ std::vector<std::shared_ptr<const DiscreteFunctionVariant>> variant_vector;
+ auto convert = [&variant_vector](auto&& df) {
+ using DF_T = std::decay_t<decltype(df)>;
+ if constexpr (is_discrete_function_v<DF_T> or std::is_same_v<DiscreteFunctionVariant, DF_T>) {
+ variant_vector.push_back(std::make_shared<DiscreteFunctionVariant>(df));
+ } else if constexpr (is_shared_ptr_v<DF_T>) {
+ using DF_Value_T = std::decay_t<typename DF_T::element_type>;
+ if constexpr (is_discrete_function_v<DF_Value_T> or std::is_same_v<DiscreteFunctionVariant, DF_Value_T>) {
+ variant_vector.push_back(std::make_shared<DiscreteFunctionVariant>(*df));
+ } else {
+ static_assert(is_false_v<DF_T>, "unexpected type");
+ }
+ } else {
+ static_assert(is_false_v<DF_T>, "unexpected type");
+ }
+ };
+
+ (convert(std::forward<DiscreteFunctionT>(input)), ...);
+ return this->build(variant_vector);
+ }
+
PolynomialReconstruction() {}
};
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index 5242b6b041dda5429a6c3509f4b0523e4924b5ac..4099960e31d6a5cbc4dc5b10c0887a8f70178de5 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -17,31 +17,6 @@
#include <MeshDataBaseForTests.hpp>
-template <typename... DiscreteFunctionT>
-std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
-build_list(DiscreteFunctionT... input)
-{
- std::vector<std::shared_ptr<const DiscreteFunctionVariant>> variant_vector;
- auto convert = [&variant_vector](auto&& df) {
- using DF_T = std::decay_t<decltype(df)>;
- if constexpr (is_discrete_function_v<DF_T> or std::is_same_v<DiscreteFunctionVariant, DF_T>) {
- variant_vector.push_back(std::make_shared<DiscreteFunctionVariant>(df));
- } else if constexpr (is_shared_ptr_v<DF_T>) {
- using DF_Value_T = std::decay_t<typename DF_T::element_type>;
- if constexpr (is_discrete_function_v<DF_Value_T> or std::is_same_v<DiscreteFunctionVariant, DF_Value_T>) {
- variant_vector.push_back(std::make_shared<DiscreteFunctionVariant>(*df));
- } else {
- static_assert(is_false_v<DF_T>, "unexpected type");
- }
- } else {
- static_assert(is_false_v<DF_T>, "unexpected type");
- }
- };
-
- (convert(input), ...);
- return variant_vector;
-}
-
// clazy:excludeall=non-pod-global-static
TEST_CASE("PolynomialReconstruction", "[scheme]")
@@ -69,7 +44,10 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
DiscreteFunctionVariant fh_v(fh);
std::shared_ptr<const DiscreteFunctionVariant> p_fh_v = std::make_shared<DiscreteFunctionVariant>(fh);
- auto reconstructions = PolynomialReconstruction{}.build(build_list(p_fh_v, fh_v, fh));
+ std::shared_ptr<DiscreteFunctionP0<const double>> p_fh =
+ std::make_shared<DiscreteFunctionP0<const double>>(fh);
+
+ auto reconstructions = PolynomialReconstruction{}.build(p_fh, p_fh_v, fh_v, fh);
auto dpk = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>();
@@ -144,7 +122,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
SECTION("R^mn data")
{
- using R32 = TinyMatrix<3, 2>;
+ using R33 = TinyMatrix<3, 3>;
for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
SECTION(named_mesh.name())
@@ -152,14 +130,16 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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 affine = [](const R1& x) -> R33 {
+ return R33{
+ +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<R32> fh{p_mesh};
+ DiscreteFunctionP0<R33> fh{p_mesh};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
@@ -178,12 +158,15 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R32 reconstructed_slope =
+ 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}));
- max_slope_error = std::max(max_slope_error, frobeniusNorm(reconstructed_slope - R32{+1.7, +2.1, //
- -0.6, -2.3, //
- +3.1, -3.6}));
+ R33 slops = R33{+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));
}
@@ -308,7 +291,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
SECTION("R^mn data")
{
- using R32 = TinyMatrix<3, 2>;
+ using R22 = TinyMatrix<2, 2>;
for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) {
SECTION(named_mesh.name())
@@ -316,14 +299,13 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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 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 xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<R32> fh{p_mesh};
+ DiscreteFunctionP0<R22> fh{p_mesh};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
@@ -342,12 +324,11 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R32 reconstructed_slope =
+ 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}));
- 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}));
+ max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R22{+1.7, +2.1, //
+ -0.6, -2.3}));
}
REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
}
@@ -355,12 +336,11 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R32 reconstructed_slope =
+ 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}));
- 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}));
+ max_y_slope_error = std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R22{+1.2, -2.2, //
+ -2.1, +1.3}));
}
REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
}
@@ -508,7 +488,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
SECTION("R^mn data")
{
- using R32 = TinyMatrix<3, 2>;
+ using R22 = TinyMatrix<2, 2>;
for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
SECTION(named_mesh.name())
@@ -516,14 +496,14 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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 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 xj = MeshDataManager::instance().getMeshData(mesh).xj();
- DiscreteFunctionP0<R32> fh{p_mesh};
+ DiscreteFunctionP0<R22> fh{p_mesh};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
@@ -542,12 +522,11 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R32 reconstructed_slope =
+ 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}));
- 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}));
+ max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R22{+1.7, 2.1, //
+ -2.3, +3.1}));
}
REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12));
}
@@ -555,12 +534,11 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R32 reconstructed_slope =
+ 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}));
- 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}));
+ max_y_slope_error = std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R22{+1.2, -2.2, //
+ 1.3, +0.8}));
}
REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12));
}
@@ -568,12 +546,11 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
{
double max_z_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const R32 reconstructed_slope =
+ 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}));
- 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}));
+ max_z_slope_error = std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R22{-1.3, -2.4, //
+ +1.4, -1.8}));
}
REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12));
}