diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index b38620923d176ddecf16215d771172c2d4cd1d6d..f268df7f20a0c0233715b9c53bfc4a1f26c908f1 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -126,20 +126,20 @@ void mlawNonLinearTVP::getModifiedMandel(const STensor3& C, const IPNonLinearTVP
     M = 0.5*(corKir + temp2);
 };
 
-void mlawNonLinearTVP::checkCommutavity(const STensor3& C, const STensor3& S, const IPNonLinearTVP *q1) const{
+void mlawNonLinearTVP::checkCommutavity(STensor3& commuteCheck,const STensor3& Ce, const STensor3& S, const IPNonLinearTVP *q1) const{
     
     const STensor3& M = q1->getConstRefToModifiedMandel();
     const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
     
-    STensor3 temp1, temp2, temp3, temp4, temp5, temp6;
-    STensorOperation::multSTensor3(C, S, temp1);
-    STensorOperation::multSTensor3(S, C, temp2);
-    STensorOperation::multSTensor3(corKir, C, temp3);
-    STensorOperation::multSTensor3(C, corKir, temp4);
+    STensor3 temp1, temp2, temp3, temp4, temp6;
+    STensorOperation::multSTensor3(Ce, S, temp1);
+    STensorOperation::multSTensor3(S, Ce, temp2);
+    STensorOperation::multSTensor3(corKir, Ce, temp3);
+    STensorOperation::multSTensor3(Ce, corKir, temp4);
     
     // Check S 
-    temp5 = temp1;
-    temp5 -= temp2;
+    commuteCheck = temp1;
+    commuteCheck -= temp2;
     
     // Check corKir (just for fun)
     temp6 = temp3;
@@ -1321,7 +1321,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
           P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
     }
 
-
+  // check Commutavity
+  static STensor3 Check;
+  checkCommutavity(Check,Ce,S,q1);
+  
   // defo energy
   // q1->getRefToElasticEnergy()=deformationEnergy(*q1);
   q1->getRefToViscousEnergyPart()=viscousEnergy(*q0,*q1)+q0->getConstRefToViscousEnergyPart();
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.h b/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
index aa906beadf2fd1cb88b648374e529ecb334de3c2..060ed56a478b7cf903d2f5d2648a015fb1ae4ee0 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
@@ -32,7 +32,7 @@ class mlawNonLinearTVP : public mlawNonLinearTVE{
         // NEW
         virtual void hardening(const IPNonLinearTVP* q0, IPNonLinearTVP* q1, const double& T) const;
         virtual void getModifiedMandel(const STensor3& C, const IPNonLinearTVP *q0_, IPNonLinearTVP *q1) const;
-        virtual void checkCommutavity(const STensor3& C, const STensor3& S, const IPNonLinearTVP *q1) const;
+        virtual void checkCommutavity(STensor3& commuteCheck, const STensor3& Ce, const STensor3& S, const IPNonLinearTVP *q1) const;
         virtual void getChabocheCoeffs(fullVector<double>& coeffs, const double& opt, const IPNonLinearTVP *q1) const;
         virtual void getPlasticPotential(const double& phiP, const double& phiE, double& Px) const;
         virtual void getYieldCoefficientTempDers(const IPNonLinearTVP *q1, const double& T, fullVector<double>& DcoeffsDT, fullVector<double>& DDcoeffsDTT) const;