diff --git a/tests/test_PolynomialReconstruction_degree_1.cpp b/tests/test_PolynomialReconstruction_degree_1.cpp
index 7c4f07266514eaf8b277af5f26751c2ec17bc8ae..6c333913e85fa811abe34af119176d3d7209f156 100644
--- a/tests/test_PolynomialReconstruction_degree_1.cpp
+++ b/tests/test_PolynomialReconstruction_degree_1.cpp
@@ -1196,24 +1196,24 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
}
}
- std::vector<PolynomialReconstructionDescriptor> descriptor_list =
- {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1,
- std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
- std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
- PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1,
- std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
- std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
- PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, 1,
- std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
- std::make_shared<NamedBoundaryDescriptor>("XMIN")}}};
-
- for (auto descriptor : descriptor_list) {
- SECTION(name(descriptor.integrationMethodType()))
- {
- SECTION("1D")
+ SECTION("1D")
+ {
+ std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1,
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+ std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1,
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+ std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, 1,
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
+ std::make_shared<NamedBoundaryDescriptor>("XMIN")}}};
+
+ using R1 = TinyVector<1>;
+
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
{
- using R1 = TinyVector<1>;
-
SECTION("R^1 data")
{
auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
@@ -1323,588 +1323,381 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
}
}
}
+ }
+ }
- SECTION("2D")
+ SECTION("2D")
+ {
+ std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1,
+ std::vector<std::shared_ptr<
+ const IBoundaryDescriptor>>{std::
+ make_shared<NamedBoundaryDescriptor>(
+ "XMAX"),
+ std::make_shared<NamedBoundaryDescriptor>(
+ "YMAX")}},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1,
+ std::vector<std::shared_ptr<
+ const IBoundaryDescriptor>>{std::
+ make_shared<NamedBoundaryDescriptor>(
+ "XMAX"),
+ std::make_shared<NamedBoundaryDescriptor>(
+ "YMAX")}},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, 1,
+ std::vector<std::shared_ptr<
+ const IBoundaryDescriptor>>{std::
+ make_shared<NamedBoundaryDescriptor>(
+ "XMAX"),
+ std::make_shared<NamedBoundaryDescriptor>(
+ "YMAX")}}};
+
+ using R2 = TinyVector<2>;
+
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
{
-#warning not done
- using R2 = TinyVector<2>;
+ SECTION("R^2 data")
+ {
+ auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+ auto& mesh = *p_mesh;
+
+ auto R2_affine = [](const R2& x) -> R2 { return R2{2.3 * (x[0] - 2), -1.3 * (x[1] - 1)}; };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0<R2> fh{p_mesh};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R2_affine(xj[cell_id]); });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
- // SECTION("R data")
- // {
- // 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 R_affine = [](const R2& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1]; };
- // auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- // DiscreteFunctionP0<double> fh{p_mesh};
-
- // parallel_for(
- // mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(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));
- // }
- // }
- // }
- // }
-
- // SECTION("R^3 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 R3_affine = [](const R2& x) -> R3 {
- // return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1], //
- // +1.4 - 0.6 * x[0] + 1.3 * x[1], //
- // -0.2 + 3.1 * x[0] - 1.1 * x[1]};
- // };
- // auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- // DiscreteFunctionP0<R3> uh{p_mesh};
-
- // parallel_for(
- // mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(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));
- // }
- // }
- // }
- // }
-
- // SECTION("R^2x2 data")
- // {
- // using R2x2 = TinyMatrix<2, 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 R2x2_affine = [](const R2& x) -> R2x2 {
- // return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1], //
- // +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]};
- // };
- // auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- // DiscreteFunctionP0<R2x2> Ah{p_mesh};
-
- // parallel_for(
- // mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(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));
- // }
- // }
- // }
- // }
-
- // 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{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_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2>>();
+
+ {
+ 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_fh[cell_id](xj[cell_id]) - R2_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 R2 reconstructed_slope =
+ 1. / 0.2 * (dpk_fh[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0.1, 0}));
+
+ max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R2{2.3, 0}));
+ }
+ 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 R2 reconstructed_slope =
+ 1 / 0.2 * (dpk_fh[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0, 0.1}));
+
+ max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - (R2{0, -1.3})));
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
+ }
+ }
+
+ SECTION("vector of R2")
+ {
+ using R4 = TinyVector<4>;
+
+ auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+ auto& mesh = *p_mesh;
+
+ auto vector_affine = [](const R2& x) -> R4 {
+ return R4{+1.7 * (x[0] - 2), //
+ -0.6 * (x[1] - 1), //
+ -2.3 * (x[0] - 2), //
+ +1.1 * (x[1] - 1)};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0Vector<R2> Vh{p_mesh, 2};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+ Vh[cell_id][0][0] = vector[0];
+ Vh[cell_id][0][1] = vector[1];
+ Vh[cell_id][1][0] = vector[2];
+ Vh[cell_id][1][1] = vector[3];
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[0] - vector[0]));
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[1] - vector[1]));
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[0] - vector[2]));
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[1] - vector[3]));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
+
+ {
+ double max_x_slope_error = 0;
+ const R4 slope{+1.7, 0, -2.3, 0};
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R2 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R2{0.1, 0} + xj[cell_id]) -
+ dpk_Vh(cell_id, 0)(xj[cell_id] - R2{0.1, 0}));
+
+ const R2 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R2{0.1, 0} + xj[cell_id]) -
+ dpk_Vh(cell_id, 1)(xj[cell_id] - R2{0.1, 0}));
+
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[0] - slope[2]));
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[1] - slope[3]));
+ }
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ double max_y_slope_error = 0;
+ const R4 slope{0, -0.6, 0, 1.1};
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R2 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R2{0, 0.1} + xj[cell_id]) -
+ dpk_Vh(cell_id, 0)(xj[cell_id] - R2{0, 0.1}));
+
+ const R2 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R2{0, 0.1} + xj[cell_id]) -
+ dpk_Vh(cell_id, 1)(xj[cell_id] - R2{0, 0.1}));
+
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[0] - slope[2]));
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[1] - slope[3]));
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
}
+ }
+ }
- SECTION("3D")
+ SECTION("3D")
+ {
+ std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, 1,
+ std::vector<std::shared_ptr<
+ const IBoundaryDescriptor>>{std::
+ make_shared<NamedBoundaryDescriptor>(
+ "XMAX"),
+ std::make_shared<NamedBoundaryDescriptor>(
+ "YMAX"),
+ std::make_shared<NamedBoundaryDescriptor>(
+ "ZMAX")}},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::element, 1,
+ std::vector<std::shared_ptr<
+ const IBoundaryDescriptor>>{std::
+ make_shared<NamedBoundaryDescriptor>(
+ "XMAX"),
+ std::make_shared<NamedBoundaryDescriptor>(
+ "YMAX"),
+ std::make_shared<NamedBoundaryDescriptor>(
+ "ZMAX")}},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, 1,
+ std::vector<std::shared_ptr<
+ const IBoundaryDescriptor>>{std::
+ make_shared<NamedBoundaryDescriptor>(
+ "XMAX"),
+ std::make_shared<NamedBoundaryDescriptor>(
+ "YMAX"),
+ std::make_shared<NamedBoundaryDescriptor>(
+ "ZMAX")}}};
+
+ using R3 = TinyVector<3>;
+
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
{
- using R3 = TinyVector<3>;
-#warning not done
-
- // 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 R_affine = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; };
- // auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- // DiscreteFunctionP0<double> fh{p_mesh};
-
- // parallel_for(
- // mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); });
-
- // auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
-
- // auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
-
- // {
- // double max_mean_error = 0;
- // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- // max_mean_error =
- // std::max(max_mean_error, std::abs(dpk_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](R3{0.1, 0, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0.1, 0, 0})) /
- // 0.2;
-
- // max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
- // }
- // 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](R3{0, 0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0.1, 0})) /
- // 0.2;
-
- // max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
- // }
- // REQUIRE(parallel::allReduceMax(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 double reconstructed_slope =
- // (dpk_fh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0, 0.1})) /
- // 0.2;
-
- // max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - 2.1));
- // }// REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
- // }
- // }
- // }
- // }
-
- // SECTION("R^3 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 R3_affine = [](const R3& x) -> R3 {
- // return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2], //
- // +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2], //
- // -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]};
- // };
- // auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- // DiscreteFunctionP0<R3> uh{p_mesh};
-
- // parallel_for(
- // mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); });
-
- // auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
-
- // auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
-
- // {
- // double max_mean_error = 0;
- // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- // max_mean_error =
- // std::max(max_mean_error, l2Norm(dpk_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](R3{0.1, 0, 0} + xj[cell_id]) -
- // dpk_uh[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
-
- // max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7,
- // -0.6, 3.1}));
- // }
- // REQUIRE(parallel::allReduceMax(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 R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
- // dpk_uh[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
-
- // max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3,
- // -1.1}));
- // }
- // REQUIRE(parallel::allReduceMax(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 R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
- // dpk_uh[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
-
- // max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R3{1.8,
- // -3.7, 1.9}));
- // }
- // REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
- // }
- // }
- // }
- // }
-
- // SECTION("R^2x2 data")
- // {
- // using R2x2 = TinyMatrix<2, 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 R2x2_affine = [](const R3& x) -> R2x2 {
- // return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 *
- // x[2],
- // //
- // +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 *
- // x[2]};
- // };
- // auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
-
- // DiscreteFunctionP0<R2x2> Ah{p_mesh};
-
- // parallel_for(
- // mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]);
- // });
-
- // descriptor.setRowWeighting(false);
- // auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
-
- // auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
-
- // {
- // double max_mean_error = 0;
- // for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- // max_mean_error =
- // std::max(max_mean_error, frobeniusNorm(dpk_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](R3{0.1, 0, 0} + xj[cell_id]) -
- // dpk_Ah[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
-
- // max_x_slope_error =
- // std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, 2.1, //
- // -2.3, +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 R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
- // dpk_Ah[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
-
- // max_y_slope_error =
- // std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, //
- // 1.3, +0.8}));
- // }
- // REQUIRE(parallel::allReduceMax(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 R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
- // dpk_Ah[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
-
- // max_z_slope_error =
- // std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R2x2{-1.3, -2.4, //
- // +1.4, -1.8}));
- // }
- // REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
- // }
- // }
- // }
- // }
-
- // 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];
- // }
- // });
-
- // descriptor.setPreconditioning(false);
- // auto reconstructions = PolynomialReconstruction{descriptor}.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(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
- // }
-
- // {
- // 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(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, 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(parallel::allReduceMax(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(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-13));
- // }
- // }
- // }
- // }
+ SECTION("R^3 data")
+ {
+ auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+ auto& mesh = *p_mesh;
+
+ auto R3_affine = [](const R3& x) -> R3 {
+ return R3{2.3 * (x[0] - 2), -1.3 * (x[1] - 1), 1.4 * (x[2] - 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] = R3_affine(xj[cell_id]); });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ max_mean_error =
+ std::max(max_mean_error, l2Norm(dpk_fh[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_fh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
+
+ max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{2.3, 0, 0}));
+ }
+ 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_fh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
+
+ max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - (R3{0, -1.3, 0})));
+ }
+ REQUIRE(parallel::allReduceMax(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 R3 reconstructed_slope =
+ 1 / 0.2 *
+ (dpk_fh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
+
+ max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - (R3{0, 0, 1.4})));
+ }
+ REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
+ }
+ }
+
+ SECTION("vector of R2")
+ {
+ using R6 = TinyVector<6>;
+
+ auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+ auto& mesh = *p_mesh;
+
+ auto vector_affine = [](const R3& x) -> R6 {
+ return R6{+1.7 * (x[0] - 2), //
+ -0.6 * (x[1] - 1), //
+ +1.2 * (x[2] - 1), //
+ -2.3 * (x[0] - 2), //
+ +1.1 * (x[1] - 1), //
+ -0.3 * (x[2] - 1)};
+ };
+ auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+
+ DiscreteFunctionP0Vector<R3> Vh{p_mesh, 2};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+ Vh[cell_id][0][0] = vector[0];
+ Vh[cell_id][0][1] = vector[1];
+ Vh[cell_id][0][2] = vector[2];
+ Vh[cell_id][1][0] = vector[3];
+ Vh[cell_id][1][1] = vector[4];
+ Vh[cell_id][1][2] = vector[5];
+ });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
+
+ auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3>>();
+
+ {
+ double max_mean_error = 0;
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ auto vector = vector_affine(xj[cell_id]);
+
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[0] - vector[0]));
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[1] - vector[1]));
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[2] - vector[2]));
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[0] - vector[3]));
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[1] - vector[4]));
+ max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[2] - vector[5]));
+ }
+ REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
+ }
+
+ {
+ double max_x_slope_error = 0;
+ const R6 slope{+1.7, 0, 0, -2.3, 0, 0};
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R3{0.1, 0, 0} + xj[cell_id]) -
+ dpk_Vh(cell_id, 0)(xj[cell_id] - R3{0.1, 0, 0}));
+
+ const R3 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R3{0.1, 0, 0} + xj[cell_id]) -
+ dpk_Vh(cell_id, 1)(xj[cell_id] - R3{0.1, 0, 0}));
+
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[2] - slope[2]));
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[0] - slope[3]));
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[1] - slope[4]));
+ max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[2] - slope[5]));
+ }
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ double max_y_slope_error = 0;
+ const R6 slope{0, -0.6, 0, 0, 1.1, 0};
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R3{0, 0.1, 0} + xj[cell_id]) -
+ dpk_Vh(cell_id, 0)(xj[cell_id] - R3{0, 0.1, 0}));
+
+ const R3 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R3{0, 0.1, 0} + xj[cell_id]) -
+ dpk_Vh(cell_id, 1)(xj[cell_id] - R3{0, 0.1, 0}));
+
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[2] - slope[2]));
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[0] - slope[3]));
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[1] - slope[4]));
+ max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[2] - slope[5]));
+ }
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
+
+ {
+ double max_z_slope_error = 0;
+ const R6 slope{0, 0, 1.2, 0, 0, -0.3};
+ for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+ const R3 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R3{0, 0, 0.1} + xj[cell_id]) -
+ dpk_Vh(cell_id, 0)(xj[cell_id] - R3{0, 0, 0.1}));
+
+ const R3 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R3{0, 0, 0.1} + xj[cell_id]) -
+ dpk_Vh(cell_id, 1)(xj[cell_id] - R3{0, 0, 0.1}));
+
+ max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
+ max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
+ max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope0[2] - slope[2]));
+ max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope1[0] - slope[3]));
+ max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope1[1] - slope[4]));
+ max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope1[2] - slope[5]));
+ }
+ REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-13));
+ }
+ }
}
}
}