diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 0f2fa8fa8792b6d2dbe8413c55cc32c13b189648..8eba43f94c5e457656dfe190586a238eff016618 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -27,6 +27,7 @@
 #include "dgFunctionIntegrator.h"
 #include "Bindings.h"
 #include "dgResidual.h"
+#include "drawContext.h"
 
 extern "C" {
   #include "lua.h"
@@ -196,6 +197,10 @@ static int luaClear (lua_State *L){
   return 0;
 }
 #endif
+static int luaRefresh (lua_State *L){
+  drawContext::global()->draw();
+  return 0;
+}
 
 int binding::readFile(const char *filename)
 {
@@ -319,6 +324,7 @@ binding::binding(){
   #ifdef HAVE_READLINE
   lua_register(L,"saveHistory",luaSave);
   lua_register(L,"clearHistory",luaClear);
+  lua_register(L,"refreshGraphics",luaRefresh);
   #endif
 
   //  lua_pushcfunction(L, luaopen_io);
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 1836721e214d16e51de15e55b158589ca65bf800..08c0e9c256917aeab1e58c06ee1ad01b610a5952 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -64,15 +64,16 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
   _order=polyOrder;
   _dimUVW = _dimXYZ = e[0]->getDim();
   // this is the biggest piece of data ... the mappings
-  int nbPsi = _fs.coefficients.size1();
-  _redistributionFluxes[0] = new fullMatrix<double> (nbPsi,_integration->size1());
-  _redistributionFluxes[1] = new fullMatrix<double> (nbPsi,_integration->size1());
-  _redistributionFluxes[2] = new fullMatrix<double> (nbPsi,_integration->size1());
-  _redistributionSource = new fullMatrix<double> (nbPsi,_integration->size1());
-  _collocation = new fullMatrix<double> (_integration->size1(),nbPsi);
-  _mapping = new fullMatrix<double> (e.size(), 10 * _integration->size1());
-  _imass = new fullMatrix<double> (nbPsi,nbPsi*e.size()); 
-  _dPsiDx = new fullMatrix<double> ( _integration->size1()*3,nbPsi*e.size());
+  int nbNodes = _fs.coefficients.size1();
+  int nbQP = _integration->size1();
+  _redistributionFluxes[0] = new fullMatrix<double> (nbNodes, nbQP);
+  _redistributionFluxes[1] = new fullMatrix<double> (nbNodes, nbQP);
+  _redistributionFluxes[2] = new fullMatrix<double> (nbNodes, nbQP);
+  _redistributionSource = new fullMatrix<double> (nbNodes, nbQP);
+  _collocation = new fullMatrix<double> (nbQP, nbNodes);
+  _mapping = new fullMatrix<double> (e.size(), 10 * nbQP);
+  _imass = new fullMatrix<double> (nbNodes,nbNodes*e.size()); 
+  _dPsiDx = new fullMatrix<double> ( nbQP*3,nbNodes*e.size());
   _elementVolume = new fullMatrix<double> (e.size(),1);
   _innerRadii = new fullMatrix<double> (e.size(),1);
   double g[256][3],f[256];
@@ -80,7 +81,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
   for (int i=0;i<_elements.size();i++){
     MElement *e = _elements[i];
     element_to_index[e] = i;
-    fullMatrix<double> imass(*_imass,nbPsi*i,nbPsi);
+    fullMatrix<double> imass(*_imass,nbNodes*i,nbNodes);
     (*_innerRadii)(i,0)=e->getInnerRadius();
     for (int j=0;j< _integration->size1() ; j++ ){
       _fs.f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
@@ -100,14 +101,14 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
       (*_mapping)(i,10*j + 6) = ijac[2][0];
       (*_mapping)(i,10*j + 7) = ijac[2][1];
       (*_mapping)(i,10*j + 8) = ijac[2][2];
-      for (int k=0;k<nbPsi;k++){ 
-        for (int l=0;l<nbPsi;l++) { 
+      for (int k=0;k<nbNodes;k++){ 
+        for (int l=0;l<nbNodes;l++) { 
           imass(k,l) += f[k]*f[l]*weight*detjac;
         }
         // (iQP*3+iXYZ , iFace*NPsi+iPsi)
-        (*_dPsiDx)(j*3  ,i*nbPsi+k) = g[k][0]*ijac[0][0]+g[k][1]*ijac[0][1]+g[k][2]*ijac[0][2];
-        (*_dPsiDx)(j*3+1,i*nbPsi+k) = g[k][0]*ijac[1][0]+g[k][1]*ijac[1][1]+g[k][2]*ijac[1][2];
-        (*_dPsiDx)(j*3+2,i*nbPsi+k) = g[k][0]*ijac[2][0]+g[k][1]*ijac[2][1]+g[k][2]*ijac[2][2];
+        (*_dPsiDx)(j*3  ,i*nbNodes+k) = g[k][0]*ijac[0][0]+g[k][1]*ijac[0][1]+g[k][2]*ijac[0][2];
+        (*_dPsiDx)(j*3+1,i*nbNodes+k) = g[k][0]*ijac[1][0]+g[k][1]*ijac[1][1]+g[k][2]*ijac[1][2];
+        (*_dPsiDx)(j*3+2,i*nbNodes+k) = g[k][0]*ijac[2][0]+g[k][1]*ijac[2][1]+g[k][2]*ijac[2][2];
       }
     }
     imass.invertInPlace();
@@ -115,20 +116,34 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
   // redistribution matrix
   // quadrature weight x parametric gradients in quadrature points
 
-  for (int j=0;j<_integration->size1();j++) {
-    _fs.df((*_integration)(j,0),
-	   (*_integration)(j,1),
-	   (*_integration)(j,2), g);
-    _fs.f((*_integration)(j,0),
-	   (*_integration)(j,1),
-	   (*_integration)(j,2), f);
-    const double weight = (*_integration)(j,3);
-    for (int k=0;k<_fs.coefficients.size1();k++){ 
-      (*_redistributionFluxes[0])(k,j) = g[k][0] * weight;
-      (*_redistributionFluxes[1])(k,j) = g[k][1] * weight;
-      (*_redistributionFluxes[2])(k,j) = g[k][2] * weight;
-      (*_redistributionSource)(k,j) = f[k] * weight;
-      (*_collocation)(j,k) = f[k];
+  _PsiDPsiDXi = fullMatrix<double> (_dimUVW*nbQP, nbNodes*nbNodes);
+  //reditribution of the diffusive jacobian dimUVW*dimUVW*nbIntegrationPoints x nbNodes*nbNodes
+  _dPsiDXDPsiDXi = fullMatrix<double> (_dimUVW*_dimUVW*nbQP, nbNodes*nbNodes);
+
+  for (int xi=0;xi<nbQP; xi++) {
+    _fs.df((*_integration)(xi,0),
+	   (*_integration)(xi,1),
+	   (*_integration)(xi,2), g);
+    _fs.f((*_integration)(xi,0),
+	   (*_integration)(xi,1),
+	   (*_integration)(xi,2), f);
+    const double weight = (*_integration)(xi,3);
+    for (int k=0;k<nbNodes; k++){ 
+      (*_redistributionFluxes[0])(k,xi) = g[k][0] * weight;
+      (*_redistributionFluxes[1])(k,xi) = g[k][1] * weight;
+      (*_redistributionFluxes[2])(k,xi) = g[k][2] * weight;
+      (*_redistributionSource)(k,xi) = f[k] * weight;
+      (*_collocation)(xi,k) = f[k];
+    }
+    for (int alpha=0; alpha<_dimUVW; alpha++) for (int beta=0; beta<_dimUVW; beta++) {
+      for (int i=0; i<nbNodes; i++) for (int j=0; j<nbNodes; j++) {
+        _dPsiDXDPsiDXi((alpha*_dimUVW+beta)*nbQP+xi, i*nbNodes+j) = g[i][alpha]*g[j][beta]*weight;
+      }
+    }
+    for (int alpha=0; alpha<_dimUVW; alpha++){
+      for (int i=0; i<nbNodes; i++) for (int j=0; j<nbNodes; j++) {
+        _dPsiDXDPsiDXi(alpha*nbQP+xi, i*nbNodes+j) = g[i][alpha]*f[j]*weight;
+      }
     }
   }
 }
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 987e86593a50aad9b08d9235ece254d93d44a126..dc22cfcfd2fad169b9272b620a7f28392be90ec3 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -73,6 +73,10 @@ class dgGroupOfElements {
   fullMatrix<double> *_imass;
   //
   fullMatrix<double> *_dPsiDx;
+  //redistributions of the convective jacobian diumUVW*nbIntegrationPoints x nbNodes*nbNodes
+  fullMatrix<double> _PsiDPsiDXi;
+  //reditribution of the diffusive jacobian dimUVW*dimUVW*nbIntegrationPoints x nbNodes*nbNodes
+  fullMatrix<double> _dPsiDXDPsiDXi;
   // dimension of the parametric space and of the real space
   // may be different if the domain is a surface in 3D (manifold)
   int _dimUVW, _dimXYZ;
@@ -91,7 +95,6 @@ protected:
   bool _multirateInnerBuffer;
   bool _multirateOuterBuffer;
 public:
-
   inline int getMultirateExponent() const {return _multirateExponent;}
   inline int getIsInnerMultirateBuffer() const {return _multirateInnerBuffer;}
   inline int getIsOuterMultirateBuffer() const {return _multirateOuterBuffer;}
@@ -109,6 +112,8 @@ 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> & getDPsiDXDPsiDXi() const {return _dPsiDXDPsiDXi;}
+  inline const fullMatrix<double> & getPsiDPsiDXi() const {return _PsiDPsiDXi;}
   inline double getElementVolume (int iElement)const {return (*_elementVolume)(iElement,0);}
   inline double getInnerRadius(int iElement)const {return (*_innerRadii)(iElement,0);}
   inline double getDetJ (int iElement, int iGaussPoint) const {return (*_mapping)(iElement, 10*iGaussPoint + 9);}
diff --git a/Solver/dgResidual.cpp b/Solver/dgResidual.cpp
index c8f34f56478c338ff6d994dfffbc9be741e099fa..0289ba87d33fcc91d9fbb0a21529ca1c137fd111 100644
--- a/Solver/dgResidual.cpp
+++ b/Solver/dgResidual.cpp
@@ -2,6 +2,7 @@
 #include "dgConservationLaw.h"
 #include "dgGroupOfElements.h"
 #include "function.h"
+#include "functionDerivator.h"
 #include "dgDofContainer.h"
 #include "MElement.h"
 
@@ -108,6 +109,153 @@ void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double
   }
 }
 
+void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullMatrix<double> &solution, fullMatrix<double> &residual, fullMatrix<double> &jacobian) 
+{
+  residual.scale(0);
+  _cacheMap->setNbEvaluationPoints(group.getNbIntegrationPoints());
+  _UVW.set(group.getIntegrationPointsMatrix());
+  // ----- 1 ----  get the solution at quadrature points
+  // ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints) 
+  fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
+  // ----- 1.2 --- multiply the solution by the collocation matrix
+  group.getCollocationMatrix().mult(solution  , solutionQP); 
+  // ----- 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(), _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 )};
+  fullMatrix<double> Source = fullMatrix<double> (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
+  // ----- 2.2 --- allocate parametric fluxes (computed in UVW coordinates) for all elements at all integration points
+  fullMatrix<double> Fuvw[3] = {fullMatrix<double> ( group.getNbIntegrationPoints(), group.getNbElements() * _nbFields),
+    fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * _nbFields),
+    fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * _nbFields)};
+  fullMatrix<double> fuvwe;
+  fullMatrix<double> source;
+  fullMatrix<double> dPsiDx,dofs;
+  int nColA = group.getDimUVW()*group.getNbIntegrationPoints();
+  int nColB = group.getDimUVW()*group.getDimUVW()*group.getNbIntegrationPoints();
+  fullMatrix<double> A (_nbFields*_nbFields, nColA*group.getNbElements());
+  fullMatrix<double> B (_nbFields*_nbFields, nColB*group.getNbElements());
+  fullMatrix<double> a, b;
+  functionDerivator *dDiffusiveFluxdGradU,*dConvectiveFluxdU;
+  if(_diffusiveFlux)
+    dDiffusiveFluxdGradU=new functionDerivator(*_diffusiveFlux,_gradientSolutionQPe,1e-6);
+  if(_convectiveFlux)
+    dConvectiveFluxdU=new functionDerivator(*_convectiveFlux,_solutionQPe,1e-6);
+  // ----- 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
+    _solutionQPe.setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
+    a.setAsProxy(A, iElement*nColA, nColA);
+    b.setAsProxy(B, iElement*nColB, nColB);
+
+    if(_gradientSolutionQPe.somethingDependOnMe()){
+      dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
+      dofs.setAsProxy(solution, _nbFields*iElement, _nbFields);
+      dPsiDx.mult(dofs, _gradientSolutionQPe.set());
+    }
+    _cacheElement.set(group.getElement(iElement));
+    if(_convectiveFlux || _diffusiveFlux) {
+      // ----- 2.3.3 --- compute fluxes in UVW coordinates
+      for (int iUVW=0;iUVW<group.getDimUVW();iUVW++) {
+        // ----- 2.3.3.1 --- get a proxy on the big local flux matrix
+        fuvwe.setAsProxy(Fuvw[iUVW], iElement*_nbFields,_nbFields);
+        // ----- 2.3.3.1 --- compute \vec{f}^{UVW} =\vec{f}^{XYZ} J^{-1} and multiply bt detJ
+        for (int iXYZ=0;iXYZ<group.getDimXYZ();iXYZ++) {
+          for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+            // get the right inv jacobian matrix dU/dX element
+            const double invJ = group.getInvJ (iElement, iPt, iUVW, iXYZ);
+            // get the detJ at this point
+            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<_nbFields;k++){
+              if(_convectiveFlux)
+                fuvwe(iPt,k) += ( (*_convectiveFlux)(iPt,iXYZ*_nbFields+k)) * factor;
+              if(_diffusiveFlux){
+                fuvwe(iPt,k) += ( (*_diffusiveFlux)(iPt,iXYZ*_nbFields+k)) * factor;
+              }
+            }
+          }
+        } 
+      }
+    }
+    if (_sourceTerm){
+      source.setAsProxy(Source, iElement*_nbFields,_nbFields);
+      for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+        const double detJ = group.getDetJ (iElement, iPt);
+        for (int k=0;k<_nbFields;k++){
+          source(iPt,k) = (*_sourceTerm)(iPt,k)*detJ;
+        }
+      }
+    }
+    int nQP = group.getNbIntegrationPoints();
+    int dim = group.getDimUVW();
+    if(_diffusiveFlux) {
+      dDiffusiveFluxdGradU->compute();
+      for (int alpha=0; alpha < dim; alpha++) for (int beta=0; beta < dim; beta++) {
+        for(int x=0;x <group.getDimXYZ();x++) for(int y=0;y <group.getDimXYZ();x++) {
+          for (int xi=0; xi <group.getNbIntegrationPoints(); xi++) {
+            const double invJx = group.getInvJ (iElement, xi, alpha, x);
+            const double invJy = group.getInvJ (iElement, xi, beta, y);
+            const double detJ = group.getDetJ (iElement, xi);
+            const double factor = invJx * invJy * detJ;
+            for (int k=0; k< _nbFields; k++) for (int l=0; l<_nbFields; l++) {
+              b(k*_nbFields+l,(alpha*dim+beta)*nQP+xi) = dDiffusiveFluxdGradU->get(k*dim+alpha,l*dim+beta,xi)*factor;
+            }
+          }
+        }
+      }
+    }
+    if(_convectiveFlux) {
+      dConvectiveFluxdU->compute();
+      for (int alpha=0; alpha < dim; alpha++) {
+        for(int x=0;x <group.getDimXYZ();x++) {
+          for (int xi=0; xi <group.getNbIntegrationPoints(); xi++){
+            const double invJ = group.getInvJ (iElement, xi, alpha, x);
+            const double detJ = group.getDetJ (iElement, xi);
+            const double factor = invJ * detJ;
+            for (int k=0; k< _nbFields; k++) for (int l=0; l<_nbFields; l++){
+              a(k*_nbFields+l,alpha*nQP+xi) = dConvectiveFluxdU->get(k*dim+alpha,l,xi)*factor;
+            }
+          }
+        }
+      }
+    }
+  }
+
+  // ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
+  if(_convectiveFlux || _diffusiveFlux){
+    for (int iUVW=0;iUVW<group.getDimUVW();iUVW++){
+      residual.gemm(group.getFluxRedistributionMatrix(iUVW),Fuvw[iUVW]);
+    }
+  }
+  if(_sourceTerm){
+    residual.gemm(group.getSourceRedistributionMatrix(),Source);
+  }
+  int nbNodes = group.getNbNodes();
+  fullMatrix<double> jacobianK (_nbFields*_nbFields,nbNodes*nbNodes);
+  if (_convectiveFlux) {
+    jacobianK.gemm(group.getPsiDPsiDXi(),B);
+  }
+  if (_convectiveFlux) {
+    jacobianK.gemm(group.getDPsiDXDPsiDXi(),A);
+  }
+  fullMatrix<double> jacobianE, jacobianKE;
+  for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
+    jacobianKE.setAsProxy(jacobianK, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
+    jacobianKE.setAsProxy(jacobian, iElement*_nbFields*nbNodes, _nbFields*nbNodes);
+    for (int k=0; k<_nbFields;k++) for (int l=0;l<_nbFields;l++) {
+      for(int i=0; i<nbNodes; i++) for(int j=0; j<nbNodes; j++) {
+        jacobianE(l*nbNodes+j, k*nbNodes+i) = jacobianKE(k*_nbFields+l, i*nbNodes+j);
+      }
+    }
+  }
+}
+
 void dgResidualVolume::computeAndMap1Group(dgGroupOfElements &group, dgDofContainer &solution, dgDofContainer &residual) 
 {
   compute1Group (group, solution.getGroupProxy(&group), residual.getGroupProxy(&group));
@@ -384,6 +532,10 @@ void dgResidual::registerBindings (binding *b)
   mb->setDescription("compute the residual for one group given a fullMatrix proxy to the solution relative to this group"); 
   mb->setArgNames("group", "solution", "residual", NULL);
 
+  mb = cb->addMethod("compute1GroupWithJacobian", &dgResidualVolume::compute1GroupWithJacobian);
+  mb->setDescription("compute the residual for one group given a fullMatrix proxy to the solution relative to this group"); 
+  mb->setArgNames("group", "solution", "residual", "jacobian", NULL);
+
   cb = b->addClass<dgResidualInterface>("dgResidualInterface");
   cb->setDescription("compute the residual for one dgGroupOfElements");
   mb = cb->addMethod("computeAndMap1Group",&dgResidualInterface::computeAndMap1Group);
diff --git a/Solver/dgResidual.h b/Solver/dgResidual.h
index ed7a1c678d60aeb55120aad7dda1c655dc5eb84e..47ab9226d2e93fb4f83377dd47c71d67754a8aa2 100644
--- a/Solver/dgResidual.h
+++ b/Solver/dgResidual.h
@@ -22,6 +22,7 @@ class dgResidualVolume {
   dgResidualVolume (const dgConservationLaw &claw);
   ~dgResidualVolume();
   void compute1Group(dgGroupOfElements &group, fullMatrix<double> &solution, fullMatrix<double> &residual);
+  void compute1GroupWithJacobian(dgGroupOfElements &group, fullMatrix<double> &solution, fullMatrix<double> &residual, fullMatrix<double> &jacobian);
   void computeAndMap1Group(dgGroupOfElements &group, dgDofContainer &solution, dgDofContainer &residual);
 };