_

#include "dgGroupOfElements.h"
#include "MElement.h"
#include "functionSpace.h"
#include "Numeric.h"
#include "MTriangle.h"
#include "MLine.h"

static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
  int npts;
  IntPt *pts;
  e->getIntegrationPoints(2*p+1, &npts, &pts);
  fullMatrix<double> *m = new fullMatrix<double> (npts, 4);
  for (int i=0;i<npts;i++){
    (*m)(i,0) = pts[i].pt[0];
    (*m)(i,1) = pts[i].pt[1];
    (*m)(i,2) = pts[i].pt[2];
    (*m)(i,3) = pts[i].weight;
  }
  return m;
}

static fullMatrix<double> * dgGetFaceIntegrationRuleOnElement (
      const functionSpace *fsFace,
      const fullMatrix<double> &intgFace,
      const functionSpace *fsElement,
      const std::vector <int> *closure) {
  int npts=intgFace.size1();
  fullMatrix<double> *m = new fullMatrix<double> (npts, 4);
  double f[256];
  for (int i=0;i<npts;i++){
    fsFace->f(intgFace(i,0),intgFace(i,1),intgFace(i,2),f);
    for(size_t j=0; j<closure->size();j++){
      int jNod=(*closure)[j];
      (*m)(i,0) += f[j] * fsElement->points(jNod,0);
      (*m)(i,1) += f[j] * fsElement->points(jNod,1);
      (*m)(i,2) += f[j] * fsElement->points(jNod,2);
      (*m)(i,3) = intgFace(i,3);
    }
  }
  return m;
}


dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOrder)
  : _elements(e), 
    _fs(*_elements[0]->getFunctionSpace(polyOrder)),
    _integration(dgGetIntegrationRule (_elements[0], polyOrder)
    )
{
  _dimUVW = _dimXYZ = e[0]->getDim();
  // this is the biggest piece of data ... the mappings
  int nbNodes = _fs.coefficients.size1();
  _redistributionFluxes[0] = new fullMatrix<double> (nbNodes,_integration->size1());
  _redistributionFluxes[1] = new fullMatrix<double> (nbNodes,_integration->size1());
  _redistributionFluxes[2] = new fullMatrix<double> (nbNodes,_integration->size1());
  _redistributionSource = new fullMatrix<double> (nbNodes,_integration->size1());
  _collocation = new fullMatrix<double> (_integration->size1(),nbNodes);
  _mapping = new fullMatrix<double> (e.size(), 10 * _integration->size1());
  _imass = new fullMatrix<double> (nbNodes,nbNodes*e.size()); 
  double g[256][3],f[256];
  for (int i=0;i<_elements.size();i++){
    MElement *e = _elements[i];
    fullMatrix<double> imass(*_imass,nbNodes*i,nbNodes);
    for (int j=0;j< _integration->size1() ; j++ ){
      _fs.f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
      double jac[3][3],ijac[3][3],detjac;
      (*_mapping)(i,10*j + 9) =  
        e->getJacobian ((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), jac);
      const double weight = (*_integration)(j,3);
      detjac=inv3x3(jac,ijac);