diff --git a/NonLinearSolver/materialLaw/mlaw.h b/NonLinearSolver/materialLaw/mlaw.h
index dfb67d87ed586c7e1e9031b6a259dbb6a72a80c3..3412224798ab340fbacf7d8d862a73a4b31e2557 100644
--- a/NonLinearSolver/materialLaw/mlaw.h
+++ b/NonLinearSolver/materialLaw/mlaw.h
@@ -399,21 +399,21 @@ class nonLocalDamageThermoMechanicalMaterialLaw : public nonLocalDamageMaterialL
   			    std::vector<STensor3>    &dStressDNonLocalVariable,
                             fullMatrix<double>       &dLocalVariableDNonLocalVariable,
                             std::vector<STensor3>    &characteristicLength,
-                            STensor3 &dPdT, STensor33 &dPdgradT,  fullVector<double> &dLocalVarDT,
+                            STensor3 &dPdT, STensor33 &dPdgradT,  fullMatrix<double> &dLocalVarDT,
                        
                             const double &T0, const double &Tn, 
                             const SVector3 &gradT0, const SVector3 &gradT,
                               
                             SVector3 &fluxT,
-                            STensor33 &dfluxTdF, std::vector<SVector3> &dfluxTdNonLocalVar,
+                            STensor33 &dfluxTdF, std::vector< std::vector<SVector3> > &dfluxTdNonLocalVar,
                             SVector3  &dfluxTdT, STensor3 &dfluxTdgradT,
                               
                             double &thermalSource, 
-                            STensor3 &dthermalSourcedF, fullVector<double> &dthermalSourcedNonLocalVar,
+                            STensor3 &dthermalSourcedF, fullMatrix<double> &dthermalSourcedNonLocalVar,
                             double &dthermalSourcedT, SVector3 &dthermalSourcedgradT,
          
                             double &mechanicalSource,
