From dc5cccb476534edfc91193a297fb2b833023876f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 12 May 2025 15:06:51 +0200
Subject: [PATCH] Continue clean-up and prepare higher level threads data
 handling

---
 src/scheme/PolynomialReconstruction.cpp       | 998 +++++++++---------
 src/scheme/PolynomialReconstruction.hpp       |   3 +
 ...aryIntegralReconstructionMatrixBuilder.hpp |   6 +-
 ...entIntegralReconstructionMatrixBuilder.hpp |   6 +-
 src/scheme/test_reconstruction.cpp            |  26 +-
 5 files changed, 543 insertions(+), 496 deletions(-)

diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 3c257c2d..464ce1a0 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -43,6 +43,508 @@ symmetrize_coordinates(const TinyVector<Dimension>& origin,
   return u - 2 * dot(u - origin, normal) * normal;
 }
 
+template <MeshConcept MeshType>
+class PolynomialReconstruction::Internal
+{
+ private:
+  using Rd = TinyVector<MeshType::Dimension>;
+
+  friend PolynomialReconstruction;
+
+  template <typename MatrixType>
+  static void
+  buildB(const CellId& cell_j_id,
+         const CellToCellStencilArray& stencil_array,
+         const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
+         SmallArray<const Rd> symmetry_normal_list,
+         ShrinkMatrixView<MatrixType>& B)
+  {
+    auto stencil_cell_list = stencil_array[cell_j_id];
+
+    size_t column_begin = 0;
+    for (size_t i_discrete_function_variant = 0; i_discrete_function_variant < discrete_function_variant_list.size();
+         ++i_discrete_function_variant) {
+      const auto& discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
+
+      std::visit(
+        [&](auto&& discrete_function) {
+          using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
+          if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+            using DataType     = std::decay_t<typename DiscreteFunctionT::data_type>;
+            const DataType& qj = discrete_function[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 DataType& qi_qj  = discrete_function[cell_i_id] - qj;
+              if constexpr (std::is_arithmetic_v<DataType>) {
+                B(index, column_begin) = qi_qj;
+              } else if constexpr (is_tiny_vector_v<DataType>) {
+                for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
+                  B(index, kB) = qi_qj[k];
+                }
+              } else if constexpr (is_tiny_matrix_v<DataType>) {
+                for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                  const size_t kB = column_begin + p * DataType::NumberOfColumns;
+                  for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                    B(index, kB + q) = qi_qj(p, q);
+                  }
+                }
+              }
+            }
+
+            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];
+              for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+                const CellId cell_i_id = ghost_cell_list[i];
+                if constexpr (std::is_arithmetic_v<DataType>) {
+                  const DataType& qi_qj  = discrete_function[cell_i_id] - qj;
+                  B(index, column_begin) = qi_qj;
+                } else if constexpr (is_tiny_vector_v<DataType>) {
+                  if constexpr (DataType::Dimension == MeshType::Dimension) {
+                    const Rd& normal = symmetry_normal_list[i_symmetry];
+
+                    const DataType& qi    = discrete_function[cell_i_id];
+                    const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
+                    for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
+                      B(index, kB) = qi_qj[k];
+                    }
+                  } else {
+                    // LCOV_EXCL_START
+                    std::stringstream error_msg;
+                    error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+                              << " using a mesh of dimension " << MeshType::Dimension;
+                    throw UnexpectedError(error_msg.str());
+                    // LCOV_EXCL_STOP
+                  }
+                } else if constexpr (is_tiny_matrix_v<DataType>) {
+                  if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
+                                (DataType::NumberOfColumns == MeshType::Dimension)) {
+                    const Rd& normal = symmetry_normal_list[i_symmetry];
+
+                    const DataType& qi    = discrete_function[cell_i_id];
+                    const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
+                    for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                      for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                        B(index, column_begin + p * DataType::NumberOfColumns + q) = qi_qj(p, q);
+                      }
+                    }
+                  } else {
+                    // LCOV_EXCL_START
+                    std::stringstream error_msg;
+                    error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
+                              << DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
+                    throw UnexpectedError(error_msg.str());
+                    // LCOV_EXCL_STOP
+                  }
+                }
+              }
+            }
+
+            if constexpr (std::is_arithmetic_v<DataType>) {
+              ++column_begin;
+            } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
+              column_begin += DataType::Dimension;
+            }
+          } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+            using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+
+            const auto qj_vector = discrete_function[cell_j_id];
+
+            if constexpr (std::is_arithmetic_v<DataType>) {
+              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];
+                for (size_t l = 0; l < qj_vector.size(); ++l) {
+                  const DataType& qj         = qj_vector[l];
+                  const DataType& qi_qj      = discrete_function[cell_i_id][l] - qj;
+                  B(index, column_begin + l) = qi_qj;
+                }
+              }
+
+              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];
+                for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+                  const CellId cell_i_id = ghost_cell_list[i];
+                  for (size_t l = 0; l < qj_vector.size(); ++l) {
+                    const DataType& qj         = qj_vector[l];
+                    const DataType& qi_qj      = discrete_function[cell_i_id][l] - qj;
+                    B(index, column_begin + l) = qi_qj;
+                  }
+                }
+              }
+            } else if constexpr (is_tiny_vector_v<DataType>) {
+              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];
+                for (size_t l = 0; l < qj_vector.size(); ++l) {
+                  const DataType& qj    = qj_vector[l];
+                  const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
+                  for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension; ++k, ++kB) {
+                    B(index, kB) = qi_qj[k];
+                  }
+                }
+              }
+
+              for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                   ++i_symmetry) {
+                if constexpr (DataType::Dimension == MeshType::Dimension) {
+                  auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+                  auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+                  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];
+
+                    for (size_t l = 0; l < qj_vector.size(); ++l) {
+                      const DataType& qj    = qj_vector[l];
+                      const DataType& qi    = discrete_function[cell_i_id][l];
+                      const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
+                      for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
+                           ++k, ++kB) {
+                        B(index, kB) = qi_qj[k];
+                      }
+                    }
+                  }
+                } else {
+                  // LCOV_EXCL_START
+                  std::stringstream error_msg;
+                  error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+                            << " using a mesh of dimension " << MeshType::Dimension;
+                  throw UnexpectedError(error_msg.str());
+                  // LCOV_EXCL_STOP
+                }
+              }
+            } else if constexpr (is_tiny_matrix_v<DataType>) {
+              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];
+                for (size_t l = 0; l < qj_vector.size(); ++l) {
+                  const DataType& qj    = qj_vector[l];
+                  const DataType& qi    = discrete_function[cell_i_id][l];
+                  const DataType& qi_qj = qi - qj;
+
+                  for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                    const size_t kB = column_begin + l * DataType::Dimension + p * DataType::NumberOfColumns;
+                    for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                      B(index, kB + q) = qi_qj(p, q);
+                    }
+                  }
+                }
+              }
+
+              for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
+                   ++i_symmetry) {
+                if constexpr ((DataType::NumberOfRows == MeshType::Dimension) and
+                              (DataType::NumberOfColumns == MeshType::Dimension)) {
+                  auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+                  auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+                  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];
+
+                    for (size_t l = 0; l < qj_vector.size(); ++l) {
+                      const DataType& qj    = qj_vector[l];
+                      const DataType& qi    = discrete_function[cell_i_id][l];
+                      const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
+
+                      for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                        const size_t kB = column_begin + l * DataType::Dimension + p * DataType::NumberOfColumns;
+                        for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                          B(index, kB + q) = qi_qj(p, q);
+                        }
+                      }
+                    }
+                  }
+                } else {
+                  // LCOV_EXCL_START
+                  std::stringstream error_msg;
+                  error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+                            << " using a mesh of dimension " << MeshType::Dimension;
+                  throw UnexpectedError(error_msg.str());
+                  // LCOV_EXCL_STOP
+                }
+              }
+            }
+
+            if constexpr (std::is_arithmetic_v<DataType>) {
+              column_begin += qj_vector.size();
+            } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
+              column_begin += qj_vector.size() * DataType::Dimension;
+            }
+
+          } else {
+            // LCOV_EXCL_START
+            throw UnexpectedError("invalid discrete function type");
+            // LCOV_EXCL_STOP
+          }
+        },
+        discrete_function_variant->discreteFunction());
+    }
+  }
+
+  template <typename MatrixType>
+  static void
+  rowWeighting(const CellId& cell_j_id,
+               const CellToCellStencilArray& stencil_array,
+               const CellValue<const Rd>& xj,
+               const SmallArray<const Rd>& symmetry_origin_list,
+               const SmallArray<const Rd>& symmetry_normal_list,
+               ShrinkMatrixView<MatrixType>& A,
+               ShrinkMatrixView<MatrixType>& B)
+  {
+    // Add row weighting (give more importance to cells that are
+    // closer to j)
+    auto stencil_cell_list = stencil_array[cell_j_id];
+
+    const Rd& Xj = 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 double weight    = 1. / l2Norm(xj[cell_i_id] - Xj);
+      for (size_t l = 0; l < A.numberOfColumns(); ++l) {
+        A(index, l) *= weight;
+      }
+      for (size_t l = 0; l < B.numberOfColumns(); ++l) {
+        B(index, l) *= weight;
+      }
+    }
+    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];
+        const double weight    = 1. / l2Norm(symmetrize_coordinates(origin, normal, xj[cell_i_id]) - Xj);
+        for (size_t l = 0; l < A.numberOfColumns(); ++l) {
+          A(index, l) *= weight;
+        }
+        for (size_t l = 0; l < B.numberOfColumns(); ++l) {
+          B(index, l) *= weight;
+        }
+      }
+    }
+  }
+
+  template <typename MatrixType>
+  static void
+  solveCollectionInPlaceWithPreconditionner(const ShrinkMatrixView<MatrixType>& A,
+                                            const SmallMatrix<double>& X,
+                                            const ShrinkMatrixView<MatrixType>& B,
+                                            const SmallVector<double>& G)
+  {
+    for (size_t l = 0; l < A.numberOfColumns(); ++l) {
+      double g = 0;
+      for (size_t i = 0; i < A.numberOfRows(); ++i) {
+        const double Ail = A(i, l);
+
+        g += Ail * Ail;
+      }
+      G[l] = std::sqrt(g);
+    }
+
+    for (size_t l = 0; l < A.numberOfColumns(); ++l) {
+      const double Gl = G[l];
+      for (size_t i = 0; i < A.numberOfRows(); ++i) {
+        A(i, l) *= Gl;
+      }
+    }
+
+    Givens::solveCollectionInPlace(A, X, B);
+
+    for (size_t l = 0; l < X.numberOfRows(); ++l) {
+      const double Gl = G[l];
+      for (size_t i = 0; i < X.numberOfColumns(); ++i) {
+        X(l, i) *= Gl;
+      }
+    }
+  }
+
+  static void
+  populateDiscreteFunctionDPkByCell(
+    const CellId& cell_j_id,
+    const size_t& degree,
+    const SmallMatrix<double>& X,
+    const SmallArray<double>& mean_j_of_ejk,
+    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
+    const std::vector<PolynomialReconstruction::MutableDiscreteFunctionDPkVariant>&
+      mutable_discrete_function_dpk_variant_list)
+  {
+    size_t column_begin = 0;
+    for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size(); ++i_dpk_variant) {
+      const auto& dpk_variant = mutable_discrete_function_dpk_variant_list[i_dpk_variant];
+
+      const auto& discrete_function_variant = discrete_function_variant_list[i_dpk_variant];
+
+      std::visit(
+        [&](auto&& dpk_function, auto&& p0_function) {
+          using DPkFunctionT = std::decay_t<decltype(dpk_function)>;
+          using P0FunctionT  = std::decay_t<decltype(p0_function)>;
+          using DataType     = std::remove_const_t<std::decay_t<typename DPkFunctionT::data_type>>;
+          using P0DataType   = std::remove_const_t<std::decay_t<typename P0FunctionT::data_type>>;
+
+          if constexpr (std::is_same_v<DataType, P0DataType>) {
+            if constexpr (is_discrete_function_P0_v<P0FunctionT>) {
+              if constexpr (is_discrete_function_dpk_scalar_v<DPkFunctionT>) {
+                auto dpk_j = dpk_function.coefficients(cell_j_id);
+                dpk_j[0]   = p0_function[cell_j_id];
+
+                if constexpr (std::is_arithmetic_v<DataType>) {
+                  if (degree > 1) {
+                    for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                      dpk_j[0] -= X(i, column_begin) * mean_j_of_ejk[i];
+                    }
+                  }
+
+                  for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                    auto& dpk_j_ip1 = dpk_j[i + 1];
+                    dpk_j_ip1       = X(i, column_begin);
+                  }
+                  ++column_begin;
+                } else if constexpr (is_tiny_vector_v<DataType>) {
+                  if (degree > 1) {
+                    for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                      auto& dpk_j_0 = dpk_j[0];
+                      for (size_t k = 0; k < DataType::Dimension; ++k) {
+                        dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
+                      }
+                    }
+                  }
+
+                  for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                    auto& dpk_j_ip1 = dpk_j[i + 1];
+                    for (size_t k = 0; k < DataType::Dimension; ++k) {
+                      dpk_j_ip1[k] = X(i, column_begin + k);
+                    }
+                  }
+                  column_begin += DataType::Dimension;
+                } else if constexpr (is_tiny_matrix_v<DataType>) {
+                  if (degree > 1) {
+                    for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                      auto& dpk_j_0 = dpk_j[0];
+                      for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+                        for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+                          dpk_j_0(k, l) -= X(i, column_begin + k * DataType::NumberOfColumns + l) * mean_j_of_ejk[i];
+                        }
+                      }
+                    }
+                  }
+
+                  for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                    auto& dpk_j_ip1 = dpk_j[i + 1];
+                    for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+                      for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+                        dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
+                      }
+                    }
+                  }
+                  column_begin += DataType::Dimension;
+                } else {
+                  // LCOV_EXCL_START
+                  throw UnexpectedError("unexpected data type");
+                  // LCOV_EXCL_STOP
+                }
+              } else {
+                // LCOV_EXCL_START
+                throw UnexpectedError("unexpected discrete dpk function type");
+                // LCOV_EXCL_STOP
+              }
+            } else if constexpr (is_discrete_function_P0_vector_v<P0FunctionT>) {
+              if constexpr (is_discrete_function_dpk_vector_v<DPkFunctionT>) {
+                auto dpk_j        = dpk_function.coefficients(cell_j_id);
+                auto cell_vector  = p0_function[cell_j_id];
+                const size_t size = X.numberOfRows() + 1;
+
+                for (size_t l = 0; l < cell_vector.size(); ++l) {
+                  const size_t component_begin = l * size;
+                  dpk_j[component_begin]       = cell_vector[l];
+                  if constexpr (std::is_arithmetic_v<DataType>) {
+                    if (degree > 1) {
+                      for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                        dpk_j[component_begin] -= X(i, column_begin) * mean_j_of_ejk[i];
+                      }
+                    }
+
+                    for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                      auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
+                      dpk_j_ip1       = X(i, column_begin);
+                    }
+                    ++column_begin;
+                  } else if constexpr (is_tiny_vector_v<DataType>) {
+                    if (degree > 1) {
+                      for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                        auto& dpk_j_0 = dpk_j[component_begin];
+                        for (size_t k = 0; k < DataType::Dimension; ++k) {
+                          dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
+                        }
+                      }
+                    }
+
+                    for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                      auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
+                      for (size_t k = 0; k < DataType::Dimension; ++k) {
+                        dpk_j_ip1[k] = X(i, column_begin + k);
+                      }
+                    }
+                    column_begin += DataType::Dimension;
+                  } else if constexpr (is_tiny_matrix_v<DataType>) {
+                    if (degree > 1) {
+                      for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                        auto& dpk_j_0 = dpk_j[component_begin];
+                        for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                          for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                            dpk_j_0(p, q) -= X(i, column_begin + p * DataType::NumberOfColumns + q) * mean_j_of_ejk[i];
+                          }
+                        }
+                      }
+                    }
+
+                    for (size_t i = 0; i < X.numberOfRows(); ++i) {
+                      auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
+                      for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                        for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                          dpk_j_ip1(p, q) = X(i, column_begin + p * DataType::NumberOfColumns + q);
+                        }
+                      }
+                    }
+                    column_begin += DataType::Dimension;
+                  } else {
+                    // LCOV_EXCL_START
+                    throw UnexpectedError("unexpected data type");
+                    // LCOV_EXCL_STOP
+                  }
+                }
+              } else {
+                // LCOV_EXCL_START
+                throw UnexpectedError("unexpected discrete dpk function type");
+                // LCOV_EXCL_STOP
+              }
+            } else {
+              // LCOV_EXCL_START
+              throw UnexpectedError("unexpected discrete function type");
+              // LCOV_EXCL_STOP
+            }
+          } else {
+            // LCOV_EXCL_START
+            throw UnexpectedError("incompatible data types");
+            // LCOV_EXCL_STOP
+          }
+        },
+        dpk_variant.mutableDiscreteFunctionDPk(), discrete_function_variant->discreteFunction());
+    }
+  }
+};
+
 size_t
 PolynomialReconstruction::_getNumberOfColumns(
   const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
@@ -246,11 +748,7 @@ PolynomialReconstruction::_build(
   SmallArray<SmallVector<double>> G_pool(Kokkos::DefaultExecutionSpace::concurrency());
   SmallArray<SmallMatrix<double>> X_pool(Kokkos::DefaultExecutionSpace::concurrency());
 
-  SmallArray<SmallArray<double>> inv_Vj_wq_detJ_ek_pool(Kokkos::DefaultExecutionSpace::concurrency());
   SmallArray<SmallArray<double>> mean_j_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
-  SmallArray<SmallArray<double>> mean_i_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
-
-  SmallArray<SmallArray<double>> inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool(Kokkos::DefaultExecutionSpace::concurrency());
 
   for (size_t i = 0; i < A_pool.size(); ++i) {
     A_pool[i] = SmallMatrix<double>(max_stencil_size, basis_dimension - 1);
@@ -258,11 +756,7 @@ PolynomialReconstruction::_build(
     G_pool[i] = SmallVector<double>(basis_dimension - 1);
     X_pool[i] = SmallMatrix<double>(basis_dimension - 1, number_of_columns);
 
-    inv_Vj_wq_detJ_ek_pool[i] = SmallArray<double>(basis_dimension);
-    mean_j_of_ejk_pool[i]     = SmallArray<double>(basis_dimension - 1);
-    mean_i_of_ejk_pool[i]     = SmallArray<double>(basis_dimension - 1);
-
-    inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[i] = SmallArray<double>(basis_dimension);
+    mean_j_of_ejk_pool[i] = SmallArray<double>(basis_dimension - 1);
   }
 
   std::shared_ptr p_cell_center_reconstruction_matrix_builder =
@@ -273,14 +767,11 @@ PolynomialReconstruction::_build(
     element_integral_reconstruction_matrix_builder_pool{A_pool.size()};
 
   for (size_t t = 0; t < element_integral_reconstruction_matrix_builder_pool.size(); ++t) {
-    SmallArray<double>& inv_Vj_wq_detJ_ek = inv_Vj_wq_detJ_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];
+    SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
 
     element_integral_reconstruction_matrix_builder_pool[t] =
       std::make_shared<ElementIntegralReconstructionMatrixBuilder<MeshType>>(*p_mesh, m_descriptor.degree(),
-                                                                             inv_Vj_wq_detJ_ek, mean_j_of_ejk,
-                                                                             mean_i_of_ejk, symmetry_origin_list,
+                                                                             mean_j_of_ejk, symmetry_origin_list,
                                                                              symmetry_normal_list, stencil_array);
   }
 
@@ -289,15 +780,11 @@ PolynomialReconstruction::_build(
 
   if constexpr (MeshType::Dimension == 2) {
     for (size_t t = 0; t < boundary_integral_reconstruction_matrix_builder_pool.size(); ++t) {
-      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];
+      SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
 
       boundary_integral_reconstruction_matrix_builder_pool[t] =
         std::make_shared<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(*p_mesh, m_descriptor.degree(),
-                                                                                inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
-                                                                                mean_j_of_ejk, mean_i_of_ejk,
-                                                                                symmetry_origin_list,
+                                                                                mean_j_of_ejk, symmetry_origin_list,
                                                                                 symmetry_normal_list, stencil_array);
     }
   }
@@ -307,239 +794,10 @@ PolynomialReconstruction::_build(
       if (cell_is_owned[cell_j_id]) {
         const int32_t t = tokens.acquire();
 
-        auto stencil_cell_list = stencil_array[cell_j_id];
-
+        ShrinkMatrixView A(A_pool[t], full_stencil_size(cell_j_id));
         ShrinkMatrixView B(B_pool[t], full_stencil_size(cell_j_id));
 
-        size_t column_begin = 0;
-        for (size_t i_discrete_function_variant = 0;
-             i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
-          const auto& discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
-
-          std::visit(
-            [&](auto&& discrete_function) {
-              using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
-              if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
-                using DataType     = std::decay_t<typename DiscreteFunctionT::data_type>;
-                const DataType& qj = discrete_function[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 DataType& qi_qj  = discrete_function[cell_i_id] - qj;
-                  if constexpr (std::is_arithmetic_v<DataType>) {
-                    B(index, column_begin) = qi_qj;
-                  } else if constexpr (is_tiny_vector_v<DataType>) {
-                    for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
-                      B(index, kB) = qi_qj[k];
-                    }
-                  } else if constexpr (is_tiny_matrix_v<DataType>) {
-                    for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
-                      const size_t kB = column_begin + p * DataType::NumberOfColumns;
-                      for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
-                        B(index, kB + q) = qi_qj(p, q);
-                      }
-                    }
-                  }
-                }
-
-                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];
-                  for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
-                    const CellId cell_i_id = ghost_cell_list[i];
-                    if constexpr (std::is_arithmetic_v<DataType>) {
-                      const DataType& qi_qj  = discrete_function[cell_i_id] - qj;
-                      B(index, column_begin) = qi_qj;
-                    } else if constexpr (is_tiny_vector_v<DataType>) {
-                      if constexpr (DataType::Dimension == MeshType::Dimension) {
-                        const Rd& normal = symmetry_normal_list[i_symmetry];
-
-                        const DataType& qi    = discrete_function[cell_i_id];
-                        const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
-                        for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
-                          B(index, kB) = qi_qj[k];
-                        }
-                      } else {
-                        // LCOV_EXCL_START
-                        std::stringstream error_msg;
-                        error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
-                                  << " using a mesh of dimension " << MeshType::Dimension;
-                        throw UnexpectedError(error_msg.str());
-                        // LCOV_EXCL_STOP
-                      }
-                    } else if constexpr (is_tiny_matrix_v<DataType>) {
-                      if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
-                                    (DataType::NumberOfColumns == MeshType::Dimension)) {
-                        const Rd& normal = symmetry_normal_list[i_symmetry];
-
-                        const DataType& qi    = discrete_function[cell_i_id];
-                        const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
-                        for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
-                          for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
-                            B(index, column_begin + p * DataType::NumberOfColumns + q) = qi_qj(p, q);
-                          }
-                        }
-                      } else {
-                        // LCOV_EXCL_START
-                        std::stringstream error_msg;
-                        error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
-                                  << DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
-                        throw UnexpectedError(error_msg.str());
-                        // LCOV_EXCL_STOP
-                      }
-                    }
-                  }
-                }
-
-                if constexpr (std::is_arithmetic_v<DataType>) {
-                  ++column_begin;
-                } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
-                  column_begin += DataType::Dimension;
-                }
-              } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
-                using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
-
-                const auto qj_vector = discrete_function[cell_j_id];
-
-                if constexpr (std::is_arithmetic_v<DataType>) {
-                  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];
-                    for (size_t l = 0; l < qj_vector.size(); ++l) {
-                      const DataType& qj         = qj_vector[l];
-                      const DataType& qi_qj      = discrete_function[cell_i_id][l] - qj;
-                      B(index, column_begin + l) = qi_qj;
-                    }
-                  }
-
-                  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];
-                    for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
-                      const CellId cell_i_id = ghost_cell_list[i];
-                      for (size_t l = 0; l < qj_vector.size(); ++l) {
-                        const DataType& qj         = qj_vector[l];
-                        const DataType& qi_qj      = discrete_function[cell_i_id][l] - qj;
-                        B(index, column_begin + l) = qi_qj;
-                      }
-                    }
-                  }
-                } else if constexpr (is_tiny_vector_v<DataType>) {
-                  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];
-                    for (size_t l = 0; l < qj_vector.size(); ++l) {
-                      const DataType& qj    = qj_vector[l];
-                      const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
-                      for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
-                           ++k, ++kB) {
-                        B(index, kB) = qi_qj[k];
-                      }
-                    }
-                  }
-
-                  for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
-                       ++i_symmetry) {
-                    if constexpr (DataType::Dimension == MeshType::Dimension) {
-                      auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
-                      auto ghost_cell_list = ghost_stencil[cell_j_id];
-
-                      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];
-
-                        for (size_t l = 0; l < qj_vector.size(); ++l) {
-                          const DataType& qj    = qj_vector[l];
-                          const DataType& qi    = discrete_function[cell_i_id][l];
-                          const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
-                          for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
-                               ++k, ++kB) {
-                            B(index, kB) = qi_qj[k];
-                          }
-                        }
-                      }
-                    } else {
-                      // LCOV_EXCL_START
-                      std::stringstream error_msg;
-                      error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
-                                << " using a mesh of dimension " << MeshType::Dimension;
-                      throw UnexpectedError(error_msg.str());
-                      // LCOV_EXCL_STOP
-                    }
-                  }
-                } else if constexpr (is_tiny_matrix_v<DataType>) {
-                  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];
-                    for (size_t l = 0; l < qj_vector.size(); ++l) {
-                      const DataType& qj    = qj_vector[l];
-                      const DataType& qi    = discrete_function[cell_i_id][l];
-                      const DataType& qi_qj = qi - qj;
-
-                      for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
-                        const size_t kB = column_begin + l * DataType::Dimension + p * DataType::NumberOfColumns;
-                        for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
-                          B(index, kB + q) = qi_qj(p, q);
-                        }
-                      }
-                    }
-                  }
-
-                  for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
-                       ++i_symmetry) {
-                    if constexpr ((DataType::NumberOfRows == MeshType::Dimension) and
-                                  (DataType::NumberOfColumns == MeshType::Dimension)) {
-                      auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
-                      auto ghost_cell_list = ghost_stencil[cell_j_id];
-
-                      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];
-
-                        for (size_t l = 0; l < qj_vector.size(); ++l) {
-                          const DataType& qj    = qj_vector[l];
-                          const DataType& qi    = discrete_function[cell_i_id][l];
-                          const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
-
-                          for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
-                            const size_t kB = column_begin + l * DataType::Dimension + p * DataType::NumberOfColumns;
-                            for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
-                              B(index, kB + q) = qi_qj(p, q);
-                            }
-                          }
-                        }
-                      }
-                    } else {
-                      // LCOV_EXCL_START
-                      std::stringstream error_msg;
-                      error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
-                                << " using a mesh of dimension " << MeshType::Dimension;
-                      throw UnexpectedError(error_msg.str());
-                      // LCOV_EXCL_STOP
-                    }
-                  }
-                }
-
-                if constexpr (std::is_arithmetic_v<DataType>) {
-                  column_begin += qj_vector.size();
-                } else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
-                  column_begin += qj_vector.size() * DataType::Dimension;
-                }
-
-              } else {
-                // LCOV_EXCL_START
-                throw UnexpectedError("invalid discrete function type");
-                // LCOV_EXCL_STOP
-              }
-            },
-            discrete_function_variant->discreteFunction());
-        }
-
-        ShrinkMatrixView A(A_pool[t], full_stencil_size(cell_j_id));
+        Internal<MeshType>::buildB(cell_j_id, stencil_array, discrete_function_variant_list, symmetry_normal_list, B);
 
         if ((m_descriptor.degree() == 1) and
             (m_descriptor.integrationMethodType() == IntegrationMethodType::cell_center)) {
@@ -562,248 +820,24 @@ PolynomialReconstruction::_build(
         }
 
         if (m_descriptor.rowWeighting()) {
-          // Add row weighting (give more importance to cells that are
-          // closer to j)
-          const Rd& Xj = 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 double weight    = 1. / l2Norm(xj[cell_i_id] - Xj);
-            for (size_t l = 0; l < basis_dimension - 1; ++l) {
-              A(index, l) *= weight;
-            }
-            for (size_t l = 0; l < number_of_columns; ++l) {
-              B(index, l) *= weight;
-            }
-          }
-          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];
-              const double weight    = 1. / l2Norm(symmetrize_coordinates(origin, normal, xj[cell_i_id]) - Xj);
-              for (size_t l = 0; l < basis_dimension - 1; ++l) {
-                A(index, l) *= weight;
-              }
-              for (size_t l = 0; l < number_of_columns; ++l) {
-                B(index, l) *= weight;
-              }
-            }
-          }
+          Internal<MeshType>::rowWeighting(cell_j_id, stencil_array, xj, symmetry_origin_list, symmetry_normal_list, A,
+                                           B);
         }
 
         const SmallMatrix<double>& X = X_pool[t];
 
         if (m_descriptor.preconditioning()) {
-          // Add column  weighting preconditioning (increase the presition)
+          // Add column  weighting preconditioning (increase the precision)
           SmallVector<double>& G = G_pool[t];
 
-          for (size_t l = 0; l < A.numberOfColumns(); ++l) {
-            double g = 0;
-            for (size_t i = 0; i < A.numberOfRows(); ++i) {
-              const double Ail = A(i, l);
-
-              g += Ail * Ail;
-            }
-            G[l] = std::sqrt(g);
-          }
-
-          for (size_t l = 0; l < A.numberOfColumns(); ++l) {
-            const double Gl = G[l];
-            for (size_t i = 0; i < A.numberOfRows(); ++i) {
-              A(i, l) *= Gl;
-            }
-          }
-
-          Givens::solveCollectionInPlace(A, X, B);
-
-          for (size_t l = 0; l < X.numberOfRows(); ++l) {
-            const double Gl = G[l];
-            for (size_t i = 0; i < X.numberOfColumns(); ++i) {
-              X(l, i) *= Gl;
-            }
-          }
+          Internal<MeshType>::solveCollectionInPlaceWithPreconditionner(A, X, B, G);
         } else {
           Givens::solveCollectionInPlace(A, X, B);
         }
 
