From 2b2cfbad4a26e5cb7d5f0aeb0706943bcc98a45f Mon Sep 17 00:00:00 2001
From: FLE_Knight <ujwalkishore.jinaga@mail.polimi.it>
Date: Sat, 6 Jul 2024 17:59:53 +0200
Subject: [PATCH] Made changes to rotation, eigenvalue alignment functions by
 adding a boolean. segmentation fault in norm+dotProduct function

---
 .../materialLaw/STensorOperations.cpp         | 12 ++++---
 .../materialLaw/STensorOperations.h           |  4 +--
 .../mlawNonLinearTVENonLinearTVP.cpp          |  5 +--
 .../mlawNonLinearTVENonLinearTVP2.cpp         | 31 ++++++++++---------
 .../mlawNonLinearTVENonLinearTVP2.h           |  2 +-
 .../materialLaw/mlawNonLinearTVM.cpp          | 30 +++---------------
 .../materialLaw/mlawNonLinearTVM.h            |  6 +---
 7 files changed, 35 insertions(+), 55 deletions(-)

diff --git a/NonLinearSolver/materialLaw/STensorOperations.cpp b/NonLinearSolver/materialLaw/STensorOperations.cpp
index 8e423bd91..059751a9b 100644
--- a/NonLinearSolver/materialLaw/STensorOperations.cpp
+++ b/NonLinearSolver/materialLaw/STensorOperations.cpp
@@ -632,7 +632,7 @@ void STensorOperation::alignEigenDecomposition_NormBased(const STensor3& A, cons
 	a3_align = eigenValReal1_permuted(2);
 }
 
-void STensorOperation::alignEigenDecomposition_NormBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex){
+void STensorOperation::alignEigenDecomposition_NormBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, bool rotate){
 	// Align A to B, where B is fixed
 
 	// Initialise 1
@@ -686,6 +686,7 @@ void STensorOperation::alignEigenDecomposition_NormBased_indexOnly(const STensor
 	normamised_norm = norm_min; // /pow(A.dotprod(),0.5);
 	Msg::Info("normamised_norm of A = %f",normamised_norm);
 	if(normamised_norm > 1.e-5){
+		rotate = true;
 		alignedIndex = eigenValReal1_permuted_Index;
 	}
 }
@@ -764,7 +765,7 @@ void STensorOperation::alignEigenDecomposition_TangentBased_indexOnly(const STen
 	alignedIndex = eigenValReal1_permuted_Index;
 }
 
-void STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex){
+void STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, bool rotate){
 	// Align A to B, where B is fixed
 
 	// Initialise 1
@@ -823,7 +824,7 @@ void STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(
 
 	// Generate permutated eigVal vector
 	static fullVector<double> eigenValReal1_permuted(3);
-	static fullVector<int> eigenValReal1_permuted_Index(3);
+	static fullVector<int> eigenValReal1_permuted_Index(3); eigenValReal1_permuted_Index(0) = 0; eigenValReal1_permuted_Index(1) = 1; eigenValReal1_permuted_Index(2) = 2;
 	for (int p=0; p<3; p++)
 		for (int q=0; q<3; q++){
 			if(q!=p){
@@ -859,7 +860,10 @@ void STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(
 		}
 
 	// Assign the aligned indices
-	alignedIndex = eigenValReal1_permuted_Index;
+	if (eigenValReal1_permuted_Index(0) != 0 || eigenValReal1_permuted_Index(1) != 1 || eigenValReal1_permuted_Index(2) != 2){
+		rotate = true;
+		alignedIndex = eigenValReal1_permuted_Index;
+	}
 }
 
 void STensorOperation::VRDecomposition(const STensor3& a, STensor3& V, STensor3&R)
diff --git a/NonLinearSolver/materialLaw/STensorOperations.h b/NonLinearSolver/materialLaw/STensorOperations.h
index 92bd167fd..a60a00393 100644
--- a/NonLinearSolver/materialLaw/STensorOperations.h
+++ b/NonLinearSolver/materialLaw/STensorOperations.h
@@ -1933,9 +1933,9 @@ namespace STensorOperation{
 
       }
 
-  void alignEigenDecomposition_NormBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex);
+  void alignEigenDecomposition_NormBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, bool rotate = false);
 
-  void alignEigenDecomposition_NormAndDotProductBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex);
+  void alignEigenDecomposition_NormAndDotProductBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, bool rotate = false);
 
   void alignEigenDecomposition_NormBased(const STensor3& A, const STensor3& B,
 		  	  	  	  	  	  	  	  	  	  	const double& b1, const double& b2, const double& b3, const STensor3& B1, const STensor3& B2, const STensor3& B3,
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
index 8492c3f2c..1149ee3e2 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
@@ -481,6 +481,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
   
   // NEW EIGEN
   static fullVector<int> alignedIndex(3); alignedIndex(0) = 0; alignedIndex(1) = 1; alignedIndex(2) = 2;
+  bool rotate = false;
   static STensor3 Ee0_rot, R;
   static STensor43 dRdEe, dRtdEe, rotationStiffness_pr ,rotationStiffness;
   STensorOperation::zero(rotationStiffness_pr); STensorOperation::zero(rotationStiffness);
@@ -515,7 +516,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
 		  mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm,&Ge_TrEeTensor);
 	  }
 	  else{
-		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,Bd_stiffnessTerm,&Ge_TrEeTensor,&rotationStiffness); // NEW EIGEN
+		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,rotate,Bd_stiffnessTerm,&Ge_TrEeTensor,&rotationStiffness); // NEW EIGEN
 	  }
       Ge_TrEeTensor_pr = Ge_TrEeTensor;
       rotationStiffness_pr = rotationStiffness;
