#ifndef OPERATOR_HH
#define OPERATOR_HH
//- Dune includes
#include <dune/common/fmatrix.hh>
#include <dune/fem/quadrature/cachingquadrature.hh>
#include <dune/fem/operator/common/operator.hh>
#include <dune/fem/operator/common/stencil.hh>
#include <dune/fem/operator/common/differentiableoperator.hh>

// FlowOperator
// ----------------
template<class DiscreteFunction, class Model >
struct FlowOperator
  : public virtual Dune::Fem::Operator< DiscreteFunction >
{
  typedef DiscreteFunction DiscreteFunctionType;
  typedef Model            ModelType;

protected:
  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
  typedef typename DiscreteFunctionType::LocalFunctionType         LocalFunctionType;
  typedef typename LocalFunctionType::RangeType                    RangeType;
  typedef typename LocalFunctionType::JacobianRangeType            JacobianRangeType;

  typedef typename DiscreteFunctionSpaceType::IteratorType         IteratorType;
  typedef typename IteratorType::Entity                            EntityType;
  typedef typename EntityType::Geometry                            GeometryType;

  typedef typename DiscreteFunctionSpaceType::DomainType DomainType;

  typedef typename DiscreteFunctionSpaceType::GridPartType  GridPartType;

  typedef typename EntityType::EntityPointer  EntityPointerType;

  typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
  typedef typename IntersectionIteratorType::Intersection IntersectionType;
  typedef typename IntersectionType::Geometry             IntersectionGeometryType;

  typedef Dune::Fem::ElementQuadrature< GridPartType, 1 > FaceQuadratureType;

  typedef Dune::Fem::CachingQuadrature< GridPartType, 0 > QuadratureType;

  static const int dimDomain = LocalFunctionType::dimDomain;
  static const int dimRange = LocalFunctionType::dimRange;

public:
  //! contructor
  FlowOperator( const ModelType &model,
          const DiscreteFunctionType &solution )
    : model_( model )
    {}

  // prepare the solution vector
  template <class Function>
  void prepare( const Function &func,
          DiscreteFunctionType &u )
    {
    }
  //! application operator
  virtual void operator() ( const DiscreteFunctionType &u,
          DiscreteFunctionType &w ) const;

protected:
  const ModelType &model () const { return model_; }

  double penalty() const { return model_.penalty(); }

  void volumetricPart( const EntityType& entity,
           const int& quadOrder,
           const GeometryType &geometry,
           const LocalFunctionType& uLocal,
           LocalFunctionType& wLocalEn ) const;

  void internal_Bnd_terms( const EntityType& entity,
         const int& quadOrder,
         const GeometryType &geometry,
         const IntersectionType &intersection,
         const DiscreteFunctionSpaceType &dfSpace,
         const LocalFunctionType& uLocal,
         const DiscreteFunctionType& u,
         LocalFunctionType& wLocalEn) const;

public:
  //eps NIPG SIPG  (epsil=0 BOmethod, epsil=1 SIPGmethod, epsil=-1 NIPGmethod)
  const int epsil = Dune::Fem::Parameter::getValue< double >( "phaseflow.DGeps", 1 );
  const bool with_upw = Dune::Fem::Parameter::getValue< bool >( "phaseflow.with_upw", false );

private:
  ModelType model_;
};

// DifferentiableFlowOperator
// ------------------------------

template< class JacobianOperator, class Model >
struct DifferentiableFlowOperator
  : public FlowOperator< typename JacobianOperator::DomainFunctionType, Model >,
    public Dune::Fem::DifferentiableOperator< JacobianOperator >
{
  typedef FlowOperator< typename JacobianOperator::DomainFunctionType, Model > BaseType;

  typedef JacobianOperator JacobianOperatorType;

  typedef typename BaseType::DiscreteFunctionType DiscreteFunctionType;
  typedef typename BaseType::ModelType            ModelType;

protected:
  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
  typedef typename DiscreteFunctionType::LocalFunctionType         LocalFunctionType;
  typedef typename LocalFunctionType::RangeType                    RangeType;
  typedef typename LocalFunctionType::JacobianRangeType            JacobianRangeType;

  typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
  typedef typename IteratorType::Entity       EntityType;
  typedef typename EntityType::Geometry       GeometryType;

  typedef typename DiscreteFunctionSpaceType :: DomainType DomainType;

  typedef typename DiscreteFunctionSpaceType::GridPartType  GridPartType;

  typedef typename BaseType::QuadratureType QuadratureType;

  typedef typename LocalFunctionType::RangeFieldType RangeFieldType;

  typedef typename EntityType::EntityPointer  EntityPointerType;

  typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
  typedef typename IntersectionIteratorType::Intersection IntersectionType;
  typedef typename IntersectionType::Geometry             IntersectionGeometryType;

  typedef Dune::Fem::ElementQuadrature< GridPartType, 1 > FaceQuadratureType;

  static const int dimDomain = LocalFunctionType::dimDomain;
  static const int dimRange = LocalFunctionType::dimRange;

  typedef typename JacobianOperator::LocalMatrixType LocalMatrixType;
  typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType BasisFunctionSetType;

  typedef typename DiscreteFunctionSpaceType :: DomainFieldType DomainFieldType;
  typedef Dune::FieldVector< DomainType , dimRange > DeoModType;
  typedef Dune::FieldMatrix< DomainFieldType, dimDomain , dimDomain > MatrixType;
  typedef Dune::DGFEntityKey<int> KeyType;
  // global coordinates
  typedef typename GeometryType::GlobalCoordinate GlobalCoordinateType;
  // local coordinates
  typedef typename GeometryType::LocalCoordinate LocalCoordinateType;
  //

public:

  DifferentiableFlowOperator( const ModelType &model, const DiscreteFunctionType &solution,DiscreteFunctionSpaceType &space)
    : BaseType( model, solution )/*, stencil_(space,space)*/
    {}

  //! method to setup the jacobian of the operator for storage in a matrix
  void jacobian( const DiscreteFunctionType &u, JacobianOperatorType &jOp ) const;

protected:

  void element_Local_Contribution( const EntityType& entity,
                                   const GeometryType &geometry,
                                   const DiscreteFunctionSpaceType &dfSpace,
                                   const BasisFunctionSetType &baseSet,
                                   const LocalFunctionType& uLocal,
                                   LocalMatrixType& jLocal ) const;

  void interface_Local_Contribution( const EntityType& entity,
                                     const GeometryType &geometry,
                                     const DiscreteFunctionSpaceType &dfSpace,
                                     const BasisFunctionSetType &baseSet,
                                     const IntersectionType& intersection,
                                     const LocalFunctionType& uLocal,
                                     const DiscreteFunctionType& u,
                                     const JacobianOperatorType &jOp,
                                     LocalMatrixType& jLocal ) const;

protected:
  using BaseType::model;
  using BaseType::penalty;
  using BaseType::epsil;
  using BaseType::with_upw;

};

