From 6a5fcd4b13b81249de3d1f752908a732dc8d3772 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Wed, 30 Apr 2025 18:50:56 +0200
Subject: [PATCH] Begin code clean-up
---
src/scheme/PolynomialReconstruction.cpp | 937 +++++++++++++++---------
src/scheme/PolynomialReconstruction.hpp | 9 +
2 files changed, 585 insertions(+), 361 deletions(-)
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 2201e5f0f..661630187 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -49,126 +49,64 @@ symmetrize_coordinates(const TinyVector<Dimension>& origin,
return u - 2 * dot(u - origin, normal) * normal;
}
-class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
+template <MeshConcept MeshType>
+class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder
{
- 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>>,
+ private:
+ using Rd = TinyVector<MeshType::Dimension>;
- 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>>>;
+ const size_t m_basis_dimension;
+ const SmallArray<const Rd> m_symmetry_origin_list;
+ const SmallArray<const Rd> m_symmetry_normal_list;
- private:
- Variant m_mutable_discrete_function_dpk;
+ const CellToCellStencilArray& m_stencil_array;
+ const CellValue<const Rd> m_xj;
public:
- PUGS_INLINE
- const Variant&
- mutableDiscreteFunctionDPk() const
- {
- return m_mutable_discrete_function_dpk;
- }
-
- template <typename DiscreteFunctionDPkT>
- PUGS_INLINE auto&&
- get() const
+ template <typename MatrixType>
+ void
+ build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A)
{
- 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());
+ 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];
+ }
+ }
}
-#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;
+ CellCenterReconstructionMatrixBuilder(const SmallArray<const Rd>& symmetry_origin_list,
+ const SmallArray<const Rd>& symmetry_normal_list,
+ const CellToCellStencilArray& stencil_array,
+ const CellValue<const Rd>& xj)
+ : m_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(1)},
+ m_symmetry_origin_list{symmetry_origin_list},
+ m_symmetry_normal_list{symmetry_normal_list},
+ m_stencil_array{stencil_array},
+ m_xj{xj}
+ {}
+
+ ~CellCenterReconstructionMatrixBuilder() = default;
};
template <typename ConformTransformationT>
@@ -179,7 +117,7 @@ _computeEjkMean(const QuadratureFormula<1>& quadrature,
const double Vi,
const size_t degree,
const size_t dimension,
- SmallArray<double>& inv_Vj_wq_detJ_ek,
+ const SmallArray<double>& inv_Vj_wq_detJ_ek,
SmallArray<double>& mean_of_ejk)
{
using Rd = TinyVector<1>;
@@ -216,7 +154,7 @@ _computeEjkMean(const QuadratureFormula<2>& quadrature,
const double Vi,
const size_t degree,
const size_t dimension,
- SmallArray<double>& inv_Vj_wq_detJ_ek,
+ const SmallArray<double>& inv_Vj_wq_detJ_ek,
SmallArray<double>& mean_of_ejk)
{
using Rd = TinyVector<2>;
@@ -269,7 +207,7 @@ _computeEjkMean(const QuadratureFormula<3>& quadrature,
const double Vi,
const size_t degree,
const size_t dimension,
- SmallArray<double>& inv_Vj_wq_detJ_ek,
+ const SmallArray<double>& inv_Vj_wq_detJ_ek,
SmallArray<double>& mean_of_ejk)
{
#warning Compute this one for all and rework design
@@ -358,7 +296,7 @@ void _computeEjkMean(const TinyVector<Dimension>& Xj,
const double Vi,
const size_t degree,
const size_t basis_dimension,
- SmallArray<double>& inv_Vj_wq_detJ_ek,
+ const SmallArray<double>& inv_Vj_wq_detJ_ek,
SmallArray<double>& mean_of_ejk);
template <typename NodeListT>
@@ -370,7 +308,7 @@ _computeEjkMean(const TinyVector<1>& Xj,
const double Vi,
const size_t degree,
const size_t basis_dimension,
- SmallArray<double>& inv_Vj_wq_detJ_ek,
+ const SmallArray<double>& inv_Vj_wq_detJ_ek,
SmallArray<double>& mean_of_ejk)
{
if (cell_type == CellType::Line) {
@@ -396,7 +334,7 @@ _computeEjkMeanInSymmetricCell(const TinyVector<1>& origin,
const double Vi,
const size_t degree,
const size_t basis_dimension,
- SmallArray<double>& inv_Vj_wq_detJ_ek,
+ const SmallArray<double>& inv_Vj_wq_detJ_ek,
SmallArray<double>& mean_of_ejk)
{
if (cell_type == CellType::Line) {
@@ -414,6 +352,308 @@ _computeEjkMeanInSymmetricCell(const TinyVector<1>& origin,
}
}
+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,
+ const 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,
+ const 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,
+ const 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,
+ const 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};
+
+ 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 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};
+
+ 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 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]]);
+
+ const CubeTransformation T{x0, x1, x2, x3, x4, x5, x6, x7};
+
+ 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 <MeshConcept MeshType>
+class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
+{
+ private:
+ using Rd = TinyVector<MeshType::Dimension>;
+
+ const size_t m_basis_dimension;
+ const size_t m_polynomial_degree;
+
+ const SmallArray<double> m_inv_Vj_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;
+
+ public:
+ template <typename MatrixType>
+ void
+ build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A)
+ {
+ const auto& stencil_cell_list = m_stencil_array[cell_j_id];
+
+ const Rd& Xj = m_xj[cell_j_id];
+
+ _computeEjkMean(Xj, m_xr, m_cell_to_node_matrix[cell_j_id], m_cell_type[cell_j_id], m_Vj[cell_j_id],
+ m_polynomial_degree, m_basis_dimension, m_inv_Vj_wq_detJ_ek, 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];
+
+ _computeEjkMean(Xj, m_xr, m_cell_to_node_matrix[cell_i_id], m_cell_type[cell_i_id], m_Vj[cell_i_id],
+ m_polynomial_degree, m_basis_dimension, m_inv_Vj_wq_detJ_ek, 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];
+
+ _computeEjkMeanInSymmetricCell(origin, normal, //
+ Xj, m_xr, m_cell_to_node_matrix[cell_i_id], m_cell_type[cell_i_id],
+ m_Vj[cell_i_id], m_polynomial_degree, m_basis_dimension, m_inv_Vj_wq_detJ_ek,
+ 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];
+ }
+ }
+ }
+ }
+
+ ElementIntegralReconstructionMatrixBuilder(const MeshType& mesh,
+ const size_t polynomial_degree,
+ const SmallArray<double>& inv_Vj_wq_detJ_ek,
+ const SmallArray<double>& mean_j_of_ejk,
+ const SmallArray<double>& mean_i_of_ejk,
+ 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_inv_Vj_wq_detJ_ek{inv_Vj_wq_detJ_ek},
+ m_mean_j_of_ejk{mean_j_of_ejk},
+ m_mean_i_of_ejk{mean_i_of_ejk},
+
+ 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()}
+ {}
+
+ ~ElementIntegralReconstructionMatrixBuilder() = default;
+};
+
template <typename ConformTransformationT>
void _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
const ConformTransformationT& T,
@@ -421,7 +661,7 @@ void _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
const double Vi,
const size_t degree,
const size_t dimension,
- SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+ const SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
SmallArray<double>& mean_of_ejk);
template <>
@@ -432,7 +672,7 @@ _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
const double Vi,
const size_t degree,
const size_t dimension,
- SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+ const SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
SmallArray<double>& mean_of_ejk)
{
using Rd = TinyVector<2>;
@@ -481,7 +721,7 @@ _computeEjkMeanByBoundary(const MeshType& mesh,
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,
+ const 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));
@@ -518,7 +758,7 @@ _computeEjkMeanByBoundaryInSymmetricCell(const TinyVector<2>& origin,
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,
+ const 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));
@@ -540,214 +780,234 @@ _computeEjkMeanByBoundaryInSymmetricCell(const TinyVector<2>& origin,
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);
- }
- }
-}
-
-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)});
- }
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ }
}
}
-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)
+template <MeshConcept MeshType>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
{
- if (cell_type == CellType::Tetrahedron) {
- const TetrahedronTransformation T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]], xr[node_list[3]]};
+ private:
+ using Rd = TinyVector<MeshType::Dimension>;
- const auto& quadrature = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{degree});
+ const MeshType& m_mesh;
+ const size_t m_basis_dimension;
+ const size_t m_polynomial_degree;
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ const SmallArray<double> m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek;
+ SmallArray<double> m_mean_j_of_ejk;
+ SmallArray<double> m_mean_i_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 ItemToItemMatrix<ItemType::cell, ItemType::face> m_cell_to_face_matrix;
+ const CellToCellStencilArray& m_stencil_array;
- const auto& quadrature = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{degree + 1});
+ const SmallArray<const Rd> m_symmetry_origin_list;
+ const SmallArray<const Rd> m_symmetry_normal_list;
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ const CellValue<const double> m_Vj;
+ const CellValue<const Rd> m_xj;
+ const NodeValue<const Rd> m_xr;
- } 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]]};
+ public:
+ template <typename MatrixType>
+ void
+ build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A)
+ {
+ const auto& stencil_cell_list = m_stencil_array[cell_j_id];
- const auto& quadrature = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{degree + 1});
+ auto face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+ auto cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
+ const Rd& Xj = m_xj[cell_j_id];
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ _computeEjkMeanByBoundary(m_mesh, Xj, cell_j_id, m_cell_to_face_matrix, face_to_node_matrix, cell_face_is_reversed,
+ m_Vj, m_polynomial_degree, m_basis_dimension, m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+ m_mean_j_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]]};
+ 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 auto& quadrature =
- QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{degree + 1});
+ _computeEjkMeanByBoundary(m_mesh, Xj, cell_i_id, m_cell_to_face_matrix, face_to_node_matrix,
+ cell_face_is_reversed, m_Vj, m_polynomial_degree, m_basis_dimension,
+ m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek, m_mean_i_of_ejk);
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_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];
- } else {
- throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
- }
-}
+ const Rd& origin = m_symmetry_origin_list[i_symmetry];
+ const Rd& normal = m_symmetry_normal_list[i_symmetry];
-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]]);
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
- const TetrahedronTransformation T{x0, x1, x2, x3};
+ _computeEjkMeanByBoundaryInSymmetricCell(origin, normal, //
+ Xj, cell_i_id, m_xr, m_cell_to_face_matrix, face_to_node_matrix,
+ cell_face_is_reversed, m_Vj, m_polynomial_degree, m_basis_dimension,
+ m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek, m_mean_i_of_ejk);
- const auto& quadrature = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{degree});
+ 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];
+ }
+ }
+ }
+ }
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh,
+ const size_t polynomial_degree,
+ const SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+ const SmallArray<double>& mean_j_of_ejk,
+ const SmallArray<double>& mean_i_of_ejk,
+ 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_inv_Vj_alpha_p_1_wq_X_prime_orth_ek{inv_Vj_alpha_p_1_wq_X_prime_orth_ek},
+ m_mean_j_of_ejk{mean_j_of_ejk},
+ m_mean_i_of_ejk{mean_i_of_ejk},
+
+ m_cell_to_face_matrix{mesh.connectivity().cellToFaceMatrix()},
+ 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()}
+ {}
+ BoundaryIntegralReconstructionMatrixBuilder() = default;
+ ~BoundaryIntegralReconstructionMatrixBuilder() = default;
+};
- } 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]]);
+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>>,
- const PrismTransformation T{x0, x1, x2, //
- x3, x4, x5};
+ 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>>,
- const auto& quadrature = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{degree + 1});
+ 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>>,
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ 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>>,
- } 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};
+ 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>>,
- const auto& quadrature = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{degree + 1});
+ 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>>>;
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ private:
+ Variant m_mutable_discrete_function_dpk;
- } 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]]);
+ public:
+ PUGS_INLINE
+ const Variant&
+ mutableDiscreteFunctionDPk() const
+ {
+ return m_mutable_discrete_function_dpk;
+ }
- const CubeTransformation T{x0, x1, x2, x3, x4, x5, x6, x7};
+ 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
- const auto& quadrature =
- QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{degree + 1});
+ return std::get<DiscreteFunctionDPkT>(this->mutableDiscreteFunctionDPk());
+ }
- _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ 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");
+ }
- } else {
- throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ 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;
+};
size_t
PolynomialReconstruction::_getNumberOfColumns(
@@ -893,7 +1153,6 @@ PolynomialReconstruction::_build(
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 full_stencil_size = [&](const CellId cell_id) {
@@ -959,7 +1218,6 @@ PolynomialReconstruction::_build(
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);
@@ -974,6 +1232,10 @@ PolynomialReconstruction::_build(
inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[i] = SmallArray<double>(basis_dimension);
}
+ std::shared_ptr p_cell_center_reconstruction_matrix_builder =
+ std::make_shared<CellCenterReconstructionMatrixBuilder<MeshType>>(symmetry_origin_list, symmetry_normal_list,
+ stencil_array, xj);
+
parallel_for(
mesh.numberOfCells(), PUGS_CLASS_LAMBDA(const CellId cell_j_id) {
if (cell_is_owned[cell_j_id]) {
@@ -1215,32 +1477,7 @@ PolynomialReconstruction::_build(
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];
- }
- }
- }
-
+ p_cell_center_reconstruction_matrix_builder->build(cell_j_id, A);
} else if ((m_descriptor.integrationMethodType() == IntegrationMethodType::element) or
(m_descriptor.integrationMethodType() == IntegrationMethodType::boundary)) {
#warning implement boundary based reconstruction for 1d and 3d
@@ -1251,6 +1488,14 @@ PolynomialReconstruction::_build(
SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
SmallArray<double>& mean_i_of_ejk = mean_i_of_ejk_pool[t];
+ BoundaryIntegralReconstructionMatrixBuilder
+ boundary_integral_reconstruction_matrix_builder(*p_mesh, m_descriptor.degree(),
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_j_of_ejk,
+ mean_i_of_ejk, symmetry_origin_list,
+ symmetry_normal_list, stencil_array);
+
+ boundary_integral_reconstruction_matrix_builder.build(cell_j_id, A);
+
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];
@@ -1297,48 +1542,18 @@ PolynomialReconstruction::_build(
throw UnexpectedError("invalid mesh dimension");
}
} else {
+
+#warning incorporate in reconstruction_matrix_builder
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];
+ ElementIntegralReconstructionMatrixBuilder
+ element_integral_reconstruction_matrix_builder(*p_mesh, m_descriptor.degree(), inv_Vj_wq_detJ_ek,
+ mean_j_of_ejk, mean_i_of_ejk, symmetry_origin_list,
+ symmetry_normal_list, stencil_array);
- _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);
-
- for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
- }
- }
- }
+ element_integral_reconstruction_matrix_builder.build(cell_j_id, A);
}
} else {
throw UnexpectedError("invalid integration strategy");
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index bf6013639..349ae4577 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -12,6 +12,15 @@ class PolynomialReconstruction
private:
class MutableDiscreteFunctionDPkVariant;
+ template <MeshConcept MeshType>
+ class CellCenterReconstructionMatrixBuilder;
+
+ template <MeshConcept MeshType>
+ class ElementIntegralReconstructionMatrixBuilder;
+
+ template <MeshConcept MeshType>
+ class BoundaryIntegralReconstructionMatrixBuilder;
+
const PolynomialReconstructionDescriptor m_descriptor;
size_t _getNumberOfColumns(
--
GitLab