#include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>

#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/SymmetryUtils.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshDataManager.hpp>
#include <scheme/DiscreteFunctionDPk.hpp>

template <MeshConcept MeshTypeT>
template <typename ConformTransformationT>
void
PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkBoundaryMean(
  const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
  const ConformTransformationT& T,
  const Rd& Xj,
  const double inv_Vi,
  SmallArray<double>& mean_of_ejk)
{
  const double velocity_perp_e1 = T.velocity()[1] * inv_Vi;

  for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
    const double wq          = quadrature.weight(i_q);
    const TinyVector<1> xi_q = quadrature.point(i_q);

    const Rd X_Xj = T(xi_q) - Xj;

    const double x_xj = X_Xj[0];
    const double y_yj = X_Xj[1];

    {
      size_t k                                   = 0;
      m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1;
      for (; k <= m_polynomial_degree; ++k) {
        m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
          x_xj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * m_antiderivative_coef[k];   // ((1. * k) / (k + 1));
      }

      for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
        const size_t begin_i_y_1 = ((i_y - 1) * (2 * m_polynomial_degree - i_y + 4)) / 2;
        for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l) {
          m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
        }
      }
    }

    for (size_t k = 1; k < m_basis_dimension; ++k) {
      mean_of_ejk[k - 1] += m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
    }
  }
}

template <MeshConcept MeshTypeT>
void
PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkMeanByBoundary(
  const Rd& Xj,
  const CellId& cell_id,
  SmallArray<double>& mean_of_ejk)
{
  const auto& quadrature =
    QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));

  const double inv_Vi = 1. / m_Vj[cell_id];

  mean_of_ejk.fill(0);
  auto cell_face_list = m_cell_to_face_matrix[cell_id];
  for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
    bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];

    const FaceId face_id = cell_face_list[i_face];
    auto face_node_list  = m_face_to_node_matrix[face_id];
    if (is_reversed) {
      const LineTransformation<2> T{m_xr[face_node_list[1]], m_xr[face_node_list[0]]};
      _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
    } else {
      const LineTransformation<2> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]]};
      _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
    }
  }
}

template <MeshConcept MeshTypeT>
void
PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
  MeshTypeT>::_computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
                                                       const Rd& normal,
                                                       const Rd& Xj,
                                                       const CellId& cell_id,
                                                       SmallArray<double>& mean_of_ejk)
{
  const auto& quadrature =
    QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));

  const double inv_Vi = 1. / m_Vj[cell_id];

  mean_of_ejk.fill(0);
  auto cell_face_list = m_cell_to_face_matrix[cell_id];
  for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
    bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];

    const FaceId face_id = cell_face_list[i_face];
    auto face_node_list  = m_face_to_node_matrix[face_id];

    const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
    const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);

    if (is_reversed) {
      const LineTransformation<2> T{x1, x0};
      _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
    } else {
      const LineTransformation<2> T{x0, x1};
      _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
    }
  }
}

template <MeshConcept MeshTypeT>
void
PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::build(
  const CellId cell_j_id,
  ShrinkMatrixView<SmallMatrix<double>>& A)
{
  if constexpr (MeshType::Dimension == 2) {
    const auto& stencil_cell_list = m_stencil_array[cell_j_id];

    const Rd& Xj = m_xj[cell_j_id];

    _computeEjkMeanByBoundary(Xj, cell_j_id, m_mean_j_of_ejk);

    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];

      _computeEjkMeanByBoundary(Xj, cell_i_id, m_mean_i_of_ejk);

      for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
        A(index, l) = m_mean_i_of_ejk[l] - m_mean_j_of_ejk[l];
      }
    }
    for (size_t i_symmetry = 0; i_symmetry < m_stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry) {
      auto& ghost_stencil  = m_stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
      auto ghost_cell_list = ghost_stencil[cell_j_id];

      const Rd& origin = m_symmetry_origin_list[i_symmetry];
      const Rd& normal = m_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];

        _computeEjkMeanByBoundaryInSymmetricCell(origin, normal, Xj, cell_i_id, m_mean_i_of_ejk);

        for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
          A(index, l) = m_mean_i_of_ejk[l] - m_mean_j_of_ejk[l];
        }
      }
    }
  } else {
    throw NotImplementedError("invalid mesh dimension");
  }
}

template <MeshConcept MeshTypeT>
PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
  MeshTypeT>::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)
  : m_mesh{mesh},
    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{m_basis_dimension},
    m_mean_j_of_ejk{m_basis_dimension - 1},
    m_mean_i_of_ejk{m_basis_dimension - 1},
    m_cell_to_face_matrix{mesh.connectivity().cellToFaceMatrix()},
    m_face_to_node_matrix{mesh.connectivity().faceToNodeMatrix()},
    m_cell_face_is_reversed{mesh.connectivity().cellFaceIsReversed()},
    m_stencil_array{stencil_array},
    m_symmetry_origin_list{symmetry_origin_list},
    m_symmetry_normal_list{symmetry_normal_list},
    m_Vj{MeshDataManager::instance().getMeshData(mesh).Vj()},
    m_xj{MeshDataManager::instance().getMeshData(mesh).xj()},
    m_xr{mesh.xr()}
{
  if constexpr (MeshType::Dimension == 2) {
    SmallArray<double> antiderivative_coef(m_polynomial_degree + 1);
    for (size_t k = 0; k < antiderivative_coef.size(); ++k) {
      antiderivative_coef[k] = ((1. * k) / (k + 1));
    }

    m_antiderivative_coef = antiderivative_coef;
  }
}

template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::build(
  const CellId,
  ShrinkMatrixView<SmallMatrix<double>>&);

template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
  Mesh<2>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
                                                        const size_t,
                                                        const SmallArray<const Rd>&,
                                                        const SmallArray<const Rd>&,
                                                        const CellToCellStencilArray&);