template< class DiscreteFunction, class Model >
void FlowOperator< DiscreteFunction, Model >
:: volumetricPart( const EntityType& entity,
                   const int& quadOrder,
                   const GeometryType &geometry,
                   const LocalFunctionType& uLocal,
                   LocalFunctionType& wLocalEn ) const
{
  QuadratureType quadrature( entity, quadOrder );
  const size_t numQuadraturePoints = quadrature.nop();
  for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
  {
    const typename QuadratureType::CoordinateType &x = quadrature.point( pt );
    const double weight = quadrature.weight( pt ) * geometry.integrationElement( x );

    RangeType vu( RangeType( 0 ) );
    uLocal.evaluate( quadrature[ pt ], vu );
    JacobianRangeType du,dun;
    uLocal.jacobian( quadrature[ pt ], du );

    // compute mass contribution (studying linear case so linearizing around zero)
    RangeType avu( RangeType( 0 ) );
    model().source( entity, quadrature[ pt ], vu, avu );
    avu *= weight;
    // add to local functional wLocal.axpy( quadrature[ pt ], avu );

    JacobianRangeType adu( 0 ), gpdu( 0 ), wpdu( 0 ), cpdu( 0 );

    model().gravity_Flux( entity, quadrature[ pt ], vu, du, gpdu );
    model().wetting_Pressure_DiffusiveFlux( entity, quadrature[ pt ], vu, du, wpdu );
    model().capillary_Pressure_DiffusiveFlux(  entity, quadrature[ pt ], vu, du, cpdu );

    adu = wpdu;
    adu += gpdu;
    adu += cpdu;
    adu *= weight;
    // add to local function
    wLocalEn.axpy( quadrature[ pt ], avu, adu );
  }
}


