diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index c5123a1321e4e2f6877c2f86cce1dcfa3c326bfc..c60fb150a339d9c4a4db8353245d24975fa1c348 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/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,65 +6145,141 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
}
}
}
+ }
+
+ if (_flag_rotation){
- if (_flag_rotation){
- const STensor3& R = ipvcur->getConstReftoRotationMatrix();
- static STensor3 RT; STensorOperation::zero(RT); STensorOperation::transposeSTensor3(R,RT);
+ const STensor3& R = ipvcur->getConstReftoRotationMatrix();
+ static STensor3 RT; STensorOperation::zero(RT); STensorOperation::transposeSTensor3(R, RT);
+ static STensor3 Finvrot, S, Srot, Prot, dSdT;
+
+ Prot = P;
+ STensorOperation::inverseSTensor3(F,Finvrot); // Note that F is Frot here
- STensorOperation::multSTensor3InPlace2nd(R,dPdT);
- STensorOperation::multSTensor3SVector3InPlace2nd(R,dqdT);
+ 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
+ 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);
+ }
- static STensor3 dHdH_rot; STensorOperation::zero(dHdH_rot);
- static STensor43 dFdF_rot; STensorOperation::zero(dFdF_rot);
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);
}
+
+ dFrot_dF = dUrot_dF; // R_def_rot = I;
- static STensor43 L_rot;
+ 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);
}
- ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
void StochTMDMNDG3DMaterialLaw::setTime(const double t,const double dtime){