#ifndef ESTIMATOR_HH
#define ESTIMATOR_HH

//- Dune-fem includes
#include <dune/fem/quadrature/caching/twistutility.hh>
#include <dune/fem/quadrature/cachingquadrature.hh>
#include <dune/fem/quadrature/intersectionquadrature.hh>
#include <dune/fem/operator/common/spaceoperatorif.hh>
#include <dune/fem/operator/matrix/blockmatrix.hh>
// Local includes
#include "../models/phaseflowmodel.hh"
#include "../algorithms/phaseflowscheme.hh"


// Estimator
// ---------
template< class DiscreteFunction, class ImplicitModel, class ExplicitModel >
class Estimator
{
  typedef Estimator< DiscreteFunction,ImplicitModel,ExplicitModel > ThisType;

public:
  typedef DiscreteFunction DiscreteFunctionType;

  typedef typename DiscreteFunctionType :: DiscreteFunctionSpaceType
  DiscreteFunctionSpaceType;
  typedef typename DiscreteFunctionType :: LocalFunctionType LocalFunctionType;

  typedef typename DiscreteFunctionSpaceType :: DomainFieldType DomainFieldType;
  typedef typename DiscreteFunctionSpaceType :: RangeFieldType RangeFieldType;
  typedef typename DiscreteFunctionSpaceType :: DomainType DomainType;
  typedef typename DiscreteFunctionSpaceType :: RangeType RangeType;
  typedef typename DiscreteFunctionSpaceType :: JacobianRangeType JacobianRangeType;
  typedef typename DiscreteFunctionSpaceType :: GridPartType GridPartType;
  typedef typename DiscreteFunctionSpaceType :: IteratorType IteratorType;

  typedef typename GridPartType :: GridType GridType;
  typedef typename GridPartType :: IndexSetType IndexSetType;
  typedef typename GridPartType :: IntersectionIteratorType IntersectionIteratorType;

  typedef typename IntersectionIteratorType :: Intersection IntersectionType;

  typedef typename GridType :: template Codim< 0 > :: Entity ElementType;
  typedef typename GridType :: template Codim< 0 > :: EntityPointer ElementPointerType;
  typedef typename ElementType::Geometry GeometryType;
  static const int dimension = GridType :: dimension;

  typedef Dune::Fem::CachingQuadrature< GridPartType, 0 > ElementQuadratureType;
  typedef Dune::Fem::CachingQuadrature< GridPartType, 1 > FaceQuadratureType;

  typedef Dune::FieldMatrix<double,dimension,dimension> JacobianInverseType;
  typedef std :: vector< double > ErrorIndicatorType;

  typedef Dune::Fem::FunctionSpace< double, double, GridType::dimensionworld, 2 > FunctionSpaceType;
  // type of the mathematical model used
  typedef InitialData< FunctionSpaceType > ProblemType;
  typedef PhaseFlowModel< FunctionSpaceType, GridPartType > ModelType;

  typedef ImplicitModel ImplicitModelType;

private:
  const DiscreteFunctionType &uh_;
  const DiscreteFunctionSpaceType &dfSpace_;
  GridPartType &gridPart_;
  const IndexSetType &indexSet_;
  GridType &grid_;
  ErrorIndicatorType indicator_;
  const int maxLevel_;
  DiscreteFunctionType un_;
  ImplicitModelType model_;

protected:
  const ImplicitModelType &model () const { return model_; }
  const DiscreteFunctionType &un () const { return un_; }

public:
  explicit Estimator ( const ImplicitModelType &model, const DiscreteFunctionType &uh )
    : uh_( uh ),
      dfSpace_( uh.space() ),
      gridPart_( dfSpace_.gridPart() ),
      indexSet_( gridPart_.indexSet() ),
      grid_( gridPart_.grid() ),
      indicator_( indexSet_.size( 0 ) ),
      maxLevel_( Dune::Fem::Parameter::getValue< int >( "phasefield.maxlevel", 5 ) ), un_( uh_ ),
      model_( model )
  {}

  void setUn( DiscreteFunctionType &u )
  {
    un_.assign( u );
  }

  //! mark all elements due to given tolerance
  bool mark ( const double tolerance ) const
  {
    int marked = 0;
    // loop over all elements
    const IteratorType end = dfSpace_.end();
    for( IteratorType it = dfSpace_.begin(); it != end; ++it )
      {
        const ElementType &entity = *it;

        const Dune::ReferenceElement< double, dimension > &refElement
        = Dune::ReferenceElements< double, dimension >::general( entity.type() );
        RangeType val;
        // evaluate the phase field at the bary center (note
        // refElement.position(0,0) is the barycenter in local coordinates)
        double markVal;
        JacobianRangeType grad;
        uh_.localFunction( entity ).jacobian( refElement.position(0,0),grad );
        markVal = grad[1].two_norm();

        if( markVal > tolerance )
        {
          // make sure grid is not overly refined...
          if ( entity.level() < maxLevel_ )
          {
            // mark entity for refinement
            grid_.mark( 1, entity );
            // grid was marked
            marked = 1;
           }
         }
         else
         {
           // mark for coarsening
           grid_.mark( -1, entity );
         }
      }
    // get global max
    marked = grid_.comm().max( marked );
    return bool( marked );
  }

};

#endif // #ifndef ESTIMATOR_HH
