#ifndef __NON_SMOOTH_MINIMIZATION_PHC_HH__
#define __NON_SMOOTH_MINIMIZATION_PHC_HH__

// include header of adaptive scheme
#include <dune/acfem/common/discretefunctionselector.hh>
#include <dune/acfem/functions/basicfunctions.hh>
#include <dune/acfem/models/basicmodels.hh>
#include <dune/acfem/models/modelexpression.hh>
#include <dune/acfem/algorithms/ellipticfemscheme.hh>
#include <dune/acfem/operators/l2projection.hh>

using namespace Dune::ACFem ;

//the projection class that does step one of the algorithm
// we now do this weakly.
// p = 1/ max{ |..|, 1} * (p + sigma * nabla u)
// with testfunctions in H^1_0
// so dirichlet 0 boundary
template<class ForwardDiscreteFunctionType, class AdjointDiscreteFunctionType, bool pContinuous, bool uContinuous>
class ProjectionHelperClass
{
public:
  //constructor
  ProjectionHelperClass(const ForwardDiscreteFunctionType& projG, const double sigma, const double tau, const double lambda_2 )
    : projG_(projG), sigma_(sigma), tau_(tau), lambda_2_(lambda_2){};

  void entitywiseProjection ( const ForwardDiscreteFunctionType & uBar, AdjointDiscreteFunctionType & p )
  {
    //the gradient model defines the weak gradient using continuous test-functions
    //from the space of p
    auto gradU = gradientModel(uBar);
    auto Dbc0 = dirichletZeroModel(p);
    auto mass = massModel(p);
    //we solve p = sigma * \nabla u + p
    //with dirichlet zero boundary
    auto model = mass - sigma_ * gradU - p + Dbc0 ;
    //Testspace = AdjointDiscreteFunctionSpaceType = space of p
    typedef EllipticFemScheme<AdjointDiscreteFunctionType, decltype(model)> SchemeType;
    //p is the returned solution
    SchemeType scheme(p, model);

    scheme.solve();

    typedef typename AdjointDiscreteFunctionType::GridPartType GridPartType;
    typedef typename GridPartType::template Codim< 0 >::IteratorType IteratorType;
    typedef typename IteratorType::Entity EntityType;

    typedef typename AdjointDiscreteFunctionType::LocalFunctionType AdjointLocalFunctionType;

    typedef typename AdjointDiscreteFunctionType::RangeType AdjointRangeType;
    const int dimDomain = AdjointLocalFunctionType::dimDomain;

    const GridPartType & gridPart = p.gridPart();

    const IteratorType end = gridPart.template end<0>();

    for(IteratorType it = gridPart.template begin<0>(); it != end; ++it)
    {
      const EntityType & entity = *it;
      AdjointLocalFunctionType pLocal = p.localFunction(entity);
      for(int i = 0; i < pLocal.numScalarDofs(); ++i)
      {
        AdjointRangeType result;
        for(int j = 0; j < dimDomain ; ++j)
        {
          result[j] = pLocal[dimDomain * i+j];
        }
        if(result.two_norm() > 1.)
        {
          //pLocal[i] accesses the dof vector - it is not an evaluation
          for(int j = 0; j < dimDomain; ++j)
          {
            pLocal[dimDomain* i+j] /= result.two_norm();
          }
        }
      }
    }
  }


  void calculateZ (const AdjointDiscreteFunctionType &p, const ForwardDiscreteFunctionType & uOld, ForwardDiscreteFunctionType& z)
  {
    auto divP  = divergence(p);
    L2Projection(1./(1.+tau_*lambda_2_) * ( tau_ * divP + uOld + tau_ * lambda_2_ * projG_), z);
  }

private:
  const ForwardDiscreteFunctionType & projG_;
  const double sigma_,tau_,lambda_2_;
};

//specialisation of PHC if the adjoint is not contiuous.
// we now assume u has a gradient
// (if it is pw constant the gradient will evaluate to zero)
template<class ForwardDiscreteFunctionType, class AdjointDiscreteFunctionType>
class ProjectionHelperClass< ForwardDiscreteFunctionType, AdjointDiscreteFunctionType, false, true>
{
public:
  //constructor
  ProjectionHelperClass(const ForwardDiscreteFunctionType& projG, const double sigma, const double tau, const double lambda_2 )
    : projG_(projG), sigma_(sigma), tau_(tau), lambda_2_(lambda_2){};

