/**************************************************************************/
//   Elementmatrizen fuer LDG Verfahren
/**************************************************************************/

// Dune includes
#include"config.h"           // file constructed by ./configure script

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

#include"dune/grid/sgrid.hh" // load sgrid definition
#include"dune/grid/common/gridinfo.hh" // definition of gridinfo

//Dune-fem includes
#include<dune/fem/space/dgspace.hh>
#include <dune/fem/quadrature/elementquadrature.hh>
#include <dune/grid/common/defaultindexsets.hh>
#include <dune/grid/common/gridpart.hh>

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

const int polord = 2;
static const int dim = 2;
static const int range = 1;
typedef FieldMatrix<double,30,30> MatrixType;

// declare function for computation of mass-matrix
template<class VecDiscreteFunctionSpc,class EntityType,class QuadType>
MatrixType massMatrix(const EntityType& e, VecDiscreteFunctionSpc& spc, 
                      QuadType quad);

// declare function for computation of gradient-matrix 
// (after partial integration this corresponds to the term with 
// divergence of vecorial basis function Psi)
template<class VecDiscreteFunctionSpc,class ScalDiscreteFunctionSpc,
         class EntityType,class QuadType>
MatrixType gradMatrix(const EntityType& e, VecDiscreteFunctionSpc& vecspc,
                      ScalDiscreteFunctionSpc& scalspc, QuadType quad);

// declare function for computation of divergence-matrix 
// (after partial integration this corresponds to the term with 
// gradient of scalar basis function Phi)
template<class VecDiscreteFunctionSpc,class ScalDiscreteFunctionSpc,
         class EntityType,class QuadType>
MatrixType divMatrix(const EntityType& e, VecDiscreteFunctionSpc& vecspc,
                     ScalDiscreteFunctionSpc& scalspc, QuadType quad);


/**************************************************************************/
// implementations of matrix computation routines
/**************************************************************************/

template<class VecDiscreteFunctionSpc,class EntityType,class QuadType>
MatrixType massMatrix(const EntityType& e,VecDiscreteFunctionSpc& spc,
                      QuadType quad)
{
  const typename VecDiscreteFunctionSpc::BaseFunctionSetType & 
      baseset = spc.getBaseFunctionSet(e);
  typedef typename VecDiscreteFunctionSpc::RangeType RangeType;
  int numDofs= baseset.numBaseFunctions();
  
  //const int size =numDofs*numDofs;
  MatrixType a(0.0);
  //double *a=new double[size];
  
  double val;
  RangeType psi;
  for(int qp=0;qp<quad.nop();qp++){
    
    // volumenelemt und quadweight
    double intel=
        quad.weight(qp)*e.geometry().integrationElement(quad.point(qp));
    
    for(int i=0;i<numDofs;i++)
    {	
      baseset.eval(i,quad,qp,psi);
      for(int j=0;j<numDofs;j++)
      {
        val=baseset.evaluateSingle(j,quad,qp,psi);
        val*=intel;
        a[i][j]+=val;
        
      }
    }
  }
  for(int i=0;i<numDofs;i++)
  {
    for(int j=0;j<numDofs;j++)
    {
      std::cout<<a[i][j]<<" ";
    }
    std::cout<<"\n";
  }
  
  return a;
}



template<class VecDiscreteFunctionSpc,class ScalDiscreteFunctionSpc, 
         class EntityType,class QuadType>
MatrixType gradMatrix(const EntityType& e,VecDiscreteFunctionSpc& vecspc, 
                      ScalDiscreteFunctionSpc& scalspc,QuadType quad)
{  
  const typename VecDiscreteFunctionSpc::BaseFunctionSetType & 
      vecbaseset = vecspc.getBaseFunctionSet(e);
  const typename ScalDiscreteFunctionSpc::BaseFunctionSetType &
      scalbaseset = scalspc.getBaseFunctionSet(e);
  typedef typename VecDiscreteFunctionSpc::JacobianRangeType JacobianRangeType;
  typedef typename ScalDiscreteFunctionSpc::RangeType RangeType;
  int vecnumDofs= vecbaseset.numBaseFunctions();
  int scalnumDofs= scalbaseset.numBaseFunctions();
  
  MatrixType a=(0.0);
  double val;
  
  JacobianRangeType helpmat(0.0);
  RangeType phi;
  for(int qp=0;qp<quad.nop();qp++){
    
    // volumenelemt und quadweight
    double intel =
        quad.weight(qp)*e.geometry().integrationElement(quad.point(qp));
    for(int j=0;j<scalnumDofs;j++)
    {
      scalbaseset.eval(j,quad,qp,phi);
      
      for(int k = 0; k < dim ; ++k)
	  helpmat[k][k]=phi[0];
      
      for(int i=0;i<vecnumDofs;i++)
      {
        val=vecbaseset.evaluateGradientSingle(i,e,quad,qp,helpmat);
        val*=intel;
        a[i][j]+=val;
      }
      
    }
  }
  for(int i=0;i<vecnumDofs;i++)
  {
    for(int j=0;j<scalnumDofs;j++)
    {
      std::cout<<a[i][j]<<" ";
    }
    std::cout<<"\n";
  }
  
  return a; 
}


