Skip to content
Snippets Groups Projects
Commit 1d505a11 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Add tests for PolynomialReconstruction of degree 2

parent 00d64dec
No related branches found
No related tags found
1 merge request!205High-order polynomial reconstruction
Show whitespace changes
Inline Side-by-side
tests/test_PolynomialReconstruction_degree_2.cpp
+ 562
− 196
View file @ 1d505a11
...@@ -477,201 +477,567 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]") ...@@ -477,201 +477,567 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
SECTION("with symmetries") SECTION("with symmetries")
{ {
#warning continue here SECTION("1D")
// SECTION("1D") {
// { std::vector<PolynomialReconstructionDescriptor> descriptor_list =
// std::vector<PolynomialReconstructionDescriptor> descriptor_list = {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
// {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree, std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
// std::vector<std::shared_ptr<const IBoundaryDescriptor>>{ std::make_shared<NamedBoundaryDescriptor>("XMIN")}},
// std::make_shared<NamedBoundaryDescriptor>("XMIN")}}, PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
// PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree, std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
// std::vector<std::shared_ptr<const IBoundaryDescriptor>>{ std::make_shared<NamedBoundaryDescriptor>("XMIN")}}};
// std::make_shared<NamedBoundaryDescriptor>("XMIN")}}}; using R1 = TinyVector<1>;
// using R1 = TinyVector<1>; auto p0 = [](const R1& x) { return +1.7 * (x[0] + 1) * (x[0] + 1) - 1.1; };
auto p1 = [](const R1& x) { return -1.2 * (x[0] + 1) * (x[0] + 1) + 1.3; };
// for (auto descriptor : descriptor_list) { auto p2 = [](const R1& x) { return +1.4 * (x[0] + 1) * (x[0] + 1) - 0.6; };
// SECTION(name(descriptor.integrationMethodType()))// {
// SECTION("R^1 data") for (auto descriptor : descriptor_list) {
// { SECTION(name(descriptor.integrationMethodType()))
// auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>(); {
SECTION("R data")
// auto& mesh = *p_mesh; {
auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
// auto R1_exact = [](const R1& x) { return R1{1.7 * (x[0] + 1)}; };
auto& mesh = *p_mesh;
// DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R1_exact));
auto R_exact = p0;
// auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));
// auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1>>();
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
// 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)); 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));
// SECTION("R1 vector data") REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
// { }
// for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
// SECTION(named_mesh.name()) SECTION("R1x1 data")
// { {
// auto p_mesh = named_mesh.mesh()->get<Mesh<1>>(); using R1x1 = TinyMatrix<1>;
// auto& mesh = *p_mesh;
auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
// std::array<std::function<R1(const R1&)>, 2> vector_exact //
// = {[](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1)}; }, auto& mesh = *p_mesh;
// [](const R1& x) -> R1 { return R1{-0.3 * (x[0] + 1)}; }};
auto R1x1_exact = [&](const R1& x) { return R1x1{p0(x)}; };
// DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R1x1_exact));
// auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
// auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R1>>();
auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1x1>>();
// double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
// REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R1x1_exact));
// } REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
// } }
// }
// } SECTION("R vector data")
// } {
// } auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
auto& mesh = *p_mesh;
// SECTION("2D")
// { std::array<std::function<double(const R1&)>, 3> vector_exact //
// std::vector<PolynomialReconstructionDescriptor> descriptor_list = = {p0, p1, p2};
// {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
// std::vector<std::shared_ptr< DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
// const IBoundaryDescriptor>>{std::
// make_shared<NamedBoundaryDescriptor>( auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
// "XMAX"),
// std::make_shared<NamedBoundaryDescriptor>( auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
// "YMAX")}},
// PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree, double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
// std::vector<std::shared_ptr< REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
// const IBoundaryDescriptor>>{std:: }
// make_shared<NamedBoundaryDescriptor>(
// "XMAX"), SECTION("R1x1 vector data")
// std::make_shared<NamedBoundaryDescriptor>( {
// "YMAX")}}}; using R1x1 = TinyMatrix<1>;
// using R2 = TinyVector<2>; auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
auto& mesh = *p_mesh;
// for (auto descriptor : descriptor_list) {
// SECTION(name(descriptor.integrationMethodType())) std::array<std::function<R1x1(const R1&)>, 3> vector_exact //
// { = {[&](const R1& x) { return R1x1{p2(x)}; }, //
// SECTION("R^2 data") [&](const R1& x) { return R1x1{p0(x)}; }, //
// { [&](const R1& x) { return R1x1{p1(x)}; }};
// auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
// auto& mesh = *p_mesh; DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
// auto R2_exact = [](const R2& x) -> R2 { return R2{2.3 * (x[0] - 2), -1.3 * (x[1] - 1)}; }; auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
// DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R2_exact)); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R1x1>>();
// auto reconstructions = PolynomialReconstruction{descriptor}.build(uh); double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
// 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("2D")
// SECTION("vector of R2") {
// { std::vector<PolynomialReconstructionDescriptor> descriptor_list =
// auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>(); {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
// auto& mesh = *p_mesh; std::vector<std::shared_ptr<
const IBoundaryDescriptor>>{std::
// std::array<std::function<R2(const R2&)>, 2> vector_exact // make_shared<NamedBoundaryDescriptor>(
// = {[](const R2& x) -> R2 { "XMAX"),
// return R2{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1)}; std::make_shared<NamedBoundaryDescriptor>(
// }, "YMAX")}},
// [](const R2& x) -> R2 { PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
// return R2{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1)}; std::vector<std::shared_ptr<
// }}; const IBoundaryDescriptor>>{std::
make_shared<NamedBoundaryDescriptor>(
// DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact); "XMAX"),
std::make_shared<NamedBoundaryDescriptor>(
// auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); "YMAX")}}};
// auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2>>(); using R2 = TinyVector<2>;
// double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact); auto p_initial_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
// REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); auto& initial_mesh = *p_initial_mesh;
// }
// } constexpr double theta = 1;
// } TinyMatrix<2> T{std::cos(theta), -std::sin(theta), //
// } std::sin(theta), std::cos(theta)};
// SECTION("3D") auto xr = initial_mesh.xr();
// {
// std::vector<PolynomialReconstructionDescriptor> descriptor_list = NodeValue<R2> new_xr{initial_mesh.connectivity()};
// {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree, parallel_for(
// std::vector<std::shared_ptr< initial_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { new_xr[node_id] = T * xr[node_id]; });
// const IBoundaryDescriptor>>{std::
// make_shared<NamedBoundaryDescriptor>( std::shared_ptr p_mesh = std::make_shared<const Mesh<2>>(initial_mesh.shared_connectivity(), new_xr);
// "XMAX"), const Mesh<2>& mesh = *p_mesh;
// std::make_shared<NamedBoundaryDescriptor>( // inverse rotation
// "YMAX"), TinyMatrix<2> inv_T{std::cos(theta), std::sin(theta), //
// std::make_shared<NamedBoundaryDescriptor>( -std::sin(theta), std::cos(theta)};
// "ZMAX")}},
// PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree, auto p0 = [&inv_T](const R2& X) {
// std::vector<std::shared_ptr< const R2 Y = inv_T * X;
// const IBoundaryDescriptor>>{std:: const double x = Y[0] - 2;
// make_shared<NamedBoundaryDescriptor>( const double y = Y[1] - 1;
// "XMAX"), return +1.7 * x * x + 2 * y * y - 1.1;
// std::make_shared<NamedBoundaryDescriptor>( };
// "YMAX"),
// std::make_shared<NamedBoundaryDescriptor>( auto p1 = [&inv_T](const R2& X) {
// "ZMAX")}}}; const R2 Y = inv_T * X;
const double x = Y[0] - 2;
// using R3 = TinyVector<3>; const double y = Y[1] - 1;
return -1.3 * x * x - 0.2 * y * y + 0.7;
// for (auto descriptor : descriptor_list) { };
// SECTION(name(descriptor.integrationMethodType()))
// { auto p2 = [&inv_T](const R2& X) {
// SECTION("R^3 data") const R2 Y = inv_T * X;
// { const double x = Y[0] - 2;
// auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); const double y = Y[1] - 1;
// auto& mesh = *p_mesh; return +2.6 * x * x - 1.4 * y * y - 1.9;
};
// auto R3_exact = [](const R3& x) -> R3 { return R3{2.3 * (x[0] - 2), -1.3 * (x[1] - 1), 1.4 * (x[2] - 1)};
// }; auto p3 = [&inv_T](const R2& X) {
const R2 Y = inv_T * X;
// DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact)); const double x = Y[0] - 2;
const double y = Y[1] - 1;
// auto reconstructions = PolynomialReconstruction{descriptor}.build(uh); return -0.6 * x * x + 1.4 * y * y + 2.3;
};
// auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
auto q0 = [&inv_T](const R2& X) {
// double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact)); const R2 Y = inv_T * X;
// REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); const double x = Y[0] - 2;
// } const double y = Y[1] - 1;
return 3 * x * y;
// SECTION("vector of R3") };
// {
// auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); auto q1 = [&inv_T](const R2& X) {
// auto& mesh = *p_mesh; const R2 Y = inv_T * X;
const double x = Y[0] - 2;
// std::array<std::function<R3(const R3&)>, 2> vector_exact // const double y = Y[1] - 1;
// = {[](const R3& x) -> R3 { return -1.3 * x * y;
// return R3{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1), +1.2 * (x[2] - 1)}; };
// },
// [](const R3& x) -> R3 { for (auto descriptor : descriptor_list) {
// return R3{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1), -0.3 * (x[2] - 1)}; SECTION(name(descriptor.integrationMethodType()))
// }}; {
SECTION("R data")
// DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact); {
auto R_exact = p0;
// auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
// auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3>>();
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
// double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
// REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12)); auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
// }
// } double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
// } REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
// } }
SECTION("R2x2 data")
{
using R2x2 = TinyMatrix<2>;
auto R2x2_exact = [&](const R2& X) -> R2x2 {
return T * TinyMatrix<2>{p0(X), q0(X), q1(X), p1(X)} * inv_T;
};
DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R2x2_exact));
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
SECTION("vector of R")
{
std::array<std::function<double(const R2&)>, 4> vector_exact //
= {p0, p1, p2, p3};
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-12));
}
SECTION("vector of R2x2")
{
using R2x2 = TinyMatrix<2>;
std::array<std::function<R2x2(const R2&)>, 2> vector_R2x2_exact //
= {[&](const R2& X) {
return T * R2x2{p0(X), q0(X), q1(X), p1(X)} * inv_T;
},
[&](const R2& X) {
return T * R2x2{p0(X), q1(X), 0, p1(X)} * inv_T;
}};
DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_R2x2_exact);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2x2>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_R2x2_exact);
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
SECTION("list")
{
using R2x2 = TinyMatrix<2>;
auto R_exact = p0;
DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
std::array<std::function<double(const R2&)>, 4> vector_exact //
= {p0, p1, p2, p3};
DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
std::array<std::function<R2x2(const R2&)>, 2> vector_R2x2_exact //
= {[&](const R2& X) {
return T * R2x2{p0(X), q0(X), q1(X), p1(X)} * inv_T;
},
[&](const R2& X) {
return T * R2x2{p2(X), q1(X), 0, p3(X)} * inv_T;
}};
DiscreteFunctionP0Vector Wh = test_only::exact_projection(mesh, degree, vector_R2x2_exact);
auto R2x2_exact = [&](const R2& X) -> R2x2 {
return T * TinyMatrix<2>{p0(X), q1(X), q0(X), p3(X)} * inv_T;
};
DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh, Vh, Wh, Ah);
auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>();
auto dpk_Vh = reconstructions[1]->get<DiscreteFunctionDPkVector<2, const double>>();
auto dpk_Wh = reconstructions[2]->get<DiscreteFunctionDPkVector<2, const R2x2>>();
auto dpk_Ah = reconstructions[3]->get<DiscreteFunctionDPk<2, const R2x2>>();
{
double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
{
double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
{
double max_error = test_only::max_reconstruction_error(mesh, dpk_Wh, vector_R2x2_exact);
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
{
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("3D")
{
using R3 = TinyVector<3>;
auto p_initial_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
auto& initial_mesh = *p_initial_mesh;
constexpr double theta = 1;
TinyMatrix<3> T{std::cos(theta), -std::sin(theta), 0,
//
std::sin(theta), std::cos(theta), 0,
//
0, 0, 1};
auto xr = initial_mesh.xr();
NodeValue<R3> new_xr{initial_mesh.connectivity()};
parallel_for(
initial_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { new_xr[node_id] = T * xr[node_id]; });
std::shared_ptr p_mesh = std::make_shared<const Mesh<3>>(initial_mesh.shared_connectivity(), new_xr);
const Mesh<3>& mesh = *p_mesh;
// inverse rotation
TinyMatrix<3> inv_T{std::cos(theta), std::sin(theta), 0,
//
-std::sin(theta), std::cos(theta), 0,
//
0, 0, 1};
auto p0 = [&inv_T](const R3& X) {
const R3 Y = inv_T * X;
const double x = Y[0] - 2;
const double y = Y[1] - 1;
const double z = Y[2] - 1;
return +1.7 * x * x + 2 * y * y + 1.3 * z * z - 1.1;
};
auto p1 = [&inv_T](const R3& X) {
const R3 Y = inv_T * X;
const double x = Y[0] - 2;
const double y = Y[1] - 1;
const double z = Y[2] - 1;
return +2.1 * x * x - 1.4 * y * y - 3.1 * z * z + 2.2;
};
auto p2 = [&inv_T](const R3& X) {
const R3 Y = inv_T * X;
const double x = Y[0] - 2;
const double y = Y[1] - 1;
const double z = Y[2] - 1;
return -1.1 * x * x - 1.2 * y * y + 1.3 * z * z - 1.7;
};
auto p3 = [&inv_T](const R3& X) {
const R3 Y = inv_T * X;
const double x = Y[0] - 2;
const double y = Y[1] - 1;
const double z = Y[2] - 1;
return 1.9 * x * x + 2.1 * y * y - 3.1 * z * z + 1.6;
};
auto p4 = [&inv_T](const R3& X) {
const R3 Y = inv_T * X;
const double x = Y[0] - 2;
const double y = Y[1] - 1;
const double z = Y[2] - 1;
return -2.4 * x * x + 3.3 * y * y - 1.7 * z * z + 2.1;
};
SECTION("3 symmetries")
{
std::vector<PolynomialReconstructionDescriptor> descriptor_list =
{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")}}};
for (auto descriptor : descriptor_list) {
SECTION("R data")
{
auto R_exact = p0;
DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
SECTION("vector of R")
{
std::array<std::function<double(const R3&)>, 4> vector_exact //
= {p0, p1, p2, p3};
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 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("1 symmetry")
{
// Matrix and their transformations are kept simple to
// derive exact solutions
std::vector<PolynomialReconstructionDescriptor> descriptor_list =
{PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree,
std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
std::make_shared<NamedBoundaryDescriptor>("XMAX")}},
PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree,
std::vector<std::shared_ptr<const IBoundaryDescriptor>>{
std::make_shared<NamedBoundaryDescriptor>("XMAX")}}};
for (auto descriptor : descriptor_list) {
SECTION(name(descriptor.integrationMethodType()))
{
SECTION("R3x3 data")
{
// Matrix and their transformations are kept simple to
// derive exact solutions
using R3x3 = TinyMatrix<3>;
auto R2x2_exact = [&](const R3& X) -> R3x3 {
return T * TinyMatrix<3>{p0(X), 0, 0, //
0, p1(X), p2(X), //
0, p3(X), p4(X)} *
inv_T;
};
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<3, const R3x3>>();
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 of R3x3")
{
using R3x3 = TinyMatrix<3>;
std::array<std::function<R3x3(const R3&)>, 2> vector_R3x3_exact //
= {[&](const R3& X) {
return T * R3x3{p0(X), 0, 0, //
0, p1(X), p2(X), //
0, p3(X), p4(X)} *
inv_T;
},
[&](const R3& X) {
return T * R3x3{p0(X), 0, 0, //
0, p2(X), p4(X), //
0, p3(X), p1(X)} *
inv_T;
}};
DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_R3x3_exact);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3x3>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_R3x3_exact);
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
SECTION("list")
{
using R3x3 = TinyMatrix<3>;
auto R_exact = p0;
DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R_exact));
std::array<std::function<double(const R3&)>, 4> vector_exact //
= {p0, p1, p2, p3};
DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
std::array<std::function<R3x3(const R3&)>, 2> vector_R3x3_exact //
= {[&](const R3& X) {
return T * R3x3{p1(X), 0, 0, //
0, p3(X), p0(X), //
0, p2(X), p4(X)} *
inv_T;
},
[&](const R3& X) {
return T * R3x3{p2(X), 0, 0, //
0, p0(X), p1(X), //
0, p3(X), p4(X)} *
inv_T;
}};
DiscreteFunctionP0Vector Wh = test_only::exact_projection(mesh, degree, vector_R3x3_exact);
auto R3x3_exact = [&](const R3& X) -> R3x3 {
return T * R3x3{p0(X), 0, 0, //
0, p1(X), p2(X), //
0, p3(X), p4(X)} *
inv_T;
};
DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R3x3_exact));
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh, Vh, Wh, Ah);
auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
auto dpk_Vh = reconstructions[1]->get<DiscreteFunctionDPkVector<3, const double>>();
auto dpk_Wh = reconstructions[2]->get<DiscreteFunctionDPkVector<3, const R3x3>>();
auto dpk_Ah = reconstructions[3]->get<DiscreteFunctionDPk<3, const R3x3>>();
{
double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R_exact));
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
{
double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
{
double max_error = test_only::max_reconstruction_error(mesh, dpk_Wh, vector_R3x3_exact);
REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
{
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-12));
}
}
}
}
}
}
} }
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment