Skip to content
Snippets Groups Projects
Commit 74a11ab2 authored by Ujwal Kishore Jinaga's avatar Ujwal Kishore Jinaga :clown:
Browse files

[MAJOR PATCH] Rotation tangents in TMDMN and DMN are now fixed. For future:... 

[MAJOR PATCH] Rotation tangents in TMDMN and DMN are now fixed. For future: dont make derivatives of F,P in terms of U,S respectively?
parent 61c52988
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
dG3D/src/dG3DMaterialLaw.cpp
+ 95
− 101
View file @ 74a11ab2
...@@ -4333,6 +4333,9 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c ...@@ -4333,6 +4333,9 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
const StochDMNDG3DIPVariable* ipv_prev = static_cast<const StochDMNDG3DIPVariable*>(ipvprev); const StochDMNDG3DIPVariable* ipv_prev = static_cast<const StochDMNDG3DIPVariable*>(ipvprev);
const STensor3& F_unrot = ipvcur->getConstRefToDeformationGradient(); 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; static STensor3 F; STensorOperation::zero(F); F = F_unrot;
if (_flag_rotation){ if (_flag_rotation){
   
...@@ -4347,7 +4350,9 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c ...@@ -4347,7 +4350,9 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
static STensor3 RT; STensorOperation::zero(RT); static STensor3 RT; STensorOperation::zero(RT);
STensorOperation::transposeSTensor3(R,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(); STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress();
...@@ -4490,8 +4495,23 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c ...@@ -4490,8 +4495,23 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
} }
   
if (_flag_rotation){ if (_flag_rotation){
const STensor3& R = ipvcur->getConstReftoRotationMatrix(); const STensor3& R = ipvcur->getConstReftoRotationMatrix();
static STensor3 RT; STensorOperation::zero(RT); STensorOperation::transposeSTensor3(R, RT); 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 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 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 static STensor43 I4_first_type; STensorOperation::zero(I4_first_type); // first type of fourth order unit tensor
...@@ -4499,31 +4519,45 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c ...@@ -4499,31 +4519,45 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
for(int n=0; n<3; n++) for(int n=0; n<3; n++)
for(int k=0; k<3; k++) for(int k=0; k<3; k++)
for(int l=0; l<3; l++) 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 STensor43 dFrot_dF, dFinvrot_dF; STensorOperation::zero(dFrot_dF); STensorOperation::zero(dFinvrot_dF);
static STensor43 dFdF_rot; STensorOperation::zero(dFdF_rot);
for(int m=0; m<3; m++) for(int m=0; m<3; m++)
for(int n=0; n<3; n++) for(int n=0; n<3; n++)
for(int k=0; k<3; k++) for(int k=0; k<3; k++){
for(int l=0; l<3; l++){ for(int l=0; l<3; l++)
dFdF_rot(m,n,k,l) = 0.;
for(int p=0; p<3; p++) 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++)
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
} }
static STensor43 L_rot; for(int m=0; m<3; m++)
STensorOperation::multSTensor43InPlace(L,dFdF_rot); for(int n=0; n<3; n++)
for(int i=0; i<3; i++)
for(int k=0; k<3; k++) for(int k=0; k<3; k++)
for(int o=0; o<3; o++) for(int l=0; l<3; l++)
for(int p=0; p<3; p++){ for(int p=0; p<3; p++)
L_rot(i,k,o,p) = 0.; for(int q=0; q<3; q++)
for(int j=0; j<3; j++) dFinvrot_dF(m,n,k,l) += dFinvrot_dFrot(m,n,p,q)*dFrot_dF(p,q,k,l);
L_rot(i,k,o,p) += R(i,j)*L(j,k,o,p);
}
// 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; L = L_rot;
} }
}
   
} }
} }
...@@ -5753,10 +5787,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5753,10 +5787,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
   
STensorOperation::multSTensor3InPlace2nd(R,U); STensorOperation::multSTensor3InPlace2nd(R,U);
STensorOperation::multSTensor3InPlace1st(U,RT); STensorOperation::multSTensor3InPlace1st(U,RT);
// STensorOperation::multSTensor3(R_def,U,F);
F = U; // R_def in the rotated frame = I; F = U; // R_def in the rotated frame = I;
   
// STensorOperation::multSTensor3InPlace1st(F,RT);
STensorOperation::multSTensor3SVector3InPlace2nd(R,H); STensorOperation::multSTensor3SVector3InPlace2nd(R,H);
} }
   
...@@ -6154,27 +6186,36 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -6154,27 +6186,36 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
   
