#ifndef ODE_DG_1D
#define ODE_DG_1D

#include <rang.hpp>

#include <utils/ArrayUtils.hpp>

#include <utils/PugsAssert.hpp>

#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshNodeBoundary.hpp>

#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>

#include <analysis/Polynomial.hpp>
#include <analysis/PolynomialBasis.hpp>

#include <mesh/ItemValueUtils.hpp>
#include <mesh/SubItemValuePerItem.hpp>

#include <scheme/DiscontinuousGalerkin1DTools.hpp>
#include <utils/Exceptions.hpp>
#include <utils/Messenger.hpp>

#include <iostream>
enum class FluxType
{
  centered,
  stable
};

template <size_t Degree>
class ODEDG1D
{
  using MeshType                    = Mesh<Connectivity1D>;
  constexpr static size_t Dimension = MeshType::Dimension;
  static_assert(Dimension == 1, "Galerkin 1D only works in dimension 1");
  using MeshDataType = MeshData<Dimension>;
  //  using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;

  const BasisType& m_basis_type;
  std::shared_ptr<const MeshType> m_mesh;
  const typename MeshType::Connectivity& m_connectivity;

  // class DirichletBoundaryCondition;

  // using BoundaryCondition = std::variant<DirichletBoundaryCondition>;

  // using BoundaryConditionList = std::vector<BoundaryCondition>;
  // DiscontinuousGalerkin1DTools<Degree>& m_dgtools;
  using Rd  = TinyVector<Degree>;
  using Rdd = TinyMatrix<Degree>;

  PolynomialBasis<Degree> m_basis;
  using R1 = TinyVector<Dimension>;
  // const BoundaryConditionList m_boundary_condition_list;
  // NodeValue<R1> m_flux;
  TinyVector<Degree + 1>
  _buildZeros(double xmin, double xmax)
  {
    TinyVector<Degree + 1> zeros(zero);
    if constexpr (Degree == 0) {
      zeros[0] = 0.5 * (xmax - xmin);
    } else {
      zeros[0]  = xmin;
      double dx = (xmax - xmin) / Degree;
      for (size_t i = 1; i <= Degree; ++i) {
        zeros[i] = xmin + i * dx;
      }
    }
    return zeros;
  }

 public:
  // void
  // _applyBC()
  // {
  //   for (const auto& boundary_condition : m_boundary_condition_list) {
  //     std::visit(
  //       [&](auto&& bc) {
  //         using T = std::decay_t<decltype(bc)>;
  //         if constexpr (std::is_same_v<DirichletBoundaryCondition, T>) {
  //           const auto& node_list  = bc.nodeList();
  //           const auto& value_list = bc.valueList();

  //           parallel_for(
  //             node_list.size(), PUGS_LAMBDA(size_t i_node) {
  //               NodeId node_id    = node_list[i_node];
  //               const auto& value = value_list[i_node];

  //               m_flux[node_id] = value;
  //             });
  //         }
  //       },
  //       boundary_condition);
  //   }
  // }

  // NodeValue<R1>
  // computeFluxes(FluxType flux_type, const PolynomialBasis<Degree> Basis, Polynomial<Degree>& U)
  // {
  //   switch (flux_type) {
  //   // case FluxType::centered: {
  //   //   return computeCenteredFlux(Basis, U);
  //   //   break;
  //   // }
  //   case FluxType::stable: {
  //     return computeStableFlux(Basis, U);
  //     break;
  //   }
  //   default:
  //     throw UnexpectedError("unknown flux type");
  //   }
  // }

  double
  computeStableFluxes(const CellValue<const Polynomial<Degree>>& U, NodeId r)
  {
    double flux;
    const NodeValue<const R1> xr    = m_mesh->xr();
    const auto& node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix();

    const auto& node_cells = node_to_cell_matrix[r];
    if (node_cells.size() == 1) {
      flux = 1.;
    } else {
      const CellId j1  = node_cells[0];
      const double xr0 = xr[r][0];
      flux             = U[j1](xr0);
      std::cout << "U[" << j1 << "]=" << U[j1] << "\n";
    }
    return flux;
  }