@@ -525,7 +526,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
 		  mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm);
 	  }
 	  else{
-	      mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,Bd_stiffnessTerm,NULL,&rotationStiffness); // NEW EIGEN
+	      mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,rotate,Bd_stiffnessTerm,NULL,&rotationStiffness); // NEW EIGEN
 	  }
 	  rotationStiffness_pr = rotationStiffness;
   }
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.cpp
index ecdf05275..6890f8141 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.cpp
@@ -482,6 +482,7 @@ void mlawNonLinearTVENonLinearTVP2::predictorCorrector_TVP_nonAssociatedFlow_non
 
   // NEW EIGEN
   static fullVector<int> alignedIndex(3); alignedIndex(0) = 0; alignedIndex(1) = 1; alignedIndex(2) = 2;
+  bool rotate = false;
   STensor3& R = q1->_R;
   STensor43& dRdEe = q1->_dRdEe;
   STensor43& dRtdEe = q1->_dRtdEe;
@@ -501,10 +502,10 @@ void mlawNonLinearTVENonLinearTVP2::predictorCorrector_TVP_nonAssociatedFlow_non
   }
   else if(_rotationCorrectionScheme == 2){
 	  Ee0_rot = q0->_Ee;
-	  STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex);
+	  // STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate);
 	  // STensorOperation::alignEigenDecomposition_EigenVectorDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex);
 	  // STensorOperation::alignEigenDecomposition_TangentBased_indexOnly(q0->_Ee,Ee,alignedIndex);
-	  // STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex);
+	  STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate);
   }
   // NEW EIGEN
 
@@ -515,7 +516,7 @@ void mlawNonLinearTVENonLinearTVP2::predictorCorrector_TVP_nonAssociatedFlow_non
 	  }
 	  else{
 		  // mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm,&Ge_TrEeTensor);
-		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,Bd_stiffnessTerm,&Ge_TrEeTensor,&rotationStiffness); // NEW EIGEN
+		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,rotate,Bd_stiffnessTerm,&Ge_TrEeTensor,&rotationStiffness); // NEW EIGEN
 	  }
       Ge_TrEeTensor_pr = Ge_TrEeTensor;
       rotationStiffness_pr = rotationStiffness;
@@ -526,7 +527,7 @@ void mlawNonLinearTVENonLinearTVP2::predictorCorrector_TVP_nonAssociatedFlow_non
 	  }
 	  else{
 		  // mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm);
-	      mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,Bd_stiffnessTerm,NULL,&rotationStiffness); // NEW EIGEN
+	      mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,rotate,Bd_stiffnessTerm,NULL,&rotationStiffness); // NEW EIGEN
 	  }
 	  rotationStiffness_pr = rotationStiffness;
   }
@@ -860,11 +861,11 @@ void mlawNonLinearTVENonLinearTVP2::predictorCorrector_TVP_nonAssociatedFlow_non
             if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                 // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
                 getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA,
-                    Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,alignedIndex,&Ge_TrEeTensor,&rotationStiffness);
+                    Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,alignedIndex,rotate,&Ge_TrEeTensor,&rotationStiffness);
             }
             else{
                 getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA,
-                    Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,alignedIndex,NULL,&rotationStiffness);
+                    Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,alignedIndex,rotate,NULL,&rotationStiffness);
             }
             PhiEq = sqrt(1.5*devPhi.dotprod());
 