const STensor3& R = ipvcur->getConstReftoRotationMatrix(); const STensor3& R = ipvcur->getConstReftoRotationMatrix();
static STensor3 RT; STensorOperation::zero(RT); STensorOperation::transposeSTensor3(R, RT); static STensor3 RT; STensorOperation::zero(RT); STensorOperation::transposeSTensor3(R, RT);
static STensor3 Finvrot, S, Srot, Prot, dSdT; static STensor3 Finvrot, Prot;
   
Prot = P; 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::multSTensor3SVector3InPlace2nd(RT,Q);
STensorOperation::multSTensor3InPlace2nd(RT,Srot);
STensorOperation::multSTensor3(Srot,R,S);
   
if(stiff){ if(stiff){
static STensor43 dFinvrot_dFrot; STensorOperation::zero(dFinvrot_dFrot); static STensor43 dFinvrot_dFrot; STensorOperation::zero(dFinvrot_dFrot);
STensorOperation::inverseSTensor3(F,Finvrot,&dFinvrot_dFrot,false); 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); static STensor3 dProt_dT;
STensorOperation::multSTensor3InPlace1st(dSdT,R); // Here dSrot_dT becomes dSdT dProt_dT = dPdT;
STensorOperation::multSTensor3(F_unrot,dSdT,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); STensorOperation::multSTensor3SVector3InPlace2nd(RT,dqdT);
   
static STensor3 I2; STensorOperation::zero(I2); I2(0,0) = I2(1,1) = I2(2,2) = 1.; 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 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 m=0; m<3; m++)
for(int n=0; n<3; n++) for(int n=0; n<3; n++)
...@@ -6183,40 +6224,18 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -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); I4_first_type(m,n,k,l) = I2(m,k)*I2(n,l);
   
static STensor3 dHrot_dH; STensorOperation::zero(dHrot_dH); static STensor3 dHrot_dH; STensorOperation::zero(dHrot_dH);
static STensor43 dUdF, dUrot_dU, dUrot_dF, dFrot_dF, dFinvrot_dF; static STensor43 dFrot_dF, dFinvrot_dF; STensorOperation::zero(dFrot_dF); STensorOperation::zero(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);
}
   
for(int m=0; m<3; m++) for(int m=0; m<3; m++)
for(int n=0; n<3; n++) for(int n=0; n<3; n++)
for(int k=0; k<3; k++){ for(int k=0; k<3; k++){
dHrot_dH(m,n) += R(m,k)*I2(k,n); dHrot_dH(m,n) += R(m,k)*I2(k,n);
for(int l=0; l<3; l++) 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++) for(int r=0; r<3; r++)
for(int s=0; s<3; s++) 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
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 m=0; m<3; m++) for(int m=0; m<3; m++)
for(int n=0; n<3; n++) for(int n=0; n<3; n++)
...@@ -6228,39 +6247,19 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -6228,39 +6247,19 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
// dPdF // dPdF
static STensor43 dSdF, dSrot_dF, L_rot; static STensor43 L_rot; STensorOperation::zero(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);
STensorOperation::multSTensor43InPlace(L,dFrot_dF); // L = dProt_dFrot and becomes dProt_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 m=0; m<3; m++)
for(int n=0; n<3; n++) for(int n=0; n<3; n++)
for(int k=0; k<3; k++) for(int k=0; k<3; k++)
for(int l=0; l<3; l++) for(int l=0; l<3; l++)
for(int p=0; p<3; p++) for(int p=0; p<3; p++)
for(int q=0; q<3; q++) 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 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);
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 // dqdF
static STensor33 dqdF_rot; STensorOperation::zero(dqdF_rot);
for(int m=0; m<3; m++) for(int m=0; m<3; m++)
for(int n=0; n<3; n++) for(int n=0; n<3; n++)
for(int k=0; k<3; k++) for(int k=0; k<3; k++)
...@@ -6274,18 +6273,13 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -6274,18 +6273,13 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
STensorOperation::multSTensor3InPlace2nd(RT,dqdgradT); STensorOperation::multSTensor3InPlace2nd(RT,dqdgradT);
   
// thermSrc, mechSrc // thermSrc, mechSrc
static STensor3 dwdf_rot, dmechSourcedf_rot;
STensorOperation::multSTensor3STensor43(dwdf,dFrot_dF,dwdf_rot); STensorOperation::multSTensor3STensor43(dwdf,dFrot_dF,dwdf_rot);
STensorOperation::multSTensor3STensor43(dmechSourcedf,dFrot_dF,dmechSourcedf_rot); STensorOperation::multSTensor3STensor43(dmechSourcedf,dFrot_dF,dmechSourcedf_rot);
// Reassign // Reassign
L = L_rot; dqdF = dqdF_rot; dwdf = dwdf_rot; dmechSourcedf = dmechSourcedf_rot; 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); ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment