diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
new file mode 100644
index 0000000000000000000000000000000000000000..3749fab61bd8f7c1e7fb23c676a4e869278ea112
--- /dev/null
+++ b/Solver/dgAlgorithm.h
@@ -0,0 +1,30 @@
+#ifndef _DG_ALGORITHM_H_
+#define _DG_ALGORITHM_H_
+
+#include "dofManager.h"
+
+class groupOfElements;
+class groupOfFaces;
+class dgConservationLaw;
+
+class dgAlgorithm () {
+ public :
+  void residualVolume ( dofManager &dof,
+			dgConservationLaw &law,
+			groupOfElements & group ,
+			fullMatrix<double> &solution,
+			fullMatrix<double> &residual);
+  void residualInterface ( dofManager &dof,
+			   dgConservationLaw &law,
+			   groupOfFaces & group ,
+			   fullMatrix<double> &solution,
+			   fullMatrix<double> &residual);
+  void residualBoundaryConditions ( dofManager &dof,
+				    dgConservationLaw &law,
+				    groupOfElements & group ,
+				    fullMatrix<double> &solution,
+				    fullMatrix<double> &residual);
+};
+
+
+#endif
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
new file mode 100644
index 0000000000000000000000000000000000000000..805806f3abca51fcb9d0aef7e6091630c4548e76
--- /dev/null
+++ b/Solver/dgConservationLaw.h
@@ -0,0 +1,46 @@
+#ifndef _DG_CONSERVATION_LAW_H_
+#define _DG_CONSERVATION_LAW_H_
+
+/*
+  \partial_t L(u) =   \nabla \cdot (\vec{f}(u,forcings))  -> convective flux f
+                    + \nabla \cdot (\vec{g}(u,\nabla u,forcings)  -> diffusive flux g
+		    + r(u,forcings)                       -> source term r
+*/
+
+class dgElement; 
+class dgFace;
+
+class dgTerm{
+ public:
+  virtual ~dgTerm () {}
+  virtual void operator () (const dgElement &, fullMatrix<double> &fcx) const = 0;
+};
+
+class dgFaceTerm{
+ public:
+  virtual ~dgFaceTerm () {}
+  virtual void operator () (const dgFace &, fullMatrix<double> &fcx) const = 0;
+};
+
+class dgConservationLaw {
+  int _nbf;
+  dgTerm *_diffusive, *_convective, *_source, *_maxConvectiveSpeed;
+  dgFaceTerm *_riemannSolver;
+public:
+  dgConservationLaw () : _diffusive(0), _convective(0), _source (0), 
+			 _riemannSolver(0),_maxConvectiveSpeed (0) {}
+  ~dgConservationLaw () {}
+
+  int nbFields() const {return nbf;}
+
+  inline const dgTerm     * convectiveFlux () {return _convective;}
+  inline const dgTerm     * diffusiveFlux  () {return _diffusive;}
+  inline const dgTerm     * sourceTerm     () {return _source;}
+  inline const dgFaceTerm * riemannSolver  () {return _riemannSolver;}
+  inline const dgTerm     * maxConvectiveSpeed () {return _maxConvectiveSpeed;}
+
+};
+
+
+
+#endif
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..899a374d0cbf309334f9e4dfe45ac3be81a50370
--- /dev/null
+++ b/Solver/dgGroupOfElements.cpp
@@ -0,0 +1,183 @@
+#include "dgGroupOfElements.h"
+#include "MElement.h"
+#include "functionSpace.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;
+}
+
+
+// _collocation(i,j) = fs(i)(point(j))
+static fullMatrix<double> * dgGetCollocation (const functionSpace &fs, 
+					      fullMatrix<double> *intg){
+  fullMatrix<double> *m = new fullMatrix<double> (fs.coefficients.size1(), 
+						  intg->size1());
+  double *sf = new double [fs.coefficients.size1()];
+  for (int i=0;i<intg->size1();i++){
+    fs.f((*intg)(i,0),(*intg)(i,1),(*intg)(i,2),sf);
+    for (int j=0;j<fs.coefficients.size1();j++){
+      (*m)(j,i) = sf[j];
+    }
+  }  
+  delete [] sf;
+  return m;
+}
+
+
+
+dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOrder)
+  : _elements(e), 
+    _fs(_elements[0]->getFunctionSpace(polyOrder)),
+    _integration(dgGetIntegrationRule (_elements[0], polyOrder)),
+    _collocation(dgGetCollocation(_fs,_integration))    
+{
+  // this is the biggest piece of data ... the mappings
+  _mapping = new fullMatrix<double> (_elements.size(), 10 * _integration->size1());
+  for (int i=0;i<_elements.size();i++){
+    MElement *e = _elements[i];
+    for (int j=0;j< _integration->size1() ; j++ ){
+      double ijac[3][3];
+      (*_mapping)(i,10*j + 9) =  
+	e->getJacobian ((*_integration)(j,0),
+			(*_integration)(j,1),
+			(*_integration)(j,2),
+			ijac);
+      (*_mapping)(i,10*j + 0) = ijac[0][0]; 
+      (*_mapping)(i,10*j + 1) = ijac[0][1]; 
+      (*_mapping)(i,10*j + 2) = ijac[0][2]; 
+      (*_mapping)(i,10*j + 3) = ijac[1][0]; 
+      (*_mapping)(i,10*j + 4) = ijac[1][1]; 
+      (*_mapping)(i,10*j + 5) = ijac[1][2]; 
+      (*_mapping)(i,10*j + 6) = ijac[2][0]; 
+      (*_mapping)(i,10*j + 7) = ijac[2][1]; 
+      (*_mapping)(i,10*j + 8) = ijac[2][2]; 
+    }
+  }
+  // redistribution matrix
+  // quadrature weight x parametric gradients in quadrature points
+
+  _redistribution[0] = new fullMatrix<double> (_fs.coefficients.size1(),_integration->size1());
+  _redistribution[1] = new fullMatrix<double> (_fs.coefficients.size1(),_integration->size1());
+  _redistribution[2] = new fullMatrix<double> (_fs.coefficients.size1(),_integration->size1());
+
+  double g[256][3];
+  
+  for (int j=0;j<_integration->size1();j++) {
+    _fs.df((*_integration)(j,0),
+	   (*_integration)(j,1),
+	   (*_integration)(j,2), g);
+    const double weight = (*_integration)(j,3);
+    for (int k=0;k<_fs.coefficients.size1();k++){ 
+      (*_redistribution[0])(k,j) = g[j][0] * weight;
+      (*_redistribution[1])(k,j) = g[j][1] * weight;
+      (*_redistribution[2])(k,j) = g[j][2] * weight;
+    }
+  }
+}
+
+dgGroupOfElements::~dgGroupOfElements(){
+  delete _integration;
+  delete _redistribution[0];
+  delete _redistribution[1];
+  delete _redistribution[2];
+  delete _mapping;
+  delete _collocation;
+}
+
+
+// dgGroupOfFaces
+void dgGroupOfFaces::dgCreateFaceElements (const std::vector<MFace> &topo_faces){
+  
+  _faces.resize(topo_faces.size());
+  // compute all closures
+  for (int i=0;i<topo_faces.size();i++){
+    // compute closures for the interpolation
+    int ithFace, sign, rot;
+    _left[i]->getFaceInfo (topo_faces[i], ithFace, sign, rot);
+    _closuresLeft.push_back(&(_fsLeft.getFaceClosure(ithFace, sign, rot)));
+    _right[i]->getFaceInfo (topo_faces[i], ithFace, sign, rot);
+    _closuresRight.push_back(&(_fsRight.getFaceClosure(ithFace, sign, rot)));        
+    // compute the face element that correspond to the geometrical closure
+    // get the vertices of the face
+    std::vector<MVertex*> _vertices;
+    std::vector<int> & geomClosure = _right[i]->getFunctionSpace().getFaceClosure(ithFace, sign, rot);
+    for (int j=0; j<geomClosure.size() ; j++){
+      int iNod = geomClosure[j];
+      MVertex *v = _left[i]->getVertex(iNod);
+      _vertices.push_back(v);
+    }
+    // triangular face
+    if (topo_faces[i]->getNumVertices() == 3){
+      switch(_vertices.size()){
+      case 3  : _faces.push_back(new MTriangle (_vertices) ); break;
+      case 6  : _faces.push_back(new MTriangle2 (_vertices) ); break;
+      case 10 : _faces.push_back(new MTriangleN (_vertices, 3) ); break;
+      case 15 : _faces.push_back(new MTriangleN (_vertices, 4) ); break;
+      case 21 : _faces.push_back(new MTriangleN (_vertices, 5) ); break;
+      default : throw;
+      }
+    }
+    // quad face 2 do
+    else throw;
+  }  
+}
+
+// dgGroupOfFaces
+void dgGroupOfFaces::dgCreateFaceElements (const std::vector<MEdge> &topo_edges){
+  
+  _faces.resize(topo_edges.size());
+  // compute all closures
+  for (int i=0;i<topo_edges.size();i++){
+    int ithEdge, sign;
+    _left[i]->getEdgeInfo (topo_edges[i], ithEdge, sign);
+    _closuresLeft.push_back(&(_fsLeft->getEdgeClosure(ithEdge, sign)));
+    _right[i]->getEdgeInfo (topo_edges[i], ithEdge, sign);
+    _closuresRight.push_back(&(_fsRight->getEdgeClosure(ithEdge, sign)));    
+    // get the vertices of the edge
+    std::vector<int> & geomClosure = _right[i]->getFunctionSpace().getEdgeClosure(ithFace, sign);
+    std::vector<MVertex*> _vertices;
+    for (int j=0; j<geomClosure.size() ; j++){
+      int iNod = geomClosure[j];
+      MVertex *v = _left[i]->getVertex(iNod);
+      _vertices.push_back(v);
+    }
+    switch(_vertices.size()){
+    case 2  : _faces.push_back(new MLine (_vertices) ); break;
+    case 3  : _faces.push_back(new MLine2 (_vertices) ); break;
+    default : _faces.push_back(new MLineN (_vertices) ); break;
+    }
+    else throw;
+  }
+}
+
+
+dgGroupOfFaces::dgGroupOfFaces (const std::vector<MFace> &topo_faces, 		  
+				const std::vector<MElement*> &l, 
+				const std::vector<MElement*> &r,
+				int pOrder) : _left(l), _right(r) 
+{
+  _fsLeft  = _left[0]->getFunctionSpace (pOrder);
+  _fsRight = _right[0]->getFunctionSpace (pOrder);
+  dgCreateFaceElements (topo_faces);
+}
+
+dgGroupOfFaces::dgGroupOfFaces (const std::vector<MEdge> &topo_edges, 		  
+				const std::vector<MElement*> &l, 
+				const std::vector<MElement*> &r,
+				int pOrder) : _left(l), _right(r) 
+{
+  _fsLeft  = _left[0]->getFunctionSpace (pOrder);
+  _fsRight = _right[0]->getFunctionSpace (pOrder);
+  dgCreateFaceElements (topo_edges);
+}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
new file mode 100644
index 0000000000000000000000000000000000000000..f8db09450bff8c7f36ff50021054726e2e52cf2c
--- /dev/null
+++ b/Solver/dgGroupOfElements.h
@@ -0,0 +1,91 @@
+#ifndef _DG_GROUP_OF_ELEMENTS_
+#define _DG_GROUP_OF_ELEMENTS_
+
+#include <vector>
+#include "fullMatrix.h"
+/*
+  A group of element contains N elements
+  that are of the same type and of the 
+  same order. All element have the same
+  number of integration points Np, the same
+  number of vertices Nv.
+
+  we DO NOT store N, Ni or Nv (matrices and vectors 
+  have the right sizes)
+*/
+class MElement;
+class functionSpace;
+
+class dgGroupOfElements {
+  // the ONLY function space that is used to 
+  // inerpolated the fields (may be different to the 
+  // one of the elements)
+  const functionSpace &_fs;
+  // N elements in the group
+  std::vector<MElement*> _elements;
+  // Ni integration points, matrix of size Ni x 4 (u,v,w,weight)
+  fullMatrix<double> *_integration;
+  // collocation matrix that maps vertices to integration points.
+  fullMatrix<double> *_collocation;
+  // inverse of the mapping and det mapping (00 01 02 10 11 12 20 21 22 det)
+  fullMatrix<double> *_mapping;
+  // redistribution of the fluxes to vertices multiplied by
+  // the gradient of shape functions (in parametric space)
+  // for both diffusive and convective fluxes
+  fullMatrix<double> *_redistribution[3];
+  // redistribution for the source term
+  fullMatrix<double> *_redistributionSource;
+  // forbid the copy 
+  //  dgGroupOfElements (const dgGroupOfElements &e, int order) {}
+  //  dgGroupOfElements & operator = (const dgGroupOfElements &e) {}
+public:
+  dgGroupOfElements (const std::vector<MElement*> &e, int);
+  virtual ~dgGroupOfElements ();
+};
+
+class dgFace {
+  int nbGaussPoints;
+  MElement *_left, *_right;
+  MElement *_face;
+  double *_normals;
+  double *_detJ;
+  int *_closureLeft,*_closureRight;
+public:
+  dgFace ();
+};
+
+class dgGroupOfFaces {
+  // Two functionSpaces for left and right elements
+  // the group has always the same types for left and right
+  const functionSpace &_fsLeft,&_fsRight, &_fsFace;
+  // N elements in the group
+  std::vector<MElement*> _left, _right, _faces;
+  // Ni integration points, matrix of size Ni x 3 (u,v,weight)
+  fullMatrix<double> *_integration;
+  // there is a finite number of combinations of orientations, senses
+  // and rotations of the faces (typically 24 for tets). Any pair
+  // is characterized by a single integer which is the combination
+  // this closure is for the interpolation that MAY BE DIFFERENT THAN THE
+  // GEOMETRICAL CLOSURE !!!
+  std::vector<std::vector<int> * > _closuresLeft; 
+  std::vector<std::vector<int> * > _closuresRight; 
+  // normals at integration points 
+  fullMatrix<double> *_normals;
+  // detJac at integration points
+  fullMatrix<double> *_detJac;
+  // collocation matrices \psi_i (GP_j) 
+  fullMatrix<double> *_collocationLeft, *_collocationRight;
+  // redistribution matrices \psi_i (GP_j) * weight_j
+  fullMatrix<double> *_redistributionLeft, *_redistributionRight;
+  
+public:
+  dgGroupOfFaces (const std::vector<MFace> &faces, 		  
+		  const std::vector<MElement*> &l, 
+		  const std::vector<MElement*> &r,
+		  int pOrder);
+  dgFace operator () (int iFace) const;
+  virtual ~dgGroupOfFaces ();
+};
+
+
+#endif