template< class DiscreteFunction, class Model >
void FlowOperator< DiscreteFunction, Model >
:: internal_Bnd_terms( const EntityType& entity,
                       const int& quadOrder,
                       const GeometryType &geometry,
                       const IntersectionType &intersection,
                       const DiscreteFunctionSpaceType &dfSpace,
                       const LocalFunctionType& uLocal,
                       const DiscreteFunctionType& u,
                       LocalFunctionType& wLocalEn ) const
{
  const double area = entity.geometry().volume();

  if ( intersection.neighbor() )
  {
    const EntityPointerType pOutside = intersection.outside(); // pointer to outside element.
    const EntityType &outside = pOutside;
    typedef typename IntersectionType::Geometry  IntersectionGeometryType;
    const IntersectionGeometryType &intersectionGeometry = intersection.geometry();

    // compute penalty factor
    const double intersectionArea = intersectionGeometry.volume();

    LocalFunctionType uOutLocal = u.localFunction( outside ); // local u on outisde element

    FaceQuadratureType quadInside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::INSIDE );
    FaceQuadratureType quadOutside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::OUTSIDE );
    const size_t numQuadraturePoints = quadInside.nop();
    //! [Compute skeleton terms: iterate over intersections]
    for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
    {
      //! [Compute skeleton terms: obtain required values on the intersection]
      // get coordinate of quadrature point on the reference element of the intersection
      const typename FaceQuadratureType::LocalCoordinateType &x = quadInside.localPoint( pt );
      const DomainType normal = intersection.integrationOuterNormal( x );
      const double weight = quadInside.weight( pt );

      RangeType vu( RangeType( 0 ) );
      uLocal.evaluate( quadInside[ pt ], vu );
      JacobianRangeType du,dun;
      uLocal.jacobian( quadInside[ pt ], du );

      RangeType value( RangeType( 0 ) );
      JacobianRangeType dvalue, advalue, gpdadvalue( 0 ), wpdadvalue( 0 ), cpdadvalue( 0 );
      JacobianRangeType dvalue0( 0 ), dvalue1( 0 );

      RangeType vuIn( RangeType( 0 ) ), vuOut( RangeType( 0 ) ), jump( RangeType( 0 ) );

      // RangeType vunIn( RangeType( 0 ) ),vunOut( RangeType( 0 ) );
      JacobianRangeType dunIn,dunOut;

      JacobianRangeType duIn, aduIn, duOut, aduOut;

      JacobianRangeType gpduIn( 0 ), wpduIn( 0 ),cpduIn( 0 ), gpduOut( 0 ), wpduOut( 0 ),cpduOut( 0 );

      uLocal.evaluate( quadInside[ pt ], vuIn );
      uLocal.jacobian( quadInside[ pt ], duIn );

      // const DomainType xGlobal = entity.geometry().global( coordinate(quadInside[ pt ] ) );

      double upwMob( 0 );
      double upwgradnonwetmob( 0 );
      double permeabilityIn( 0 ), permeabilityOut( 0 );
      model_.permeability_K( entity,quadInside[ pt ], permeabilityIn);
      model_.permeability_K( outside,quadOutside[ pt ], permeabilityOut);

      if ( with_upw == true  )
      {
        model_.upwindingMobility(entity,quadInside[ pt ],outside,quadOutside[ pt ],normal,vuIn,vuOut,duIn,duOut,upwMob,upwgradnonwetmob);
      }

      double betapress( 0 );
      double betasat( 0 );
      model_.penality_pressure(entity,quadInside[ pt ],outside,quadOutside[ pt ], betapress);
      model_.penality_saturation(entity,quadInside[ pt ],outside,quadOutside[ pt ], betasat);
      betapress *= intersectionArea / std::min( area, outside.geometry().volume() );
      betasat *= intersectionArea / std::min( area, outside.geometry().volume() );

      model().gravity_Flux( entity, quadInside[ pt ], vuIn, duIn, gpduIn, upwMob,with_upw );
      model().wetting_Pressure_DiffusiveFlux( entity, quadInside[ pt ], vuIn, duIn, wpduIn,upwMob,with_upw );
      model().capillary_Pressure_DiffusiveFlux( entity, quadInside[ pt ], vuIn, duIn, cpduIn , upwMob,with_upw );

      aduIn = wpduIn;
      aduIn += gpduIn;
      aduIn += cpduIn;

      uOutLocal.evaluate( quadOutside[ pt ], vuOut );
      uOutLocal.jacobian( quadOutside[ pt ], duOut );

      model().wetting_Pressure_DiffusiveFlux( entity, quadInside[ pt ], vuIn, dvalue0, wpdadvalue, upwMob, with_upw );
      model().gravity_Flux( outside, quadOutside[ pt ], vuOut, duOut, gpduOut, upwMob, with_upw );
      model().wetting_Pressure_DiffusiveFlux( outside, quadOutside[ pt ], vuOut, duOut, wpduOut, upwMob, with_upw );
      model().capillary_Pressure_DiffusiveFlux( entity, quadInside[ pt ], vuIn, dvalue1, cpdadvalue , upwMob, with_upw );
      model().capillary_Pressure_DiffusiveFlux( outside, quadOutside[ pt ], vuOut, duOut, cpduOut , upwMob, with_upw );

      advalue = wpdadvalue;
      // advalue+=cpdadvalue;

      aduOut = wpduOut;
      aduOut += gpduOut;
      aduOut += cpduOut;

      //! [Compute skeleton terms: obtain required values on the intersection]

      //! [Compute skeleton terms: compute factors for axpy method]
      // penalty term : beta [u] [phi] = beta (u+ - u-)(phi+ - phi-)=beta (u+ - u-)phi+

      RangeType cap_pressIn( RangeType( 0 ) ),cap_pressOut( RangeType( 0 ) );
      model_.cap_Press( entity, quadInside[ pt ], vuIn, cap_pressIn );
      model_.cap_Press( entity, quadOutside[ pt ], vuOut, cap_pressOut );

      jump[0] = vuIn[0] - vuOut[0];
      jump[1] = vuIn[1] - vuOut[1];

      value[0] = jump[0];
      value [1] = jump[1];//[sw]

      value[0] *= betapress * intersectionGeometry.integrationElement( x );
      value[1] *= betasat * intersectionGeometry.integrationElement( x );

      aduIn*=permeabilityOut;
      aduOut*=permeabilityIn;
      aduIn += aduOut;
      aduIn /= -(permeabilityIn+permeabilityOut);
      aduIn.umv(normal,value);
      //
      for (int r=0;r<dimRange;++r)
      {
        for (int d=0;d<dimDomain;++d)
        {
          dvalue0[0][d] = -( permeabilityOut / ( permeabilityIn+permeabilityOut ) ) *epsil* normal[d] * jump[0];
          dvalue0[1][d] = -( permeabilityOut / ( permeabilityIn+permeabilityOut ) ) *epsil* normal[d] * jump[0];

          dvalue1[0][d] = -( permeabilityOut / ( permeabilityIn+permeabilityOut ) ) *epsil* normal[d] * jump[1];
          dvalue1[1][d] = -( permeabilityOut / ( permeabilityIn+permeabilityOut ) ) *epsil* normal[d] * jump[1];
        }
      }
      value *= weight;
      advalue *= weight;
      wLocalEn.axpy( quadInside[ pt ], value, advalue );
    }
  }
  else if( intersection.boundary() && model().isDirichletIntersection( intersection ) )
  {
    typedef typename IntersectionType::Geometry  IntersectionGeometryType;
    const IntersectionGeometryType &intersectionGeometry = intersection.geometry();

    const double intersectionArea = intersectionGeometry.volume();

    FaceQuadratureType quadInside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::INSIDE );
    const size_t numQuadraturePoints = quadInside.nop();

    for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
    {
      const typename FaceQuadratureType::LocalCoordinateType &x = quadInside.localPoint( pt );
      const DomainType normal = intersection.integrationOuterNormal( x );
      const double weight = quadInside.weight( pt );

      RangeType value( RangeType( 0 ) );
      JacobianRangeType dvalue( 0 ), dvalue0( 0 ), dvalue1( 0 ),advalue( 0 ),gpdadvalue( 0 ),wpdadvalue( 0 ),cpdadvalue( 0 );

      double upwMob( 0 );
      // double upwgradnonwetmob( 0 );
      RangeType vuIn( RangeType( 0 ) ),jump( RangeType( 0 ) );
      JacobianRangeType duIn,dunIn, aduIn( 0 );
      JacobianRangeType gpduIn( 0 ), wpduIn( 0 ),cpduIn( 0 );

      uLocal.evaluate( quadInside[ pt ], vuIn );
      uLocal.jacobian( quadInside[ pt ], duIn );

      if (model().isRHS())
      {
        model().g( RangeType( 0 ), entity, quadInside.point(pt), vuIn );
      }
      else
      {
        model().gravity_Flux( entity, quadInside[ pt ], vuIn, duIn, gpduIn, upwMob );
        model().wetting_Pressure_DiffusiveFlux( entity, quadInside[ pt ], vuIn, duIn, wpduIn);
        model().capillary_Pressure_DiffusiveFlux( entity, quadInside[ pt ], vuIn, duIn, cpduIn);

        aduIn = wpduIn;
        aduIn += gpduIn;
        aduIn += cpduIn;
      }
      jump[0] = vuIn[0];
      jump[1] = vuIn[1];

      value = jump;

      double  betapress( 0 );
      double  betasat( 0 );
      model_.penality_pressure( entity, quadInside[ pt ], entity, quadInside[ pt ], betapress);
      model_.penality_saturation( entity, quadInside[ pt ], entity, quadInside[ pt ], betasat);
      betapress*= intersectionArea / area;
      betasat*= intersectionArea / area;

      value[0] *= betapress * intersectionGeometry.integrationElement( x );
      value[1] *= betasat * intersectionGeometry.integrationElement( x );

      aduIn.umv(normal,value);

      for (int r=0;r<dimRange;++r)
        for (int d=0;d<dimDomain;++d)
        {
          dvalue0[0][d] = -1 *epsil* normal[d] * jump[0];
          dvalue0[1][d] = -1 *epsil* normal[d] * jump[0];

          dvalue1[0][d] = -1 *epsil* normal[d] * jump[1];
          dvalue1[1][d] = -1 *epsil* normal[d] * jump[1];
        }

      if (model().isRHS())
      {
        model().g( RangeType( 0 ), entity, quadInside.point( pt ), vuIn );

        model_.wetting_Pressure_DiffusiveFlux( entity, quadInside[ pt ], vuIn, dvalue0, wpdadvalue, upwMob, false, true );
        model_.capillary_Pressure_DiffusiveFlux( entity, quadInside[ pt ], vuIn, dvalue1, cpdadvalue, upwMob, false, true );

        advalue=wpdadvalue;
        // advalue+=cpdadvalue;
      }
      else
      {
        //    model_.Gravity_Flux( vuIn, duIn, entity, quadInside[ pt ], vuIn, dvalue0, gpdadvalue );
        model_.wetting_Pressure_DiffusiveFlux( entity, quadInside[ pt ], vuIn, dvalue0, wpdadvalue );
        model_.capillary_Pressure_DiffusiveFlux( entity, quadInside[ pt ], vuIn, dvalue1, cpdadvalue );
        advalue=wpdadvalue;
        //      advalue+=cpdadvalue;
      }

      value *= weight;
      advalue *= weight;
      wLocalEn.axpy( quadInside[ pt ], value, advalue );
    }
  }



  else if( intersection.boundary() && model().isNeumanIntersection( intersection ) )
  {
    typedef typename IntersectionType::Geometry  IntersectionGeometryType;
    // const IntersectionGeometryType &intersectionGeometry = intersection.geometry();

    // const double intersectionArea = intersectionGeometry.volume();

    FaceQuadratureType quadInside( dfSpace.gridPart(), intersection, quadOrder, FaceQuadratureType::INSIDE );
    const size_t numQuadraturePoints = quadInside.nop();

    for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
    {
      const typename FaceQuadratureType::LocalCoordinateType &x = quadInside.localPoint( pt );
      const DomainType normal = intersection.integrationOuterNormal( x );
      double faceVol = normal.two_norm();

      const double weight = quadInside.weight( pt )*faceVol;

      RangeType value( RangeType( 0 ) );
      JacobianRangeType dvalue( 0 ), dvalue0( 0 ), dvalue1( 0 ),advalue( 0 ),wpdadvalue( 0 ),cpdadvalue( 0 );

      RangeType vuIn( RangeType( 0 ) ), jump(RangeType( 0 ) );
      JacobianRangeType duIn,dunIn, aduIn( 0 );
      JacobianRangeType wpduIn( 0 ),cpduIn( 0 );

      uLocal.jacobian( quadInside[ pt ], duIn );

      if (model().isRHS())
      {
        model().gN( RangeType( 0 ), entity, quadInside.point( pt ), vuIn );
      }
      else
      {
        vuIn = 0;
      }

      jump[0] = -1* vuIn[0];
      jump[1] = -1*vuIn[1];

      value = jump;

      value *= weight;
      advalue = 0;
      wLocalEn.axpy( quadInside[ pt ], value, advalue );
    }
  }

  else
  {
  }
}


