/**************************************************************************
**  demo program performing a l2-projection on Lagrange or DG-spaces
**  only SGRID supported
**  Andreas Dedner, Bernard Haasdonk 30.5.2006
**************************************************************************/

const int POLORDER = 1;

#ifndef DIM_OF_WORLD 
static const int dimw = 2;
#else
static const int dimw = DIM_OF_WORLD;
#endif

#ifndef DIM
static const int dimp = 2;
#else
static const int dimp = DIM;
#endif

#include <iostream>
#include <dune/config.h>

#define HAVE_ALBERTA 0
#include <dune/common/stdstreams.cc>

using namespace Dune;

#include <dune/grid/sgrid.hh>
typedef SGrid  < dimp, dimw > GridType;

// include continuous lagrange functions
#include <dune/fem/space/lagrangespace/lagrange.hh>
// include discontinuous lagrange functions
#include <dune/fem/space/dgspace.hh>

#include <dune/fem/discretefunction/dfadapt.hh>
#include <dune/fem/quadrature/elementquadrature.hh>

#include <dune/grid/common/defaultindexsets.hh>
#include <dune/grid/common/gridpart.hh>
#include <dune/grid/common/referenceelements.hh>

// #if HAVE_GRAPE
#include <dune/grid/io/visual/grapedatadisplay.hh>
// #endif

// polynom approximation order of quadratures, 
// at least polynom order of basis functions 
const int polOrd = POLORDER;

//! the index set we are using 
typedef LeafGridPart < GridType > GridPartType; 

//! define the function space, \f[ \R^2 \rightarrow \R \f]
// see dune/common/functionspace.hh
typedef FunctionSpace < double , double, dimp , 1 > FunctionSpaceType;

//! define the function space our unkown belong to 
//typedef DiscontinuousGalerkinSpace<FunctionSpaceType, GridPartType, 
//	polOrd,CachingStorage> DiscreteFunctionSpaceType;

typedef LagrangeDiscreteFunctionSpace 
< FunctionSpaceType, GridPartType, polOrd >        DiscreteFunctionSpaceType;

//! define the type of discrete function we are using , see
typedef DFAdapt < DiscreteFunctionSpaceType > DiscreteFunctionType;  

//! the exact solution to the problem for EOC calculation 
class ExactSolution : public Function < FunctionSpaceType , ExactSolution > 
{
  typedef FunctionSpaceType::RangeType RangeType;
  typedef FunctionSpaceType::RangeFieldType RangeFieldType;
  typedef FunctionSpaceType::DomainType DomainType;
public:
  ExactSolution (FunctionSpaceType &f) : Function < FunctionSpaceType , ExactSolution > ( f ) {}
 
  void evaluate (const DomainType & x , RangeType & ret)  const
        {
          ret = (4096/9.); // maximum of function is 1
          for(int i=0; i<DomainType::dimension; i++)
              ret *= x[i]*(1.0 -x[i])*(0.5-x[i]);
          ret += 1.0;
          ret *= 0.5;
        }
  void evaluate (const DomainType & x , RangeFieldType time , RangeType & ret) const
  {
    evaluate ( x , ret );
  }
};
 
// ********************************************************************
template <class FunctionType, int polOrd>
class L2LagrangeProjection
{
  typedef typename DiscreteFunctionType::FunctionSpaceType FunctionSpaceType;

 public:
  static void project (const FunctionType &f, DiscreteFunctionType &discFunc) {
    typedef typename DiscreteFunctionType::Traits::DiscreteFunctionSpaceType 
      FunctionSpaceType;
    typedef typename FunctionSpaceType::Traits::GridType GridType;
    typedef typename FunctionSpaceType::Traits::IteratorType Iterator;
    typedef typename FunctionSpaceType::Traits::RangeType RangeType;
    typedef typename FunctionSpaceType::Traits::DomainType DomainType;

    const FunctionSpaceType& space =  discFunc.getFunctionSpace();

    discFunc.clear();

    typedef typename DiscreteFunctionType::LocalFunctionType LocalFuncType;

    typename FunctionSpaceType::RangeType ret (0.0);

    Iterator it = space.begin();
    Iterator endit = space.end();
    
    for( ; it != endit ; ++it) {
      LocalFuncType lf = discFunc.localFunction(*it);
      const typename GridType::template Codim<0>::Entity::Geometry& 
          itGeom = (*it).geometry();
      
      // assert 2D and poldeg 1 => Lagrange points are corners of element
      assert (polOrd==1 && dimp==2);
      
      for (int lP=0;lP < 4;lP ++)
      {
        RangeType ret;
        // cycle through the 4 points (0,0),(1,0),(1,1),(0,1)
        // and evaluate function
        DomainType lpoint(0.0);
        if (lP==2 || lP==1)  
            lpoint[0] = 1.0;
        if (lP==2 || lP==3)  
            lpoint[1] = 1.0;
	f.evaluate(itGeom.global(lpoint), ret);
        lf[lP] = ret ;           
      }   
    }
  }
};

