diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 6616301877473e41ef3ea38909faf0da931c3072..25323391004f4004fd594444ae10d2ba55bf34dc 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -25,6 +25,9 @@
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
 
+#include <scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp>
+#include <scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp>
+
 template <size_t Dimension>
 PUGS_INLINE auto
 symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimension>& u)
@@ -49,611 +52,6 @@ symmetrize_coordinates(const TinyVector<Dimension>& origin,
   return u - 2 * dot(u - origin, normal) * normal;
 }
 
-template <MeshConcept MeshType>
-class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder
-{
- private:
-  using Rd = TinyVector<MeshType::Dimension>;
-
-  const size_t m_basis_dimension;
-  const SmallArray<const Rd> m_symmetry_origin_list;
-  const SmallArray<const Rd> m_symmetry_normal_list;
-
-  const CellToCellStencilArray& m_stencil_array;
-  const CellValue<const Rd> m_xj;
-
- 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];
-
-    const Rd& Xj = m_xj[cell_j_id];
-    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];
-      const Rd Xi_Xj         = m_xj[cell_i_id] - Xj;
-      for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
-        A(index, l) = Xi_Xj[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];
-        const Rd Xi_Xj         = symmetrize_coordinates(origin, normal, m_xj[cell_i_id]) - Xj;
-        for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
-          A(index, l) = Xi_Xj[l];
-        }
-      }
-    }
-  }
-
-  CellCenterReconstructionMatrixBuilder(const SmallArray<const Rd>& symmetry_origin_list,
-                                        const SmallArray<const Rd>& symmetry_normal_list,
-                                        const CellToCellStencilArray& stencil_array,
-                                        const CellValue<const Rd>& xj)
-    : m_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(1)},
-      m_symmetry_origin_list{symmetry_origin_list},
-      m_symmetry_normal_list{symmetry_normal_list},
-      m_stencil_array{stencil_array},
-      m_xj{xj}
-  {}
-
-  ~CellCenterReconstructionMatrixBuilder() = default;
-};
-
-template <typename ConformTransformationT>
-void
-_computeEjkMean(const QuadratureFormula<1>& quadrature,
-                const ConformTransformationT& T,
-                const TinyVector<1>& Xj,
-                const double Vi,
-                const size_t degree,
-                const size_t dimension,
-                const SmallArray<double>& inv_Vj_wq_detJ_ek,
-                SmallArray<double>& mean_of_ejk)
-{
-  using Rd = TinyVector<1>;
-  mean_of_ejk.fill(0);
-
-  for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
-    const double wq   = quadrature.weight(i_q);
-    const Rd& xi_q    = quadrature.point(i_q);
-    const double detT = T.jacobianDeterminant();
-
-    const Rd X_Xj = T(xi_q) - Xj;
-
-    const double x_xj = X_Xj[0];
-
-    {
-      size_t k               = 0;
-      inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
-      for (; k <= degree; ++k) {
-        inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
-      }
-    }
-
-    for (size_t k = 1; k < dimension; ++k) {
-      mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
-    }
-  }
-}
-
-template <typename ConformTransformationT>
-void
-_computeEjkMean(const QuadratureFormula<2>& quadrature,
-                const ConformTransformationT& T,
-                const TinyVector<2>& Xj,
-                const double Vi,
-                const size_t degree,
-                const size_t dimension,
-                const SmallArray<double>& inv_Vj_wq_detJ_ek,
-                SmallArray<double>& mean_of_ejk)
-{
-  using Rd = TinyVector<2>;
-
-  mean_of_ejk.fill(0);
-
-  for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
-    const double wq   = quadrature.weight(i_q);
-    const Rd& xi_q    = quadrature.point(i_q);
-    const double detT = [&] {
-      if constexpr (std::is_same_v<TriangleTransformation<2>, std::decay_t<decltype(T)>>) {
-        return T.jacobianDeterminant();
-      } else {
-        return T.jacobianDeterminant(xi_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_wq_detJ_ek[k++] = wq * detT / Vi;
-      for (; k <= degree; ++k) {
-        inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
-      }
-
-      for (size_t i_y = 1; i_y <= degree; ++i_y) {
-#warning store y_row index and size
-        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_wq_detJ_ek[k++] = y_yj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
-        }
-      }
-    }
-
-    for (size_t k = 1; k < dimension; ++k) {
-      mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
-    }
-  }
-}
-
-template <typename ConformTransformationT>
-void
-_computeEjkMean(const QuadratureFormula<3>& quadrature,
-                const ConformTransformationT& T,
-                const TinyVector<3>& Xj,
-                const double Vi,
-                const size_t degree,
-                const size_t dimension,
-                const SmallArray<double>& inv_Vj_wq_detJ_ek,
-                SmallArray<double>& mean_of_ejk)
-{
-#warning Compute this one for all and rework design
-  SmallArray<size_t> m_yz_row_index((degree + 2) * (degree + 1) / 2 + 1);
-  SmallArray<size_t> m_z_triangle_index(degree + 1);
-
-  {
-    size_t i_z  = 0;
-    size_t i_yz = 0;
-
-    m_yz_row_index[i_yz++] = 0;
-    for (ssize_t n = degree + 1; n >= 1; --n) {
-      m_z_triangle_index[i_z++] = i_yz - 1;
-      for (ssize_t i = n; i >= 1; --i) {
-        m_yz_row_index[i_yz] = m_yz_row_index[i_yz - 1] + i;
-        ++i_yz;
-      }
-    }
-  }
-
-  SmallArray<size_t> m_yz_row_size{m_yz_row_index.size() - 1};
-  for (size_t i = 0; i < m_yz_row_size.size(); ++i) {
-    m_yz_row_size[i] = m_yz_row_index[i + 1] - m_yz_row_index[i];
-  }
-
-  using Rd = TinyVector<3>;
-  mean_of_ejk.fill(0);
-
-  for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
-    const double wq   = quadrature.weight(i_q);
-    const Rd& xi_q    = quadrature.point(i_q);
-    const double detT = [&] {
-      if constexpr (std::is_same_v<TetrahedronTransformation, std::decay_t<decltype(T)>>) {
-        return T.jacobianDeterminant();
-      } else {
-        return T.jacobianDeterminant(xi_q);
-      }
-    }();
-
-    const Rd X_Xj = T(xi_q) - Xj;
-
-    const double x_xj = X_Xj[0];
-    const double y_yj = X_Xj[1];
-    const double z_zj = X_Xj[2];
-
-    {
-      size_t k               = 0;
-      inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
-      for (; k <= degree; ++k) {
-        inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
-      }
-
-      for (size_t i_y = 1; i_y <= degree; ++i_y) {
-        const size_t begin_i_y_1 = m_yz_row_index[i_y - 1];
-        const size_t nb_monoms   = m_yz_row_size[i_y];
-        for (size_t l = 0; l < nb_monoms; ++l) {
-          inv_Vj_wq_detJ_ek[k++] = y_yj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
-        }
-      }
-
-      for (size_t i_z = 1; i_z <= degree; ++i_z) {
-        const size_t nb_y      = m_yz_row_size[m_z_triangle_index[i_z]];
-        const size_t index_z   = m_z_triangle_index[i_z];
-        const size_t index_z_1 = m_z_triangle_index[i_z - 1];
-        for (size_t i_y = 0; i_y < nb_y; ++i_y) {
-          const size_t begin_i_yz_1 = m_yz_row_index[index_z_1 + i_y];
-          const size_t nb_monoms    = m_yz_row_size[index_z + i_y];
-          for (size_t l = 0; l < nb_monoms; ++l) {
-            inv_Vj_wq_detJ_ek[k++] = z_zj * inv_Vj_wq_detJ_ek[begin_i_yz_1 + l];
-          }
-        }
-      }
-    }
-
-    for (size_t k = 1; k < dimension; ++k) {
-      mean_of_ejk[k - 1] += inv_Vj_wq_detJ_ek[k];
-    }
-  }
-}
-
-template <typename NodeListT, size_t Dimension>
-void _computeEjkMean(const TinyVector<Dimension>& Xj,
-                     const NodeValue<const TinyVector<Dimension>>& xr,
-                     const NodeListT& node_list,
-                     const CellType& cell_type,
-                     const double Vi,
-                     const size_t degree,
-                     const size_t basis_dimension,
-                     const SmallArray<double>& inv_Vj_wq_detJ_ek,
-                     SmallArray<double>& mean_of_ejk);
-
-template <typename NodeListT>
-void
-_computeEjkMean(const TinyVector<1>& Xj,
-                const NodeValue<const TinyVector<1>>& xr,
-                const NodeListT& node_list,
-                const CellType& cell_type,
-                const double Vi,
-                const size_t degree,
-                const size_t basis_dimension,
-                const SmallArray<double>& inv_Vj_wq_detJ_ek,
-                SmallArray<double>& mean_of_ejk)
-{
-  if (cell_type == CellType::Line) {
-    const LineTransformation<1> T{xr[node_list[0]], xr[node_list[1]]};
-
-    const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{degree});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-
-  } else {
-    throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
-  }
-}
-
-template <typename NodeListT>
-void
-_computeEjkMeanInSymmetricCell(const TinyVector<1>& origin,
-                               const TinyVector<1>& normal,
-                               const TinyVector<1>& Xj,
-                               const NodeValue<const TinyVector<1>>& xr,
-                               const NodeListT& node_list,
-                               const CellType& cell_type,
-                               const double Vi,
-                               const size_t degree,
-                               const size_t basis_dimension,
-                               const SmallArray<double>& inv_Vj_wq_detJ_ek,
-                               SmallArray<double>& mean_of_ejk)
-{
-  if (cell_type == CellType::Line) {
-    const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
-    const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
-
-    const LineTransformation<1> T{x0, x1};
-
-    const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{degree});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-
-  } else {
-    throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
-  }
-}
-
-template <typename NodeListT>
-void
-_computeEjkMean(const TinyVector<2>& Xj,
-                const NodeValue<const TinyVector<2>>& xr,
-                const NodeListT& node_list,
-                const CellType& cell_type,
-                const double Vi,
-                const size_t degree,
-                const size_t basis_dimension,
-                const SmallArray<double>& inv_Vj_wq_detJ_ek,
-                SmallArray<double>& mean_of_ejk)
-{
-  switch (cell_type) {
-  case CellType::Triangle: {
-    const TriangleTransformation<2> T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]]};
-    const auto& quadrature = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{degree});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-    break;
-  }
-  case CellType::Quadrangle: {
-    const SquareTransformation<2> T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]], xr[node_list[3]]};
-    const auto& quadrature =
-      QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor{degree + 1});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-    break;
-  }
-  default: {
-    throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
-  }
-  }
-}
-
-template <typename NodeListT>
-void
-_computeEjkMeanInSymmetricCell(const TinyVector<2>& origin,
-                               const TinyVector<2>& normal,
-                               const TinyVector<2>& Xj,
-                               const NodeValue<const TinyVector<2>>& xr,
-                               const NodeListT& node_list,
-                               const CellType& cell_type,
-                               const double Vi,
-                               const size_t degree,
-                               const size_t basis_dimension,
-                               const SmallArray<double>& inv_Vj_wq_detJ_ek,
-                               SmallArray<double>& mean_of_ejk)
-{
-  switch (cell_type) {
-  case CellType::Triangle: {
-    const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
-    const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
-    const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
-
-    const TriangleTransformation<2> T{x0, x1, x2};
-    const auto& quadrature = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{degree});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-    break;
-  }
-  case CellType::Quadrangle: {
-    const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
-    const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
-    const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
-    const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
-
-    const SquareTransformation<2> T{x0, x1, x2, x3};
-    const auto& quadrature =
-      QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor{degree + 1});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-    break;
-  }
-  default: {
-    throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
-  }
-  }
-}
-
-template <typename NodeListT>
-void
-_computeEjkMean(const TinyVector<3>& Xj,
-                const NodeValue<const TinyVector<3>>& xr,
-                const NodeListT& node_list,
-                const CellType& cell_type,
-                const double Vi,
-                const size_t degree,
-                const size_t basis_dimension,
-                const SmallArray<double>& inv_Vj_wq_detJ_ek,
-                SmallArray<double>& mean_of_ejk)
-{
-  if (cell_type == CellType::Tetrahedron) {
-    const TetrahedronTransformation T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]], xr[node_list[3]]};
-
-    const auto& quadrature = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{degree});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-
-  } else if (cell_type == CellType::Prism) {
-    const PrismTransformation T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]],   //
-                                xr[node_list[3]], xr[node_list[4]], xr[node_list[5]]};
-
-    const auto& quadrature = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{degree + 1});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-
-  } else if (cell_type == CellType::Pyramid) {
-    const PyramidTransformation T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]], xr[node_list[3]],
-                                  xr[node_list[4]]};
-
-    const auto& quadrature = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{degree + 1});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-
-  } else if (cell_type == CellType::Hexahedron) {
-    const CubeTransformation T{xr[node_list[0]], xr[node_list[1]], xr[node_list[2]], xr[node_list[3]],
-                               xr[node_list[4]], xr[node_list[5]], xr[node_list[6]], xr[node_list[7]]};
-
-    const auto& quadrature =
-      QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{degree + 1});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-
-  } else {
-    throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
-  }
-}
-
-template <typename NodeListT>
-void
-_computeEjkMeanInSymmetricCell(const TinyVector<3>& origin,
-                               const TinyVector<3>& normal,
-                               const TinyVector<3>& Xj,
-                               const NodeValue<const TinyVector<3>>& xr,
-                               const NodeListT& node_list,
-                               const CellType& cell_type,
-                               const double Vi,
-                               const size_t degree,
-                               const size_t basis_dimension,
-                               const SmallArray<double>& inv_Vj_wq_detJ_ek,
-                               SmallArray<double>& mean_of_ejk)
-{
-  if (cell_type == CellType::Tetrahedron) {
-    const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
-    const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
-    const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
-    const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
-
-    const TetrahedronTransformation T{x0, x1, x2, x3};
-
-    const auto& quadrature = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{degree});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-
-  } else if (cell_type == CellType::Prism) {
-    const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
-    const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
-    const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
-    const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[4]]);
-    const auto x4 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
-    const auto x5 = symmetrize_coordinates(origin, normal, xr[node_list[5]]);
-
-    const PrismTransformation T{x0, x1, x2,   //
-                                x3, x4, x5};
-
-    const auto& quadrature = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{degree + 1});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-
-  } else if (cell_type == CellType::Pyramid) {
-    const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
-    const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
-    const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
-    const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
-    const auto x4 = symmetrize_coordinates(origin, normal, xr[node_list[4]]);
-    const PyramidTransformation T{x0, x1, x2, x3, x4};
-
-    const auto& quadrature = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{degree + 1});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-
-  } else if (cell_type == CellType::Hexahedron) {
-    const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
-    const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
-    const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
-    const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
-    const auto x4 = symmetrize_coordinates(origin, normal, xr[node_list[7]]);
-    const auto x5 = symmetrize_coordinates(origin, normal, xr[node_list[6]]);
-    const auto x6 = symmetrize_coordinates(origin, normal, xr[node_list[5]]);
-    const auto x7 = symmetrize_coordinates(origin, normal, xr[node_list[4]]);
-
-    const CubeTransformation T{x0, x1, x2, x3, x4, x5, x6, x7};
-
-    const auto& quadrature =
-      QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{degree + 1});
-
-    _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
-
-  } else {
-    throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
-  }
-}
-
-template <MeshConcept MeshType>
-class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
-{
- private:
-  using Rd = TinyVector<MeshType::Dimension>;
-
-  const size_t m_basis_dimension;
-  const size_t m_polynomial_degree;
-
-  const SmallArray<double> m_inv_Vj_wq_detJ_ek;
-  SmallArray<double> m_mean_j_of_ejk;
-  SmallArray<double> m_mean_i_of_ejk;
-
-  const ItemToItemMatrix<ItemType::cell, ItemType::node> m_cell_to_node_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 CellType> m_cell_type;
-  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];
-
-    const Rd& Xj = m_xj[cell_j_id];
-
-    _computeEjkMean(Xj, m_xr, m_cell_to_node_matrix[cell_j_id], m_cell_type[cell_j_id], m_Vj[cell_j_id],
-                    m_polynomial_degree, m_basis_dimension, m_inv_Vj_wq_detJ_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];
-
-      _computeEjkMean(Xj, m_xr, m_cell_to_node_matrix[cell_i_id], m_cell_type[cell_i_id], m_Vj[cell_i_id],
-                      m_polynomial_degree, m_basis_dimension, m_inv_Vj_wq_detJ_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];
-
-        _computeEjkMeanInSymmetricCell(origin, normal,   //
-                                       Xj, m_xr, m_cell_to_node_matrix[cell_i_id], m_cell_type[cell_i_id],
-                                       m_Vj[cell_i_id], m_polynomial_degree, m_basis_dimension, m_inv_Vj_wq_detJ_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];
-        }
-      }
-    }
-  }
-
-  ElementIntegralReconstructionMatrixBuilder(const MeshType& mesh,
-                                             const size_t polynomial_degree,
-                                             const SmallArray<double>& inv_Vj_wq_detJ_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_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(
-        polynomial_degree)},
-      m_polynomial_degree{polynomial_degree},
-
-      m_inv_Vj_wq_detJ_ek{inv_Vj_wq_detJ_ek},
-      m_mean_j_of_ejk{mean_j_of_ejk},
-      m_mean_i_of_ejk{mean_i_of_ejk},
-
-      m_cell_to_node_matrix{mesh.connectivity().cellToNodeMatrix()},
-      m_stencil_array{stencil_array},
-      m_symmetry_origin_list{symmetry_origin_list},
-      m_symmetry_normal_list{symmetry_normal_list},
-      m_cell_type{mesh.connectivity().cellType()},
-      m_Vj{MeshDataManager::instance().getMeshData(mesh).Vj()},
-      m_xj{MeshDataManager::instance().getMeshData(mesh).xj()},
-      m_xr{mesh.xr()}
-  {}
-
-  ~ElementIntegralReconstructionMatrixBuilder() = default;
-};
-
 template <typename ConformTransformationT>
 void _computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
                              const ConformTransformationT& T,
