diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3159d9471f1d59fd8e573c4f85576156c243c3df..46126f9121cfafeb028c9d6f32a2ddf826c9fcf1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -741,6 +741,7 @@ target_link_libraries(
PugsMesh
PugsAlgebra
PugsScheme
+ PugsSchemeReconstructionUtils
PugsUtils
PugsOutput
PugsLanguageUtils
@@ -778,6 +779,7 @@ target_link_libraries(
PugsLanguageModules
PugsLanguageUtils
PugsScheme
+ PugsSchemeReconstructionUtils
PugsDev
PugsAnalysis
PugsAlgebra
@@ -820,6 +822,7 @@ install(TARGETS
PugsMesh
PugsOutput
PugsScheme
+ PugsSchemeReconstructionUtils
kokkos
Catch2
diff --git a/src/geometry/SquareTransformation.hpp b/src/geometry/SquareTransformation.hpp
index d6474f33e4ecbaa4bcfba525dadc9fd38e95962d..b0e970c6eda673b48c5cf5c0d0a1ad7e73f8ef74 100644
--- a/src/geometry/SquareTransformation.hpp
+++ b/src/geometry/SquareTransformation.hpp
@@ -84,6 +84,25 @@ class SquareTransformation
return m_coefficients * X + m_shift;
}
+ PUGS_INLINE
+ TinyVector<Dimension>
+ areaVariation(const TinyVector<2>& X) const
+ {
+ const auto& T = m_coefficients;
+ const double& x = X[0];
+ const double& y = X[1];
+
+ const TinyVector<Dimension> dxT{T(0, 0) + T(0, 2) * y, //
+ T(1, 0) + T(1, 2) * y, //
+ T(2, 0) + T(2, 2) * y};
+
+ const TinyVector<Dimension> dyT{T(0, 1) + T(0, 2) * x, //
+ T(1, 1) + T(1, 2) * x, //
+ T(2, 1) + T(2, 2) * x};
+
+ return crossProduct(dxT, dyT);
+ }
+
PUGS_INLINE double
areaVariationNorm(const TinyVector<2>& X) const
{
diff --git a/src/geometry/SymmetryUtils.hpp b/src/geometry/SymmetryUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6524a0791aef79b5440110d7180a472865c12d89
--- /dev/null
+++ b/src/geometry/SymmetryUtils.hpp
@@ -0,0 +1,32 @@
+#ifndef SYMMETRY_UTILS_HPP
+#define SYMMETRY_UTILS_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <utils/PugsMacros.hpp>
+
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimension>& u)
+{
+ return u - 2 * dot(u, normal) * normal;
+}
+
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_matrix(const TinyVector<Dimension>& normal, const TinyMatrix<Dimension>& A)
+{
+ const TinyMatrix S = TinyMatrix<Dimension>{identity} - 2 * tensorProduct(normal, normal);
+ return S * A * S;
+}
+
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_coordinates(const TinyVector<Dimension>& origin,
+ const TinyVector<Dimension>& normal,
+ const TinyVector<Dimension>& u)
+{
+ return u - 2 * dot(u - origin, normal) * normal;
+}
+
+#endif // SYMMETRY_UTILS_HPP
diff --git a/src/geometry/TriangleTransformation.hpp b/src/geometry/TriangleTransformation.hpp
index df8444d894a7ed56c0fe3825ffbf1cc14a742531..c49db368684d8c96aeb5061041828799e8509fac 100644
--- a/src/geometry/TriangleTransformation.hpp
+++ b/src/geometry/TriangleTransformation.hpp
@@ -66,6 +66,7 @@ class TriangleTransformation
private:
TinyMatrix<Dimension, 2> m_jacobian;
TinyVector<Dimension> m_shift;
+ TinyVector<Dimension> m_area_variation;
double m_area_variation_norm;
public:
@@ -76,6 +77,21 @@ class TriangleTransformation
return m_jacobian * x + m_shift;
}
+ PUGS_INLINE
+ TinyVector<Dimension>
+ areaVariation() const
+ {
+ return m_area_variation;
+ }
+
+ PUGS_INLINE
+ TinyVector<Dimension>
+ areaVariation(const TinyVector<2>&) const
+ {
+ return m_area_variation;
+ }
+
+ PUGS_INLINE
double
areaVariationNorm() const
{
@@ -92,7 +108,8 @@ class TriangleTransformation
m_shift = a;
- m_area_variation_norm = l2Norm(crossProduct(b - a, c - a));
+ m_area_variation = crossProduct(b - a, c - a);
+ m_area_variation_norm = l2Norm(m_area_variation);
}
~TriangleTransformation() = default;
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 8385e55ccbc7ebcfedfd6880d7def1df5a381847..c03b412af73bb25d0b71eeb060e33484c02eb824 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -1,5 +1,7 @@
# ------------------- Source files --------------------
+add_subdirectory(reconstruction_utils)
+
add_library(
PugsScheme
AcousticSolver.cpp
@@ -30,3 +32,9 @@ target_link_libraries(
PugsScheme
${HIGHFIVE_TARGET}
)
+
+# Additional dependencies
+add_dependencies(
+ PugsScheme
+ PugsSchemeReconstructionUtils
+)
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index cf4784a285a85435ddd650d44ceac9d80ff06032..f3d3768d9dd4a880386c66c20ca69fd5847be8e2 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -2,19 +2,8 @@
#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 <geometry/SymmetryUtils.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshFlatFaceBoundary.hpp>
@@ -24,702 +13,532 @@
#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
+#include <scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp>
+#include <scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp>
+#include <scheme/reconstruction_utils/MutableDiscreteFunctionDPkVariant.hpp>
-template <size_t Dimension>
-PUGS_INLINE auto
-symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimension>& u)
-{
- return u - 2 * dot(u, normal) * normal;
-}
-
-template <size_t Dimension>
-PUGS_INLINE auto
-symmetrize_matrix(const TinyVector<Dimension>& normal, const TinyMatrix<Dimension>& A)
-{
- const TinyMatrix S = TinyMatrix<Dimension>{identity} - 2 * tensorProduct(normal, normal);
- return S * A * S;
-}
-
-template <size_t Dimension>
-PUGS_INLINE auto
-symmetrize_coordinates(const TinyVector<Dimension>& origin,
- const TinyVector<Dimension>& normal,
- const TinyVector<Dimension>& u)
-{
- return u - 2 * dot(u - origin, normal) * normal;
-}
-
-class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
+template <MeshConcept MeshType>
+class PolynomialReconstruction::Internal
{
- public:
- using Variant = std::variant<DiscreteFunctionDPk<1, double>,
- DiscreteFunctionDPk<1, TinyVector<1>>,
- DiscreteFunctionDPk<1, TinyVector<2>>,
- DiscreteFunctionDPk<1, TinyVector<3>>,
- DiscreteFunctionDPk<1, TinyMatrix<1>>,
- DiscreteFunctionDPk<1, TinyMatrix<2>>,
- DiscreteFunctionDPk<1, TinyMatrix<3>>,
-
- DiscreteFunctionDPk<2, double>,
- DiscreteFunctionDPk<2, TinyVector<1>>,
- DiscreteFunctionDPk<2, TinyVector<2>>,
- DiscreteFunctionDPk<2, TinyVector<3>>,
- DiscreteFunctionDPk<2, TinyMatrix<1>>,
- DiscreteFunctionDPk<2, TinyMatrix<2>>,
- DiscreteFunctionDPk<2, TinyMatrix<3>>,
-
- DiscreteFunctionDPk<3, double>,
- DiscreteFunctionDPk<3, TinyVector<1>>,
- DiscreteFunctionDPk<3, TinyVector<2>>,
- DiscreteFunctionDPk<3, TinyVector<3>>,
- DiscreteFunctionDPk<3, TinyMatrix<1>>,
- DiscreteFunctionDPk<3, TinyMatrix<2>>,
- DiscreteFunctionDPk<3, TinyMatrix<3>>,
-
- DiscreteFunctionDPkVector<1, double>,
- DiscreteFunctionDPkVector<1, TinyVector<1>>,
- DiscreteFunctionDPkVector<1, TinyVector<2>>,
- DiscreteFunctionDPkVector<1, TinyVector<3>>,
- DiscreteFunctionDPkVector<1, TinyMatrix<1>>,
- DiscreteFunctionDPkVector<1, TinyMatrix<2>>,
- DiscreteFunctionDPkVector<1, TinyMatrix<3>>,
-
- DiscreteFunctionDPkVector<2, double>,
- DiscreteFunctionDPkVector<2, TinyVector<1>>,
- DiscreteFunctionDPkVector<2, TinyVector<2>>,
- DiscreteFunctionDPkVector<2, TinyVector<3>>,
- DiscreteFunctionDPkVector<2, TinyMatrix<1>>,
- DiscreteFunctionDPkVector<2, TinyMatrix<2>>,
- DiscreteFunctionDPkVector<2, TinyMatrix<3>>,
-
- DiscreteFunctionDPkVector<3, double>,
- DiscreteFunctionDPkVector<3, TinyVector<1>>,
- DiscreteFunctionDPkVector<3, TinyVector<2>>,
- DiscreteFunctionDPkVector<3, TinyVector<3>>,
- DiscreteFunctionDPkVector<3, TinyMatrix<1>>,
- DiscreteFunctionDPkVector<3, TinyMatrix<2>>,
- DiscreteFunctionDPkVector<3, TinyMatrix<3>>>;
-
private:
- Variant m_mutable_discrete_function_dpk;
+ using Rd = TinyVector<MeshType::Dimension>;
- public:
- PUGS_INLINE
- const Variant&
- mutableDiscreteFunctionDPk() const
- {
- return m_mutable_discrete_function_dpk;
- }
+ friend PolynomialReconstruction;
- template <typename DiscreteFunctionDPkT>
- PUGS_INLINE auto&&
- get() const
+ template <typename MatrixType>
+ static void
+ buildB(const CellId& cell_j_id,
+ const CellToCellStencilArray& stencil_array,
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
+ SmallArray<const Rd> symmetry_normal_list,
+ ShrinkMatrixView<MatrixType>& B)
{
- static_assert(is_discrete_function_dpk_v<DiscreteFunctionDPkT>, "invalid template argument");
-
-#ifndef NDEBUG
- if (not std::holds_alternative<DiscreteFunctionDPkT>(this->m_mutable_discrete_function_dpk)) {
- std::ostringstream error_msg;
- error_msg << "invalid discrete function type\n";
- error_msg << "- required " << rang::fgB::red << demangle<DiscreteFunctionDPkT>() << rang::fg::reset << '\n';
- error_msg << "- contains " << rang::fgB::yellow
- << std::visit([](auto&& f) -> std::string { return demangle<decltype(f)>(); },
- this->m_mutable_discrete_function_dpk)
- << rang::fg::reset;
- throw NormalError(error_msg.str());
- }
-#endif // NDEBUG
-
- return std::get<DiscreteFunctionDPkT>(this->mutableDiscreteFunctionDPk());
- }
+ auto stencil_cell_list = stencil_array[cell_j_id];
+
+ size_t column_begin = 0;
+ for (size_t i_discrete_function_variant = 0; i_discrete_function_variant < discrete_function_variant_list.size();
+ ++i_discrete_function_variant) {
+ const auto& discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
+
+ std::visit(
+ [&](auto&& discrete_function) {
+ using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+ if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+ using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+ const DataType& qj = discrete_function[cell_j_id];
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const DataType& qi_qj = discrete_function[cell_i_id] - qj;
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ B(index, column_begin) = qi_qj;
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
+ B(index, kB) = qi_qj[k];
+ }
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+ const size_t kB = column_begin + p * DataType::NumberOfColumns;
+ for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+ B(index, kB + q) = qi_qj(p, q);
+ }
+ }
+ }
+ }
- template <size_t Dimension, typename DataType>
- MutableDiscreteFunctionDPkVariant(const DiscreteFunctionDPk<Dimension, DataType>& discrete_function_dpk)
- : m_mutable_discrete_function_dpk{discrete_function_dpk}
- {
- static_assert(std::is_same_v<DataType, double> or //
- std::is_same_v<DataType, TinyVector<1, double>> or //
- std::is_same_v<DataType, TinyVector<2, double>> or //
- std::is_same_v<DataType, TinyVector<3, double>> or //
- std::is_same_v<DataType, TinyMatrix<1, 1, double>> or //
- std::is_same_v<DataType, TinyMatrix<2, 2, double>> or //
- std::is_same_v<DataType, TinyMatrix<3, 3, double>>,
- "DiscreteFunctionDPk with this DataType is not allowed in variant");
- }
+ for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+ ++i_symmetry) {
+ auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ const DataType& qi_qj = discrete_function[cell_i_id] - qj;
+ B(index, column_begin) = qi_qj;
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ if constexpr (DataType::Dimension == MeshType::Dimension) {
+ const Rd& normal = symmetry_normal_list[i_symmetry];
- template <size_t Dimension, typename DataType>
- MutableDiscreteFunctionDPkVariant(const DiscreteFunctionDPkVector<Dimension, DataType>& discrete_function_dpk)
- : m_mutable_discrete_function_dpk{discrete_function_dpk}
- {
- static_assert(std::is_same_v<DataType, double> or //
- std::is_same_v<DataType, TinyVector<1, double>> or //
- std::is_same_v<DataType, TinyVector<2, double>> or //
- std::is_same_v<DataType, TinyVector<3, double>> or //
- std::is_same_v<DataType, TinyMatrix<1, 1, double>> or //
- std::is_same_v<DataType, TinyMatrix<2, 2, double>> or //
- std::is_same_v<DataType, TinyMatrix<3, 3, double>>,
- "DiscreteFunctionDPkVector with this DataType is not allowed in variant");
- }
+ const DataType& qi = discrete_function[cell_i_id];
+ const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
+ for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
+ B(index, kB) = qi_qj[k];
+ }
+ } else {
+ // LCOV_EXCL_START
+ std::stringstream error_msg;
+ error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+ << " using a mesh of dimension " << MeshType::Dimension;
+ throw UnexpectedError(error_msg.str());
+ // LCOV_EXCL_STOP
+ }
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
+ (DataType::NumberOfColumns == MeshType::Dimension)) {
+ const Rd& normal = symmetry_normal_list[i_symmetry];
+
+ const DataType& qi = discrete_function[cell_i_id];
+ const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
+ for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+ for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+ B(index, column_begin + p * DataType::NumberOfColumns + q) = qi_qj(p, q);
+ }
+ }
+ } else {
+ // LCOV_EXCL_START
+ std::stringstream error_msg;
+ error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
+ << DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
+ throw UnexpectedError(error_msg.str());
+ // LCOV_EXCL_STOP
+ }
+ }
+ }
+ }
- MutableDiscreteFunctionDPkVariant& operator=(MutableDiscreteFunctionDPkVariant&&) = default;
- MutableDiscreteFunctionDPkVariant& operator=(const MutableDiscreteFunctionDPkVariant&) = default;
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ ++column_begin;
+ } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
+ column_begin += DataType::Dimension;
+ }
+ } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+ using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
- MutableDiscreteFunctionDPkVariant(const MutableDiscreteFunctionDPkVariant&) = default;
- MutableDiscreteFunctionDPkVariant(MutableDiscreteFunctionDPkVariant&&) = default;
+ const auto qj_vector = discrete_function[cell_j_id];
- MutableDiscreteFunctionDPkVariant() = delete;
- ~MutableDiscreteFunctionDPkVariant() = default;
-};
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ for (size_t l = 0; l < qj_vector.size(); ++l) {
+ const DataType& qj = qj_vector[l];
+ const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
+ B(index, column_begin + l) = qi_qj;
+ }
+ }
-template <typename ConformTransformationT>
-void
-_computeEjkMean(const QuadratureFormula<1>& quadrature,
- const ConformTransformationT& T,
- const TinyVector<1>& Xj,
- const double Vi,
- const size_t degree,
- const size_t dimension,
- SmallArray<double>& inv_Vj_wq_detJ_ek,
- SmallArray<double>& mean_of_ejk)
-{
- using Rd = TinyVector<1>;
- mean_of_ejk.fill(0);
+ for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+ ++i_symmetry) {
+ auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+ for (size_t l = 0; l < qj_vector.size(); ++l) {
+ const DataType& qj = qj_vector[l];
+ const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
+ B(index, column_begin + l) = qi_qj;
+ }
+ }
+ }
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ for (size_t l = 0; l < qj_vector.size(); ++l) {
+ const DataType& qj = qj_vector[l];
+ const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
+ for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension; ++k, ++kB) {
+ B(index, kB) = qi_qj[k];
+ }
+ }
+ }
- 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 = T.jacobianDeterminant();
+ for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+ ++i_symmetry) {
+ if constexpr (DataType::Dimension == MeshType::Dimension) {
+ auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
- const Rd X_Xj = T(xi_q) - Xj;
+ const Rd& normal = symmetry_normal_list[i_symmetry];
- const double x_xj = X_Xj[0];
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
- {
- 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 l = 0; l < qj_vector.size(); ++l) {
+ const DataType& qj = qj_vector[l];
+ const DataType& qi = discrete_function[cell_i_id][l];
+ const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
+ for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
+ ++k, ++kB) {
+ B(index, kB) = qi_qj[k];
+ }
+ }
+ }
+ } else {
+ // LCOV_EXCL_START
+ std::stringstream error_msg;
+ error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+ << " using a mesh of dimension " << MeshType::Dimension;
+ throw UnexpectedError(error_msg.str());
+ // LCOV_EXCL_STOP
+ }
+ }
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ for (size_t l = 0; l < qj_vector.size(); ++l) {
+ const DataType& qj = qj_vector[l];
+ const DataType& qi = discrete_function[cell_i_id][l];
+ const DataType& qi_qj = qi - qj;
+
+ for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+ const size_t kB = column_begin + l * DataType::Dimension + p * DataType::NumberOfColumns;
+ for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+ B(index, kB + q) = qi_qj(p, q);
+ }
+ }
+ }
+ }
- for (size_t k = 1; k < dimension; ++k) {
- mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
- }
- }
-}
+ for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+ ++i_symmetry) {
+ if constexpr ((DataType::NumberOfRows == MeshType::Dimension) and
+ (DataType::NumberOfColumns == MeshType::Dimension)) {
+ auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
-template <typename ConformTransformationT>
-void
-_computeEjkMean(const QuadratureFormula<2>& quadrature,
- const ConformTransformationT& T,
- const TinyVector<2>& Xj,
- const double Vi,
- const size_t degree,
- const size_t dimension,
- SmallArray<double>& inv_Vj_wq_detJ_ek,
- SmallArray<double>& mean_of_ejk)
-{
- using Rd = TinyVector<2>;
-
- 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<TriangleTransformation<2>, std::decay_t<decltype(T)>>) {
- return T.jacobianDeterminant();
- } else {
- return T.jacobianDeterminant(xi_q);
- }
- }();
+ const Rd& normal = symmetry_normal_list[i_symmetry];
- const Rd X_Xj = T(xi_q) - Xj;
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
- const double x_xj = X_Xj[0];
- const double y_yj = X_Xj[1];
+ for (size_t l = 0; l < qj_vector.size(); ++l) {
+ const DataType& qj = qj_vector[l];
+ const DataType& qi = discrete_function[cell_i_id][l];
+ const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
- {
- 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 p = 0; p < DataType::NumberOfRows; ++p) {
+ const size_t kB = column_begin + l * DataType::Dimension + p * DataType::NumberOfColumns;
+ for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+ B(index, kB + q) = qi_qj(p, q);
+ }
+ }
+ }
+ }
+ } else {
+ // LCOV_EXCL_START
+ std::stringstream error_msg;
+ error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+ << " using a mesh of dimension " << MeshType::Dimension;
+ throw UnexpectedError(error_msg.str());
+ // LCOV_EXCL_STOP
+ }
+ }
+ }
- 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];
- }
- }
- }
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ column_begin += qj_vector.size();
+ } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
+ column_begin += qj_vector.size() * DataType::Dimension;
+ }
- for (size_t k = 1; k < dimension; ++k) {
- mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("invalid discrete function type");
+ // LCOV_EXCL_STOP
+ }
+ },
+ discrete_function_variant->discreteFunction());
}
}
-}
-template <typename ConformTransformationT>
-void
-_computeEjkMean(const QuadratureFormula<3>& quadrature,
- const ConformTransformationT& T,
- const TinyVector<3>& Xj,
- const double Vi,
- const size_t degree,
- const size_t dimension,
- SmallArray<double>& inv_Vj_wq_detJ_ek,
- SmallArray<double>& mean_of_ejk)
-{
- using Rd = TinyVector<3>;
- 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);
+ template <typename MatrixType>
+ static void
+ rowWeighting(const CellId& cell_j_id,
+ const CellToCellStencilArray& stencil_array,
+ const CellValue<const Rd>& xj,
+ const SmallArray<const Rd>& symmetry_origin_list,
+ const SmallArray<const Rd>& symmetry_normal_list,
+ ShrinkMatrixView<MatrixType>& A,
+ ShrinkMatrixView<MatrixType>& B)
+ {
+ // Add row weighting (give more importance to cells that are
+ // closer to j)
+ auto stencil_cell_list = stencil_array[cell_j_id];
+
+ const Rd& Xj = xj[cell_j_id];
+
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const double weight = 1. / l2Norm(xj[cell_i_id] - Xj);
+ for (size_t l = 0; l < A.numberOfColumns(); ++l) {
+ A(index, l) *= weight;
}
- }();
-
- 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 l = 0; l < B.numberOfColumns(); ++l) {
+ B(index, l) *= weight;
}
-
- 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_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry) {
+ auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+ const Rd& origin = symmetry_origin_list[i_symmetry];
+ const Rd& normal = symmetry_normal_list[i_symmetry];
+
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+ const double weight = 1. / l2Norm(symmetrize_coordinates(origin, normal, xj[cell_i_id]) - Xj);
+ for (size_t l = 0; l < A.numberOfColumns(); ++l) {
+ A(index, l) *= weight;
}
- }
-
- for (size_t i_z = 1; i_z <= degree; ++i_z) {
- const size_t begin_i_z_1 = ((i_z - 1) * (3 * (degree + 2) * (degree + 3) + (i_z - (3 * degree + 8)) * i_z)) / 6;
- 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 l = 0; l < B.numberOfColumns(); ++l) {
+ B(index, l) *= weight;
}
}
}
-
- for (size_t k = 1; k < dimension; ++k) {
- mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
- }
}
-}
-
-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
-_computeEjkMeanInSymmetricCell(const TinyVector<1>& origin,
- const TinyVector<1>& normal,
- 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 auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
- const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
- const LineTransformation<1> T{x0, x1};
-
- 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 ConformTransformationT>
-void _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
- const ConformTransformationT& T,
- const TinyVector<2>& Xj,
- const double Vi,
- const size_t degree,
- const size_t dimension,
- SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
- SmallArray<double>& mean_of_ejk);
-
-template <>
-void
-_computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
- const LineTransformation<2>& T,
- const TinyVector<2>& Xj,
- const double Vi,
- const size_t degree,
- const size_t dimension,
- SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
- SmallArray<double>& mean_of_ejk)
-{
- using Rd = TinyVector<2>;
-
- const double velocity_perp_e1 = T.velocity()[1];
-
- for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
- const double wq = quadrature.weight(i_q);
- const TinyVector<1> xi_q = quadrature.point(i_q);
-
- const Rd X_Xj = T(xi_q) - Xj;
-
- const double x_xj = X_Xj[0];
- const double y_yj = X_Xj[1];
+ template <typename MatrixType>
+ static void
+ solveCollectionInPlaceWithPreconditionner(const ShrinkMatrixView<MatrixType>& A,
+ const SmallMatrix<double>& X,
+ const ShrinkMatrixView<MatrixType>& B,
+ const SmallVector<double>& G)
+ {
+ for (size_t l = 0; l < A.numberOfColumns(); ++l) {
+ double g = 0;
+ for (size_t i = 0; i < A.numberOfRows(); ++i) {
+ const double Ail = A(i, l);
- {
- size_t k = 0;
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1 / Vi;
- for (; k <= degree; ++k) {
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
- x_xj * inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * ((1. * k) / (k + 1));
+ g += Ail * Ail;
}
+ G[l] = std::sqrt(g);
+ }
- 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_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
- }
+ for (size_t l = 0; l < A.numberOfColumns(); ++l) {
+ const double Gl = G[l];
+ for (size_t i = 0; i < A.numberOfRows(); ++i) {
+ A(i, l) *= Gl;
}
}
- for (size_t k = 1; k < dimension; ++k) {
- mean_of_ejk[k - 1] += inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
- }
- }
-}
+ Givens::solveCollectionInPlace(A, X, B);
-template <MeshConcept MeshType>
-void
-_computeEjkMeanByBoundary(const MeshType& mesh,
- const TinyVector<2>& Xj,
- const CellId& cell_id,
- const auto& cell_to_face_matrix,
- const auto& face_to_node_matrix,
- const auto& cell_face_is_reversed,
- const CellValue<const double>& Vi,
- const size_t degree,
- const size_t basis_dimension,
- SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
- SmallArray<double>& mean_of_ejk)
-{
- const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(degree + 1));
-
- auto xr = mesh.xr();
- mean_of_ejk.fill(0);
- auto cell_face_list = cell_to_face_matrix[cell_id];
- for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- bool is_reversed = cell_face_is_reversed[cell_id][i_face];
-
- const FaceId face_id = cell_face_list[i_face];
- auto face_node_list = face_to_node_matrix[face_id];
- if (is_reversed) {
- const LineTransformation<2> T{xr[face_node_list[1]], xr[face_node_list[0]]};
- _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
- } else {
- const LineTransformation<2> T{xr[face_node_list[0]], xr[face_node_list[1]]};
- _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
- }
- }
-}
-
-void
-_computeEjkMeanByBoundaryInSymmetricCell(const TinyVector<2>& origin,
- const TinyVector<2>& normal,
- const TinyVector<2>& Xj,
- const CellId& cell_id,
- const NodeValue<const TinyVector<2>>& xr,
- const auto& cell_to_face_matrix,
- const auto& face_to_node_matrix,
- const auto& cell_face_is_reversed,
- const CellValue<const double>& Vi,
- const size_t degree,
- const size_t basis_dimension,
- SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
- SmallArray<double>& mean_of_ejk)
-{
- const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(degree + 1));
-
- mean_of_ejk.fill(0);
- auto cell_face_list = cell_to_face_matrix[cell_id];
- for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- bool is_reversed = cell_face_is_reversed[cell_id][i_face];
-
- const FaceId face_id = cell_face_list[i_face];
- auto face_node_list = face_to_node_matrix[face_id];
-
- const auto x0 = symmetrize_coordinates(origin, normal, xr[face_node_list[1]]);
- const auto x1 = symmetrize_coordinates(origin, normal, xr[face_node_list[0]]);
-
- if (is_reversed) {
- const LineTransformation<2> T{x1, x0};
- _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
- } else {
- const LineTransformation<2> T{x0, x1};
- _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ for (size_t l = 0; l < X.numberOfRows(); ++l) {
+ const double Gl = G[l];
+ for (size_t i = 0; i < X.numberOfColumns(); ++i) {
+ X(l, i) *= Gl;
+ }
}
}
-}
-
-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
-_computeEjkMeanInSymmetricCell(const TinyVector<2>& origin,
- const TinyVector<2>& normal,
- 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 auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
- const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
- const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
-
- const TriangleTransformation<2> T{x0, x1, x2};
- 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 auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
- const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
- const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
- const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
-
- const SquareTransformation<2> T{x0, x1, x2, x3};
- 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)});
- }
-}
-
-template <typename NodeListT>
-void
-_computeEjkMeanInSymmetricCell(const TinyVector<3>& origin,
- const TinyVector<3>& normal,
- 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 auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
- const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
- const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
- const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
-
- const TetrahedronTransformation T{x0, x1, x2, x3};
-
- 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 auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
- const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
- const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
- const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[4]]);
- const auto x4 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
- const auto x5 = symmetrize_coordinates(origin, normal, xr[node_list[5]]);
-
- const PrismTransformation T{x0, x1, x2, //
- x3, x4, x5};
+ template <typename ReconstructionMatrixBuilderType>
+ static void
+ populateDiscreteFunctionDPkByCell(
+ const CellId& cell_j_id,
+ const size_t& degree,
+ const SmallMatrix<double>& X,
+ const ReconstructionMatrixBuilderType& reconstruction_matrix_builder,
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
+ const std::vector<PolynomialReconstruction::MutableDiscreteFunctionDPkVariant>&
+ mutable_discrete_function_dpk_variant_list)
+ {
+ size_t column_begin = 0;
+ for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size(); ++i_dpk_variant) {
+ const auto& dpk_variant = mutable_discrete_function_dpk_variant_list[i_dpk_variant];
- const auto& quadrature = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{degree + 1});
+ const auto& discrete_function_variant = discrete_function_variant_list[i_dpk_variant];
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ std::visit(
+ [&](auto&& dpk_function, auto&& p0_function) {
+ using DPkFunctionT = std::decay_t<decltype(dpk_function)>;
+ using P0FunctionT = std::decay_t<decltype(p0_function)>;
+ using DataType = std::remove_const_t<std::decay_t<typename DPkFunctionT::data_type>>;
+ using P0DataType = std::remove_const_t<std::decay_t<typename P0FunctionT::data_type>>;
- } else if (cell_type == CellType::Pyramid) {
- const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
- const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
- const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
- const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
- const auto x4 = symmetrize_coordinates(origin, normal, xr[node_list[4]]);
- const PyramidTransformation T{x0, x1, x2, x3, x4};
+ if constexpr (std::is_same_v<DataType, P0DataType>) {
+ if constexpr (is_discrete_function_P0_v<P0FunctionT>) {
+ if constexpr (is_discrete_function_dpk_scalar_v<DPkFunctionT>) {
+ auto dpk_j = dpk_function.coefficients(cell_j_id);
+ dpk_j[0] = p0_function[cell_j_id];
- const auto& quadrature = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{degree + 1});
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ if (degree > 1) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ dpk_j[0] -= X(i, column_begin) * mean_j_of_ejk[i];
+ }
+ }
+ }
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_ip1 = dpk_j[i + 1];
+ dpk_j_ip1 = X(i, column_begin);
+ }
+ ++column_begin;
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ if (degree > 1) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_0 = dpk_j[0];
+ for (size_t k = 0; k < DataType::Dimension; ++k) {
+ dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
+ }
+ }
+ }
+ }
- } else if (cell_type == CellType::Hexahedron) {
- const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
- const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
- const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
- const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
- const auto x4 = symmetrize_coordinates(origin, normal, xr[node_list[7]]);
- const auto x5 = symmetrize_coordinates(origin, normal, xr[node_list[6]]);
- const auto x6 = symmetrize_coordinates(origin, normal, xr[node_list[5]]);
- const auto x7 = symmetrize_coordinates(origin, normal, xr[node_list[4]]);
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_ip1 = dpk_j[i + 1];
+ for (size_t k = 0; k < DataType::Dimension; ++k) {
+ dpk_j_ip1[k] = X(i, column_begin + k);
+ }
+ }
+ column_begin += DataType::Dimension;
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ if (degree > 1) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_0 = dpk_j[0];
+ for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+ for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+ dpk_j_0(k, l) -= X(i, column_begin + k * DataType::NumberOfColumns + l) * mean_j_of_ejk[i];
+ }
+ }
+ }
+ }
+ }
- const CubeTransformation T{x0, x1, x2, x3, x4, x5, x6, x7};
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_ip1 = dpk_j[i + 1];
+ for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+ for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+ dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
+ }
+ }
+ }
+ column_begin += DataType::Dimension;
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected data type");
+ // LCOV_EXCL_STOP
+ }
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected discrete dpk function type");
+ // LCOV_EXCL_STOP
+ }
+ } else if constexpr (is_discrete_function_P0_vector_v<P0FunctionT>) {
+ if constexpr (is_discrete_function_dpk_vector_v<DPkFunctionT>) {
+ auto dpk_j = dpk_function.coefficients(cell_j_id);
+ auto cell_vector = p0_function[cell_j_id];
+ const size_t size = X.numberOfRows() + 1;
+
+ for (size_t l = 0; l < cell_vector.size(); ++l) {
+ const size_t component_begin = l * size;
+ dpk_j[component_begin] = cell_vector[l];
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ if (degree > 1) {
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ dpk_j[component_begin] -= X(i, column_begin) * mean_j_of_ejk[i];
+ }
+ }
+ }
- const auto& quadrature =
- QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{degree + 1});
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
+ dpk_j_ip1 = X(i, column_begin);
+ }
+ ++column_begin;
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ if (degree > 1) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ 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];
+ }
+ }
+ }
+ }
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
+ for (size_t k = 0; k < DataType::Dimension; ++k) {
+ dpk_j_ip1[k] = X(i, column_begin + k);
+ }
+ }
+ column_begin += DataType::Dimension;
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ if constexpr (ReconstructionMatrixBuilderType::handles_high_degrees) {
+ if (degree > 1) {
+ const auto& mean_j_of_ejk = reconstruction_matrix_builder.meanjOfEjk();
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_0 = dpk_j[component_begin];
+ for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+ for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+ dpk_j_0(p, q) -=
+ X(i, column_begin + p * DataType::NumberOfColumns + q) * mean_j_of_ejk[i];
+ }
+ }
+ }
+ }
+ }
- } else {
- throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ for (size_t i = 0; i < X.numberOfRows(); ++i) {
+ auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
+ for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+ for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+ dpk_j_ip1(p, q) = X(i, column_begin + p * DataType::NumberOfColumns + q);
+ }
+ }
+ }
+ column_begin += DataType::Dimension;
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected data type");
+ // LCOV_EXCL_STOP
+ }
+ }
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected discrete dpk function type");
+ // LCOV_EXCL_STOP
+ }
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("unexpected discrete function type");
+ // LCOV_EXCL_STOP
+ }
+ } else {
+ // LCOV_EXCL_START
+ throw UnexpectedError("incompatible data types");
+ // LCOV_EXCL_STOP
+ }
+ },
+ dpk_variant.mutableDiscreteFunctionDPk(), discrete_function_variant->discreteFunction());
+ }
}
-}
+};
size_t
PolynomialReconstruction::_getNumberOfColumns(
@@ -836,12 +655,14 @@ PolynomialReconstruction::_checkDataAndSymmetriesCompatibility(
}
}
-template <MeshConcept MeshType>
+template <typename ReconstructionMatrixBuilderType, MeshConcept MeshType>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::_build(
const std::shared_ptr<const MeshType>& p_mesh,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
+ static_assert(std::is_same_v<MeshType, typename ReconstructionMatrixBuilderType::MeshType>);
+
const MeshType& mesh = *p_mesh;
using Rd = TinyVector<MeshType::Dimension>;
@@ -863,10 +684,8 @@ PolynomialReconstruction::_build(
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
- auto cell_is_owned = mesh.connectivity().cellIsOwned();
- auto cell_type = mesh.connectivity().cellType();
- auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- auto cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+ auto cell_is_owned = mesh.connectivity().cellIsOwned();
+ auto cell_type = mesh.connectivity().cellType();
auto full_stencil_size = [&](const CellId cell_id) {
auto stencil_cell_list = stencil_array[cell_id];
@@ -926,24 +745,19 @@ PolynomialReconstruction::_build(
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());
- SmallArray<SmallArray<double>> mean_j_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
- SmallArray<SmallArray<double>> mean_i_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
-
- SmallArray<SmallArray<double>> inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool(Kokkos::DefaultExecutionSpace::concurrency());
- SmallArray<SmallArray<double>> mean_l_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
-
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] = 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);
- mean_j_of_ejk_pool[i] = SmallArray<double>(basis_dimension - 1);
- mean_i_of_ejk_pool[i] = SmallArray<double>(basis_dimension - 1);
+ SmallArray<std::shared_ptr<ReconstructionMatrixBuilderType>> reconstruction_matrix_builder_pool(A_pool.size());
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[i] = SmallArray<double>(basis_dimension);
+ for (size_t t = 0; t < reconstruction_matrix_builder_pool.size(); ++t) {
+ reconstruction_matrix_builder_pool[t] =
+ std::make_shared<ReconstructionMatrixBuilderType>(*p_mesh, m_descriptor.degree(), symmetry_origin_list,
+ symmetry_normal_list, stencil_array);
}
parallel_for(
@@ -951,544 +765,34 @@ PolynomialReconstruction::_build(
if (cell_is_owned[cell_j_id]) {
const int32_t t = tokens.acquire();
- auto stencil_cell_list = stencil_array[cell_j_id];
-
- ShrinkMatrixView B(B_pool[t], full_stencil_size(cell_j_id));
-
- size_t column_begin = 0;
- for (size_t i_discrete_function_variant = 0;
- i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
- const auto& discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
-
- std::visit(
- [&](auto&& discrete_function) {
- using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
- if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
- using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
- const DataType& qj = discrete_function[cell_j_id];
- size_t index = 0;
- for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = stencil_cell_list[i];
- const DataType& qi_qj = discrete_function[cell_i_id] - qj;
- if constexpr (std::is_arithmetic_v<DataType>) {
- B(index, column_begin) = qi_qj;
- } else if constexpr (is_tiny_vector_v<DataType>) {
- for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
- B(index, kB) = qi_qj[k];
- }
- } else if constexpr (is_tiny_matrix_v<DataType>) {
- for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
- }
- }
- }
- }
-
- for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
- ++i_symmetry) {
- auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
- auto ghost_cell_list = ghost_stencil[cell_j_id];
- for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = ghost_cell_list[i];
- if constexpr (std::is_arithmetic_v<DataType>) {
- const DataType& qi_qj = discrete_function[cell_i_id] - qj;
- B(index, column_begin) = qi_qj;
- } else if constexpr (is_tiny_vector_v<DataType>) {
- if constexpr (DataType::Dimension == MeshType::Dimension) {
- const Rd& normal = symmetry_normal_list[i_symmetry];
-
- const DataType& qi = discrete_function[cell_i_id];
- const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
- for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
- B(index, kB) = qi_qj[k];
- }
- } else {
- // LCOV_EXCL_START
- std::stringstream error_msg;
- error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
- << " using a mesh of dimension " << MeshType::Dimension;
- throw UnexpectedError(error_msg.str());
- // LCOV_EXCL_STOP
- }
- } else if constexpr (is_tiny_matrix_v<DataType>) {
- if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
- (DataType::NumberOfColumns == MeshType::Dimension)) {
- const Rd& normal = symmetry_normal_list[i_symmetry];
- const DataType& qi = discrete_function[cell_i_id];
- const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
- for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
- }
- }
- } else {
- // LCOV_EXCL_START
- std::stringstream error_msg;
- error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
- << DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
- throw UnexpectedError(error_msg.str());
- // LCOV_EXCL_STOP
- }
- }
- }
- }
-
- if constexpr (std::is_arithmetic_v<DataType>) {
- ++column_begin;
- } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
- column_begin += DataType::Dimension;
- }
- } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
- using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
-
- const auto qj_vector = discrete_function[cell_j_id];
-
- if constexpr (std::is_arithmetic_v<DataType>) {
- size_t index = 0;
- for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = stencil_cell_list[i];
- for (size_t l = 0; l < qj_vector.size(); ++l) {
- const DataType& qj = qj_vector[l];
- const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
- B(index, column_begin + l) = qi_qj;
- }
- }
-
- for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
- ++i_symmetry) {
- auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
- auto ghost_cell_list = ghost_stencil[cell_j_id];
- for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = ghost_cell_list[i];
- for (size_t l = 0; l < qj_vector.size(); ++l) {
- const DataType& qj = qj_vector[l];
- const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
- B(index, column_begin + l) = qi_qj;
- }
- }
- }
- } else if constexpr (is_tiny_vector_v<DataType>) {
- size_t index = 0;
- for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = stencil_cell_list[i];
- for (size_t l = 0; l < qj_vector.size(); ++l) {
- const DataType& qj = qj_vector[l];
- const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
- for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
- ++k, ++kB) {
- B(index, kB) = qi_qj[k];
- }
- }
- }
-
- for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
- ++i_symmetry) {
- if constexpr (DataType::Dimension == MeshType::Dimension) {
- auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
- auto ghost_cell_list = ghost_stencil[cell_j_id];
-
- const Rd& normal = symmetry_normal_list[i_symmetry];
-
- for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = ghost_cell_list[i];
-
- for (size_t l = 0; l < qj_vector.size(); ++l) {
- const DataType& qj = qj_vector[l];
- const DataType& qi = discrete_function[cell_i_id][l];
- const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
- for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
- ++k, ++kB) {
- B(index, kB) = qi_qj[k];
- }
- }
- }
- } else {
- // LCOV_EXCL_START
- std::stringstream error_msg;
- error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
- << " using a mesh of dimension " << MeshType::Dimension;
- throw UnexpectedError(error_msg.str());
- // LCOV_EXCL_STOP
- }
- }
- } else if constexpr (is_tiny_matrix_v<DataType>) {
- for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
- ++i_symmetry) {
- auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
- auto ghost_cell_list = ghost_stencil[cell_j_id];
- if (ghost_cell_list.size() > 0) {
- throw NotImplementedError("TinyMatrix symmetries for reconstruction of DiscreteFunctionP0Vector");
- }
- }
- }
-
- if constexpr (std::is_arithmetic_v<DataType>) {
- column_begin += qj_vector.size();
- } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
- column_begin += qj_vector.size() * DataType::Dimension;
- }
-
- } else {
- // LCOV_EXCL_START
- throw UnexpectedError("invalid discrete function type");
- // LCOV_EXCL_STOP
- }
- },
- discrete_function_variant->discreteFunction());
- }
-
ShrinkMatrixView A(A_pool[t], full_stencil_size(cell_j_id));
+ ShrinkMatrixView B(B_pool[t], full_stencil_size(cell_j_id));
- if ((m_descriptor.degree() == 1) and
- (m_descriptor.integrationMethodType() == IntegrationMethodType::cell_center)) {
- const Rd& Xj = xj[cell_j_id];
- size_t index = 0;
- for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = stencil_cell_list[i];
- const Rd Xi_Xj = xj[cell_i_id] - Xj;
- for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(index, l) = Xi_Xj[l];
- }
- }
- for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
- ++i_symmetry) {
- auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
- auto ghost_cell_list = ghost_stencil[cell_j_id];
-
- const Rd& origin = symmetry_origin_list[i_symmetry];
- const Rd& normal = symmetry_normal_list[i_symmetry];
-
- for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = ghost_cell_list[i];
- const Rd Xi_Xj = symmetrize_coordinates(origin, normal, xj[cell_i_id]) - Xj;
- for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(index, l) = Xi_Xj[l];
- }
- }
- }
-
- } else if ((m_descriptor.integrationMethodType() == IntegrationMethodType::element) or
- (m_descriptor.integrationMethodType() == IntegrationMethodType::boundary)) {
- if ((m_descriptor.integrationMethodType() == IntegrationMethodType::boundary) and
- (MeshType::Dimension == 2)) {
- if constexpr (MeshType::Dimension == 2) {
- SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek = inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[t];
- SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
- SmallArray<double>& mean_i_of_ejk = mean_i_of_ejk_pool[t];
-
- auto face_to_node_matrix = p_mesh->connectivity().faceToNodeMatrix();
- auto cell_face_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
- const Rd& Xj = xj[cell_j_id];
-
- _computeEjkMeanByBoundary(*p_mesh, Xj, cell_j_id, cell_to_face_matrix, face_to_node_matrix,
- cell_face_is_reversed, Vj, m_descriptor.degree(), basis_dimension,
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_j_of_ejk);
-
- size_t index = 0;
- for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = stencil_cell_list[i];
-
- _computeEjkMeanByBoundary(*p_mesh, Xj, cell_i_id, cell_to_face_matrix, face_to_node_matrix,
- cell_face_is_reversed, Vj, m_descriptor.degree(), basis_dimension,
- inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_i_of_ejk);
-
- for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
- }
- }
- for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
- ++i_symmetry) {
- auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
- auto ghost_cell_list = ghost_stencil[cell_j_id];
-
- const Rd& origin = symmetry_origin_list[i_symmetry];
- const Rd& normal = symmetry_normal_list[i_symmetry];
-
- for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = ghost_cell_list[i];
-
- _computeEjkMeanByBoundaryInSymmetricCell(origin, normal, //
- Xj, cell_i_id, xr, cell_to_face_matrix, face_to_node_matrix,
- cell_face_is_reversed, Vj, m_descriptor.degree(),
- basis_dimension, inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
- mean_i_of_ejk);
-
- for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
- }
- }
- }
- } else {
- throw UnexpectedError("invalid mesh dimension");
- }
- } else {
- SmallArray<double>& inv_Vj_wq_detJ_ek = inv_Vj_wq_detJ_ek_pool[t];
- SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
- SmallArray<double>& mean_i_of_ejk = mean_i_of_ejk_pool[t];
-
- const Rd& Xj = xj[cell_j_id];
-
- _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);
-
- size_t index = 0;
- for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = stencil_cell_list[i];
-
- _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(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
- }
- }
- for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
- ++i_symmetry) {
- auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
- auto ghost_cell_list = ghost_stencil[cell_j_id];
-
- const Rd& origin = symmetry_origin_list[i_symmetry];
- const Rd& normal = symmetry_normal_list[i_symmetry];
-
- for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = ghost_cell_list[i];
-
- _computeEjkMeanInSymmetricCell(origin, normal, //
- 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);
+ Internal<MeshType>::buildB(cell_j_id, stencil_array, discrete_function_variant_list, symmetry_normal_list, B);
- for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
- }
- }
- }
- }
- } else {
- throw UnexpectedError("invalid integration strategy");
- }
+ ReconstructionMatrixBuilderType& reconstruction_matrix_builder = *reconstruction_matrix_builder_pool[t];
+ reconstruction_matrix_builder.build(cell_j_id, A);
if (m_descriptor.rowWeighting()) {
- // Add row weighting (give more importance to cells that are
- // closer to j)
- const Rd& Xj = xj[cell_j_id];
-
- size_t index = 0;
- for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = stencil_cell_list[i];
- const double weight = 1. / l2Norm(xj[cell_i_id] - Xj);
- for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(index, l) *= weight;
- }
- for (size_t l = 0; l < number_of_columns; ++l) {
- B(index, l) *= weight;
- }
- }
- for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
- ++i_symmetry) {
- auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
- auto ghost_cell_list = ghost_stencil[cell_j_id];
-
- const Rd& origin = symmetry_origin_list[i_symmetry];
- const Rd& normal = symmetry_normal_list[i_symmetry];
-
- for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
- const CellId cell_i_id = ghost_cell_list[i];
- const double weight = 1. / l2Norm(symmetrize_coordinates(origin, normal, xj[cell_i_id]) - Xj);
- for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(index, l) *= weight;
- }
- for (size_t l = 0; l < number_of_columns; ++l) {
- B(index, l) *= weight;
- }
- }
- }
+ Internal<MeshType>::rowWeighting(cell_j_id, stencil_array, xj, symmetry_origin_list, symmetry_normal_list, A,
+ B);
}
const SmallMatrix<double>& X = X_pool[t];
if (m_descriptor.preconditioning()) {
- // Add column weighting preconditioning (increase the presition)
+ // Add column weighting preconditioning (increase the precision)
SmallVector<double>& G = G_pool[t];
- for (size_t l = 0; l < A.numberOfColumns(); ++l) {
- double g = 0;
- for (size_t i = 0; i < A.numberOfRows(); ++i) {
- const double Ail = A(i, l);
-
- g += Ail * Ail;
- }
- G[l] = std::sqrt(g);
- }
-
- for (size_t l = 0; l < A.numberOfColumns(); ++l) {
- const double Gl = G[l];
- for (size_t i = 0; i < A.numberOfRows(); ++i) {
- A(i, l) *= Gl;
- }
- }
-
- Givens::solveCollectionInPlace(A, X, B);
-
- for (size_t l = 0; l < X.numberOfRows(); ++l) {
- const double Gl = G[l];
- for (size_t i = 0; i < X.numberOfColumns(); ++i) {
- X(l, i) *= Gl;
- }
- }
+ Internal<MeshType>::solveCollectionInPlaceWithPreconditionner(A, X, B, G);
} else {
Givens::solveCollectionInPlace(A, X, B);
}
- column_begin = 0;
- for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size();
- ++i_dpk_variant) {
- const auto& dpk_variant = mutable_discrete_function_dpk_variant_list[i_dpk_variant];
-
- const auto& discrete_function_variant = discrete_function_variant_list[i_dpk_variant];
-
- std::visit(
- [&](auto&& dpk_function, auto&& p0_function) {
- using DPkFunctionT = std::decay_t<decltype(dpk_function)>;
- using P0FunctionT = std::decay_t<decltype(p0_function)>;
- using DataType = std::remove_const_t<std::decay_t<typename DPkFunctionT::data_type>>;
- using P0DataType = std::remove_const_t<std::decay_t<typename P0FunctionT::data_type>>;
-
- if constexpr (std::is_same_v<DataType, P0DataType>) {
- if constexpr (is_discrete_function_P0_v<P0FunctionT>) {
- if constexpr (is_discrete_function_dpk_scalar_v<DPkFunctionT>) {
- auto dpk_j = dpk_function.coefficients(cell_j_id);
- dpk_j[0] = p0_function[cell_j_id];
-
- if constexpr (std::is_arithmetic_v<DataType>) {
- 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) {
- dpk_j[0] -= X(i, column_begin) * mean_j_of_ejk[i];
- }
- }
-
- for (size_t i = 0; i < basis_dimension - 1; ++i) {
- auto& dpk_j_ip1 = dpk_j[i + 1];
- dpk_j_ip1 = X(i, column_begin);
- }
- ++column_begin;
- } else if constexpr (is_tiny_vector_v<DataType>) {
- 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];
- for (size_t k = 0; k < DataType::Dimension; ++k) {
- dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
- }
- }
- }
-
- for (size_t i = 0; i < basis_dimension - 1; ++i) {
- auto& dpk_j_ip1 = dpk_j[i + 1];
- for (size_t k = 0; k < DataType::Dimension; ++k) {
- dpk_j_ip1[k] = X(i, column_begin + k);
- }
- }
- column_begin += DataType::Dimension;
- } else if constexpr (is_tiny_matrix_v<DataType>) {
- 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];
- for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- dpk_j_0(k, l) -=
- X(i, column_begin + k * DataType::NumberOfColumns + l) * mean_j_of_ejk[i];
- }
- }
- }
- }
-
- for (size_t i = 0; i < basis_dimension - 1; ++i) {
- auto& dpk_j_ip1 = dpk_j[i + 1];
- for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
- }
- }
- }
- column_begin += DataType::Dimension;
- } else {
- // LCOV_EXCL_START
- throw UnexpectedError("unexpected data type");
- // LCOV_EXCL_STOP
- }
- } else {
- // LCOV_EXCL_START
- throw UnexpectedError("unexpected discrete dpk function type");
- // LCOV_EXCL_STOP
- }
- } else if constexpr (is_discrete_function_P0_vector_v<P0FunctionT>) {
- if constexpr (is_discrete_function_dpk_vector_v<DPkFunctionT>) {
- auto dpk_j = dpk_function.coefficients(cell_j_id);
- auto cell_vector = p0_function[cell_j_id];
- const size_t size = basis_dimension;
-
- for (size_t l = 0; l < cell_vector.size(); ++l) {
- const size_t component_begin = l * size;
- dpk_j[component_begin] = cell_vector[l];
- if constexpr (std::is_arithmetic_v<DataType>) {
- 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) {
- dpk_j[component_begin] -= X(i, column_begin) * mean_j_of_ejk[i];
- }
- }
-
- for (size_t i = 0; i < basis_dimension - 1; ++i) {
- auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
- dpk_j_ip1 = X(i, column_begin);
- }
- ++column_begin;
- } else if constexpr (is_tiny_vector_v<DataType>) {
- 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[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];
- }
- }
- }
-
- for (size_t i = 0; i < basis_dimension - 1; ++i) {
- auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
- for (size_t k = 0; k < DataType::Dimension; ++k) {
- dpk_j_ip1[k] = X(i, column_begin + k);
- }
- }
- column_begin += DataType::Dimension;
- } else {
- // LCOV_EXCL_START
- throw UnexpectedError("unexpected data type");
- // LCOV_EXCL_STOP
- }
- }
- } else {
- // LCOV_EXCL_START
- throw UnexpectedError("unexpected discrete dpk function type");
- // LCOV_EXCL_STOP
- }
- } else {
- // LCOV_EXCL_START
- throw UnexpectedError("unexpected discrete function type");
- // LCOV_EXCL_STOP
- }
- } else {
- // LCOV_EXCL_START
- throw UnexpectedError("incompatible data types");
- // LCOV_EXCL_STOP
- }
- },
- dpk_variant.mutableDiscreteFunctionDPk(), discrete_function_variant->discreteFunction());
- }
+ Internal<MeshType>::template populateDiscreteFunctionDPkByCell(cell_j_id, m_descriptor.degree(), X,
+ reconstruction_matrix_builder,
+ discrete_function_variant_list,
+ mutable_discrete_function_dpk_variant_list);
tokens.release(t);
}
@@ -1509,6 +813,30 @@ PolynomialReconstruction::_build(
return discrete_function_dpk_variant_list;
}
+template <MeshConcept MeshType>
+[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
+PolynomialReconstruction::_build(
+ const std::shared_ptr<const MeshType>& p_mesh,
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+ switch (m_descriptor.integrationMethodType()) {
+ case IntegrationMethodType::cell_center: {
+ return this->_build<CellCenterReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
+ }
+ case IntegrationMethodType::boundary: {
+ return this->_build<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
+ }
+ case IntegrationMethodType::element: {
+ return this->_build<ElementIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("invalid reconstruction matrix builder type");
+ }
+ // LCOV_EXCL_STOP
+ }
+}
+
std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::build(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index bf6013639f732fd52f9c145ba4982701c0e4bcfc..4f76844711ada4641b69234c5ac787c472247a37 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -10,8 +10,20 @@ class DiscreteFunctionVariant;
class PolynomialReconstruction
{
private:
+ template <MeshConcept MeshType>
+ class Internal;
+
class MutableDiscreteFunctionDPkVariant;
+ template <MeshConcept MeshType>
+ class CellCenterReconstructionMatrixBuilder;
+
+ template <MeshConcept MeshType>
+ class ElementIntegralReconstructionMatrixBuilder;
+
+ template <MeshConcept MeshType>
+ class BoundaryIntegralReconstructionMatrixBuilder;
+
const PolynomialReconstructionDescriptor m_descriptor;
size_t _getNumberOfColumns(
@@ -26,6 +38,11 @@ class PolynomialReconstruction
const std::shared_ptr<const MeshType>& p_mesh,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+ template <typename ReconstructionMatrixBuilderType, MeshConcept MeshType>
+ [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build(
+ const std::shared_ptr<const MeshType>& p_mesh,
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;
+
template <MeshConcept MeshType>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build(
const std::shared_ptr<const MeshType>& p_mesh,
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..58b8fcd5ea34a64bd8d5808c41ea95079c2fcf84
--- /dev/null
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
@@ -0,0 +1,553 @@
+#include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureFormula.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <geometry/SquareTransformation.hpp>
+#include <geometry/SymmetryUtils.hpp>
+#include <geometry/TriangleTransformation.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+
+template <MeshConcept MeshTypeT>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::SpecificDimensionalData
+{
+ private:
+ const ItemToItemMatrix<ItemType::cell, ItemType::face> m_cell_to_face_matrix;
+ const ItemToItemMatrix<ItemType::face, ItemType::node> m_face_to_node_matrix;
+ const FaceValuePerCell<const bool> m_cell_face_is_reversed;
+
+ public:
+ PUGS_INLINE const auto&
+ cellToFaceMatrix() const noexcept
+ {
+ return m_cell_to_face_matrix;
+ }
+
+ PUGS_INLINE
+ const auto&
+ faceToNodeMatrix() const noexcept
+ {
+ return m_face_to_node_matrix;
+ }
+
+ PUGS_INLINE
+ const auto&
+ cellFaceIsReversed() const noexcept
+ {
+ return m_cell_face_is_reversed;
+ }
+
+ SpecificDimensionalData(const Connectivity<MeshType::Dimension>& connectivity)
+ : m_cell_to_face_matrix{connectivity.cellToFaceMatrix()},
+ m_face_to_node_matrix{connectivity.faceToNodeMatrix()},
+ m_cell_face_is_reversed{connectivity.cellFaceIsReversed()}
+ {}
+
+ ~SpecificDimensionalData() = default;
+};
+
+template <>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<1>>::SpecificDimensionalData
+{
+ private:
+ const ItemToItemMatrix<ItemType::cell, ItemType::node> m_cell_to_node_matrix;
+
+ public:
+ PUGS_INLINE
+ const auto&
+ cellToNodeMatrix() const noexcept
+ {
+ return m_cell_to_node_matrix;
+ }
+
+ SpecificDimensionalData(const Connectivity<1>& connectivity) : m_cell_to_node_matrix{connectivity.cellToNodeMatrix()}
+ {}
+
+ ~SpecificDimensionalData() = default;
+};
+
+template <>
+template <typename ConformTransformationT>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<3>>::_computeEjkBoundaryMean(
+ const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
+ const ConformTransformationT& T,
+ const Rd& Xj,
+ const double inv_Vi,
+ SmallArray<double>& mean_of_ejk)
+{
+ for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+ const double wq = quadrature.weight(i_q);
+ const TinyVector<2> xi_q = quadrature.point(i_q);
+
+ const double area_variation_e1 = T.areaVariation(xi_q)[0] * inv_Vi;
+
+ 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;
+ m_tmp_array[k++] = x_xj * wq * area_variation_e1;
+ for (; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = x_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
+ const size_t begin_i_y_1 = m_yz_row_index[i_y - 1];
+ const size_t nb_monoms = m_yz_row_size[i_y];
+ for (size_t l = 0; l < nb_monoms; ++l, ++k) {
+ m_tmp_array[k] = y_yj * m_tmp_array[begin_i_y_1 + l];
+ }
+ }
+
+ for (size_t i_z = 1; i_z <= m_polynomial_degree; ++i_z) {
+ const size_t nb_y = m_yz_row_size[m_z_triangle_index[i_z]];
+ const size_t index_z = m_z_triangle_index[i_z];
+ const size_t index_z_1 = m_z_triangle_index[i_z - 1];
+ for (size_t i_y = 0; i_y < nb_y; ++i_y) {
+ const size_t begin_i_yz_1 = m_yz_row_index[index_z_1 + i_y];
+ const size_t nb_monoms = m_yz_row_size[index_z + i_y];
+ for (size_t l = 0; l < nb_monoms; ++l, ++k) {
+ m_tmp_array[k] = z_zj * m_tmp_array[begin_i_yz_1 + l];
+ }
+ }
+ }
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] += m_tmp_array[k];
+ }
+ }
+}
+
+template <>
+template <typename ConformTransformationT>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::_computeEjkBoundaryMean(
+ const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
+ const ConformTransformationT& T,
+ const Rd& Xj,
+ const double inv_Vi,
+ SmallArray<double>& mean_of_ejk)
+{
+ const double velocity_perp_e1 = T.velocity()[1] * inv_Vi;
+
+ for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+ const double wq = quadrature.weight(i_q);
+ const TinyVector<1> xi_q = quadrature.point(i_q);
+
+ const Rd X_Xj = T(xi_q) - Xj;
+
+ const double x_xj = X_Xj[0];
+ const double y_yj = X_Xj[1];
+
+ {
+ size_t k = 0;
+ m_tmp_array[k++] = x_xj * wq * velocity_perp_e1;
+ for (; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = x_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
+ const size_t begin_i_y_1 = m_y_row_index[i_y - 1];
+ for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l) {
+ m_tmp_array[k++] = y_yj * m_tmp_array[begin_i_y_1 + l];
+ }
+ }
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] += m_tmp_array[k];
+ }
+ }
+}
+
+template <MeshConcept MeshTypeT>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkMeanByBoundary(
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ const double inv_Vi = 1. / m_Vj[cell_id];
+
+ const auto& face_is_reversed = m_dimensional_data_ptr->cellFaceIsReversed()[cell_id];
+ const auto& cell_face_list = m_dimensional_data_ptr->cellToFaceMatrix()[cell_id];
+ const auto& face_to_node_matrix = m_dimensional_data_ptr->faceToNodeMatrix();
+
+ mean_of_ejk.fill(0);
+
+ if constexpr (MeshType::Dimension == 2) {
+ const auto& quadrature =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+ if (face_is_reversed[i_face]) {
+ const LineTransformation<2> T{m_xr[face_node_list[1]], m_xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const LineTransformation<2> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ }
+ } else {
+ static_assert(MeshType::Dimension == 3);
+
+ const auto& triangle_quadrature =
+ QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(m_polynomial_degree + 1));
+ const auto& square_quadrature =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+ switch (face_node_list.size()) {
+ case 3: {
+ if (face_is_reversed[i_face]) {
+ const TriangleTransformation<3> T{m_xr[face_node_list[2]], m_xr[face_node_list[1]], m_xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const TriangleTransformation<3> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]], m_xr[face_node_list[2]]};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ case 4: {
+ if (face_is_reversed[i_face]) {
+ const SquareTransformation<3> T{m_xr[face_node_list[3]], m_xr[face_node_list[2]], m_xr[face_node_list[1]],
+ m_xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const SquareTransformation<3> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]], m_xr[face_node_list[2]],
+ m_xr[face_node_list[3]]};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ default: {
+ throw NotImplementedError("invalid face type");
+ }
+ }
+ }
+ }
+}
+
+template <MeshConcept MeshTypeT>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ MeshTypeT>::_computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
+ const Rd& normal,
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ const double inv_Vi = 1. / m_Vj[cell_id];
+
+ const auto& face_is_reversed = m_dimensional_data_ptr->cellFaceIsReversed()[cell_id];
+ const auto& cell_face_list = m_dimensional_data_ptr->cellToFaceMatrix()[cell_id];
+ const auto& face_to_node_matrix = m_dimensional_data_ptr->faceToNodeMatrix();
+
+ mean_of_ejk.fill(0);
+ if constexpr (MeshType::Dimension == 2) {
+ const auto& quadrature =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
+
+ if (face_is_reversed[i_face]) {
+ const LineTransformation<2> T{x1, x0};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const LineTransformation<2> T{x0, x1};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ }
+ } else {
+ static_assert(MeshType::Dimension == 3);
+
+ const auto& triangle_quadrature =
+ QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(m_polynomial_degree + 1));
+ const auto& square_quadrature =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+ switch (face_node_list.size()) {
+ case 3: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[2]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
+
+ if (face_is_reversed[i_face]) {
+ const TriangleTransformation<3> T{x2, x1, x0};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const TriangleTransformation<3> T{x0, x1, x2};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ case 4: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[3]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[2]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
+
+ if (face_is_reversed[i_face]) {
+ const SquareTransformation<3> T{x3, x2, x1, x0};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const SquareTransformation<3> T{x0, x1, x2, x3};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ default: {
+ throw NotImplementedError("invalid face type");
+ }
+ }
+ }
+ }
+}
+
+template <>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<1>>::_computeEjkMeanByBoundary(
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ const double inv_Vi = 1. / m_Vj[cell_id];
+
+ const auto& cell_node_list = m_dimensional_data_ptr->cellToNodeMatrix()[cell_id];
+
+ const double xr1_xj = (m_xr[cell_node_list[1]] - Xj)[0];
+
+ m_tmp_array[0] = xr1_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr1_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] = m_tmp_array[k];
+ }
+
+ const double xr0_xj = (m_xr[cell_node_list[0]] - Xj)[0];
+
+ m_tmp_array[0] = xr0_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr0_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] -= m_tmp_array[k];
+ }
+}
+
+template <>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<1>>::_computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
+ const Rd& normal,
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ const double inv_Vi = 1. / m_Vj[cell_id];
+
+ const auto& cell_node_list = m_dimensional_data_ptr->cellToNodeMatrix()[cell_id];
+
+ const double xr1_xj = (symmetrize_coordinates(origin, normal, m_xr[cell_node_list[0]]) - Xj)[0];
+
+ m_tmp_array[0] = xr1_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr1_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] = m_tmp_array[k];
+ }
+
+ const double xr0_xj = (symmetrize_coordinates(origin, normal, m_xr[cell_node_list[1]]) - Xj)[0];
+
+ m_tmp_array[0] = xr0_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr0_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] -= m_tmp_array[k];
+ }
+}
+
+template <MeshConcept MeshTypeT>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::build(
+ const CellId cell_j_id,
+ ShrinkMatrixView<SmallMatrix<double>>& A)
+{
+ if constexpr (MeshType::Dimension <= 3) {
+ const auto& stencil_cell_list = m_stencil_array[cell_j_id];
+
+ const Rd& Xj = m_xj[cell_j_id];
+
+ _computeEjkMeanByBoundary(Xj, cell_j_id, m_mean_j_of_ejk);
+
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = stencil_cell_list[i];
+
+ _computeEjkMeanByBoundary(Xj, cell_i_id, m_mean_i_of_ejk);
+
+ for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
+ A(index, l) = m_mean_i_of_ejk[l] - m_mean_j_of_ejk[l];
+ }
+ }
+ for (size_t i_symmetry = 0; i_symmetry < m_stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry) {
+ auto& ghost_stencil = m_stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+ const Rd& origin = m_symmetry_origin_list[i_symmetry];
+ const Rd& normal = m_symmetry_normal_list[i_symmetry];
+
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+
+ _computeEjkMeanByBoundaryInSymmetricCell(origin, normal, Xj, cell_i_id, m_mean_i_of_ejk);
+
+ for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
+ A(index, l) = m_mean_i_of_ejk[l] - m_mean_j_of_ejk[l];
+ }
+ }
+ }
+ } else {
+ throw NotImplementedError("invalid mesh dimension");
+ }
+}
+
+template <MeshConcept MeshTypeT>
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ MeshTypeT>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh,
+ const size_t polynomial_degree,
+ const SmallArray<const Rd>& symmetry_origin_list,
+ const SmallArray<const Rd>& symmetry_normal_list,
+ const CellToCellStencilArray& stencil_array)
+ : m_mesh{mesh},
+ m_basis_dimension{
+ DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(polynomial_degree)},
+ m_polynomial_degree{polynomial_degree},
+ m_tmp_array{m_basis_dimension},
+ m_mean_j_of_ejk{m_basis_dimension - 1},
+ m_mean_i_of_ejk{m_basis_dimension - 1},
+
+ m_dimensional_data_ptr{std::make_shared<SpecificDimensionalData>(mesh.connectivity())},
+
+ m_stencil_array{stencil_array},
+ m_symmetry_origin_list{symmetry_origin_list},
+ m_symmetry_normal_list{symmetry_normal_list},
+ m_Vj{MeshDataManager::instance().getMeshData(mesh).Vj()},
+ m_xj{MeshDataManager::instance().getMeshData(mesh).xj()},
+ m_xr{mesh.xr()}
+{
+ SmallArray<double> antiderivative_correction_coef(m_polynomial_degree + 1);
+ for (size_t k = 0; k < antiderivative_correction_coef.size(); ++k) {
+ // The antiderivative of x^k is k/(k+1) times the antiderivative of x^(k-1)
+ antiderivative_correction_coef[k] = ((1. * k) / (k + 1));
+ }
+
+ m_antiderivative_correction_coef = antiderivative_correction_coef;
+
+ if constexpr (MeshType::Dimension == 2) {
+ SmallArray<size_t> y_row_index(m_polynomial_degree + 1);
+
+ size_t i_y = 0;
+
+ y_row_index[i_y++] = 0;
+ for (ssize_t n = m_polynomial_degree + 1; n > 1; --n, ++i_y) {
+ y_row_index[i_y] = y_row_index[i_y - 1] + n;
+ }
+
+ m_y_row_index = y_row_index;
+
+ } else if constexpr (MeshType::Dimension == 3) {
+ SmallArray<size_t> yz_row_index((m_polynomial_degree + 2) * (m_polynomial_degree + 1) / 2 + 1);
+ SmallArray<size_t> z_triangle_index(m_polynomial_degree + 1);
+
+ {
+ size_t i_z = 0;
+ size_t i_yz = 0;
+
+ yz_row_index[i_yz++] = 0;
+ for (ssize_t n = m_polynomial_degree + 1; n >= 1; --n) {
+ z_triangle_index[i_z++] = i_yz - 1;
+ for (ssize_t i = n; i >= 1; --i) {
+ yz_row_index[i_yz] = yz_row_index[i_yz - 1] + i;
+ ++i_yz;
+ }
+ }
+ }
+
+ SmallArray<size_t> yz_row_size{yz_row_index.size() - 1};
+ for (size_t i = 0; i < yz_row_size.size(); ++i) {
+ yz_row_size[i] = yz_row_index[i + 1] - yz_row_index[i];
+ }
+
+ m_yz_row_index = yz_row_index;
+ m_z_triangle_index = z_triangle_index;
+ m_yz_row_size = yz_row_size;
+ }
+}
+
+template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<1>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<3>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<1>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
+
+template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<2>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
+
+template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<3>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..60707189b909e9d5f17ae2852be655cbfb877351
--- /dev/null
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
@@ -0,0 +1,89 @@
+#ifndef BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
+#define BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
+
+#include <algebra/ShrinkMatrixView.hpp>
+#include <algebra/SmallMatrix.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/StencilArray.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/SmallArray.hpp>
+
+template <size_t Dimension>
+class QuadratureFormula;
+
+template <MeshConcept MeshTypeT>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
+{
+ public:
+ using MeshType = MeshTypeT;
+
+ constexpr static bool handles_high_degrees = true;
+
+ private:
+ using Rd = TinyVector<MeshType::Dimension>;
+
+ const MeshType& m_mesh;
+ const size_t m_basis_dimension;
+ const size_t m_polynomial_degree;
+
+ const SmallArray<double> m_tmp_array;
+ SmallArray<double> m_mean_j_of_ejk;
+ SmallArray<double> m_mean_i_of_ejk;
+
+ class SpecificDimensionalData;
+ std::shared_ptr<SpecificDimensionalData> m_dimensional_data_ptr;
+
+ const CellToCellStencilArray& m_stencil_array;
+
+ const SmallArray<const Rd> m_symmetry_origin_list;
+ const SmallArray<const Rd> m_symmetry_normal_list;
+
+ const CellValue<const double> m_Vj;
+ const CellValue<const Rd> m_xj;
+ const NodeValue<const Rd> m_xr;
+
+ // 2D
+ SmallArray<const size_t> m_y_row_index;
+
+ // 3D
+ SmallArray<const size_t> m_yz_row_index;
+ SmallArray<const size_t> m_z_triangle_index;
+ SmallArray<const size_t> m_yz_row_size;
+
+ SmallArray<const double> m_antiderivative_correction_coef;
+
+ template <typename ConformTransformationT>
+ void _computeEjkBoundaryMean(const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
+ const ConformTransformationT& T,
+ const Rd& Xj,
+ const double inv_Vi,
+ SmallArray<double>& mean_of_ejk);
+
+ void _computeEjkMeanByBoundary(const Rd& Xj, const CellId& cell_id, SmallArray<double>& mean_of_ejk);
+
+ void _computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
+ const Rd& normal,
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk);
+
+ public:
+ PUGS_INLINE
+ SmallArray<const double>
+ meanjOfEjk() const
+ {
+ return m_mean_j_of_ejk;
+ }
+
+ void build(const CellId cell_j_id, ShrinkMatrixView<SmallMatrix<double>>& A);
+
+ BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh,
+ const size_t polynomial_degree,
+ const SmallArray<const Rd>& symmetry_origin_list,
+ const SmallArray<const Rd>& symmetry_normal_list,
+ const CellToCellStencilArray& stencil_array);
+
+ ~BoundaryIntegralReconstructionMatrixBuilder() = default;
+};
+
+#endif // BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
diff --git a/src/scheme/reconstruction_utils/CMakeLists.txt b/src/scheme/reconstruction_utils/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1d4ebcb48c3d47c9c19421151257deb8a31bf9a9
--- /dev/null
+++ b/src/scheme/reconstruction_utils/CMakeLists.txt
@@ -0,0 +1,17 @@
+# ------------------- Source files --------------------
+
+add_library(PugsSchemeReconstructionUtils
+ BoundaryIntegralReconstructionMatrixBuilder.cpp
+ CellCenterReconstructionMatrixBuilder.cpp
+ ElementIntegralReconstructionMatrixBuilder.cpp
+)
+
+add_dependencies(
+ PugsUtils
+ PugsMesh
+)
+
+target_link_libraries(
+ PugsSchemeReconstructionUtils
+ ${HIGHFIVE_TARGET}
+)
diff --git a/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..94d1ef9a79ac1fb8f7e3cc0ece4c8eb21672bef7
--- /dev/null
+++ b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.cpp
@@ -0,0 +1,92 @@
+#include <scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp>
+
+#include <geometry/SymmetryUtils.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+
+template <MeshConcept MeshTypeT>
+void
+PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<MeshTypeT>::build(
+ const CellId cell_j_id,
+ ShrinkMatrixView<SmallMatrix<double>>& A)
+{
+ const auto& stencil_cell_list = m_stencil_array[cell_j_id];
+
+ const Rd& Xj = m_xj[cell_j_id];
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const Rd Xi_Xj = m_xj[cell_i_id] - Xj;
+ for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
+ A(index, l) = Xi_Xj[l];
+ }
+ }
+ for (size_t i_symmetry = 0; i_symmetry < m_stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry) {
+ auto& ghost_stencil = m_stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+ const Rd& origin = m_symmetry_origin_list[i_symmetry];
+ const Rd& normal = m_symmetry_normal_list[i_symmetry];
+
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+ const Rd Xi_Xj = symmetrize_coordinates(origin, normal, m_xj[cell_i_id]) - Xj;
+ for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
+ A(index, l) = Xi_Xj[l];
+ }
+ }
+ }
+}
+
+template <MeshConcept MeshTypeT>
+PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<MeshTypeT>::CellCenterReconstructionMatrixBuilder(
+ const MeshType& mesh,
+ const size_t polynomial_degree,
+ const SmallArray<const Rd>& symmetry_origin_list,
+ const SmallArray<const Rd>& symmetry_normal_list,
+ const CellToCellStencilArray& stencil_array)
+ : m_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(
+ polynomial_degree)},
+ m_symmetry_origin_list{symmetry_origin_list},
+ m_symmetry_normal_list{symmetry_normal_list},
+ m_stencil_array{stencil_array},
+ m_xj{MeshDataManager::instance().getMeshData(mesh).xj()}
+{
+ if (polynomial_degree != 1) {
+ throw NormalError("cell center based reconstruction is only valid for first order");
+ }
+}
+
+template void PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<Mesh<1>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template void PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<Mesh<2>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template void PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<Mesh<3>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<
+ Mesh<1>>::CellCenterReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
+
+template PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<
+ Mesh<2>>::CellCenterReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
+
+template PolynomialReconstruction::CellCenterReconstructionMatrixBuilder<
+ Mesh<3>>::CellCenterReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
diff --git a/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a8205a625dc0d62e35a33469fe2864719f8147a6
--- /dev/null
+++ b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp
@@ -0,0 +1,41 @@
+#ifndef CELL_CENTER_RECONSTRUCTION_MATRIX_BUILDER_HPP
+#define CELL_CENTER_RECONSTRUCTION_MATRIX_BUILDER_HPP
+
+#include <algebra/ShrinkMatrixView.hpp>
+#include <algebra/SmallMatrix.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/StencilArray.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/SmallArray.hpp>
+
+template <MeshConcept MeshTypeT>
+class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder
+{
+ public:
+ using MeshType = MeshTypeT;
+
+ constexpr static bool handles_high_degrees = false;
+
+ private:
+ using Rd = TinyVector<MeshType::Dimension>;
+
+ const size_t m_basis_dimension;
+ const SmallArray<const Rd> m_symmetry_origin_list;
+ const SmallArray<const Rd> m_symmetry_normal_list;
+
+ const CellToCellStencilArray& m_stencil_array;
+ const CellValue<const Rd> m_xj;
+
+ public:
+ void build(const CellId cell_j_id, ShrinkMatrixView<SmallMatrix<double>>& A);
+
+ CellCenterReconstructionMatrixBuilder(const MeshType& mesh,
+ const size_t polynomial_degree,
+ const SmallArray<const Rd>& symmetry_origin_list,
+ const SmallArray<const Rd>& symmetry_normal_list,
+ const CellToCellStencilArray& stencil_array);
+
+ ~CellCenterReconstructionMatrixBuilder() = default;
+};
+
+#endif // CELL_CENTER_RECONSTRUCTION_MATRIX_BUILDER_HPP
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..09e410811ed203610926bd605be3618e4219f396
--- /dev/null
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
@@ -0,0 +1,501 @@
+#include <scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.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/SymmetryUtils.hpp>
+#include <geometry/TetrahedronTransformation.hpp>
+#include <geometry/TriangleTransformation.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+
+template <MeshConcept MeshTypeT>
+template <typename ConformTransformationT>
+void
+PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkMean(
+ const QuadratureFormula<MeshType::Dimension>& quadrature,
+ const ConformTransformationT& T,
+ const Rd& Xj,
+ const double Vi,
+ SmallArray<double>& mean_of_ejk) noexcept(NO_ASSERT)
+{
+ 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 Rd X_Xj = T(xi_q) - Xj;
+
+ if constexpr (MeshType::Dimension == 1) {
+ const double detT = T.jacobianDeterminant();
+
+ const double x_xj = X_Xj[0];
+
+ {
+ m_wq_detJ_ek[0] = wq * detT;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_wq_detJ_ek[k] = x_xj * m_wq_detJ_ek[k - 1];
+ }
+ }
+
+ } else if constexpr (MeshType::Dimension == 2) {
+ const double detT = [&] {
+ if constexpr (std::is_same_v<TriangleTransformation<2>, std::decay_t<decltype(T)>>) {
+ return T.jacobianDeterminant();
+ } else {
+ return T.jacobianDeterminant(xi_q);
+ }
+ }();
+
+ const double x_xj = X_Xj[0];
+ const double y_yj = X_Xj[1];
+
+ {
+ size_t k = 0;
+ m_wq_detJ_ek[k++] = wq * detT;
+ for (; k <= m_polynomial_degree; ++k) {
+ m_wq_detJ_ek[k] = x_xj * m_wq_detJ_ek[k - 1];
+ }
+
+ for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
+ const size_t begin_i_y_1 = m_y_row_index[i_y - 1];
+ for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l, ++k) {
+ m_wq_detJ_ek[k] = y_yj * m_wq_detJ_ek[begin_i_y_1 + l];
+ }
+ }
+ }
+
+ } else if constexpr (MeshType::Dimension == 3) {
+ static_assert(MeshType::Dimension == 3);
+
+ const double detT = T.jacobianDeterminant(xi_q);
+
+ 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;
+ m_wq_detJ_ek[k++] = wq * detT;
+ for (; k <= m_polynomial_degree; ++k) {
+ m_wq_detJ_ek[k] = x_xj * m_wq_detJ_ek[k - 1];
+ }
+
+ for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
+ const size_t begin_i_y_1 = m_yz_row_index[i_y - 1];
+ const size_t nb_monoms = m_yz_row_size[i_y];
+ for (size_t l = 0; l < nb_monoms; ++l, ++k) {
+ m_wq_detJ_ek[k] = y_yj * m_wq_detJ_ek[begin_i_y_1 + l];
+ }
+ }
+
+ for (size_t i_z = 1; i_z <= m_polynomial_degree; ++i_z) {
+ const size_t nb_y = m_yz_row_size[m_z_triangle_index[i_z]];
+ const size_t index_z = m_z_triangle_index[i_z];
+ const size_t index_z_1 = m_z_triangle_index[i_z - 1];
+ for (size_t i_y = 0; i_y < nb_y; ++i_y) {
+ const size_t begin_i_yz_1 = m_yz_row_index[index_z_1 + i_y];
+ const size_t nb_monoms = m_yz_row_size[index_z + i_y];
+ for (size_t l = 0; l < nb_monoms; ++l, ++k) {
+ m_wq_detJ_ek[k] = z_zj * m_wq_detJ_ek[begin_i_yz_1 + l];
+ }
+ }
+ }
+ }
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] += m_wq_detJ_ek[k];
+ }
+ }
+
+ const double inv_Vi = 1. / Vi;
+ for (size_t k = 0; k < mean_of_ejk.size(); ++k) {
+ mean_of_ejk[k] *= inv_Vi;
+ }
+}
+
+template <MeshConcept MeshTypeT>
+void
+PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkMean(
+ const Rd& Xj,
+ const CellId& cell_i_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ const CellType cell_type = m_cell_type[cell_i_id];
+ const auto node_list = m_cell_to_node_matrix[cell_i_id];
+ const double Vi = m_Vj[cell_i_id];
+
+ if constexpr (MeshType::Dimension == 1) {
+ if (m_cell_type[cell_i_id] == CellType::Line) {
+ const LineTransformation<1> T{m_xr[node_list[0]], m_xr[node_list[1]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+ } else {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+ } else if constexpr (MeshType::Dimension == 2) {
+ switch (cell_type) {
+ case CellType::Triangle: {
+ const TriangleTransformation<2> T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]]};
+ const auto& quadrature =
+ QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{m_polynomial_degree});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+ break;
+ }
+ case CellType::Quadrangle: {
+ const SquareTransformation<2> T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]], m_xr[node_list[3]]};
+ const auto& quadrature =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree + 1});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+ break;
+ }
+ default: {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+ }
+ } else {
+ static_assert(MeshType::Dimension == 3);
+
+ switch (cell_type) {
+ case CellType::Tetrahedron: {
+ const TetrahedronTransformation T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]], m_xr[node_list[3]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{m_polynomial_degree});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+ break;
+ }
+ case CellType::Prism: {
+ const PrismTransformation T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]], //
+ m_xr[node_list[3]], m_xr[node_list[4]], m_xr[node_list[5]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{m_polynomial_degree + 1});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+ break;
+ }
+ case CellType::Pyramid: {
+ const PyramidTransformation T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]], m_xr[node_list[3]],
+ m_xr[node_list[4]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{m_polynomial_degree + 1});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+ break;
+ }
+ case CellType::Hexahedron: {
+ const CubeTransformation T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]], m_xr[node_list[3]],
+ m_xr[node_list[4]], m_xr[node_list[5]], m_xr[node_list[6]], m_xr[node_list[7]]};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree + 1});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+ break;
+ }
+ default: {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+ }
+ }
+}
+
+template <MeshConcept MeshTypeT>
+void
+PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkMeanInSymmetricCell(
+ const Rd& origin,
+ const Rd& normal,
+ const Rd& Xj,
+ const CellId& cell_i_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ if constexpr (MeshType::Dimension == 1) {
+ auto node_list = m_cell_to_node_matrix[cell_i_id];
+ const CellType cell_type = m_cell_type[cell_i_id];
+ const double Vi = m_Vj[cell_i_id];
+
+ if (cell_type == CellType::Line) {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+
+ const LineTransformation<1> T{x0, x1};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+ } else {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+ } else if constexpr (MeshType::Dimension == 2) {
+ auto node_list = m_cell_to_node_matrix[cell_i_id];
+ const CellType cell_type = m_cell_type[cell_i_id];
+ const double Vi = m_Vj[cell_i_id];
+
+ switch (cell_type) {
+ case CellType::Triangle: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+
+ const TriangleTransformation<2> T{x0, x1, x2};
+ const auto& quadrature =
+ QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{m_polynomial_degree});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+ break;
+ }
+ case CellType::Quadrangle: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[3]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+
+ const SquareTransformation<2> T{x0, x1, x2, x3};
+ const auto& quadrature =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree + 1});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+ break;
+ }
+ default: {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+ }
+ } else {
+ static_assert(MeshType::Dimension == 3);
+ auto node_list = m_cell_to_node_matrix[cell_i_id];
+ const CellType cell_type = m_cell_type[cell_i_id];
+ const double Vi = m_Vj[cell_i_id];
+ switch (cell_type) {
+ case CellType::Tetrahedron: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, m_xr[node_list[3]]);
+
+ const TetrahedronTransformation T{x0, x1, x2, x3};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{m_polynomial_degree});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+ break;
+ }
+ case CellType::Prism: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, m_xr[node_list[4]]);
+ const auto x4 = symmetrize_coordinates(origin, normal, m_xr[node_list[3]]);
+ const auto x5 = symmetrize_coordinates(origin, normal, m_xr[node_list[5]]);
+
+ const PrismTransformation T{x0, x1, x2, //
+ x3, x4, x5};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{m_polynomial_degree + 1});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+ break;
+ }
+ case CellType::Pyramid: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[3]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+ const auto x4 = symmetrize_coordinates(origin, normal, m_xr[node_list[4]]);
+ const PyramidTransformation T{x0, x1, x2, x3, x4};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{m_polynomial_degree + 1});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+ break;
+ }
+ case CellType::Hexahedron: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[3]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+ const auto x4 = symmetrize_coordinates(origin, normal, m_xr[node_list[7]]);
+ const auto x5 = symmetrize_coordinates(origin, normal, m_xr[node_list[6]]);
+ const auto x6 = symmetrize_coordinates(origin, normal, m_xr[node_list[5]]);
+ const auto x7 = symmetrize_coordinates(origin, normal, m_xr[node_list[4]]);
+
+ const CubeTransformation T{x0, x1, x2, x3, x4, x5, x6, x7};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree + 1});
+
+ this->_computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+ break;
+ }
+ default: {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+ }
+ }
+}
+
+template <MeshConcept MeshTypeT>
+void
+PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<MeshTypeT>::build(
+ const CellId cell_j_id,
+ ShrinkMatrixView<SmallMatrix<double>>& A)
+{
+ const auto& stencil_cell_list = m_stencil_array[cell_j_id];
+
+ const Rd& Xj = m_xj[cell_j_id];
+
+ this->_computeEjkMean(Xj, cell_j_id, m_mean_j_of_ejk);
+
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = stencil_cell_list[i];
+
+ this->_computeEjkMean(Xj, cell_i_id, m_mean_i_of_ejk);
+
+ for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
+ A(index, l) = m_mean_i_of_ejk[l] - m_mean_j_of_ejk[l];
+ }
+ }
+
+ for (size_t i_symmetry = 0; i_symmetry < m_stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry) {
+ auto& ghost_stencil = m_stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+ const Rd& origin = m_symmetry_origin_list[i_symmetry];
+ const Rd& normal = m_symmetry_normal_list[i_symmetry];
+
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+
+ this->_computeEjkMeanInSymmetricCell(origin, normal, Xj, cell_i_id, m_mean_i_of_ejk);
+
+ for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
+ A(index, l) = m_mean_i_of_ejk[l] - m_mean_j_of_ejk[l];
+ }
+ }
+ }
+}
+
+template <MeshConcept MeshTypeT>
+PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<
+ MeshTypeT>::ElementIntegralReconstructionMatrixBuilder(const MeshType& mesh,
+ const size_t polynomial_degree,
+ const SmallArray<const Rd>& symmetry_origin_list,
+ const SmallArray<const Rd>& symmetry_normal_list,
+ const CellToCellStencilArray& stencil_array)
+ : m_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(
+ polynomial_degree)},
+ m_polynomial_degree{polynomial_degree},
+
+ m_wq_detJ_ek{m_basis_dimension},
+ m_mean_j_of_ejk{m_basis_dimension - 1},
+ m_mean_i_of_ejk{m_basis_dimension - 1},
+
+ m_cell_to_node_matrix{mesh.connectivity().cellToNodeMatrix()},
+ m_stencil_array{stencil_array},
+ m_symmetry_origin_list{symmetry_origin_list},
+ m_symmetry_normal_list{symmetry_normal_list},
+ m_cell_type{mesh.connectivity().cellType()},
+ m_Vj{MeshDataManager::instance().getMeshData(mesh).Vj()},
+ m_xj{MeshDataManager::instance().getMeshData(mesh).xj()},
+ m_xr{mesh.xr()}
+{
+ if constexpr (MeshType::Dimension == 2) {
+ SmallArray<size_t> y_row_index(m_polynomial_degree + 1);
+
+ size_t i_y = 0;
+
+ y_row_index[i_y++] = 0;
+ for (ssize_t n = m_polynomial_degree + 1; n > 1; --n, ++i_y) {
+ y_row_index[i_y] = y_row_index[i_y - 1] + n;
+ }
+
+ m_y_row_index = y_row_index;
+
+ } else if constexpr (MeshType::Dimension == 3) {
+ SmallArray<size_t> yz_row_index((m_polynomial_degree + 2) * (m_polynomial_degree + 1) / 2 + 1);
+ SmallArray<size_t> z_triangle_index(m_polynomial_degree + 1);
+
+ {
+ size_t i_z = 0;
+ size_t i_yz = 0;
+
+ yz_row_index[i_yz++] = 0;
+ for (ssize_t n = m_polynomial_degree + 1; n >= 1; --n) {
+ z_triangle_index[i_z++] = i_yz - 1;
+ for (ssize_t i = n; i >= 1; --i) {
+ yz_row_index[i_yz] = yz_row_index[i_yz - 1] + i;
+ ++i_yz;
+ }
+ }
+ }
+
+ SmallArray<size_t> yz_row_size{yz_row_index.size() - 1};
+ for (size_t i = 0; i < yz_row_size.size(); ++i) {
+ yz_row_size[i] = yz_row_index[i + 1] - yz_row_index[i];
+ }
+
+ m_yz_row_index = yz_row_index;
+ m_z_triangle_index = z_triangle_index;
+ m_yz_row_size = yz_row_size;
+ }
+}
+
+template void PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<Mesh<1>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template void PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<Mesh<2>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template void PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<Mesh<3>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<
+ Mesh<1>>::ElementIntegralReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
+
+template PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<
+ Mesh<2>>::ElementIntegralReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
+
+template PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<
+ Mesh<3>>::ElementIntegralReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7a3e8f7fc04d6d4ba21d4a36ed02ba24dd0d67d0
--- /dev/null
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
@@ -0,0 +1,88 @@
+#ifndef ELEMENT_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
+#define ELEMENT_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
+
+#include <algebra/ShrinkMatrixView.hpp>
+#include <algebra/SmallMatrix.hpp>
+#include <mesh/CellType.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/StencilArray.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/SmallArray.hpp>
+
+template <size_t Dimension>
+class QuadratureFormula;
+
+template <MeshConcept MeshTypeT>
+class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
+{
+ public:
+ using MeshType = MeshTypeT;
+
+ constexpr static bool handles_high_degrees = true;
+
+ private:
+ using Rd = TinyVector<MeshType::Dimension>;
+
+ const size_t m_basis_dimension;
+ const size_t m_polynomial_degree;
+
+ const SmallArray<double> m_wq_detJ_ek;
+ SmallArray<double> m_mean_j_of_ejk;
+ SmallArray<double> m_mean_i_of_ejk;
+
+ const ItemToItemMatrix<ItemType::cell, ItemType::node> m_cell_to_node_matrix;
+ const CellToCellStencilArray& m_stencil_array;
+
+ const SmallArray<const Rd> m_symmetry_origin_list;
+ const SmallArray<const Rd> m_symmetry_normal_list;
+
+ const CellValue<const CellType> m_cell_type;
+ const CellValue<const double> m_Vj;
+ const CellValue<const Rd> m_xj;
+ const NodeValue<const Rd> m_xr;
+
+ // 2D
+ SmallArray<const size_t> m_y_row_index;
+
+ // 3D
+ SmallArray<const size_t> m_yz_row_index;
+ SmallArray<const size_t> m_z_triangle_index;
+ SmallArray<const size_t> m_yz_row_size;
+
+ template <typename ConformTransformationT>
+ void _computeEjkMean(const QuadratureFormula<MeshType::Dimension>& quadrature,
+ const ConformTransformationT& T,
+ const Rd& Xj,
+ const double Vi,
+ SmallArray<double>& mean_of_ejk) noexcept(NO_ASSERT);
+
+ void _computeEjkMean(const TinyVector<MeshType::Dimension>& Xj,
+ const CellId& cell_i_id,
+ SmallArray<double>& mean_of_ejk);
+
+ void _computeEjkMeanInSymmetricCell(const Rd& origin,
+ const Rd& normal,
+ const Rd& Xj,
+ const CellId& cell_i_id,
+ SmallArray<double>& mean_of_ejk);
+
+ public:
+ PUGS_INLINE
+ SmallArray<const double>
+ meanjOfEjk() const
+ {
+ return m_mean_j_of_ejk;
+ }
+
+ void build(const CellId cell_j_id, ShrinkMatrixView<SmallMatrix<double>>& A);
+
+ ElementIntegralReconstructionMatrixBuilder(const MeshType& mesh,
+ const size_t polynomial_degree,
+ const SmallArray<const Rd>& symmetry_origin_list,
+ const SmallArray<const Rd>& symmetry_normal_list,
+ const CellToCellStencilArray& stencil_array);
+
+ ~ElementIntegralReconstructionMatrixBuilder() = default;
+};
+
+#endif // ELEMENT_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
diff --git a/src/scheme/reconstruction_utils/MutableDiscreteFunctionDPkVariant.hpp b/src/scheme/reconstruction_utils/MutableDiscreteFunctionDPkVariant.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..761f3c1129b00609b1be9a23b5d256db4b48c084
--- /dev/null
+++ b/src/scheme/reconstruction_utils/MutableDiscreteFunctionDPkVariant.hpp
@@ -0,0 +1,130 @@
+#ifndef MUTABLE_DISCRETE_FUNCTION_D_PK_VARIANT_HPP
+#define MUTABLE_DISCRETE_FUNCTION_D_PK_VARIANT_HPP
+
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+
+#include <variant>
+
+class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
+{
+ public:
+ using Variant = std::variant<DiscreteFunctionDPk<1, double>,
+ DiscreteFunctionDPk<1, TinyVector<1>>,
+ DiscreteFunctionDPk<1, TinyVector<2>>,
+ DiscreteFunctionDPk<1, TinyVector<3>>,
+ DiscreteFunctionDPk<1, TinyMatrix<1>>,
+ DiscreteFunctionDPk<1, TinyMatrix<2>>,
+ DiscreteFunctionDPk<1, TinyMatrix<3>>,
+
+ DiscreteFunctionDPk<2, double>,
+ DiscreteFunctionDPk<2, TinyVector<1>>,
+ DiscreteFunctionDPk<2, TinyVector<2>>,
+ DiscreteFunctionDPk<2, TinyVector<3>>,
+ DiscreteFunctionDPk<2, TinyMatrix<1>>,
+ DiscreteFunctionDPk<2, TinyMatrix<2>>,
+ DiscreteFunctionDPk<2, TinyMatrix<3>>,
+
+ DiscreteFunctionDPk<3, double>,
+ DiscreteFunctionDPk<3, TinyVector<1>>,
+ DiscreteFunctionDPk<3, TinyVector<2>>,
+ DiscreteFunctionDPk<3, TinyVector<3>>,
+ DiscreteFunctionDPk<3, TinyMatrix<1>>,
+ DiscreteFunctionDPk<3, TinyMatrix<2>>,
+ DiscreteFunctionDPk<3, TinyMatrix<3>>,
+
+ DiscreteFunctionDPkVector<1, double>,
+ DiscreteFunctionDPkVector<1, TinyVector<1>>,
+ DiscreteFunctionDPkVector<1, TinyVector<2>>,
+ DiscreteFunctionDPkVector<1, TinyVector<3>>,
+ DiscreteFunctionDPkVector<1, TinyMatrix<1>>,
+ DiscreteFunctionDPkVector<1, TinyMatrix<2>>,
+ DiscreteFunctionDPkVector<1, TinyMatrix<3>>,
+
+ DiscreteFunctionDPkVector<2, double>,
+ DiscreteFunctionDPkVector<2, TinyVector<1>>,
+ DiscreteFunctionDPkVector<2, TinyVector<2>>,
+ DiscreteFunctionDPkVector<2, TinyVector<3>>,
+ DiscreteFunctionDPkVector<2, TinyMatrix<1>>,
+ DiscreteFunctionDPkVector<2, TinyMatrix<2>>,
+ DiscreteFunctionDPkVector<2, TinyMatrix<3>>,
+
+ DiscreteFunctionDPkVector<3, double>,
+ DiscreteFunctionDPkVector<3, TinyVector<1>>,
+ DiscreteFunctionDPkVector<3, TinyVector<2>>,
+ DiscreteFunctionDPkVector<3, TinyVector<3>>,
+ DiscreteFunctionDPkVector<3, TinyMatrix<1>>,
+ DiscreteFunctionDPkVector<3, TinyMatrix<2>>,
+ DiscreteFunctionDPkVector<3, TinyMatrix<3>>>;
+
+ private:
+ Variant m_mutable_discrete_function_dpk;
+
+ public:
+ PUGS_INLINE
+ const Variant&
+ mutableDiscreteFunctionDPk() const
+ {
+ return m_mutable_discrete_function_dpk;
+ }
+
+ template <typename DiscreteFunctionDPkT>
+ PUGS_INLINE auto&&
+ get() const
+ {
+ static_assert(is_discrete_function_dpk_v<DiscreteFunctionDPkT>, "invalid template argument");
+#ifndef NDEBUG
+ if (not std::holds_alternative<DiscreteFunctionDPkT>(this->m_mutable_discrete_function_dpk)) {
+ std::ostringstream error_msg;
+ error_msg << "invalid discrete function type\n";
+ error_msg << "- required " << rang::fgB::red << demangle<DiscreteFunctionDPkT>() << rang::fg::reset << '\n';
+ error_msg << "- contains " << rang::fgB::yellow
+ << std::visit([](auto&& f) -> std::string { return demangle<decltype(f)>(); },
+ this->m_mutable_discrete_function_dpk)
+ << rang::fg::reset;
+ throw NormalError(error_msg.str());
+ }
+#endif // NDEBUG
+
+ return std::get<DiscreteFunctionDPkT>(this->mutableDiscreteFunctionDPk());
+ }
+
+ template <size_t Dimension, typename DataType>
+ MutableDiscreteFunctionDPkVariant(const DiscreteFunctionDPk<Dimension, DataType>& discrete_function_dpk)
+ : m_mutable_discrete_function_dpk{discrete_function_dpk}
+ {
+ static_assert(std::is_same_v<DataType, double> or //
+ std::is_same_v<DataType, TinyVector<1, double>> or //
+ std::is_same_v<DataType, TinyVector<2, double>> or //
+ std::is_same_v<DataType, TinyVector<3, double>> or //
+ std::is_same_v<DataType, TinyMatrix<1, 1, double>> or //
+ std::is_same_v<DataType, TinyMatrix<2, 2, double>> or //
+ std::is_same_v<DataType, TinyMatrix<3, 3, double>>,
+ "DiscreteFunctionDPk with this DataType is not allowed in variant");
+ }
+
+ template <size_t Dimension, typename DataType>
+ MutableDiscreteFunctionDPkVariant(const DiscreteFunctionDPkVector<Dimension, DataType>& discrete_function_dpk)
+ : m_mutable_discrete_function_dpk{discrete_function_dpk}
+ {
+ static_assert(std::is_same_v<DataType, double> or //
+ std::is_same_v<DataType, TinyVector<1, double>> or //
+ std::is_same_v<DataType, TinyVector<2, double>> or //
+ std::is_same_v<DataType, TinyVector<3, double>> or //
+ std::is_same_v<DataType, TinyMatrix<1, 1, double>> or //
+ std::is_same_v<DataType, TinyMatrix<2, 2, double>> or //
+ std::is_same_v<DataType, TinyMatrix<3, 3, double>>,
+ "DiscreteFunctionDPkVector with this DataType is not allowed in variant");
+ }
+
+ MutableDiscreteFunctionDPkVariant& operator=(MutableDiscreteFunctionDPkVariant&&) = default;
+ MutableDiscreteFunctionDPkVariant& operator=(const MutableDiscreteFunctionDPkVariant&) = default;
+
+ MutableDiscreteFunctionDPkVariant(const MutableDiscreteFunctionDPkVariant&) = default;
+ MutableDiscreteFunctionDPkVariant(MutableDiscreteFunctionDPkVariant&&) = default;
+
+ MutableDiscreteFunctionDPkVariant() = delete;
+ ~MutableDiscreteFunctionDPkVariant() = default;
+};
+
+#endif // MUTABLE_DISCRETE_FUNCTION_D_PK_VARIANT_HPP
diff --git a/src/utils/SignalManager.cpp b/src/utils/SignalManager.cpp
index 4b2046a98796fc72cef406f77e7498930089f533..5054622c0d9959c089ddfa3f3a323820614c5806 100644
--- a/src/utils/SignalManager.cpp
+++ b/src/utils/SignalManager.cpp
@@ -52,17 +52,27 @@ SignalManager::signalName(int signal)
void
SignalManager::pauseForDebug(int signal)
{
- if (std::string(PUGS_BUILD_TYPE) != "Release") {
- if (s_pause_on_error) {
- // Each failing process must write
- std::cerr.clear();
+ if (s_pause_on_error) {
+ // Each failing process must write
+ std::cerr.clear();
+ if (std::string(PUGS_BUILD_TYPE) != "Release") {
+ char hostname[HOST_NAME_MAX + 1];
+ gethostname(hostname, HOST_NAME_MAX + 1);
std::cerr << "\n======================================\n"
- << rang::style::reset << rang::fg::reset << rang::style::bold << "to attach gdb to this process run\n"
+ << rang::style::reset << rang::fg::reset << rang::style::bold << "Process paused on host \""
+ << rang::fg::yellow << hostname << rang::fg::reset << "\"\n"
+ << "to attach gdb to this process run\n"
<< "\tgdb -pid " << rang::fg::red << getpid() << rang::fg::reset << '\n'
<< "else press Control-C to exit\n"
<< rang::style::reset << "======================================\n"
<< std::flush;
pause();
+ } else {
+ std::cerr << '\n'
+ << rang::style::bold
+ << "Pausing is useless for Release version.\n"
+ "To attach debugger use Debug built type."
+ << rang::style::reset << '\n';
}
}
std::exit(signal);
@@ -73,56 +83,52 @@ SignalManager::handler(int signal)
{
static std::mutex mutex;
- if (mutex.try_lock()) {
- std::signal(SIGTERM, SIG_DFL);
- std::signal(SIGINT, SIG_DFL);
- std::signal(SIGABRT, SIG_DFL);
-
- // Each failing process must write
- std::cerr.clear();
+ std::lock_guard<std::mutex> lock(mutex);
- std::cerr << BacktraceManager{} << '\n';
+ std::signal(SIGINT, SIG_BLOCK);
- std::cerr << "\n *** " << rang::style::reset << rang::fg::reset << rang::style::bold << "Signal " << rang::fgB::red
- << signalName(signal) << rang::fg::reset << " caught" << rang::style::reset << " ***\n\n";
+ // Each failing process must write
+ std::cerr.clear();
- std::exception_ptr eptr = std::current_exception();
- try {
- if (eptr) {
- std::rethrow_exception(eptr);
- } else {
- std::ostringstream error_msg;
- error_msg << "received " << signalName(signal);
- std::cerr << ASTExecutionStack::getInstance().errorMessageAt(error_msg.str()) << '\n';
- }
- }
- catch (const IBacktraceError& backtrace_error) {
- auto source_location = backtrace_error.sourceLocation();
- std::cerr << rang::fgB::cyan << source_location.file_name() << ':' << source_location.line() << ':'
- << source_location.column() << ':' << rang::fg::reset << rang::fgB::yellow
- << " threw the following exception" << rang::fg::reset << "\n\n";
- std::cerr << ASTExecutionStack::getInstance().errorMessageAt(backtrace_error.what()) << '\n';
- }
- catch (const ParseError& parse_error) {
- auto p = parse_error.positions().front();
- std::cerr << rang::style::bold << p.source << ':' << p.line << ':' << p.column << ':' << rang::style::reset
- << rang::fgB::red << " error: " << rang::fg::reset << rang::style::bold << parse_error.what()
- << rang::style::reset << '\n';
- }
- catch (const IExitError& exit_error) {
- std::cerr << ASTExecutionStack::getInstance().errorMessageAt(exit_error.what()) << '\n';
- }
- catch (const AssertError& assert_error) {
- std::cerr << assert_error << '\n';
- }
- catch (...) {
- std::cerr << "Unknown exception!\n";
- }
+ std::cerr << BacktraceManager{} << '\n';
- SignalManager::pauseForDebug(signal);
+ std::cerr << "\n *** " << rang::style::reset << rang::fg::reset << rang::style::bold << "Signal " << rang::fgB::red
+ << signalName(signal) << rang::fg::reset << " caught" << rang::style::reset << " ***\n\n";
- mutex.unlock();
+ std::exception_ptr eptr = std::current_exception();
+ try {
+ if (eptr) {
+ std::rethrow_exception(eptr);
+ } else {
+ std::ostringstream error_msg;
+ error_msg << "received " << signalName(signal);
+ std::cerr << ASTExecutionStack::getInstance().errorMessageAt(error_msg.str()) << '\n';
+ }
+ }
+ catch (const IBacktraceError& backtrace_error) {
+ auto source_location = backtrace_error.sourceLocation();
+ std::cerr << rang::fgB::cyan << source_location.file_name() << ':' << source_location.line() << ':'
+ << source_location.column() << ':' << rang::fg::reset << rang::fgB::yellow
+ << " threw the following exception" << rang::fg::reset << "\n\n";
+ std::cerr << ASTExecutionStack::getInstance().errorMessageAt(backtrace_error.what()) << '\n';
}
+ catch (const ParseError& parse_error) {
+ auto p = parse_error.positions().front();
+ std::cerr << rang::style::bold << p.source << ':' << p.line << ':' << p.column << ':' << rang::style::reset
+ << rang::fgB::red << " error: " << rang::fg::reset << rang::style::bold << parse_error.what()
+ << rang::style::reset << '\n';
+ }
+ catch (const IExitError& exit_error) {
+ std::cerr << ASTExecutionStack::getInstance().errorMessageAt(exit_error.what()) << '\n';
+ }
+ catch (const AssertError& assert_error) {
+ std::cerr << assert_error << '\n';
+ }
+ catch (...) {
+ std::cerr << "Unknown exception!\n";
+ }
+
+ SignalManager::pauseForDebug(signal);
}
void
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 52d3a91ecb06c63b44bb6b391eb1c51bde95ed83..8771bb4cc0bae2587948acbe5eb6e49c737690b4 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -285,6 +285,7 @@ add_executable (mpi_unit_tests
test_Partitioner.cpp
test_PolynomialReconstruction_degree_1.cpp
test_PolynomialReconstruction_degree_2.cpp
+ test_PolynomialReconstruction_degree_3.cpp
test_PolynomialReconstructionDescriptor.cpp
test_RandomEngine.cpp
test_StencilBuilder_cell2cell.cpp
@@ -320,6 +321,7 @@ target_link_libraries (unit_tests
PugsAlgebra
PugsAnalysis
PugsScheme
+ PugsSchemeReconstructionUtils
PugsOutput
PugsUtils
PugsCheckpointing
@@ -350,6 +352,7 @@ target_link_libraries (mpi_unit_tests
PugsUtils
PugsLanguageUtils
PugsScheme
+ PugsSchemeReconstructionUtils
PugsOutput
PugsUtils
PugsCheckpointing
diff --git a/tests/test_PolynomialReconstruction_degree_2.cpp b/tests/test_PolynomialReconstruction_degree_2.cpp
index 435084af00385cda388588fb21860fe32680c73e..db13612c1d8c5efad94a6f890251e366826d014b 100644
--- a/tests/test_PolynomialReconstruction_degree_2.cpp
+++ b/tests/test_PolynomialReconstruction_degree_2.cpp
@@ -118,7 +118,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -139,7 +139,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -172,7 +172,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
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));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -219,7 +219,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
{
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));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
{
@@ -231,7 +231,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
{
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));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -477,201 +477,567 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
SECTION("with symmetries")
{
-#warning continue here
- // SECTION("1D")
- // {
- // std::vector<PolynomialReconstructionDescriptor> descriptor_list =
- // {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::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::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));
- // }
- // }
- // }
- // }
+ SECTION("1D")
+ {
+ std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ {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>;
+
+ 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; };
+ auto p2 = [](const R1& x) { return +1.4 * (x[0] + 1) * (x[0] + 1) - 0.6; };
+
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
+ {
+ SECTION("R data")
+ {
+ auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ auto& mesh = *p_mesh;
+
+ auto R_exact = p0;
+
+ 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-12));
+ }
+
+ SECTION("R1x1 data")
+ {
+ using R1x1 = TinyMatrix<1>;
+
+ auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ auto& mesh = *p_mesh;
+
+ auto R1x1_exact = [&](const R1& x) { return R1x1{p0(x)}; };
+
+ DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree, std::function(R1x1_exact));
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
+
+ auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1x1>>();
+
+ 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;
+
+ std::array<std::function<double(const R1&)>, 3> vector_exact //
+ = {p0, p1, p2};
+
+ 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-12));
+ }
+
+ SECTION("R1x1 vector data")
+ {
+ using R1x1 = TinyMatrix<1>;
+
+ auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
+
+ std::array<std::function<R1x1(const R1&)>, 3> vector_exact //
+ = {[&](const R1& x) { return R1x1{p2(x)}; }, //
+ [&](const R1& x) { return R1x1{p0(x)}; }, //
+ [&](const R1& x) { return R1x1{p1(x)}; }};
+
+ 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 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));
+ }
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ 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")}},
+ 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>;
+
+ auto p_initial_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+ 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)};
+
+ auto xr = initial_mesh.xr();
+
+ NodeValue<R2> 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<2>>(initial_mesh.shared_connectivity(), new_xr);
+ const Mesh<2>& mesh = *p_mesh;
+ // inverse rotation
+ TinyMatrix<2> inv_T{std::cos(theta), std::sin(theta), //
+ -std::sin(theta), std::cos(theta)};
+
+ auto p0 = [&inv_T](const R2& X) {
+ const R2 Y = inv_T * X;
+ const double x = Y[0] - 2;
+ const double y = Y[1] - 1;
+ return +1.7 * x * x + 2 * y * y - 1.1;
+ };
+
+ auto p1 = [&inv_T](const R2& X) {
+ const R2 Y = inv_T * X;
+ const double x = Y[0] - 2;
+ const double y = Y[1] - 1;
+ return -1.3 * x * x - 0.2 * y * y + 0.7;
+ };
+
+ auto p2 = [&inv_T](const R2& X) {
+ const R2 Y = inv_T * X;
+ const double x = Y[0] - 2;
+ const double y = Y[1] - 1;
+ return +2.6 * x * x - 1.4 * y * y - 1.9;
+ };
+
+ auto p3 = [&inv_T](const R2& X) {
+ const R2 Y = inv_T * X;
+ const double x = Y[0] - 2;
+ const double y = Y[1] - 1;
+ return -0.6 * x * x + 1.4 * y * y + 2.3;
+ };
+
+ auto q0 = [&inv_T](const R2& X) {
+ const R2 Y = inv_T * X;
+ const double x = Y[0] - 2;
+ const double y = Y[1] - 1;
+ return 3 * x * y;
+ };
+
+ auto q1 = [&inv_T](const R2& X) {
+ const R2 Y = inv_T * X;
+ const double x = Y[0] - 2;
+ const double y = Y[1] - 1;
+ return -1.3 * x * y;
+ };
+
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
+ {
+ 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<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));
+ }
+ }
+ }
+ }
+ }
+ }
}
}
diff --git a/tests/test_PolynomialReconstruction_degree_3.cpp b/tests/test_PolynomialReconstruction_degree_3.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e20fa89a3bdcba29f61a771deacba273c536f28b
--- /dev/null
+++ b/tests/test_PolynomialReconstruction_degree_3.cpp
@@ -0,0 +1,766 @@
+#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 <mesh/CartesianMeshBuilder.hpp>
+
+#include <DiscreteFunctionDPkForTests.hpp>
+#include <MeshDataBaseForTests.hpp>
+
+#include <NbGhostLayersTester.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PolynomialReconstruction_degree_3", "[scheme]")
+{
+ constexpr size_t degree = 3;
+
+ constexpr size_t nb_ghost_layers = 3;
+ NbGhostLayersTester t{nb_ghost_layers};
+
+ SECTION("without symmetries")
+ {
+ std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ {PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree},
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree}};
+
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
+ {
+ SECTION("1D")
+ {
+ using R1 = TinyVector<1>;
+
+ auto p0 = [](const R1& X) {
+ const double x = X[0];
+ return +2.3 + (1.4 + (1.7 - 2.3 * x) * x) * x;
+ };
+ auto p1 = [](const R1& X) {
+ const double x = X[0];
+ return -1.2 - (2.3 - (1.1 + 2.1 * x) * x) * x;
+ };
+ auto p2 = [](const R1& X) {
+ const double x = X[0];
+ return +2.1 + (2.1 + (3.1 - 1.7 * x) * x) * x;
+ };
+ auto p3 = [](const R1& X) {
+ const double x = X[0];
+ return -1.7 + (1.4 + (1.6 - 3.1 * x) * x) * x;
+ };
+ auto p4 = [](const R1& X) {
+ const double x = X[0];
+ return +1.1 - (1.3 - (2.1 + 1.5 * x) * x) * x;
+ };
+ auto p5 = [](const R1& X) {
+ const double x = X[0];
+ return +1.9 - (2.1 + (1.6 - 2.7 * x) * x) * x;
+ };
+ auto p6 = [](const R1& X) {
+ const double x = X[0];
+ return -0.7 + (1.4 + (2.1 + 1.1 * x) * x) * x;
+ };
+ auto p7 = [](const R1& X) {
+ const double x = X[0];
+ return -1.4 - (1.2 + (1.5 - 2.1 * x) * x) * x;
+ };
+ auto p8 = [](const R1& X) {
+ const double x = X[0];
+ return -2.1 + (1.1 - (1.7 + 1.2 * x) * x) * x;
+ };
+ auto p9 = [](const R1& X) {
+ const double x = X[0];
+ return +1.8 - (3.1 + (2.1 - 2.4 * x) * x) * x;
+ };
+
+ 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 = p0;
+
+ 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-12));
+ }
+ }
+ }
+
+ 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{p2(x), p4(x), p1(x)}; };
+
+ 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-12));
+ }
+ }
+ }
+
+ 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{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));
+
+ 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-12));
+ }
+ }
+ }
+
+ 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 = {p1, p7, p9};
+
+ 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-12));
+ }
+ }
+ }
+
+ 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&)>, 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);
+
+ 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-12));
+ }
+ }
+ }
+
+ 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 = p0;
+
+ auto R3_exact = [&](const R1& x) -> R3 { return R3{p9(x), p4(x), p7(x)}; };
+
+ 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 = {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));
+ 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-12));
+ }
+
+ {
+ 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-12));
+ }
+
+ {
+ 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-12));
+ }
+
+ {
+ 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-12));
+ }
+ }
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ using R2 = TinyVector<2>;
+ auto p0 = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double y2 = y * y;
+ const double x3 = x2 * x;
+ const double x2y = x2 * y;
+ const double xy2 = x * y2;
+ const double y3 = y * y2;
+
+ return +2.3 //
+ + 1.7 * x - 1.3 * y //
+ + 1.2 * x2 + 1.3 * xy - 3.2 * y2 //
+ - 1.3 * x3 + 2.1 * x2y - 1.6 * xy2 + 2.1 * y3;
+ };
+
+ auto p1 = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double y2 = y * y;
+ const double x3 = x2 * x;
+ const double x2y = x2 * y;
+ const double xy2 = x * y2;
+ const double y3 = y * y2;
+
+ return +1.4 //
+ + 1.6 * x - 2.1 * y //
+ + 1.3 * x2 + 2.6 * xy - 1.4 * y2 //
+ - 1.2 * x3 - 1.7 * x2y + 2.1 * xy2 - 2.2 * y3;
+ };
+
+ auto p2 = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double y2 = y * y;
+ const double x3 = x2 * x;
+ const double x2y = x2 * y;
+ const double xy2 = x * y2;
+ const double y3 = y * y2;
+
+ return -1.2 //
+ + 2.3 * x - 1.6 * y //
+ - 1.2 * x2 + 2.4 * xy - 1.9 * y2 //
+ + 0.9 * x3 + 1.5 * x2y - 2.3 * xy2 - 1.6 * y3;
+ };
+
+ auto p3 = [](const R2& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double x2 = x * x;
+ const double xy = x * y;
+ const double y2 = y * y;
+ const double x3 = x2 * x;
+ const double x2y = x2 * y;
+ const double xy2 = x * y2;
+ const double y3 = y * y2;
+
+ return +2.4 //
+ + 2.5 * x + 1.4 * y //
+ - 2.7 * x2 + 1.9 * xy - 2.2 * y2 //
+ - 1.3 * x3 + 2.3 * x2y - 1.4 * xy2 + 2.2 * y3;
+ };
+
+ 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 = p0;
+
+ DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree + 1, 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-12));
+ }
+ }
+ }
+
+ 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{p1(x), p2(x), p3(x)}; };
+
+ DiscreteFunctionP0 uh = test_only::exact_projection(mesh, degree + 1, 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-12));
+ }
+ }
+ }
+
+ 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{p0(x), p1(x), //
+ p2(x), p3(x)};
+ };
+
+ DiscreteFunctionP0 Ah = test_only::exact_projection(mesh, degree + 1, 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-12));
+ }
+ }
+ }
+
+ 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 = {p0, p1, p2, p3};
+
+ DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree + 1, 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("3D")
+ {
+ using R3 = TinyVector<3>;
+
+ auto p = [](const R3& X, const std::array<double, 20>& a) -> double {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ const double xy = x * y;
+ const double xz = x * z;
+ const double yz = y * z;
+ const double x2 = x * x;
+ const double y2 = y * y;
+ const double z2 = z * z;
+ const double xyz = x * y * z;
+ const double x2y = x2 * y;
+ const double x2z = x2 * z;
+ const double xy2 = x * y2;
+ const double y2z = y2 * z;
+ const double xz2 = x * z2;
+ const double yz2 = y * z2;
+ const double x3 = x2 * x;
+ const double y3 = y2 * y;
+ const double z3 = z2 * z;
+
+ return a[0] + a[1] * x + a[2] * y + a[3] * z //
+ + a[4] * x2 + a[5] * y2 + a[6] * z2 //
+ + a[7] * xy + a[8] * xz + a[9] * yz //
+ + a[10] * x3 + a[11] * y3 + a[12] * z3 //
+ + a[13] * x2y + a[14] * x2z + a[15] * xy2 //
+ + a[16] * y2z + a[17] * xz2 + a[18] * yz2 //
+ + a[19] * xyz //
+ ;
+ };
+
+ constexpr std::array<double, 20> a0 = {+2.3, +1.7, -1.3, +2.1, +1.7, +1.4, +1.7, -2.3, +1.6, -1.9,
+ +1.2, -2.1, -1.1, -1.7, -1.3, +0.9, -0.7, +1.5, -0.7, +2.8};
+
+ auto p0 = [&p, &a0](const R3& X) -> double { return p(X, a0); };
+
+ constexpr std::array<double, 20> a1 = {-1.3, +2.2, +0.1, -2.5, +0.2, -2.3, -1.4, +0.9, +0.2, -0.3,
+ +2.4, -1.2, +1.7, -2.2, +0.6, +1.9, +1.0, -0.8, +2.4, +2.4};
+
+ auto p1 = [&p, &a1](const R3& X) -> double { return p(X, a1); };
+
+ constexpr std::array<double, 20> a2 = {+1.9, -1.2, -0.4, -1.2, -0.8, +1.4, +0.5, -1.6, +1.1, -0.7,
+ +0.6, +2.3, -1.8, -1.9, -0.3, -2.4, -1.7, +0.2, -2.4, +1.9};
+
+ auto p2 = [&p, &a2](const R3& X) -> double { return p(X, a2); };
+
+ constexpr std::array<double, 20> a3 = {+0.8, +0.5, +1.3, -2.3, +0.9, -0.4, -2.0, +1.8, +0.5, +0.7,
+ +1.0, -0.4, +1.1, +1.8, -0.4, +1.1, -0.0, +1.4, +1.9, -2.2};
+
+ auto p3 = [&p, &a3](const R3& X) -> double { return p(X, a3); };
+
+ auto p_mesh = CartesianMeshBuilder{TinyVector<3>{-0.5, -0.5, -0.5}, TinyVector<3>{3.5, 3.5, 3.5},
+ TinyVector<3, size_t>{4, 4, 4}}
+ .mesh()
+ ->get<Mesh<3>>();
+ const auto& mesh = *p_mesh;
+
+ SECTION("R data")
+ {
+ auto R_exact = p0;
+
+ DiscreteFunctionP0 fh = test_only::exact_projection(mesh, degree + 3, 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")
+ {
+ 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));
+
+ 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>;
+
+ 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));
+
+ 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")
+ {
+ std::array<std::function<double(const R3&)>, 4> vector_exact = {p0, p1, p2, p3};
+
+ 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("with symmetries")
+ {
+ SECTION("1D")
+ {
+ std::vector<PolynomialReconstructionDescriptor> descriptor_list =
+ {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>;
+
+ auto p0 = [](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1) * (x[0] + 1) * (x[0] + 1)}; };
+ auto p1 = [](const R1& x) -> R1 { return R1{-1.2 * (x[0] + 1) * (x[0] + 1) * (x[0] + 1)}; };
+ auto p2 = [](const R1& x) -> R1 { return R1{+1.4 * (x[0] + 1) * (x[0] + 1) * (x[0] + 1)}; };
+
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
+ {
+ SECTION("R1 data")
+ {
+ auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ auto& mesh = *p_mesh;
+
+ auto R1_exact = p0;
+
+ 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")
+ {
+ auto p_mesh = MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+ auto& mesh = *p_mesh;
+
+ std::array<std::function<R1(const R1&)>, 3> vector_exact //
+ = {p0, p1, p2};
+
+ 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::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>;
+
+ auto p_mesh = CartesianMeshBuilder{TinyVector<2>{0, 0}, TinyVector<2>{2, 1}, TinyVector<2, size_t>{3, 3}}
+ .mesh()
+ ->get<Mesh<2>>();
+
+ auto& mesh = *p_mesh;
+
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
+ {
+ SECTION("R^2 data")
+ {
+ auto R2_exact = [](const R2& x) -> R2 {
+ return R2{+2.3 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2), //
+ -1.3 * (x[1] - 1) * (x[1] - 1) * (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")
+ {
+ std::array<std::function<R2(const R2&)>, 2> vector_exact //
+ = {[](const R2& x) -> R2 {
+ return R2{+1.7 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2), //
+ -0.6 * (x[1] - 1) * (x[1] - 1) * (x[1] - 1)};
+ },
+ [](const R2& x) -> R2 {
+ return R2{-2.3 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2), //
+ +1.1 * (x[1] - 1) * (x[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::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>;
+
+ auto p_mesh = CartesianMeshBuilder{TinyVector<3>{0, 0, 0}, TinyVector<3>{2, 1, 1}, TinyVector<3, size_t>{3, 3, 3}}
+ .mesh()
+ ->get<Mesh<3>>();
+
+ auto& mesh = *p_mesh;
+
+ for (auto descriptor : descriptor_list) {
+ SECTION(name(descriptor.integrationMethodType()))
+ {
+ SECTION("R^3 data")
+ {
+ auto R3_exact = [](const R3& x) -> R3 {
+ return R3{+2.3 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2), //
+ -1.3 * (x[1] - 1) * (x[1] - 1) * (x[1] - 1), //
+ +1.4 * (x[2] - 1) * (x[2] - 1) * (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")
+ {
+ std::array<std::function<R3(const R3&)>, 2> vector_exact //
+ = {[](const R3& x) -> R3 {
+ return R3{+1.7 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2), //
+ -0.6 * (x[1] - 1) * (x[1] - 1) * (x[1] - 1), //
+ +1.2 * (x[2] - 1) * (x[2] - 1) * (x[2] - 1)};
+ },
+ [](const R3& x) -> R3 {
+ return R3{-2.3 * (x[0] - 2) * (x[0] - 2) * (x[0] - 2), //
+ +1.1 * (x[1] - 1) * (x[1] - 1) * (x[1] - 1), //
+ -0.3 * (x[2] - 1) * (x[2] - 1) * (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));
+ }
+ }
+ }
+ }
+ }
+}