#ifndef PHASEFLOW_MODEL_HH
#define PHASEFLOW_MODEL_HH

#include <cassert>
#include <cmath>
#include <iostream>

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

#include "../algorithms/temporalprobleminterface.hh"
#include "../models/model.hh"

#include "../models/physicalparam.hh"


using namespace Dune::Fem;

template <class FunctionSpace>
class InitialData : public TemporalProblemInterface < FunctionSpace >
{
  typedef TemporalProblemInterface < FunctionSpace >  BaseType;
public:
  typedef typename BaseType :: RangeType            RangeType;
  typedef typename BaseType :: DomainType           DomainType;
  typedef typename BaseType :: JacobianRangeType    JacobianRangeType;

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

  // get time function from base class
  using BaseType :: time ;

  InitialData( const Dune::Fem::TimeProviderBase &timeProvider )
    : BaseType( timeProvider )
  {}

  //! the right hand side data (default = 0)
  virtual void f( const DomainType& x,
                  RangeType& phi ) const
  {
    phi[0] = 0;
    phi[1] = 0;
  }

  //! the exact solution
  virtual void u( const DomainType& x,
                  RangeType& phi ) const
  {
    if ( dimDomain == 3)
    {
      phi[0] = (1-x[2])*9810;
      phi[1] = 0;
    }
    else if ( dimDomain==2)
    {
      phi[0] = (0.65-x[1])*9810;
      phi[1] = 0;
    }
  }

  //! the jacobian of the exact solution
  virtual void uJacobian( const DomainType& x,
                          JacobianRangeType& ret ) const
  {
    // no exact solution
    ret = JacobianRangeType(0);
  }
};

template< class FunctionSpace, class GridPart >
struct PhaseFlowModel : protected DiffusionModel<FunctionSpace,GridPart>, protected PhysicalParam<FunctionSpace,GridPart>
{
  typedef DiffusionModel<FunctionSpace,GridPart> BaseType;

  typedef PhysicalParam<FunctionSpace,GridPart> BaseTypebis;

  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 typename BaseType::ProblemType InitialFunctionType;

  typedef Dune::Fem::TimeProviderBase TimeProviderType;

  typedef  std::vector< double > GravitationVectorType;

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

  PhaseFlowModel( const ProblemType& problem,
                  const GridPart &gridPart,
                  const bool implicit )
    : BaseType( problem,gridPart ),
      BaseTypebis( problem,gridPart ),
      implicit_( implicit ),
      timeStepFactor_( 0 )
  {
    // get theta for theta scheme (only used for non-linear source)
    double theta = Dune::Fem::Parameter::getValue< double >("phaseflow.theta", 1.0 );

    if (implicit)
      timeStepFactor_ = theta ;
    else
      timeStepFactor_ = -( 1.0 - theta );
  }

  //Effective saturation s_e= 1-\frac{1-s_w-s_rw}{1-s_rw -s_rn}
  template< class Entity, class Point >
  void effective_sat ( const Entity &entity,
		       const Point &x,
		       const RangeType &value,
		       double &effsat ) const
  {
    return BaseTypebis::effective_sat(entity, x, value, effsat);
  }

  //Capillary pressure function p_c=p_d (1-s_n)^{-\frac{1}{inhom_theta_}}
  template< class Entity, class Point >
  void cap_Press( const Entity &entity,
                  const Point &x,
                  const RangeType &value,
                  RangeType &cap_press ) const
  {
    return BaseTypebis::cap_Press( entity, x, value, cap_press );
  }

  //First derivative of the cappilary pressure
  template< class Entity, class Point >
  void first_der_cap_Press( const Entity &entity,
                            const Point &x,
                            const RangeType &value,
                            RangeType &first_der_cap_press ) const
  {
    return BaseTypebis::first_der_cap_Press( entity, x, value, first_der_cap_press );
  }

  //Second derivative of the cappilary pressure
  template< class Entity, class Point >
  void second_der_cap_Press( const Entity &entity,
                             const Point &x,
                             const RangeType &value,
                             RangeType &sec_der_cap_press ) const
  {
    return BaseTypebis::second_der_cap_Press( entity, x, value, sec_der_cap_press );
  }

  //gradient of the cappilary pressure
  template< class Entity, class Point >
  void grad_cap_Press ( const Entity &entity,
                        const Point &x,
                        const RangeType &value,
                        const JacobianRangeType & gradient,
                        JacobianRangeType &grad_cap_press ) const
  {
    return BaseTypebis::grad_cap_Press( entity, x, value, gradient, grad_cap_press );
  }

