diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index dd64c2a427de4710f32cdc0752d877c3a5d605e5..6b0ed3fbddcbcc351c12a6fa8d6bac4925d07bf7 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -4333,6 +4333,9 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
const StochDMNDG3DIPVariable* ipv_prev = static_cast<const StochDMNDG3DIPVariable*>(ipvprev);
const STensor3& F_unrot = ipvcur->getConstRefToDeformationGradient();
+ 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;
if (_flag_rotation){
@@ -4346,8 +4349,10 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
}
static STensor3 RT; STensorOperation::zero(RT);
STensorOperation::transposeSTensor3(R,RT);
-
- STensorOperation::multSTensor3InPlace2nd(RT,F);
+
+ STensorOperation::multSTensor3InPlace2nd(R,U);
+ STensorOperation::multSTensor3InPlace1st(U,RT);
+ F = U; // R_def in the rotated frame = I;
}
STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress();
@@ -4490,39 +4495,68 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
}
if (_flag_rotation){
- const STensor3& R = ipvcur->getConstReftoRotationMatrix();
- static STensor3 RT; STensorOperation::zero(RT); STensorOperation::transposeSTensor3(R,RT);
- static STensor3 I2; STensorOperation::zero(I2); I2(0,0) = I2(1,1) = I2(2,2) = 1.;
- 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 STensor43 dFdF_rot; STensorOperation::zero(dFdF_rot);
+ const STensor3& R = ipvcur->getConstReftoRotationMatrix();
+ static STensor3 RT; STensorOperation::zero(RT); STensorOperation::transposeSTensor3(R, RT);
+ static STensor3 Finvrot, Prot;
+
+ Prot = P;
+ STensorOperation::zero(P);
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++){
- dFdF_rot(m,n,k,l) = 0.;
- for(int p=0; p<3; p++)
- dFdF_rot(m,n,k,l) += RT(m,p)*I4_first_type(p,n,k,l);
- }
-
- static STensor43 L_rot;
- STensorOperation::multSTensor43InPlace(L,dFdF_rot);
- for(int i=0; i<3; i++)
- for(int k=0; k<3; k++)
- for(int o=0; o<3; o++)
- for(int p=0; p<3; p++){
- L_rot(i,k,o,p) = 0.;
- for(int j=0; j<3; j++)
- L_rot(i,k,o,p) += R(i,j)*L(j,k,o,p);
+ for(int p=0; p<3; p++)
+ for(int q=0; q<3; q++)
+ for(int r=0; r<3; r++)
+ P(m,n) += R_def(m,p)*RT(p,q)*Prot(q,r)*R(r,n);
+
+ if(stiff){
+ static STensor43 dFinvrot_dFrot; STensorOperation::zero(dFinvrot_dFrot);
+ STensorOperation::inverseSTensor3(F,Finvrot,&dFinvrot_dFrot,false);
+
+ static STensor3 I2; STensorOperation::zero(I2); I2(0,0) = I2(1,1) = I2(2,2) = 1.;
+ 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);
+
+ static STensor43 dFrot_dF, dFinvrot_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++)
+ for(int q=0; q<3; q++)
+ for(int r=0; r<3; r++)
+ dFrot_dF(m,n,k,l) += R(m,p)*R_def(q,p)*I4_first_type(q,r,k,l)*RT(r,n); // R_def_T is used here
}
-
- L = 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 L_rot; STensorOperation::zero(L_rot);
+ 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++)
+ for(int q=0; q<3; q++)
+ for(int r=0; r<3; r++)
+ L_rot(m,n,k,l) += R_def(m,p)*RT(p,q)*L(q,r,k,l)*R(r,n);
+
+ // Reassign
+ L = L_rot;
+ }
}
}
@@ -5753,10 +5787,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
STensorOperation::multSTensor3InPlace2nd(R,U);
STensorOperation::multSTensor3InPlace1st(U,RT);
- // STensorOperation::multSTensor3(R_def,U,F);
F = U; // R_def in the rotated frame = I;
- // STensorOperation::multSTensor3InPlace1st(F,RT);
STensorOperation::multSTensor3SVector3InPlace2nd(R,H);
}
@@ -6154,27 +6186,36 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
const STensor3& R = ipvcur->getConstReftoRotationMatrix();
static STensor3 RT; STensorOperation::zero(RT); STensorOperation::transposeSTensor3(R, RT);
- static STensor3 Finvrot, S, Srot, Prot, dSdT;
+ static STensor3 Finvrot, Prot;
Prot = P;
- STensorOperation::inverseSTensor3(F,Finvrot); // Note that F is Frot here
+ STensorOperation::zero(P);
+ for(int m=0; m<3; m++)
+ for(int n=0; n<3; n++)
+ for(int p=0; p<3; p++)
+ for(int q=0; q<3; q++)
+ for(int r=0; r<3; r++)
+ P(m,n) += R_def(m,p)*RT(p,q)*Prot(q,r)*R(r,n);
- STensorOperation::multSTensor3(Finvrot,Prot,Srot);
- STensorOperation::multSTensor3InPlace2nd(RT,Srot);
- STensorOperation::multSTensor3(Srot,R,S);
+ STensorOperation::multSTensor3SVector3InPlace2nd(RT,Q);
if(stiff){
static STensor43 dFinvrot_dFrot; STensorOperation::zero(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);
-
+
+ static STensor3 dProt_dT;
+ dProt_dT = dPdT;
+ STensorOperation::zero(dPdT);
+ for(int m=0; m<3; m++)
+ for(int n=0; n<3; n++)
+ for(int p=0; p<3; p++)
+ for(int q=0; q<3; q++)
+ for(int r=0; r<3; r++)
+ dPdT(m,n) += R_def(m,p)*RT(p,q)*dProt_dT(q,r)*R(r,n);
+
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++)
@@ -6183,40 +6224,18 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
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) += R_def(p,m)*I4_first_type(p,n,k,l); // R_def_T is used here
- 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 STensor43 dFrot_dF, dFinvrot_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++){
dHrot_dH(m,n) += R(m,k)*I2(k,n);
- for(int l=0; l<3; l++)
- for(int r=0; r<3; r++)
- for(int s=0; s<3; s++)
- dUrot_dF(m,n,k,l) += dUrot_dU(m,n,r,s)*dUdF(r,s,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++)
- dFrot_dF(m,n,k,l) += R_def(m,p)*dUrot_dF(p,n,k,l);*/
- dFrot_dF = dUrot_dF; // R_def_rot = I;
+ for(int q=0; q<3; q++)
+ for(int r=0; r<3; r++)
+ dFrot_dF(m,n,k,l) += R(m,p)*R_def(q,p)*I4_first_type(q,r,k,l)*RT(r,n); // R_def_T is used here
+ }
for(int m=0; m<3; m++)
for(int n=0; n<3; n++)
@@ -6228,39 +6247,19 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
// dPdF
- static STensor43 dSdF, dSrot_dF, L_rot;
- static STensor33 dqdF_rot;
- static STensor3 dwdf_rot, dmechSourcedf_rot;
-
- STensorOperation::zero(L_rot);
- STensorOperation::zero(dqdF_rot);
- STensorOperation::zero(dSdF);
- STensorOperation::zero(dSrot_dF);
-
+ static STensor43 L_rot; STensorOperation::zero(L_rot);
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);
+ for(int r=0; r<3; r++)
+ L_rot(m,n,k,l) += R_def(m,p)*RT(p,q)*L(q,r,k,l)*R(r,n);
// dqdF
+ static STensor33 dqdF_rot; STensorOperation::zero(dqdF_rot);
for(int m=0; m<3; m++)
for(int n=0; n<3; n++)
for(int k=0; k<3; k++)
@@ -6274,18 +6273,13 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
STensorOperation::multSTensor3InPlace2nd(RT,dqdgradT);
// thermSrc, mechSrc
+ static STensor3 dwdf_rot, dmechSourcedf_rot;
STensorOperation::multSTensor3STensor43(dwdf,dFrot_dF,dwdf_rot);
STensorOperation::multSTensor3STensor43(dmechSourcedf,dFrot_dF,dmechSourcedf_rot);
// Reassign
L = L_rot; dqdF = dqdF_rot; dwdf = dwdf_rot; dmechSourcedf = dmechSourcedf_rot;
}
-
- STensorOperation::multSTensor3(F_unrot,S,P);
-
- // STensorOperation::multSTensor3(Prot,R,P);
- STensorOperation::multSTensor3SVector3InPlace2nd(RT,Q);
-
}
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}