diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 1a2b1e072734154516b01001cb9831bfaff8ea73..5101d5da271894ed113fc0f57f5b19effe5f29ac 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -1452,7 +1452,7 @@ PolynomialReconstruction::_build(
if (m_descriptor.degree() > 1) {
auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
for (size_t i = 0; i < basis_dimension - 1; ++i) {
- auto& dpk_j_0 = dpk_j[0];
+ auto& dpk_j_0 = dpk_j[component_begin];
for (size_t k = 0; k < DataType::Dimension; ++k) {
dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
}
diff --git a/tests/test_PolynomialReconstruction_degree_2.cpp b/tests/test_PolynomialReconstruction_degree_2.cpp
index dce627cabab012e398a271b88aafc5131e4338d5..435084af00385cda388588fb21860fe32680c73e 100644
--- a/tests/test_PolynomialReconstruction_degree_2.cpp
+++ b/tests/test_PolynomialReconstruction_degree_2.cpp
@@ -38,6 +38,17 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
{
using R1 = TinyVector<1>;
+ auto p0 = [](const R1& x) { return +2.3 + 1.7 * x[0] - 2.3 * x[0] * x[0]; };
+ auto p1 = [](const R1& x) { return -1.7 + 2.1 * x[0] + 1.2 * x[0] * x[0]; };
+ auto p2 = [](const R1& x) { return +1.4 - 0.6 * x[0] - 2.0 * x[0] * x[0]; };
+ auto p3 = [](const R1& x) { return +2.4 - 2.3 * x[0] + 1.1 * x[0] * x[0]; };
+ auto p4 = [](const R1& x) { return -0.2 + 3.1 * x[0] - 0.7 * x[0] * x[0]; };
+ auto p5 = [](const R1& x) { return -3.2 - 3.6 * x[0] + 0.1 * x[0] * x[0]; };
+ auto p6 = [](const R1& x) { return -4.1 + 3.1 * x[0] - 0.2 * x[0] * x[0]; };
+ auto p7 = [](const R1& x) { return +0.8 + 2.9 * x[0] + 4.1 * x[0] * x[0]; };
+ auto p8 = [](const R1& x) { return -1.6 + 2.3 * x[0] - 1.7 * x[0] * x[0]; };
+ auto p9 = [](const R1& x) { return +2.3 + 1.7 * x[0] - 1.4 * x[0] * x[0]; };
+
SECTION("R data")
{
for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
@@ -46,7 +57,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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] - 1.4 * x[0] * x[0]; };
+ auto R_exact = p0;
DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));
@@ -70,11 +81,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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] - x[0] * x[0], //
- +1.4 - 0.6 * x[0] + 2 * x[0] * x[0], //
- -0.2 + 3.1 * x[0] + 1.4 * x[0] * x[0]};
- };
+ auto R3_exact = [&](const R1& x) -> R3 { return R3{p2(x), p4(x), p1(x)}; };
DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
@@ -98,18 +105,10 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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] - 2.3 * x[0] * x[0], //
- -1.7 + 2.1 * x[0] + 1.2 * x[0] * x[0], //
- +1.4 - 0.6 * x[0] - 2.0 * x[0] * x[0], //
- //
- +2.4 - 2.3 * x[0] + 1.1 * x[0] * x[0], //
- -0.2 + 3.1 * x[0] - 0.7 * x[0] * x[0], //
- -3.2 - 3.6 * x[0] + 0.1 * x[0] * x[0], //
- //
- -4.1 + 3.1 * x[0] - 0.2 * x[0] * x[0], //
- +0.8 + 2.9 * x[0] + 4.1 * x[0] * x[0], //
- -1.6 + 2.3 * x[0] - 1.7 * x[0] * x[0]};
+ auto R3x3_exact = [&](const R1& x) -> R3x3 {
+ return R3x3{p1(x), p2(x), p3(x), //
+ p4(x), p5(x), p6(x), //
+ p7(x), p8(x), p9(x)};
};
DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R3x3_exact));
@@ -131,10 +130,8 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
{
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] + 1.2 * x[0] * x[0]; },
- [](const R1& x) -> double { return -1.7 + 2.1 * x[0] + 2.1 * x[0] * x[0]; },
- [](const R1& x) -> double { return +1.4 - 0.6 * x[0] - 1.3 * x[0] * x[0]; }};
+
+ std::array<std::function<double(const R1&)>, 3> vector_exact = {p1, p7, p9};
DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
@@ -157,13 +154,16 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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] + 0.8 * x[0] * x[0], //
- -1.7 + 2.1 * x[0] - 0.7 * x[0] * x[0], //
- +1.4 - 0.6 * x[0] + 1.9 * x[0] * 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]}; }};
+ std::array<std::function<R3(const R1&)>, 3> vector_exact //
+ = {[&](const R1& x) -> R3 {
+ return R3{p1(x), p2(x), p3(x)};
+ },
+ [&](const R1& x) -> R3 {
+ return R3{p5(x), p7(x), p0(x)};
+ },
+ [&](const R1& x) -> R3 {
+ return R3{p9(x), p8(x), p4(x)};
+ }};
DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
@@ -188,34 +188,17 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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 R_exact = p0;
- auto R3_exact = [](const R1& x) -> R3 {
- return R3{+2.3 + 1.7 * x[0] + 2.3 * x[0] * x[0], //
- +1.4 - 0.6 * x[0] - 1.5 * x[0] * x[0], //
- -0.2 + 3.1 * x[0] + 2.9 * x[0] * x[0]};
- };
+ auto R3_exact = [&](const R1& x) -> R3 { return R3{p9(x), p4(x), p7(x)}; };
- auto R3x3_exact = [](const R1& x) -> R3x3 {
- return R3x3{
- +2.3 + 1.7 * x[0] - 0.3 * x[0] * x[0],
- -1.7 + 2.1 * x[0] + 1.4 * x[0] * x[0],
- +1.4 - 0.6 * x[0] - 2.3 * x[0] * x[0],
- //
- +2.4 - 2.3 * x[0] + 1.8 * x[0] * x[0],
- -0.2 + 3.1 * x[0] - 1.7 * x[0] * x[0],
- -3.2 - 3.6 * x[0] + 0.7 * x[0] * x[0],
- //
- -4.1 + 3.1 * x[0] - 1.9 * x[0] * x[0],
- +0.8 + 2.9 * x[0] + 2.2 * x[0] * x[0],
- -1.6 + 2.3 * x[0] - 1.3 * x[0] * x[0],
- };
+ auto R3x3_exact = [&](const R1& x) -> R3x3 {
+ return R3x3{p2(x), p1(x), p0(x), //
+ p3(x), p2(x), p4(x), //
+ p6(x), p5(x), p9(x)};
};
- std::array<std::function<double(const R1&)>, 3> vector_exact =
- {[](const R1& x) -> double { return +2.3 + 1.7 * x[0] + 2.2 * x[0] * x[0]; },
- [](const R1& x) -> double { return -1.7 + 2.1 * x[0] - 1.9 * x[0] * x[0]; },
- [](const R1& x) -> double { return +1.4 - 0.6 * x[0] + 3.1 * x[0] * x[0]; }};
+ std::array<std::function<double(const R1&)>, 3> vector_exact = {p1, p8, p7};
DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));
DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
@@ -258,6 +241,19 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
SECTION("2D")
{
using R2 = TinyVector<2>;
+ auto p0 = [](const R2& x) {
+ return +2.3 + 1.7 * x[0] - 1.3 * x[1] + 1.2 * x[0] * x[0] + 1.3 * x[0] * x[1] - 3.2 * x[1] * x[1];
+ };
+
+ auto p1 = [](const R2& x) {
+ return +2.3 + 1.7 * x[0] - 2.2 * x[1] - 2.1 * x[0] * x[0] - 2.3 * x[0] * x[1] - 3.2 * x[1] * x[1];
+ };
+ auto p2 = [](const R2& x) {
+ return +1.4 - 0.6 * x[0] + 1.3 * x[1] + 2.3 * x[0] * x[0] - 1.3 * x[0] * x[1] + 1.2 * x[1] * x[1];
+ };
+ auto p3 = [](const R2& x) {
+ return -0.2 + 3.1 * x[0] - 1.1 * x[1] - 2.1 * x[0] * x[0] + 1.3 * x[0] * x[1] - 1.1 * x[1] * x[1];
+ };
SECTION("R data")
{
@@ -267,10 +263,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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] //
- + 1.2 * x[0] * x[0] + 1.3 * x[0] * x[1] - 3.2 * x[1] * x[1];
- };
+ auto R_exact = p0;
DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));
@@ -294,14 +287,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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] //
- - 2.1 * x[0] * x[0] - 2.3 * x[0] * x[1] - 3.2 * x[1] * x[1], //
- +1.4 - 0.6 * x[0] + 1.3 * x[1] //
- + 2.3 * x[0] * x[0] - 1.3 * x[0] * x[1] + 1.2 * x[1] * x[1], //
- -0.2 + 3.1 * x[0] - 1.1 * x[1] //
- - 2.1 * x[0] * x[0] + 1.3 * x[0] * x[1] - 1.1 * x[1] * x[1]};
- };
+ auto R3_exact = [&](const R2& x) -> R3 { return R3{p1(x), p2(x), p3(x)}; };
DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
@@ -325,15 +311,9 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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] //
- - 2.1 * x[0] * x[0] + 1.3 * x[0] * x[1] + 1.2 * x[1] * x[1], //
- -1.7 + 2.1 * x[0] - 2.2 * x[1] //
- - 1.2 * x[0] * x[0] + 2.1 * x[0] * x[1] - 1.3 * x[1] * x[1], //
- +1.4 - 0.6 * x[0] - 2.1 * x[1] //
- - 1.1 * x[0] * x[0] - 2.3 * x[0] * x[1] + 2.1 * x[1] * x[1],
- +2.4 - 2.3 * x[0] + 1.3 * x[1] //
- + 2.7 * x[0] * x[0] + 2.1 * x[0] * x[1] - 2.7 * x[1] * x[1]};
+ auto R2x2_exact = [&](const R2& x) -> R2x2 {
+ return R2x2{p0(x), p1(x), //
+ p2(x), p3(x)};
};
DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
@@ -355,19 +335,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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] + 1.2 * x[0] * x[0] - 2.1 * x[0] * x[1] + 3.1 * x[1] * x[1];
- },
- [](const R2& x) -> double {
- return -1.7 + 2.1 * x[0] - 2.2 * x[1] - 0.7 * x[0] * x[0] + 2.2 * x[0] * x[1] - 1.6 * x[1] * x[1];
- },
- [](const R2& x) -> double {
- return +1.4 - 0.6 * x[0] - 2.1 * x[1] + 2.3 * x[0] * x[0] + 2.3 * x[0] * x[1] - 2.9 * x[1] * x[1];
- },
- [](const R2& x) -> double {
- return +2.4 - 2.3 * x[0] + 1.3 * x[1] - 2.7 * x[0] * x[0] - 1.2 * x[0] * x[1] - 0.7 * x[1] * x[1];
- }};
+ std::array<std::function<double(const R2&)>, 4> vector_exact = {p0, p1, p2, p3};
DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
@@ -385,6 +353,30 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
{
using R3 = TinyVector<3>;
+ auto p0 = [](const R3& x) {
+ return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2] //
+ + 1.7 * x[0] * x[0] + 1.4 * x[1] * x[1] + 1.7 * x[2] * x[2] //
+ - 2.3 * x[0] * x[1] + 1.6 * x[0] * x[2] - 1.9 * x[1] * x[2];
+ };
+
+ auto p1 = [](const R3& x) {
+ return +2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2] //
+ + 1.7 * x[0] * x[0] - 2.4 * x[1] * x[1] - 2.3 * x[2] * x[2] //
+ - 2.1 * x[0] * x[1] + 2.6 * x[0] * x[2] + 1.6 * x[1] * x[2];
+ };
+
+ auto p2 = [](const R3& x) {
+ return +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2] //
+ + 3.1 * x[0] * x[0] - 1.1 * x[1] * x[1] + 1.7 * x[2] * x[2] //
+ - 2.3 * x[0] * x[1] - 2.6 * x[0] * x[2] - 1.9 * x[1] * x[2];
+ };
+
+ auto p3 = [](const R3& x) {
+ return -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2] //
+ - 1.5 * x[0] * x[0] + 1.4 * x[1] * x[1] - 1.2 * x[2] * x[2] //
+ - 1.7 * x[0] * x[1] - 1.3 * x[0] * x[2] + 2.1 * x[1] * x[2];
+ };
+
SECTION("R data")
{
for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
@@ -393,11 +385,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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] //
- + 1.7 * x[0] * x[0] + 1.4 * x[1] * x[1] + 1.7 * x[2] * x[2] //
- - 2.3 * x[0] * x[1] + 1.6 * x[0] * x[2] - 1.9 * x[1] * x[2];
- };
+ auto R_exact = p0;
DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R_exact));
@@ -419,17 +407,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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.7 * x[0] * x[0] - 2.4 * x[1] * x[1] - 2.3 * x[2] * x[2] //
- - 2.1 * x[0] * x[1] + 2.6 * x[0] * x[2] + 1.6 * x[1] * x[2], //
- +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2] //
- + 3.1 * x[0] * x[0] - 1.1 * x[1] * x[1] + 1.7 * x[2] * x[2] //
- - 2.3 * x[0] * x[1] - 2.6 * x[0] * x[2] - 1.9 * x[1] * x[2], //
- -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2] //
- - 1.5 * x[0] * x[0] + 1.4 * x[1] * x[1] - 1.2 * x[2] * x[2] //
- - 1.7 * x[0] * x[1] - 1.3 * x[0] * x[2] + 2.1 * x[1] * x[2]};
- };
+ auto R3_exact = [&](const R3& x) -> R3 { return R3{p1(x), p2(x), p3(x)}; };
DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree, std::function(R3_exact));
@@ -453,22 +431,9 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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 * x[0] * x[0] + 1.4 * x[1] * x[1] + 1.7 * x[2] * x[2] //
- - 1.3 * x[0] * x[1] + 1.6 * x[0] * x[2] - 1.9 * x[1] * x[2], //
- -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2] //
- + 3.7 * x[0] * x[0] + 1.3 * x[1] * x[1] + 1.6 * x[2] * x[2] //
- - 2.1 * x[0] * x[1] - 1.5 * x[0] * x[2] - 1.7 * x[1] * x[2], //
- //
- +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2] //
- - 2.1 * x[0] * x[0] + 1.7 * x[1] * x[1] + 1.8 * x[2] * x[2] //
- - 1.4 * x[0] * x[1] + 1.3 * x[0] * x[2] - 2.9 * x[1] * x[2], //
- -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2] //
- + 1.6 * x[0] * x[0] + 2.1 * x[1] * x[1] - 2.1 * x[2] * x[2] //
- - 1.1 * x[0] * x[1] - 1.3 * x[0] * x[2] + 1.6 * x[1] * x[2], //
- };
+ auto R2x2_exact = [&](const R3& x) -> R2x2 {
+ return R2x2{p0(x), p1(x), //
+ p2(x), p3(x)};
};
DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree, std::function(R2x2_exact));
@@ -490,12 +455,8 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
{
auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
auto& mesh = *p_mesh;
-#warning continue here
- 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]; }};
+
+ std::array<std::function<double(const R3&)>, 4> vector_exact = {p0, p1, p2, p3};
DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
@@ -516,6 +477,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
SECTION("with symmetries")
{
+#warning continue here
// SECTION("1D")
// {
// std::vector<PolynomialReconstructionDescriptor> descriptor_list =
@@ -529,8 +491,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
// using R1 = TinyVector<1>;
// for (auto descriptor : descriptor_list) {
- // SECTION(name(descriptor.integrationMethodType()))
- // {
+ // SECTION(name(descriptor.integrationMethodType()))// {
// SECTION("R^1 data")
// {
// auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();