#ifndef ALGORITHM_HH
#define ALGORITHM_HH

#include <cassert>
#include <cmath>

// iostream includes
#include <iostream>

// include grid part
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/alugrid/dgf.hh>
// include output
#include <dune/fem/io/file/dataoutput.hh>
#include <dune/fem/space/common/adaptmanager.hh>
#include <dune/fem/space/common/restrictprolongfunction.hh>
//#include <dune/fem/function/subfunction.hh>

#include "../dune/twophaseDG/models/phaseflowmodel.hh"
#include "../dune/twophaseDG/algorithms/phaseflowscheme.hh"

template <class HGridType>
class Algorithm
{
public:

  typedef Dune::Fem::FunctionSpace< double, double, HGridType::dimensionworld, 2 > FunctionSpaceType;
  typedef Dune::Fem::AdaptiveLeafGridPart< HGridType, Dune::InteriorBorder_Partition > GridPartType;

  typedef InitialData< FunctionSpaceType > ProblemType;
  typedef PhaseFlowModel< FunctionSpaceType, GridPartType > ModelType;
  typedef PhaseFlowScheme< ModelType, ModelType > SchemeType;

  typedef Dune::tuple< const typename SchemeType::DiscreteFunctionType * > IOTupleType;
  typedef Dune::Fem::DataOutput< HGridType, IOTupleType > DataOutputType;


public:
  //
  Algorithm( HGridType & grid )
    : grid_( grid )
  {
  }

  void compute ()
  {
    // create time provider
    Dune::Fem::GridTimeProvider< HGridType > timeProvider( grid_ );
    // we want to solve the problem on the leaf elements of the grid
    GridPartType gridPart( grid_ );
    // type of the mathematical model used
    ProblemType problem( timeProvider ) ;
    // implicit model for left hand side
    ModelType implicitModel( problem, gridPart, true );
    // explicit model for right hand side
    ModelType explicitModel( problem, gridPart, false );
    // create scheme
    SchemeType scheme( gridPart, implicitModel, explicitModel );
    //! input/output tuple and setup datawritter
    IOTupleType ioTuple( &(scheme.solution())) ; // tuple with pointers
    DataOutputType dataOutput( grid_, ioTuple );

    const bool local_adapt = Dune::Fem::Parameter::getValue< bool >("phaseflow.local_adapt", false );
    const double endTime  = Dune::Fem::Parameter::getValue< double >( "phaseflow.endtime", 3.0 );
    const double dtreducefactor = Dune::Fem::Parameter::getValue< double >("phaseflow.reducetimestepfactor", 1 );
    double timeStep = Dune::Fem::Parameter::getValue< double >( "phaseflow.timestep", 0.00125 );
    double tolerance = Dune::Fem::Parameter::getValue< double >("phaseflow.tolerance", 0.5 );
    int step=1;

    timeStep *= pow( dtreducefactor,step );

    // initialize with fixed time step
    timeProvider.init( timeStep ) ;
    scheme.initialize();
    // write initial solve
    dataOutput.write( timeProvider );

    // time loop, increment with fixed time step
    for( ; timeProvider.time() < endTime; timeProvider.next( timeStep ) )
      {
        std::cout << "t=" << timeProvider.time() << std::endl;

        if (local_adapt)
        {
          // Grid adaptation only WORKS with alugrid_CUBE
          // mark element for adaptation
          scheme.mark( tolerance );
          // adapt grid
          scheme.adapt();
          scheme.prepare();
          scheme.solve( true );
        }
        else
        {
          scheme.prepare();
          scheme.solve( true );
        }

        dataOutput.write( timeProvider );

      }
    dataOutput.write( timeProvider );

  }

protected:
  HGridType& grid_; // reference to grid, i.e. the hierarchical grid
};
#endif // #ifndef ALGORITHM_HH