-        column_begin = 0;
-        for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size();
-             ++i_dpk_variant) {
-          const auto& dpk_variant = mutable_discrete_function_dpk_variant_list[i_dpk_variant];
-
-          const auto& discrete_function_variant = discrete_function_variant_list[i_dpk_variant];
-
-          std::visit(
-            [&](auto&& dpk_function, auto&& p0_function) {
-              using DPkFunctionT = std::decay_t<decltype(dpk_function)>;
-              using P0FunctionT  = std::decay_t<decltype(p0_function)>;
-              using DataType     = std::remove_const_t<std::decay_t<typename DPkFunctionT::data_type>>;
-              using P0DataType   = std::remove_const_t<std::decay_t<typename P0FunctionT::data_type>>;
-
-              if constexpr (std::is_same_v<DataType, P0DataType>) {
-                if constexpr (is_discrete_function_P0_v<P0FunctionT>) {
-                  if constexpr (is_discrete_function_dpk_scalar_v<DPkFunctionT>) {
-                    auto dpk_j = dpk_function.coefficients(cell_j_id);
-                    dpk_j[0]   = p0_function[cell_j_id];
-
-                    if constexpr (std::is_arithmetic_v<DataType>) {
-                      if (m_descriptor.degree() > 1) {
-                        auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
-                        for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                          dpk_j[0] -= X(i, column_begin) * mean_j_of_ejk[i];
-                        }
-                      }
-
-                      for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                        auto& dpk_j_ip1 = dpk_j[i + 1];
-                        dpk_j_ip1       = X(i, column_begin);
-                      }
-                      ++column_begin;
-                    } else if constexpr (is_tiny_vector_v<DataType>) {
-                      if (m_descriptor.degree() > 1) {
-                        auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
-                        for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                          auto& dpk_j_0 = dpk_j[0];
-                          for (size_t k = 0; k < DataType::Dimension; ++k) {
-                            dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
-                          }
-                        }
-                      }
-
-                      for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                        auto& dpk_j_ip1 = dpk_j[i + 1];
-                        for (size_t k = 0; k < DataType::Dimension; ++k) {
-                          dpk_j_ip1[k] = X(i, column_begin + k);
-                        }
-                      }
-                      column_begin += DataType::Dimension;
-                    } else if constexpr (is_tiny_matrix_v<DataType>) {
-                      if (m_descriptor.degree() > 1) {
-                        auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
-                        for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                          auto& dpk_j_0 = dpk_j[0];
-                          for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
-                            for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
-                              dpk_j_0(k, l) -=
-                                X(i, column_begin + k * DataType::NumberOfColumns + l) * mean_j_of_ejk[i];
-                            }
-                          }
-                        }
-                      }
-
-                      for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                        auto& dpk_j_ip1 = dpk_j[i + 1];
-                        for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
-                          for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
-                            dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
-                          }
-                        }
-                      }
-                      column_begin += DataType::Dimension;
-                    } else {
-                      // LCOV_EXCL_START
-                      throw UnexpectedError("unexpected data type");
-                      // LCOV_EXCL_STOP
-                    }
-                  } else {
-                    // LCOV_EXCL_START
-                    throw UnexpectedError("unexpected discrete dpk function type");
-                    // LCOV_EXCL_STOP
-                  }
-                } else if constexpr (is_discrete_function_P0_vector_v<P0FunctionT>) {
-                  if constexpr (is_discrete_function_dpk_vector_v<DPkFunctionT>) {
-                    auto dpk_j        = dpk_function.coefficients(cell_j_id);
-                    auto cell_vector  = p0_function[cell_j_id];
-                    const size_t size = basis_dimension;
-
-                    for (size_t l = 0; l < cell_vector.size(); ++l) {
-                      const size_t component_begin = l * size;
-                      dpk_j[component_begin]       = cell_vector[l];
-                      if constexpr (std::is_arithmetic_v<DataType>) {
-                        if (m_descriptor.degree() > 1) {
-                          auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
-                          for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                            dpk_j[component_begin] -= X(i, column_begin) * mean_j_of_ejk[i];
-                          }
-                        }
-
-                        for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                          auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
-                          dpk_j_ip1       = X(i, column_begin);
-                        }
-                        ++column_begin;
-                      } else if constexpr (is_tiny_vector_v<DataType>) {
-                        if (m_descriptor.degree() > 1) {
-                          auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
-                          for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                            auto& dpk_j_0 = dpk_j[component_begin];
-                            for (size_t k = 0; k < DataType::Dimension; ++k) {
-                              dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
-                            }
-                          }
-                        }
-
-                        for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                          auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
-                          for (size_t k = 0; k < DataType::Dimension; ++k) {
-                            dpk_j_ip1[k] = X(i, column_begin + k);
-                          }
-                        }
-                        column_begin += DataType::Dimension;
-                      } else if constexpr (is_tiny_matrix_v<DataType>) {
-                        if (m_descriptor.degree() > 1) {
-                          auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
-                          for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                            auto& dpk_j_0 = dpk_j[component_begin];
-                            for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
-                              for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
-                                dpk_j_0(p, q) -=
-                                  X(i, column_begin + p * DataType::NumberOfColumns + q) * mean_j_of_ejk[i];
-                              }
-                            }
-                          }
-                        }
-
-                        for (size_t i = 0; i < basis_dimension - 1; ++i) {
-                          auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
-                          for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
-                            for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
-                              dpk_j_ip1(p, q) = X(i, column_begin + p * DataType::NumberOfColumns + q);
-                            }
-                          }
-                        }
-                        column_begin += DataType::Dimension;
-                      } else {
-                        // LCOV_EXCL_START
-                        throw UnexpectedError("unexpected data type");
-                        // LCOV_EXCL_STOP
-                      }
-                    }
-                  } else {
-                    // LCOV_EXCL_START
-                    throw UnexpectedError("unexpected discrete dpk function type");
-                    // LCOV_EXCL_STOP
-                  }
-                } else {
-                  // LCOV_EXCL_START
-                  throw UnexpectedError("unexpected discrete function type");
-                  // LCOV_EXCL_STOP
-                }
-              } else {
-                // LCOV_EXCL_START
-                throw UnexpectedError("incompatible data types");
-                // LCOV_EXCL_STOP
-              }
-            },
-            dpk_variant.mutableDiscreteFunctionDPk(), discrete_function_variant->discreteFunction());
-        }
+        Internal<MeshType>::populateDiscreteFunctionDPkByCell(cell_j_id, m_descriptor.degree(), X,
+                                                              mean_j_of_ejk_pool[t], discrete_function_variant_list,
+                                                              mutable_discrete_function_dpk_variant_list);
 
         tokens.release(t);
       }
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index 349ae457..fabd3e8e 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -10,6 +10,9 @@ class DiscreteFunctionVariant;
 class PolynomialReconstruction
 {
  private:
+  template <MeshConcept MeshType>
+  class Internal;
+
   class MutableDiscreteFunctionDPkVariant;
 
   template <MeshConcept MeshType>
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
index 350618f5..d22ecbb1 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
@@ -188,9 +188,7 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
 
   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)
