Begin BoundaryIntegralReconstructionMatrixBuilder cleaning

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;
 };