@@ -2196,7 +2197,7 @@ void mlawNonLinearTVENonLinearTVP2::getIterated_DPhi(const double& T0, const dou
                                             double& Ke, double& Ge, STensor43& Ge_Tensor,
                                             double& ptilde, STensor3& devPhi,
                                             STensor3& Phi, STensor3& N, STensor3& expGN, STensor43& dexpAdA,
-                                            STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv, fullVector<int>& alignedIndex,
+                                            STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv, fullVector<int>& alignedIndex, bool rotate,
                                             STensor43* Ge_TrEeTensor, STensor43* rotationStiffness) const{
 
   // This function calculates Phi iteratively.
@@ -2247,10 +2248,10 @@ void mlawNonLinearTVENonLinearTVP2::getIterated_DPhi(const double& T0, const dou
   }
   else if(_rotationCorrectionScheme == 2){
 	  Ee0_rot = q0->_Ee;
-	  STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex);
+	  // STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate);
 	  // STensorOperation::alignEigenDecomposition_EigenVectorDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex);
 	  // STensorOperation::alignEigenDecomposition_TangentBased_indexOnly(q0->_Ee,Ee,alignedIndex);
-	  // STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex);
+	  STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate);
   }
   // NEW EIGEN
 
@@ -2259,7 +2260,7 @@ void mlawNonLinearTVENonLinearTVP2::getIterated_DPhi(const double& T0, const dou
 		  mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm);
 	  }
 	  else{
-		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,true,alignedIndex,Bd_stiffnessTerm,NULL,rotationStiffness); // NEW EIGEN
+		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,true,alignedIndex,rotate,Bd_stiffnessTerm,NULL,rotationStiffness); // NEW EIGEN
 	  }
   }
   else{ // TC asymmetry
@@ -2267,7 +2268,7 @@ void mlawNonLinearTVENonLinearTVP2::getIterated_DPhi(const double& T0, const dou
 		  mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm,Ge_TrEeTensor);
 	  }
 	  else{
-		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,true,alignedIndex,Bd_stiffnessTerm,Ge_TrEeTensor,rotationStiffness); // NEW EIGEN
+		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,true,alignedIndex,rotate,Bd_stiffnessTerm,Ge_TrEeTensor,rotationStiffness); // NEW EIGEN
 	  }
   }
 
@@ -2493,10 +2494,10 @@ void mlawNonLinearTVENonLinearTVP2::getIterated_DPhi(const double& T0, const dou
       // EIGEN NEW
       if(_rotationCorrectionScheme == 2){
     	  Ee0_rot = q0->_Ee;
-    	  STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex);
+    	  // STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate);
     	  // STensorOperation::alignEigenDecomposition_EigenVectorDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex);
     	  // STensorOperation::alignEigenDecomposition_TangentBased_indexOnly(q0->_Ee,Ee,alignedIndex);
-    	  // STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex);
+    	  STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate);
       }
 
       // STensorOperation::getEigenDecomposition(Ee,e1,e2,e3,E1,E2,E3);
