PolynomialReconstruction.hpp

#ifndef POLYNOMIAL_RECONSTRUCTION_HPP
#define POLYNOMIAL_RECONSTRUCTION_HPP

#include <algebra/ShrinkMatrixView.hpp>
#include <algebra/ShrinkVectorView.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/SmallVector.hpp>
#include <mesh/MeshTraits.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp>

#warning MOVE TO .cpp WHEN FINISHED
#include <algebra/Givens.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/LineTransformation.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/StencilManager.hpp>

class DiscreteFunctionDPkVariant;
class DiscreteFunctionVariant;

class PolynomialReconstruction
{
 private:
  class MutableDiscreteFunctionDPkVariant;

  template <MeshConcept MeshType>
  [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build(
    const size_t degree,
    const std::shared_ptr<const MeshType>& p_mesh,
    const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const;

 public:
  template <MeshConcept MeshType, typename DataType>
  PUGS_INLINE DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>>
  build(const MeshType& mesh, const DiscreteFunctionP0<DataType> p0_function) const
  {
    using Rd = TinyVector<MeshType::Dimension>;

    const size_t degree    = 1;
    const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity());

    auto xj = MeshDataManager::instance().getMeshData(mesh).xj();

    auto cell_is_owned = mesh.connectivity().cellIsOwned();

    const size_t max_stencil_size = [&]() {
      size_t max_size = 0;
      for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
        auto stencil_cell_list = stencil[cell_id];
        if (cell_is_owned[cell_id] and stencil_cell_list.size() > max_size) {
          max_size = stencil_cell_list.size();
        }
      }
      return max_size;
    }();

    Kokkos::Experimental::UniqueToken<Kokkos::DefaultExecutionSpace::execution_space,
                                      Kokkos::Experimental::UniqueTokenScope::Global>
      tokens;

    DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>> dpk{mesh.meshVariant(), degree};

    if constexpr (is_polygonal_mesh_v<MeshType>) {
      if constexpr (std::is_arithmetic_v<DataType>) {
        SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
        SmallArray<SmallVector<double>> b_pool(Kokkos::DefaultExecutionSpace::concurrency());
        for (size_t i = 0; i < A_pool.size(); ++i) {