From 1cd6b6f3bc5dbfb05aa0f83fad3b0f47142e5c1b Mon Sep 17 00:00:00 2001
From: FLE_Knight <ujwalkishore.jinaga@mail.polimi.it>
Date: Tue, 7 Feb 2023 15:38:26 +0100
Subject: [PATCH] restoredtoFixedDCeDF_doesntWork_why?
---
.../materialLaw/.mlawNonLinearTVP.cpp.swp | Bin 16384 -> 0 bytes
.../materialLaw/mlawNonLinearTVP.cpp | 1228 ++++++++++++++---
.../materialLaw/mlawNonLinearTVP.h | 18 +-
dG3D/src/dG3DMaterialLaw.cpp | 3 +-
dG3D/src/dG3DMaterialLaw.h | 4 +-
5 files changed, 1039 insertions(+), 214 deletions(-)
delete mode 100644 NonLinearSolver/materialLaw/.mlawNonLinearTVP.cpp.swp
diff --git a/NonLinearSolver/materialLaw/.mlawNonLinearTVP.cpp.swp b/NonLinearSolver/materialLaw/.mlawNonLinearTVP.cpp.swp
deleted file mode 100644
index 0a449d84b681569ab414d2a0a02038080d9f5ca3..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001
literal 16384
zcmYc?2=nw+u+TGPU|?VnU|@JD`Y^eDY7>K#6&FKkR(WC$NSpw6^Ks2A&PXjvC7`Yj
zY@mK}uCY&MQc+@2W@@p%Uw)oXW?pJyQE+}vS!$7fZemGlQD$O}Phz=#Zcbu3Ts|Z$
zKrgwVfPfXF(xV|T8Ui>&fS1A8$j|^Jt*oS|AS@ILVvgd`5Eu=C(GVC7fzc2c4S~@R
z7!85Z5Eu=C5fTC=1u_iv3=9lRQ2$0jX+|{qKNKWG#iOA#Ox-`IygyWYAC!j4|AoqX
zLB)BY5>Sfa4^%!9DjyD&fKm*<q4I%H`3+DBD8=v#DxU(CXM}nPCjS#Ep9qytMU(#l
zl{bgVA48M>4wbir%Kt@^{|1$}gUXvRLl`joze43<`b*H{zd+?(A=(($qRD@T%EQ9_
z6Po-ds5~rwd|4m{!R-GCl`nz_Gi*SU{{WTGfXXYeLKrap@1gSf5MhQcH2HUE@;}k!
z-$LcHp!ze|APkuKZ=mw9_`8HA{~9U}8)GtKN7erdDi1S%Hk$lPsC+F%8v`o`gaNbv
z1ynv4BFxZ)CjT5NUksIhk0$>NDi1S1i4&p_W<E5!7=$3g4ENFGp~9ooXb6mkz-S1J
zhQMeDjE2By2#kinXb6mkz-S22ECdo$7#O@67#PljIw}GT46y$HUw#ILd;AOxxA_?u
zBKR2?LirgO6!{q#<oFpF?(#7(+~Q+kxX#DGaE*_FVLKlK!!$kyh6+9ghH^dzhB7_|
z21Y&xhF)F<22NfEh7=wK200!E23Z~khIiZy3_H0Q7}~fQ7=pMN7-YB^7=*YP82Gsv
z82Gpu7(Q?@FihuSV3@?kz!1;Hz~IBhz~Ifrz`(@C!0?!pfnh5r1H%?h28KdT28Lix
z1_n(|1_m`w28O#F3=I1?7#JpTFfjCTFfjOVFfeFyFfjaPXJGij&cN`Moq^#CI|IWi
zb_Rx}><kPw><kRP><kQl*ccd|voSENXJcTP&c?u?&&I&O$;QBNiIst&h?RlCmz9A*
zij{$3D+>cdKMMmx6bl2xZ)OID%ghW6vzQqeCNMKFbTKn9R4_9zlrl3gfco{m%nS@H
z%nS?%nHU(VnHU)KnHU&!m>3wenHU&uFfuS)XJlZwzz7NRbD(giV(jYcEBLvFGcf4u
zD>!?&I(sX4x+!=D1cx|=xGMO$y1E1_lw_nT6es4UDx{>QW#(m;Waj57Bo-@V<|%;I
zTLk2nrxrzKrskx0R2HNbDd_6zDujDDhS=LPFw|NzC@3f>WTq)-#FwPzx|QZ7$CqUm
zCo3r0D)@!^_-HDmq~@fSq$;4uTBC9yO05+XFwOEv!fldA5|TL{N%%~`V+gVtc#H_H
zG{9|0aHWB@0>V!qwNRrNKmmxnV8zO+4qZq?vl<?dMX4pFMR^KZB^jB;kf5yv2MQ>V
zA;AM?DInPf4sKh8;-X|-JCsn?wM)*)&r8+Ngct<U2#P06g>a9-t;$RzEjsn}6-rWb
z3sQ>`OG=AC$*3STFC{fEIaMJ)O#$S<{Nl`#O1SYLr@)=+k)(j`dXFTc+>Vqu;6@|6
z97!G{8G;fIURP&l=B4H)mSiR?WF!`)q~>MjrNb>jbw8%_iE};vVE2GaW4hg10k6~H
z4pt~GDM~Fa#_n!NDn(5V!IcI?xjnej0G7ewh9ieO$ZVJb4CjM$3tkU^%b>&@g-TGN
zgYJoVk05x6djuiUc6>-;WlnxkU}<JvYNcCZa!Gy>TpgwqL|<l_MsX3Mu&M@oCOJQ^
zxCC7CA!l4IP!3bDO{vUF%*{-WPfjc@vB9CzPD8CYzqBYh6`J$(^}(j0m7`Wx`QQwb
zUu3I+q8zP!aK)(*qi|6J+k@<CWW7+=B39F>DM0Ool|Nwb!!sr*-$FuA4@C?mg+NOb
zR2A@~i`76-y+9xZgRDjr8K|bf$_s>*XvGDFlaR^*EJhLG9PDKQnpx;U<$*t(5c!({
zsWA6YfL0tN=5SQEfU+OBOhtqddd|iOB-AVok8p@lM1&F6JPUCaEW_%>!}2T6ysC#q
z29i}V61$ZZT1E-bfTnP?HKJ;uoU3Vto;O@USzlA38q_Y*0ZXD~P_R6_W(Jk{uo@U1
zI-qo_V5^XqTCM=LJ29^$I5{yVu?W-<0yU#F4E3Oe0CFUQYCKZ)kZuhm>WH&PUtb}q
z64d-lEX^s=QOL{3QJ*t_V-y?~xCtzAs9~UIsi~t7?-2wF5lA_U+W@c(TI{-N<dx>?
zD7fk<<dy0u6lLV=D3s*q=qQ9h6P}KOV@^RvqK<-hNQQHMo{oZZfsR6MVo6bEMN(xz
zYEemPQBoo(N^}$oQj1D-K#ks_+{7H$;*v~IJ)fT!P?TDhnO|D$oS&DLnXairj90<2
zu8@?UpM%3Lh=P>-(xjYJsKpQ=u%42PRJ?AmBHj&9D-gi}u?8v(4I30OSU{l&Iv3~=
zpVKn)N)$k`f?|TJ4hk0(TPOnHs6!C}#Uld)gK}nGa!zSVs)7<))CcHgC?Uz_!PbZS
zrY2`V)`+7j2o6cjE6y+S2Uo9|`FX`4X$F0L2vTs)FQ_caOwTB>QV32>RdDtXi1hSx
zS4hjuNmVGxNGwrE&d*IPRw&Ob$xr~b!xd68i%W_!lS)CU2Wp06X-P(Yk(Gi?C}=^v
zf_G+dMt)JMf|r7>f}4+PymubRd^;Tl4QEXSBLgENs5Y0>;^d;t0+2yg3V!){x{!tL
z3L&6$kgpq-S)81&3-VTCUS@K!f(Ej<YffTu2~<r0n62QGSgr|iC)lj~<is2Wm&Dw}
z^i&1Uypq(Sw8Z381&!kT{5(An(1g0eSz8++Z>8W1Z4rQs1E~P@|8;m77$)#S*8jui
zC(rXUFr4ORU^vCkz>vVtz@X00z@W^}z_5;wfx(QAf#DY~1H(6728MUM3=D^P85s8S
zGBE7pWnkFA%fK*?mw{m#)XaX6xm4`mjv6%@0;3@?8UmvsFd71*Aut*OqaiRF0;3@?
z8Un*51Zu4npi^TC3eedgeSHNtPahwJ5D(8_2FUbVaY<%c8fYd7rl`<R*DgJ^Bq%j4
zB;T{BD77rLs5mnzC)G7CwJ5z(LsP+4!9Wi*ww9VxoC+F|ftXE%Hjt76P!BXS8D<#B
z9Pp?Oh|mD_n!%R)CMFl<L%Qf1nqZND#F7lR{G6Qp^31$+A5f3gGcUO)H8(Y{B&4z+
z733^ig#u(nzKMB>>8V9lRxX~w!JYw*A)fwz@veTZLGF>7njm+B+z+0Cq>;S=K90d5
zp3X2k74nN1NC;z4V5Ov{<s*j_Xz-OrVdU!O4|8ZWXs`}5Qo&;fh?rB*wpDOX%*{>I
NO35$r4)(KV006(xkMRHi
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index ea29b71df..937e8e80d 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -23,8 +23,6 @@ mlawNonLinearTVP::mlawNonLinearTVP(const int num,const double E,const double nu,
// by default, no temperature dependence
_temFunc_Sy0 = new constantScalarFunction(1.);
_temFunc_H = new constantScalarFunction(1.);
- _temFunc_Hb = new constantScalarFunction(1.);
- _temFunc_visc = new constantScalarFunction(1.);
};
@@ -39,17 +37,9 @@ mlawNonLinearTVP::mlawNonLinearTVP(const mlawNonLinearTVP& src): mlawNonLinearTV
_temFunc_Sy0 = src._temFunc_Sy0->clone();
_temFunc_H = NULL; // temperature dependence of hardening stress
- if (src._temFunc_H != NULL)
+ if (src._temFunc_H != NULL){
_temFunc_H = src._temFunc_H->clone();
-
- _temFunc_Hb = NULL; // temperature dependence of kinematic hardening
- if (src._temFunc_Hb != NULL)
- _temFunc_Hb = src._temFunc_Hb->clone();
-
- _temFunc_visc = NULL; // temperature dependence of viscosity
- if (src._temFunc_visc != NULL)
- _temFunc_visc = src._temFunc_visc->clone();
-
+ }
};
mlawNonLinearTVP& mlawNonLinearTVP::operator=(const materialLaw& source){
@@ -64,15 +54,9 @@ mlawNonLinearTVP& mlawNonLinearTVP::operator=(const materialLaw& source){
if (src->_temFunc_Sy0!= NULL)
_temFunc_Sy0 = src->_temFunc_Sy0->clone();
if(_temFunc_H != NULL) delete _temFunc_H; // temperature dependence of hardening stress
- if (src->_temFunc_H != NULL)
+ if (src->_temFunc_H != NULL){
_temFunc_H = src->_temFunc_H->clone();
- if(_temFunc_Hb != NULL) delete _temFunc_Hb; // temperature dependence of kinematic hardening
- if (src->_temFunc_Hb != NULL)
- _temFunc_Hb = src->_temFunc_Hb->clone();
- if(_temFunc_visc != NULL) delete _temFunc_visc; // temperature dependence of viscosity
- if (src->_temFunc_visc != NULL)
- _temFunc_visc = src->_temFunc_visc->clone();
-
+ }
}
return *this;
};
@@ -80,8 +64,6 @@ mlawNonLinearTVP& mlawNonLinearTVP::operator=(const materialLaw& source){
mlawNonLinearTVP::~mlawNonLinearTVP(){
if (_temFunc_Sy0 != NULL) delete _temFunc_Sy0; _temFunc_Sy0 = NULL;
if (_temFunc_H != NULL) delete _temFunc_H; _temFunc_H= NULL;
- if (_temFunc_Hb != NULL) delete _temFunc_Hb; _temFunc_Hb= NULL;
- if (_temFunc_visc != NULL) delete _temFunc_visc; _temFunc_visc= NULL;
};
// CHECK IF IPSTATE NEEDS the same definition as in mlawPowerYieldHyper --- WHAT??
@@ -126,20 +108,20 @@ void mlawNonLinearTVP::getModifiedMandel(const STensor3& C, const IPNonLinearTVP
M = 0.5*(corKir + temp2);
};
-void mlawNonLinearTVP::checkCommutavity(STensor3& commuteCheck,const STensor3& Ce, const STensor3& S, const IPNonLinearTVP *q1) const{
+void mlawNonLinearTVP::checkCommutavity(const STensor3& C, const STensor3& S, const IPNonLinearTVP *q1) const{
const STensor3& M = q1->getConstRefToModifiedMandel();
const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
- STensor3 temp1, temp2, temp3, temp4, temp6;
- STensorOperation::multSTensor3(Ce, S, temp1);
- STensorOperation::multSTensor3(S, Ce, temp2);
- STensorOperation::multSTensor3(corKir, Ce, temp3);
- STensorOperation::multSTensor3(Ce, corKir, temp4);
+ STensor3 temp1, temp2, temp3, temp4, temp5, temp6;
+ STensorOperation::multSTensor3(C, S, temp1);
+ STensorOperation::multSTensor3(S, C, temp2);
+ STensorOperation::multSTensor3(corKir, C, temp3);
+ STensorOperation::multSTensor3(C, corKir, temp4);
// Check S
- commuteCheck = temp1;
- commuteCheck -= temp2;
+ temp5 = temp1;
+ temp5 -= temp2;
// Check corKir (just for fun)
temp6 = temp3;
@@ -155,7 +137,7 @@ void mlawNonLinearTVP::getChabocheCoeffs(fullVector<double>& coeffs, const doubl
coeffs.resize(2);
- coeffs(0) = 1; //DEBUGGING
+ coeffs(0) = 0; // 1;
coeffs(1) = 0;
// add dependency on p later
@@ -188,16 +170,13 @@ void mlawNonLinearTVP::getYieldCoefficientTempDers(const IPNonLinearTVP *q1, con
double DDmDTT = (ddsigtdTT*sigc - sigt*ddsigcdTT)/(sigc*sigc) - 2*(dsigtdT*sigc- sigt*dsigcdT)/(sigc*sigc*sigc)*dsigcdT;
double Da1Dm = 3./sigc*(_n*pow(m,_n-1.)/(m+1.) - (pow(m,_n)-1.)/(m+1.)/(m+1.));
double DDa1Dmm = 3./sigc*(_n*(_n-1)*pow(m,_n-2.)/(m+1.) - _n*(pow(m,_n)-1.)/(m+1.)/(m+1.) - _n*pow(m,_n-1)/(m+1.)/(m+1.) + 2.*(pow(m,_n)-1.)/(m+1.)/(m+1.)/(m+1.));
- double term = _n*pow(m,_n) + _n*pow(m,_n-1.) - pow(m,_n);
- double dterm = pow(_n,2)*pow(m,_n-1.) + _n*(_n-1)*pow(m,_n-2.) - _n*pow(m,_n-1);
DcoeffsDT(2) = -_n*pow(sigc,-_n-1.)*dsigcdT;
DcoeffsDT(1) = Da1Dm*DmDT - 3.*(pow(m,_n)-1.)/(m+1.)/(sigc*sigc)*dsigcdT;
DcoeffsDT(0) = ((_n*pow(m,_n-1)+1.)/(m+1) - (pow(m,_n)+m)/(m+1.)/(m+1.))*DmDT;
DDcoeffsDTT(2) = -_n*pow(sigc,-_n-1.)*ddsigcdTT -_n*(-_n-1)*pow(sigc,-_n-2.)*dsigcdT*dsigcdT;
- DDcoeffsDTT(1) = -1/sigc*dsigcdT*DcoeffsDT(1) + 3/sigc*( (dterm/(m+1.)/(m+1.) -2*term/(m+1.)/(m+1.)/(m+1.))*DmDT*DmDT + term/(m+1.)/(m+1.)*DDmDTT
- + 1/sigc*( (1/sigc*pow(dsigtdT,2) - ddsigcdTT) * (pow(m,_n)-1.)/(m+1.) + dsigtdT * term/(m+1.)/(m+1.) ) );
+ DDcoeffsDTT(1) = DDa1Dmm*DmDT*DmDT + Da1Dm*DDmDTT + 3.*2.*(pow(m,_n)-1.)/(m+1.)/(sigc*sigc*sigc)*dsigcdT*dsigcdT - 3.*(pow(m,_n)-1.)/(m+1.)/(sigc*sigc)*ddsigcdTT;
DDcoeffsDTT(0) = ((_n*pow(m,_n-1)+1.)/(m+1) - (pow(m,_n)+m)/(m+1.)/(m+1.))*DDmDTT
+ ((_n*(_n-1)*pow(m,_n-2))/(m+1) + 2*(pow(m,_n)+m)/(m+1.)/(m+1.)/(m+1.))*DmDT*DmDT;
@@ -368,13 +347,13 @@ void mlawNonLinearTVP::getFreeEnergyTVM(const IPNonLinearTVP *q0, IPNonLinearTVP
};
void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double T0, const double& T,
- const STensor3& DphiPDF, std::vector<double>* ddRdTT, std::vector<STensor3>* ddRdFdT) const{
+ const STensor3& dPhiPdF, std::vector<double>* ddRdTT, std::vector<STensor3>* ddRdFdT) const{
const double& Gamma_0 = q0->_Gamma;
const double& Gamma = q1->_Gamma;
const double& DgammaDT = q1->_DgammaDT;
const double& dGammaDT_0 = q0->_DGammaDT;
- const double& dGammaDT = q1->_DGammaDT; // = 0. for debugging //DEBUGGING
+ const double& dGammaDT = q1->_DGammaDT;
const STensor3& DgammaDF = q1->_DgammaDF;
const STensor3& DGammaDF = q1->_DGammaDF;
const double& gamma0 = q0->_epspbarre;
@@ -393,10 +372,8 @@ void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNo
STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
STensorOperation::decomposeDevTr(q0->_DModMandelDT,dDevMeDT,dTrMeDT_0);
- trMe /= 3; trX1 /= 3; dTrMeDT /= 3; dTrXdT /= 3;
- trMe_0 /= 3; trX1_0 /= 3; dTrMeDT_0 /= 3; dTrXdT_0 /= 3;
-
- // dTrMeDT = 0.; dTrXdT = 0.; // for debugging only //DEBUGGING
+ trMe *= 1/3; trX1 *= 1/3; dTrMeDT *= 1/3; dTrXdT *= 1/3;
+ trMe_0 *= 1/3; trX1_0 *= 1/3; dTrMeDT_0 *= 1/3; dTrXdT_0 *= 1/3;
// get Dgamma
double Dgamma = gamma1 - gamma0;
@@ -443,14 +420,14 @@ void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNo
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
- dRdF[0](i,j) = -Da(1)*(trMe-trX1)*DgammaDF(i,j)*sigc - a(1)*(DphiPDF(i,j)*sigc + (trMe-trX1)*Hc*DgammaDF(i,j));
- dRdF[1](i,j) = -Da(0)*DgammaDF(i,j)*sigc - a(0)*Hc*DgammaDF(i,j);
+ dRdF[0](i,j) = -Da(1)*(trMe-trX1)*DgammaDF(i,j)*sigc - a(1)*(dPhiPdF(i,j)*sigc + (trMe-trX1)*Hc);
+ dRdF[1](i,j) = -Da(0)*DgammaDF(i,j)*sigc - a(0)*Hc;
}
if (Gamma>0 and etaOverDt>0){
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
- dRdF[2](i,j) = - _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(eta*DGammaDF(i,j)/this->getTimeStep())*sigc - pow(eta*Gamma/this->getTimeStep(),_p) * Hc*DgammaDF(i,j);
+ dRdF[2](i,j) = - _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(eta*DGammaDF(i,j)/this->getTimeStep())*sigc - pow(eta*Gamma/this->getTimeStep(),_p) * Hc;
}
}
@@ -469,12 +446,12 @@ void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNo
}
if (delT>1e-10){
- double DphiPDT = dTrMeDT- dTrXdT; double DphiPDT_0 = dTrMeDT_0 - dTrXdT_0;
- double DDphiPDTT = (DphiPDT-DphiPDT_0)/delT; // (DphiPDT - (trMe-trMe_0+trX1_0-trX1)/delT)/delT; // 2*DphiPDT/delT - 2*(trMe-trMe_0+trX1_0-trX1)/pow(delT,2); // DEBUGGING
+ double DphiPDT = dTrMeDT- dTrXdT;
+ double DDphiPDTT = 2*DphiPDT/delT - 2*(trMe-trMe_0+trX1_0-trX1)/pow(delT,2);
ddRdTT->at(0) -= a(1) * DDphiPDTT * sigc;
if (Gamma>0 and etaOverDt>0){
- double DDGammaDTT = (dGammaDT-dGammaDT_0)/delT; // 2*dGammaDT/delT - 2*(Gamma-Gamma_0)/pow(delT,2);
+ double DDGammaDTT = 2*dGammaDT/delT - 2*(Gamma-Gamma_0)/pow(delT,2);
ddRdTT->at(2) -= _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(eta*DDGammaDTT)/this->getTimeStep() * sigc;
}
}
@@ -497,7 +474,7 @@ void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNo
};
void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double T0, const double& T,
- const STensor3& dFpdT, const STensor43& dFpdF, const STensor3& DphiPDF,
+ const STensor3& dFpdT, const STensor43& dFpdF, const STensor3& dPhiPdF,
double *Wm, STensor3 *dWmdF, double *dWmdT) const{
const double& Gamma = q1->_Gamma;
@@ -512,7 +489,7 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
const STensor3& Me = q1->_ModMandel;
const STensor3& X1 = q1->_backsig;
const STensor3& X0 = q0->_backsig;
- const STensor3& dXdT = q1->_DbackSigDT; // Just for debugging = 0. //DEBUGGING
+ const STensor3& dXdT = q1->_DbackSigDT;
const STensor3& dXdT_0 = q0->_DbackSigDT;
const STensor43& dXdF_0 = q0->_DbackSigDF;
const STensor43& dXdF = q1->_DbackSigDF;
@@ -546,8 +523,7 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
STensorOperation::inverseSTensor3(Fp1,Fpinv);
Fp_dot = Fp1;
Fp_dot -= Fp0;
- Fp_dot *= 1/this->getTimeStep();
- STensorOperation::multSTensor3(Fp_dot,Fpinv,Dp);
+ STensorOperation::multSTensor3(Dp,Fp_dot,Fpinv);
// get alpha
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
@@ -566,12 +542,12 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
std::vector<double> ddIsoHardForcedTT; ddIsoHardForcedTT.resize(3,0.);
std::vector<STensor3> ddIsoHardForcedFdT; ddIsoHardForcedFdT.resize(3,0.);
- getIsotropicHardeningForce(q0,q1,T0,T,DphiPDF,&ddIsoHardForcedTT,&ddIsoHardForcedFdT);
+ getIsotropicHardeningForce(q0,q1,T0,T,dPhiPdF,&ddIsoHardForcedTT,&ddIsoHardForcedFdT);
const std::vector<double>& IsoHardForce = q1->_IsoHardForce;
const std::vector<double>& dIsoHardForcedT = q1->_dIsoHardForcedT;
const std::vector<STensor3>& dIsoHardForcedF = q1->_dIsoHardForcedF;
- double R(0.), dRdT(0.), ddRdTT(0.);
+ double R, dRdT, ddRdTT;
STensor3 dRdF, ddRdFdT; STensorOperation::zero(dRdF); STensorOperation::zero(ddRdFdT);
for (int i=0; i<3; i++){
@@ -579,7 +555,7 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
dRdT += dIsoHardForcedT[i];
ddRdTT += ddIsoHardForcedTT[i];
dRdF += dIsoHardForcedF[i];
- ddRdFdT += ddIsoHardForcedFdT[i]; //DEBUGGING
+ ddRdFdT += ddIsoHardForcedFdT[i];
}
R -= T*dRdT;
@@ -598,7 +574,6 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
// Some Generic Conversion Tensors
static STensor43 dDpdFp, dFpinvdFp;
- static STensor3 dFpinvdT;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
@@ -607,7 +582,7 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
dFpinvdFp(i,j,k,l) = 0.;
for (int m=0; m<3; m++)
for (int s=0; s<3; s++)
- dFpinvdFp(i,j,k,l) -= Fpinv(i,m)*_I4(m,s,k,l)*Fpinv(s,j);
+ dFpinvdFp(i,j,k,l) -= Fpinv(i,m)*_I4(m,j,k,s)*Fpinv(s,l);
}
for (int i=0; i<3; i++)
@@ -616,17 +591,9 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
for (int l=0; l<3; l++){
dDpdFp(i,j,k,l) = 0.;
for (int m=0; m<3; m++)
- dDpdFp(i,j,k,l) += 1/this->getTimeStep()*_I4(i,m,k,l)*Fpinv(m,j) + Fp_dot(i,m)*dFpinvdFp(m,j,k,l);
+ dDpdFp(i,j,k,l) += 1/this->getTimeStep()*_I4(i,j,k,m)*Fpinv(m,l) + Fp_dot(i,m)*dFpinvdFp(m,j,k,l);
}
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- dFpinvdT(i,j) = 0.;
- for (int m=0; m<3; m++)
- for (int s=0; s<3; s++)
- dFpinvdT(i,j) -= Fpinv(i,m)*dFpdT(m,s)*Fpinv(s,j);
- }
-
// dDpdF and dDpdT
static STensor43 dDpdF;
// STensorOperation::zero(dDpdF);
@@ -636,9 +603,10 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dDpdT(i,j) = 0.;
- for (int m=0; m<3; m++){
- dDpdT(i,j) += dFpdT(i,m)*Fpinv(m,j)*1/this->getTimeStep() + Fp_dot(i,m)*dFpinvdT(m,j); // dDpdFp(i,j,k,l)*dFpdT(k,l); // maybe use Dp = Fp_dot Fpinv
- }}
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++)
+ dDpdT(i,j) += dDpdFp(i,j,k,l)*dFpdT(k,l);
+ }
// dXdF is available
// const STensor43 dXdF = q1->_DbackSigDF;
@@ -652,10 +620,10 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
if (delT>1e-10){
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
- ddXdTT(i,j) += (dXdT(i,j)-dXdT_0(i,j))/delT; // 2*dXdT(i,j)/(T-T0) - 2*(X1(i,j)-X0(i,j))/pow((T-T0),2); // = 0. in DEBUGGING
+ ddXdTT(i,j) += 2*dXdT(i,j)/(T-T0) - 2*(X1(i,j)-X0(i,j))/pow((T-T0),2);
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
- ddXdTF(i,j,k,l) += (dXdF(i,j,k,l)-dXdF_0(i,j,k,l))/delT; // = 0. in DEBUGGING
+ ddXdTF(i,j,k,l) += (dXdF(i,j,k,l)-dXdF_0(i,j,k,l))/(T-T0);
}
}
}
@@ -689,8 +657,8 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
dmechSourceTVPdF(i,j) = - (dRdF(i,j) - T*ddRdFdT(i,j))*Dgamma/this->getTimeStep() - R*DgammaDF(i,j)/this->getTimeStep();
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
- dmechSourceTVPdF(i,j) += dMedF(i,j,k,l)*Dp(i,j) + Me(i,j)*dDpdF(i,j,k,l)
- - (dXdF(i,j,k,l) - T*ddXdTF(i,j,k,l))*alpha1(i,j) - tempX(i,j)*dAlphadF(i,j,k,l)/this->getTimeStep();
+ dmechSourceTVPdF(i,j) += dMedF(i,j,k,l)*Dp(k,l) + Me(i,j)*dDpdF(i,j,k,l)
+ - (dXdF(i,j,k,l) - T*ddXdTF(i,j,k,l))*alpha1(k,l) - tempX(i,j)*dAlphadF(i,j,k,l)/this->getTimeStep();
}
}
}
@@ -1077,24 +1045,30 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
DlnDCe = DlnDCepr;
// update A, B
+ // double Ge, Ke;
double Ke, Ge = 0.; double Kde, Gde = 0.; double KTsum, GTsum = 0.;
double DKe, DGe = 0.; double DKde, DGde = 0.; double DKDTsum, DGDTsum = 0.;
ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T);
getTVEdCorKirDT(q0,q1,T0,T);
+ // NEW
static STensor3 Me, Mepr, Kepr; // Ke is corKir
getModifiedMandel(Ce,q0,q1);
Me = q1->_ModMandel;
Mepr = Me;
Kepr = q1-> _kirchhoff;
+ // const STensor3& dKeprDT = q1->getConstRefToDcorKirDT();
+ // NEW
+ // NEW
static STensor3 devXn;
static double trXn;
STensorOperation::decomposeDevTr(q0->_backsig,devXn,trXn); // needed for chaboche model
+ // NEW
static STensor3 PhiPr;
- PhiPr = q1->_ModMandel;
+ PhiPr = q1->_ModMandel; // _kirchhoff; //NEW
PhiPr -= q1->_backsig;
static STensor3 devPhipr,devPhi; // effective dev stress predictor
@@ -1113,7 +1087,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
// double Gamma = 0.; // flow rule parameter
- double& Gamma = q1->_Gamma; // flow rule parameter // still initialised to zero
+ double& Gamma = q1->_Gamma; // flow rule parameter
double PhiEq = PhiEqpr; // current effective deviator
@@ -1122,14 +1096,13 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static fullVector<double> a(3), Da(3); // yield coefficients and derivatives in respect to plastic deformation
this->getYieldCoefficients(q1,a);
- double Hb(0.), dHb(0.), dHbdT(0.), ddHbddT(0.);
+ double Hb(0.), dHbdT(0.), ddHbddT(0.);
if (q1->_ipKinematic != NULL){
Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
- dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter
dHbdT = Hb * q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
ddHbddT = Hb * q1->_ipKinematic->getDDRDDT();
}
-
+
//a.print("a init");
static STensor3 devN; // dev part of yield normal
@@ -1149,15 +1122,17 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
double dDgammaDGamma = 0.;
double Dgamma = 0.; // eqplastic strain
- double Cx(0.), Px(0.);
+// NEW
static fullVector<double> m(2);
getChabocheCoeffs(m,0,q1);
- getPlasticPotential(ptilde, PhiEq, Px);
-
- Cx = m(0)*(-q0->_epspbarre + 0.) + m(1)*Px - 1/Hb * dHbdT * (T-T0) + 1; // CHECK -> CHANGED m(0)*(q0->_epsbarre + Dgamma) and no del
+ double Cx(0.), Px(0.);
+ getPlasticPotential(ptilde, PhiEq, Px);
+ Cx = m(0)*(q0->_epspbarre + Dgamma) + m(1)*Px - 1/Hb * dHbdT + 1;
double Gt = Ge + pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
double Kt = Ke + pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
+// NEW
+
if (q1->dissipationIsBlocked()){
q1->getRefToDissipationActive() = false;
@@ -1169,6 +1144,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
int ite = 0;
int maxite = 100; // maximal number of iters
+
//Msg::Error("plasticity occurs f = %e",f);
//double f0 = fabs(f);
@@ -1179,19 +1155,22 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
_viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT);
double etaOverDt = eta/this->getTimeStep();
+ // NEW
double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A);
+ // double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A); // -> the same
+ // NEW
dDgammaDGamma = kk*(A+Gamma*dAdGamma); // mistake in the paper (VD 2016)
this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*ptilde -Da(0);
-
+ // NEW
double dCxdDgamma = m(0);
double Stemp = STensorOperation::doubledot(devPhi,devXn);
- double He = (PhiEq*6*Gamma*kk*kk*Hb - Stemp*3.*1./PhiEq)/(2*Cx*Cx*u) * dCxdDgamma; // CHECK -> CHANGE
- double Hp = ptilde*(2*_b*Gamma*kk*kk*Hb - trXn)/(3*Cx*Cx*v) * dCxdDgamma; // CHECK
+ double He = (PhiEq*6*Gamma*kk*kk*Hb - Stemp*3.*1./PhiEq)/(2*Cx*Cx*u) * dCxdDgamma; // CHANGE THIS DERIVATIVE -> tensor subtracting scalar
+ double Hp = ptilde*(2*_b*Gamma*kk*kk*Hb - trXn)/(3*Cx*Cx*v) * dCxdDgamma;
dfdDgamma += a(2)*_n*pow(PhiEq,(_n-1))*He - a(1)*Hp;
-
+ // NEW
if (Gamma>0 and etaOverDt>0)
dfdDgamma -= _p*pow(etaOverDt,_p-1.)*Deta/this->getTimeStep()*pow(Gamma,_p); // THIS term is absent in the paper!! WHAT??
@@ -1199,7 +1178,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
if (Gamma>0 and etaOverDt>0)
DfDGamma -=pow(etaOverDt,_p)*_p*pow(Gamma,_p-1.);
- double dGamma = -f/DfDGamma; // WHAT
+ double dGamma = -f/DfDGamma; // WHAT
if (Gamma + dGamma <=0.){
Gamma /= 2.; // WHAT
@@ -1211,11 +1190,14 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
u = 1.+6.*Gt*Gamma;
v = 1.+2.*_b*Kt*Gamma;
-
+ // NEW
+ // PhiEq = PhiEqpr/u; // -> original
+ // ptilde = ptildepr/v; // -> original
devPhi = (devPhipr + devXn * (1-1/Cx));
devPhi *= 1./u;
PhiEq = sqrt(1.5*devPhi.dotprod());
ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v;
+ // NEW
A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
Dgamma = kk*Gamma*A;
@@ -1228,7 +1210,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
//a.print("a update");
f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
- double viscoTerm = etaOverDt*Gamma;
+ double viscoTerm = etaOverDt*Gamma; // I THINK temperature dependency of viscosity should be present in eta, but this relation remains the same
if (Gamma>0 and etaOverDt>0) f-= pow(viscoTerm,_p);
ite++;
@@ -1247,9 +1229,11 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
q1->_DgammaDt = Dgamma/this->getTimeStep();
// update effective stress tensor
+ // NEW
devPhi = (devPhipr + devXn * (1-1/Cx));
- devPhi *= 1./u;
- ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v;
+ devPhi *= 1./u; //NEW
+ ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v; //NEW
+ // NEW
// update normal
devN = devPhi;
@@ -1282,39 +1266,61 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// update A, B
ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T); //,false);
getTVEdCorKirDT(q0,q1,T0,T); // update dCorKirdT
- // DKDTsum and DGDTsum = sum of bulk/shear moduli derivatives wrt T
+ // DKDTsum and DGDTsum = sum of bulk/shear moduli derivatives wrt T
// backstress
static STensor3 DB; // increment
DB = (N); // increment
- DB *= (pow(kk,1)*Hb*Gamma); // pow(kk,2) CHANGE
+ DB *= (kk*Hb*Gamma); // pow(kk,2) CHANGE
+ // NEW
DB += q0->_backsig;
DB *= 1/Cx;
DB -= q0->_backsig;
+ // NEW
q1->_backsig += DB; // update
+
+ /*
+ // corotationaal Kirchhoff stress
+ q1->_kirchhoff = devPhi;
+ q1->_kirchhoff += q1->_backsig;
+
+ q1->_kirchhoff(0,0) += (ptilde);
+ q1->_kirchhoff(1,1) += (ptilde);
+ q1->_kirchhoff(2,2) += (ptilde);
+ */
+ // NEW
q1->_ModMandel = devPhi;
q1->_ModMandel += q1->_backsig;
q1->_ModMandel(0,0) += (ptilde);
q1->_ModMandel(1,1) += (ptilde);
q1->_ModMandel(2,2) += (ptilde);
-
+
+ // NEW
+
}
else{
q1->getRefToDissipationActive() = false;
}
}
-
- // update Stress
+
+
const STensor3& KS = q1->_kirchhoff;
- getModifiedMandel(Ce, q0, q1); // update Mandel
+ // second Piola Kirchhoff stress
+ // static STensor3 S;
+ // STensorOperation::multSTensor3STensor43(KS,DlnDCe,S);
+
+ // NEW
+ getModifiedMandel(Ce, q0, q1); // update Mandel // this step is probably wrong -> Mandel is already being updated above in the NR loop. Should corKir come from Mandel then?
+ // WHAT ???-
const STensor3& MS = q1->_ModMandel;
- static STensor3 S, Ceinv; // get 2nd PK
+ // second Piola Kirchhoff stress
+ static STensor3 S, Ceinv;
STensorOperation::inverseSTensor3(Ce,Ceinv);
- STensorOperation::multSTensor3(Ceinv,q1->_ModMandel,S);
-
- // first PK
+ STensorOperation::multSTensor3(Ceinv,q1->_ModMandel,S);
+ // NEW
+
for(int i=0; i<3; i++)
for(int j=0; j<3; j++){
P(i,j) = 0.;
@@ -1323,18 +1329,18 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
}
- // check Commutavity
- static STensor3 Check;
- checkCommutavity(Check,Ce,S,q1);
-
+
// defo energy
// q1->getRefToElasticEnergy()=deformationEnergy(*q1);
q1->getRefToViscousEnergyPart()=viscousEnergy(*q0,*q1)+q0->getConstRefToViscousEnergyPart();
q1->getRefToPlasticEnergy() = q0->plasticEnergy();
-
if (Gamma > 0){
+ // NEW
+ // double dotKSN = dot(KS,N);
+ //q1->getRefToPlasticEnergy() += Gamma*dotKSN;
double dotMSN = dot(MS,N);
q1->getRefToPlasticEnergy() += Gamma*dotMSN; // = Dp:Me
+ // NEW
}
@@ -1355,13 +1361,56 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
STensorOperation::multSTensor3SecondTranspose(Finv,Finv,Cinv);
STensorOperation::multSTensor3SVector3(Cinv,gradT,fluxT);
fluxT *= (-KThConT*J);
-
- // ThermSrc and MechSrc after dPdF and dPdT - But need the following tensors for the mechSrcTVP function call
- static STensor3 DphiPDF;
- STensorOperation::zero(dFpdF);
- STensorOperation::zero(dFpdT);
+ // thermalEnergy -> assuming stiff = false
+ double CpT;
+ if (stiff){
+ double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
+ getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
+ CpT = -T*DDpsiTVMdTT;
+ }
+ else{
+ getCp(CpT,T);
+ }
+
+ q1->_thermalEnergy = CpT*T;
+
+ // thermalSource
+ if (this->getTimeStep() > 0.){
+ thermalSource = -CpT*(T-T0)/this->getTimeStep();
+ }
+ else
+ thermalSource = 0.;
+
+ // NEW - mechSourceTVP
+ double& Wm = q1->getRefToMechSrcTVP();
+ double& Wm_TVE = q1->getRefToMechSrcTVE();
+
+ // Have to do this to update DcorKirDT in IP, and mechSourceTVE and derivatives
+ mechanicalSource = 0.;
+ mlawNonLinearTVE::getMechSourceTVE(q0,q1,T0,T,DKDTsum,DGDTsum,_I4,&Wm_TVE);
+ // STensor3 temp101 = q1->_DcorKirDT;
+ mechanicalSource += Wm_TVE;
+
+ STensor3 P_TVE(0.);
+ SVector3 fluxT_TVE(0.);
+ double thermalSource_TVE(0.);
+ double mechanicalSource_TVE(0.);
+ /* mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(q0->getConstRefToFe(), q1->getConstRefToFe(), P_TVE, q0, q1, Tangent, dFedF, dFedT,
+ T0, T, gradT0, gradT, fluxT_TVE, dPdT, dfluxTdgradT, dfluxTdT, dfluxTdF,
+ thermalSource_TVE, dthermalSourcedT, dthermalSourcedF,
+ mechanicalSource, dmechanicalSourcedT, dmechanicalSourceF,
+ stiff);*/
+
+
+ getMechSourceTVP(q0,q1,T0,T,_I,_I4,_I,&Wm);
+ mechanicalSource += Wm;
+
+ // freeEnergy and elastic energy
+ double& psiTVM = q1->getRefToFreeEnergyTVM();
+ getFreeEnergyTVM(q0,q1,T0,T,&psiTVM,NULL);
+ q1->_elasticEnergy = mlawNonLinearTVE::freeEnergyMechanical(*q1,T); // deformationEnergy(*q1,T);
// I didnt make this
if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
@@ -1405,6 +1454,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
//2. get DCeprDCepr and DCeinvprDCepr
static STensor43 DCeprDCepr, DCeinvprDCepr;
DCeprDCepr = _I4;
+
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
@@ -1551,7 +1601,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
temp8 *= (ptildepr + 1/3*trXn*(1-1/Cx));
temp8 *= 1/pow(v,2);
- DphiPDCepr += temp7;
+ DphiPDCepr -= temp7;
DphiPDCepr -= temp8;
// DdevphiDCepr
@@ -1563,15 +1613,15 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
STensorOperation::zero(temp13);
STensorOperation::zero(temp14);
- temp12 = devXn;
- temp12 *= 1/pow(Cx,2);
+ temp13 = devXn;
+ temp13 *= 1/pow(Cx,2);
STensorOperation::prod(temp12, dCxdCepr, 1., temp10);
temp10 *= 1/u;
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
- temp13(i,j) = devPhipr(i,j) + (1-1/Cx)*devXn(i,j);
+ temp13(i,j) = devPhipr(i,j)+ (1-1/Cx)*devXn(i,j);
temp14(i,j) = 6*DGDCepr(i,j)*( Ge + pow(kk,1)*Hb/(2*Cx) ) - 6*Gamma*(pow(kk,1)*Hb/(2*pow(Cx,2))*dCxdCepr(i,j)); // pow(kk,2) CHANGE
}
}
@@ -1579,7 +1629,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
STensorOperation::prod(temp13, temp14, 1., temp11);
temp11 *= 1/pow(u,2);
- DdevphiDCepr += temp10;
+ DdevphiDCepr -= temp10;
DdevphiDCepr -= temp11;
//9. get DtrNDCepr DdevNdCepr
@@ -1734,18 +1784,17 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// Done below -> need dCedF
STensor43& dXdF = q1->getRefToDbackStressdF();
- static STensor3 dTrXdF;
-
+ static STensor3 dTrXdF, DphiPDF;
STensorOperation::multSTensor43(dXdCepr,CeprToF,dXdF);
- // DphiPDF - use this for the pressure term in R
+ // DphiPDCepr - use this for the pressure term in R
// for the pressure term in R
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
DphiPDF(i,j) = 0.;
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
- DphiPDF(i,j) += DphiPDCepr(k,l)*CeprToF(k,l,i,j);
+ DphiPDF(i,j) = DphiPDCepr(k,l)*CeprToF(k,l,i,j);
}
}
@@ -1759,8 +1808,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
}
else{
STensorOperation::zero(DgammaDF);
- STensorOperation::zero(DGammaDF);
- // STensorOperation::zero(dFpdF);
+ STensorOperation::zero(dFpdF);
}
// 16. Everything else and Tangent
@@ -1925,15 +1973,15 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static double temp18;
static STensor3 temp19;
- dPhiPdT = dPhipprDT/v + 1/3*trXn/pow(Cx,2)*dCxdT/v; // no correction to this
+ dPhiPdT = dPhipprDT/v - 1/3*trXn/pow(Cx,2)*dCxdT/v; // no correction to this
temp18 = (ptildepr+1/3*trXn*(1-1/Cx))*(2*_b*Gamma*dKxdT)/pow(v,2); // Gamma derivatives will change this
dPhiPdT -= temp18;
- temp19 = (devPhipr + devXn*(1-1/Cx))*(6*Gamma*dGxdT);
+ temp19 = (devPhipr+devXn*(1-1/Cx))*(6*Gamma*dGxdT);
temp19 *= 1/pow(u,2);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
- dDevPhiDT(i,j) = (dDevPhiprDT(i,j) + devXn(i,j)/pow(Cx,2)*dCxdT)/u - temp19(i,j);
+ dDevPhiDT(i,j) = (dDevPhiprDT(i,j) - devXn(i,j)/pow(Cx,2)*dCxdT)/u - temp19(i,j);
}
}
@@ -1947,6 +1995,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static double DgammaDT;
static double DtrNdT;
static STensor3 DdevNdT;
+ static STensor3 dFpDT;
double eta(0.),Deta(0.),DetaDT(0.);
if (_viscosity != NULL)
@@ -1960,13 +2009,8 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
fullVector<double> DaDT(3), DDaDTT(3);
getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
- double a1, a2, a3;
- double a4, a5, a6;
- a1 = a(0); a2 = a(1); a3 = a(2);
- a4 = DaDT(0); a5 = DaDT(1); a6 = DaDT(2);
-
dAdT = (4.*_b*_b*ptilde/(A*3.))*dPhiPdT + (6./A)*PhiEq*dPhiEdT;
- dfdT = DaDT(2)*pow(PhiEq,_n) + a(2)*_n*pow(PhiEq,_n-1.)*dPhiEdT - DaDT(1)*ptilde - a(1)*dPhiPdT - DaDT(0);
+ dfdT = DaDT(2)*pow(PhiEq,_n) + a(2)*_n*pow(PhiEq,_n-1)*dPhiEdT - DaDT(1)*ptilde - a(1)*dPhiPdT - DaDT(0);
if (this->getTimeStep()>0){
dfdT -= (DetaDT*Gamma/this->getTimeStep())*_p*pow((eta*Gamma/this->getTimeStep()),(_p-1));
}
@@ -1999,12 +2043,12 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
DdevNdT = dDevPhiDT;
DdevNdT *= 3;
- // 10. get dFpdT
+ // 10. get dFpDT
static STensor3 temp20;
STensorOperation::zero(temp20);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
- temp20(i,j) += N(i,j)*dGammaDT + Gamma*DdevNdT(i,j) + 1/3*Gamma*_I(i,j)*DtrNdT;
+ temp20(i,j) += N(i,j)*dGammaDT + Gamma*DdevNdT(i,j) + Gamma/3.*_I(i,j)*DtrNdT;
static STensor43 EprFp0;
for (int i=0; i<3; i++)
@@ -2019,10 +2063,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
- dFpdT(i,j) = 0.;
+ dFpDT(i,j) = 0.;
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
- dFpdT(i,j) += EprFp0(i,j,k,l)*temp20(k,l);
+ dFpDT(i,j) += EprFp0(i,j,k,l)*temp20(k,l);
}
}
}
@@ -2048,16 +2092,12 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
else{
// elastic
STensorOperation::zero(DgammaDT);
- // STensorOperation::zero(dFpdT);
+ STensorOperation::zero(dFpDT);
STensorOperation::zero(DtrNdT);
STensorOperation::zero(DdevNdT);
STensorOperation::zero(dGammaDT);
}
- // 11.1
- q1->_DGammaDT = dGammaDT;
- q1->_DgammaDT = DgammaDT;
-
// 12. get dKcorDT
// done above
@@ -2069,7 +2109,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
for (int j=0; j<3; j++)
for (int p=0; p<3; p++)
for (int q=0; q<3; q++){
- DinvFpdT(i,j) -= Fpinv(i,p)*dFpdT(p,q)*Fpinv(q,j);
+ DinvFpdT(i,j) -= Fpinv(i,p)*dFpDT(p,q)*Fpinv(q,j);
}
STensorOperation::zero(dFedT);
@@ -2131,9 +2171,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
dPdT(i,j) += (dFedT(i,k)*S(k,l)*Fpinv(j,l) + Fe(i,k)*dSdT(k,l)*Fpinv(j,l) + Fe(i,k)*S(k,l)*DinvFpdT(j,l));
}
}
- }
- if(stiff){
// TVP - flux, thermalEnergy, thermalSource
if (!_thermalEstimationPreviousConfig){
@@ -2169,65 +2207,22 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
}
}
}
- }
-
- // thermSrc and MechSrc
-
- // thermalEnergy
- double CpT;
- if (stiff){
- double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
- getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
- // CpT = -T*DDpsiTVMdTT;
- getCp(CpT,T);
- }
- else{
- getCp(CpT,T);
- }
-
- q1->_thermalEnergy = CpT*T;
-
- // thermalSource
- if (this->getTimeStep() > 0.){
- thermalSource = -CpT*(T-T0)/this->getTimeStep();
- }
- else
- thermalSource = 0.;
-
- // mechSourceTVP
- double& Wm = q1->getRefToMechSrcTVP();
- double& Wm_TVE = q1->getRefToMechSrcTVE();
-
- mechanicalSource = 0.;
- mlawNonLinearTVE::getMechSourceTVE(q0,q1,T0,T,DKDTsum,DGDTsum,_I4,&Wm_TVE); // mechSourceTVE
- mechanicalSource += Wm_TVE;
- getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,&Wm);
- mechanicalSource += Wm;
-
- // freeEnergy and elastic energy
- double& psiTVM = q1->getRefToFreeEnergyTVM();
- getFreeEnergyTVM(q0,q1,T0,T,&psiTVM,NULL);
-
- q1->_elasticEnergy = mlawNonLinearTVE::freeEnergyMechanical(*q1,T); // deformationEnergy(*q1,T);
-
-if (stiff){
-
const double& DDpsiTVMdTT_0 = q0->getConstRefToDDFreeEnergyTVMdTT();
double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
- // double CpT = -T*DDpsiTVMdTT;
- // double CpT_0 = -T0*DDpsiTVMdTT_0;
+ double CpT = -T*DDpsiTVMdTT;
+ double CpT_0 = -T0*DDpsiTVMdTT_0;
+
+
+ // thermal energy
+ q1->_thermalEnergy = CpT*T;
// thermal source derivatives
- double DCpDT(0.);
- mlawNonLinearTVE::getCp(CpT,T,&DCpDT);
static STensor3 DCpDF;
STensorOperation::zero(DCpDF); // CHANGE for DCpDF
-
if (this->getTimeStep() > 0){
- dthermalSourcedT = -DCpDT*(T-T0)/this->getTimeStep() - CpT/this->getTimeStep();
- // dthermalSourcedT = -CpT/this->getTimeStep() - (CpT-CpT_0)/this->getTimeStep(); // thermalSource = -CpT*(T-T0)/this->getTimeStep();
+ dthermalSourcedT = -CpT/this->getTimeStep() - (CpT-CpT_0)/this->getTimeStep();
for(int i = 0; i< 3; i++){
for(int j = 0; j< 3; j++){
dthermalSourcedF(i,j) = -DCpDF(i,j)*(T-T0)/this->getTimeStep();
@@ -2269,7 +2264,7 @@ if (stiff){
double& dWmdT = q1->getRefTodMechSrcTVPdT();
STensor3& dWmdF = q1->getRefTodMechSrcTVPdF();
- getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,NULL,&dWmdF,&dWmdT);
+ getMechSourceTVP(q0,q1,T0,T,dFpDT,dFpdF,DphiPDF,NULL,&dWmdF,&dWmdT);
dmechanicalSourceF += dWmdF;
dmechanicalSourcedT += dWmdT;
@@ -2541,4 +2536,849 @@ void mlawNonLinearTVP::predictorCorrector_TVP_AssociatedFlow(const STensor3& F,
if (stiff){
// FILL THIS
}
-};
\ No newline at end of file
+};
+
+/*
+void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double T0, const double& T,
+ const STensor3& dFpdT, const STensor43& dFpdF, const STensor3& dPhiPdF,
+ double *Wm, STensor3 *dWmdF, double *dWmdT) const{
+
+ const double& DgammaDT = q1->_DgammaDT;
+ const double& dGammaDT = q1->_DGammaDT;
+ const STensor3& DgammaDF = q1->_DgammaDF;
+ const STensor3& DGammaDF = q1->_DGammaDF;
+ const double& gamma0 = q0->_epspbarre;
+ const double& gamma1 = q1->_epspbarre;
+ const STensor3& Fp0 = q0->_Fp;
+ const STensor3& Fp1 = q1->_Fp;
+ const STensor3& Me = q1->_ModMandel;
+ const STensor3& X1 = q1->_backsig;
+ const STensor3& X0 = q0->_backsig;
+ const STensor3& dXdT = q1->_DbackSigDT;
+ const STensor3& dXdT_0 = q0->_DbackSigDT;
+ const STensor43& dXdF_0 = q0->_DbackSigDF;
+ const STensor43& dXdF = q1->_DbackSigDF;
+ const STensor3& dMeDT = q1->_DModMandelDT;
+ const STensor3& dMeDT_0 = q0->_DModMandelDT;
+ const STensor43& dMedF = q1->_DModMandelDF;
+ const STensor43& dMedF_0 = q0->_DModMandelDF;
+
+ static STensor3 devMe, devX1, dDevXdT, dDevMeDT;
+ static double trMe, trX1, dTrXdT, dTrMeDT;
+
+ STensorOperation::decomposeDevTr(q1->_ModMandel,devMe,trMe);
+ STensorOperation::decomposeDevTr(q1->_backsig,devX1,trX1);
+ STensorOperation::decomposeDevTr(q1->_DbackSigDT,dDevXdT,dTrXdT);
+ STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
+
+ // get Dgamma
+ double Dgamma = gamma1 - gamma0;
+
+ // get Hb
+ double Hb(0.), dHb(0), dHbdT(0.), ddHbddT(0.);
+ if (q1->_ipKinematic != NULL){
+ Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
+ dHb = q1->_ipKinematic->getDDR();
+ dHbdT = q1->_ipKinematic->getDRDT();
+ ddHbddT = q1->_ipKinematic->getDDRDDT();
+ }
+
+ // get Dp
+ static STensor3 Fpinv, Fp_dot, Dp;
+ STensorOperation::inverseSTensor3(Fp1,Fpinv);
+ Fp_dot = Fp1;
+ Fp_dot -= Fp0;
+ STensorOperation::multSTensor3(Dp,Fp_dot,Fpinv);
+
+ // get alpha
+ double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
+ static STensor3 alpha0, alpha1;
+ alpha0 = X0; alpha1 = X1;
+ alpha0 *= 1/(pow(kk,2)*Hb); alpha1 *= 1/(pow(kk,2)*Hb); // how to get alpha0-current or previous Hb??? WHAT?
+ alpha1 -= alpha0;
+ alpha1 *= 1/this->getTimeStep();
+
+ // get TempX - convenient backStress tensor
+ STensor3 tempX;
+ tempX = X1;
+ tempX -= T*dXdT;
+
+ // get R
+ getIsotropicHardeningForce(q0,q1,T0,T,Gamma,dFpdT,dFpdF,dPhiPdF,ddRdTT,ddRdFdT);
+
+ double eta(0.),Deta(0.),DetaDT(0.);
+ fullVector<double> a(3), Da(3), DaDT(3), DDaDTT(3);
+ getYieldCoefficients(q1,a);
+ getYieldCoefficientDerivatives(q1,q1->_nup,Da);
+ getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
+
+ if (_viscosity != NULL)
+ _viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT);
+ double etaOverDt = eta/this->getTimeStep();
+ double viscoTerm = etaOverDt*Gamma;
+
+ const double& R0 = 0;
+
+ // double& R = q1->getRefToIsotropicHardeningForce(); //change
+ const double& dRdT_0 = 0;
+ // double& dRdT = q1->_dIsoHardForcedT;
+ double R(0.);
+ double dRdT(0.);
+ R = -a(1)/a(2) * (trMe-trX1);
+ dRdT = (-DaDT(1)/a(2) + a(1)/pow(a(2),2) *DaDT(2))*(trMe-trX1) + (-a(1)/a(2))*(dTrMeDT-dTrXdT);
+ R -= a(0)/a(2);
+ dRdT += (-DaDT(0)/a(2) + a(0)/pow(a(2),2)*DaDT(2));
+ if (Gamma>0 and etaOverDt>0){
+ R -= 1/a(2)*pow(eta*Gamma/this->getTimeStep(),_p);
+ dRdT += (1/pow(a(2),2)*DaDT(2)*pow(eta*Gamma/this->getTimeStep(),_p)
+ - 1/a(2)*_p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(DetaDT*Gamma/this->getTimeStep()
+ + eta*dGammaDT/this->getTimeStep()) );
+ }
+
+ R -= T*dRdT;
+
+ // mechSourceTVP
+ if(Wm!=NULL){
+ double mechSourceTVP = STensorOperation::doubledot(Me,Dp);
+ mechSourceTVP -= STensorOperation::doubledot(tempX,alpha1);
+ mechSourceTVP -= (R*Dgamma)/this->getTimeStep();
+ *Wm = mechSourceTVP;
+ // Second Term TBD - need dXdT, dRdT -> added above
+ }
+
+ // Some Generic Conversion Tensors
+ static STensor43 dDpdFp, dFpinvdFp;
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++){
+ dFpinvdFp(i,j,k,l) = 0.;
+ for (int m=0; m<3; m++)
+ for (int s=0; s<3; s++)
+ dFpinvdFp(i,j,k,l) -= Fpinv(i,m)*_I4(m,j,k,s)*Fpinv(s,l);
+ }
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++){
+ dDpdFp(i,j,k,l) = 0.;
+ for (int m=0; m<3; m++)
+ dDpdFp(i,j,k,l) += 1/this->getTimeStep()*_I4(i,j,k,m)*Fpinv(m,l) + Fp_dot(i,m)*dFpinvdFp(m,j,k,l);
+ }
+
+ // dDpdF and dDpdT
+ static STensor43 dDpdF;
+ // STensorOperation::zero(dDpdF);
+ STensorOperation::multSTensor43(dDpdFp, dFpdF, dDpdF);
+
+ static STensor3 dDpdT;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dDpdT(i,j) = 0.;
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++)
+ dDpdT(i,j) += dDpdFp(i,j,k,l)*dFpdT(k,l);
+ }
+
+ // dXdF is available
+ // const STensor43 dXdF = q1->_DbackSigDF;
+
+ // ddXdTF -> make it -- CHANGE!!
+ static STensor43 ddXdTF; // ddGammaDTF and ddQDTF are very difficult to estimate
+ static STensor3 ddXdTT; // ddGammaDTT and ddQDTT are very difficult to estimate
+ STensorOperation::zero(ddXdTF);
+ STensorOperation::zero(ddXdTT);
+ double delT = T-T0;
+ if (delT>1e-10){
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ ddXdTT(i,j) += 2*dXdT(i,j)/(T-T0) - 2*(X1(i,j)-X0(i,j))/pow((T-T0),2);
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++){
+ ddXdTF(i,j,k,l) += (dXdF(i,j,k,l)-dXdF_0(i,j,k,l))/(T-T0);
+ }
+ }
+ }
+
+ // dAlphadF - from dXdF
+ static STensor43 dAlphadF;
+ STensorOperation::zero(dAlphadF);
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++){
+ dAlphadF(i,j,k,l) += dXdF(i,j,k,l)*1/(pow(kk,2)*Hb) + X1(i,j)*DgammaDF(k,l)*1/(pow(kk*Hb,2))*dHb;
+ }
+
+ // dAlphadT - from dXdT
+ static STensor3 dAlphadT;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dAlphadT(i,j) = dXdT(i,j)*1/(pow(kk,2)*Hb) + X1(i,j)*DgammaDT*1/(pow(kk*Hb,2))*dHbdT;
+ }
+
+
+ // dRdF
+ // static STensor3 dRdF;
+ const STensor3& dRdF_0 = 0.; //q0->_dIsoHardForcedF;
+ // STensor3& dRdF = q1->_dIsoHardForcedF;
+ static STensor3 dRdF;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dRdF(i,j) = (-Da(1)/a(2) + a(1)/pow(a(2),2) *Da(2))*(trMe-trX1)*DgammaDF(i,j) + (-a(1)/a(2))*dPhiPdF(i,j) + (-Da(0)/a(2) + a(0)/pow(a(2),2)*Da(2))*DgammaDF(i,j);
+ }
+
+ if (Gamma>0 and etaOverDt>0){
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dRdF += (1/pow(a(2),2)*Da(2)*pow(eta*Gamma/this->getTimeStep(),_p)*DgammaDF(i,j)
+ - 1/a(2)*_p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(eta*DGammaDF(i,j)/this->getTimeStep()));
+ }
+ }
+
+ // ddRdTT, ddRdTF -> make all these -- CHANGE!!
+ static STensor3 ddRdTF;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ ddRdTF(i,j) = (dRdF(i,j)-dRdF_0(i,j))/(T-T0);
+ }
+ // STensorOperation::zero(ddRdTF); // ddGammaDTF is very difficult to estimate
+
+ static double ddRdTT;
+ ddRdTT = 2*dRdT/(T-T0) -2*(R-R0)/pow((T-T0),2); // 0.; // ddGammaDTT is very difficult to estimate
+
+
+ // dMedF is available
+ // const STensor43 dMedF = q1->_DModMandelDF;
+
+ if(dWmdF!=NULL){
+ STensor3 dmechSourceTVPdF;
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dmechSourceTVPdF(i,j) = - (dRdF(i,j) - T*ddRdTF(i,j))*Dgamma/this->getTimeStep() - R*DgammaDF(i,j)/this->getTimeStep();
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++){
+ dmechSourceTVPdF(i,j) += dMedF(i,j,k,l)*Dp(k,l) + Me(i,j)*dDpdF(i,j,k,l)
+ - (dXdF(i,j,k,l) - T*ddXdTF(i,j,k,l))*alpha1(k,l) - tempX(i,j)*dAlphadF(i,j,k,l)/this->getTimeStep();
+ }
+ }
+
+ *dWmdF = dmechSourceTVPdF;
+ }
+
+ if(dWmdT!=NULL){
+ double dmechSourceTVPdT(0.);
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dmechSourceTVPdT += dMeDT(i,j)*Dp(i,j) + Me(i,j)*dDpdT(i,j)
+ + T*ddXdTT(i,j)*alpha1(i,j) - tempX(i,j)*dAlphadT(i,j)/this->getTimeStep()
+ + T*ddRdTT*Dgamma/this->getTimeStep() - R*DgammaDT/this->getTimeStep() ;
+ }
+ *dWmdT = dmechSourceTVPdT;
+ }
+}; */
+
+/*void mlawNonLinearTVP::getFreeEnergyTVM(const IPNonLinearTVP *q0, const IPNonLinearTVP *q1, const double& T0, const double& T,
+ double *psiTVM, double *DpsiTVMdT, double *DDpsiTVMdTT) const{
+
+ // get some Things
+ const double& Gamma = q1->_Gamma;
+ const double& DgammaDT = q1->_DgammaDT;
+ const double& dGammaDT = q1->_DGammaDT;
+ const STensor3& Me = q1->_ModMandel;
+ const STensor3& X = q1->_backsig;
+ const STensor3& dXdT = q1->_DbackSigDT;
+ const STensor3& ddMedTT(0.); // CHANGE!!
+ const STensor3& ddXdTT(0.); // CHANGE!!
+
+ STensor3 devMe, devX, dDevXdT, dDevMeDT, ddDevXdTT, ddDevMeDTT;
+ double trMe, trX, dTrXdT, dTrMeDT, ddTrXdTT, ddTrMedTT;
+ STensorOperation::decomposeDevTr(q1->_ModMandel,devMe,trMe);
+ STensorOperation::decomposeDevTr(q1->_backsig,devX,trX);
+ STensorOperation::decomposeDevTr(q1->_DbackSigDT,dDevXdT,dTrXdT);
+ STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
+ STensorOperation::decomposeDevTr(ddXdTT,ddDevXdTT,ddTrXdTT);
+ STensorOperation::decomposeDevTr(ddMedTT,ddDevMeDTT,ddTrMedTT);
+
+ // Start
+ double psiVE, psiVP, psiT = 0.;
+ double DpsiVEdT, DpsiVPdT, DpsiTdT = 0.;
+ double DDpsiVEdTT, DDpsiVPdTT, DpsiTdTT = 0.;
+
+ double CpT, DCpDT, DDCpDTT;
+ mlawNonLinearTVE::getCp(CpT,T,&DCpDT,&DDCpDTT);
+
+ // psiT
+ psiT = CpT*((T-_Tinitial)-T*log(T/_Tinitial));
+ DpsiTdT = -CpT*log(T/_Tinitial) - psiT/CpT*DCpDT;
+ DpsiTdTT = -CpT/T - psiT/CpT*DDCpDTT;
+
+ //
+ // psiVP -> make IP for this CHANGE!!!
+
+ // get Isotropic Hardening stuff
+ double eta(0.),Deta(0.),DetaDT(0.);
+ fullVector<double> a(3), Da(3), DaDT(3), DDaDTT(3);
+ getYieldCoefficients(q1,a);
+ getYieldCoefficientDerivatives(q1,q1->_nup,Da);
+ getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
+
+ if (_viscosity != NULL)
+ _viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT);
+ double etaOverDt = eta/this->getTimeStep();
+ double viscoTerm = etaOverDt*Gamma;
+
+ double R1(0.), R2(0.), R3(0.), dR1dT(0.), dR2dT(0.), dR3dT(0.), ddR1dTT(0.), ddR2dTT(0.), ddR3dTT(0.);
+ R1 = -a(1)/a(2) * (trMe-trX);
+ dR1dT = (-DaDT(1)/a(2) + a(1)/pow(a(2),2) *DaDT(2))*(trMe-trX) + (-a(1)/a(2))*(dTrMeDT-dTrXdT);
+ R2 = -a(0)/a(2);
+ dR2dT = (-DaDT(0)/a(2) + a(0)/pow(a(2),2)*DaDT(2));
+ if (Gamma>0 and etaOverDt>0){
+ R3 -= 1/a(2)*pow(eta*Gamma/this->getTimeStep(),_p);
+ dR3dT += (1/pow(a(2),2)*DaDT(2)*pow(eta*Gamma/this->getTimeStep(),_p)
+ - 1/a(2)*_p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(DetaDT*Gamma/this->getTimeStep()
+ + eta*dGammaDT/this->getTimeStep()) );
+ }
+
+ // get Kinematic Hardening stuff
+ double Hb(0.), dHbdT(0.), ddHbddT(0.);
+ if (q1->_ipKinematic != NULL){
+ Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
+ dHbdT = q1->_ipKinematic->getDRDT();
+ ddHbddT = q1->_ipKinematic->getDDRDDT();
+ }
+
+ double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
+ static STensor3 alpha;
+ alpha = X;
+ alpha *= 1/(pow(kk,2)*Hb);
+
+ static STensor3 dAlphadT;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dAlphadT(i,j) = dXdT(i,j)*1/(pow(kk,2)*Hb) + X(i,j)*DgammaDT*1/(pow(kk*Hb,2))*dHbdT;
+ }
+
+
+ // finally
+ STensorOperation::doubleContractionSTensor3(X,alpha,psiVP);
+ psiVP *= 0.5;
+ psiVP += 0.5*( 1/_HR[0]*R1*a(1)/a(2) + 1/_HR[1]*pow(R2,2) + 1/_HR[2]*pow(R3,2) );
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ DpsiVPdT += 0.5*(dXdT(i,j)*alpha(i,j)+X(i,j)*dAlphadT(i,j));
+ DDpsiVPdTT += 0.5*(0.);
+ }
+ DpsiVPdT += 0.5*( 1/_HR[0]*(dR1dT*a(1)/a(2) + R1*(DaDT(1)/a(2) - a(1)/pow(a(2),2) *DaDT(2))) + 2/_HR[1]*R2*dR2dT + 2/_HR[2]*R3*dR3dT );
+
+ DDpsiVPdTT += 0.; // CHANGE!! 2nd derivative of R and X are difficult
+
+ //
+ // psiVE
+
+ double KT, GT, cteT, dKdT, ddKdTT, dGdT, ddGdTT, DcteDT, DDcteDTT; getK(KT,T,&dKdT,&ddKdTT); getG(GT,T,&dGdT,&ddGdTT); getAlpha(cteT,T,&DcteDT,&DDcteDTT);
+
+ const STensor3& E = q1->getConstRefToElasticStrain();
+ static STensor3 devKeinf, DdevKeinfDT, DDdevKeinfDTT, devEe;
+ static double trKeinf, DtrKeinfDT, DDtrKeinfDTT, trEe;
+ STensorOperation::decomposeDevTr(E,devEe,trEe);
+ mlawNonLinearTVE::corKirInfinity(devEe,trEe,T,devKeinf,trKeinf);
+
+ DtrKeinfDT = trKeinf*dKdT/KT - 3*KT*(DcteDT*(T-_Tinitial) + cteT);
+ DDtrKeinfDTT = trKeinf*ddKdTT/KT - 6*dKdT*(DcteDT*(T-_Tinitial) + cteT)- 3*KT*(DDcteDTT*(T-_Tinitial) + 2*DcteDT);
+ DdevKeinfDT = devKeinf;
+ DdevKeinfDT *= dGdT/GT;
+ DDdevKeinfDTT = devKeinf;
+ DDdevKeinfDTT *= ddGdTT/GT;
+
+ const std::vector<STensor3>& devOi = q1->getConstRefToDevViscoElasticOverStress();
+ const std::vector<double>& trOi = q1->getConstRefToTrViscoElasticOverStress();
+ const std::vector<STensor3>& devDOiDT = q1->getConstRefToDevDOiDT();
+ const std::vector<double>& trDOiDT = q1->getConstRefToTrDOiDT();
+ const std::vector<STensor3>& devDDOiDTT = q1->getConstRefToDevDDOiDTT();
+ const std::vector<double>& trDDOiDTT = q1->getConstRefToTrDDOiDTT();
+
+ psiVE = mlawNonLinearTVE::freeEnergyMechanical(*q1,T);
+
+ DpsiVEdT = 1/9 * (1/KT * trKeinf * DtrKeinfDT - 1/(2*pow(KT,2)) * dKdT * pow(trKeinf,2));
+ DDpsiVEdTT = 1/9 * (1/KT * (pow(DtrKeinfDT,2)+trKeinf*DDtrKeinfDTT) - 1/pow(KT,2)*dKdT*trKeinf*DtrKeinfDT
+ - 1/(2*pow(KT,2)) *(ddKdTT*pow(trKeinf,2) + 2*dKdT*trKeinf*DtrKeinfDT)
+ + 1/pow(KT,3)*pow(dKdT,2)*pow(trKeinf,2));
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ DpsiVEdT += 2/GT * devKeinf(i,j)*DdevKeinfDT(i,j) - 1/pow(GT,2) * dGdT * devKeinf(i,j)*devKeinf(i,j);
+ DDpsiVEdTT += 2/GT*(DdevKeinfDT(i,j)*DdevKeinfDT(i,j) + devKeinf(i,j)*DDdevKeinfDTT(i,j) - dGdT*devKeinf(i,j)*DdevKeinfDT(i,j))
+ - 1/pow(GT,2) * (ddGdTT*devKeinf(i,j)*devKeinf(i,j) + 2*dGdT*devKeinf(i,j)*DdevKeinfDT(i,j));
+ }
+ if ((_Ki.size() > 0) or (_Gi.size() > 0)){
+ for (int k=0; k<_Gi.size(); k++){
+ DpsiVEdT += 1/9 * 1/_Ki[k] * trOi[k] * trDOiDT[k];
+ DDpsiVEdTT += 1/9 * 1/_Ki[k] * (trOi[k]*trDDOiDTT[k] + pow(trDOiDT[k],2));
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ DpsiVEdT += 2/_Gi[k] * devOi[k](i,j) * devDOiDT[k](i,j);
+ DDpsiVEdTT += 2/_Gi[k] * (devOi[k](i,j) * devDDOiDTT[k](i,j) + pow(devDOiDT[k](i,j),2));
+ }
+ }
+ }
+
+ // FINALLY
+ if(psiTVM!=NULL){
+ *psiTVM = psiT + psiVE + psiVP;
+ }
+ if(DpsiTVMdT!=NULL){
+ *DpsiTVMdT = DpsiTdT + DpsiVEdT + DpsiVPdT;
+ }
+ if(DDpsiTVMdTT!=NULL){
+ *DDpsiTVMdTT = DpsiTdTT + DDpsiVEdTT + DDpsiVPdTT;
+ }
+
+};*/
+
+/* OLD PLasticity Correction step
+
+ * // compute elastic predictor
+ STensor3& Fp1 = q1->_Fp;
+ const STensor3& Fp0 = q0->_Fp;
+
+ Fp1 = Fp0; // plastic deformation tensor
+ q1->_epspbarre = q0->_epspbarre; // plastic equivalent strain
+ q1->_epspCompression = q0->_epspCompression;
+ q1->_epspTraction = q0->_epspTraction;
+ q1->_epspShear = q0->_epspShear;
+ q1->_backsig = q0->_backsig; // backstress tensor
+ q1->_DbackSigDT = q0->_DbackSigDT; // DbackSigDT
+ q1->_DgammaDt = 0.;
+
+
+ static STensor3 Fpinv, Ce, Fepr;
+ STensorOperation::inverseSTensor3(Fp1,Fpinv);
+ STensorOperation::multSTensor3(F,Fpinv,Fepr);
+ STensorOperation::multSTensor3FirstTranspose(Fepr,Fepr,Ce);
+
+ // NEW
+ static STensor3 Cepr, Ceinvpr;
+ Cepr = Ce;
+ STensorOperation::inverseSTensor3(Cepr,Ceinvpr);
+ // NEW
+
+ static STensor3 invFp0; // plastic predictor
+ invFp0= Fpinv;
+ STensor3& Fe = q1->_Fe;
+ Fe = Fepr;
+
+ static STensor43 DlnDCepr, DlnDCe;
+ static STensor63 DDlnDDCe;
+ static STensor3 expGN;
+ static STensor43 dexpAdA; // estimation of dexpA/dA
+
+ STensor3& Ee = q1->_Ee;
+ STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCepr,&DDlnDDCe);
+ Ee *= 0.5;
+ DlnDCe = DlnDCepr;
+
+ // update A, B
+ // double Ge, Ke;
+ double Ke, Ge = 0.; double Kde, Gde = 0.; double KTsum, GTsum = 0.;
+ double DKe, DGe = 0.; double DKde, DGde = 0.; double DKDTsum, DGDTsum = 0.;
+ ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T);
+ getTVEdCorKirDT(q0,q1,T0,T);
+
+ // NEW
+ static STensor3 Me, Mepr, Kepr; // Ke is corKir
+
+ getModifiedMandel(Ce,q0,q1);
+ Me = q1->_ModMandel;
+ Mepr = Me;
+ Kepr = q1-> _kirchhoff;
+ // const STensor3& dKeprDT = q1->getConstRefToDcorKirDT();
+ // NEW
+
+ // NEW
+ static STensor3 devXn;
+ static double trXn;
+ STensorOperation::decomposeDevTr(q0->_backsig,devXn,trXn); // needed for chaboche model
+ // NEW
+
+ static STensor3 PhiPr;
+ PhiPr = q1->_ModMandel; // _kirchhoff; //NEW
+ PhiPr -= q1->_backsig;
+
+ static STensor3 devPhipr,devPhi; // effective dev stress predictor
+ static double ptildepr,ptilde;
+ STensorOperation::decomposeDevTr(PhiPr,devPhipr,ptildepr);
+ ptildepr/= 3.;
+
+ ptilde = ptildepr; // current effective pressure -> first invariant
+ devPhi = devPhipr;
+
+ double PhiEqpr2 = 1.5*devPhipr.dotprod();
+ double PhiEqpr = sqrt(PhiEqpr2); // -> second invariant
+
+ // plastic poisson ratio
+ q1->_nup = (9.-2.*_b)/(18.+2.*_b);
+ double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
+
+ // double Gamma = 0.; // flow rule parameter
+ double& Gamma = q1->_Gamma; // flow rule parameter
+ double PhiEq = PhiEqpr; // current effective deviator
+
+
+ // hardening
+ this->hardening(q0,q1,T);
+ static fullVector<double> a(3), Da(3); // yield coefficients and derivatives in respect to plastic deformation
+ this->getYieldCoefficients(q1,a);
+
+ double Hb(0.), dHbdT(0.), ddHbddT(0.);
+ if (q1->_ipKinematic != NULL){
+ Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
+ dHbdT = Hb * q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
+ ddHbddT = Hb * q1->_ipKinematic->getDDRDDT();
+ }
+
+ //a.print("a init");
+
+ static STensor3 devN; // dev part of yield normal
+ double trN = 0.; // trace part of yield normal
+ static STensor3 N; // yield normal
+ STensorOperation::zero(devN);
+ STensorOperation::zero(N);
+
+ double f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
+
+ double DfDGamma = 0.;
+ double dfdDgamma = 0.;
+
+ // NEW
+ // Change u and v inside iterations
+ double u = 1.;
+ double v = 1.;
+ // NEW
+
+ double A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
+
+
+ double dDgammaDGamma = 0.;
+ double Dgamma = 0.; // eqplastic strain
+
+// NEW
+ double Gt = 0.; // Ge + kk*Hb/2.;
+ double Kt = 0.; // Ke + kk*Hb/3.;
+ static fullVector<double> m(2);
+ getChabocheCoeffs(m,0,q1);
+ double Cx(0.), Px(0.);
+ getPlasticPotential(ptilde, PhiEq, Px);
+ Cx = m(0)*(q0->_epspbarre + Dgamma) + m(1)*Px - 1/Hb * dHbdT + 1;
+ Gt = Ge + pow(kk,2) * Hb/(2*Cx);
+ Kt = Ke + pow(kk,2) * Hb/(3*Cx);
+// NEW
+
+ if (q1->dissipationIsBlocked()){
+ q1->getRefToDissipationActive() = false;
+ }
+ else{
+ if (f>_tol){
+ q1->getRefToDissipationActive() = true;
+ // plasticity
+ int ite = 0;
+ int maxite = 100; // maximal number of iters
+
+
+ //Msg::Error("plasticity occurs f = %e",f);
+
+ //double f0 = fabs(f);
+
+ while (fabs(f) >_tol or ite <1){
+ double eta(0.),Deta(0.),DetaDT(0.);
+ if (_viscosity != NULL)
+ _viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT);
+ double etaOverDt = eta/this->getTimeStep();
+
+ // NEW
+ // DEFINE THESE INSIDE PLASTIC ITERATIONS -> Because dgamma is available here
+ // need to calculate m(0) and m(1), preferably like yield coeffs (a(i)), they are functions of gamma
+ // make sure that Px is already multiplied to m(1) or something
+ // put these inside the iterations? Because dgamma isnt available here -> may have to put it in IP?
+ // add temp funcs for dHbdT
+
+ // Cx = m(0)*(q0->_epspbarre + Dgamma) + m(1)*Px - 1/Hb * dHbdT + 1;
+ // Gt = Ge + pow(kk,2) * Hb/(2*Cx);
+ // Kt = Ke + pow(kk,2) * Hb/(3*Cx);
+ // NEW
+
+ // NEW
+ double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A);
+ // double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A); // -> the same
+ // NEW
+
+ dDgammaDGamma = kk*(A+Gamma*dAdGamma); // mistake in the paper (VD 2016)
+
+ this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
+ dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*ptilde -Da(0);
+ // NEW
+ double dCxdDgamma = m(0);
+ double Stemp = STensorOperation::doubledot(devPhi,devXn);
+ double He = (PhiEq*6*Gamma*kk*kk*Hb - Stemp*3.*1./PhiEq)/(2*Cx*Cx*u) * dCxdDgamma; // CHANGE THIS DERIVATIVE -> tensor subtracting scalar
+ double Hp = ptilde*(2*_b*Gamma*kk*kk*Hb - trXn)/(3*Cx*Cx*v) * dCxdDgamma;
+ dfdDgamma += a(2)*_n*pow(PhiEq,(_n-1))*He - a(1)*Hp;
+ // NEW
+ if (Gamma>0 and etaOverDt>0)
+ dfdDgamma -= _p*pow(etaOverDt,_p-1.)*Deta/this->getTimeStep()*pow(Gamma,_p); // THIS term is absent in the paper!! WHAT??
+
+ DfDGamma = dfdDgamma*dDgammaDGamma - (_n*a(2)*6.*Gt)*pow(PhiEq,_n)/u + a(1)*ptilde*2.*_b*Kt/v;
+ if (Gamma>0 and etaOverDt>0)
+ DfDGamma -=pow(etaOverDt,_p)*_p*pow(Gamma,_p-1.);
+
+ double dGamma = -f/DfDGamma; // WHAT
+
+ if (Gamma + dGamma <=0.){
+ Gamma /= 2.; // WHAT
+ }
+ else
+ Gamma += dGamma;
+
+ //Msg::Error("Gamma = %e",Gamma);
+
+ u = 1.+6.*Gt*Gamma;
+ v = 1.+2.*_b*Kt*Gamma;
+ // NEW
+ // PhiEq = PhiEqpr/u; // -> original
+ // ptilde = ptildepr/v; // -> original
+ devPhi = (devPhipr + devXn * (1-1/Cx));
+ devPhi *= 1./u;
+ PhiEq = sqrt(1.5*devPhi.dotprod());
+ ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v;
+ // NEW
+
+ A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
+ Dgamma = kk*Gamma*A;
+
+ //Msg::Error("it = %d, u=%e, v=%e, Dgamma=%e",ite,u,v,Dgamma);
+
+ updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
+ hardening(q0,q1,T);
+ getYieldCoefficients(q1,a);
+ //a.print("a update");
+
+ f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
+ double viscoTerm = etaOverDt*Gamma; // I THINK temperature dependency of viscosity should be present in eta, but this relation remains the same
+ if (Gamma>0 and etaOverDt>0) f-= pow(viscoTerm,_p);
+
+ ite++;
+ //if (ite> maxite-5)
+ //Msg::Error("it = %d, DfDGamma = %e error = %e dGamma = %e, Gamma = %e",ite,DfDGamma,f,dGamma,Gamma);
+
+ if (fabs(f) <_tol) break;
+
+ if(ite > maxite){
+ Msg::Error("No convergence for plastic correction in mlawNonLinearTVP nonAssociatedFlow Maxwell iter = %d, f = %e!!",ite,f);
+ P(0,0) = P(1,1) = P(2,2) = sqrt(-1.);
+ return;
+ }
+ };
+
+ q1->_DgammaDt = Dgamma/this->getTimeStep();
+
+ // update effective stress tensor
+ // NEW
+ devPhi = (devPhipr + devXn * (1-1/Cx));
+ devPhi *= 1./u; //NEW
+ ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v; //NEW
+ // NEW
+
+ // update normal
+ devN = devPhi;
+ devN *= 3.;
+ trN = 2.*_b*ptilde;
+ N = devN;
+ N(0,0) += trN/3.;
+ N(1,1) += trN/3.;
+ N(2,2) += trN/3.;
+
+ // estimate exp(GammaN)
+ // static STensor3 expGN;
+ static STensor3 GammaN;
+ GammaN = N;
+ GammaN *= Gamma;
+ STensorOperation::expSTensor3(GammaN,_order,expGN,&dexpAdA);
+
+ // update plastic deformation tensor
+ STensorOperation::multSTensor3(expGN,Fp0,Fp1);
+ // update IP
+ updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma); // Why again???? WHAT??
+
+ // update elastic deformation tensor, corotational stress
+ STensorOperation::inverseSTensor3(Fp1,Fpinv);
+ STensorOperation::multSTensor3(F,Fpinv,Fe);
+ STensorOperation::multSTensor3FirstTranspose(Fe,Fe,Ce);
+ STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCe,&DDlnDDCe);
+ Ee *= 0.5; // Here We have assumed that, Ce commutes with Q (i.e. with M)
+
+ // update A, B
+ ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T); //,false);
+ getTVEdCorKirDT(q0,q1,T0,T); // update dCorKirdT
+ // DKDTsum and DGDTsum = sum of bulk/shear moduli derivatives wrt T
+
+ // backstress
+ static STensor3 DB; // increment
+ DB = (N); // increment
+ DB *= (kk*Hb*Gamma);
+ // NEW
+ DB += q0->_backsig;
+ DB *= 1/Cx;
+ DB -= q0->_backsig;
+ // NEW
+ q1->_backsig += DB; // update
+
+ /*
+ // corotationaal Kirchhoff stress
+ q1->_kirchhoff = devPhi;
+ q1->_kirchhoff += q1->_backsig;
+
+ q1->_kirchhoff(0,0) += (ptilde);
+ q1->_kirchhoff(1,1) += (ptilde);
+ q1->_kirchhoff(2,2) += (ptilde);
+
+ */
+ /*
+ // NEW
+ q1->_ModMandel = devPhi;
+ q1->_ModMandel += q1->_backsig;
+
+ q1->_ModMandel(0,0) += (ptilde);
+ q1->_ModMandel(1,1) += (ptilde);
+ q1->_ModMandel(2,2) += (ptilde);
+
+ // NEW
+
+ }
+ else{
+ q1->getRefToDissipationActive() = false;
+ }
+ }
+
+
+ const STensor3& KS = q1->_kirchhoff;
+ // second Piola Kirchhoff stress
+ // static STensor3 S;
+ // STensorOperation::multSTensor3STensor43(KS,DlnDCe,S);
+
+ // NEW
+ getModifiedMandel(Ce, q0, q1); // update Mandel // this step is probably wrong -> Mandel is already being updated above in the NR loop. Should corKir come from Mandel then?
+ // WHAT ???
+ const STensor3& MS = q1->_ModMandel;
+ // second Piola Kirchhoff stress
+ static STensor3 S, Ceinv;
+ STensorOperation::inverseSTensor3(Ce,Ceinv);
+ STensorOperation::multSTensor3(Ceinv,MS,S);
+ // NEW
+
+ for(int i=0; i<3; i++)
+ for(int j=0; j<3; j++){
+ P(i,j) = 0.;
+ for(int k=0; k<3; k++)
+ for(int l=0; l<3; l++)
+ P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
+ }
+
+
+ // defo energy
+ // q1->getRefToElasticEnergy()=deformationEnergy(*q1);
+ q1->getRefToViscousEnergyPart()=viscousEnergy(*q0,*q1)+q0->getConstRefToViscousEnergyPart();
+ q1->getRefToPlasticEnergy() = q0->plasticEnergy();
+ if (Gamma > 0){
+ // NEW
+ // double dotKSN = dot(KS,N);
+ //q1->getRefToPlasticEnergy() += Gamma*dotKSN;
+ double dotMSN = dot(MS,N);
+ q1->getRefToPlasticEnergy() += Gamma*dotMSN; // = Dp:Me
+ // NEW
+ }
+
+
+ // fluxT
+ double KThConT, DKThConDT; mlawNonLinearTVE::getKTHCon(KThConT,T,&DKThConDT);
+ double J = 1.;
+ STensor3 Finv(0.);
+ if (_thermalEstimationPreviousConfig){ // ADD _thermalEstimationPreviousConfig
+ STensorOperation::inverseSTensor3(F0,Finv);
+ J = STensorOperation::determinantSTensor3(F0);
+ }
+ else{
+ STensorOperation::inverseSTensor3(F,Finv);
+ J = STensorOperation::determinantSTensor3(F);
+ }
+
+ static STensor3 Cinv;
+ STensorOperation::multSTensor3SecondTranspose(Finv,Finv,Cinv);
+ STensorOperation::multSTensor3SVector3(Cinv,gradT,fluxT);
+ fluxT *= (-KThConT*J);
+
+ // thermalEnergy -> assuming stiff = false
+ double CpT;
+ if (stiff){
+ double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
+ getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
+ CpT = -T*DDpsiTVMdTT;
+ }
+ else{
+ getCp(CpT,T);
+ }
+
+ q1->_thermalEnergy = CpT*T;
+
+ // thermalSource
+ if (this->getTimeStep() > 0.){
+ thermalSource = -CpT*(T-T0)/this->getTimeStep();
+ }
+ else
+ thermalSource = 0.;
+
+ // NEW - mechSourceTVP
+ double& Wm = q1->getRefToMechSrcTVP();
+ double& Wm_TVE = q1->getRefToMechSrcTVE();
+
+ // Have to do this to update DcorKirDT in IP, and mechSourceTVE and derivatives
+ mechanicalSource = 0.;
+ mlawNonLinearTVE::getMechSourceTVE(q0,q1,T0,T,DKDTsum,DGDTsum,_I4,&Wm_TVE);
+ // STensor3 temp101 = q1->_DcorKirDT;
+ mechanicalSource += Wm_TVE;
+
+ STensor3 P_TVE(0.);
+ SVector3 fluxT_TVE(0.);
+ double thermalSource_TVE(0.);
+ double mechanicalSource_TVE(0.);
+ /* mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(q0->getConstRefToFe(), q1->getConstRefToFe(), P_TVE, q0, q1, Tangent, dFedF, dFedT,
+ T0, T, gradT0, gradT, fluxT_TVE, dPdT, dfluxTdgradT, dfluxTdT, dfluxTdF,
+ thermalSource_TVE, dthermalSourcedT, dthermalSourcedF,
+ mechanicalSource, dmechanicalSourcedT, dmechanicalSourceF,
+ stiff);
+
+
+ getMechSourceTVP(q0,q1,T0,T,_I,_I4,_I,&Wm);
+ mechanicalSource += Wm;
+
+ // freeEnergy and elastic energy
+ double& psiTVM = q1->getRefToFreeEnergyTVM();
+ getFreeEnergyTVM(q0,q1,T0,T,&psiTVM,NULL);
+
+ q1->_elasticEnergy = mlawNonLinearTVE::freeEnergyMechanical(*q1,T); // deformationEnergy(*q1,T);
+
+ // I didnt make this
+ if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
+ q1->getRefToIrreversibleEnergy() = q1->defoEnergy();
+ }
+ else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
+ (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
+ q1->getRefToIrreversibleEnergy() = q1->plasticEnergy();
+ }
+ else{
+ q1->getRefToIrreversibleEnergy() = 0.;
+ }
+
+ */
\ No newline at end of file
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.h b/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
index 060ed56a4..eaef7b8ec 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
@@ -22,8 +22,6 @@ class mlawNonLinearTVP : public mlawNonLinearTVE{
scalarFunction* _temFunc_Sy0; // temperature dependence of initial yield stress
scalarFunction* _temFunc_H; // temperature dependence of hardening stress
- scalarFunction* _temFunc_visc; // temperature dependence of viscosity law
- scalarFunction* _temFunc_Hb; // temperature dependence of kinematic hardening modulus
std::vector<double> _HR; // isotropic hardening moduli for R
@@ -32,7 +30,7 @@ class mlawNonLinearTVP : public mlawNonLinearTVE{
// NEW
virtual void hardening(const IPNonLinearTVP* q0, IPNonLinearTVP* q1, const double& T) const;
virtual void getModifiedMandel(const STensor3& C, const IPNonLinearTVP *q0_, IPNonLinearTVP *q1) const;
- virtual void checkCommutavity(STensor3& commuteCheck, const STensor3& Ce, const STensor3& S, const IPNonLinearTVP *q1) const;
+ virtual void checkCommutavity(const STensor3& C, const STensor3& S, const IPNonLinearTVP *q1) const;
virtual void getChabocheCoeffs(fullVector<double>& coeffs, const double& opt, const IPNonLinearTVP *q1) const;
virtual void getPlasticPotential(const double& phiP, const double& phiE, double& Px) const;
virtual void getYieldCoefficientTempDers(const IPNonLinearTVP *q1, const double& T, fullVector<double>& DcoeffsDT, fullVector<double>& DDcoeffsDTT) const;
@@ -45,7 +43,6 @@ class mlawNonLinearTVP : public mlawNonLinearTVE{
virtual void getFreeEnergyTVM(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double& T0, const double& T,
double *psiTVM = NULL, double *DpsiTVMdT = NULL, double *DDpsiTVMdTT = NULL) const;
-
// NEW
virtual void predictorCorrector_TVP_nonAssociatedFlow(const STensor3& F0, const STensor3& F, const IPNonLinearTVP *q0_, IPNonLinearTVP *q1,
@@ -162,19 +159,6 @@ class mlawNonLinearTVP : public mlawNonLinearTVE{
void setTemperatureFunction_Hardening(const scalarFunction& Tfunc){
if (_temFunc_H != NULL) delete _temFunc_H;
_temFunc_H = Tfunc.clone();
- // _compression->setTemperatureFunction(Tfunc);
- // _traction->setTemperatureFunction(Tfunc); // - this is difficult, because of multiple parameters
- }
- void setTemperatureFunction_KinematicHardening(const scalarFunction& Tfunc){
- if (_temFunc_Hb != NULL) delete _temFunc_Hb;
- _temFunc_Hb = Tfunc.clone();
- // _kinematic->setTemperatureFunction_K(Tfunc);
- // - this is not a pure virtual parent class, on the other hand the viscosity/kine harde type isnt set here but in the .py file.
- }
- void setTemperatureFunction_Viscosity(const scalarFunction& Tfunc){
- if (_temFunc_visc != NULL) delete _temFunc_visc;
- _temFunc_visc = Tfunc.clone();
- // _viscosity->setTemperatureFunction(Tfunc);
}
void setTaylorQuineyFactor(const double f){_TaylorQuineyFactor = f;};
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index be9f036c6..9b814e932 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -8575,6 +8575,7 @@ void NonLinearTVPDG3DMaterialLaw::setTaylorQuineyFactor(const double f){
void NonLinearTVPDG3DMaterialLaw::setTemperatureFunction_InitialYieldStress(const scalarFunction& Tfunc){
_viscoLaw.setTemperatureFunction_InitialYieldStress(Tfunc);
};
+/*
void NonLinearTVPDG3DMaterialLaw::setTemperatureFunction_Hardening(const scalarFunction& Tfunc){
_viscoLaw.setTemperatureFunction_Hardening(Tfunc);
};
@@ -8584,7 +8585,7 @@ void NonLinearTVPDG3DMaterialLaw::setTemperatureFunction_KinematicHardening(cons
void NonLinearTVPDG3DMaterialLaw::setTemperatureFunction_Viscosity(const scalarFunction& Tfunc){
_viscoLaw.setTemperatureFunction_Viscosity(Tfunc);
};
-
+*/
// Added extra argument in the NonLinearTVEDG3DIPVariable contructor - Tinitial (added from LinearThermoMech)
void NonLinearTVPDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 234321640..a70797517 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -3000,9 +3000,9 @@ class NonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{ // public Material
void setTaylorQuineyFactor(const double f);
void setTemperatureFunction_InitialYieldStress(const scalarFunction& Tfunc);
- void setTemperatureFunction_Hardening(const scalarFunction& Tfunc);
+ /*void setTemperatureFunction_Hardening(const scalarFunction& Tfunc);
void setTemperatureFunction_KinematicHardening(const scalarFunction& Tfunc);
- void setTemperatureFunction_Viscosity(const scalarFunction& Tfunc);
+ void setTemperatureFunction_Viscosity(const scalarFunction& Tfunc);*/
#ifndef SWIG
NonLinearTVPDG3DMaterialLaw(const NonLinearTVPDG3DMaterialLaw &source);
--
GitLab