#include "config.h"

#include <iostream>

#if !HAVE_DUNE_FOAMGRID
#error "Need dune-foamgrid for the definition of the convergence test"
#endif

#if !HAVE_DUNE_LOCALFUNCTIONS
#error "Need dune-localfunctions for the definition of the convergence test"
#endif

#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/test/testsuite.hh>
#include <dune/curvedgrid/curvedgrid.hh>
#include <dune/foamgrid/foamgrid.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/localfunctions/lagrange/pk.hh>

#if WRITE_FILES
#include <dune/vtk/vtkwriter.hh>
#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
#endif

// 1 .. Sphere
// 2 .. Ellipsoid
// 3 .. Torus
// 4 .. Cylinder
#ifndef GEOMETRY_TYPE
#define GEOMETRY_TYPE 1
#endif

#if GEOMETRY_TYPE == 1
#include <dune/curvedgrid/geometries/sphere.hh>
#elif GEOMETRY_TYPE == 2
#include <dune/curvedgrid/geometries/ellipsoid.hh>
#elif GEOMETRY_TYPE == 3
#include <dune/curvedgrid/geometries/torus.hh>
#elif GEOMETRY_TYPE == 4
#include <dune/curvedgrid/geometries/cylinder.hh>
#else
#error "Unknown GEOMETRY_TYPE. Must bei either 1, 2, 3, or 4"
#endif

// polynomial order of geometry
#ifndef GEOMETRY_ORDER
#define GEOMETRY_ORDER 3
#endif

#include "convergence.hh"

const int order = GEOMETRY_ORDER;
const int quad_order = order+6;
const int num_levels = 5;

int main (int argc, char** argv)
{
  using namespace Dune;
  MPIHelper::instance(argc, argv);

  using HostGrid = FoamGrid<2,3>;
#if GEOMETRY_TYPE == 1
  std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( DUNE_GRID_PATH "sphere.msh");
  auto gridFct = sphereGridFunction<HostGrid>(1.0);
#elif GEOMETRY_TYPE == 2
  std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( DUNE_GRID_PATH "ellipsoid3.msh");
  auto gridFct = ellipsoidGridFunction<HostGrid>(1.0, 0.75, 1.25);
#elif GEOMETRY_TYPE == 3
  std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( DUNE_GRID_PATH "torus2.msh");
  auto gridFct = torusGridFunction<HostGrid>(2.0, 1.0);
#elif GEOMETRY_TYPE == 4
  std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( DUNE_GRID_PATH "cylinder.msh");
  auto gridFct = cylinderGridFunction<HostGrid>(1.0);
#endif

  CurvedGrid grid(*hostGrid, gridFct, order);
  using Grid = decltype(grid);
  std::vector<typename Grid::ctype> inf_errors, L2_errors, normal_errors, curvature_errors, edge_lengths;
  for (int i = 0; i < num_levels; ++i) {
    if (i > 0)
      grid.globalRefine(1);

#if WRITE_FILES
    using DC = LagrangeDataCollector<typename Grid::LeafGridView, order>;
    VtkUnstructuredGridWriter<typename Grid::LeafGridView, DC> writer(grid.leafGridView());
    writer.write("level_" + std::to_string(i) + ".vtu");
#endif

    std::cout << "level = " << i << std::endl;
    inf_errors.push_back( inf_error(grid, gridFct.f(), quad_order) );
    L2_errors.push_back( L2_error(grid, gridFct.f(), quad_order) );
    normal_errors.push_back( normal_error(grid, gridFct.f(), quad_order) );
    curvature_errors.push_back( curvature_error(grid, gridFct.f(), quad_order) );
    edge_lengths.push_back( edge_length(grid) );
  }

  using std::log;
  using std::abs;

  // calculate the experimental order of convergence
  std::vector<typename Grid::ctype> eocInf(num_levels, 0), eocL2(num_levels, 0), eocNormal(num_levels, 0), eocCurvature(num_levels, 0);
  for (int i = 1; i < num_levels; ++i) {
    eocInf[i]       = log(inf_errors[i]/inf_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
    eocL2[i]        = log(L2_errors[i]/L2_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
    eocNormal[i]    = log(normal_errors[i]/normal_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
    eocCurvature[i] = log(curvature_errors[i]/curvature_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
  }

  TestSuite test;
  test.check(abs(eocInf.back()    - (order+1)) < 0.5, "|d|_L^inf != order+1");
  test.check(abs(eocL2.back()     - (order+1)) < 0.5, "|d|_L^2 != order+1");
  test.check(abs(eocNormal.back() - (order+0)) < 0.5, "|n - nh|_L^2 != order");
  test.check(abs(eocCurvature.back()   - (order-1)) < 0.5, "|H - Hh|_L^2 != order-1");


  auto print_line = [](auto i, const auto& data, auto accept) {
    std::cout << i;
    for (const auto& d : data) {
      std::cout << '\t' << "| ";
      if (accept(d))
        std::cout << d;
      else
        std::cout << '\t';
    }
    std::cout << std::endl;
  };

  auto print_break = [] {
    std::cout.width(8 + 9*16);
    std::cout.fill('-');
    std::cout << '-' << std::endl;
    std::cout.fill(' ');
  };

  std::cout.setf(std::ios::scientific);
  print_line("level", std::vector<std::string>{"h       ","err_inf", "eoc_inf", "err_L2", "eoc_L2", "err_norm", "eoc_norm", "err_curv", "eoc_curv"},
    [](auto&&) { return true; });
  print_break();

  for (int i = 0; i < num_levels; ++i)
    print_line(i, std::vector<typename Grid::ctype>{edge_lengths[i], inf_errors[i], eocInf[i], L2_errors[i], eocL2[i], normal_errors[i], eocNormal[i], curvature_errors[i], eocCurvature[i]},
      [](auto&&) { return true; });

  return test.exit();
}