#ifndef MODEL_HH
#define MODEL_HH

#include <cassert>
#include <cmath>

#include <dune/fem/solver/timeprovider.hh>
#include <dune/fem/io/parameter.hh>

#include "../algorithms/probleminterface.hh"


template< class FunctionSpace, class GridPart >
struct DiffusionModel
{
  typedef FunctionSpace FunctionSpaceType;
  typedef GridPart GridPartType;

  typedef typename FunctionSpaceType::DomainType DomainType;
  typedef typename FunctionSpaceType::RangeType RangeType;
  typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;

  typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
  typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;

  typedef TemporalProblemInterface< FunctionSpaceType > ProblemType ;

  typedef Dune::Fem::TimeProviderBase TimeProviderType;

  enum { dimRange  = FunctionSpaceType :: dimRange };
  enum { dimDomain = FunctionSpaceType :: dimDomain };

protected:
  enum FunctionId { rhs, bnd_Dir, bnd_Neum };
  template <FunctionId id>
  class FunctionWrapper;
public:
  typedef Dune::Fem::GridFunctionAdapter< FunctionWrapper<rhs>, GridPartType > RightHandSideType;
  typedef Dune::Fem::GridFunctionAdapter< FunctionWrapper<bnd_Dir>, GridPartType > DirichletBoundaryType;
  typedef Dune::Fem::GridFunctionAdapter< FunctionWrapper<bnd_Neum>, GridPartType > NeumanBoundaryType;

