///////////////////////////////////////////////////////////////////////////////
// 
// Flussmatrizen fuer LDG Verfahrrn
// code by M. Kraenkel 11.7.06
//
//////////////////////////////////////////////////////////////////////////////
//Flussanteil von M_sigma_u
template<class VecDiscreteFunctionSpc,class ScalDiscreteFunctionSpc,class EntityType>
void flux_gradMatrix(const EntityType& e,
		     VecDiscreteFunctionSpc& vecspc,
		     ScalDiscreteFunctionSpc& scalspc,
		     GlobalMatrixType& ret)
{ 
  const typename VecDiscreteFunctionSpc::BaseFunctionSetType &vecbaseset = vecspc.getBaseFunctionSet(e);
  const typename ScalDiscreteFunctionSpc::BaseFunctionSetType &scalbaseset = scalspc.getBaseFunctionSet(e);
  typedef typename VecDiscreteFunctionSpc::GridType GridType;
  enum{polOrd= VecDiscreteFunctionSpc::polOrd} ;
  typedef typename VecDiscreteFunctionSpc::RangeType VecRangeType;
  typedef typename VecDiscreteFunctionSpc::DomainType DomainType;
  typedef typename ScalDiscreteFunctionSpc::RangeType RangeType;
  typedef typename EntityType::IntersectionIterator IntersectionIterator;
  typedef typename EntityType::EntityPointer EntityPointerType;
  int vecnumDofs= vecbaseset.numBaseFunctions();
  int scalnumDofs= scalbaseset.numBaseFunctions();

  IntersectionIterator endnit = e.iend();
  VecRangeType psi(0.0);
  RangeType enVal,neighVal, enflux,neighflux;
  const RangeType zero(0.0);


  for(IntersectionIterator nit =e.ibegin();nit !=endnit;++nit)
    { 
     
      if (nit.leafNeighbor())
	{
	  EntityPointerType neighEp=nit.outside();
	  EntityType& nb=  *neighEp;
	  
	  const typename ScalDiscreteFunctionSpc::BaseFunctionSetType &scalbasesetNeigh = scalspc.getBaseFunctionSet(nb);
	  ElementQuadrature<GridType,1> innerQuad(nit,2*polOrd,ElementQuadrature<GridType,1>::INSIDE);
	  ElementQuadrature<GridType,1> outerQuad(nit,2*polOrd,ElementQuadrature<GridType,1>::OUTSIDE);
	  
	  const int quadNop = innerQuad.nop();
	   
	  for (int l = 0; l < quadNop ; ++l)
	     {
	       DomainType normal=nit.integrationOuterNormal(innerQuad.localPoint(l));
	       DomainType unitnormal= nit.unitOuterNormal(innerQuad.localPoint(l));
		 for(int j=0;j<vecnumDofs;++j)
		   {
		    int enrow = vecspc.mapToGlobal(e,j);
		    vecbaseset.eval(j,innerQuad,l,psi);
		    double psi_dot_normal = unitnormal*psi;
	    
		    for(int k;k<scalnumDofs;++k)
		      {
			int encol = scalspc.mapToGlobal(e,j);
			int neighcol = scalspc.mapToGlobal(nb,j);
			
			scalbaseset.eval(k,innerQuad,l,enVal);
			scalbasesetNeigh.eval(k,outerQuad,l,neighVal);
		  
			uflux(normal,enVal,zero,enflux);
			uflux(normal,zero,neighVal,neighflux);
			
			enflux *= psi_dot_normal;
			neighflux *= psi_dot_normal;
			
			ret.add(enrow,encol,enflux[0]);
			ret.add(enrow,neighcol,neighflux[0]);
		      }
		  }
	    }
	}
    }
}