  void entitywiseProjection ( const ForwardDiscreteFunctionType & uBar, AdjointDiscreteFunctionType & p )
  {
    // p is not continuous - so u has to be at least of order 1
    //so we have a function and basically need an L2 projection onto
    // th space of p with p= 0 or p*n =0 on the boundary
    //possible solution:
    // l2 projection and then project

    auto gradU = gradient(uBar);

    L2Projection(sigma_ * gradU +p, p);


    typedef typename AdjointDiscreteFunctionType::GridPartType GridPartType;
    typedef typename GridPartType::template Codim< 0 >::IteratorType IteratorType;
    typedef typename IteratorType::Entity EntityType;

    typedef typename AdjointDiscreteFunctionType::LocalFunctionType AdjointLocalFunctionType;

    typedef typename AdjointDiscreteFunctionType::RangeType AdjointRangeType;
    const int dimDomain = AdjointLocalFunctionType::dimDomain;


    const GridPartType & gridPart = p.gridPart();

    const IteratorType end = gridPart.template end<0>();

    for(IteratorType it = gridPart.template begin<0>(); it != end; ++it)
    {
      const EntityType & entity = *it;
      AdjointLocalFunctionType pLocal = p.localFunction(entity);

      for(int i = 0; i < pLocal.numScalarDofs(); ++i)
      {
        AdjointRangeType result;
        for(int j = 0; j < dimDomain ; ++j)
        {
          result[j] = pLocal[dimDomain * i+j];
        }
        //DG dofs have to be scaled by sqrt(2)
        if(result.two_norm() > 1./std::sqrt(2.))
        {
          //pLocal[i] accesses the dof vector - it is not an evaluation
          for(int j = 0; j < dimDomain; ++j)
          {
            //DG dofs have to be scaled by sqrt(2)
            pLocal[dimDomain * i+j] /= std::sqrt(2.) * result.two_norm();
          }
        }
      }

      //the condition that p does not have outflow
      // p * n = 0
      //this is done by replacing p on boundary elements by
      // p_new = p - (p*n) n
      if(entity.hasBoundaryIntersections())
      {
        for(auto iit = gridPart.ibegin(entity); iit!=gridPart.iend(entity); ++ iit)
        {
           const auto intersection = *iit;
           if (intersection.boundary())
           {
            auto normal = intersection.centerUnitOuterNormal();
            double scalar = pLocal[0] *normal[0];
            for(int k = 1; k < dimDomain; ++k)
              scalar += pLocal[k] * normal[k];
            for(int k = 0; k < dimDomain; ++k)
              pLocal[k] -= scalar * normal[k];
           }
        }
      }
    }
  }

  void calculateZ (const AdjointDiscreteFunctionType &p, const ForwardDiscreteFunctionType & uOld, ForwardDiscreteFunctionType& z)
  {
    //get the mass model of the Forward space
    auto U_Phi = massModel(projG_.space());

    //calculate z
    auto weakDiv_P = weakDivergenceModel(p);
    auto projModel = U_Phi -1./(1.+tau_*lambda_2_) * ( tau_ * weakDiv_P + uOld + tau_ * lambda_2_ * projG_);

    typedef
    EllipticFemScheme<ForwardDiscreteFunctionType, decltype(projModel)>
    SchemeType;

    SchemeType scheme(z, projModel);

    scheme.solve();
  }

private:
  const ForwardDiscreteFunctionType & projG_;
  const double sigma_,tau_,lambda_2_;
};

//specialisation of PHC if u is not contiuous.
// we now assume p has a divergence
// (if it is pw constant the gradient will evaluate to zero)
template<class ForwardDiscreteFunctionType, class AdjointDiscreteFunctionType>
class ProjectionHelperClass< ForwardDiscreteFunctionType, AdjointDiscreteFunctionType, true, false>
{
public:
  //constructor
  ProjectionHelperClass(const ForwardDiscreteFunctionType& projG, const double sigma, const double tau, const double lambda_2 )
    : projG_(projG), sigma_(sigma), tau_(tau), lambda_2_(lambda_2){};

