#ifndef BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP
#define BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP

#include <algebra/ShrinkMatrixView.hpp>
#include <algebra/SmallMatrix.hpp>
#include <analysis/QuadratureFormula.hpp>
#include <geometry/LineTransformation.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/StencilArray.hpp>
#include <mesh/SubItemValuePerItem.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp>

template <MeshConcept MeshTypeT>
class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
{
 public:
  using MeshType = MeshTypeT;

  constexpr static bool handles_high_degrees = true;

 private:
  using Rd = TinyVector<MeshType::Dimension>;

  const MeshType& m_mesh;
  const size_t m_basis_dimension;
  const size_t m_polynomial_degree;

  const SmallArray<double> m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek;
  SmallArray<double> m_mean_j_of_ejk;
  SmallArray<double> m_mean_i_of_ejk;

  const ItemToItemMatrix<ItemType::cell, ItemType::face> m_cell_to_face_matrix;
  const ItemToItemMatrix<ItemType::face, ItemType::node> m_face_to_node_matrix;
  const FaceValuePerCell<const bool> m_cell_face_is_reversed;

  const CellToCellStencilArray& m_stencil_array;

  const SmallArray<const Rd> m_symmetry_origin_list;
  const SmallArray<const Rd> m_symmetry_normal_list;

  const CellValue<const double> m_Vj;
  const CellValue<const Rd> m_xj;
  const NodeValue<const Rd> m_xr;

  SmallArray<const double> m_antiderivative_coef;

  template <typename ConformTransformationT>
  void _computeEjkBoundaryMean(const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
                               const ConformTransformationT& T,
                               const Rd& Xj,
                               const double inv_Vi,
                               SmallArray<double>& mean_of_ejk);

  void _computeEjkMeanByBoundary(const Rd& Xj, const CellId& cell_id, SmallArray<double>& mean_of_ejk);

  void _computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
                                                const Rd& normal,
                                                const Rd& Xj,
                                                const CellId& cell_id,
                                                SmallArray<double>& mean_of_ejk);

 public:
  PUGS_INLINE
  SmallArray<const double>
  meanjOfEjk() const
  {
    return m_mean_j_of_ejk;
  }

  void build(const CellId cell_j_id, ShrinkMatrixView<SmallMatrix<double>>& A);

  BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh,
                                              const size_t polynomial_degree,
                                              const SmallArray<const Rd>& symmetry_origin_list,
                                              const SmallArray<const Rd>& symmetry_normal_list,
                                              const CellToCellStencilArray& stencil_array);

  BoundaryIntegralReconstructionMatrixBuilder()  = default;
  ~BoundaryIntegralReconstructionMatrixBuilder() = default;
};

#endif   // BOUNDARY_INTEGRAL_RECONSTRUCTION_MATRIX_BUILDER_HPP