From ffe899397ff5ac711608fc73292e0b5fa78df817 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com> Date: Wed, 7 May 2025 17:38:17 +0200 Subject: [PATCH] Begin BoundaryIntegralReconstructionMatrixBuilder cleaning --- ...aryIntegralReconstructionMatrixBuilder.hpp | 32 +++++++++++++------ 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp index a4a282fd..350618f5 100644 --- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp +++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp @@ -41,14 +41,16 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder const CellValue<const Rd> m_xj; const NodeValue<const Rd> m_xr; + SmallArray<const double> m_antiderivative_coef; + void _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature, const LineTransformation<2>& T, const Rd& Xj, - const double Vi, + const double inv_Vi, SmallArray<double>& mean_of_ejk) { - const double velocity_perp_e1 = T.velocity()[1]; + 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); @@ -61,10 +63,10 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder { size_t k = 0; - m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1 / Vi; + 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] * ((1. * k) / (k + 1)); + 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) { @@ -87,6 +89,8 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder 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) { @@ -96,10 +100,10 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder 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, m_Vj[cell_id], mean_of_ejk); + _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, m_Vj[cell_id], mean_of_ejk); + _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk); } } } @@ -114,6 +118,8 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder 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) { @@ -127,10 +133,10 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder if (is_reversed) { const LineTransformation<2> T{x1, x0}; - _computeEjkBoundaryMean(quadrature, T, Xj, m_Vj[cell_id], mean_of_ejk); + _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk); } else { const LineTransformation<2> T{x0, x1}; - _computeEjkBoundaryMean(quadrature, T, Xj, m_Vj[cell_id], mean_of_ejk); + _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk); } } } @@ -204,7 +210,15 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder 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; + } + BoundaryIntegralReconstructionMatrixBuilder() = default; ~BoundaryIntegralReconstructionMatrixBuilder() = default; }; -- GitLab