  //Non wetting mobility
  template< class Entity, class Point >
  void nonWet_mobility( const Entity &entity,
                        const Point &x,
                        const RangeType &value,
                        double &nwmob ) const
  {
    return BaseTypebis::nonWet_mobility( entity, x, value, nwmob );
  }

  //Wetting mobility
  template< class Entity, class Point >
  void wet_mobility( const Entity &entity,
                     const Point &x,
                     const RangeType &value,
                     double &wmob ) const
  {
    return BaseTypebis::wet_mobility( entity, x, value, wmob );
  }

  //Total mobility
  template< class Entity, class Point >
  void tot_mobility ( const Entity &entity,
                      const Point &x,
                      const RangeType &value,
                      double &totmob ) const
  {
    return BaseTypebis::tot_mobility( entity, x, value, totmob );
  }

  template< class EntityInside, class EntityOutside, class PointInside, class PointOutside , class DomainTypenormal>
  void upwindingMobility( const EntityInside &entityIn,
                          const PointInside &xIn,
                          const EntityOutside &entityOut,
                          const PointOutside &xOut,
                          const DomainTypenormal normal,
                          const RangeType& u0En,
                          const RangeType& u0OutEn,
                          const JacobianRangeType& u0EnJac1,
                          const JacobianRangeType& u0OutJac1,
                          double &upwMob,
                          double &upwgradnonwetmob) const
  {

    return BaseTypebis::upwindingMobility (entityIn, xIn, entityOut, xOut, normal, u0En, u0OutEn, u0EnJac1, u0OutJac1, upwMob, upwgradnonwetmob);
  }

  //grad of the Non wetting mobility
  template< class Entity, class Point >
  void grad_nonWet_mobility( const Entity &entity,
                             const Point &x,
                             const RangeType &value,
                             double &gradnonwetmob ) const
  {
    return BaseTypebis::grad_nonWet_mobility ( entity, x, value, gradnonwetmob );
  }

  //gradient of the wet_mobility
  template< class Entity, class Point >
  void grad_wet_mobility( const Entity &entity,
                          const Point &x,
                          const RangeType &value,
                          double &gradwetmob ) const
  {
    return BaseTypebis::grad_wet_mobility( entity, x, value, gradwetmob );
  }

  //gradient of the tot mob
  template< class Entity, class Point >
  void grad_tot_mobility( const Entity &entity,
                          const Point &x,
                          const RangeType &value,
                          double &gradtotmob ) const
  {
    return BaseTypebis::grad_tot_mobility( entity, x, value, gradtotmob );
  }

  template< class Entity, class Point >
  void source( const Entity &entity,
              const Point &x,
              const RangeType &value,
              RangeType &flux ) const
  {
    double PHI(0);
    porosity_PHI(entity, x, PHI);
    // time discretization
    flux[0] = 0;
    flux[1] = 1*PHI*value[1];
    if ( !implicit_ )
      {
        const DomainType xGlobal = entity.geometry().global( coordinate( x ) );
        // evaluate right hand side
        RangeType rhs ;
        problem_.f( xGlobal, rhs );
        rhs  *= timeProvider_.deltaT();
        flux += rhs ;
      }
  }

  template< class Entity, class Point >
  void linSource ( const RangeType& valueUn,
		   const Entity &entity,
		   const Point &x,
		   const RangeType &value,
		   RangeType &flux ) const
  {
    double PHI( 0 );
    porosity_PHI( entity, x, PHI);
    flux[0] = 0;
    flux[1] = 1*PHI*value[1];
  }

  template< class Entity, class Point >
  void gravity_Flux( const Entity &entity,
                     const Point &x,
                     const RangeType &value,
                     const JacobianRangeType &gradient,
                     JacobianRangeType &flux,
                     const double lambda_n_upw=1,
                     const bool &useupw=false,
                     const bool &buildRhs=false ) const
  {
    double K_perm( 0 );
    double lambda_n( 0 ), lambda_w( 0 );
    GravitationVectorType modif_grad ( dimDomain, 0 ), gravi( dimDomain, 0 );
    JacobianRangeType gvector( 0 );
    //
    flux = 0;
    gravi = gravity_ ;
    permeability_K( entity, x, K_perm );
    nonWet_mobility( entity, x, value, lambda_n );
    wet_mobility( entity, x, value, lambda_w );

    for (int d=0; d<dimDomain; d++)
      {
        for (int r=0; r<dimRange; r++)
          {
            gvector[r][d] = gravi[d] ;
          }
      }
    //
    flux[0].axpy( -1*K_perm*(lambda_n*density_n+lambda_w*density_w), gvector[1] );
    flux[0] *= timeStepFactor_;
    //
    if ( useupw )
    {
      flux[1].axpy( -1*K_perm*lambda_n_upw*density_n, gvector[1] );
      flux[1] *= timeStepFactor_*timeProvider_.deltaT();
    }
    else
    {
      flux[1].axpy(-1*K_perm*lambda_n*density_n, gvector[1]);
      flux[1] *= timeStepFactor_*timeProvider_.deltaT();
    }
  }

  //! return the diffusive flux
  template< class Entity, class Point >
  void wetting_Pressure_DiffusiveFlux( const Entity &entity,
                                       const Point &x,
                                       const RangeType &value,
                                       const JacobianRangeType &gradient,
                                       JacobianRangeType &flux,
                                       const double lambda_n_upw=1,
                                       const bool &useupw=false,
                                       const bool &buildRhs=false ) const
  {
    double K_perm(0);
    double lambda_n(0), lambda_t(0);
    //
    flux = 0;
    permeability_K(entity, x, K_perm);
    tot_mobility ( entity, x, value, lambda_t);
    nonWet_mobility ( entity, x, value, lambda_n);
    //
    if ( buildRhs == false )
    {
      flux[0].axpy(1*K_perm*lambda_t, gradient[0]);
      flux[0]*=timeStepFactor_;
      //
      if ( useupw )
      {
        flux[1].axpy(1*K_perm*lambda_n_upw, gradient[0]);
        flux[1] *= timeStepFactor_*timeProvider_.deltaT();
      }
      else
      {
        flux[1].axpy(1*K_perm*lambda_n, gradient[0]);
        flux[1] *= timeStepFactor_*timeProvider_.deltaT();
      }
    }
    else
    {
      flux[0].axpy(1*K_perm*lambda_t,gradient[0]);
      //
      flux[1].axpy(1*K_perm*lambda_n,gradient[0]);
      flux[1] *= 1*timeProvider_.deltaT();
    }
  }

  //! return the diffusive flux
  template< class Entity, class Point >
  void capillary_Pressure_DiffusiveFlux( const Entity &entity,
                                         const Point &x,
                                         const RangeType &value,
                                         const JacobianRangeType &gradient,
                                         JacobianRangeType &flux,
                                         const double lambda_n_upw=1,
                                         const bool &useupw=false,
                                         const bool &buildRhs=false ) const
  {
    double K_perm( 0 );
    double lambda_n( 0 );
    JacobianRangeType grad_cap_press( 0 );
    //
    flux = 0;
    permeability_K( entity, x, K_perm );
    nonWet_mobility( entity, x, value, lambda_n );
    grad_cap_Press( entity, x, value, gradient, grad_cap_press );
    //
    if ( buildRhs == false )
    {
      flux[0].axpy( 1*K_perm*lambda_n,grad_cap_press[1] );
      flux[0] *= timeStepFactor_;
      //
      flux[1].axpy( 1*K_perm*lambda_n,grad_cap_press[1] );
      flux[1] *= timeStepFactor_*timeProvider_.deltaT();
    }
    else
    {
      flux[0].axpy( 1*K_perm*lambda_n,grad_cap_press[1] );
      flux[0]*=1;
      //
      flux[1].axpy( 1*K_perm*lambda_n,grad_cap_press[1] );
      flux[1] *= 1*timeProvider_.deltaT();
    }
  }

  //! return the diffusive flux
  template< class Entity, class Point >
  void linGravity_Flux( const RangeType &valueUn,
                        const JacobianRangeType &gradientUn,
                        const Entity &entity,
                        const Point &x,
                        const RangeType &value,
                        const JacobianRangeType &gradient,
                        JacobianRangeType &flux,
                        const double lambda_n_upw = 1,
                        const double gradnonwetmob_upw = 1,
                        const bool useupw=false ) const
  {
    double K_perm( 0 );
    double gradwetmob( 0 ), gradnonwetmob( 0 ),gradtotmob( 0 );
    GravitationVectorType modif_grad ( dimDomain, 0 ), gravi( dimDomain, 0 );
    JacobianRangeType gvector( 0 );
    //
    flux = 0;
    permeability_K( entity, x, K_perm );
    gravi = gravity_ ;

    for (int d=0; d<dimDomain; d++)
      {
        for (int r=0; r<dimRange; r++)
          {
            gvector[r][d] = gravi[d];
          }
      }
    //
    grad_wet_mobility( entity, x, valueUn, gradwetmob );
    grad_nonWet_mobility ( entity, x, valueUn, gradnonwetmob );
    grad_tot_mobility ( entity, x, valueUn, gradtotmob );
    //
    flux[0].axpy(-1*K_perm*(gradnonwetmob*density_n + gradwetmob*density_w)*value[1], gvector[1]);
    flux[0] *= timeStepFactor_;
    //
    if ( useupw )
    {
      flux[1].axpy(-1*K_perm*gradnonwetmob_upw*density_n*value[1], gvector[1]);
      flux[1] *= timeStepFactor_*timeProvider_.deltaT();
    }
    else{
      flux[1].axpy(-1*K_perm*gradnonwetmob*density_n*value[1], gvector[1]);
      flux[1] *= timeStepFactor_*timeProvider_.deltaT();
    }
  }

  //! return the diffusive flux
  template< class Entity, class Point >
  void linWetting_Pressure_DiffusiveFlux( const RangeType &valueUn,
                                          const JacobianRangeType &gradientUn,
                                          const Entity &entity,
                                          const Point &x,
                                          const RangeType &value,
                                          const JacobianRangeType &gradient,
                                          JacobianRangeType &flux,
                                          const double lambda_n_upw = 1,
                                          const double gradnonwetmob_upw = 1,
                                          const bool useupw=false ) const
  {
    double K_perm( 0 );
    double lambda_n( 0 ), lambda_t( 0 );
    double gradwetmob( 0 ), gradnonwetmob( 0 ),gradtotmob( 0 );
    //
    flux = 0;
    permeability_K( entity, x, K_perm );
    tot_mobility( entity, x, valueUn, lambda_t);
    nonWet_mobility( entity, x, valueUn, lambda_n);
    grad_wet_mobility( entity, x, valueUn, gradwetmob );
    grad_nonWet_mobility( entity, x, valueUn, gradnonwetmob );
    grad_tot_mobility( entity, x, valueUn, gradtotmob );
    //
    flux[0].axpy( 1*K_perm*lambda_t, gradient[0] );
    flux[0].axpy( 1*K_perm*gradtotmob*value[1], gradientUn[0] );
    flux[0] *= timeStepFactor_;
    //
    if ( useupw )
     {
      flux[1] = gradient[0];
      flux[1] *= 1*K_perm*lambda_n_upw;
      flux[1].axpy( 1* K_perm*gradnonwetmob_upw*value[1], gradientUn[0] );
      flux[1] *= timeStepFactor_*timeProvider_.deltaT();
    }
    else
    {
      flux[1] = gradient[0];
      flux[1] *= 1*K_perm*lambda_n;
      flux[1].axpy(1* K_perm*gradnonwetmob*value[1], gradientUn[0]);
      flux[1] *= timeStepFactor_*timeProvider_.deltaT();
    }
  }

  //! return the diffusive flux
  template< class Entity, class Point >
  void linCapillary_Pressure_DiffusiveFlux( const RangeType &valueUn,
                                            const JacobianRangeType &gradientUn,
                                            const Entity &entity,
                                            const Point &x,
                                            const RangeType &value,
                                            const JacobianRangeType &gradient,
                                            JacobianRangeType &flux ) const
  {
    double K_perm( 0 );
    double lambda_n( 0 );
    double gradnonwetmob( 0 );
    RangeType dercap_press( 0 ), secdercap_press( 0 );
    JacobianRangeType grad_cap_press( 0 );
    //
    flux = 0;
    permeability_K( entity, x, K_perm );
    nonWet_mobility( entity, x, valueUn, lambda_n );
    grad_cap_Press( entity, x, valueUn, gradientUn, grad_cap_press );
    first_der_cap_Press( entity, x, valueUn, dercap_press );
    second_der_cap_Press( entity, x, valueUn, secdercap_press );
    grad_nonWet_mobility( entity, x, valueUn, gradnonwetmob );
    //
    flux[0] = grad_cap_press[1];
    flux[0] *= (1* K_perm*gradnonwetmob*value[1]);
    flux[0].axpy(1* K_perm*lambda_n*dercap_press[1],gradient[1]);
    flux[0].axpy(1* K_perm*lambda_n*value[1]*secdercap_press[1],gradientUn[1]);
    flux[0] *= timeStepFactor_;
    //
    flux[1] = grad_cap_press[1];
    flux[1] *= (1* K_perm*gradnonwetmob*value[1]);
    flux[1].axpy(1* K_perm*lambda_n*dercap_press[1],gradient[1]);
    flux[1].axpy(1* K_perm*lambda_n*value[1]*secdercap_press[1],gradientUn[1]);
    flux[1] *= timeStepFactor_*timeProvider_.deltaT();
  }