template< class JacobianOperator, class Model >
void DifferentiableFlowOperator< JacobianOperator, Model >
::element_Local_Contribution(const EntityType& entity,
                             const GeometryType &geometry,
                             const DiscreteFunctionSpaceType &dfSpace,
                             const BasisFunctionSetType &baseSet,
                             const LocalFunctionType& uLocal,
                             LocalMatrixType& jLocal) const
{
  QuadratureType quadrature( entity, 2*dfSpace.order()+1 );
  const int blockSize = dfSpace.localBlockSize; // is equal to 1 for scalar functions
  std::vector< typename LocalFunctionType::RangeType > phi( dfSpace.blockMapper().maxNumDofs()*blockSize );
  std::vector< typename LocalFunctionType::JacobianRangeType > dphi( dfSpace.blockMapper().maxNumDofs()*blockSize );
  const unsigned int numBasisFunctions = baseSet.size();

  const size_t numQuadPoints = quadrature.nop();
  for( size_t pt = 0; pt < numQuadPoints; ++pt )
  {

    const typename QuadratureType::CoordinateType &x = quadrature.point( pt );
    const double weight = quadrature.weight( pt ) * geometry.integrationElement( x );

    // evaluate all basis functions at given quadrature point
    baseSet.evaluateAll( quadrature[ pt ], phi );

    // evaluate jacobians of all basis functions at given quadrature point
    baseSet.jacobianAll( quadrature[ pt ], dphi );

    //! [model for non-linear time integration]
    // get value for linearization and for u^n
    RangeType u0( RangeType( 0 ) );
    uLocal.evaluate( quadrature[ pt ], u0 );
    JacobianRangeType jacU0;
    uLocal.jacobian( quadrature[ pt ], jacU0 );

    RangeType aphi( RangeType( 0 ) );
    JacobianRangeType adphi( 0 ), lgpdphi( 0 ),lwpdphi( 0 ),cpdphi( 0 );
    for( unsigned int localCol = 0; localCol < numBasisFunctions; ++localCol )
    {
      // if mass terms or right hand side is present
      model().linSource( u0, entity, quadrature[ pt ], phi[ localCol ], aphi );
      // if gradient term is present
      model().linGravity_Flux( u0, jacU0, entity, quadrature[ pt ], phi[ localCol ], dphi[ localCol ], lgpdphi );
      model().linWetting_Pressure_DiffusiveFlux( u0, jacU0, entity, quadrature[ pt ], phi[ localCol ], dphi[ localCol ], lwpdphi );
      model().linCapillary_Pressure_DiffusiveFlux( u0, jacU0, entity, quadrature[ pt ], phi[ localCol ], dphi[ localCol ], cpdphi );
      adphi=lwpdphi;
      adphi+=lgpdphi;
      adphi+=cpdphi;
      // get column object and call axpy method
      jLocal.column( localCol ).axpy( phi, dphi, aphi, adphi, weight );
    }
  }
}



