diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 4fa7cc1d3e3f669eff35969405869a8173ae9721..69d24b9e588587e83aab48ef194e8ed20883fca8 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -7,7 +7,11 @@
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp>
+#include <geometry/CubeTransformation.hpp>
+#include <geometry/PrismTransformation.hpp>
+#include <geometry/PyramidTransformation.hpp>
#include <geometry/SquareTransformation.hpp>
+#include <geometry/TetrahedronTransformation.hpp>
#include <geometry/TriangleTransformation.hpp>
class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
@@ -220,6 +224,16 @@ PolynomialReconstruction::_build(
mean_i_of_ejk_pool[i] = SmallArray<double>(basis_dimension - 1);
}
+#warning remove after rework
+ const bool FORCE_INTEGRATE = []() {
+ auto value = std::getenv("INTEGRATE");
+ if (value == nullptr) {
+ return false;
+ } else {
+ return true;
+ }
+ }();
+
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
if (cell_is_owned[cell_j_id]) {
@@ -286,17 +300,7 @@ PolynomialReconstruction::_build(
ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
-#warning remove after rework
- std::string ENV_SWITCH = []() {
- auto value = std::getenv("INTEGRATE");
- if (value == nullptr) {
- return std::string{""};
- } else {
- return std::string{value};
- }
- }();
-
- if ((m_descriptor.degree() == 1) and ((MeshType::Dimension == 3) or (ENV_SWITCH == ""))) {
+ if ((m_descriptor.degree() == 1) and (not FORCE_INTEGRATE)) {
const Rd& Xj = xj[cell_j_id];
for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
const CellId cell_i_id = stencil_cell_list[i];
@@ -352,7 +356,7 @@ PolynomialReconstruction::_build(
compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
} else {
- throw NotImplementedError("unexpected cell type");
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_j_id])});
}
for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
@@ -370,7 +374,7 @@ PolynomialReconstruction::_build(
mean_i_of_ejk);
} else {
- throw NotImplementedError("unexpected cell type");
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_i_id])});
}
for (size_t l = 0; l < basis_dimension - 1; ++l) {
@@ -443,7 +447,7 @@ PolynomialReconstruction::_build(
compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
} else {
- throw NotImplementedError("unexpected cell type");
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_j_id])});
}
for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
@@ -470,13 +474,171 @@ PolynomialReconstruction::_build(
mean_i_of_ejk);
} else {
- throw NotImplementedError("unexpected cell type");
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_i_id])});
+ }
+
+ for (size_t l = 0; l < basis_dimension - 1; ++l) {
+ A(i, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+ }
+ }
+
+ } else if constexpr (MeshType::Dimension == 3) {
+ auto compute_mean_ejk = [&inv_Vj_wq_detJ_ek](const auto& quadrature, const auto& T, const Rd& Xj,
+ const double Vi, const size_t degree, const size_t dimension,
+ SmallArray<double>& mean_of_ejk) {
+ mean_of_ejk.fill(0);
+
+ for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+ const double wq = quadrature.weight(i_q);
+ const Rd& xi_q = quadrature.point(i_q);
+ const double detT = [&] {
+ if constexpr (std::is_same_v<TetrahedronTransformation, std::decay_t<decltype(T)>>) {
+ return T.jacobianDeterminant();
+ } else {
+ return T.jacobianDeterminant(xi_q);
+ }
+ }();
+
+ const Rd X_Xj = T(xi_q) - Xj;
+
+ const double x_xj = X_Xj[0];
+ const double y_yj = X_Xj[1];
+ const double z_zj = X_Xj[2];
+
+ {
+ size_t k = 0;
+ inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
+ for (; k <= degree; ++k) {
+ inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
+ }
+
+ for (size_t i_y = 1; i_y <= degree; ++i_y) {
+ const size_t begin_i_y_1 = ((i_y - 1) * (2 * degree - i_y + 4)) / 2;
+ for (size_t l = 0; l <= degree - i_y; ++l) {
+ inv_Vj_wq_detJ_ek[k++] = y_yj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
+ }
+ }
+
+ for (size_t i_z = 1; i_z <= degree; ++i_z) {
+ const size_t begin_i_z_1 = i_z * (i_z + 1) * (i_z + 2) / 6 - 1;
+ for (size_t i_y = 0; i_y <= degree - i_z; ++i_y) {
+ const size_t begin_i_y_1 = (i_y * (2 * degree - i_y + 3)) / 2 + begin_i_z_1;
+ for (size_t l = 0; l <= degree - i_y - i_z; ++l) {
+ inv_Vj_wq_detJ_ek[k++] = z_zj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
+ }
+ }
+ }
+ }
+
+ for (size_t k = 1; k < dimension; ++k) {
+ mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
+ }
+ }
+ };
+
+ const Rd& Xj = xj[cell_j_id];
+
+ if (cell_type[cell_j_id] == CellType::Tetrahedron) {
+ auto cell_node = cell_to_node_matrix[cell_j_id];
+
+ const TetrahedronTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], xr[cell_node[3]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{m_descriptor.degree()});
+
+ compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+
+ } else if (cell_type[cell_j_id] == CellType::Prism) {
+ auto cell_node = cell_to_node_matrix[cell_j_id];
+
+ const PrismTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], //
+ xr[cell_node[3]], xr[cell_node[4]], xr[cell_node[5]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{m_descriptor.degree() + 1});
+
+ compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+
+ } else if (cell_type[cell_j_id] == CellType::Pyramid) {
+ auto cell_node = cell_to_node_matrix[cell_j_id];
+
+ const PyramidTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], xr[cell_node[3]],
+ xr[cell_node[4]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{m_descriptor.degree() + 1});
+
+ compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+
+ } else if (cell_type[cell_j_id] == CellType::Hexahedron) {
+ auto cell_node = cell_to_node_matrix[cell_j_id];
+
+ const CubeTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], xr[cell_node[3]],
+ xr[cell_node[4]], xr[cell_node[5]], xr[cell_node[6]], xr[cell_node[7]]};
+
+ const auto& quadrature = QuadratureManager::instance().getCubeFormula(
+ GaussLegendreQuadratureDescriptor{m_descriptor.degree() + 1});
+
+ compute_mean_ejk(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension, mean_j_of_ejk);
+
+ } else {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_j_id])});
+ }
+
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+
+ auto cell_node = cell_to_node_matrix[cell_i_id];
+
+ if (cell_type[cell_i_id] == CellType::Tetrahedron) {
+ const TetrahedronTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]],
+ xr[cell_node[3]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{m_descriptor.degree()});
+
+ compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ mean_i_of_ejk);
+
+ } else if (cell_type[cell_i_id] == CellType::Prism) {
+ const PrismTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], //
+ xr[cell_node[3]], xr[cell_node[4]], xr[cell_node[5]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{m_descriptor.degree() + 1});
+
+ compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ mean_i_of_ejk);
+
+ } else if (cell_type[cell_i_id] == CellType::Pyramid) {
+ const PyramidTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], xr[cell_node[3]],
+ xr[cell_node[4]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{m_descriptor.degree() + 1});
+
+ compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ mean_i_of_ejk);
+
+ } else if (cell_type[cell_i_id] == CellType::Hexahedron) {
+ const CubeTransformation T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], xr[cell_node[3]],
+ xr[cell_node[4]], xr[cell_node[5]], xr[cell_node[6]], xr[cell_node[7]]};
+
+ const auto& quadrature = QuadratureManager::instance().getCubeFormula(
+ GaussLegendreQuadratureDescriptor{m_descriptor.degree() + 1});
+
+ compute_mean_ejk(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
+ mean_i_of_ejk);
+
+ } else {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_i_id])});
}
for (size_t l = 0; l < basis_dimension - 1; ++l) {
A(i, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
}
}
+
} else {
throw NotImplementedError("unexpected dimension");
}
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index 0e9398fab9398ef286ccb7e3932cd6ca4fa31bf8..47a538f0d1f4e32d8367616739b878ddbe728cb7 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -692,7 +692,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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-14));
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
@@ -703,7 +703,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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));
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
@@ -714,7 +714,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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-14));
+ REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -761,7 +761,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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));
+ REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-12));
}
{
@@ -772,7 +772,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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));
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
}
{
@@ -783,7 +783,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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-13));
+ REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -847,7 +847,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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-13));
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
}
{
@@ -860,7 +860,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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-13));
+ REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -934,7 +934,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
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));
+ REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
}
{