Skip to content
Snippets Groups Projects
Select Git revision
  • 82d750676e323b9fa68ee460b04188ab7142979f
  • develop default protected
  • origin/stage/bouguettaia
  • feature/kinetic-schemes
  • feature/reconstruction
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/variational-hydro
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • feature/escobar-smoother
  • feature/hypoelasticity-clean
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

ParallelChecker.hpp

Blame
  • PolynomialReconstruction.cpp 39.95 KiB
    #include <scheme/PolynomialReconstruction.hpp>
    
    #include <algebra/Givens.hpp>
    #include <algebra/ShrinkMatrixView.hpp>
    #include <algebra/ShrinkVectorView.hpp>
    #include <algebra/SmallMatrix.hpp>
    #include <mesh/MeshData.hpp>
    #include <mesh/MeshDataManager.hpp>
    #include <mesh/MeshFlatFaceBoundary.hpp>
    #include <mesh/NamedBoundaryDescriptor.hpp>
    #include <mesh/StencilDescriptor.hpp>
    #include <mesh/StencilManager.hpp>
    #include <scheme/DiscreteFunctionDPkVariant.hpp>
    #include <scheme/DiscreteFunctionUtils.hpp>
    #include <scheme/DiscreteFunctionVariant.hpp>
    #include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
    #include <scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp>
    #include <scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp>
    #include <scheme/reconstruction_utils/MutableDiscreteFunctionDPkVariant.hpp>
    
    #warning put in a file in geometry
    template <size_t Dimension>
    PUGS_INLINE auto
    symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimension>& u)
    {
      return u - 2 * dot(u, normal) * normal;
    }
    
    template <size_t Dimension>
    PUGS_INLINE auto
    symmetrize_matrix(const TinyVector<Dimension>& normal, const TinyMatrix<Dimension>& A)
    {
      const TinyMatrix S = TinyMatrix<Dimension>{identity} - 2 * tensorProduct(normal, normal);
      return S * A * S;
    }
    
    template <size_t Dimension>
    PUGS_INLINE auto
    symmetrize_coordinates(const TinyVector<Dimension>& origin,
                           const TinyVector<Dimension>& normal,
                           const TinyVector<Dimension>& u)
    {
      return u - 2 * dot(u - origin, normal) * normal;
    }
    
    size_t
    PolynomialReconstruction::_getNumberOfColumns(
      const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
    {
      size_t number_of_columns = 0;
      for (auto discrete_function_variant_p : discrete_function_variant_list) {
        number_of_columns += std::visit(
          [](auto&& discrete_function) -> size_t {
            using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
            if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
              using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
              if constexpr (std::is_arithmetic_v<data_type>) {
                return 1;
              } else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
                return data_type::Dimension;
              } else {
                // LCOV_EXCL_START
                throw UnexpectedError("unexpected data type " + demangle<data_type>());
                // LCOV_EXCL_STOP
              }
            } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
              using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
              if constexpr (std::is_arithmetic_v<data_type>) {
                return discrete_function.size();
              } else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
                return discrete_function.size() * data_type::Dimension;
              } else {
                // LCOV_EXCL_START
                throw UnexpectedError("unexpected data type " + demangle<data_type>());
                // LCOV_EXCL_STOP
              }
            } else {
              // LCOV_EXCL_START
              throw UnexpectedError("unexpected discrete function type");
              // LCOV_EXCL_STOP
            }
          },
          discrete_function_variant_p->discreteFunction());
      }
      return number_of_columns;
    }
    
    template <MeshConcept MeshType>
    std::vector<PolynomialReconstruction::MutableDiscreteFunctionDPkVariant>
    PolynomialReconstruction::_createMutableDiscreteFunctionDPKVariantList(
      const std::shared_ptr<const MeshType>& p_mesh,
      const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
    {
      std::vector<MutableDiscreteFunctionDPkVariant> mutable_discrete_function_dpk_variant_list;
      for (size_t i_discrete_function_variant = 0; i_discrete_function_variant < discrete_function_variant_list.size();
           ++i_discrete_function_variant) {
        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::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
              mutable_discrete_function_dpk_variant_list.push_back(
                DiscreteFunctionDPk<MeshType::Dimension, DataType>(p_mesh, m_descriptor.degree()));
            } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
              using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
              mutable_discrete_function_dpk_variant_list.push_back(
                DiscreteFunctionDPkVector<MeshType::Dimension, DataType>(p_mesh, m_descriptor.degree(),
                                                                         discrete_function.size()));
            } else {
              // LCOV_EXCL_START
              throw UnexpectedError("unexpected discrete function type");
              // LCOV_EXCL_STOP
            }
          },
          discrete_function_variant->discreteFunction());
      }
    
      return mutable_discrete_function_dpk_variant_list;
    }
    
    template <MeshConcept MeshType>
    void
    PolynomialReconstruction::_checkDataAndSymmetriesCompatibility(
      const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
    {
      for (auto&& discrete_function_variant : discrete_function_variant_list) {
        std::visit(
          [&](auto&& discrete_function) {
            using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
            if constexpr (is_discrete_function_P0_v<DiscreteFunctionT> or
                          is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
              using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
              if constexpr (is_tiny_vector_v<DataType>) {
                if constexpr (DataType::Dimension != MeshType::Dimension) {
                  std::stringstream error_msg;
                  error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
                            << " using a mesh of dimension " << MeshType::Dimension;
                  throw NormalError(error_msg.str());
                }
              } else if constexpr (is_tiny_matrix_v<DataType>) {
                if constexpr ((DataType::NumberOfRows != MeshType::Dimension) or
                              (DataType::NumberOfColumns != MeshType::Dimension)) {
                  std::stringstream error_msg;
                  error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
                            << DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
                  throw NormalError(error_msg.str());
                }
              }
            } else {
              // LCOV_EXCL_START
              throw UnexpectedError("invalid discrete function type");
              // LCOV_EXCL_STOP
            }
          },
          discrete_function_variant->discreteFunction());
      }
    }
    
    template <MeshConcept MeshType>
    [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
    PolynomialReconstruction::_build(
      const std::shared_ptr<const MeshType>& p_mesh,
      const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
    {
      const MeshType& mesh = *p_mesh;
    
      using Rd = TinyVector<MeshType::Dimension>;
    
      if (m_descriptor.symmetryBoundaryDescriptorList().size() > 0) {
        this->_checkDataAndSymmetriesCompatibility<MeshType>(discrete_function_variant_list);
      }
    
      const size_t number_of_columns = this->_getNumberOfColumns(discrete_function_variant_list);
    
      const size_t basis_dimension =
        DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(m_descriptor.degree());
    
      const auto& stencil_array =
        StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), m_descriptor.stencilDescriptor(),
                                                             m_descriptor.symmetryBoundaryDescriptorList());
    
      auto xr = mesh.xr();
      auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
      auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
    
      auto cell_is_owned       = mesh.connectivity().cellIsOwned();
      auto cell_type           = mesh.connectivity().cellType();
      auto cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
    
      auto full_stencil_size = [&](const CellId cell_id) {
        auto stencil_cell_list = stencil_array[cell_id];
        size_t stencil_size    = stencil_cell_list.size();
        for (size_t i = 0; i < m_descriptor.symmetryBoundaryDescriptorList().size(); ++i) {
          auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i].stencilArray();
          stencil_size += ghost_stencil[cell_id].size();
        }
    
        return stencil_size;
      };
    
      const size_t max_stencil_size = [&]() {
        size_t max_size = 0;
        for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
          const size_t stencil_size = full_stencil_size(cell_id);
          if (cell_is_owned[cell_id] and stencil_size > max_size) {
            max_size = stencil_size;
          }
        }
        return max_size;
      }();
    
      SmallArray<const Rd> symmetry_normal_list = [&] {
        SmallArray<Rd> normal_list(m_descriptor.symmetryBoundaryDescriptorList().size());
        size_t i_symmetry_boundary = 0;
        for (auto p_boundary_descriptor : m_descriptor.symmetryBoundaryDescriptorList()) {
          const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor;
    
          auto symmetry_boundary             = getMeshFlatFaceBoundary(mesh, boundary_descriptor);
          normal_list[i_symmetry_boundary++] = symmetry_boundary.outgoingNormal();
        }
        return normal_list;
      }();
    
      SmallArray<const Rd> symmetry_origin_list = [&] {
        SmallArray<Rd> origin_list(m_descriptor.symmetryBoundaryDescriptorList().size());
        size_t i_symmetry_boundary = 0;
        for (auto p_boundary_descriptor : m_descriptor.symmetryBoundaryDescriptorList()) {
          const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor;
    
          auto symmetry_boundary             = getMeshFlatFaceBoundary(mesh, boundary_descriptor);
          origin_list[i_symmetry_boundary++] = symmetry_boundary.origin();
        }
        return origin_list;
      }();
    
      Kokkos::Experimental::UniqueToken<Kokkos::DefaultExecutionSpace::execution_space,
                                        Kokkos::Experimental::UniqueTokenScope::Global>
        tokens;
    
      auto mutable_discrete_function_dpk_variant_list =
        this->_createMutableDiscreteFunctionDPKVariantList(p_mesh, discrete_function_variant_list);
    
      SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
      SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency());
      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);
        B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns);
        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);
      }
    
      std::shared_ptr p_cell_center_reconstruction_matrix_builder =
        std::make_shared<CellCenterReconstructionMatrixBuilder<MeshType>>(symmetry_origin_list, symmetry_normal_list,
                                                                          stencil_array, xj);
    
      parallel_for(
        mesh.numberOfCells(), PUGS_CLASS_LAMBDA(const CellId cell_j_id) {
          if (cell_is_owned[cell_j_id]) {
            const int32_t t = tokens.acquire();
    
            auto stencil_cell_list = stencil_array[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));
    
            if ((m_descriptor.degree() == 1) and
                (m_descriptor.integrationMethodType() == IntegrationMethodType::cell_center)) {
              p_cell_center_reconstruction_matrix_builder->build(cell_j_id, A);
            } else if ((m_descriptor.integrationMethodType() == IntegrationMethodType::element) or
                       (m_descriptor.integrationMethodType() == IntegrationMethodType::boundary)) {
    #warning implement boundary based reconstruction for 1d and 3d
              if ((m_descriptor.integrationMethodType() == IntegrationMethodType::boundary) and
                  (MeshType::Dimension == 2)) {
                if constexpr (MeshType::Dimension == 2) {
    #warning incorporate in reconstruction_matrix_builder
                  SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek = inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[t];
                  SmallArray<double>& mean_j_of_ejk                       = mean_j_of_ejk_pool[t];
                  SmallArray<double>& mean_i_of_ejk                       = mean_i_of_ejk_pool[t];
    
                  BoundaryIntegralReconstructionMatrixBuilder
                    boundary_integral_reconstruction_matrix_builder(*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,
                                                                    symmetry_normal_list, stencil_array);
    
                  boundary_integral_reconstruction_matrix_builder.build(cell_j_id, A);
    
                } else {
                  throw NotImplementedError("invalid mesh dimension");
                }
              } else {
    
    #warning incorporate in reconstruction_matrix_builder
                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];
    
                ElementIntegralReconstructionMatrixBuilder
                  element_integral_reconstruction_matrix_builder(*p_mesh, m_descriptor.degree(), inv_Vj_wq_detJ_ek,
                                                                 mean_j_of_ejk, mean_i_of_ejk, symmetry_origin_list,
                                                                 symmetry_normal_list, stencil_array);
    
                element_integral_reconstruction_matrix_builder.build(cell_j_id, A);
              }
            } else {
              throw UnexpectedError("invalid integration strategy");
            }
    
            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;
                  }
                }
              }
            }
    
            const SmallMatrix<double>& X = X_pool[t];
    
            if (m_descriptor.preconditioning()) {
              // Add column  weighting preconditioning (increase the presition)
              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;
                }
              }
            } 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());
            }
    
            tokens.release(t);
          }
        });
    
      std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> discrete_function_dpk_variant_list;
    
      for (auto discrete_function_dpk_variant_p : mutable_discrete_function_dpk_variant_list) {
        std::visit(
          [&](auto&& mutable_function_dpk) {
            synchronize(mutable_function_dpk.cellArrays());
            discrete_function_dpk_variant_list.push_back(
              std::make_shared<DiscreteFunctionDPkVariant>(mutable_function_dpk));
          },
          discrete_function_dpk_variant_p.mutableDiscreteFunctionDPk());
      }
    
      return discrete_function_dpk_variant_list;
    }
    
    std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
    PolynomialReconstruction::build(
      const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
    {
      if (not hasSameMesh(discrete_function_variant_list)) {
        throw NormalError("cannot reconstruct functions living of different meshes simultaneously");
      }
    
      auto mesh_v = getCommonMesh(discrete_function_variant_list);
    
      return std::visit([&](auto&& p_mesh) { return this->_build(p_mesh, discrete_function_variant_list); },
                        mesh_v->variant());
    }