template< class JacobianOperator, class Model >
void DifferentiableFlowOperator< JacobianOperator, Model >
::interface_Local_Contribution(const EntityType& entity,
                               const GeometryType &geometry,
                               const DiscreteFunctionSpaceType &dfSpace,
                               const BasisFunctionSetType &baseSet,
                               const IntersectionType& intersection,
                               const LocalFunctionType& uLocal,
                               const DiscreteFunctionType& u,
                               const JacobianOperatorType &jOp,
                               LocalMatrixType& jLocal) const
{
  double area = geometry.volume();

  const GridPartType& gridPart = dfSpace.gridPart();
  const int blockSize = dfSpace.localBlockSize; // is equal to 1 for scalar functions
  std::vector< typename LocalFunctionType::RangeType > phi( dfSpace.blockMapper().maxNumDofs()*blockSize );
  std::vector< typename LocalFunctionType::JacobianRangeType > dphi( dfSpace.blockMapper().maxNumDofs()*blockSize );
  std::vector< typename LocalFunctionType::RangeType > phiNb( dfSpace.blockMapper().maxNumDofs()*blockSize );
  std::vector< typename LocalFunctionType::JacobianRangeType > dphiNb( dfSpace.blockMapper().maxNumDofs()*blockSize );

  const unsigned int numBasisFunctions = baseSet.size();


  if( intersection.neighbor() )
  {
    EntityPointerType ep = intersection.outside();
    const EntityType& neighbor = ep;
    typedef typename IntersectionType::Geometry  IntersectionGeometryType;
    const IntersectionGeometryType &intersectionGeometry = intersection.geometry();

    const LocalFunctionType uOutLocal = u.localFunction( neighbor );

    //! [Assemble skeleton terms: get contributions on off diagonal block]
    // get local matrix for face entries
    LocalMatrixType localOpNb = jOp.localMatrix( neighbor, entity );
    // get neighbor's base function set
    const BasisFunctionSetType &baseSetNb = localOpNb.domainBasisFunctionSet();
    //! [Assemble skeleton terms: get contributions on off diagonal block]

    // compute penalty factor
    const double intersectionArea = intersectionGeometry.volume();

    // here we assume that the intersection is conforming
    FaceQuadratureType faceQuadInside(gridPart, intersection, 2*dfSpace.order() + 1,
              FaceQuadratureType::INSIDE);
    FaceQuadratureType faceQuadOutside(gridPart, intersection, 2*dfSpace.order() + 1,
               FaceQuadratureType::OUTSIDE);

    const size_t numFaceQuadPoints = faceQuadInside.nop();
    for( size_t pt = 0; pt < numFaceQuadPoints; ++pt )
    {
      const typename FaceQuadratureType::LocalCoordinateType &x = faceQuadInside.localPoint( pt );
      DomainType normal = intersection.integrationOuterNormal( x );
      double faceVol = normal.two_norm();
      normal /= faceVol; // make it into a unit normal

      const double quadWeight = faceQuadInside.weight( pt );
      const double weight = quadWeight * faceVol;

      //! [Assemble skeleton terms: obtain values om quadrature point]
      RangeType u0En( RangeType( 0 ) );
      JacobianRangeType u0EnJac;
      uLocal.evaluate( faceQuadInside[ pt ], u0En );
      uLocal.jacobian( faceQuadInside[ pt ], u0EnJac );

      RangeType u0OutEn( RangeType( 0 ) );
      JacobianRangeType u0OutEnJac;
      uOutLocal.evaluate( faceQuadOutside[ pt ], u0OutEn );
      uOutLocal.jacobian( faceQuadOutside[ pt ], u0OutEnJac );

      /////////////////////////////////////////////////////////////
      // evaluate basis function of face inside E^- (entity)
      /////////////////////////////////////////////////////////////

      // evaluate all basis functions for quadrature point pt
      baseSet.evaluateAll( faceQuadInside[ pt ], phi );

      // evaluate the jacobians of all basis functions
      baseSet.jacobianAll( faceQuadInside[ pt ], dphi );

      /////////////////////////////////////////////////////////////
      // evaluate basis function of face inside E^+ (neighbor)
      /////////////////////////////////////////////////////////////

      // evaluate all basis functions for quadrature point pt on neighbor
      baseSetNb.evaluateAll( faceQuadOutside[ pt ], phiNb );

      // evaluate the jacobians of all basis functions on neighbor
      baseSetNb.jacobianAll( faceQuadOutside[ pt ], dphiNb );

      // RangeType cap_pressEn( RangeType( 0 ) ), cap_pressNb( RangeType( 0 ) );
      JacobianRangeType dvalue0En( 0 ), dvalue0Nb( 0 ), dvalue1En( 0 ), dvalue1Nb( 0 ), advalueEn( 0 ),advalueNb( 0 );

      for( unsigned int i = 0; i < numBasisFunctions; ++i )
      {
        JacobianRangeType lgpdphiEn( 0 ), lwpdphiEn( 0 ),cpdphiEn( 0 ), lgpdphiNb( 0 ), lwpdphiNb( 0 ),cpdphiNb( 0 ),gpdadvalueEn( 0 ), wpdadvalueEn( 0 ), cpdadvalueEn( 0 ), gpdadvalueNb( 0 ), wpdadvalueNb( 0 ), cpdadvalueNb( 0 );
        JacobianRangeType adphiEn = dphi[ i ];
        JacobianRangeType adphiNb = dphiNb[ i ];
        double upwMob( 0 );
        double upwgradnonwetmob( 0 );

        double permeabilityIn, permeabilityOut;
        model().permeability_K( entity, faceQuadInside[ pt ], permeabilityIn );
        model().permeability_K( neighbor, faceQuadOutside[ pt ], permeabilityOut );

        // const DomainType xGlobal = entity.geometry().global( coordinate( faceQuadInside[ pt ] ) );

        if ( with_upw == true )
        {
          model().upwindingMobility( entity, faceQuadInside[ pt ], neighbor, faceQuadOutside[ pt ], normal, u0En, u0OutEn, u0EnJac, u0OutEnJac, upwMob, upwgradnonwetmob );
        }

        model().linGravity_Flux( u0En, u0EnJac, entity,   faceQuadInside[ pt ], phi[i], adphiEn, lgpdphiEn, upwMob, upwgradnonwetmob, with_upw );
        model().linWetting_Pressure_DiffusiveFlux( u0En, u0EnJac, entity, faceQuadInside[ pt ], phi[i], adphiEn, lwpdphiEn, upwMob, upwgradnonwetmob, with_upw );
        model().linCapillary_Pressure_DiffusiveFlux( u0En, u0EnJac, entity, faceQuadInside[ pt ], phi[i], adphiEn, cpdphiEn);

        dphi[ i ] = lwpdphiEn;
        dphi[ i ] += lgpdphiEn;
        dphi[ i ] += cpdphiEn;

        model().linGravity_Flux(u0OutEn, u0OutEnJac, neighbor, faceQuadOutside[ pt ], phiNb[i], adphiNb, lgpdphiNb , upwMob, upwgradnonwetmob, with_upw );
        model().linWetting_Pressure_DiffusiveFlux( u0OutEn, u0OutEnJac,  neighbor, faceQuadOutside[ pt ], phiNb[i], adphiNb, lwpdphiNb , upwMob, upwgradnonwetmob, with_upw );
        model().linWetting_Pressure_DiffusiveFlux( u0En, u0EnJac, entity, faceQuadInside[ pt ], phi[i], dvalue0En, wpdadvalueEn, upwMob, upwgradnonwetmob, with_upw );
        model().linWetting_Pressure_DiffusiveFlux( u0OutEn, u0OutEnJac, neighbor, faceQuadOutside[ pt ], phiNb[i], dvalue0Nb, wpdadvalueNb, upwMob, upwgradnonwetmob, with_upw );
        model().linCapillary_Pressure_DiffusiveFlux(  u0OutEn, u0OutEnJac,neighbor, faceQuadOutside[ pt ], phiNb[i], adphiNb, cpdphiNb );
        model().linCapillary_Pressure_DiffusiveFlux(  u0OutEn, u0OutEnJac,  neighbor, faceQuadOutside[ pt ], phiNb[i], dvalue1Nb, cpdadvalueNb );
        model().linCapillary_Pressure_DiffusiveFlux( u0En, u0EnJac,  entity,   faceQuadInside[ pt ], phi[i], dvalue1En, cpdadvalueEn );

        dphiNb[ i ] = lwpdphiNb;
        dphiNb[ i ] += lgpdphiNb;
        dphiNb[ i ] += cpdphiNb;

        // model().linGravity_Flux( u0En, u0EnJac, u0En, u0EnJac, entity,   faceQuadInside[ pt ], phi[i], dvalue0En, gpdadvalueEn);
        // model().linGravity_Flux(  u0OutEn, u0OutEnJac, u0En, u0EnJac, neighbor,   faceQuadOutside[ pt ], phiNb[i], dvalue0Nb, gpdadvalueNb);

        advalueEn = wpdadvalueEn;
        //   advalueEn += cpdadvalueEn;
        advalueNb = wpdadvalueNb;
        //   advalueNb += cpdadvalueNb;
      }

      //! [Assemble skeleton terms: compute factors for axpy method]
      for( unsigned int localCol = 0; localCol < numBasisFunctions; ++localCol )
      {
        RangeType valueEn( RangeType( 0 ) ), valueNb( RangeType( 0 ) );
        JacobianRangeType dvalueEn( 0 ), dvalueNb( 0 );
        double Kin( 0 ),Kout( 0 );
        model().permeability_K( entity, faceQuadInside[ pt ], Kin);
        model().permeability_K( neighbor, faceQuadOutside[ pt ], Kout);

        //  -{ A grad u } * [ phi_en ]
        dphi[localCol].usmv( -( Kout / ( Kin + Kout ) ), normal, valueEn );

        //  -{ A grad u } * [ phi_en ]
        dphiNb[localCol].usmv( -( Kin /( Kin + Kout ) ), normal, valueNb );

        double  betapress( 0 );
        double  betasat( 0 );
        model().penality_pressure( entity, faceQuadInside[ pt ], neighbor, faceQuadOutside[ pt ], betapress );
        model().penality_saturation( entity, faceQuadInside[ pt ], neighbor, faceQuadOutside[ pt ], betasat );

        // const DomainType xGlobal = entity.geometry().global( coordinate( faceQuadInside[ pt ] ) );

        betapress*= intersectionArea / std::min( area, neighbor.geometry().volume() );
        betasat*= intersectionArea / std::min( area, neighbor.geometry().volume() );

        valueEn[0] += betapress*phi[localCol][0];
        valueEn[1] += betasat*phi[localCol][1];
        valueNb[0] -= betapress*phiNb[localCol][0];
        valueNb[1] -= betasat*phiNb[localCol][1];

        // here we need a diadic product of u x n
        for ( int r = 0; r < dimRange; ++r )
          for ( int d = 0; d < dimDomain; ++d )
          {
            dvalue0En[r][d] = -(Kout/(Kin+Kout)) *epsil* normal[d] * phi[localCol][0];
            dvalue1En[r][d] = -(Kout/(Kin+Kout)) *epsil* normal[d] * phi[localCol][1];
            //
            dvalue0Nb[r][d] = (Kin/(Kin+Kout)) *epsil* normal[d] * phiNb[localCol][0];
            dvalue1Nb[r][d] = (Kin/(Kin+Kout)) *epsil* normal[d] * phiNb[localCol][1];
          }

        jLocal.column( localCol ).axpy( phi, dphi, valueEn, advalueEn, weight );
        localOpNb.column( localCol ).axpy( phi, dphi, valueNb, advalueNb, weight );
      }
      //! [Assemble skeleton terms: compute factors for axpy method]
    }
  }
  else if( intersection.boundary() )
  {
    if ( ! model().isDirichletIntersection( intersection ) )
    {
    }

    else
    {
      typedef typename IntersectionType::Geometry  IntersectionGeometryType;
      const IntersectionGeometryType &intersectionGeometry = intersection.geometry();

      // compute penalty factor
      const double intersectionArea = intersectionGeometry.volume();

      // here we assume that the intersection is conforming
      FaceQuadratureType faceQuadInside(gridPart, intersection, 2*dfSpace.order() + 1,
                FaceQuadratureType::INSIDE);

      const size_t numFaceQuadPoints = faceQuadInside.nop();
      for( size_t pt = 0; pt < numFaceQuadPoints; ++pt )
      {
        const typename FaceQuadratureType::LocalCoordinateType &x = faceQuadInside.localPoint( pt );
        DomainType normal = intersection.integrationOuterNormal( x );
        double faceVol = normal.two_norm();
        normal /= faceVol; // make it into a unit normal

        const double quadWeight = faceQuadInside.weight( pt );
        const double weight = quadWeight * faceVol;

        RangeType u0En( RangeType( 0 ) );
        JacobianRangeType u0EnJac;
        uLocal.evaluate( faceQuadInside[ pt ], u0En );
        uLocal.jacobian( faceQuadInside[ pt ], u0EnJac );

        /////////////////////////////////////////////////////////////
        // evaluate basis function of face inside E^- (entity)
        /////////////////////////////////////////////////////////////

        // evaluate all basis functions for quadrature point pt
        baseSet.evaluateAll( faceQuadInside[ pt ], phi );

        // evaluate the jacobians of all basis functions
        baseSet.jacobianAll( faceQuadInside[ pt ], dphi );

        // RangeType cap_pressEn( RangeType( 0 ) );
        JacobianRangeType dvalue0En( 0 ),dvalue1En( 0 ), advalueEn( 0 );

          for( unsigned int i = 0; i < numBasisFunctions; ++i )
          {
          JacobianRangeType adphiEn = dphi[ i ];
          JacobianRangeType lgpdphiEn( 0 ), lwpdphiEn( 0 ),cpdphiEn( 0 );
          JacobianRangeType gpdadvalueEn( 0 ), wpdadvalueEn( 0 ), cpdadvalueEn( 0 );

          model().linGravity_Flux( u0En, u0EnJac, entity,   faceQuadInside[ pt ], phi[i], adphiEn, lgpdphiEn );

          model().linWetting_Pressure_DiffusiveFlux( u0En, u0EnJac, entity,   faceQuadInside[ pt ], phi[i], adphiEn, lwpdphiEn );
          model().linCapillary_Pressure_DiffusiveFlux( u0En, u0EnJac, entity,   faceQuadInside[ pt ], phi[i], adphiEn, cpdphiEn );

          dphi[ i ] = lwpdphiEn;
          dphi[ i ] += lgpdphiEn;
          dphi[ i ] += cpdphiEn;
          model().linWetting_Pressure_DiffusiveFlux( u0En, u0EnJac, entity,   faceQuadInside[ pt ], phi[i], dvalue0En, wpdadvalueEn );
          model().linCapillary_Pressure_DiffusiveFlux( u0En, u0EnJac, entity,   faceQuadInside[ pt ], phi[i], dvalue1En, cpdadvalueEn );

          advalueEn = wpdadvalueEn;
          //                    advalueEn+=cpdadvalueEn;
        }

        for( unsigned int localCol = 0; localCol < numBasisFunctions; ++localCol )
        {
          RangeType valueEn( RangeType( 0 ) );
          JacobianRangeType dvalueEn( 0 );

          dphi[localCol].usmv( -1.0 , normal, valueEn );

          double  betapress( 0 );
          double  betasat( 0 );
          model().penality_pressure(entity,faceQuadInside[ pt ],entity,faceQuadInside[ pt ], betapress);
          model().penality_saturation(entity,faceQuadInside[ pt ],entity,faceQuadInside[ pt ], betasat);
          betapress *=  intersectionArea / area;
          betasat *=  intersectionArea / area;

          valueEn[0] += betapress*phi[localCol][0];
          valueEn[1] += betasat*phi[localCol][1];

          for ( int r=0; r< dimRange; ++r )
            for ( int d=0; d< dimDomain; ++d )
            {
              dvalue0En[r][d] = - 1.0 * epsil * normal[d] * phi[localCol][0];
              dvalue1En[r][d] = - 1.0 * epsil * normal[d] * phi[localCol][1];
            }

          jLocal.column( localCol ).axpy( phi, dphi, valueEn, advalueEn, weight );
        }
      }

    }
  }

}