  void entitywiseProjection ( const ForwardDiscreteFunctionType & uBar, AdjointDiscreteFunctionType & p )
  {
    //the gradient model defines the weak gradient using continuous test-functions
    //from the space of p
    auto gradU = gradientModel(uBar);
    auto Dbc0 = dirichletZeroModel(p);
    auto mass = massModel(p);
    //we solve p = sigma * \nabla u + p
    //with dirichlet zero boundary
    auto model = mass - sigma_ * gradU - p + Dbc0 ;
    //Testspace = AdjointDiscreteFunctionSpaceType = space of p
    typedef EllipticFemScheme<AdjointDiscreteFunctionType, decltype(model)> SchemeType;
    //p is the returned solution
    SchemeType scheme(p, model);

    scheme.solve();

    //now project p to be feasible, i.e. |p(x)|_2 \leq 1 \forall x \in \Omega

    typedef typename AdjointDiscreteFunctionType::GridPartType GridPartType;
    typedef typename GridPartType::template Codim< 0 >::IteratorType IteratorType;
    typedef typename IteratorType::Entity EntityType;

    typedef typename AdjointDiscreteFunctionType::LocalFunctionType AdjointLocalFunctionType;

    typedef typename AdjointDiscreteFunctionType::RangeType AdjointRangeType;
    const int dimDomain = AdjointLocalFunctionType::dimDomain;

    const GridPartType & gridPart = p.gridPart();

    const IteratorType end = gridPart.template end<0>();
    //iterate over all elements
    for(IteratorType it = gridPart.template begin<0>(); it != end; ++it)
    {
      const EntityType & entity = *it;
      AdjointLocalFunctionType pLocal = p.localFunction(entity);
      //iterate over all dofs
      for(int i = 0; i < pLocal.numScalarDofs(); ++i)
      {
        AdjointRangeType result;
        for(int j = 0; j < dimDomain ; ++j)
        {
          result[j] = pLocal[dimDomain * i+j];
        }
        //if |p(x)|_2> 1, then restrict to 1
        if(result.two_norm() > 1)
        {
          //pLocal[i] accesses the dof vector - it is not an evaluation
          for(int j = 0; j < dimDomain; ++j)
          {
            pLocal[dimDomain* i+j] /= result.two_norm();
          }
        }
      }
    }
  }


  void calculateZ (const AdjointDiscreteFunctionType &p, const ForwardDiscreteFunctionType & uOld, ForwardDiscreteFunctionType& z)
  {
    auto divP  = divergence(p);
    L2Projection(1./(1.+tau_*lambda_2_) * ( tau_ * divP + uOld + tau_ * lambda_2_ * projG_), z);
  }

private:
  const ForwardDiscreteFunctionType & projG_;
  const double sigma_,tau_,lambda_2_;
};

/* Specialization of PHC if neither u nor p is continuous
    This just throws an error. We cannot handle this correctly as 
    we need one of the spaces to be continuous. This may change,
    once jumps are implemented in ACFem
*/
template<class ForwardDiscreteFunctionType, class AdjointDiscreteFunctionType>
class ProjectionHelperClass< ForwardDiscreteFunctionType, AdjointDiscreteFunctionType, false, false>
{
public:
  //constructor
  ProjectionHelperClass(const ForwardDiscreteFunctionType& projG, const double sigma, const double tau, const double lambda_2 )
    : projG_(projG), sigma_(sigma), tau_(tau), lambda_2_(lambda_2){ std::cerr<< "\n \n Either p or u should be continuous!!!! \n\n terminating..." << std::endl; exit(1); };

  void entitywiseProjection ( const ForwardDiscreteFunctionType & uBar, AdjointDiscreteFunctionType & p )
  {}

  void calculateZ (const AdjointDiscreteFunctionType &p, const ForwardDiscreteFunctionType & u, ForwardDiscreteFunctionType& z)
  {}

private:
  const ForwardDiscreteFunctionType & projG_;
  const double sigma_,tau_,lambda_2_;
};



#endif //__NON_SMOOTH_MINIMIZATION_PHC_HH__