-                            STensor3 &dmechanicalSourceF, fullVector<double> &dmechanicalSourcedNonLocalVar,
+                            STensor3 &dmechanicalSourceF, fullMatrix<double> &dmechanicalSourcedNonLocalVar,
                             double &dmechanicalSourcedT, SVector3 &dmechanicalSourcedgradT,
 			    const bool stiff,
                             STensor43* elasticTangent = NULL,
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVEVPwithFailure.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVEVPwithFailure.cpp
index 0a8f03f38fcfa7ff1906c59c88c29242e3b2048e..90da9050b3d1237246bd46c284d1315604b6b33a 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVEVPwithFailure.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVEVPwithFailure.cpp
@@ -243,21 +243,21 @@ void mlawNonLocalDamagaNonLinearTVEVP::constitutive(const STensor3& F0, const ST
                               std::vector<STensor3>    &characteristicLength,
 
                               STensor3& dPdT, STensor33& dPdgradT,
-                              fullVector<double>& dLocalPlasticStrainDT,
+                              fullMatrix<double>& dLocalPlasticStrainDT,
 
                               const double& T0, const double& T,
                               const SVector3& gradT0, const SVector3& gradT,
 
                               SVector3& fluxT,
-                              STensor33& dfluxTdF, std::vector<SVector3>& dfluxTdNonLocalPlasticStrain,
+                              STensor33& dfluxTdF, std::vector< std::vector<SVector3> >& dfluxTdNonLocalPlasticStrain,
                               SVector3& dfluxTdT, STensor3& dfluxTdgradT,
 
                               double& thermalSource,
-                              STensor3& dthermalSourcedF, fullVector<double>& dthermalSourcedNonLocalPlasticStrain,
+                              STensor3& dthermalSourcedF, fullMatrix<double>& dthermalSourcedNonLocalPlasticStrain,
                               double& dthermalSourcedT, SVector3& dthermalSourcedgradT,
 
                               double& mechanicalSource,
-                              STensor3& dmechanicalSourcedF, fullVector<double>& dmechanicalSourcedNonLocalPlasticStrain,
+                              STensor3& dmechanicalSourcedF, fullMatrix<double>& dmechanicalSourcedNonLocalPlasticStrain,
                               double& dmechanicalSourcedT, SVector3& dmechanicalSourcedgradT,
 
                               const bool stiff,
@@ -268,7 +268,7 @@ void mlawNonLocalDamagaNonLinearTVEVP::constitutive(const STensor3& F0, const ST
 
   if(elasticTangent != NULL) Msg::Error("mlawNonLocalDamagaNonLinearTVEVP elasticTangent not defined");
   if(_tangentByPerturbation) Msg::Error("mlawNonLocalDamagaNonLinearTVEVP numerical tangent not implemented. Switching to algorithmic tangent.");
-  if(dTangent && dCalgdeps!=NULL)
+  if(dTangentdeps && dCalgdeps!=NULL)
   {
     Msg::Error("mlawNonLocalDamagaNonLinearTVEVP::constitutive: dCalgdeps has to be checked");
   }
@@ -398,19 +398,19 @@ void mlawNonLocalDamagaNonLinearTVEVP::constitutive(const STensor3& F0, const ST
 
 
     // partial p/partial F
-    dLocalPlasticStrainDStrain = DgammadF;
+    dLocalPlasticStrainDStrain[0] = DgammadF;
 
-    dLocalPlasticStrainDT = DgammadT;
+    dLocalPlasticStrainDT(0,0) = DgammadT;
     // -hat{P} partial D/partial tilde p
     // if (q1->getNonLocalToLocal()){   // CHECK what this does
       // STensorOperation::zero(dStressDNonLocalPlasticStrain);
     // }
     // else{
-    dStressDNonLocalPlasticStrain = Peff*(-1*q1->getDDamageDp());
+    dStressDNonLocalPlasticStrain[0] = Peff*(-1*q1->getDDamageDp());
     // }
 
     // partial p partial tilde p (0 if no MFH)
-    dLocalPlasticStrainDNonLocalPlasticStrain = 0.;
+    dLocalPlasticStrainDNonLocalPlasticStrain(0,0) = 0.;
 
     dPdT=dPeffdT;
     dPdT*=(1.-D);
@@ -421,9 +421,9 @@ void mlawNonLocalDamagaNonLinearTVEVP::constitutive(const STensor3& F0, const ST
     // we assume dDammage dT=0 // CHECK AND REFORMULATE
 
     STensorOperation::zero(dPdgradT);
-    STensorOperation::zero(dfluxTdNonLocalPlasticStrain);
-    dthermalSourcedNonLocalPlasticStrain=0.;
-    dmechanicalSourcedNonLocalPlasticStrain=0.;
+    STensorOperation::zero(dfluxTdNonLocalPlasticStrain[0][0]);
+    dthermalSourcedNonLocalPlasticStrain(0,0)=0.;
+    dmechanicalSourcedNonLocalPlasticStrain(0,0)=0.;
 
     // Plastic power stuff
     // STensor3& DplasticPowerDF = q1->getRefToDPlasticPowerDF(); // CHECK if exists or ADD
@@ -521,11 +521,8 @@ mlawLocalDamageNonLinearTVEVPWithFailure& mlawLocalDamageNonLinearTVEVPWithFailu
   const mlawLocalDamageNonLinearTVEVPWithFailure* src = dynamic_cast<const mlawLocalDamageNonLinearTVEVPWithFailure*>(&source);
   if(src != NULL)
   {
-    if(damLaw != NULL)
-    {
-      for (int i=0; i< damLaw.size(); i++){
-        if (damLaw[i]!=0) delete damLaw[i];
-      }
+    for (int i=0; i< damLaw.size(); i++){
+        if (damLaw[i]!=NULL) delete damLaw[i];
     }
     damLaw.clear();
     for (int i=0; i< src->damLaw.size(); i++){
@@ -828,12 +825,9 @@ mlawNonLocalDamageNonLinearTVEVPWithFailure& mlawNonLocalDamageNonLinearTVEVPWit
     _blockHardeningAfterFailure=src->_blockHardeningAfterFailure;
     _blockDamageAfterFailure=src->_blockDamageAfterFailure;
      
-    if(cLLaw != NULL)
-    {
-      for (int i=0; i< cLLaw.size(); i++){
-        if (cLLaw[i]!=0) delete cLLaw[i];
-        if (damLaw[i]!=0) delete damLaw[i];
-      }
+    for (int i=0; i< cLLaw.size(); i++){
+      if (cLLaw[i]!=0) delete cLLaw[i];
+      if (damLaw[i]!=0) delete damLaw[i];
     }
     cLLaw.clear();
     damLaw.clear();
@@ -946,24 +940,24 @@ void mlawNonLocalDamageNonLinearTVEVPWithFailure::predictorCorrector(const STens
                               std::vector<STensor3>    &characteristicLength,
 
                               STensor3& dPdT, STensor33& dPdgradT,
-                              fullVector<double>& dLocalVariableDT,
+                              fullMatrix<double>& dLocalVariableDT,
 
                               const double& T0, const double& T,
                               const SVector3& gradT0, const SVector3& gradT,
 
                               SVector3& fluxT,
                               STensor33& dfluxTdF,
-                              std::vector<SVector3>& dfluxTdNonLocalVariable,
+                              std::vector< std::vector<SVector3> >& dfluxTdNonLocalVariable,
                               SVector3& dfluxTdT, STensor3& dfluxTdgradT,
 
                               double& thermalSource,
                               STensor3& dthermalSourcedF,
-                              fullVector<double>& dthermalSourcedNonLocalVariable,
+                              fullMatrix<double>& dthermalSourcedNonLocalVariable,
                               double& dthermalSourcedT, SVector3& dthermalSourcedgradT,
 
                               double& mechanicalSource,
                               STensor3& dmechanicalSourcedF,
-                              fullVector<double>& dmechanicalSourcedNonLocalVariable,
+                              fullMatrix<double>& dmechanicalSourcedNonLocalVariable,
                               double& dmechanicalSourcedT, SVector3& dmechanicalSourcedgradT,
 
                               const bool stiff,
@@ -1131,8 +1125,8 @@ void mlawNonLocalDamageNonLinearTVEVPWithFailure::predictorCorrector(const STens
     // dLocalVariableDT = dgFdT;
     // IS THIS NEEDED ?? CHECK
     dLocalVariableDT.setAll(0.); // CHECK
-    dLocalVariableDT(0) = DgammadT;
-    dLocalVariableDT(1) = dgFdT;
+    dLocalVariableDT(0,0) = DgammadT;
+    dLocalVariableDT(1,0) = dgFdT;
 
 
     dPdT=dPeffdT;
@@ -1144,8 +1138,8 @@ void mlawNonLocalDamageNonLinearTVEVPWithFailure::predictorCorrector(const STens
     // we assume dDammage dT=0 // CHECK and REFORMULATE
 
     STensorOperation::zero(dPdgradT);
-    STensorOperation::zero(dfluxTdNonLocalVariable[0]);
-    STensorOperation::zero(dfluxTdNonLocalVariable[1]);
+    STensorOperation::zero(dfluxTdNonLocalVariable[0][0]);
+    STensorOperation::zero(dfluxTdNonLocalVariable[0][1]);
     dthermalSourcedNonLocalVariable.setAll(0.);
     dmechanicalSourcedNonLocalVariable.setAll(0.);
 
@@ -1261,24 +1255,24 @@ void mlawNonLocalDamageNonLinearTVEVPWithFailure::constitutive(const STensor3& F
                               std::vector<STensor3>    &characteristicLength,
 
                               STensor3& dPdT, STensor33& dPdgradT,
-                              fullVector<double>& dLocalVariableDT,
+                              fullMatrix<double>& dLocalVariableDT,
 
                               const double& T0, const double& T,
                               const SVector3& gradT0, const SVector3& gradT,
 
                               SVector3& fluxT,
                               STensor33& dfluxTdF,
-                              std::vector<SVector3>& dfluxTdNonLocalVariable,
+                              std::vector< std::vector<SVector3> >& dfluxTdNonLocalVariable,
                               SVector3& dfluxTdT, STensor3& dfluxTdgradT,
 
                               double& thermalSource,
                               STensor3& dthermalSourcedF,
-                              fullVector<double>& dthermalSourcedNonLocalVariable,
+                              fullMatrix<double>& dthermalSourcedNonLocalVariable,
                               double& dthermalSourcedT, SVector3& dthermalSourcedgradT,
 
                               double& mechanicalSource,
                               STensor3& dmechanicalSourcedF,
-                              fullVector<double>& dmechanicalSourcedNonLocalVariable,
+                              fullMatrix<double>& dmechanicalSourcedNonLocalVariable,
                               double& dmechanicalSourcedT, SVector3& dmechanicalSourcedgradT,
 
                               const bool stiff,            // if true compute the tangents
@@ -1288,7 +1282,7 @@ void mlawNonLocalDamageNonLinearTVEVPWithFailure::constitutive(const STensor3& F
 {
 
   if (!_tangentByPerturbation){
-  	predictorCorrector(F0,F1,P,q0_,q1_,Tangent,dLocalVariableDStrain,dStressDNonLocalVariable,dLocalVariableDNonLocalVariable,dPdT,dPdgradT,dLocalVariableDT,
+  	predictorCorrector(F0,F1,P,q0_,q1_,Tangent,nonLocalVariable0,nonLocalVariable,localVariable,  dLocalVariableDStrain,dStressDNonLocalVariable,dLocalVariableDNonLocalVariable,characteristicLength,dPdT,dPdgradT,dLocalVariableDT,
                             T0,T,gradT0,gradT,fluxT,dfluxTdF,dfluxTdNonLocalVariable,dfluxTdT,dfluxTdgradT,
                  	    thermalSource,dthermalSourcedF,dthermalSourcedNonLocalVariable,dthermalSourcedT,dthermalSourcedgradT,
           		    mechanicalSource,dmechanicalSourcedF,dmechanicalSourcedNonLocalVariable,dmechanicalSourcedT,dmechanicalSourcedgradT,
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVEVPwithFailure.h b/NonLinearSolver/materialLaw/mlawNonLinearTVEVPwithFailure.h
index eaee72d406ead8d64ba6bc9f8d81c6dabb919c51..72fff9cad6a00c077518951d0af7c0a9ba824085 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVEVPwithFailure.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVEVPwithFailure.h
@@ -112,21 +112,21 @@ class mlawNonLocalDamagaNonLinearTVEVP : public mlawNonLinearTVENonLinearTVP, pu
                               std::vector<STensor3>    &characteristicLength,
  
                               STensor3& dPdT, STensor33& dPdgradT,
-                              fullVector<double>& dLocalPlasticStrainDT,
+                              fullMatrix<double>& dLocalPlasticStrainDT,
 
                               const double& T0, const double& T,
                               const SVector3& gradT0, const SVector3& gradT,
 
                               SVector3& fluxT,
-                              STensor33& dfluxTdF,  std::vector<STensor3>& dfluxTdNonLocalPlasticStrain,
+                              STensor33& dfluxTdF,  std::vector< std::vector<SVector3> >& dfluxTdNonLocalPlasticStrain,
                               SVector3& dfluxTdT, STensor3& dfluxTdgradT,
 
                               double& thermalSource,
-                              STensor3& dthermalSourcedF, fullVector<double>& dthermalSourcedNonLocalPlasticStrain,
+                              STensor3& dthermalSourcedF, fullMatrix<double>& dthermalSourcedNonLocalPlasticStrain,
                               double& dthermalSourcedT, SVector3& dthermalSourcedgradT,
 
                               double& mechanicalSource,
-                              STensor3& dmechanicalSourcedF, fullVector<double>& dmechanicalSourcedNonLocalPlasticStrain,
+                              STensor3& dmechanicalSourcedF, fullMatrix<double>& dmechanicalSourcedNonLocalPlasticStrain,
                               double& dmechanicalSourcedT, SVector3& dmechanicalSourcedgradT,
 
                               const bool stiff,
@@ -266,24 +266,24 @@ class mlawNonLocalDamageNonLinearTVEVPWithFailure : public mlawNonLinearTVEVPwit
                               std::vector<STensor3>    &characteristicLength,
                               
                               STensor3& dPdT, STensor33& dPdgradT,
-                              fullVector<double>& dLocalVariableDT,
+                              fullMatrix<double>& dLocalVariableDT,
 
                               const double& T0, const double& T,
                               const SVector3& gradT0, const SVector3& gradT,
 
                               SVector3& fluxT,
                               STensor33& dfluxTdF,
-                              std::vector<SVector3>& dfluxTdNonLocalVariable,
+                              std::vector< std::vector<SVector3> >& dfluxTdNonLocalVariable,
                               SVector3& dfluxTdT, STensor3& dfluxTdgradT,
 
                               double& thermalSource,
                               STensor3& dthermalSourcedF,
-                              fullVector<double>& dthermalSourcedNonLocalVariable,
+                              fullMatrix<double>& dthermalSourcedNonLocalVariable,
                               double& dthermalSourcedT, SVector3& dthermalSourcedgradT,
 
                               double& mechanicalSource,
                               STensor3& dmechanicalSourcedF,
-                              fullVector<double>& dmechanicalSourcedNonLocalVariable,
+                              fullMatrix<double>& dmechanicalSourcedNonLocalVariable,
                               double& dmechanicalSourcedT, SVector3& dmechanicalSourcedgradT,
 
                               const bool stiff,
@@ -314,24 +314,24 @@ class mlawNonLocalDamageNonLinearTVEVPWithFailure : public mlawNonLinearTVEVPwit
 			      std::vector<STensor3>    &characteristicLength,
  
                               STensor3& dPdT, STensor33& dPdgradT,
-                              fullVector<double>& dLocalVariableDT,
+                              fullMatrix<double>& dLocalVariableDT,
 
                               const double& T0, const double& T,
                               const SVector3& gradT0, const SVector3& gradT,
 
                               SVector3& fluxT,
                               STensor33& dfluxTdF,
-                              std::vector<SVector3>& dfluxTdNonLocalVariable,
+                              std::vector< std::vector<SVector3> >& dfluxTdNonLocalVariable,
                               SVector3& dfluxTdT, STensor3& dfluxTdgradT,
 
                               double& thermalSource,
                               STensor3& dthermalSourcedF,
-                              fullVector<double>& dthermalSourcedNonLocalVariable,
+                              fullMatrix<double>& dthermalSourcedNonLocalVariable,
                               double& dthermalSourcedT, SVector3& dthermalSourcedgradT,
 
                               double& mechanicalSource,
                               STensor3& dmechanicalSourcedF,
-                              fullVector<double>& dmechanicalSourcedNonLocalVariable,
+                              fullMatrix<double>& dmechanicalSourcedNonLocalVariable,
                               double& dmechanicalSourcedT, SVector3& dmechanicalSourcedgradT,
 
                               const bool stiff,            // if true compute the tangents
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2FullyCoupledThermoMechanics.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2FullyCoupledThermoMechanics.cpp
index 03a6c1698a31423a9d91dcca239c8956a101ad09..a84af4072bd0fb5f036bf71bf2dcb1f15a1c5083 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2FullyCoupledThermoMechanics.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2FullyCoupledThermoMechanics.cpp
@@ -170,21 +170,21 @@ void mlawNonLocalDamageJ2FullyCoupledThermoMechanics::constitutive(
   			      std::vector<STensor3>    &dStressDNonLocalPlasticStrain,
                               fullMatrix<double>       &dLocalPlasticStrainDNonLocalPlasticStrain,
                               std::vector<STensor3>    &characteristicLength,
-                              STensor3 &dPdT, STensor33 &dPdgradT,  fullVector<double> &dLocalPlasticStrainDT,
+                              STensor3 &dPdT, STensor33 &dPdgradT,  fullMatrix<double> &dLocalPlasticStrainDT,
                               
                               const double &T0, const double &Tn, 
                               const SVector3 &gradT0, const SVector3 &gradT,
                               
                               SVector3 &fluxT,
-                              STensor33 &dfluxTdF, std::vector<SVector3> &dfluxTdNonLocalPlasticStrain,
+                              STensor33 &dfluxTdF, std::vector< std::vector<SVector3> > &dfluxTdNonLocalPlasticStrain,
                               SVector3  &dfluxTdT, STensor3 &dfluxTdgradT,
                               
                               double &thermalSource, 
-                              STensor3 &dthermalSourcedF, fullVector<double> &dthermalSourcedNonLocalPlasticStrain,
+                              STensor3 &dthermalSourcedF, fullMatrix<double> &dthermalSourcedNonLocalPlasticStrain,
                               double &dthermalSourcedT, SVector3 &dthermalSourcedgradT,
           
                               double &mechanicalSource,
-                              STensor3 &dmechanicalSourceF, fullVector<double> &dmechanicalSourcedNonLocalPlasticStrain,
+                              STensor3 &dmechanicalSourceF, fullMatrix<double> &dmechanicalSourcedNonLocalPlasticStrain,
                               double &dmechanicalSourcedT, SVector3 &dmechanicalSourcedgradT,
                               const bool stiff,
                               STensor43* elasticTangent,
@@ -278,21 +278,21 @@ void mlawNonLocalDamageJ2FullyCoupledThermoMechanics::predictorCorrector(
                               fullMatrix<double>       &dLocalPlasticStrainDNonLocalPlasticStrain,
                               std::vector<STensor3>    &characteristicLength,
                               const bool stiff,
-                              STensor3 &dPdT, STensor33 &dPdgradT,  fullVector<double> &dLocalPlasticStrainDT,
+                              STensor3 &dPdT, STensor33 &dPdgradT,  fullMatrix<double> &dLocalPlasticStrainDT,
                               
                               const double &T0, const double &Tn, 
                               const SVector3 &gradT0, const SVector3 &gradT,
                               
                               SVector3 &fluxT,
-                              STensor33 &dfluxTdF, std::vector<SVector3> &dfluxTdNonLocalPlasticStrain,
+                              STensor33 &dfluxTdF, std::vector< std::vector<SVector3> > &dfluxTdNonLocalPlasticStrain,
                               SVector3  &dfluxTdT, STensor3 &dfluxTdgradT,
                               
                               double &thermalSource, 
-                              STensor3 &dthermalSourcedF, fullVector<double> &dthermalSourcedNonLocalPlasticStrain,
+                              STensor3 &dthermalSourcedF, fullMatrix<double> &dthermalSourcedNonLocalPlasticStrain,
                               double &dthermalSourcedT, SVector3 &dthermalSourcedgradT,
           
                               double &mechanicalSource,
-                              STensor3 &dmechanicalSourcedF, fullVector<double> &dmechanicalSourcedNonLocalPlasticStrain,
+                              STensor3 &dmechanicalSourcedF, fullMatrix<double> &dmechanicalSourcedNonLocalPlasticStrain,
                               double &dmechanicalSourcedT, SVector3 &dmechanicalSourcedgradT,
                               STensor43* elasticTangent  
                            ) const
@@ -313,7 +313,7 @@ void mlawNonLocalDamageJ2FullyCoupledThermoMechanics::predictorCorrector(
   
   
   mlawJ2FullyCoupledThermoMechanics::predictorCorector(F0,Fn,Peff,ipvprev,ipvcur,Tangent,dFpdF,dFpdT,  
-                      dFedF,dFedT,dLocalPlasticStrainDStrain[0],dLocalPlasticStrainDT[0],T0,Tn,gradT0,gradT,fluxT,dPeffdT,dfluxTdgradT,dfluxTdT,dfluxTdF,
+                      dFedF,dFedT,dLocalPlasticStrainDStrain[0],dLocalPlasticStrainDT(0,0),T0,Tn,gradT0,gradT,fluxT,dPeffdT,dfluxTdgradT,dfluxTdT,dfluxTdF,
                       thermalSource,dthermalSourcedT,dthermalSourcedF,
                       mechanicalSource,dmechanicalSourcedT,dmechanicalSourcedF,stiff);
 
@@ -550,7 +550,7 @@ void mlawNonLocalDamageJ2FullyCoupledThermoMechanics::predictorCorrector(
 
     
     STensorOperation::zero(dPdgradT); 
-    STensorOperation::zero(dfluxTdNonLocalPlasticStrain[0]);
+    STensorOperation::zero(dfluxTdNonLocalPlasticStrain[0][0]);
     dthermalSourcedNonLocalPlasticStrain(0,0)=0.;
     dmechanicalSourcedNonLocalPlasticStrain(0,0)=0.;
 
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2FullyCoupledThermoMechanics.h b/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2FullyCoupledThermoMechanics.h
index f5f491324b65d8a098044fa3d523aedf6eebda95..469e9c7c010ad4aafe5486be9670e03e84eb879d 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2FullyCoupledThermoMechanics.h
+++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2FullyCoupledThermoMechanics.h
@@ -102,21 +102,21 @@ class mlawNonLocalDamageJ2FullyCoupledThermoMechanics : public mlawJ2FullyCouple
   			      std::vector<STensor3>    &dStressDNonLocalPlasticStrain,
                               fullMatrix<double>       &dLocalPlasticStrainDNonLocalPlasticStrain,
                               std::vector<STensor3>    &characteristicLength,
-                              STensor3 &dPdT, STensor33 &dPdgradT,  fullVector<double> &dLocalVarDT,
+                              STensor3 &dPdT, STensor33 &dPdgradT,  fullMatrix<double> &dLocalVarDT,
                        
                               const double &T0, const double &Tn, 
                               const SVector3 &gradT0, const SVector3 &gradT,
                               
                               SVector3 &fluxT,
-                              STensor33 &dfluxTdF, std::vector<SVector3> &dfluxTdNonLocalVar,
+                              STensor33 &dfluxTdF, std::vector< std::vector<SVector3> > &dfluxTdNonLocalVar,
                               SVector3  &dfluxTdT, STensor3 &dfluxTdgradT,
                               
                               double &thermalSource, 
-                              STensor3 &dthermalSourcedF, fullVector<double> &dthermalSourcedNonLocalVar,
+                              STensor3 &dthermalSourcedF, fullMatrix<double> &dthermalSourcedNonLocalVar,
                               double &dthermalSourcedT, SVector3 &dthermalSourcedgradT,
           
                               double &mechanicalSource,
-                              STensor3 &dmechanicalSourceF, fullVector<double> &dmechanicalSourcedNonLocalVar,
+                              STensor3 &dmechanicalSourceF, fullMatrix<double> &dmechanicalSourcedNonLocalVar,
                               double &dmechanicalSourcedT, SVector3 &dmechanicalSourcedgradT,
                               const bool stiff,
                               STensor43* elasticTangent = NULL,
@@ -140,21 +140,21 @@ class mlawNonLocalDamageJ2FullyCoupledThermoMechanics : public mlawJ2FullyCouple
                               fullMatrix<double>       &dLocalPlasticStrainDNonLocalPlasticStrain,
                               std::vector<STensor3>    &characteristicLength,
                               const bool stiff,
-                              STensor3 &dPdT, STensor33 &dPdgradT,  fullVector<double> &dLocalPlasticStrainDT,
+                              STensor3 &dPdT, STensor33 &dPdgradT,  fullMatrix<double> &dLocalPlasticStrainDT,
                            
                               const double &T0, const double &Tn, 
                               const SVector3 &gradT0, const SVector3 &gradT,
                               
                               SVector3 &fluxT,
-                              STensor33 &dfluxTdF, std::vector<SVector3> &dfluxTdNonLocalPlasticStrain,
+                              STensor33 &dfluxTdF, std::vector< std::vector<SVector3> > &dfluxTdNonLocalPlasticStrain,
                               SVector3  &dfluxTdT, STensor3 &dfluxTdgradT,
                               
                               double &thermalSource, 
-                              STensor3 &dthermalSourcedF, fullVector<double> &dthermalSourcedNonLocalPlasticStrain,
+                              STensor3 &dthermalSourcedF, fullMatrix<double> &dthermalSourcedNonLocalPlasticStrain,
                               double &dthermalSourcedT, SVector3 &dthermalSourcedgradT,
           
                               double &mechanicalSource,
-                              STensor3 &dmechanicalSourceF, fullVector<double> &dmechanicalSourcedNonLocalPlasticStrain,
+                              STensor3 &dmechanicalSourceF, fullMatrix<double> &dmechanicalSourcedNonLocalPlasticStrain,
                               double &dmechanicalSourcedT, SVector3 &dmechanicalSourcedgradT,
                               STensor43* elasticTangent 
                            ) const;
diff --git a/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp b/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp
index 2baa463d0f22c6d1f3affea8b7e302e193c61b81..46da8d4f678287d1127b0400867994f5227026f0 100644
--- a/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp
+++ b/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp
@@ -940,24 +940,25 @@ void NonLocalDamageNonLinearTVEVPDG3DMaterialLaw::stress(IPVariable* ipv, const
 
   /* compute stress */
   _nldNonLinearTVEVP->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
-
-                        ipvcur->getRefToDLocalVariableDStrain()[0],ipvcur->getRefToDStressDNonLocalVariable()[0],
-                        ipvcur->getRefToDLocalVariableDNonLocalVariable()(0,0),
-
-                        ipvcur->getRefTodPdField()[0], ipvcur->getRefTodPdGradField()[0],  ipvcur->getRefTodLocalVariableDExtraDofDiffusionField()(0,0),
+                        ipvprev->getConstRefToNonLocalVariable(), ipvcur->getConstRefToNonLocalVariable(), ipvcur->getRefToLocalVariable(),
+                        ipvcur->getRefToDLocalVariableDStrain(),ipvcur->getRefToDStressDNonLocalVariable(),
+                        ipvcur->getRefToDLocalVariableDNonLocalVariable(),
+                        ipvcur->getRefToCharacteristicLengthMatrix(),
+                        
+                        ipvcur->getRefTodPdField()[0], ipvcur->getRefTodPdGradField()[0],  ipvcur->getRefTodLocalVariableDExtraDofDiffusionField(),
                         ipvprev->getConstRefToField(0),ipvcur->getConstRefToField(0),
                         ipvprev->getConstRefToGradField()[0],ipvcur->getConstRefToGradField()[0],
 
                         ipvcur->getRefToFlux()[0],
-                        ipvcur->getRefTodFluxdF()[0], ipvcur->getRefTodFluxdNonLocalVariable()[0][0],
+                        ipvcur->getRefTodFluxdF()[0], ipvcur->getRefTodFluxdNonLocalVariable(),
                         ipvcur->getRefTodFluxdField()[0][0], ipvcur->getRefTodFluxdGradField()[0][0],
 
                         ipvcur->getRefToFieldSource()(0),
-                        ipvcur->getRefTodFieldSourcedF()[0], ipvcur->getRefTodFieldSourcedNonLocalVariable()(0,0),
+                        ipvcur->getRefTodFieldSourcedF()[0], ipvcur->getRefTodFieldSourcedNonLocalVariable(),
                         ipvcur->getRefTodFieldSourcedField()(0,0), ipvcur->getRefTodFieldSourcedGradField()[0][0],
 
                         ipvcur->getRefToMechanicalSource()(0),
-                        ipvcur->getRefTodMechanicalSourcedF()[0],ipvcur->getRefTodMechanicalSourcedNonLocalVariable()(0,0),
+                        ipvcur->getRefTodMechanicalSourcedF()[0],ipvcur->getRefTodMechanicalSourcedNonLocalVariable(),
                         ipvcur->getRefTodMechanicalSourcedField()(0,0), ipvcur->getRefTodMechanicalSourcedGradField()[0][0],
 
                         stiff,&elasticL,dTangent,ipvcur->getPtrTodCmdFm());
@@ -1584,9 +1585,11 @@ void NonLocalDamageNonLinearTVEVPDG3DMaterialLawWithFailure::stress(IPVariable*
 
   /* compute stress */
   _nldNonLinearTVEVP->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
-
+  
+                        ipvprev->getConstRefToNonLocalVariable(), ipvcur->getConstRefToNonLocalVariable(), ipvcur->getRefToLocalVariable(),
                         ipvcur->getRefToDLocalVariableDStrain(),ipvcur->getRefToDStressDNonLocalVariable(),
                         ipvcur->getRefToDLocalVariableDNonLocalVariable(),
+                        ipvcur->getRefToCharacteristicLengthMatrix(),
 
                         ipvcur->getRefTodPdField()[0], ipvcur->getRefTodPdGradField()[0],  ipvcur->getRefTodLocalVariableDExtraDofDiffusionField(),
 
@@ -1594,7 +1597,7 @@ void NonLocalDamageNonLinearTVEVPDG3DMaterialLawWithFailure::stress(IPVariable*
                         ipvprev->getConstRefToGradField()[0],ipvcur->getConstRefToGradField()[0],
 
                         ipvcur->getRefToFlux()[0],
-                        ipvcur->getRefTodFluxdF()[0], ipvcur->getRefTodFluxdNonLocalVariable()[0],
+                        ipvcur->getRefTodFluxdF()[0], ipvcur->getRefTodFluxdNonLocalVariable(),
                         ipvcur->getRefTodFluxdField()[0][0], ipvcur->getRefTodFluxdGradField()[0][0],
 
                         ipvcur->getRefToFieldSource()(0),
diff --git a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
index 5c49dc603213ae945702235fea47ac37a8d478ed..04add7ee1bdb85f3071fb283da30c3e85485a47c 100644
--- a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
+++ b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
@@ -2509,19 +2509,18 @@ void NonLocalDamageJ2FullyCoupledDG3DMaterialLaw::stress(IPVariable* ipv, const
                               ipvprev->getConstRefToNonLocalVariable(), ipvcur->getConstRefToNonLocalVariable(), ipvcur->getRefToLocalVariable(),
                               ipvcur->getRefToDLocalVariableDStrain(), ipvcur->getRefToDStressDNonLocalVariable(),
                               ipvcur->getRefToDLocalVariableDNonLocalVariable(), ipvcur->getRefToCharacteristicLengthMatrix(),
-                              stiff,
-                              ipvcur->getRefTodPdField()[0], ipvcur->getRefTodPdGradField()[0],  ipvcur->getRefTodLocalVariableDExtraDofDiffusionField()(0,0),
+                              ipvcur->getRefTodPdField()[0], ipvcur->getRefTodPdGradField()[0],  ipvcur->getRefTodLocalVariableDExtraDofDiffusionField(),
                               ipvprev->getConstRefToField(0),ipvcur->getConstRefToField(0),
                               ipvprev->getConstRefToGradField()[0],ipvcur->getConstRefToGradField()[0],
                               ipvcur->getRefToFlux()[0],
-                              ipvcur->getRefTodFluxdF()[0], ipvcur->getRefTodFluxdNonLocalVariable()[0][0],
+                              ipvcur->getRefTodFluxdF()[0], ipvcur->getRefTodFluxdNonLocalVariable(),
                               ipvcur->getRefTodFluxdField()[0][0], ipvcur->getRefTodFluxdGradField()[0][0],
                               ipvcur->getRefToFieldSource()(0), 
-                              ipvcur->getRefTodFieldSourcedF()[0], ipvcur->getRefTodFieldSourcedNonLocalVariable()(0,0),
+                              ipvcur->getRefTodFieldSourcedF()[0], ipvcur->getRefTodFieldSourcedNonLocalVariable(),
                               ipvcur->getRefTodFieldSourcedField()(0,0), ipvcur->getRefTodFieldSourcedGradField()[0][0],  
                               ipvcur->getRefToMechanicalSource()(0),
-                              ipvcur->getRefTodMechanicalSourcedF()[0],ipvcur->getRefTodMechanicalSourcedNonLocalVariable()(0,0),
-                              ipvcur->getRefTodMechanicalSourcedField()(0,0), ipvcur->getRefTodMechanicalSourcedGradField()[0][0]);
+                              ipvcur->getRefTodMechanicalSourcedF()[0],ipvcur->getRefTodMechanicalSourcedNonLocalVariable(),
+                              ipvcur->getRefTodMechanicalSourcedField()(0,0), ipvcur->getRefTodMechanicalSourcedGradField()[0][0],stiff);
                             
                               
   ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);