[NEX PATCH] Rotation implementation in StochTMDMN.. has been updated. Feature yet to be debugged.

dG3D/src/dG3DMaterialLaw.cpp

@@ -5733,6 +5733,9 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
   const SVector3& H_unrot = ipvcur->getConstRefToGradT();
   const double T    = ipvcur->getConstRefToTemperature();
+  static STensor3 U, R_def; // R_def is the Rotation in F = R_def . U
+  STensorOperation::RUDecomposition(F_unrot,U,R_def);
+  
   static STensor3 F; STensorOperation::zero(F); F = F_unrot;
   static SVector3 H; STensorOperation::zero(H); H = H_unrot;
   if (_flag_rotation){
@@ -5748,8 +5751,10 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
     static STensor3 RT; STensorOperation::zero(RT);
     STensorOperation::transposeSTensor3(R,RT);
-    STensorOperation::multSTensor3InPlace2nd(RT,F);
-    STensorOperation::multSTensor3SVector3InPlace2nd(RT,H);
+    STensorOperation::multSTensor3InPlace2nd(R,U);
+    STensorOperation::multSTensor3InPlace1st(U,RT);
+    F = U; // R_def in the rotated frame = I;
+    STensorOperation::multSTensor3SVector3InPlace2nd(R,H);
   }
   STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress(); STensorOperation::zero(P);
@@ -6055,12 +6060,6 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
       }  
     }
-    if (_flag_rotation){
-      const STensor3& R = ipvcur->getConstReftoRotationMatrix();
-      STensorOperation::multSTensor3InPlace2nd(R,P);
-      STensorOperation::multSTensor3SVector3InPlace2nd(R,Q);
-    }
-
     if(stiff){
       
       NTVC.mult(_I, dRdF); dRdT = NTV_dPQdT;
@@ -6146,66 +6145,142 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
             }
         }
       }