  //! constructor taking problem reference
  DiffusionModel( const ProblemType& problem, const GridPart &gridPart )
    : problem_( problem ),
      gridPart_( gridPart ),
      rhs_( problem_ ),
      bnd_Dirich_( problem_ ),
      bnd_Neum_( problem_ ),
      timeProvider_( problem_.timeProvider() ),
      penalty_( Dune::Fem::Parameter::getValue<double>( "dg.penalty", 0.0 ) ),
      density_w( Dune::Fem::Parameter::getValue< double >( "phaseflow.densitywetting" , 1 ) ),
      density_n( Dune::Fem::Parameter::getValue< double >( "phaseflow.densitynonwetting" , 1 ) ),
      coord_x_Dirich_bnd0( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_x_Dirich_bnd0" , 0 ) ),
      coord_x_Dirich_bnd1( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_x_Dirich_bnd1" , 0 ) ),
      coord_y_Dirich_bnd0( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_y_Dirich_bnd0" , 0 ) ),
      coord_y_Dirich_bnd1( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_y_Dirich_bnd1" , 0 ) ),
      coord_z_Dirich_bnd0( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_z_Dirich_bnd0" , 0 ) ),
      coord_z_Dirich_bnd1( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_z_Dirich_bnd1" , 0 ) ),
      coord_x_Neum_bnd0( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_x_Neum_bnd0" , 0 ) ),
      coord_x_Neum_bnd1( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_x_Neum_bnd1" , 0 ) ),
      coord_y_Neum_bnd0( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_y_Neum_bnd0" , 0 ) ),
      coord_y_Neum_bnd1( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_y_Neum_bnd1" , 0 ) ),
      coord_z_Neum_bnd0( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_z_Neum_bnd0" , 0 ) ),
      coord_z_Neum_bnd1( Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_z_Neum_bnd1" , 0 ) )

  {
  }

  //                   const Point &x,
  //                   const RangeType &value,
  //                   RangeType &flux ) const
  //  {
  //    const DomainType xGlobal = entity.geometry().global( coordinate( x ) );
  //    RangeType m;
  //    problem_.m(xGlobal,m);
  //    for (unsigned int i=0;i<flux.size();++i)
  //      flux[i] = m[i]*value[i];
  //  }
  //  //! return the diffusive flux
  //  template< class Entity, class Point >
  //  void diffusiveFlux ( const Entity &entity,
  //                       const Point &x,
  //                       const RangeType &value,
  //                       const JacobianRangeType &gradient,
  //                       JacobianRangeType &flux ) const
  //  {
  //    linDiffusiveFlux( value, gradient, entity, x, value, gradient, flux );
  //  }
  //  // linearization of diffusiveFlux
  //  template< class Entity, class Point >
  //  void linDiffusiveFlux ( const RangeType& uBar,
  //                          const JacobianRangeType& gradientBar,
  //                          const Entity &entity,
  //                          const Point &x,
  //                          const RangeType &value,
  //                          const JacobianRangeType &gradient,
  //                          JacobianRangeType &flux ) const
  //  {
  //    // the flux is simply the identity
  //    flux = gradient;
  //  }

  bool hasDirichletBoundary () const
  {
    return problem_.hasDirichletBoundary() ;
  }

  //! return true if given intersection belongs to the Dirichlet boundary -
  //! we test here if the center is a dirichlet point
  template <class Intersection>
  bool isDirichletIntersection( const Intersection& inter ) const
  {
    return isDirichletPoint( inter.geometry().center() );
  }

  //! return true if given point belongs to the Dirichlet boundary
  bool isDirichletPoint( const DomainType& x ) const
  {
    bool test( false );

    if ( dimDomain == 3)
    {
      if( x[0] == coord_x_Dirich_bnd0 || x[1] == coord_y_Dirich_bnd0 )
      {
        test = true;
      }
      else if( x[0] == coord_x_Dirich_bnd1 || x[1] == coord_y_Dirich_bnd1 )
      {
        test = true;
      }
      return test ;
    }
    else if ( dimDomain == 2 )
    {
      if( x[0] == coord_x_Dirich_bnd0 )
      {
        test = true;
      }
      else if( x[0] == coord_x_Dirich_bnd1 )
      {
        test = true;
      }
      return test ;
    }
  }

  template< class Entity, class Point >
  void g( const RangeType& uBar,
          const Entity &entity,
          const Point &x,
          RangeType &u ) const
  {
    const DomainType xGlobal = entity.geometry().global( Dune::Fem::coordinate( x ) );
    problem_.g( xGlobal, u );
  }

  template< class Entity, class Point >
  void gN( const RangeType& uBar,
           const Entity &entity,
           const Point &x,
           RangeType &u ) const
  {
    const DomainType xGlobal = entity.geometry().global( Dune::Fem::coordinate( x ) );
    u = RangeType(0);

    if ( dimDomain == 3)
     {
      if ( ( xGlobal[0] >= coord_x_Neum_bnd0 && xGlobal[0] <= coord_x_Neum_bnd1 ) && ( xGlobal[1] >= coord_y_Neum_bnd0 && xGlobal[1] <= coord_y_Neum_bnd1 ) && xGlobal[2] >= coord_z_Neum_bnd1 )
        {
          u [0] = -0.25/density_n;
          u [1] = -0.25*timeProvider_.deltaT()/density_n;
        }
      }
    else if ( dimDomain == 2 )
    {
      if ( ( xGlobal[0] >= coord_x_Neum_bnd0 && xGlobal[0] <= coord_x_Neum_bnd1 ) && xGlobal[1] >= coord_y_Neum_bnd1 )
        {
          u [0] = -0.075/density_n;
          u [1] = -0.075*timeProvider_.deltaT()/density_n;
        }
    }
  }

  // return Fem :: Function for Dirichlet boundary values
  DirichletBoundaryType dirichletBoundary( ) const
  {
    return DirichletBoundaryType( "boundary function", bnd_Dirich_, gridPart_, 5 );
  }

  bool hasNeumanBoundary () const
  {
    return problem_.hasNeumanBoundary() ;
  }

  //! return true if given intersection belongs to the Dirichlet boundary -
  //! we test here if the center is a dirichlet point
  template <class Intersection>
  bool isNeumanIntersection( const Intersection& inter ) const
  {
    return isNeumanPoint( inter.geometry().center() );
  }

  //! return true if given point belongs to the Neumann boundary (default is true)
  bool isNeumanPoint( const DomainType& x ) const
  {
    bool test( false );

    if ( dimDomain==3 )
     {
      if ( ( x[0] >= coord_x_Neum_bnd0 && x[0] <= coord_x_Neum_bnd1 ) && (x[1] >= coord_y_Neum_bnd0 && x[1] <= coord_y_Neum_bnd1 ) && x[2] >= coord_z_Neum_bnd1 )
      {
        test=true;
      }
      else if( ( x[0] <= coord_x_Neum_bnd0 || x[0] >= coord_x_Neum_bnd1 ) && ( x[1] <= coord_y_Neum_bnd0 || x[1] >= coord_y_Neum_bnd1 ) && x[2] >= coord_z_Neum_bnd1 )
      {
        test=true;
      }
      else if( x[2] <= coord_y_Neum_bnd0 )
      {
        test=true;
      }
      return test ;
    }
    else if ( dimDomain==2)
     {
      if ( ( x[0] >= coord_x_Neum_bnd0 && x[0] <= coord_x_Neum_bnd1 ) && x[1] >= coord_y_Neum_bnd1 )
      {
        test=true;
      }
      else if( ( x[0] <= coord_x_Neum_bnd0 || x[0] >= coord_x_Neum_bnd1 ) && x[1] >= coord_y_Neum_bnd1 )
      {
        test=true;
      }
      else if( x[1] <= coord_y_Neum_bnd0 )
      {
        test=true;
      }
      return test ;
    }

  }

  // return Fem :: Function for Neumann boundary values
  NeumanBoundaryType neumanBoundary( ) const
  {
    return NeumanBoundaryType( "boundary function", bnd_Neum_, gridPart_, 5 );
  }

  // return Fem :: Function for right hand side
  RightHandSideType rightHandSide(  ) const
  {
    return RightHandSideType( "right hand side", rhs_, gridPart_, 5 );
  }

  //! return reference to Problem's time provider
  const TimeProviderType & timeProvider() const
  {
    return timeProvider_;
  }
  //penalty parameter
  double penalty() const
  {
    return penalty_;
  }

protected:
  template <FunctionId id>
  class FunctionWrapper : public Dune::Fem::Function< FunctionSpaceType, FunctionWrapper< id > >
  {
    const ProblemInterface<FunctionSpaceType>& impl_;
  public:
    FunctionWrapper( const ProblemInterface<FunctionSpaceType>& impl )
      : impl_( impl ) {}

    //! evaluate function
    void evaluate( const DomainType& x, RangeType& ret ) const
    {
      if( id == rhs )
      {
        // call right hand side of implementation
        impl_.f( x, ret );
      }
      else if( id == bnd_Dir )
      {
        // call Dirichlet boudary data of implementation
        impl_.g( x, ret );
      }
      else if( id == bnd_Neum )
      {
          // call Neumann boudary data of implementation
          impl_.gN( x, ret );
        }
      else
      {
        DUNE_THROW( Dune::NotImplemented,"FunctionId not implemented" );
      }
    }
  };

  const ProblemType& problem_;
  const GridPart &gridPart_;
  FunctionWrapper<rhs> rhs_;
  FunctionWrapper<bnd_Dir> bnd_Dirich_;
  FunctionWrapper<bnd_Neum> bnd_Neum_;
  const TimeProviderType &timeProvider_;

  double penalty_;
  const double density_w;
  const double density_n;
  const double coord_x_Dirich_bnd0,coord_x_Dirich_bnd1,coord_y_Dirich_bnd0,coord_y_Dirich_bnd1,coord_z_Dirich_bnd0,coord_z_Dirich_bnd1;
  const double coord_x_Neum_bnd0,coord_x_Neum_bnd1,coord_y_Neum_bnd0,coord_y_Neum_bnd1,coord_z_Neum_bnd0,coord_z_Neum_bnd1;

};

#endif // #ifndef MODEL_HH
