#include <catch2/catch_approx.hpp> #include <catch2/catch_test_macros.hpp> #include <catch2/matchers/catch_matchers_all.hpp> #include <utils/PugsAssert.hpp> #include <mesh/Mesh.hpp> #include <mesh/NamedBoundaryDescriptor.hpp> #include <scheme/DiscreteFunctionDPkVariant.hpp> #include <scheme/DiscreteFunctionP0.hpp> #include <scheme/DiscreteFunctionVariant.hpp> #include <scheme/PolynomialReconstruction.hpp> #include <DiscreteFunctionDPkForTests.hpp> #include <MeshDataBaseForTests.hpp> // clazy:excludeall=non-pod-global-static TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") { constexpr size_t degree = 1; SECTION("without symmetries") { std::vector<PolynomialReconstructionDescriptor> descriptor_list = {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree}, PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree}, PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree}}; for (auto descriptor : descriptor_list) { SECTION(name(descriptor.integrationMethodType())) { SECTION("1D") { using R1 = TinyVector<1>; 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 R_exact = [](const R1& x) { return 2.3 + 1.7 * x[0]; }; DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(fh); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } } } SECTION("R^3 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 R3_exact = [](const R1& x) -> R3 { return R3{+2.3 + 1.7 * x[0], // +1.4 - 0.6 * x[0], // -0.2 + 3.1 * x[0]}; }; DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(uh); auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } } } SECTION("R^3x3 data") { using R3x3 = TinyMatrix<3, 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 R3x3_exact = [](const R1& x) -> R3x3 { return R3x3{ +2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0], // +2.4 - 2.3 * x[0], -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0], -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0], }; }; DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R3x3_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah); auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R3x3_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } } } SECTION("R vector 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; std::array<std::function<double(const R1&)>, 3> vector_exact = {[](const R1& x) -> double { return +2.3 + 1.7 * x[0]; }, [](const R1& x) -> double { return -1.7 + 2.1 * x[0]; }, [](const R1& x) -> double { return +1.4 - 0.6 * x[0]; }}; DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact); auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } } } SECTION("R3 vector data") { using R3 = TinyVector<3>; for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) { SECTION(named_mesh.name()) { auto p_mesh = named_mesh.mesh()->get<Mesh<1>>(); auto& mesh = *p_mesh; std::array<std::function<R3(const R1&)>, 2> vector_exact = {[](const R1& x) -> R3 { return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]}; }, [](const R1& x) -> R3 { return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]}; }}; DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact); auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R3>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } } } SECTION("list of various types") { using R3x3 = TinyMatrix<3>; 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 R_exact = [](const R1& x) { return 2.3 + 1.7 * x[0]; }; auto R3_exact = [](const R1& x) -> R3 { return R3{+2.3 + 1.7 * x[0], // +1.4 - 0.6 * x[0], // -0.2 + 3.1 * x[0]}; }; auto R3x3_exact = [](const R1& x) -> R3x3 { return R3x3{ +2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0], // +2.4 - 2.3 * x[0], -0.2 + 3.1 * x[0], -3.2 - 3.6 * x[0], -4.1 + 3.1 * x[0], +0.8 + 2.9 * x[0], -1.6 + 2.3 * x[0], }; }; std::array<std::function<double(const R1&)>, 3> vector_exact = {[](const R1& x) -> double { return +2.3 + 1.7 * x[0]; }, [](const R1& x) -> double { return -1.7 + 2.1 * x[0]; }, [](const R1& x) -> double { return +1.4 - 0.6 * x[0]; }}; DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact)); DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact)); DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R3x3_exact)); DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact); auto reconstructions = PolynomialReconstruction{descriptor}.build(std::make_shared<DiscreteFunctionVariant>(fh), uh, std::make_shared<DiscreteFunctionP0<R3x3>>(Ah), DiscreteFunctionVariant(Vh)); { auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } { auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } { auto dpk_Ah = reconstructions[2]->get<DiscreteFunctionDPk<1, const R3x3>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R3x3_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } { auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } } } } } SECTION("2D") { using R2 = TinyVector<2>; 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_exact = [](const R2& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1]; }; DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(fh); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } } } 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_exact = [](const R2& x) -> R3 { return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1], // +1.4 - 0.6 * x[0] + 1.3 * x[1], // -0.2 + 3.1 * x[0] - 1.1 * x[1]}; }; DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(uh); auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R3>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } } } 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_exact = [](const R2& x) -> R2x2 { return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1], -1.7 + 2.1 * x[0] - 2.2 * x[1], // +1.4 - 0.6 * x[0] - 2.1 * x[1], +2.4 - 2.3 * x[0] + 1.3 * x[1]}; }; DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah); auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R2x2_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } } } SECTION("vector data") { for (auto named_mesh : MeshDataBaseForTests::get().all2DMeshes()) { SECTION(named_mesh.name()) { auto p_mesh = named_mesh.mesh()->get<Mesh<2>>(); auto& mesh = *p_mesh; std::array<std::function<double(const R2&)>, 4> vector_exact = {[](const R2& x) -> double { return +2.3 + 1.7 * x[0] + 1.2 * x[1]; }, [](const R2& x) -> double { return -1.7 + 2.1 * x[0] - 2.2 * x[1]; }, [](const R2& x) -> double { return +1.4 - 0.6 * x[0] - 2.1 * x[1]; }, [](const R2& x) -> double { return +2.4 - 2.3 * x[0] + 1.3 * x[1]; }}; DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact); auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14)); } } } } 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 R_exact = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; }; DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(fh); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-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_exact = [](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]}; }; DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(uh); auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-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_exact = [](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]}; }; DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact)); descriptor.setRowWeighting(false); auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah); auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R2x2_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); } } } SECTION("vector 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; std::array<std::function<double(const R3&)>, 4> vector_exact = {[](const R3& x) -> double { return +2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2]; }, [](const R3& x) -> double { return -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2]; }, [](const R3& x) -> double { return +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2]; }, [](const R3& x) -> double { return -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]; }}; DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact); descriptor.setPreconditioning(false); auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-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{descriptor}.build(f1, f2), "error: cannot reconstruct functions living of different meshes simultaneously"); } } } } SECTION("with symmetries") { SECTION("errors") { SECTION("1D") { auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>(); PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree, std::vector<std::shared_ptr<const IBoundaryDescriptor>>{ std::make_shared<NamedBoundaryDescriptor>("XMIN")}}; REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<2>>{p_mesh}), "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 1"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<3>>{p_mesh}), "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 1"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyVector<2>>{p_mesh, 1}), "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 1"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyVector<3>>{p_mesh, 1}), "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 1"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<2>>{p_mesh}), "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 1"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<3>>{p_mesh}), "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 1"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyMatrix<2>>{p_mesh, 1}), "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 1"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyMatrix<3>>{p_mesh, 1}), "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 1"); } SECTION("2D") { auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>(); PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree, std::vector<std::shared_ptr<const IBoundaryDescriptor>>{ std::make_shared<NamedBoundaryDescriptor>("XMIN")}}; REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<1>>{p_mesh}), "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 2"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<3>>{p_mesh}), "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 2"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyVector<1>>{p_mesh, 1}), "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 2"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyVector<3>>{p_mesh, 1}), "error: cannot symmetrize vectors of dimension 3 using a mesh of dimension 2"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<1>>{p_mesh}), "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 2"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<3>>{p_mesh}), "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 2"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyMatrix<1>>{p_mesh, 1}), "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 2"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyMatrix<3>>{p_mesh, 1}), "error: cannot symmetrize matrices of dimensions 3x3 using a mesh of dimension 2"); } SECTION("3D") { auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree, std::vector<std::shared_ptr<const IBoundaryDescriptor>>{ std::make_shared<NamedBoundaryDescriptor>("XMIN")}}; REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<1>>{p_mesh}), "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 3"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyVector<2>>{p_mesh}), "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 3"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyVector<1>>{p_mesh, 1}), "error: cannot symmetrize vectors of dimension 1 using a mesh of dimension 3"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyVector<2>>{p_mesh, 1}), "error: cannot symmetrize vectors of dimension 2 using a mesh of dimension 3"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<1>>{p_mesh}), "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 3"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build(DiscreteFunctionP0<TinyMatrix<2>>{p_mesh}), "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 3"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyMatrix<1>>{p_mesh, 1}), "error: cannot symmetrize matrices of dimensions 1x1 using a mesh of dimension 3"); REQUIRE_THROWS_WITH(PolynomialReconstruction{descriptor}.build( DiscreteFunctionP0Vector<TinyMatrix<2>>{p_mesh, 1}), "error: cannot symmetrize matrices of dimensions 2x2 using a mesh of dimension 3"); } } SECTION("1D") { std::vector<PolynomialReconstructionDescriptor> descriptor_list = {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree, std::vector<std::shared_ptr<const IBoundaryDescriptor>>{ std::make_shared<NamedBoundaryDescriptor>("XMIN")}}, PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree, std::vector<std::shared_ptr<const IBoundaryDescriptor>>{ std::make_shared<NamedBoundaryDescriptor>("XMIN")}}, PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree, 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())) { SECTION("R^1 data") { auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>(); auto& mesh = *p_mesh; auto R1_exact = [](const R1& x) { return R1{1.7 * (x[0] + 1)}; }; DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R1_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(fh); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R1_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); } SECTION("R1 vector 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; std::array<std::function<R1(const R1&)>, 2> vector_exact // = {[](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1)}; }, [](const R1& x) -> R1 { return R1{-0.3 * (x[0] + 1)}; }}; DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact); auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R1>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); } } } } } } SECTION("2D") { std::vector<PolynomialReconstructionDescriptor> descriptor_list = {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree, std::vector<std::shared_ptr< const IBoundaryDescriptor>>{std:: make_shared<NamedBoundaryDescriptor>( "XMAX"), std::make_shared<NamedBoundaryDescriptor>( "YMAX")}}, PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree, std::vector<std::shared_ptr< const IBoundaryDescriptor>>{std:: make_shared<NamedBoundaryDescriptor>( "XMAX"), std::make_shared<NamedBoundaryDescriptor>( "YMAX")}}, PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree, 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())) { SECTION("R^2 data") { auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>(); auto& mesh = *p_mesh; auto R2_exact = [](const R2& x) -> R2 { return R2{2.3 * (x[0] - 2), -1.3 * (x[1] - 1)}; }; DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R2_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(uh); auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R2_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); } SECTION("vector of R2") { auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>(); auto& mesh = *p_mesh; std::array<std::function<R2(const R2&)>, 2> vector_exact // = {[](const R2& x) -> R2 { return R2{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1)}; }, [](const R2& x) -> R2 { return R2{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1)}; }}; DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact); auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); } } } } SECTION("3D") { std::vector<PolynomialReconstructionDescriptor> descriptor_list = {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree, 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, degree, 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, degree, 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())) { SECTION("R^3 data") { auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); auto& mesh = *p_mesh; auto R3_exact = [](const R3& x) -> R3 { return R3{2.3 * (x[0] - 2), -1.3 * (x[1] - 1), 1.4 * (x[2] - 1)}; }; DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact)); auto reconstructions = PolynomialReconstruction{descriptor}.build(uh); auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact)); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); } SECTION("vector of R3") { auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); auto& mesh = *p_mesh; std::array<std::function<R3(const R3&)>, 2> vector_exact // = {[](const R3& x) -> R3 { return R3{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1), +1.2 * (x[2] - 1)}; }, [](const R3& x) -> R3 { return R3{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1), -0.3 * (x[2] - 1)}; }}; DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact); auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3>>(); double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact); REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); } } } } } }