diff --git a/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9f3ee95247f8bfbd870dfbaf62983b7b90a16aa4
--- /dev/null
+++ b/src/scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp
@@ -0,0 +1,71 @@
+#ifndef CELL_CENTER_RECONSTRUCTION_MATRIX_BUILDER_HPP
+#define CELL_CENTER_RECONSTRUCTION_MATRIX_BUILDER_HPP
+
+#include <algebra/ShrinkMatrixView.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/StencilArray.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/SmallArray.hpp>
+
+template <MeshConcept MeshType>
+class PolynomialReconstruction::CellCenterReconstructionMatrixBuilder
+{
+ private:
+  using Rd = TinyVector<MeshType::Dimension>;
+
+  const size_t m_basis_dimension;
+  const SmallArray<const Rd> m_symmetry_origin_list;
+  const SmallArray<const Rd> m_symmetry_normal_list;
+
+  const CellToCellStencilArray& m_stencil_array;
+  const CellValue<const Rd> m_xj;
+
+ 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];
+
+    const Rd& Xj = m_xj[cell_j_id];
+    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];
+      const Rd Xi_Xj         = m_xj[cell_i_id] - Xj;
+      for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
+        A(index, l) = Xi_Xj[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];
+        const Rd Xi_Xj         = symmetrize_coordinates(origin, normal, m_xj[cell_i_id]) - Xj;
+        for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
+          A(index, l) = Xi_Xj[l];
+        }
+      }
+    }
+  }
+
+  CellCenterReconstructionMatrixBuilder(const SmallArray<const Rd>& symmetry_origin_list,
+                                        const SmallArray<const Rd>& symmetry_normal_list,
+                                        const CellToCellStencilArray& stencil_array,
+                                        const CellValue<const Rd>& xj)
+    : m_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(1)},
+      m_symmetry_origin_list{symmetry_origin_list},
+      m_symmetry_normal_list{symmetry_normal_list},
+      m_stencil_array{stencil_array},
+      m_xj{xj}
+  {}
+
+  ~CellCenterReconstructionMatrixBuilder() = default;
+};
+
+#endif   // CELL_CENTER_RECONSTRUCTION_MATRIX_BUILDER_HPP
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8477f57ea818faf8a0585d4267727e5f2ecc0f8e
--- /dev/null
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
@@ -0,0 +1,534 @@
+#ifndef ELEMENT_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
+#define ELEMENT_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
+
+#include <algebra/ShrinkMatrixView.hpp>
+#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>
+#include <geometry/PyramidTransformation.hpp>
+#include <geometry/SquareTransformation.hpp>
+#include <geometry/TetrahedronTransformation.hpp>
+#include <geometry/TriangleTransformation.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::ElementIntegralReconstructionMatrixBuilder
+{
+ private:
+  using Rd = TinyVector<MeshType::Dimension>;
+
+  const size_t m_basis_dimension;
+  const size_t m_polynomial_degree;
+
+  const SmallArray<double> m_inv_Vj_wq_detJ_ek;
+  SmallArray<double> m_mean_j_of_ejk;
+  SmallArray<double> m_mean_i_of_ejk;
+
+  const ItemToItemMatrix<ItemType::cell, ItemType::node> m_cell_to_node_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 CellType> m_cell_type;
+  const CellValue<const double> m_Vj;
+  const CellValue<const Rd> m_xj;
+  const NodeValue<const Rd> m_xr;
+
+  template <typename ConformTransformationT>
+  void
+  _computeEjkMean(const QuadratureFormula<1>& quadrature,
+                  const ConformTransformationT& T,
+                  const TinyVector<1>& Xj,
+                  const double Vi,
+                  SmallArray<double>& mean_of_ejk)
+  {
+    mean_of_ejk.fill(0);
+
+    for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+      const double wq   = quadrature.weight(i_q);
+      const Rd& xi_q    = quadrature.point(i_q);
+      const double detT = T.jacobianDeterminant();
+
+      const Rd X_Xj = T(xi_q) - Xj;
+
+      const double x_xj = X_Xj[0];
+
+      {
+        size_t k                 = 0;
+        m_inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
+        for (; k <= m_polynomial_degree; ++k) {
+          m_inv_Vj_wq_detJ_ek[k] = x_xj * m_inv_Vj_wq_detJ_ek[k - 1];
+        }
+      }
+
+      for (size_t k = 1; k < m_basis_dimension; ++k) {
+        mean_of_ejk[k - 1] += m_inv_Vj_wq_detJ_ek[k];
+      }
+    }
+  }
+
+  template <typename ConformTransformationT>
+  void
+  _computeEjkMean(const QuadratureFormula<2>& quadrature,
+                  const ConformTransformationT& T,
+                  const TinyVector<2>& Xj,
+                  const double Vi,
+                  SmallArray<double>& mean_of_ejk)
+  {
+    mean_of_ejk.fill(0);
+
+    for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+      const double wq   = quadrature.weight(i_q);
+      const Rd& xi_q    = quadrature.point(i_q);
+      const double detT = [&] {
+        if constexpr (std::is_same_v<TriangleTransformation<2>, std::decay_t<decltype(T)>>) {
+          return T.jacobianDeterminant();
+        } else {
+          return T.jacobianDeterminant(xi_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_wq_detJ_ek[k++] = wq * detT / Vi;
+        for (; k <= m_polynomial_degree; ++k) {
+          m_inv_Vj_wq_detJ_ek[k] = x_xj * m_inv_Vj_wq_detJ_ek[k - 1];
+        }
+
+        for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
+#warning store y_row index and size
+          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_wq_detJ_ek[k++] = y_yj * m_inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
+          }
+        }
+      }
+
+      for (size_t k = 1; k < m_basis_dimension; ++k) {
+        mean_of_ejk[k - 1] += m_inv_Vj_wq_detJ_ek[k];
+      }
+    }
+  }
+
+  template <typename ConformTransformationT>
+  void
+  _computeEjkMean(const QuadratureFormula<3>& quadrature,
+                  const ConformTransformationT& T,
+                  const TinyVector<3>& Xj,
+                  const double Vi,
+                  SmallArray<double>& mean_of_ejk)
+  {
+#warning Compute this one 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);
+
+    {
+      size_t i_z  = 0;
+      size_t i_yz = 0;
+
+      m_yz_row_index[i_yz++] = 0;
+      for (ssize_t n = m_polynomial_degree + 1; n >= 1; --n) {
+        m_z_triangle_index[i_z++] = i_yz - 1;
+        for (ssize_t i = n; i >= 1; --i) {
+          m_yz_row_index[i_yz] = m_yz_row_index[i_yz - 1] + i;
+          ++i_yz;
+        }
+      }
+    }
+
+    SmallArray<size_t> m_yz_row_size{m_yz_row_index.size() - 1};
+    for (size_t i = 0; i < m_yz_row_size.size(); ++i) {
+      m_yz_row_size[i] = m_yz_row_index[i + 1] - m_yz_row_index[i];
+    }
+
+    mean_of_ejk.fill(0);
+
+    for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+      const double wq   = quadrature.weight(i_q);
+      const Rd& xi_q    = quadrature.point(i_q);
+      const double detT = [&] {
+        if constexpr (std::is_same_v<TetrahedronTransformation, std::decay_t<decltype(T)>>) {
+          return T.jacobianDeterminant();
+        } else {
+          return T.jacobianDeterminant(xi_q);
+        }
+      }();
+
+      const Rd X_Xj = T(xi_q) - Xj;
+
+      const double x_xj = X_Xj[0];
+      const double y_yj = X_Xj[1];
+      const double z_zj = X_Xj[2];
+
+      {
+        size_t k                 = 0;
+        m_inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
+        for (; k <= m_polynomial_degree; ++k) {
+          m_inv_Vj_wq_detJ_ek[k] = x_xj * m_inv_Vj_wq_detJ_ek[k - 1];
+        }
+
+        for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
+          const size_t begin_i_y_1 = m_yz_row_index[i_y - 1];
+          const size_t nb_monoms   = m_yz_row_size[i_y];
+          for (size_t l = 0; l < nb_monoms; ++l) {
+            m_inv_Vj_wq_detJ_ek[k++] = y_yj * m_inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
+          }
+        }
+
+        for (size_t i_z = 1; i_z <= m_polynomial_degree; ++i_z) {
+          const size_t nb_y      = m_yz_row_size[m_z_triangle_index[i_z]];
+          const size_t index_z   = m_z_triangle_index[i_z];
+          const size_t index_z_1 = m_z_triangle_index[i_z - 1];
+          for (size_t i_y = 0; i_y < nb_y; ++i_y) {
+            const size_t begin_i_yz_1 = m_yz_row_index[index_z_1 + i_y];
+            const size_t nb_monoms    = m_yz_row_size[index_z + i_y];
+            for (size_t l = 0; l < nb_monoms; ++l) {
+              m_inv_Vj_wq_detJ_ek[k++] = z_zj * m_inv_Vj_wq_detJ_ek[begin_i_yz_1 + l];
+            }
+          }
+        }
+      }
+
+      for (size_t k = 1; k < m_basis_dimension; ++k) {
+        mean_of_ejk[k - 1] += m_inv_Vj_wq_detJ_ek[k];
+      }
+    }
+  }
+
+  void
+  _computeEjkMean(const TinyVector<1>& Xj, const CellId& cell_i_id, SmallArray<double>& mean_of_ejk)
+  {
+    const CellType cell_type = m_cell_type[cell_i_id];
+    const auto node_list     = m_cell_to_node_matrix[cell_i_id];
+    const double Vi          = m_Vj[cell_i_id];
+
+    if (m_cell_type[cell_i_id] == CellType::Line) {
+      const LineTransformation<1> T{m_xr[node_list[0]], m_xr[node_list[1]]};
+
+      const auto& quadrature =
+        QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+    } else {
+      throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+    }
+  }
+
+  void
+  _computeEjkMean(const TinyVector<2>& Xj, const CellId& cell_i_id, SmallArray<double>& mean_of_ejk)
+  {
+    const CellType cell_type = m_cell_type[cell_i_id];
+    const auto node_list     = m_cell_to_node_matrix[cell_i_id];
+    const double Vi          = m_Vj[cell_i_id];
+
+    switch (cell_type) {
+    case CellType::Triangle: {
+      const TriangleTransformation<2> T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]]};
+      const auto& quadrature =
+        QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{m_polynomial_degree});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+      break;
+    }
+    case CellType::Quadrangle: {
+      const SquareTransformation<2> T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]], m_xr[node_list[3]]};
+      const auto& quadrature =
+        QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree + 1});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+      break;
+    }
+    default: {
+      throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+    }
+    }
+  }
+
+  void
+  _computeEjkMean(const TinyVector<3>& Xj, const CellId& cell_i_id, SmallArray<double>& mean_of_ejk)
+  {
+    const CellType cell_type = m_cell_type[cell_i_id];
+    const auto node_list     = m_cell_to_node_matrix[cell_i_id];
+    const double Vi          = m_Vj[cell_i_id];
+
+    switch (cell_type) {
+    case CellType::Tetrahedron: {
+      const TetrahedronTransformation T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]], m_xr[node_list[3]]};
+
+      const auto& quadrature =
+        QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{m_polynomial_degree});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+      break;
+    }
+    case CellType::Prism: {
+      const PrismTransformation T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]],   //
+                                  m_xr[node_list[3]], m_xr[node_list[4]], m_xr[node_list[5]]};
+
+      const auto& quadrature =
+        QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{m_polynomial_degree + 1});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+      break;
+    }
+    case CellType::Pyramid: {
+      const PyramidTransformation T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]], m_xr[node_list[3]],
+                                    m_xr[node_list[4]]};
+
+      const auto& quadrature =
+        QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{m_polynomial_degree + 1});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+      break;
+    }
+    case CellType::Hexahedron: {
+      const CubeTransformation T{m_xr[node_list[0]], m_xr[node_list[1]], m_xr[node_list[2]], m_xr[node_list[3]],
+                                 m_xr[node_list[4]], m_xr[node_list[5]], m_xr[node_list[6]], m_xr[node_list[7]]};
+
+      const auto& quadrature =
+        QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree + 1});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+      break;
+    }
+    default: {
+      throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+    }
+    }
+  }
+
+  void
+  _computeEjkMeanInSymmetricCell(const TinyVector<1>& origin,
+                                 const TinyVector<1>& normal,
+                                 const TinyVector<1>& Xj,
+                                 const CellId& cell_i_id,
+                                 SmallArray<double>& mean_of_ejk)
+  {
+    auto node_list           = m_cell_to_node_matrix[cell_i_id];
+    const CellType cell_type = m_cell_type[cell_i_id];
+    const double Vi          = m_Vj[cell_i_id];
+
+    if (cell_type == CellType::Line) {
+      const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+      const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+
+      const LineTransformation<1> T{x0, x1};
+
+      const auto& quadrature =
+        QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+    } else {
+      throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+    }
+  }
+
+  void
+  _computeEjkMeanInSymmetricCell(const TinyVector<2>& origin,
+                                 const TinyVector<2>& normal,
+                                 const TinyVector<2>& Xj,
+                                 const CellId& cell_i_id,
+                                 SmallArray<double>& mean_of_ejk)
+  {
+    auto node_list           = m_cell_to_node_matrix[cell_i_id];
+    const CellType cell_type = m_cell_type[cell_i_id];
+    const double Vi          = m_Vj[cell_i_id];
+
+    switch (cell_type) {
+    case CellType::Triangle: {
+      const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+      const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+      const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+
+      const TriangleTransformation<2> T{x0, x1, x2};
+      const auto& quadrature =
+        QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{m_polynomial_degree});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+      break;
+    }
+    case CellType::Quadrangle: {
+      const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[3]]);
+      const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+      const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+      const auto x3 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+
+      const SquareTransformation<2> T{x0, x1, x2, x3};
+      const auto& quadrature =
+        QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree + 1});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+      break;
+    }
+    default: {
+      throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+    }
+    }
+  }
+
+  void
+  _computeEjkMeanInSymmetricCell(const TinyVector<3>& origin,
+                                 const TinyVector<3>& normal,
+                                 const TinyVector<3>& Xj,
+                                 const CellId cell_i_id,
+                                 SmallArray<double>& mean_of_ejk)
+  {
+    auto node_list           = m_cell_to_node_matrix[cell_i_id];
+    const CellType cell_type = m_cell_type[cell_i_id];
+    const double Vi          = m_Vj[cell_i_id];
+
+    if (cell_type == CellType::Tetrahedron) {
+      const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+      const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+      const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+      const auto x3 = symmetrize_coordinates(origin, normal, m_xr[node_list[3]]);
+
+      const TetrahedronTransformation T{x0, x1, x2, x3};
+
+      const auto& quadrature =
+        QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{m_polynomial_degree});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+    } else if (cell_type == CellType::Prism) {
+      const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+      const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+      const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+      const auto x3 = symmetrize_coordinates(origin, normal, m_xr[node_list[4]]);
+      const auto x4 = symmetrize_coordinates(origin, normal, m_xr[node_list[3]]);
+      const auto x5 = symmetrize_coordinates(origin, normal, m_xr[node_list[5]]);
+
+      const PrismTransformation T{x0, x1, x2,   //
+                                  x3, x4, x5};
+
+      const auto& quadrature =
+        QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{m_polynomial_degree + 1});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+    } else if (cell_type == CellType::Pyramid) {
+      const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[3]]);
+      const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+      const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+      const auto x3 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+      const auto x4 = symmetrize_coordinates(origin, normal, m_xr[node_list[4]]);
+      const PyramidTransformation T{x0, x1, x2, x3, x4};
+
+      const auto& quadrature =
+        QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{m_polynomial_degree + 1});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+    } else if (cell_type == CellType::Hexahedron) {
+      const auto x0 = symmetrize_coordinates(origin, normal, m_xr[node_list[3]]);
+      const auto x1 = symmetrize_coordinates(origin, normal, m_xr[node_list[2]]);
+      const auto x2 = symmetrize_coordinates(origin, normal, m_xr[node_list[1]]);
+      const auto x3 = symmetrize_coordinates(origin, normal, m_xr[node_list[0]]);
+      const auto x4 = symmetrize_coordinates(origin, normal, m_xr[node_list[7]]);
+      const auto x5 = symmetrize_coordinates(origin, normal, m_xr[node_list[6]]);
+      const auto x6 = symmetrize_coordinates(origin, normal, m_xr[node_list[5]]);
+      const auto x7 = symmetrize_coordinates(origin, normal, m_xr[node_list[4]]);
+
+      const CubeTransformation T{x0, x1, x2, x3, x4, x5, x6, x7};
+
+      const auto& quadrature =
+        QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{m_polynomial_degree + 1});
+
+      _computeEjkMean(quadrature, T, Xj, Vi, mean_of_ejk);
+
+    } else {
+      throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+    }
+  }
+
+ 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];
+
+    const Rd& Xj = m_xj[cell_j_id];
+
+    _computeEjkMean(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];
+
+      _computeEjkMean(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];
+
+        _computeEjkMeanInSymmetricCell(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];
+        }
+      }
+    }
+  }
+
+  ElementIntegralReconstructionMatrixBuilder(const MeshType& mesh,
+                                             const size_t polynomial_degree,
+                                             const SmallArray<double>& inv_Vj_wq_detJ_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_basis_dimension{DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(
+        polynomial_degree)},
+      m_polynomial_degree{polynomial_degree},
+
+      m_inv_Vj_wq_detJ_ek{inv_Vj_wq_detJ_ek},
+      m_mean_j_of_ejk{mean_j_of_ejk},
+      m_mean_i_of_ejk{mean_i_of_ejk},
+
+      m_cell_to_node_matrix{mesh.connectivity().cellToNodeMatrix()},
+      m_stencil_array{stencil_array},
+      m_symmetry_origin_list{symmetry_origin_list},
+      m_symmetry_normal_list{symmetry_normal_list},
+      m_cell_type{mesh.connectivity().cellType()},
+      m_Vj{MeshDataManager::instance().getMeshData(mesh).Vj()},
+      m_xj{MeshDataManager::instance().getMeshData(mesh).xj()},
+      m_xr{mesh.xr()}
+  {}
+
+  ~ElementIntegralReconstructionMatrixBuilder() = default;
+};
+
+#endif   // ELEMENT_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP