#ifndef PHYSICAL_PARAM_HH
#define PHYSICAL_PARAM_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 PhysicalParam
{
  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 ProblemInterface< FunctionSpaceType > ProblemType ;
  typedef  std::vector< double > GravitationVectorType;

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

public:

  //! constructor taking problem reference
  PhysicalParam( const ProblemType& problem, const GridPart &gridPart )
    : problem_( problem ),
      gridPart_(gridPart),
      gravity_( dimDomain, 0 ),
      eps_(Dune::Fem::Parameter::getValue< double >( "phaseflow.eps", 0.001 ) ),
      srn_(Dune::Fem::Parameter::getValue< double >( "phaseflow.residualnonwetting", 0 ) ),
      srw_(Dune::Fem::Parameter::getValue< double >( "phaseflow.residualwetting", 0.15 ) ),
      srwlens_(Dune::Fem::Parameter::getValue< double >( "phaseflow.lensresidualwetting", 0.15 ) ),
      visc_n_(Dune::Fem::Parameter::getValue< double >( "phaseflow.vicositynonwetting", 0.0005 ) ),
      visc_w_(Dune::Fem::Parameter::getValue< double >( "phaseflow.vicositywetting", 0.002 ) ),
      inhom_theta_(Dune::Fem::Parameter::getValue< double >( "phaseflow.inhomogeneitytheta", 2 ) ),
      lens_inhom_theta_(Dune::Fem::Parameter::getValue< double >( "phaseflow.lensinhomogeneitytheta", 2 ) ),
      gravity_value_(Dune::Fem::Parameter::getValue< double >( "phaseflow.gravity_value", -9.81 ) ),
      entry_press_(Dune::Fem::Parameter::getValue< double >( "phaseflow.entrypressure", 5e4 ) ),
      lens_entry_press_(Dune::Fem::Parameter::getValue< double >( "phaseflow.lensentrypressure", 5e4 ) ),
      abs_perm1(Dune::Fem::Parameter::getValue< double >( "phaseflow.domain1absolutepermeability", 6.64e-11 ) ),
      lens_abs_perm(Dune::Fem::Parameter::getValue< double >( "phaseflow.lensabsolutepermeability", 6.64e-11 ) ),
      porosity(Dune::Fem::Parameter::getValue< double >( "phaseflow.porosity", 1 ) ),
      lens_porosity(Dune::Fem::Parameter::getValue< double >( "phaseflow.lensporosity", 1 ) ),
      density_w(Dune::Fem::Parameter::getValue< double >( "phaseflow.densitywetting", 1 ) ),
      density_n(Dune::Fem::Parameter::getValue< double >( "phaseflow.densitynonwetting", 1 ) ),
      secondlens_entry_press_(Dune::Fem::Parameter::getValue< double >( "phaseflow.secondlensentrypressure", 5e4 ) ),
      secondlens_abs_perm(Dune::Fem::Parameter::getValue< double >( "phaseflow.secondlensabsolutepermeability", 6.64e-11 ) ),
      secondlens_porosity(Dune::Fem::Parameter::getValue< double >( "phaseflow.secondlensporosity", 1 ) ),
      srw_secondlens_(Dune::Fem::Parameter::getValue< double >( "phaseflow.secondlensresidualwetting", 0.15 ) ),
      secondlens_inhom_theta_(Dune::Fem::Parameter::getValue< double >( "phaseflow.secondlensinhomogeneitytheta", 2 ) ),
      penalty_press_( Dune::Fem::Parameter::getValue< double >( "phaseflow.penaltypress", 1 ) ),
      penalty_sat_(Dune::Fem::Parameter::getValue< double >( "phaseflow.penaltysat", 1 ) ),
      /* coordinate first lens*/
      coord_x_lens0(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_x_lens0", 1 ) ),
      coord_x_lens1(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_x_lens1", 1 ) ),
      coord_y_lens0(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_y_lens0", 1 ) ),
      coord_y_lens1(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_y_lens1", 1 ) ),
      coord_z_lens0(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_z_lens0", 0 ) ),
      coord_z_lens1(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_z_lens1", 0 ) ),
      /* coordinate second lens*/
      coord_x_secondlens0(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_x_secondlens0", 0 ) ),
      coord_x_secondlens1(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_x_secondlens1", 0 ) ),
      coord_y_secondlens0(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_y_secondlens0", 0 ) ),
      coord_y_secondlens1(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_y_secondlens1", 0 ) ),
      coord_z_secondlens0(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_z_secondlens0", 0 ) ),
      coord_z_secondlens1(Dune::Fem::Parameter::getValue< double >( "phaseflow.coord_z_secondlens1", 0 ) )
  {
    if ( dimDomain == 3 )
    {
      gravity_[dimDomain-3] = 0;
      gravity_[dimDomain-2] = 0;
      gravity_[dimDomain-1] = gravity_value_;
    }
    else if ( dimDomain == 2 )
    {
      gravity_[dimDomain-2] = 0;
      gravity_[dimDomain-1] = gravity_value_;
    }
  }

  //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
  {
    double tetha_in( 0 );
    double irred_water_sat( 0 );
    //
    inhomegeneity_tetha( entity, x, tetha_in );
    irreduct_water_satu( entity, x, irred_water_sat );
    //
    effsat = value[1]-srn_;
    effsat /= (1-irred_water_sat-srn_);
  }

  //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
  {
    double tetha_in( 0 );
    double entry_press( 0 );
    double effsat( 0 );
    //
    inhomegeneity_tetha( entity, x, tetha_in );
    entry_Pressure( entity, x, entry_press );
    effective_sat( entity, x, value, effsat );
    //
    cap_press[1] = entry_press*(std::pow((1-effsat),-1/tetha_in));
  }

  //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
  {
    double tetha_in( 0 );
    double irred_water_sat( 0 );
    double entry_press( 0 );
    double effsat( 0 );
    RangeType cappress;
    //
    irreduct_water_satu( entity, x, irred_water_sat );
    inhomegeneity_tetha( entity, x, tetha_in );
    entry_Pressure( entity, x, entry_press );
    effective_sat( entity, x, value, effsat );
    cap_Press( entity, x, value, cappress );
    //
    first_der_cap_press[0] = 0;
    first_der_cap_press[1] = (1/tetha_in)*std::pow(entry_press,-1*tetha_in)*(std::pow(cappress[1],(1+tetha_in)));
    first_der_cap_press[1] /= (1-irred_water_sat-srn_);
  }

  //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
  {
    double tetha_in( 0 );
    double entry_press( 0 );
    double irred_water_sat( 0 );
    double effsat( 0 );
    RangeType first_der_cap_pressure,cappress;
    //
    irreduct_water_satu( entity, x, irred_water_sat );
    entry_Pressure( entity, x, entry_press );
    inhomegeneity_tetha( entity, x, tetha_in );
    effective_sat( entity, x, value, effsat );
    first_der_cap_Press( entity, x, value, first_der_cap_pressure );
    cap_Press( entity, x, value, cappress );
    //
    sec_der_cap_press[0] = 0;
    sec_der_cap_press[1] = (1/tetha_in)*std::pow(entry_press,-1*tetha_in)*(1+tetha_in)*(std::pow(cappress[1],tetha_in))*first_der_cap_pressure[1];
    sec_der_cap_press[1] /= (1-irred_water_sat-srn_);
  }

  //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
  {
    double tetha_in( 0 );
    double effsat( 0 );
    RangeType firs_der_cap( 0 );
    //
    inhomegeneity_tetha( entity, x, tetha_in );
    effective_sat( entity, x, value, effsat );
    first_der_cap_Press( entity, x, value, firs_der_cap );
    //
    grad_cap_press[0] = 0;
    grad_cap_press[1] = gradient[1];
    grad_cap_press[1] *= 1*firs_der_cap[1];
  }

  //Non wetting mobility
  template< class Entity, class Point >
  void nonWet_mobility( const Entity &entity,
                        const Point &x,
                        const RangeType &value,
                        double &nwmob ) const
  {
    double tetha_in( 0 );
    double effsat( 0 );
    double Krn;
    //
    inhomegeneity_tetha( entity, x, tetha_in );
    effective_sat ( entity, x, value, effsat );
    //
    Krn = std::pow(effsat,2.)*(1-std::pow((1-effsat),(2.+tetha_in)/tetha_in));
    nwmob = Krn/visc_n_;
  }

  //Wetting mobility
  template< class Entity, class Point >
  void wet_mobility( const Entity &entity,
                     const Point &x,
                     const RangeType &value,
                     double &wmob ) const
  {
    double tetha_in( 0 );
    double effsat;
    double Krw;
    //
    inhomegeneity_tetha( entity, x, tetha_in );
    effective_sat( entity, x, value, effsat );
    //
    Krw = (std::pow((1-effsat),(2.+3*tetha_in)/tetha_in));
    wmob = Krw/visc_w_;
  }

  //Total mobility
  template< class Entity, class Point >
  void tot_mobility( const Entity &entity,
                     const Point &x,
                     const RangeType &value,
                     double &totmob ) const
  {
    double nwmob( 0 ),wmob( 0 );
    //
    wet_mobility ( entity, x, value, wmob );
    nonWet_mobility ( entity, x, value, nwmob );
    //
    totmob = nwmob;
    totmob+=wmob;
  }

  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
  {
    double permeabilityIn ( 0 ), permeabilityOut( 0 );
    double tem( 0 );
    double nwmobIn( 0 ), nwmobOut( 0 ), gradnonwetmobIn( 0 ), gradnonwetmobOut( 0 );
    JacobianRangeType upw;
    RangeType first_der_cap_pressure;
    //
    permeability_K( entityIn,xIn, permeabilityIn );
    permeability_K( entityOut,xOut, permeabilityOut );
    first_der_cap_Press( entityIn, xIn, u0En, first_der_cap_pressure );
    //
    upw = u0EnJac1;
    upw[1] *= first_der_cap_pressure[1];
    upw[0] += upw[1];
    upw[0][0] += 9.81*density_n;
    upw[0][1] += 9.81*density_n;
    upw[0] *= -1*permeabilityIn;

    tem += normal[0] * upw[0][0];
    tem += normal[1] * upw[0][1];

    nonWet_mobility( entityIn, xIn, u0En, nwmobIn );
    nonWet_mobility( entityOut, xOut, u0OutEn, nwmobOut );
    grad_nonWet_mobility( entityIn,xIn,u0En,gradnonwetmobIn );
    grad_nonWet_mobility( entityOut,xOut,u0OutEn,gradnonwetmobOut );

    if ( tem <= 0 && nwmobOut > 0 )
    {
      upwMob = nwmobOut;
      upwgradnonwetmob = gradnonwetmobOut;
    }
    else if ( tem >= 0 && nwmobOut > 0 )
    {
      upwMob = nwmobIn;
      upwgradnonwetmob = gradnonwetmobIn;
    }
  }

  //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
  {
    double tetha_in( 0 );
    double irred_water_sat( 0 );
    double effsat;
    //
    irreduct_water_satu( entity, x, irred_water_sat );
    inhomegeneity_tetha( entity, x, tetha_in );
    effective_sat( entity, x, value, effsat );
    //
    gradnonwetmob = 1*(2*(effsat)*(1-std::pow(1-effsat,(2.+tetha_in)/tetha_in)))+(1*(2.+tetha_in)/tetha_in)*std::pow(effsat,2.)*(std::pow((1-effsat),(2.)/tetha_in));
    gradnonwetmob /= visc_n_;
    gradnonwetmob /= (1-irred_water_sat-srn_);
  }

  //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
  {
    double tetha_in( 0 );
    double irred_water_sat( 0 );
    double effsat;
    //
    irreduct_water_satu( entity, x, irred_water_sat );
    inhomegeneity_tetha( entity, x, tetha_in );
    effective_sat ( entity, x, value, effsat );
    //
    gradwetmob = 0;
    gradwetmob = (-1*(2.+3*tetha_in)/tetha_in)*std::pow((1-effsat),((2.+3*tetha_in)/tetha_in)-1);
    gradwetmob /= visc_w_ ;
    gradwetmob /= (1-irred_water_sat-srn_);
  }

  //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
  {
    double gradwetmob( 0 ), gradnonwetmob( 0 );
    //
    grad_wet_mobility( entity, x, value, gradwetmob );
    grad_nonWet_mobility( entity, x, value, gradnonwetmob );
    //
    gradtotmob = gradwetmob;
    gradtotmob += gradnonwetmob;
  }

  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
  {
    betapress = penalty_press_;
  }

  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
  {
    betasat = penalty_sat_;
  }

  template< class Entity, class Point >
  void porosity_PHI( const Entity &entity,
                     const Point &x,
                     double & phi_poro) const
  {
    const DomainType xGlobal = entity.geometry().global( coordinate( x ) );

    if ( dimDomain == 3 )
    {
      if ( ( xGlobal[0] >= coord_x_lens0 && xGlobal[0] <= coord_x_lens1 ) && ( xGlobal[1] >= coord_y_lens0 && xGlobal[1] <= coord_y_lens1 ) && ( xGlobal[2] >= coord_z_lens0 && xGlobal[2] <= coord_z_lens1 ) )
      phi_poro = lens_porosity;
      else if( ( xGlobal[0] >= coord_x_secondlens0 && xGlobal[0] <= coord_x_secondlens1 ) && ( xGlobal[1] >= coord_y_secondlens0 && xGlobal[1] <= coord_y_secondlens1 ) && ( xGlobal[2] >= coord_z_secondlens0 && xGlobal[2] <= coord_z_secondlens1 ) )
      phi_poro = secondlens_porosity;
      else
      phi_poro = porosity;
    }
    else if ( dimDomain == 2 )
    {
      if( ( xGlobal[0] >= coord_x_lens0 && xGlobal[0] <= coord_x_lens1 ) && ( xGlobal[1] >= coord_y_lens0 && xGlobal[1] <= coord_y_lens1 ) )
      phi_poro = lens_porosity;
      else if ( (xGlobal[0] >= coord_x_secondlens0 && xGlobal[0] <= coord_x_secondlens1) && (xGlobal[1] >= coord_y_secondlens0 && xGlobal[1] <= coord_y_secondlens1 ) )
      phi_poro = secondlens_porosity;
      else
      phi_poro = porosity;
    }
  }

  template< class Entity, class Point >
  void entry_Pressure( const Entity &entity,
                       const Point &x,
                       double & entry_press ) const
  {
    const DomainType xGlobal = entity.geometry().global( coordinate( x ) );

    if ( dimDomain == 3)
    {
      if( (xGlobal[0] >= coord_x_lens0 && xGlobal[0] <= coord_x_lens1 ) && ( xGlobal[1] >= coord_y_lens0 && xGlobal[1] <= coord_y_lens1 ) && ( xGlobal[2] >= coord_z_lens0 && xGlobal[2] <= coord_z_lens1 ) )
      entry_press =lens_entry_press_;
      else if((xGlobal[0]>=coord_x_secondlens0 && xGlobal[0]<=coord_x_secondlens1) && ( xGlobal[1] >= coord_y_secondlens0 && xGlobal[1] <= coord_y_secondlens1 ) && ( xGlobal[2] >= coord_z_secondlens0 && xGlobal[2] <= coord_z_secondlens1 ) )
      entry_press=secondlens_entry_press_;
      else
      entry_press =entry_press_;
    }
    else if ( dimDomain == 2 )
     {
      if ( ( xGlobal[0] >= coord_x_lens0 && xGlobal[0] <= coord_x_lens1 ) && ( xGlobal[1] >= coord_y_lens0 && xGlobal[1] <= coord_y_lens1 ) )
      entry_press = lens_entry_press_;
      else if ( ( xGlobal[0] >= coord_x_secondlens0 && xGlobal[0] <= coord_x_secondlens1 ) && ( xGlobal[1] >= coord_y_secondlens0 && xGlobal[1] <= coord_y_secondlens1 ) )
      entry_press = secondlens_entry_press_;
      else
      entry_press =entry_press_;
    }
  }

  template< class Entity, class Point >
  void inhomegeneity_tetha( const Entity &entity,
                            const Point &x,
                            double & tet_inhom) const
  {
    const DomainType xGlobal = entity.geometry().global( coordinate( x ) );

    if ( dimDomain == 3 )
    {
      if ( (xGlobal[0] >= coord_x_lens0 && xGlobal[0] <= coord_x_lens1 ) && ( xGlobal[1] >= coord_y_lens0 && xGlobal[1] <= coord_y_lens1 ) && ( xGlobal[2] >= coord_z_lens0 && xGlobal[2] <= coord_z_lens1 ) )
      tet_inhom = lens_inhom_theta_;
      else if ( ( xGlobal[0] >= coord_x_secondlens0 && xGlobal[0] <= coord_x_secondlens1 ) && ( xGlobal[1] >= coord_y_secondlens0 && xGlobal[1] <= coord_y_secondlens1 ) && ( xGlobal[2] >= coord_z_secondlens0 && xGlobal[2] <= coord_z_secondlens1 ) )
      tet_inhom = secondlens_inhom_theta_;
      else
      tet_inhom = inhom_theta_;
    }
    else if ( dimDomain == 2)
     {
      if ( (xGlobal[0] >= coord_x_lens0 && xGlobal[0]<=coord_x_lens1) && (xGlobal[1]>=coord_y_lens0 && xGlobal[1]<=coord_y_lens1 ) )
      tet_inhom = lens_inhom_theta_;
      else if ( ( xGlobal[0] >= coord_x_secondlens0 && xGlobal[0] <= coord_x_secondlens1) && (xGlobal[1] >= coord_y_secondlens0 && xGlobal[1] <= coord_y_secondlens1 ) )
      tet_inhom = secondlens_inhom_theta_;
      else
      tet_inhom = inhom_theta_;
    }
  }

  template< class Entity, class Point >
  void irreduct_water_satu( const Entity &entity,
                            const Point &x,
                            double & swr) const
  {
    const DomainType xGlobal = entity.geometry().global( coordinate( x ) );

    if ( dimDomain == 3 )
     {
      if ( ( xGlobal[0] >= coord_x_lens0 && xGlobal[0] <= coord_x_lens1 ) && ( xGlobal[1] >= coord_y_lens0 && xGlobal[1] <= coord_y_lens1 ) && ( xGlobal[2] >= coord_z_lens0 && xGlobal[2] <= coord_z_lens1 ) )
      swr = srwlens_;
      else if((xGlobal[0]>=coord_x_secondlens0 && xGlobal[0]<=coord_x_secondlens1) && (xGlobal[1]>=coord_y_secondlens0 && xGlobal[1]<=coord_y_secondlens1 ) && (xGlobal[2]>=coord_z_secondlens0 && xGlobal[2]<=coord_z_secondlens1 ))
      swr =srw_secondlens_;
      else
      swr =srw_;
    }
    else if ( dimDomain==2)
     {
      if ( ( xGlobal[0] >= coord_x_lens0 && xGlobal[0] <= coord_x_lens1) && ( xGlobal[1] >= coord_y_lens0 && xGlobal[1] <= coord_y_lens1 ) )
      swr = srwlens_;
      else if ( ( xGlobal[0] >= coord_x_secondlens0 && xGlobal[0] <= coord_x_secondlens1 ) && ( xGlobal[1] >= coord_y_secondlens0 && xGlobal[1] <= coord_y_secondlens1 ) )
      swr = srw_secondlens_;
      else
      swr = srw_;
    }
  }

  template< class Entity, class Point >
  void permeability_K( const Entity &entity,
                       const Point &x,
                       double & K ) const
  {
    const DomainType xGlobal = entity.geometry().global( coordinate( x ) );

    if ( dimDomain == 3 )
     {
      if ( ( xGlobal[0] >= coord_x_lens0 && xGlobal[0] <= coord_x_lens1 ) && ( xGlobal[1] >= coord_y_lens0 && xGlobal[1] <= coord_y_lens1 ) && ( xGlobal[2] >= coord_z_lens0 && xGlobal[2] <= coord_z_lens1 ) )
      K = lens_abs_perm;
      else if((xGlobal[0]>=coord_x_secondlens0 && xGlobal[0]<=coord_x_secondlens1) && (xGlobal[1]>=coord_y_secondlens0 && xGlobal[1]<=coord_y_secondlens1 ) && (xGlobal[2]>=coord_z_secondlens0 && xGlobal[2]<=coord_z_secondlens1 ))
      K = secondlens_abs_perm;
      else
      K = abs_perm1;
    }
    else if ( dimDomain == 2 )
    {
      if ( ( xGlobal[0] >= coord_x_lens0 && xGlobal[0] <= coord_x_lens1 ) && ( xGlobal[1] >= coord_y_lens0 && xGlobal[1] <= coord_y_lens1 ) )
      K = lens_abs_perm;
      else if ( ( xGlobal[0] >= coord_x_secondlens0 && xGlobal[0] <= coord_x_secondlens1 ) && ( xGlobal[1] >= coord_y_secondlens0 && xGlobal[1] <= coord_y_secondlens1 ) )
      K = secondlens_abs_perm;
      else
      K = abs_perm1;
    }
  }

protected:

  const ProblemType& problem_;
  const GridPart &gridPart_;

  GravitationVectorType gravity_;

  const double eps_, srn_, srw_,srwlens_, visc_n_, visc_w_, inhom_theta_,lens_inhom_theta_,gravity_value_;
  const double entry_press_, lens_entry_press_, abs_perm1, lens_abs_perm, porosity, lens_porosity, density_w, density_n;
  const double  secondlens_entry_press_, secondlens_abs_perm, secondlens_porosity, srw_secondlens_,secondlens_inhom_theta_;

  const double penalty_press_,penalty_sat_;

  const double coord_x_lens0,coord_x_lens1,coord_y_lens0,coord_y_lens1,coord_z_lens0,coord_z_lens1;
  const double coord_x_secondlens0,coord_x_secondlens1,coord_y_secondlens0,coord_y_secondlens1,coord_z_secondlens0,coord_z_secondlens1;
};
#endif // #ifndef PHYSICAL_PARAM_HH
