#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#include <array>
#include <cmath>
#include <functional>
#include <memory>
#include <vector>

#include <dune/common/fvector.hh>
#include <dune/common/parallel/mpihelper.hh>

#ifdef HAVE_DUNE_UGGRID
#include <dune/grid/uggrid.hh>
#else
#include <dune/grid/yaspgrid.hh>
#endif

#include <amdis/ProblemStat.hpp>
#include <amdis/ProblemStatTraits.hpp>

#include "Tests.hpp"

using namespace AMDiS;

template <class BasisCreator>
auto makeGrid(bool simplex, int globalRefines = 0)
{
  using Grid = typename BasisCreator::GlobalBasis::GridView::Grid;
  static constexpr int d = Grid::dimension; // problem dimension
  using DomainType = typename Dune::FieldVector<double, d>;

  // constants
  DomainType lowerLeft; lowerLeft = 0.0;    // lower left grid corner
  DomainType upperRight; upperRight = 1.0;  // upper right grid corner
  std::array<unsigned int, d> s; s.fill(1); // number of elements on each axis

  // make grid
  std::unique_ptr<Grid> gridPtr;
  if (simplex)
  {
    gridPtr = std::unique_ptr<Grid>{
      Dune::StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, s)};
  }
  else
  {
    gridPtr = std::unique_ptr<Grid>{
      Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, s)};
  }
  gridPtr->globalRefine(globalRefines);

  return gridPtr;
}

template <class BasisCreator, class Fcts>
auto makeProblem(typename BasisCreator::GlobalBasis::GridView::Grid& grid, Fcts const& funcs)
{
  using Problem = ProblemStat<BasisCreator>;

  // make problem
  Problem prob("test", grid);
  prob.initialize(INIT_ALL);

  auto& globalBasis = *prob.globalBasis();
  auto localView = globalBasis.localView();

  // interpolate given function to initial grid
  int k = 0;
  AMDiS::forEachLeafNode_(localView.tree(), [&](auto const& node, auto tp)
  {
    interpolate(globalBasis, tp, prob.solution(tp).coefficients(), funcs[k]);
    k++;
  });

  return prob;
}

template <class Problem, class Fcts>
double calcError(Problem const& prob, Fcts const& funcs)
{
  auto& globalBasis = *prob.globalBasis();
  auto localView = globalBasis.localView();
  auto sol = prob.solution().coefficients();
  std::vector<double> ref;
  ref.resize(globalBasis.size());
  double error = 0;
  double maxError = 0;
  int k = 0;

  // interpolate given function onto reference vector
  AMDiS::forEachLeafNode_(localView.tree(), [&](auto const& node, auto tp)
  {
    interpolate(globalBasis, tp, ref, funcs[k]);
    k++;
  });

  // compare the solution with the reference
  for (std::size_t i = 0; i < ref.size(); ++i)
  {
    error = std::abs(ref[i]-sol[i]);
    maxError = std::max(maxError, error);
  }
  return maxError;
}

template <class BasisCreator>
  using Fcts = std::vector<std::function<double(
      Dune::FieldVector<double, BasisCreator::GlobalBasis::GridView::Grid::dimension>)> >;

/// Test data transfer for the case where no grid changes are made
template <class BasisCreator>
bool unchanged_test(Fcts<BasisCreator> const& funcs, bool simplex = true)
{
  using Grid = typename BasisCreator::GlobalBasis::GridView::Grid;
  static constexpr int d = Grid::dimension; // problem dimension

  auto gridPtr = makeGrid<BasisCreator>(simplex, (d > 2 ? 3 : 5));
  auto prob = makeProblem<BasisCreator, Fcts<BasisCreator>>(*gridPtr, funcs);
  // mark a single element for coarsening -> no adaptation is done
  auto e = *gridPtr->leafGridView().template begin<0>();
  gridPtr->mark(-1, e);
  AdaptInfo adaptInfo("adapt");
  prob.adaptGrid(adaptInfo);
  auto error = calcError(prob, funcs);

  return error < AMDIS_TEST_TOL;
}

/// Test data transfer for the case of global coarsening
template <class BasisCreator>
bool coarsen_test(Fcts<BasisCreator> const& funcs, bool simplex = true)
{
  using Grid = typename BasisCreator::GlobalBasis::GridView::Grid;
  static constexpr int d = Grid::dimension; // problem dimension

  auto gridPtr = makeGrid<BasisCreator>(simplex, (d > 2 ? 2 : 4));
  auto prob = makeProblem<BasisCreator, Fcts<BasisCreator>>(*gridPtr, funcs);
  prob.globalCoarsen(1);
  auto error = calcError(prob, funcs);

  return error < AMDIS_TEST_TOL;
}

/// Test data transfer for the case of global refinement
template <class BasisCreator>
bool refine_test(Fcts<BasisCreator> const& funcs, bool simplex = true)
{
  using Grid = typename BasisCreator::GlobalBasis::GridView::Grid;
  static constexpr int d = Grid::dimension; // problem dimension

  auto gridPtr = makeGrid<BasisCreator>(simplex, (d > 2 ? 1 : 3));
  auto prob = makeProblem<BasisCreator, Fcts<BasisCreator>>(*gridPtr, funcs);
  prob.globalRefine(1);
  auto error = calcError(prob, funcs);

  return error < AMDIS_TEST_TOL;
}

template <class Grid>
  using Lagrange3 = LagrangeBasis<typename Grid::LeafGridView, 3>;
template <class Grid>
  using TaylorHood = TaylorHoodBasis<typename Grid::LeafGridView>;

  constexpr double pi = 3.141592653589793238463;