  void
  globalSolve(CellValue<Polynomial<Degree>>& U)
  {
    const NodeValue<const R1> xr    = m_mesh->xr();
    const auto& cell_to_node_matrix = m_mesh->connectivity().cellToNodeMatrix();
    CellValue<TinyVector<Degree + 1>> moments(m_connectivity);
    for (CellId j = 0; j < m_mesh->numberOfCells(); j++) {
      const auto& cell_nodes             = cell_to_node_matrix[j];
      const NodeId r1                    = cell_nodes[0];
      const NodeId r2                    = cell_nodes[1];
      const TinyVector<Dimension> xr1    = xr[r1];
      const TinyVector<Dimension> xr2    = xr[r2];
      const TinyVector<Degree + 1> zeros = _buildZeros(xr1[0], xr2[0]);
      std::cout << "x1 " << xr1[0] << " x2 " << xr2[0] << " zeros " << zeros << "\n";
      m_basis.build(m_basis_type, 0.5 * (xr1[0] + xr2[0]), zeros);
      std::cout << " B[" << j << "]=" << m_basis.elements()[0] << "\n";
      const double fluxg = computeStableFluxes(U, r1);
      std::cout << j << " flux " << fluxg << "\n";
      moments[j] = buildMoments(m_basis, fluxg, xr1[0]);
      std::cout << " moments[" << j << "]=" << moments[j] << "\n";

      U[j] = elementarySolve(moments[j], m_basis, xr1[0], xr2[0]);
      std::cout << " U[" << j << "]=" << U[j] << "\n";
    }
  }

  // void
  // integrate() const
  // {}
  PUGS_INLINE constexpr TinyVector<Degree + 1>
  buildMoments(PolynomialBasis<Degree> basis, double fluxg, double xr1)
  {
    TinyVector<Degree + 1> moments(zero);
    for (size_t i = 0; i <= Degree; i++) {
      // moments[i] = (basis.elements()[i](xr2) - basis.elements()[i](xr1)) * fluxg;
      moments[i] = -basis.elements()[i](xr1) * fluxg;
    }
    return moments;
  }

  PUGS_INLINE constexpr Polynomial<Degree>
  elementarySolve(const TinyVector<Degree + 1> moments,
                  const PolynomialBasis<Degree> Basis,
                  const double& xinf,
                  const double& xsup) const
  {
    Polynomial<Degree> U;
    U.coefficients() = zero;
    TinyMatrix<Degree + 1> M(zero);
    for (size_t i = 0; i <= Degree; ++i) {
      for (size_t j = 0; j <= Degree; ++j) {
        Polynomial<2 * Degree> P = Basis.elements()[j] * Basis.elements()[i];
        P += derivative(Basis.elements()[i]) * Basis.elements()[j];
        double fluxd = Basis.elements()[i](xsup) * Basis.elements()[j](xsup);
        M(i, j)      = integrate(P, xinf, xsup) - fluxd;
      }
    }
    std::cout << " M " << M << "\n";
    TinyMatrix inv_M                    = ::inverse(M);
    TinyVector<Degree + 1> coefficients = inv_M * moments;
    for (size_t i = 0; i <= Degree; ++i) {
      U += coefficients[i] * Basis.elements()[i];
    }
    std::cout << "U(" << xsup << ")=" << U(xsup) << "\n";
    return U;
  }

 public:
  // TODO add a flux manager to constructor
  // TODO add a basis to constructor
  ODEDG1D(const BasisType& basis_type, std::shared_ptr<const IMesh> i_mesh)   //, const BoundaryConditionList& bc_list)
  : m_basis_type(basis_type),
    m_mesh(std::dynamic_pointer_cast<const MeshType>(i_mesh)),
    m_connectivity(m_mesh->connectivity())   //,
  // m_boundary_condition_list(bc_list)
  {
    ;
  }
};

#endif   // ODE_DG_1D