diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 138f94b0f84621e1ff0839e2365c5990be533051..848b14e91ff9c0c3bd454a64ad715a300c2f6f3b 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -60,7 +60,7 @@ void fullMatrix<std::complex<double> >::mult(const fullMatrix<std::complex<doubl
 }
 
 template<> 
-void fullMatrix<double>::gemm(fullMatrix<double> &a, fullMatrix<double> &b, 
+void fullMatrix<double>::gemm(const fullMatrix<double> &a, const fullMatrix<double> &b, 
                               double alpha, double beta)
 {
   int M = size1(), N = size2(), K = a.size2();
@@ -70,8 +70,8 @@ void fullMatrix<double>::gemm(fullMatrix<double> &a, fullMatrix<double> &b,
 }
 
 template<> 
-void fullMatrix<std::complex<double> >::gemm(fullMatrix<std::complex<double> > &a, 
-                                             fullMatrix<std::complex<double> > &b, 
+void fullMatrix<std::complex<double> >::gemm(const fullMatrix<std::complex<double> > &a, 
+                                             const fullMatrix<std::complex<double> > &b, 
                                              std::complex<double> alpha, 
                                              std::complex<double> beta)
 {
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index e3387a26f5112ec055c14322dd784a6a4701baa7..88328d151e2d9ba0c317ec145c60cfc5a6ca8bcd 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -143,7 +143,7 @@ class fullMatrix
   }
 #endif
   ;
-  void gemm(fullMatrix<scalar> &a, fullMatrix<scalar> &b, 
+  void gemm(const fullMatrix<scalar> &a, const fullMatrix<scalar> &b, 
             scalar alpha=1., scalar beta=1.)
 #if !defined(HAVE_BLAS)
   {
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index ecd48affd130e24251dbf0e17562471150107678..aad2bad669a4a6f6b7069257509a48e0c6dedf3b 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -11,6 +11,7 @@ set(SRC
   SElement.cpp
   eigenSolver.cpp
   dgGroupOfElements.cpp
+  dgAlgorithm.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h) 
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index b86b8279f930b163c95cdf5242475d674fd0c2f6..f0626de940661e85437fd0acb483b860decbe3d1 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -1,45 +1,48 @@
 #include "dgAlgorithm.h"
+#include "dgGroupOfElements.h"
+#include "dgConservationLaw.h"
 
 /*
   compute 
     \int \vec{f} \cdot \grad \phi dv   
 */
 
-void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe useless here)
+void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
 				   const dgConservationLaw &claw,   // the conservation law
-				   const groupOfElements & group ) // the residual
+				   const dgGroupOfElements & group ) // the residual
 { 
   // ----- 1 ----  get the solution at quadrature points
   // ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints) 
-  fullMatrix<double> solutionQP (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints());
+  int nbFields = group.getNbFields();
+  fullMatrix<double> solutionQP (group.getNbElements() * nbFields, group.getNbIntegrationPoints());
   // ----- 1.2 --- multiply the solution by the collocation matrix
   group.getSolution().mult ( group.getCollocationMatrix() , solutionQP); 
   // ----- 1.3 --- if the conservation law is diffusive, compute the gradients too
   fullMatrix<double> *gradientSolutionQP= 0;
   if (claw.diffusiveFlux()){
-     gradientSolutionQP = new  fullMatrix<double> (group.getNbElements() * group.getNbFields() * 3, group.getNbIntegrationPoints());
+     gradientSolutionQP = new  fullMatrix<double> (group.getNbElements() * nbFields * 3, group.getNbIntegrationPoints());
      group.getGradientOfSolution().mult ( group.getCollocationMatrix() , *gradientSolutionQP); 
   }
   
 
   // ----- 2 ----  compute all fluxes (diffusive and convective) at integration points
   // ----- 2.1 --- allocate elementary fluxes (compued in XYZ coordinates)
-  fullMatrix<double> fConv[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
-				 fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
-				 fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() )};
-  fullMatrix<double> fDiff[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
-				 fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
-				 fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() )};
+  fullMatrix<double> fConv[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
+				 fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
+				 fullMatrix<double>( group.getNbIntegrationPoints(), nbFields )};
+  fullMatrix<double> fDiff[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
+				 fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
+				 fullMatrix<double>( group.getNbIntegrationPoints(), nbFields )};
   // ----- 2.2 --- allocate parametric fluxes (compued in UVW coordinates) for all elements at all integration points
-  fullMatrix<double> Fuvw[3] = {fullMatrix<double> (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints()),
-				fullMatrix<double> (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints()),
-				fullMatrix<double> (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints())};
+  fullMatrix<double> Fuvw[3] = {fullMatrix<double> (group.getNbElements() * nbFields, group.getNbIntegrationPoints()),
+				fullMatrix<double> (group.getNbElements() * nbFields, group.getNbIntegrationPoints()),
+				fullMatrix<double> (group.getNbElements() * nbFields, group.getNbIntegrationPoints())};
   // ----- 2.3 --- iterate on elements
   for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {        
     // ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
     fullMatrix<double> solutionQPe (solutionQP, iElement*claw.nbFields(),claw.nbFields() );
     fullMatrix<double> *gradSolutionQPe;
-    if (claw.diffusiveFlux()) gradSolutionQPe = new fullMatrix<double>(*gradSolutionQP, 3*iElement*claw.nbFields(),3*claw.nbFields() );      
+    if (claw.diffusiveFlux()) gradSolutionQPe = new fullMatrix<double>(*gradientSolutionQP, 3*iElement*claw.nbFields(),3*claw.nbFields() );      
     else gradSolutionQPe = new fullMatrix<double>;
     dgElement DGE( group.getElement(iElement), solutionQPe, *gradSolutionQPe, group.getIntegrationPointsMatrix());
     // ----- 2.3.2 --- compute fluxes in XYZ coordinates
@@ -59,7 +62,7 @@ void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe us
 	  const double detJ = group.getDetJ (iElement, iPt);
 	  const double factor = invJ * detJ;
 	  // compute fluxes in the reference coordinate system at this integration point
-	  for (int k=0;k<group.nbFields();k++) {
+	  for (int k=0;k<nbFields;k++) {
 	    fuvwe(iPt,k) += ( fConv[iXYZ](iPt,k) + fDiff[iXYZ](iPt,k)) * factor;
 	  }
 	}
@@ -68,5 +71,5 @@ void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe us
   }
   // ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
   for (int iUVW=0;iUVW<group.getDimUVW();iUVW++)
-    group.getResidual().dgemm(group.getRedistributionMatrix(iUVW), Fuvw[iUVW]);
+    group.getResidual().gemm(group.getFluxRedistributionMatrix(iUVW), Fuvw[iUVW]);
 }
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 3749fab61bd8f7c1e7fb23c676a4e869278ea112..dc487b640ec98f775f4e49e4efa164ae2c64a2be 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -3,27 +3,21 @@
 
 #include "dofManager.h"
 
