#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <iostream>

#if HAVE_DUNE_SPGRID
#include <dune/grid/spgrid.hh>
#endif

#include <amdis/AMDiS.hpp>
#include <amdis/Integrate.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/Operators.hpp>
#include <amdis/common/Literals.hpp>

using namespace AMDiS;
using namespace Dune::Indices;

// 1 component with polynomial degree 2
using Param   = YaspGridBasis<AMDIS_DIM, 2>;
using ElliptProblem = ProblemStat<Param>;

template <class SetBoundary>
void run(SetBoundary setBoundary)
{
  ElliptProblem prob("ellipt");
  prob.initialize(INIT_ALL);
  setBoundary(prob.boundaryManager());

  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
  prob.addMatrixOperator(opL, 0, 0);

  auto f = [](auto const& x) {
    double r2 = dot(x,x);
    double ux = std::exp(-10.0 * r2);
    return -(400.0 * r2 - 20.0 * x.size()) * ux;
  };
  auto opForce = makeOperator(tag::test{}, f, 6);
  prob.addVectorOperator(opForce, 0);

  // set boundary condition
  auto g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); };
  prob.addDirichletBC(BoundaryType{1}, 0, 0, g);

  AdaptInfo adaptInfo("adapt");

  prob.assemble(adaptInfo);
  prob.solve(adaptInfo);

  double errorL2 = integrate(sqr(g - prob.solution(0)), prob.gridView(), 6);
  msg("error_L2 = {}", errorL2);
}

void run_periodic()
{
#if HAVE_DUNE_SPGRID
  Dune::SPCube<double,2> cube({0.0,0.0},{1.0,1.0});
  Dune::SPDomain<double,2> domain({cube}, Dune::SPTopology<2>(0b01));
  Dune::SPGrid<double,2> grid(domain, Dune::SPMultiIndex<2>({2,2}));
  using Grid = Dune::SPGrid<double,2>;
#else
  Dune::YaspGrid<2> grid({1.0,1.0},{2,2},std::bitset<2>("10"),0);
  using Grid = Dune::YaspGrid<2>;
#endif

  using Traits = LagrangeBasis<typename Grid::LeafGridView, 2>;
  ProblemStat<Traits> prob("ellipt", grid);
  prob.initialize(INIT_ALL);
  prob.boundaryManager().setBoxBoundary({-1,-1,1,1});

  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
  prob.addMatrixOperator(opL, 0, 0);

  auto opForce = makeOperator(tag::test{}, 1.0, 6);
  prob.addVectorOperator(opForce, 0);

  // set boundary condition
  auto g = [](auto const& x) {
    return std::sin(0.5*M_PI*x[1])*std::sin(2*M_PI * (0.25 + x[0]))
         + std::cos(0.5*M_PI*x[1])*std::sin(2*M_PI * (-0.25 + x[0]));
  };
  prob.addDirichletBC(BoundaryType{1}, 0, 0, g);

  typename ElliptProblem::WorldMatrix A{{1.0,0.0}, {0.0,1.0}};
  typename ElliptProblem::WorldVector b{1.0, 0.0};
  prob.addPeriodicBC(BoundaryType{-1}, A, b);

  AdaptInfo adaptInfo("adapt");

  prob.assemble(adaptInfo);
  prob.solve(adaptInfo);
  prob.writeFiles(adaptInfo, true);
}


int main(int argc, char** argv)
{
  AMDiS::init(argc, argv);

  auto b = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8 || x[0] > 1.0-1.e-8 || x[1] > 1.0-1.e-8; };

  auto setBoundary1 = [](auto& boundaryManager)
  {
    boundaryManager.setBoxBoundary({1,1,1,1});
  };
  run(setBoundary1);

  auto setBoundary2 = [b](auto& boundaryManager)
  {
    boundaryManager.setIndicator([b](auto const& x) { return b(x) ? 1 : 0; });
  };
  run(setBoundary2);

  auto setBoundary3 = [b](auto& boundaryManager)
  {
    boundaryManager.setPredicate(b, 1);
  };
  run(setBoundary3);

  run_periodic();

  AMDiS::finalize();
  return 0;
}