//Flussanteil von M_u_sigma
template<class VecDiscreteFunctionSpc,class ScalDiscreteFunctionSpc,class EntityType>
void flux_divMatrix(const EntityType& e,VecDiscreteFunctionSpc& vecspc,ScalDiscreteFunctionSpc& scalspc,GlobalMatrixType& ret)
{ 
  const typename VecDiscreteFunctionSpc::BaseFunctionSetType &vecbaseset = vecspc.getBaseFunctionSet(e);
  const typename ScalDiscreteFunctionSpc::BaseFunctionSetType &scalbaseset = scalspc.getBaseFunctionSet(e);
  typedef typename VecDiscreteFunctionSpc::GridType GridType;
  enum{polOrd= VecDiscreteFunctionSpc::polOrd} ;
  typedef typename VecDiscreteFunctionSpc::RangeType VecRangeType;
  typedef typename VecDiscreteFunctionSpc::DomainType DomainType;
  typedef typename ScalDiscreteFunctionSpc::RangeType RangeType;
  typedef typename EntityType::IntersectionIterator IntersectionIterator;
  typedef typename EntityType::EntityPointer EntityPointerType;


 
  int vecnumDofs= vecbaseset.numBaseFunctions();
  int scalnumDofs= scalbaseset.numBaseFunctions();
  IntersectionIterator endnit = e.iend();
  RangeType phi(0.0);
  VecRangeType enVal,neighVal;
  RangeType enflux,neighflux;
  const RangeType uzero(0.0);
  const VecRangeType sigmazero(0.0);
  for(IntersectionIterator nit =e.ibegin();nit !=endnit;++nit)
    { 
     
     
      {
	//if (nit.leafNeighbor())
	// {
	    EntityPointerType neighEp=nit.outside();
	    EntityType& nb=  *neighEp;

	    const typename VecDiscreteFunctionSpc::BaseFunctionSetType &vecbasesetNeigh = vecspc.getBaseFunctionSet(nb);
	    // }
	ElementQuadrature<GridType,1> innerQuad(nit,2*polOrd,ElementQuadrature<GridType,1>::INSIDE);
	ElementQuadrature<GridType,1> outerQuad(nit,2*polOrd,ElementQuadrature<GridType,1>::OUTSIDE);
	 
	const int quadNop = innerQuad.nop();
	
	for (int l = 0; l < quadNop ; ++l)
	  {
	    DomainType normal=nit.integrationOuterNormal(innerQuad.localPoint(l));
	    for(int j=0;j<scalnumDofs;++j)
	      {
		int enrow = scalspc.mapToGlobal(e,j);
				    
		scalbaseset.eval(j,innerQuad,l,phi);
		    
		for(int k;k<vecnumDofs;++k)
		  {
		    int encol = vecspc.mapToGlobal(e,j);
		    vecbaseset.eval(k,innerQuad,l,enVal);
		    sigmaflux(normal,enVal,sigmazero,uzero,uzero,enflux);
		    ret.add(enrow,encol,enflux[0]*innerQuad.weight(l));
		    
		    if (nit.leafNeighbor()){
		      int neighcol = vecspc.mapToGlobal(nb,j);
		      vecbasesetNeigh.eval(k,outerQuad,l,neighVal);
		      sigmaflux(normal,sigmazero,neighVal,uzero,uzero,neighflux);
		      ret.add(enrow,neighcol,neighflux[0]*innerQuad.weight(l));
		    }
		  }
	      }
	  }
      }
    }
}

//M_u_u(besitzt nur Flussanteil)
template<class ScalDiscreteFunctionSpc,class EntityType>
void stabMatrix(const EntityType& e,ScalDiscreteFunctionSpc& scalspc,GlobalMatrixType& ret)
{ 
 
  const typename ScalDiscreteFunctionSpc::BaseFunctionSetType &scalbaseset = scalspc.getBaseFunctionSet(e);
 
  typedef typename ScalDiscreteFunctionSpc::RangeType RangeType;
  typedef typename ScalDiscreteFunctionSpc::GridType GridType;
  enum{polOrd= ScalDiscreteFunctionSpc::polOrd} ;
  typedef typename ScalDiscreteFunctionSpc::DomainType DomainType;
  typedef typename ScalDiscreteFunctionSpc::RangeType RangeType;
  typedef typename ScalDiscreteFunctionSpc::JacobianRangeType JacobianRangeType;
  typedef typename EntityType::IntersectionIterator IntersectionIterator;
  typedef typename EntityType::EntityPointer EntityPointerType;
  int scalnumDofs= scalbaseset.numBaseFunctions();
  IntersectionIterator endnit = e.iend();
  RangeType enflux,neighflux;
  RangeType phi(0.0);
  RangeType enVal,neighVal;
  const JacobianRangeType sigmazero(0.0);
  const RangeType uzero(0.0); 
  for(IntersectionIterator nit =e.ibegin();nit !=endnit;++nit)
    { 
     
    
	{
	  EntityPointerType neighEp=nit.outside();
	  EntityType& nb=  *neighEp;
	  const typename ScalDiscreteFunctionSpc::BaseFunctionSetType &scalbasesetNeigh = scalspc.getBaseFunctionSet(nb);
	  ElementQuadrature<GridType,1> innerQuad(nit,2*polOrd,ElementQuadrature<GridType,1>::INSIDE);
	  ElementQuadrature<GridType,1> outerQuad(nit,2*polOrd,ElementQuadrature<GridType,1>::OUTSIDE);
	  const int quadNop = innerQuad.nop();
	
	  for (int l = 0; l < quadNop ; ++l)
	    {
	      DomainType normal=nit.integrationOuterNormal(innerQuad.localPoint(l));
		for(int j=0;j<scalnumDofs;++j)
		  {
		    int row = scalspc.mapToGlobal(e,j);
		      
		   
		    scalbaseset.eval(j,innerQuad,l,phi);
		    
		    for(int k;k<scalnumDofs;++k)
		      {
			int encol = scalspc.mapToGlobal(e,j); 
			scalbaseset.eval(k,innerQuad,l,enVal);
			sigmaflux(normal,sigmazero[0],sigmazero[0],enVal,uzero,enflux);
			ret.add(row,encol,enflux[0]*innerQuad.weight(l));
			if (nit.leafNeighbor())
			  {
			    int neighcol = scalspc.mapToGlobal(nb,j);
			    scalbasesetNeigh.eval(k,outerQuad,l,neighVal);
			    sigmaflux(normal,sigmazero[0],sigmazero[0],uzero,neighVal,neighflux);
			    ret.add(row,neighcol,neighflux[0]*innerQuad.weight(l));
			  }
		      }
		  }
	    }
	}
    }
}

#include "fluxes.cc"