diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 25323391004f4004fd594444ae10d2ba55bf34dc..2e6e1fabbf51ddb49d54ea8fdc1d9b072ff1075c 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -25,6 +25,7 @@
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
 
+#include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
 #include <scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp>
 #include <scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp>
 
@@ -52,240 +53,6 @@ symmetrize_coordinates(const TinyVector<Dimension>& origin,
   return u - 2 * dot(u - origin, normal) * normal;
 }
 
-template <typename ConformTransformationT>
-void _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
-                             const ConformTransformationT& T,
-                             const TinyVector<2>& Xj,
-                             const double Vi,
-                             const size_t degree,
-                             const size_t dimension,
-                             const SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
-                             SmallArray<double>& mean_of_ejk);
-
-template <>
-void
-_computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
-                        const LineTransformation<2>& T,
-                        const TinyVector<2>& Xj,
-                        const double Vi,
-                        const size_t degree,
-                        const size_t dimension,
-                        const SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
-                        SmallArray<double>& mean_of_ejk)
-{
-  using Rd = TinyVector<2>;
-
-  const double velocity_perp_e1 = T.velocity()[1];
-
-  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;
-      inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1 / Vi;
-      for (; k <= degree; ++k) {
-        inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
-          x_xj * inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * ((1. * k) / (k + 1));
-      }
-
-      for (size_t i_y = 1; i_y <= degree; ++i_y) {
-        const size_t begin_i_y_1 = ((i_y - 1) * (2 * degree - i_y + 4)) / 2;
-        for (size_t l = 0; l <= degree - i_y; ++l) {
-          inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
-        }
-      }
-    }
-
-    for (size_t k = 1; k < dimension; ++k) {
-      mean_of_ejk[k - 1] += inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
-    }
-  }
-}
-
-template <MeshConcept MeshType>
-void
-_computeEjkMeanByBoundary(const MeshType& mesh,
-                          const TinyVector<2>& Xj,
-                          const CellId& cell_id,
-                          const auto& cell_to_face_matrix,
-                          const auto& face_to_node_matrix,
-                          const auto& cell_face_is_reversed,
-                          const CellValue<const double>& Vi,
-                          const size_t degree,
-                          const size_t basis_dimension,
-                          const SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
-                          SmallArray<double>& mean_of_ejk)
-{
-  const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(degree + 1));
-
-  auto xr = mesh.xr();
-  mean_of_ejk.fill(0);
-  auto cell_face_list = cell_to_face_matrix[cell_id];
-  for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
-    bool is_reversed = cell_face_is_reversed[cell_id][i_face];
-
-    const FaceId face_id = cell_face_list[i_face];
-    auto face_node_list  = face_to_node_matrix[face_id];
-    if (is_reversed) {
-      const LineTransformation<2> T{xr[face_node_list[1]], xr[face_node_list[0]]};
-      _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
-                              inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
-    } else {
-      const LineTransformation<2> T{xr[face_node_list[0]], xr[face_node_list[1]]};
-      _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
-                              inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
-    }
-  }
-}
-
-void
-_computeEjkMeanByBoundaryInSymmetricCell(const TinyVector<2>& origin,
-                                         const TinyVector<2>& normal,
-                                         const TinyVector<2>& Xj,
-                                         const CellId& cell_id,
-                                         const NodeValue<const TinyVector<2>>& xr,
-                                         const auto& cell_to_face_matrix,
-                                         const auto& face_to_node_matrix,
-                                         const auto& cell_face_is_reversed,
-                                         const CellValue<const double>& Vi,
-                                         const size_t degree,
-                                         const size_t basis_dimension,
-                                         const SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
-                                         SmallArray<double>& mean_of_ejk)
-{
-  const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(degree + 1));
-
-  mean_of_ejk.fill(0);
-  auto cell_face_list = cell_to_face_matrix[cell_id];
-  for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
-    bool is_reversed = cell_face_is_reversed[cell_id][i_face];
-
-    const FaceId face_id = cell_face_list[i_face];
-    auto face_node_list  = face_to_node_matrix[face_id];
-
-    const auto x0 = symmetrize_coordinates(origin, normal, xr[face_node_list[1]]);
-    const auto x1 = symmetrize_coordinates(origin, normal, xr[face_node_list[0]]);
-
-    if (is_reversed) {
-      const LineTransformation<2> T{x1, x0};
-      _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
-                              inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
-    } else {
-      const LineTransformation<2> T{x0, x1};
-      _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
-                              inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
-    }
-  }
-}
-
-template <MeshConcept MeshType>
-class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
-{
- private:
-  using Rd = TinyVector<MeshType::Dimension>;
-
-  const MeshType& m_mesh;
-  const size_t m_basis_dimension;
-  const size_t m_polynomial_degree;
-
-  const SmallArray<double> m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek;
-  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 CellToCellStencilArray& m_stencil_array;
-
-  const SmallArray<const Rd> m_symmetry_origin_list;
-  const SmallArray<const Rd> m_symmetry_normal_list;
-
-  const CellValue<const double> m_Vj;
-  const CellValue<const Rd> m_xj;
-  const NodeValue<const Rd> m_xr;
-
- public:
-  template <typename MatrixType>
-  void
-  build(const CellId cell_j_id, ShrinkMatrixView<MatrixType>& A)
-  {
-    const auto& stencil_cell_list = m_stencil_array[cell_j_id];
-
-    auto face_to_node_matrix   = m_mesh.connectivity().faceToNodeMatrix();
-    auto cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
-    const Rd& Xj               = m_xj[cell_j_id];
-
-    _computeEjkMeanByBoundary(m_mesh, Xj, cell_j_id, m_cell_to_face_matrix, face_to_node_matrix, cell_face_is_reversed,
-                              m_Vj, m_polynomial_degree, m_basis_dimension, m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
-                              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(m_mesh, Xj, cell_i_id, m_cell_to_face_matrix, face_to_node_matrix,
-                                cell_face_is_reversed, m_Vj, m_polynomial_degree, m_basis_dimension,
-                                m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek, 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_xr, m_cell_to_face_matrix, face_to_node_matrix,
-                                                 cell_face_is_reversed, m_Vj, m_polynomial_degree, m_basis_dimension,
-                                                 m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek, 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];
-        }
-      }
-    }
-  }
-
-  BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh,
-                                              const size_t polynomial_degree,
-                                              const SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
-                                              const SmallArray<double>& mean_j_of_ejk,
-                                              const SmallArray<double>& mean_i_of_ejk,
-                                              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{inv_Vj_alpha_p_1_wq_X_prime_orth_ek},
-      m_mean_j_of_ejk{mean_j_of_ejk},
-      m_mean_i_of_ejk{mean_i_of_ejk},
-
-      m_cell_to_face_matrix{mesh.connectivity().cellToFaceMatrix()},
-      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()}
-  {}
-  BoundaryIntegralReconstructionMatrixBuilder()  = default;
-  ~BoundaryIntegralReconstructionMatrixBuilder() = default;
-};
-
 class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
 {
  public:
@@ -882,6 +649,7 @@ PolynomialReconstruction::_build(
           if ((m_descriptor.integrationMethodType() == IntegrationMethodType::boundary) and
               (MeshType::Dimension == 2)) {
             if constexpr (MeshType::Dimension == 2) {
+#warning incorporate in reconstruction_matrix_builder
               SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek = inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[t];
               SmallArray<double>& mean_j_of_ejk                       = mean_j_of_ejk_pool[t];
               SmallArray<double>& mean_i_of_ejk                       = mean_i_of_ejk_pool[t];
@@ -894,50 +662,8 @@ PolynomialReconstruction::_build(
 
               boundary_integral_reconstruction_matrix_builder.build(cell_j_id, A);
 
-              auto face_to_node_matrix   = p_mesh->connectivity().faceToNodeMatrix();
-              auto cell_face_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
-              const Rd& Xj               = xj[cell_j_id];
-
-              _computeEjkMeanByBoundary(*p_mesh, Xj, cell_j_id, cell_to_face_matrix, face_to_node_matrix,
-                                        cell_face_is_reversed, Vj, m_descriptor.degree(), basis_dimension,
-                                        inv_Vj_alpha_p_1_wq_X_prime_orth_ek, 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(*p_mesh, Xj, cell_i_id, cell_to_face_matrix, face_to_node_matrix,
-                                          cell_face_is_reversed, Vj, m_descriptor.degree(), basis_dimension,
-                                          inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_i_of_ejk);
-
-                for (size_t l = 0; l < basis_dimension - 1; ++l) {
-                  A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
-                }
-              }
-              for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
-                   ++i_symmetry) {
-                auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
-                auto ghost_cell_list = ghost_stencil[cell_j_id];
-
-                const Rd& origin = symmetry_origin_list[i_symmetry];
-                const Rd& normal = 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, xr, cell_to_face_matrix, face_to_node_matrix,
-                                                           cell_face_is_reversed, Vj, m_descriptor.degree(),
-                                                           basis_dimension, inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
-                                                           mean_i_of_ejk);
-
-                  for (size_t l = 0; l < basis_dimension - 1; ++l) {
-                    A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
-                  }
-                }
-              }
             } else {
-              throw UnexpectedError("invalid mesh dimension");
+              throw NotImplementedError("invalid mesh dimension");
             }
           } else {
 
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a4a282fde6c970df2a1ed74db91d5e0baeef8169
--- /dev/null
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
@@ -0,0 +1,212 @@
+#ifndef BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
+#define BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
+
+#include <algebra/ShrinkMatrixView.hpp>
+#include <analysis/GaussLegendreQuadratureDescriptor.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/StencilArray.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/SmallArray.hpp>
+
+template <MeshConcept MeshType>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
+{
+ private:
+  using Rd = TinyVector<MeshType::Dimension>;
+
+  const MeshType& m_mesh;
+  const size_t m_basis_dimension;
+  const size_t m_polynomial_degree;
+
+  const SmallArray<double> m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek;
+  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;
+
+  const CellToCellStencilArray& m_stencil_array;
+
+  const SmallArray<const Rd> m_symmetry_origin_list;
+  const SmallArray<const Rd> m_symmetry_normal_list;
+
+  const CellValue<const double> m_Vj;
+  const CellValue<const Rd> m_xj;
+  const NodeValue<const Rd> m_xr;
+
+  void
+  _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
+                          const LineTransformation<2>& T,
+                          const Rd& Xj,
+                          const double Vi,
+                          SmallArray<double>& mean_of_ejk)
+  {
+    const double velocity_perp_e1 = T.velocity()[1];
+
+    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 / Vi;
+        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));
+        }
+
+        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));
+
+    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, m_Vj[cell_id], 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);
+      }
+    }
+  }
+
+  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));
+
+    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, m_Vj[cell_id], mean_of_ejk);
+      } else {
+        const LineTransformation<2> T{x0, x1};
+        _computeEjkBoundaryMean(quadrature, T, Xj, m_Vj[cell_id], mean_of_ejk);
+      }
+    }
+  }
+
+ public:
+  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");
+    }
+  }
+
+  BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh,
+                                              const size_t polynomial_degree,
+                                              const SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+                                              const SmallArray<double>& mean_j_of_ejk,
+                                              const SmallArray<double>& mean_i_of_ejk,
+                                              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{inv_Vj_alpha_p_1_wq_X_prime_orth_ek},
+      m_mean_j_of_ejk{mean_j_of_ejk},
+      m_mean_i_of_ejk{mean_i_of_ejk},
+      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()}
+  {}
+  BoundaryIntegralReconstructionMatrixBuilder()  = default;
+  ~BoundaryIntegralReconstructionMatrixBuilder() = default;
+};
+
+#endif   // BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
index 8477f57ea818faf8a0585d4267727e5f2ecc0f8e..8e1ce86607ff8ed3f913307bca652701e14a3265 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
@@ -134,7 +134,7 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
                   const double Vi,
                   SmallArray<double>& mean_of_ejk)
   {
-#warning Compute this one for all and rework design
+#warning Compute this once for all and rework design
     SmallArray<size_t> m_yz_row_index((m_polynomial_degree + 2) * (m_polynomial_degree + 1) / 2 + 1);
     SmallArray<size_t> m_z_triangle_index(m_polynomial_degree + 1);