#ifndef FLOW_FEMSCHEME_HH
#define FLOW_FEMSCHEME_HH

// lagrange interpolation
#include <dune/fem/operator/lagrangeinterpolation.hh>
//#include <dune/fem/solver/newtoninverseoperator.hh>
// local includes
#include "femscheme.hh"
#include "../estimator/estimator.hh"
// l2-projection
#include <dune/fem/operator/projection/l2projection.hh>
// Newton method
#include "../operator/newtoninvop.hh"

// Scheme
//----------

template < class ImplicitModel, class ExplicitModel >
struct PhaseFlowScheme : public FemScheme<ImplicitModel>
{
  typedef FemScheme<ImplicitModel> BaseType;
  typedef typename BaseType::GridType GridType;
  typedef typename BaseType::GridPartType GridPartType;
  typedef typename BaseType::ModelType ImplicitModelType;
  typedef ExplicitModel ExplicitModelType;
  typedef typename BaseType::FunctionSpaceType FunctionSpaceType;
  typedef typename BaseType::DiscreteFunctionType DiscreteFunctionType;
  typedef Dune::Fem::NewtonInverseOperator< typename BaseType::LinearOperatorType, typename BaseType::LinearInverseOperatorType > NonLinearInverseOperatorType;

  //! type of restriction/prolongation projection for adaptive simulations
  typedef Dune::Fem::RestrictProlongDefault< DiscreteFunctionType >  RestrictionProlongationType;
  //! type of adaptation manager handling adapation and DoF compression
  typedef Dune::Fem::AdaptationManager< GridType, RestrictionProlongationType > AdaptationManagerType;
  //! Estimator
  typedef Estimator< DiscreteFunctionType,ImplicitModel, ExplicitModel > EstimatorType ;
  //!limiter
  const bool with_limiter  = Dune::Fem::Parameter::getValue< bool >( "phaseflow.with_limiter", false );

  PhaseFlowScheme( GridPartType &gridPart,
                   const ImplicitModelType& implicitModel,
                   const ExplicitModelType& explicitModel )
    : BaseType(gridPart, implicitModel),
      explicitModel_(explicitModel),
      explicitOperator_( explicitModel_, solution_, discreteSpace_),
      estimator_(implicitModel, solution_ ),
      // restriction/prolongation operator
      restrictProlong_( solution_ ),
      // adaptation manager
      adaptationManager_( gridPart_.grid(), restrictProlong_ )
  {
  }

  void prepare()
  {
    estimator_.setUn(solution_);

    explicitOperator_( solution_, rhs_ );
    // apply constraints, e.g. Dirichlet contraints, to the result
    explicitOperator_.prepare( solution_, rhs_ );
  }

  void initialize ()
  {
    Dune::Fem::DGL2ProjectionImpl::project( explicitModel_.initialFunction(), solution_ );
  }

  void solve ( bool assemble )
  {
    NonLinearInverseOperatorType invOp( implicitOperator_ );

    invOp( rhs_, solution_ );
  }

  //! mark elements for adaptation
  bool mark ( const double tolerance )
  {
    return estimator_.mark( tolerance );
  }

  //! do the adaptation for a given marking
  void adapt()
  {
    // apply adaptation and load balancing
    adaptationManager_.adapt();
  }

private:
  using BaseType::gridPart_;
  using BaseType::discreteSpace_;
  using BaseType::solution_;
  using BaseType::implicitModel_;
  using BaseType::rhs_;
  const ExplicitModelType &explicitModel_;

  typename BaseType::FlowOperatorType explicitOperator_; // the operator for the rhs
  using BaseType::implicitOperator_;                         // the operator for the lhs

  EstimatorType  estimator_;  // residual error estimator
  RestrictionProlongationType restrictProlong_ ; // local restriction/prolongation object
  AdaptationManagerType  adaptationManager_ ;    // adaptation manager handling adaptation
};

#endif // end #if FLOW_FEMSCHEME_HH
