diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 6073da9c8ca1547a7aa36479b098b43999191c6f..e29e8b8da42f6d897ed96398be3a2a0bd6c82691 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -828,6 +828,7 @@ PolynomialReconstruction::_build(
       return this->_build<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh,
                                                                                  discrete_function_variant_list);
     }
+    std::clog << "falling back to element integration in 3d\n";
     [[fallthrough]];
   }
   case IntegrationMethodType::element: {
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
index 2a481517f6aa70a8c4203e07445a7337e16f4912..0e3c01a7d4857ad39ff3be8555141be6b47c4caf 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
@@ -1,13 +1,53 @@
 #include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
 
 #include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureFormula.hpp>
 #include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
 #include <geometry/SymmetryUtils.hpp>
-#include <mesh/Connectivity.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshDataManager.hpp>
 #include <scheme/DiscreteFunctionDPk.hpp>
 
+template <>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::ConnectivityData
+{
+ private:
+  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;
+
+ public:
+  PUGS_INLINE
+  const auto&
+  cellToFaceMatrix() const noexcept
+  {
+    return m_cell_to_face_matrix;
+  }
+
+  PUGS_INLINE
+  const auto&
+  faceToNodeMatrix() const noexcept
+  {
+    return m_face_to_node_matrix;
+  }
+
+  PUGS_INLINE
+  const auto&
+  cellFaceIsReversed() const noexcept
+  {
+    return m_cell_face_is_reversed;
+  }
+
+  ConnectivityData(const Connectivity<2>& connectivity)
+    : m_cell_to_face_matrix{connectivity.cellToFaceMatrix()},
+      m_face_to_node_matrix{connectivity.faceToNodeMatrix()},
+      m_cell_face_is_reversed{connectivity.cellFaceIsReversed()}
+  {}
+
+  ~ConnectivityData() = default;
+};
+
 template <MeshConcept MeshTypeT>
 template <typename ConformTransformationT>
 void
@@ -63,14 +103,15 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
 
   const double inv_Vi = 1. / m_Vj[cell_id];
 
+  const auto& face_is_reversed    = m_dimensional_data->cellFaceIsReversed()[cell_id];
+  const auto& cell_face_list      = m_dimensional_data->cellToFaceMatrix()[cell_id];
+  const auto& face_to_node_matrix = m_dimensional_data->faceToNodeMatrix();
+
   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) {
+    auto face_node_list  = face_to_node_matrix[face_id];
+    if (face_is_reversed[i_face]) {
       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 {
@@ -94,18 +135,19 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
 
   const double inv_Vi = 1. / m_Vj[cell_id];
 
+  const auto& face_is_reversed    = m_dimensional_data->cellFaceIsReversed()[cell_id];
+  const auto& cell_face_list      = m_dimensional_data->cellToFaceMatrix()[cell_id];
+  const auto& face_to_node_matrix = m_dimensional_data->faceToNodeMatrix();
+
   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];
+    auto face_node_list  = 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) {
+    if (face_is_reversed[i_face]) {
       const LineTransformation<2> T{x1, x0};
       _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
     } else {
@@ -174,9 +216,9 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
     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_dimensional_data{std::make_shared<ConnectivityData>(mesh.connectivity())},
+
     m_stencil_array{stencil_array},
     m_symmetry_origin_list{symmetry_origin_list},
     m_symmetry_normal_list{symmetry_normal_list},
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
index df24d76a73097d13f6379681780d90aed6f49e1d..cb638fa2bef710f8ffe92257c546911b158f4b42 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
@@ -3,14 +3,14 @@
 
 #include <algebra/ShrinkMatrixView.hpp>
 #include <algebra/SmallMatrix.hpp>
-#include <analysis/QuadratureFormula.hpp>
-#include <geometry/LineTransformation.hpp>
 #include <mesh/ItemValue.hpp>
 #include <mesh/StencilArray.hpp>
-#include <mesh/SubItemValuePerItem.hpp>
 #include <scheme/PolynomialReconstruction.hpp>
 #include <utils/SmallArray.hpp>
 
+template <size_t Dimension>
+class QuadratureFormula;
+
 template <MeshConcept MeshTypeT>
 class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
 {
@@ -30,9 +30,8 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
   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;
+  class ConnectivityData;
+  std::shared_ptr<ConnectivityData> m_dimensional_data;
 
   const CellToCellStencilArray& m_stencil_array;
 
@@ -76,7 +75,6 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
                                               const SmallArray<const Rd>& symmetry_normal_list,
                                               const CellToCellStencilArray& stencil_array);
 
-  BoundaryIntegralReconstructionMatrixBuilder()  = default;
   ~BoundaryIntegralReconstructionMatrixBuilder() = default;
 };
 
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
index faededf698a291fbfb90addca1ac6fdcf15320a5..18306a8ecc635277f5e003d10faf53dfc41042a7 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
@@ -2,6 +2,8 @@
 
 #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>
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
index 1f4573438260e321dbd1461d2566f854236b25f4..7a3e8f7fc04d6d4ba21d4a36ed02ba24dd0d67d0 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
@@ -3,14 +3,15 @@
 
 #include <algebra/ShrinkMatrixView.hpp>
 #include <algebra/SmallMatrix.hpp>
-#include <analysis/QuadratureFormula.hpp>
-#include <analysis/QuadratureManager.hpp>
 #include <mesh/CellType.hpp>
 #include <mesh/ItemValue.hpp>
 #include <mesh/StencilArray.hpp>
 #include <scheme/PolynomialReconstruction.hpp>
 #include <utils/SmallArray.hpp>
 
+template <size_t Dimension>
+class QuadratureFormula;
+
 template <MeshConcept MeshTypeT>
 class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
 {