Slightly redesign of BoundaryIntegralReconstructionMatrixBuilder

src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp

+#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&);

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,
+  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)
-  {
-    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);
+                               SmallArray<double>& mean_of_ejk);
 
-      const Rd X_Xj = T(xi_q) - Xj;
+  void _computeEjkMeanByBoundary(const Rd& Xj, const CellId& cell_id, SmallArray<double>& mean_of_ejk);
 
-      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];
-      }
-    }
-  }
-
-  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
-  _computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
+  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);
-      }
-    }
-  }
+                                                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;

src/scheme/reconstruction_utils/CMakeLists.txt

 # ------------------- Source files --------------------
 
 add_library(PugsSchemeReconstructionUtils
+  BoundaryIntegralReconstructionMatrixBuilder.cpp
   ElementIntegralReconstructionMatrixBuilder.cpp
 )