diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index f9e55f19c2073b4e814ee3740603ae4d924f2322..c20fbc897c060fa182eb30683b109160938d68c6 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -1,18 +1,26 @@
#include <scheme/PolynomialReconstruction.hpp>
-#include <scheme/DiscreteFunctionUtils.hpp>
-#include <scheme/DiscreteFunctionVariant.hpp>
-
+#include <algebra/Givens.hpp>
+#include <algebra/ShrinkMatrixView.hpp>
+#include <algebra/ShrinkVectorView.hpp>
+#include <algebra/SmallMatrix.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/CubeTransformation.hpp>
+#include <geometry/LineTransformation.hpp>
#include <geometry/PrismTransformation.hpp>
#include <geometry/PyramidTransformation.hpp>
#include <geometry/SquareTransformation.hpp>
#include <geometry/TetrahedronTransformation.hpp>
#include <geometry/TriangleTransformation.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/StencilManager.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
{
@@ -261,6 +269,126 @@ _computeEjkMean(const QuadratureFormula<3>& quadrature,
}
}
+template <typename NodeListT, size_t Dimension>
+void
+_computeEjkMean(const TinyVector<Dimension>& Xj,
+ const NodeValue<const TinyVector<Dimension>>& xr,
+ const NodeListT& node_list,
+ const CellType& cell_type,
+ const double Vi,
+ const size_t degree,
+ const size_t basis_dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk)
+{}
+
+template <typename NodeListT>
+void
+_computeEjkMean(const TinyVector<1>& Xj,
+ const NodeValue<const TinyVector<1>>& xr,
+ const NodeListT& node_list,
+ const CellType& cell_type,
+ const double Vi,
+ const size_t degree,
+ const size_t basis_dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ if (cell_type == CellType::Line) {
+ const LineTransformation<1> T{xr[node_list[0]], xr[node_list[1]]};
+
+ const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{degree});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+
+ } else {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+}
+
+template <typename NodeListT>
+void
+_computeEjkMean(const TinyVector<2>& Xj,
+ const NodeValue<const TinyVector<2>>& xr,
+ const NodeListT& node_list,
+ const CellType& cell_type,
+ const double Vi,
+ const size_t degree,
+ const size_t basis_dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ switch (cell_type) {
+ case CellType::Triangle: {
+ const TriangleTransformation<2> T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]]};
+ const auto& quadrature = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{degree});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ break;
+ }
+ case CellType::Quadrangle: {
+ const SquareTransformation<2> T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]], xr[node_list[3]]};
+ const auto& quadrature =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor{degree + 1});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ break;
+ }
+ default: {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+ }
+}
+
+template <typename NodeListT>
+void
+_computeEjkMean(const TinyVector<3>& Xj,
+ const NodeValue<const TinyVector<3>>& xr,
+ const NodeListT& node_list,
+ const CellType& cell_type,
+ const double Vi,
+ const size_t degree,
+ const size_t basis_dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ if (cell_type == CellType::Tetrahedron) {
+ const TetrahedronTransformation T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]], xr[node_list[3]]};
+
+ const auto& quadrature = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{degree});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+
+ } else if (cell_type == CellType::Prism) {
+ const PrismTransformation T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]], //
+ xr[node_list[3]], xr[node_list[4]], xr[node_list[5]]};
+
+ const auto& quadrature = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{degree + 1});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+
+ } else if (cell_type == CellType::Pyramid) {
+ const PyramidTransformation T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]], xr[node_list[3]],
+ xr[node_list[4]]};
+
+ const auto& quadrature = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{degree + 1});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+
+ } else if (cell_type == CellType::Hexahedron) {
+ const CubeTransformation T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]], xr[node_list[3]],
+ xr[node_list[4]], xr[node_list[5]], xr[node_list[6]], xr[node_list[7]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{degree + 1});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+
+ } else {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+}
+
size_t
PolynomialReconstruction::_getNumberOfColumns(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
@@ -374,7 +502,7 @@ PolynomialReconstruction::_build(
SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency());
- SmallArray<SmallArray<double>> G_pool(Kokkos::DefaultExecutionSpace::concurrency());
+ SmallArray<SmallVector<double>> G_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallMatrix<double>> X_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallArray<double>> inv_Vj_wq_detJ_ek_pool(Kokkos::DefaultExecutionSpace::concurrency());
@@ -384,7 +512,7 @@ PolynomialReconstruction::_build(
for (size_t i = 0; i < A_pool.size(); ++i) {
A_pool[i] = SmallMatrix<double>(max_stencil_size, basis_dimension - 1);
B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns);
- G_pool[i] = SmallArray<double>(basis_dimension - 1);
+ G_pool[i] = SmallVector<double>(basis_dimension - 1);
X_pool[i] = SmallMatrix<double>(basis_dimension - 1, number_of_columns);
inv_Vj_wq_detJ_ek_pool[i] = SmallArray<double>(basis_dimension);
@@ -482,217 +610,20 @@ PolynomialReconstruction::_build(
SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
SmallArray<double>& mean_i_of_ejk = mean_i_of_ejk_pool[t];
- if constexpr (MeshType::Dimension == 1) {
- const Rd& Xj = xj[cell_j_id];
-
- if (cell_type[cell_j_id] == CellType::Line) {
- auto cell_node = cell_to_node_matrix[cell_j_id];
-
- const LineTransformation<1> T{xr[cell_node[0]], xr[cell_node[1]]};
-
- const auto& quadrature =
- QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{m_descriptor.degree()});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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::Line) {
- const LineTransformation<1> T{xr[cell_node[0]], xr[cell_node[1]]};
-
- const auto& quadrature =
- QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor{m_descriptor.degree()});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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 if constexpr (MeshType::Dimension == 2) {
- const Rd& Xj = xj[cell_j_id];
-
- if (cell_type[cell_j_id] == CellType::Triangle) {
- auto cell_node = cell_to_node_matrix[cell_j_id];
-
- const TriangleTransformation<2> T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]]};
-
- const auto& quadrature =
- QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{m_descriptor.degree()});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, mean_j_of_ejk);
-
- } else if (cell_type[cell_j_id] == CellType::Quadrangle) {
- auto cell_node = cell_to_node_matrix[cell_j_id];
-
- const SquareTransformation<2> T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], xr[cell_node[3]]};
-
- const auto& quadrature = QuadratureManager::instance().getSquareFormula(
- GaussLegendreQuadratureDescriptor{m_descriptor.degree() + 1});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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::Triangle) {
- const TriangleTransformation<2> T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]]};
-
- const auto& quadrature =
- QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{m_descriptor.degree()});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, mean_i_of_ejk);
-
- } else if (cell_type[cell_i_id] == CellType::Quadrangle) {
- const SquareTransformation<2> T{xr[cell_node[0]], xr[cell_node[1]], xr[cell_node[2]], xr[cell_node[3]]};
-
- const auto& quadrature = QuadratureManager::instance().getSquareFormula(
- GaussLegendreQuadratureDescriptor{m_descriptor.degree() + 1});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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 if constexpr (MeshType::Dimension == 3) {
- 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()});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_j_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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()});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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});
-
- _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, 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 Rd& Xj = xj[cell_j_id];
- const auto& quadrature = QuadratureManager::instance().getCubeFormula(
- GaussLegendreQuadratureDescriptor{m_descriptor.degree() + 1});
+ _computeEjkMean(Xj, xr, cell_to_node_matrix[cell_j_id], cell_type[cell_j_id], Vj[cell_j_id],
+ m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek, mean_j_of_ejk);
- _computeEjkMean(quadrature, T, Xj, Vj[cell_i_id], m_descriptor.degree(), basis_dimension,
- inv_Vj_wq_detJ_ek, mean_i_of_ejk);
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
- } else {
- throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type[cell_i_id])});
- }
+ _computeEjkMean(Xj, xr, cell_to_node_matrix[cell_i_id], cell_type[cell_i_id], Vj[cell_i_id],
+ m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek, mean_i_of_ejk);
- for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(i, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
- }
+ 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");
}
}
@@ -717,7 +648,7 @@ PolynomialReconstruction::_build(
if (m_descriptor.preconditioning()) {
// Add column weighting preconditioning (increase the presition)
- SmallArray<double>& G = G_pool[t];
+ SmallVector<double>& G = G_pool[t];
for (size_t l = 0; l < basis_dimension - 1; ++l) {
double g = 0;
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index f523148728712bdf6cb29af45686d3db9b0a0bae..b2b855271d1ed4140ebb51d1b21314602caa0dc6 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -1,24 +1,9 @@
#ifndef POLYNOMIAL_RECONSTRUCTION_HPP
#define POLYNOMIAL_RECONSTRUCTION_HPP
-#include <algebra/ShrinkMatrixView.hpp>
-#include <algebra/ShrinkVectorView.hpp>
-#include <algebra/SmallMatrix.hpp>
-#include <algebra/SmallVector.hpp>
#include <mesh/MeshTraits.hpp>
-#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/PolynomialReconstructionDescriptor.hpp>
-#warning MOVE TO .cpp WHEN FINISHED
-#include <algebra/Givens.hpp>
-#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
-#include <analysis/QuadratureFormula.hpp>
-#include <analysis/QuadratureManager.hpp>
-#include <geometry/LineTransformation.hpp>
-#include <mesh/MeshData.hpp>
-#include <mesh/MeshDataManager.hpp>
-#include <mesh/StencilManager.hpp>
-
class DiscreteFunctionDPkVariant;
class DiscreteFunctionVariant;
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index 47a538f0d1f4e32d8367616739b878ddbe728cb7..b7c7019f57e6227de20f746900d87416e0a400ff 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -11,6 +11,7 @@
#include <algebra/SmallVector.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshDataManager.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/PolynomialReconstruction.hpp>