From 16ccf08d6d72831cbf610da17e3bfd4b80054ec8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 19 May 2025 14:37:14 +0200
Subject: [PATCH] Prepare boundary integration in 1d
---
src/scheme/PolynomialReconstruction.cpp | 1 +
...aryIntegralReconstructionMatrixBuilder.cpp | 70 +++++++++++++++----
...aryIntegralReconstructionMatrixBuilder.hpp | 12 ++--
...entIntegralReconstructionMatrixBuilder.cpp | 2 +
...entIntegralReconstructionMatrixBuilder.hpp | 5 +-
5 files changed, 67 insertions(+), 23 deletions(-)
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 6073da9c8..e29e8b8da 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -828,6 +828,7 @@ PolynomialReconstruction::_build(
return this->_build<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh,
discrete_function_variant_list);
}
+ std::clog << "falling back to element integration in 3d\n";
[[fallthrough]];
}
case IntegrationMethodType::element: {
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
index 2a481517f..0e3c01a7d 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
@@ -1,13 +1,53 @@
#include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
#include <geometry/SymmetryUtils.hpp>
-#include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshDataManager.hpp>
#include <scheme/DiscreteFunctionDPk.hpp>
+template <>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::ConnectivityData
+{
+ 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;
+ }
+
+ ConnectivityData(const Connectivity<2>& connectivity)
+ : m_cell_to_face_matrix{connectivity.cellToFaceMatrix()},
+ m_face_to_node_matrix{connectivity.faceToNodeMatrix()},
+ m_cell_face_is_reversed{connectivity.cellFaceIsReversed()}
+ {}
+
+ ~ConnectivityData() = default;
+};
+
template <MeshConcept MeshTypeT>
template <typename ConformTransformationT>
void
@@ -63,14 +103,15 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
const double inv_Vi = 1. / m_Vj[cell_id];
+ const auto& face_is_reversed = m_dimensional_data->cellFaceIsReversed()[cell_id];
+ const auto& cell_face_list = m_dimensional_data->cellToFaceMatrix()[cell_id];
+ const auto& face_to_node_matrix = m_dimensional_data->faceToNodeMatrix();
+
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) {
+ 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 {
@@ -94,18 +135,19 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
const double inv_Vi = 1. / m_Vj[cell_id];
+ const auto& face_is_reversed = m_dimensional_data->cellFaceIsReversed()[cell_id];
+ const auto& cell_face_list = m_dimensional_data->cellToFaceMatrix()[cell_id];
+ const auto& face_to_node_matrix = m_dimensional_data->faceToNodeMatrix();
+
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];
+ 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 (is_reversed) {
+ if (face_is_reversed[i_face]) {
const LineTransformation<2> T{x1, x0};
_computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
} else {
@@ -174,9 +216,9 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
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_dimensional_data{std::make_shared<ConnectivityData>(mesh.connectivity())},
+
m_stencil_array{stencil_array},
m_symmetry_origin_list{symmetry_origin_list},
m_symmetry_normal_list{symmetry_normal_list},
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
index df24d76a7..cb638fa2b 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
@@ -3,14 +3,14 @@
#include <algebra/ShrinkMatrixView.hpp>
#include <algebra/SmallMatrix.hpp>
-#include <analysis/QuadratureFormula.hpp>
-#include <geometry/LineTransformation.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/StencilArray.hpp>
-#include <mesh/SubItemValuePerItem.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp>
+template <size_t Dimension>
+class QuadratureFormula;
+
template <MeshConcept MeshTypeT>
class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
{
@@ -30,9 +30,8 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
SmallArray<double> m_mean_j_of_ejk;
SmallArray<double> m_mean_i_of_ejk;
- 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;
+ class ConnectivityData;
+ std::shared_ptr<ConnectivityData> m_dimensional_data;
const CellToCellStencilArray& m_stencil_array;
@@ -76,7 +75,6 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
const SmallArray<const Rd>& symmetry_normal_list,
const CellToCellStencilArray& stencil_array);
- BoundaryIntegralReconstructionMatrixBuilder() = default;
~BoundaryIntegralReconstructionMatrixBuilder() = default;
};
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
index faededf69..18306a8ec 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
@@ -2,6 +2,8 @@
#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>
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
index 1f4573438..7a3e8f7fc 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
@@ -3,14 +3,15 @@
#include <algebra/ShrinkMatrixView.hpp>
#include <algebra/SmallMatrix.hpp>
-#include <analysis/QuadratureFormula.hpp>
-#include <analysis/QuadratureManager.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
{
--
GitLab