diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2a481517f6aa70a8c4203e07445a7337e16f4912
--- /dev/null
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
@@ -0,0 +1,206 @@
+#include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/SymmetryUtils.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::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::_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_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1;
+ for (; k <= m_polynomial_degree; ++k) {
+ m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
+ x_xj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * m_antiderivative_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 = ((i_y - 1) * (2 * m_polynomial_degree - i_y + 4)) / 2;
+ for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l) {
+ m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
+ }
+ }
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] += m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
+ }
+ }
+}
+
+template <MeshConcept MeshTypeT>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkMeanByBoundary(
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ const auto& quadrature =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ const double inv_Vi = 1. / m_Vj[cell_id];
+
+ mean_of_ejk.fill(0);
+ auto cell_face_list = m_cell_to_face_matrix[cell_id];
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];
+
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = m_face_to_node_matrix[face_id];
+ if (is_reversed) {
+ 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);
+ }
+ }
+}
+
+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 auto& quadrature =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ const double inv_Vi = 1. / m_Vj[cell_id];
+
+ mean_of_ejk.fill(0);
+ auto cell_face_list = m_cell_to_face_matrix[cell_id];
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];
+
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = m_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 (is_reversed) {
+ 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);
+ }
+ }
+}
+
+template <MeshConcept MeshTypeT>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::build(
+ const CellId cell_j_id,
+ ShrinkMatrixView<SmallMatrix<double>>& A)
+{
+ if constexpr (MeshType::Dimension == 2) {
+ 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_inv_Vj_alpha_p_1_wq_X_prime_orth_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_face_matrix{mesh.connectivity().cellToFaceMatrix()},
+ m_face_to_node_matrix{mesh.connectivity().faceToNodeMatrix()},
+ m_cell_face_is_reversed{mesh.connectivity().cellFaceIsReversed()},
+ 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()}
+{
+ if constexpr (MeshType::Dimension == 2) {
+ SmallArray<double> antiderivative_coef(m_polynomial_degree + 1);
+ for (size_t k = 0; k < antiderivative_coef.size(); ++k) {
+ antiderivative_coef[k] = ((1. * k) / (k + 1));
+ }
+
+ m_antiderivative_coef = antiderivative_coef;
+ }
+}
+
+template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<2>>::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
index 105c4640e4d534c75c0f956b39d8347846f7ac48..df24d76a73097d13f6379681780d90aed6f49e1d 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
@@ -2,15 +2,12 @@
#define BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
#include <algebra/ShrinkMatrixView.hpp>
-#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <algebra/SmallMatrix.hpp>
#include <analysis/QuadratureFormula.hpp>
-#include <analysis/QuadratureManager.hpp>
#include <geometry/LineTransformation.hpp>
-#include <mesh/Connectivity.hpp>
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshDataManager.hpp>
+#include <mesh/ItemValue.hpp>
#include <mesh/StencilArray.hpp>
-#include <scheme/DiscreteFunctionDPk.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp>
@@ -48,103 +45,20 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
SmallArray<const double> m_antiderivative_coef;
- void
- _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
- const LineTransformation<2>& 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_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1;
- for (; k <= m_polynomial_degree; ++k) {
- m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
- x_xj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * m_antiderivative_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 = ((i_y - 1) * (2 * m_polynomial_degree - i_y + 4)) / 2;
- for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l) {
- m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
- }
- }
- }
-
- for (size_t k = 1; k < m_basis_dimension; ++k) {
- mean_of_ejk[k - 1] += m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
- }
- }
- }
+ 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)
- {
- const auto& quadrature =
- QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
-
- const double inv_Vi = 1. / m_Vj[cell_id];
-
- mean_of_ejk.fill(0);
- auto cell_face_list = m_cell_to_face_matrix[cell_id];
- for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];
-
- const FaceId face_id = cell_face_list[i_face];
- auto face_node_list = m_face_to_node_matrix[face_id];
- if (is_reversed) {
- 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);
- }
- }
- }
+ 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)
- {
- const auto& quadrature =
- QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
-
- const double inv_Vi = 1. / m_Vj[cell_id];
-
- mean_of_ejk.fill(0);
- auto cell_face_list = m_cell_to_face_matrix[cell_id];
- for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];
-
- const FaceId face_id = cell_face_list[i_face];
- auto face_node_list = m_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 (is_reversed) {
- 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);
- }
- }
- }
+ void _computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
+ const Rd& normal,
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk);
public:
PUGS_INLINE
@@ -154,79 +68,13 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
return m_mean_j_of_ejk;
}
- template <typename MatrixType>
- void
- build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A)
- {
- if constexpr (MeshType::Dimension == 2) {
- 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");
- }
- }
+ 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)
- : 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{m_basis_dimension},
- m_mean_j_of_ejk{m_basis_dimension - 1},
- m_mean_i_of_ejk{m_basis_dimension - 1},
- m_cell_to_face_matrix{mesh.connectivity().cellToFaceMatrix()},
- m_face_to_node_matrix{mesh.connectivity().faceToNodeMatrix()},
- m_cell_face_is_reversed{mesh.connectivity().cellFaceIsReversed()},
- 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_coef(m_polynomial_degree + 1);
- for (size_t k = 0; k < antiderivative_coef.size(); ++k) {
- antiderivative_coef[k] = ((1. * k) / (k + 1));
- }
-
- m_antiderivative_coef = antiderivative_coef;
- }
+ const CellToCellStencilArray& stencil_array);
BoundaryIntegralReconstructionMatrixBuilder() = default;
~BoundaryIntegralReconstructionMatrixBuilder() = default;
diff --git a/src/scheme/reconstruction_utils/CMakeLists.txt b/src/scheme/reconstruction_utils/CMakeLists.txt
index 66aaaa241b28ab025cec63f8bd8d9c1719b0f6f6..fa1395b2dd9615050ea557e6b83f96a3de994ca8 100644
--- a/src/scheme/reconstruction_utils/CMakeLists.txt
+++ b/src/scheme/reconstruction_utils/CMakeLists.txt
@@ -1,6 +1,7 @@
# ------------------- Source files --------------------
add_library(PugsSchemeReconstructionUtils
+ BoundaryIntegralReconstructionMatrixBuilder.cpp
ElementIntegralReconstructionMatrixBuilder.cpp
)