@@ -2508,7 +2509,7 @@ void mlawNonLinearTVENonLinearTVP2::getIterated_DPhi(const double& T0, const dou
     		  mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm);
     	  }
     	  else{
-    		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,true,alignedIndex,Bd_stiffnessTerm,NULL,rotationStiffness); // NEW EIGEN
+    		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,true,alignedIndex,rotate,Bd_stiffnessTerm,NULL,rotationStiffness); // NEW EIGEN
     	  }
       }
       else{ // TC asymmetry
@@ -2516,7 +2517,7 @@ void mlawNonLinearTVENonLinearTVP2::getIterated_DPhi(const double& T0, const dou
     		  mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm,Ge_TrEeTensor);
     	  }
     	  else{
-    		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,true,alignedIndex,Bd_stiffnessTerm,Ge_TrEeTensor,rotationStiffness); // NEW EIGEN
+    		  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,true,alignedIndex,rotate,Bd_stiffnessTerm,Ge_TrEeTensor,rotationStiffness); // NEW EIGEN
     	  }
       }
 
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.h b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.h
index c867ffa2c..a7f40ccf6 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.h
@@ -31,7 +31,7 @@ class mlawNonLinearTVENonLinearTVP2 : public mlawNonLinearTVP{
                                             double& Ke, double& Ge, STensor43& Ge_Tensor,
                                             double& ptilde, STensor3& devPhi,
                                             STensor3& Phi, STensor3& N, STensor3& expGN, STensor43& dexpAdA,
-                                            STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv, fullVector<int>& alignedIndex,
+                                            STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv, fullVector<int>& alignedIndex, bool rotate = false,
                                             STensor43* Ge_TrEeTensor = NULL, STensor43* rotationStiffness = NULL) const;
 
         virtual void getDho(const double& Gamma, const STensor3& Cepr, const STensor3& Ceinvpr, const STensor3& KS,
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
index f2891929a..049251b86 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
@@ -2759,35 +2759,13 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
     const STensor3& DcorKirDT = q1->getRefToDcorKirDT();
     const STensor3& DDcorKirDT = q1->getRefToDDcorKirDTT(); // for mech source derivatives
     const STensor43& DDcorKirDTDEe = q1->getRefToDDcorKirDTDE();
-
-    // Add corKirinf terms in IP
-    /*
-    static STensor3 corKirDevInf, devEe;
-    double corKirTrInf, trEe;
-    STensorOperation::decomposeDevTr(E,devEe,trEe);
-
-    const STensor3& DEe = q1->getConstRefToDE();
-    const STensor3& devDE = q1->getConstRefToDevDE();
-    const double& trDE = q1->getConstRefToTrDE();
-
-    corKirInfinity(q1, devEe, trEe, T, corKirDevInf, corKirTrInf); // MISTAKE HERE -> corKirTrInf is the pressure  
-    
-    // Add extrabranch to corKir_inf
-    double tr_sigExtra(0.);
-    static STensor3 dev_sigExtra;
-    STensorOperation::decomposeDevTr(sigExtra,dev_sigExtra,tr_sigExtra);
-    corKirDevInf += dev_sigExtra;
-    corKirTrInf += 1/3*tr_sigExtra; // LHS is pressure
-    
-    // q1->_trCorKirinf_TVE = 3*corKirTrInf; //CHANGE
-    // q1->_devCorKirinf_TVE = corKirDevInf; //CHANGE 
-    */
-    const STensor3& DEe = q1->getConstRefToDE();
     
     // mechanicalSource
     mechanicalSource = 0.;
     double mechanicalSource_cap = 0.; // in case of mullins, undamaged mechanicalSource
 
+    const STensor3& DEe = q1->getConstRefToDE();
+
     // 1) ThermoElastic Heat
     if (this->getTimeStep() > 0){
         mechanicalSource += (STensorOperation::doubledot(DcorKirDT,DEe)*T/this->getTimeStep());
@@ -3134,7 +3112,7 @@ void mlawNonLinearTVM::get_Ee0s_Ai0s_ders(const STensor3& Ce, const STensor3& Ee
 
 void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ce, const STensor3& Ee, const STensor3& Ee0,
           	  							const IPNonLinearTVM *q0, IPNonLinearTVM *q1, double& Ke, double& Ge, const double T0, const double T1,
-										const bool stiff, const fullVector<int>& alignedIndex, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2, STensor43* rotationStiffness) const{
+										const bool stiff, const fullVector<int>& alignedIndex, const bool rotate, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2, STensor43* rotationStiffness) const{
 
 	// Caution!! -> For _rotationCorrectionScheme == 0 and 1,  Ee0 is already rotated. _R and its derivatives are in the IP
 
@@ -3179,7 +3157,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ce, co
     	// mlawNonLinearTVM::getEe0s2(Ce,Ee,devEe0,devEe0s,dDevEe0sdEe);
     	q1->_A_rot = q0->_A;
     	// mlawNonLinearTVM::getAi0s2(Ce,Ee,q0->_A,q1->_A_rot,dAi0sdEe);
-    	mlawNonLinearTVM::get_Ee0s_Ai0s_ders(Ce,Ee,devEe0,q0->_A,alignedIndex,devEe0s,dDevEe0sdEe,q1->_A_rot,dAi0sdEe);
+    	if(rotate){mlawNonLinearTVM::get_Ee0s_Ai0s_ders(Ce,Ee,devEe0,q0->_A,alignedIndex,devEe0s,dDevEe0sdEe,q1->_A_rot,dAi0sdEe);}
         // devEe0s.print("after");
     	// Msg::Info("alignedIndex = %d, %d, %d",alignedIndex(0),alignedIndex(1),alignedIndex(2)); // FLE
     }
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
index 7d73767c0..d32f1ffef 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
@@ -11,10 +11,6 @@
 #ifndef MLAWNONLINEARTVM_H_
 #define MLAWNONLINEARTVM_H_
 
-#include "mlaw.h"
-#include "STensor3.h"
-#include "STensor43.h"
-#include "STensor63.h"
 #include "ipNonLinearTVM.h"
 #include "viscosityLaw.h"
 #include "mlawHyperelasticFailureCriterion.h"
@@ -141,7 +137,7 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
           
     void ThermoViscoElasticPredictor_forTVP(const STensor3& Ce, const STensor3& Ee, const STensor3& Ee0,
           const IPNonLinearTVM *q0, IPNonLinearTVM *q1,  double& Ke, double& Ge, const double T0, const double T1,
-		  const bool stiff,  const fullVector<int>& alignedIndex, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2 = NULL, STensor43* rotationStiffness = NULL) const;
+		  const bool stiff,  const fullVector<int>& alignedIndex, const bool rotate, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2 = NULL, STensor43* rotationStiffness = NULL) const;
               
     void isotropicHookTensor(const double K, const double G, STensor43& L) const;
                                   
-- 
GitLab