-class groupOfElements;
-class groupOfFaces;
+class dgGroupOfElements;
+class dgGroupOfFaces;
 class dgConservationLaw;
 
-class dgAlgorithm () {
+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);
+  void residualVolume ( /*dofManager &dof,*/
+			const dgConservationLaw &law,
+			const dgGroupOfElements & group);
+  void residualInterface ( /*dofManager &dof,*/
+			   const dgConservationLaw &law,
+			   const dgGroupOfFaces & group);
+  void residualBoundaryConditions ( /*dofManager &dof,*/
+				    const dgConservationLaw &law,
+				    const dgGroupOfElements & group);
 };
 
 
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 8c5bc9ff33b5d55664715e70846d4d92dd159982..1625b34e25d3a4b6f9502c6f27c2569d37a47c38 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -31,13 +31,13 @@ public:
 			 _riemannSolver(0),_maxConvectiveSpeed (0) {}
   ~dgConservationLaw () {}
 
-  int nbFields() const {return nbf;}
+  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;}
+  inline const dgTerm     * convectiveFlux () const {return _convective;}
+  inline const dgTerm     * diffusiveFlux  () const {return _diffusive;}
+  inline const dgTerm     * sourceTerm     () const {return _source;}
+  inline const dgFaceTerm * riemannSolver  () const {return _riemannSolver;}
+  inline const dgTerm     * maxConvectiveSpeed () const {return _maxConvectiveSpeed;}
 
 };
 
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 21eaa519df9ee1dcf08a5fc7b1176050e1c8d0f2..c8042a0862059f678d22450a6aaf04ae5024ab9c 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -19,14 +19,14 @@ class MEdge;
 class functionSpace;
 
 class dgElement {
-  MElement *_element;
+  const MElement *_element;
   // solution at points
   const fullMatrix<double> &_solution, &_integration, &_gradients;
 public:
-  dgElement (MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &integ)
+  dgElement (const MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &integ)
     : _element(e), _solution(sol), _integration(integ), _gradients(sol)
   {}
-  dgElement (MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &grads, const fullMatrix<double> &integ)
+  dgElement (const MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &grads, const fullMatrix<double> &integ)
     : _element(e), _solution(sol), _integration(integ), _gradients(grads)
   {}
 };
@@ -76,11 +76,11 @@ public:
   inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
   inline const fullMatrix<double> & getFluxRedistributionMatrix (int i) const {return *_redistributionFluxes[i];}
   inline const fullMatrix<double> & getSourceRedistributionMatrix () const {return *_redistributionSource;}
-  inline const fullMatrix<double> & getSolution () const {return *_solution;}
-  inline const fullMatrix<double> & getGradientOfSolution () const {return *_gradSolution;}
+  inline fullMatrix<double> & getSolution () const {return *_solution;}
+  inline fullMatrix<double> & getGradientOfSolution () const {return *_gradSolution;}
   // get a proxy on the solution for element iElement
   inline fullMatrix<double> getSolution (int iElement) const {return fullMatrix<double>(*_solution, iElement*_nbFields, _nbFields);}
-  inline const fullMatrix<double> & getResidual () const {return *_solution;}
+  inline fullMatrix<double> & getResidual () const {return *_solution;}
   // get a proxy on the residual for element iElement
   inline fullMatrix<double> getResidual (int iElement) const {return fullMatrix<double>(*_residual, iElement*_nbFields, _nbFields);}
   inline double getDetJ (int iElement, int iGaussPoint) const {return (*_mapping)(iElement, 10*iGaussPoint + 9);}