@@ -198,9 +196,9 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
       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_inv_Vj_alpha_p_1_wq_X_prime_orth_ek{m_basis_dimension},
       m_mean_j_of_ejk{mean_j_of_ejk},
-      m_mean_i_of_ejk{mean_i_of_ejk},
+      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()},
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
index efb92d0d..1e5b7d71 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
@@ -502,9 +502,7 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
 
   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)
@@ -512,9 +510,9 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
         polynomial_degree)},
       m_polynomial_degree{polynomial_degree},
 
-      m_wq_detJ_ek{inv_Vj_wq_detJ_ek},
+      m_wq_detJ_ek{m_basis_dimension},
       m_mean_j_of_ejk{mean_j_of_ejk},
-      m_mean_i_of_ejk{mean_i_of_ejk},
+      m_mean_i_of_ejk{m_basis_dimension - 1},
 
       m_cell_to_node_matrix{mesh.connectivity().cellToNodeMatrix()},
       m_stencil_array{stencil_array},
diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
index ec106f49..c0bfa81f 100644
--- a/src/scheme/test_reconstruction.cpp
+++ b/src/scheme/test_reconstruction.cpp
@@ -6,14 +6,28 @@
 void
 test_reconstruction(std::shared_ptr<const DiscreteFunctionVariant> df, size_t degree, size_t number)
 {
-  PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::boundary, degree};
-  PolynomialReconstruction rec_builder{descriptor};
+  {
+    PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::boundary, degree};
+    PolynomialReconstruction rec_builder{descriptor};
 
-  Timer t;
+    Timer t;
 
-  for (size_t i = 0; i < number; ++i) {
-    auto rec = rec_builder.build(df);
+    for (size_t i = 0; i < number; ++i) {
+      auto rec = rec_builder.build(df);
+    }
+
+    std::cout << "*** Elapsed time for " << number << " boundary reconstructions: " << t.seconds() << "s ***\n";
   }
+  {
+    PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree};
+    PolynomialReconstruction rec_builder{descriptor};
+
+    Timer t;
 
-  std::cout << "*** Elapsed time for " << number << " reconstructions: " << t.seconds() << "s ***\n";
+    for (size_t i = 0; i < number; ++i) {
+      auto rec = rec_builder.build(df);
+    }
+
+    std::cout << "*** Elapsed time for " << number << " element reconstructions: " << t.seconds() << "s ***\n";
+  }
 }
-- 
GitLab