// calculates || u-u_h ||_L2
template <class DiscreteFunctionType>
class L2Error
{
  typedef typename DiscreteFunctionType::FunctionSpaceType FunctionSpaceType;

public:
  template <int polOrd, class FunctionType>
  double norm (FunctionType &f, DiscreteFunctionType &discFunc,
      double time)
  {
    const typename DiscreteFunctionType::FunctionSpaceType
        & space = discFunc.getFunctionSpace();

    typedef typename FunctionSpaceType::GridType GridType;
    typedef typename FunctionSpaceType::IteratorType IteratorType;
    typedef typename DiscreteFunctionType::LocalFunctionType LocalFuncType;
    typedef typename FunctionSpaceType::RangeType RangeType;

    RangeType ret (0.0);
    RangeType phi (0.0);

    double sum = 0.0;

    IteratorType it    = space.begin();
    IteratorType endit = space.end();

    for(; it != endit ; ++it)
    {
      const int codim = 0;
      ElementQuadrature<GridType,codim> quad(*it, polOrd);
      LocalFuncType lf = discFunc.localFunction(*it);
      for(int qP = 0; qP < quad.nop(); qP++)
      {
        double det = (*it).geometry().integrationElement(quad.point(qP));
        f.evaluate((*it).geometry().global(quad.point(qP)),time, ret);
        lf.evaluate((*it),quad,qP,phi);
        sum += det * quad.weight(qP) * SQR(ret[0] - phi[0]);
      }
    }
    return sqrt(sum);
  }
};

// ********************************************************************
double algorithm (GridType& grid, bool visualize )
{
   GridPartType part ( grid);
   DiscreteFunctionSpaceType linFuncSpace ( part );
   DiscreteFunctionType solution ( "sol", linFuncSpace );
   solution.clear();
      
   ExactSolution f ( linFuncSpace ); 
   L2Error < DiscreteFunctionType > l2err;
       
   //! perform l2-projection
   L2LagrangeProjection<ExactSolution, polOrd>::
     project(f, solution);

   // calculation L2 error 
   // pol ord for calculation the error chould by higher than 
   // pol for evaluation the basefunctions 
   double error = l2err.norm<polOrd + 4> (f ,solution, 0.0);
   std::cout << "\nL2 Error : " << error << "\n\n";
  
// #if HAVE_GRAPE
   // if Grape was found, then display last solution 
   if(visualize)
   {
     GrapeDataDisplay < GridType > grape(grid); 
     grape.dataDisplay( solution );
   }
// #endif
   
   return error;
}


//**************************************************
//
//  main programm, run algorithm twice to calc EOC 
//
//**************************************************
int main (int argc, char **argv)
{
  if(argc != 2)
  {
    fprintf(stderr,"usage: %s <maxlevel> \n",argv[0]);
    exit(1);
  }
  int ml = atoi( argv[1] );
  double error[2];
  
  int n[dimp];
  double h[dimp];
  for(int i=0; i<dimp; i++)  { n[i] = 2; h[i] = 1.0; }
  
  GridType grid ((int *) &n, (double *) &h );
  
  grid.globalRefine(ml);
  
  bool visualize = false;
  error[0] = algorithm ( grid , visualize);
  grid.globalRefine(1);
  visualize = true;
  error[1] = algorithm ( grid , visualize);
  double eoc = log( error[0]/error[1]) / M_LN2; 
  std::cout << "EOC = " << eoc << " \n";

  return 0;
}