// Implementation of FlowOperator
// ----------------------------------

template< class DiscreteFunction, class Model >
void FlowOperator< DiscreteFunction, Model >
::operator() ( const DiscreteFunctionType &u, DiscreteFunctionType &w ) const
{
  // clear destination
  w.clear();

  // get discrete function space
  const DiscreteFunctionSpaceType &dfSpace = w.space();

  // iterate over grid
  const IteratorType end = dfSpace.end();
  for( IteratorType it = dfSpace.begin(); it != end; ++it )
  {
    // get entity (here element)
    const EntityType &entity = *it;
    // get elements geometry
    const GeometryType &geometry = entity.geometry();

    // get local representation of the discrete functions
    const LocalFunctionType uLocal = u.localFunction( entity );
    LocalFunctionType wLocal = w.localFunction( entity );

    // obtain quadrature order
    const int quadOrder = uLocal.order() + wLocal.order();

    //Computing local contibution from elements
    volumetricPart( entity, quadOrder, geometry, uLocal, wLocal );

    if ( ! dfSpace.continuous() )
    {
      //const double area = entity.geometry().volume();
      const IntersectionIteratorType iitend = dfSpace.gridPart().iend( entity );
      for( IntersectionIteratorType iit = dfSpace.gridPart().ibegin( entity ); iit != iitend; ++iit ) // looping over intersections
      {
        const IntersectionType &intersection = *iit;
        //Computing local contibution from interfaces and boundaries
        internal_Bnd_terms( entity, quadOrder, geometry, intersection, dfSpace, uLocal, u, wLocal );
      }
    }
  }

  // communicate data (in parallel runs)
  w.communicate();
}