template<class VecDiscreteFunctionSpc,class ScalDiscreteFunctionSpc,
         class EntityType,class QuadType>
MatrixType divMatrix(const EntityType& e,VecDiscreteFunctionSpc& vecspc, 
                     ScalDiscreteFunctionSpc& scalspc,QuadType quad){
  
  const typename VecDiscreteFunctionSpc::BaseFunctionSetType &
      vecbaseset = vecspc.getBaseFunctionSet(e);
  const typename ScalDiscreteFunctionSpc::BaseFunctionSetType &
      scalbaseset = scalspc.getBaseFunctionSet(e);
  typedef typename VecDiscreteFunctionSpc::RangeType RangeType;


  int vecnumDofs= vecbaseset.numBaseFunctions();
  int scalnumDofs= scalbaseset.numBaseFunctions();
  
  MatrixType a(0.0);
  double val;
  
  
  RangeType phi;
  for(int qp=0;qp<quad.nop();qp++){
    
    // volumenelemt und quadweight
    double intel=quad.weight(qp)*e.geometry().integrationElement(quad.point(qp));
    for(int j=0;j<vecnumDofs;j++)
    {
      vecbaseset.eval(j,quad,qp,phi);
      
      for(int i=0;i<scalnumDofs;i++)
      {
        val=scalbaseset.evaluateGradientSingle(i,e,quad,qp,phi[0]);	 
        val*=intel;       
	    a[i][j]+=val;
      }
    }
  }
  for(int i=0;i<scalnumDofs;i++)
  {
    for(int j=0;j<vecnumDofs;j++)
    {
      std::cout<<a[i][j]<<" ";
    }
    std::cout<<"\n";
  }
  return a;
}

/**************************************************************************/
// main program
/**************************************************************************/
int main()
{ 
  
  typedef Dune::SGrid<dim,dim> GridType;
  typedef Dune::LeafGridPart < GridType > GridPartType; 
 

  typedef GridType::Codim<0>::Entity EntityType;
  
  typedef GridType::Codim<0>::LeafIterator ElementLeafIterator;

  
  typedef FunctionSpace< double , double, dim , range > ScalFunctionSpaceType;
  typedef FunctionSpace< double , double, dim , dim > VecFunctionSpaceType;


  typedef DiscontinuousGalerkinSpace<ScalFunctionSpaceType, 
    GridPartType,polord> ScalDiscreteFunctionSpaceType;
  typedef DiscontinuousGalerkinSpace<VecFunctionSpaceType, 
    GridPartType,polord> VecDiscreteFunctionSpaceType;
  
  typedef  ScalDiscreteFunctionSpaceType::IteratorType ScalIteratorType;
  typedef  VecDiscreteFunctionSpaceType::IteratorType VecIteratorType;

  Dune::FieldVector<int,dim> N(10);
  Dune::FieldVector<GridType::ctype,dim> L(0.0);
  Dune::FieldVector<GridType::ctype,dim> H(1.0);
  GridType grid(N,L,H);
  GridPartType gridpart(grid);
  
  MatrixType massmat;
  MatrixType gradmat;
  MatrixType divmat;

  ScalDiscreteFunctionSpaceType scalspace(gridpart);
  VecDiscreteFunctionSpaceType vecspace(gridpart);
 
  VecIteratorType vecit =vecspace.begin();
  ElementQuadrature<GridType,0> quad(*vecit, 2*polord);
  std::cout<<"MassMatrix:\n";
  massmat = massMatrix(*vecit ,vecspace,quad);
  std::cout<<"GradientMatrix:\n";
  gradmat = gradMatrix(*vecit ,vecspace,scalspace,quad);
  std::cout<<"DivergenceMatrix:\n";
  divmat = divMatrix(*vecit ,vecspace,scalspace,quad);
  
  return 0;
}