+    }
       
     if (_flag_rotation){
+
       const STensor3& R = ipvcur->getConstReftoRotationMatrix();
       static STensor3 RT; STensorOperation::zero(RT); STensorOperation::transposeSTensor3(R, RT);
-        STensorOperation::multSTensor3InPlace2nd(R,dPdT);
-        STensorOperation::multSTensor3SVector3InPlace2nd(R,dqdT);
+      static STensor3 Finvrot, S, Srot, Prot, dSdT;
+
+      Prot = P;
+      STensorOperation::inverseSTensor3(F,Finvrot); // Note that F is Frot here
+      STensorOperation::multSTensor3(Finvrot,Prot,Srot);
+      STensorOperation::multSTensor3InPlace2nd(RT,Srot);
+      STensorOperation::multSTensor3(Srot,R,S);
+
+      if(stiff){
+        static STensor43 dFinvrot_dFrot;
+        STensorOperation::inverseSTensor3(F,Finvrot,&dFinvrot_dFrot,false);
+        STensorOperation::multSTensor3(Finvrot,dPdT,dSdT); // Since dPdT is dProt_dT initially, dSdT is dSrot_dT is initially
+        STensorOperation::multSTensor3InPlace2nd(RT,dSdT);
+        STensorOperation::multSTensor3InPlace1st(dSdT,R); // Here dSrot_dT becomes dSdT
+        STensorOperation::multSTensor3(F_unrot,dSdT,dPdT);
+  
+        STensorOperation::multSTensor3SVector3InPlace2nd(RT,dqdT);
         static STensor3 I2; STensorOperation::zero(I2); I2(0,0) = I2(1,1) = I2(2,2) = 1.;
+        const STensor43 I4_symm(1.,1.); // symmetric type of fourth order unit tensor
         static STensor43 I4_first_type; STensorOperation::zero(I4_first_type); // first type of fourth order unit tensor
         for(int m=0; m<3; m++)
           for(int n=0; n<3; n++)
             for(int k=0; k<3; k++)
               for(int l=0; l<3; l++)
-                I4_first_type(m,n,k,l) = I2(m,k)*I2(n,l); // = dFrotdFrot
-        static STensor3 dHdH_rot; STensorOperation::zero(dHdH_rot);
-        static STensor43 dFdF_rot; STensorOperation::zero(dFdF_rot);
+                I4_first_type(m,n,k,l) = I2(m,k)*I2(n,l);
+
+        static STensor3 dHrot_dH; STensorOperation::zero(dHrot_dH);
+        static STensor43 dUdF, dUrot_dU, dUrot_dF, dFrot_dF, dFinvrot_dF; 
+        STensorOperation::zero(dUdF);
+        STensorOperation::zero(dUrot_dU);
+        STensorOperation::zero(dUrot_dF);
+        STensorOperation::zero(dFrot_dF);
+        STensorOperation::zero(dFinvrot_dF);
+
+        for(int m=0; m<3; m++)
+          for(int n=0; n<3; n++)
+            for(int k=0; k<3; k++)
+              for(int l=0; l<3; l++)
+                for(int p=0; p<3; p++){
+                  dUdF(m,n,k,l) += I4_first_type(m,p,k,l)*R_def(n,p);
+                  for(int q=0; q<3; q++)
+                    dUrot_dU(m,n,k,l) += R(m,p)*I4_symm(p,q,k,l)*RT(q,n);
+                }
         for(int m=0; m<3; m++)
           for(int n=0; n<3; n++)
             for(int k=0; k<3; k++){
-              dHdH_rot(m,n) += RT(m,k)*I2(k,n);
-              for(int l=0; l<3; l++){
-                dFdF_rot(m,n,k,l) = 0.;
+              dHrot_dH(m,n) += R(m,k)*I2(k,n);
+              for(int l=0; l<3; l++)
                 for(int p=0; p<3; p++)
-                  dFdF_rot(m,n,k,l) += RT(m,p)*I4_first_type(p,n,k,l);
-              }
+                  for(int q=0; q<3; q++)
+                    for(int r=0; r<3; r++)
+                      for(int s=0; s<3; s++)
+                        dUrot_dF(m,n,k,l) += I4_symm(m,n,p,q)*dUrot_dU(p,q,r,s)*dUdF(r,s,k,l);
             }
                         
-        static STensor43 L_rot;
+        dFrot_dF = dUrot_dF; // R_def_rot = I;
+        
+        for(int m=0; m<3; m++)
+          for(int n=0; n<3; n++)
+            for(int k=0; k<3; k++)
+              for(int l=0; l<3; l++)
+                for(int p=0; p<3; p++)
+                  for(int q=0; q<3; q++)
+                    dFinvrot_dF(m,n,k,l) += dFinvrot_dFrot(m,n,p,q)*dFrot_dF(p,q,k,l);
+                
+        
+        // dPdF
+        static STensor43 dSdF, dSrot_dF, L_rot;
         static STensor33 dqdF_rot;
         static STensor3 dqdgradT_rot, dwdf_rot, dmechSourcedf_rot;
         
-        STensorOperation::multSTensor43InPlace(L,dFdF_rot);
-        STensorOperation::multSTensor33STensor43InPlace1st(dqdF,dFdF_rot);
-        STensorOperation::multSTensor3InPlace1st(dqdgradT,dHdH_rot);
-        STensorOperation::multSTensor3STensor43(dwdf,dFdF_rot,dwdf_rot);
-        STensorOperation::multSTensor3STensor43(dmechSourcedf,dFdF_rot,dmechSourcedf_rot);
-        for(int i=0; i<3; i++)
-          for(int k=0; k<3; k++){
-            dqdgradT_rot(i,k) = 0.;
-            for(int o=0; o<3; o++){
-              dqdF_rot(i,k,o) = 0.;
-              dqdgradT_rot(i,k) += R(i,o)*dqdgradT(o,k);
-              for(int p=0; p<3; p++){
-                L_rot(i,k,o,p) = 0.;
-                dqdF_rot(i,k,o) += R(i,p)*dqdF(p,k,o);
-                for(int j=0; j<3; j++)
-                  L_rot(i,k,o,p) += R(i,j)*L(j,k,o,p);
-              }
-            }
-          }
+        STensorOperation::zero(L_rot);
+        STensorOperation::zero(dqdF_rot);
+        STensorOperation::zero(dSdF);
+        STensorOperation::zero(dSrot_dF);
+        
+        STensorOperation::multSTensor43InPlace(L,dFrot_dF); // L = dProt_dFrot and becomes dProt_dF
+        for(int m=0; m<3; m++)
+          for(int n=0; n<3; n++)
+            for(int k=0; k<3; k++)
+              for(int l=0; l<3; l++)
+                for(int p=0; p<3; p++)
+                  dSrot_dF(m,n,k,l) += dFinvrot_dF(m,p,k,l)*Prot(p,n) + Finvrot(m,p)*L(p,n,k,l);
+        
+        for(int m=0; m<3; m++)
+          for(int n=0; n<3; n++)
+            for(int k=0; k<3; k++)
+              for(int l=0; l<3; l++)
+                for(int p=0; p<3; p++)
+                  for(int q=0; q<3; q++)
+                    dSdF(m,n,k,l) += RT(m,p)*dSrot_dF(p,q,k,l)*R(q,n);
+        
+        for(int m=0; m<3; m++)
+          for(int n=0; n<3; n++)
+            for(int k=0; k<3; k++)
+              for(int l=0; l<3; l++)
+                for(int p=0; p<3; p++)
+                  L_rot(m,n,k,l) += I4_first_type(m,p,k,l)*S(p,n) + F_unrot(m,p)*dSdF(p,n,k,l);
         
+        // dqdF  
+        for(int m=0; m<3; m++)
+          for(int n=0; n<3; n++)
+            for(int k=0; k<3; k++)
+              for(int l=0; l<3; l++)
+                for(int p=0; p<3; p++)
+                  for(int q=0; q<3; q++)
+                    dqdF_rot(m,n,k) += RT(m,l)*dqdF(l,p,q)*dFrot_dF(p,q,n,k);
+
+        // dqdH           
+        STensorOperation::multSTensor3InPlace1st(dqdgradT,dHrot_dH);
+        STensorOperation::multSTensor3InPlace2nd(RT,dqdgradT);
+
+        // thermSrc, mechSrc
+        STensorOperation::multSTensor3STensor43(dwdf,dFrot_dF,dwdf_rot);
+        STensorOperation::multSTensor3STensor43(dmechSourcedf,dFrot_dF,dmechSourcedf_rot);
+          
+        // Reassign
         L = L_rot; dqdF = dqdF_rot; dqdgradT = dqdgradT_rot; dwdf = dwdf_rot; dmechSourcedf = dmechSourcedf_rot; 
       }
     
-    }
+      STensorOperation::multSTensor3(F_unrot,S,P);
+      STensorOperation::multSTensor3SVector3InPlace2nd(RT,Q);
+
     }
     ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
   }
+} 
 void StochTMDMNDG3DMaterialLaw::setTime(const double t,const double dtime){
   dG3DMaterialLaw::setTime(t,dtime);