// Implementation of DifferentiableFlowOperator
// ------------------------------------------------

template< class JacobianOperator, class Model >
void DifferentiableFlowOperator< JacobianOperator, Model >
::jacobian ( const DiscreteFunctionType &u, JacobianOperator &jOp ) const
{
  typedef typename JacobianOperator::LocalMatrixType LocalMatrixType;
  typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType BasisFunctionSetType;

  const DiscreteFunctionSpaceType &dfSpace = u.space();
  const GridPartType& gridPart = dfSpace.gridPart();
  Dune::Fem::DiagonalAndNeighborStencil<DiscreteFunctionSpaceType,DiscreteFunctionSpaceType> stencil( dfSpace, dfSpace );
  jOp.reserve(stencil);
  jOp.clear();

  const IteratorType end = dfSpace.end();
  for( IteratorType it = dfSpace.begin(); it != end; ++it )
  {
    const EntityType &entity = *it;
    const GeometryType &geometry = entity.geometry();

    const LocalFunctionType uLocal = u.localFunction( entity );

    LocalMatrixType jLocal = jOp.localMatrix( entity, entity );

    const BasisFunctionSetType &baseSet = jLocal.domainBasisFunctionSet();

    element_Local_Contribution( entity, geometry, dfSpace, baseSet, uLocal, jLocal );

    if ( dfSpace.continuous() )
      continue;

    const IntersectionIteratorType endiit = gridPart.iend( entity );
    for ( IntersectionIteratorType iit = gridPart.ibegin( entity ); iit != endiit ; ++ iit )
    {
      const IntersectionType& intersection = *iit ;
      //Computing local contibution from interfaces and boundaries
      interface_Local_Contribution( entity, geometry, dfSpace, baseSet, intersection, uLocal, u, jOp, jLocal );
    } //! [Assembling the local matrix]
  } // end grid traversal

  jOp.communicate();
}
#endif // #ifndef OPERATOR_HH
