#ifndef FEMSCHEME_HH
#define FEMSCHEME_HH

// iostream includes
#include <iostream>

// include discrete function space
#include <dune/fem/space/lagrange.hh>

// adaptation ...
#include <dune/fem/function/adaptivefunction.hh>
#include <dune/fem/space/common/adaptmanager.hh>

// include discrete function
#include <dune/fem/function/blockvectorfunction.hh>

// include linear operators
#include <dune/fem/operator/linear/spoperator.hh>
#include <dune/fem/solver/diagonalpreconditioner.hh>

// #include <dune/fem/operator/linear/istloperator.hh>
// #include <dune/fem/solver/istlsolver.hh>
#include <dune/fem/solver/cginverseoperator.hh>
#include <dune/fem/solver/oemsolver.hh>

// lagrange interpolation
#include <dune/fem/operator/lagrangeinterpolation.hh>

// DG discrete space
#include <dune/fem/space/discontinuousgalerkin.hh>

// include norms
#include <dune/fem/misc/l2norm.hh>
#include <dune/fem/misc/h1norm.hh>

// include parameter handling
#include <dune/fem/io/parameter.hh>

// local includes
#include "probleminterface.hh"
// model
#include "../models/model.hh"
// operator
#include "../operator/operator.hh"

//----------

template < class Model>
class FemScheme
{
public:
  //! type of the mathematical model
  typedef Model ModelType ;

  //! grid view (e.g. leaf grid view) provided in the template argument list
  typedef typename ModelType::GridPartType GridPartType;

  //! type of underyling hierarchical grid needed for data output
  typedef typename GridPartType::GridType GridType;

  //! type of function space (scalar functions, \f$ f: \Omega -> R \f$)
  typedef typename ModelType :: FunctionSpaceType   FunctionSpaceType;

  //! choose type of discrete function space and the polynomial order POLORDER
#if USE_LAG_DG
  typedef Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, POLORDER > DiscreteFunctionSpaceType;
#else
  typedef Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, POLORDER >         DiscreteFunctionSpaceType;
#endif

  // choose type of discrete function, Matrix implementation and solver implementation
// #if HAVE_DUNE_ISTL && WANT_ISTL
//   typedef Dune::Fem::ISTLBlockVectorDiscreteFunction< DiscreteFunctionSpaceType >       DiscreteFunctionType;
//   typedef Dune::Fem::ISTLLinearOperator< DiscreteFunctionType, DiscreteFunctionType >   LinearOperatorType;
//   typedef Dune::Fem::ISTLGMResOp< DiscreteFunctionType, LinearOperatorType >            LinearInverseOperatorType;
// #else
  typedef Dune::Fem::AdaptiveDiscreteFunction< DiscreteFunctionSpaceType >                  DiscreteFunctionType;
  typedef Dune::Fem::SparseRowLinearOperator< DiscreteFunctionType, DiscreteFunctionType >  LinearOperatorType;
  //  typedef Dune::Fem::OEMBICGSTABOp< DiscreteFunctionType, LinearOperatorType >          LinearInverseOperatorType;
  typedef Dune::Fem::GMRESOp< DiscreteFunctionType, LinearOperatorType >                    LinearInverseOperatorType;
// #endif

  // define Jacobian operator
  typedef DifferentiableFlowOperator< LinearOperatorType, ModelType > FlowOperatorType;


  FemScheme( GridPartType &gridPart,
             const ModelType& implicitModel )
    : implicitModel_( implicitModel ),
      gridPart_( gridPart ),
      discreteSpace_( gridPart_ ),
      solution_( "solution", discreteSpace_ ),
      rhs_( "rhs", discreteSpace_ ),
      // the  operator (implicit)
      implicitOperator_( implicitModel_, solution_,discreteSpace_ ),
      // create linear operator (domainSpace,rangeSpace)
      linearOperator_( "assempled  operator", discreteSpace_, discreteSpace_ ),
      // tolerance for iterative solver
      solverEps_( Dune::Fem::Parameter::getValue< double >( "phaseflow.solvereps", 1e-8 ) )
  {
    // set all DoF to zero
    solution_.clear();
  }

  DiscreteFunctionType &solution()
  {
    return solution_;
  }
  const DiscreteFunctionType &solution() const
  {
    return solution_;
  }

  //! setup the right hand side
  void prepare()
  {
    // set boundary values for solution
    implicitOperator_.prepare( implicitModel_.dirichletBoundary(), solution_ );
    // apply constraints, e.g. Dirichlet contraints, to the result
    implicitOperator_.prepare( solution_, rhs_ );
  }

  //! solve the system
  void solve ( bool assemble )
  {
    if( assemble )
      {
        // assemble linear operator (i.e. setup matrix)
        implicitOperator_.jacobian( solution_ , linearOperator_ );
      }
    // inverse operator using linear operator
    LinearInverseOperatorType invOp( implicitOperator_, solverEps_, solverEps_ );
    // solve system
    invOp( rhs_, solution_ );
  }

protected:
  const ModelType& implicitModel_;   // the mathematical model

  GridPartType  &gridPart_;         // grid part(view)

  DiscreteFunctionSpaceType discreteSpace_; // discrete function space

  DiscreteFunctionType solution_;   // the unknown
  DiscreteFunctionType rhs_;        // the right hand side

  FlowOperatorType   implicitOperator_; // the implicit operator
  LinearOperatorType linearOperator_;  // the linear operator (i.e. jacobian of the implicit)

  const double solverEps_ ; // eps for linear solver
};

#endif // end #if FEMSCHEME_HH
