#include <config.h>

#include <array>
#include <bitset>
#include <cmath>
#include <string>
#include <vector>

#include <dune/common/fvector.hh>

#include <dune/functions/functionspacebases/interpolate.hh>

#include <dune/grid/yaspgrid.hh>

#include <amdis/linearalgebra/istl/Communication.hpp>
#include <amdis/AdaptiveGrid.hpp>
#include <amdis/Environment.hpp>
#include <amdis/ProblemStatTraits.hpp>

#include "Tests.hpp"

using namespace AMDiS;

template <class Basis>
bool test(Basis& basis, std::string const& prefix)
{
  using GV = typename Basis::GridView;
  using Grid = typename GV::Grid;
  using Coord = typename Grid::template Codim<0>::Geometry::GlobalCoordinate;
  using T = double;
  static constexpr int d = Grid::dimension;

  bool passed = true;

#if HAVE_MPI
  // Interpolate test function f onto dofs and reference
  auto f = [&](const Coord& x) -> T { return std::pow(x[0],3) + x[1] + (d == 3 ? x[2] : 0.0); };
  std::vector<T> ref(basis.size());
  Dune::Functions::interpolate(basis, ref, f);
  std::vector<T> dofs = ref;

  // Exchange dofs
  basis.communicator().get().copyOwnerToAll(dofs, dofs);

  // Compare post-exchange data with reference
  for (std::size_t i = 0; i < dofs.size(); ++i)
    if (std::abs(dofs[i]-ref[i]) > AMDIS_TEST_TOL)
      passed = false;
#endif

  return passed;
}

int main(int argc, char** argv)
{
  using namespace Dune::Functions::BasisFactory;

  Environment env(argc, argv);

  if (Environment::mpiSize() == 1)
    return 0;

  using YaspGrid2d = Dune::YaspGrid< 2, Dune::EquidistantOffsetCoordinates<double, 2> >;
  using YaspGrid3d = Dune::YaspGrid< 3, Dune::EquidistantOffsetCoordinates<double, 3> >;

  // constants
  using Coord2d = Dune::FieldVector<double, 2>;
  using Coord3d = Dune::FieldVector<double, 3>;
  Coord2d lower2d(1.0);
  Coord2d upper2d(2.0);
  std::array<int, 2> nElements2d{16,3};
  Coord3d lower3d(1.0);
  Coord3d upper3d(2.0);
  std::array<int, 3> nElements3d{16,3,3};
  int ovlp = 1;

  {
    YaspGrid2d hostGrid(lower2d, upper2d, nElements2d, std::bitset<2>(), ovlp);
    AdaptiveGrid<YaspGrid2d> grid(hostGrid);

    auto l1 = LagrangeBasis<YaspGrid2d, 1>::create(grid.leafGridView());
    AMDIS_TEST(test(l1, "Yasp2d_L1_Ovlp"));

    auto th = TaylorHoodBasis<YaspGrid2d, 1>::create(grid.leafGridView());
    AMDIS_TEST(test(th, "Yasp2d_TH_Ovlp"));
  }
  {
    YaspGrid3d hostGrid(lower3d, upper3d, nElements3d, std::bitset<3>(), ovlp);
    AdaptiveGrid<YaspGrid3d> grid(hostGrid);

    auto l1 = LagrangeBasis<YaspGrid3d, 1>::create(grid.leafGridView());
    AMDIS_TEST(test(l1, "Yasp3d_L1_Ovlp"));

    auto th = TaylorHoodBasis<YaspGrid3d, 1>::create(grid.leafGridView());
    AMDIS_TEST(test(th, "Yasp3d_TH_Ovlp"));
  }

  return report_errors();
}