  bool hasDirichletBoundary () const
  {
    return BaseType::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
  {
    bool test;
    test = isDirichletPoint( inter.geometry().center() );
    return test;
  }

  //! return true if given point belongs to the Dirichlet boundary (default is true)
  bool isDirichletPoint( const DomainType& x ) const
  {
    return BaseType::isDirichletPoint( x ) ;
  }

  bool hasNeumanBoundary () const
  {
    return BaseType::hasNeumanBoundary() ;
  }

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

  //! return true if given point belongs to the Neuman boundary (default is true)
  bool isNeumanPoint( const DomainType& x ) const
  {
    return BaseType::isNeumanPoint( x );
  }

  //! return true if RHS
  bool isRHS() const
  {
    return !implicit_ ;
  }

  template< class Entity, class Point >
  void permeability_K( const Entity &entity,
                       const Point &x,
                       double & K_perm ) const
  {
    return BaseTypebis::permeability_K(entity,x,K_perm);
  }

  template< class Entity, class Point >
  void porosity_PHI( const Entity &entity,
                     const Point &x,
                     double & phi_poro ) const
  {
    return BaseTypebis::porosity_PHI( entity, x, phi_poro );
  }

  template< class Entity, class Point, class Entityneigh, class Pointneigh >
  void penality_pressure( const Entity &entity,
                          const Point &x1,
                          const Entityneigh &entityneigh,
                          const Pointneigh &x2,
                          double & betapress ) const
  {
    return BaseTypebis::penality_pressure( entity, x1, entityneigh, x2, betapress );
  }

  template< class Entity, class Point,class Entityneigh, class Pointneigh >
  void penality_saturation( const Entity &entity,
                            const Point &x1,
                            const Entityneigh &entityneigh,
                            const Pointneigh &x2,
                            double & betasat ) const
  {
    return BaseTypebis::penality_saturation( entity, x1, entityneigh, x2, betasat );
  }

  template< class Entity, class Point >
  void entry_Pressure( const Entity &entity,
                       const Point &x,
                       double & entry_press ) const
  {
    return BaseTypebis::entry_Pressure( entity, x, entry_press );
  }

  template< class Entity, class Point >
  void inhomegeneity_tetha( const Entity &entity,
                            const Point &x,
                            double & tet_inhom ) const
  {
    return BaseTypebis::inhomegeneity_tetha( entity, x, tet_inhom );
  }

  template< class Entity, class Point >
  void irreduct_water_satu( const Entity &entity,
                            const Point &x,
                            double & swr ) const
  {
    return BaseTypebis::irreduct_water_satu( entity, x, swr );
  }

  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( coordinate( x ) );

    BaseType::g( uBar, entity, x, u);
  }

  template< class Entity, class Point >
  void gN( const RangeType& uBar,
           const Entity &entity,
           const Point &x,
           RangeType &u ) const
  {
    BaseType::gN( uBar, entity, x, u);
  }

  // return Fem::Function for Dirichlet boundary values
  typename BaseType::DirichletBoundaryType dirichletBoundary() const
  {
    return BaseType::dirichletBoundary();
  }

  // return Fem::Function for Neumann boundary values
  typename BaseType::NeumanBoundaryType neumanBoundary() const
  {
    return BaseType::NeumanBoundaryType();
  }

  //! return reference to Problem's time provider
  const TimeProviderType & timeProvider() const
  {
    return  BaseType::timeProvider();
  }

  const InitialFunctionType &initialFunction() const
  {
    return problem_;
  }

protected:
  using BaseType::problem_;
  using BaseTypebis::gravity_;

  using BaseType::timeProvider_;
  bool implicit_;
  double timeStepFactor_;

  using BaseTypebis::density_w;
  using BaseTypebis::density_n;
};
#endif // #ifndef PHASEFLOW_MODEL_HH
