diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index 44ae251b26e11b5160339404eace358e6a74b83e..80f399b0c7e9d2b2af610117c4986ebcf499b279 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -25,6 +25,7 @@ mlawNonLinearTVP::mlawNonLinearTVP(const int num,const double E,const double nu,
_temFunc_H = new constantScalarFunction(1.);
_temFunc_Hb = new constantScalarFunction(1.);
_temFunc_visc = new constantScalarFunction(1.);
+
};
// , mlawPowerYieldHyper(src)
@@ -41,6 +42,7 @@ mlawNonLinearTVP::mlawNonLinearTVP(const mlawNonLinearTVP& src): mlawNonLinearTV
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();
@@ -48,7 +50,6 @@ mlawNonLinearTVP::mlawNonLinearTVP(const mlawNonLinearTVP& src): mlawNonLinearTV
_temFunc_visc = NULL; // temperature dependence of viscosity
if (src._temFunc_visc != NULL)
_temFunc_visc = src._temFunc_visc->clone();
-
};
mlawNonLinearTVP& mlawNonLinearTVP::operator=(const materialLaw& source){
@@ -79,13 +80,11 @@ 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;
+ if (_temFunc_visc != NULL) delete _temFunc_visc; _temFunc_visc= NULL;
};
// CHECK IF IPSTATE NEEDS the same definition as in mlawPowerYieldHyper --- WHAT??
-// NEW
-
void mlawNonLinearTVP::setIsotropicHardeningCoefficients(const double HR1, const double HR2, const double HR3){
_HR.clear();
@@ -149,10 +148,9 @@ void mlawNonLinearTVP::checkCommutavity(STensor3& commuteCheck,const STensor3& C
void mlawNonLinearTVP::getChabocheCoeffs(fullVector<double>& coeffs, const double& opt, const IPNonLinearTVP *q1) const{
double p = q1->_epspbarre; // eq.plastic strain
-
coeffs.resize(2);
- coeffs(0) = 1; //DEBUGGING
+ coeffs(0) = 1; // 1; //DEBUGGING
coeffs(1) = 0;
// add dependency on p later
@@ -365,13 +363,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;
+ const double& dGammaDT = 0.; // q1->_DGammaDT; = 0. for debugging //DEBUGGING
const STensor3& DgammaDF = q1->_DgammaDF;
const STensor3& DGammaDF = q1->_DGammaDF;
const double& gamma0 = q0->_epspbarre;
@@ -390,8 +388,10 @@ void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNo
STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
STensorOperation::decomposeDevTr(q0->_DModMandelDT,dDevMeDT,dTrMeDT_0);
- 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;
+ 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
// get Dgamma
double Dgamma = gamma1 - gamma0;
@@ -438,14 +438,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);
- dRdF[1](i,j) = -Da(0)*DgammaDF(i,j)*sigc - a(0)*Hc;
+ 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);
}
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;
+ 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);
}
}
@@ -465,11 +465,11 @@ void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNo
if (delT>1e-10){
double DphiPDT = dTrMeDT- dTrXdT;
- double DDphiPDTT = 2*DphiPDT/delT - 2*(trMe-trMe_0+trX1_0-trX1)/pow(delT,2);
+ double DDphiPDTT = 0.; // (DphiPDT - (trMe-trMe_0+trX1_0-trX1)/delT)/delT; // 2*DphiPDT/delT - 2*(trMe-trMe_0+trX1_0-trX1)/pow(delT,2); // DEBUGGING
ddRdTT->at(0) -= a(1) * DDphiPDTT * sigc;
if (Gamma>0 and etaOverDt>0){
- double DDGammaDTT = 2*dGammaDT/delT - 2*(Gamma-Gamma_0)/pow(delT,2);
+ double DDGammaDTT = 0.; // 2*dGammaDT/delT - 2*(Gamma-Gamma_0)/pow(delT,2); // DEBUGGING
ddRdTT->at(2) -= _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(eta*DDGammaDTT)/this->getTimeStep() * sigc;
}
}
@@ -492,7 +492,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;
@@ -507,7 +507,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;
+ const STensor3& dXdT(0.); // = q1->_DbackSigDT; // Just for debugging = 0. //DEBUGGING
const STensor3& dXdT_0 = q0->_DbackSigDT;
const STensor43& dXdF_0 = q0->_DbackSigDF;
const STensor43& dXdF = q1->_DbackSigDF;
@@ -541,7 +541,8 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
STensorOperation::inverseSTensor3(Fp1,Fpinv);
Fp_dot = Fp1;
Fp_dot -= Fp0;
- STensorOperation::multSTensor3(Dp,Fp_dot,Fpinv);
+ Fp_dot *= 1/this->getTimeStep();
+ STensorOperation::multSTensor3(Fp_dot,Fpinv,Dp);
// get alpha
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
@@ -560,7 +561,7 @@ 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;
@@ -573,7 +574,7 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
dRdT += dIsoHardForcedT[i];
ddRdTT += ddIsoHardForcedTT[i];
dRdF += dIsoHardForcedF[i];
- ddRdFdT += ddIsoHardForcedFdT[i];
+ ddRdFdT += 0.; // ddIsoHardForcedFdT[i]; //DEBUGGING
}
R -= T*dRdT;
@@ -592,6 +593,7 @@ 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++)
@@ -600,7 +602,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,j,k,s)*Fpinv(s,l);
+ dFpinvdFp(i,j,k,l) -= Fpinv(i,m)*_I4(m,s,k,l)*Fpinv(s,j);
}
for (int i=0; i<3; i++)
@@ -609,9 +611,17 @@ 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,j,k,m)*Fpinv(m,l) + Fp_dot(i,m)*dFpinvdFp(m,j,k,l);
+ 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);
}
+ 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);
@@ -621,10 +631,9 @@ 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 k=0; k<3; k++)
- for (int l=0; l<3; l++)
- dDpdT(i,j) += dDpdFp(i,j,k,l)*dFpdT(k,l);
- }
+ 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
+ }}
// dXdF is available
// const STensor43 dXdF = q1->_DbackSigDF;
@@ -638,10 +647,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) += 2*dXdT(i,j)/(T-T0) - 2*(X1(i,j)-X0(i,j))/pow((T-T0),2);
+ // ddXdTT(i,j) += 2*dXdT(i,j)/(T-T0) - 2*(X1(i,j)-X0(i,j))/pow((T-T0),2); //DEBUGGING
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);
+ // ddXdTF(i,j,k,l) += (dXdF(i,j,k,l)-dXdF_0(i,j,k,l))/(T-T0); //DEBUGGING
}
}
}
@@ -675,8 +684,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(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();
+ 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();
}
}
}
@@ -1022,10 +1031,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
double &dmechanicalSourcedT,
STensor3 &dmechanicalSourceF) const{
- // Compute elastic predictor first
+ // compute elastic predictor
STensor3& Fp1 = q1->_Fp;
const STensor3& Fp0 = q0->_Fp;
-
+
// Initialise predictor values
Fp1 = Fp0; // plastic deformation tensor
q1->_epspbarre = q0->_epspbarre; // plastic equivalent strain
@@ -1036,7 +1045,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
q1->_DbackSigDT = q0->_DbackSigDT; // DbackSigDT
q1->_DgammaDt = 0.;
- // Get hardening stuff
+ // Get hardening parameters
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);
@@ -1044,13 +1053,12 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
double Hb(0.), dHb(0.), dHbdT(0.), ddHbddT(0.);
if (q1->_ipKinematic != NULL){
Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
- dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter derivative (dHb/dgamma ??)
+ dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter derivative (dHb/dgamma)
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");
+ //a.print("a init");
// Get Cepr, Ceinvpr
static STensor3 Fpinv, Ce, Fepr;
@@ -1061,7 +1069,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static STensor3 Cepr, Ceinvpr;
Cepr = Ce;
STensorOperation::inverseSTensor3(Cepr,Ceinvpr);
-
+
static STensor3 invFp0; // plastic predictor
invFp0= Fpinv;
STensor3& Fe = q1->_Fe;
@@ -1093,8 +1101,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// Get Xn
static STensor3 devXn;
- static double trXn;
+ static double trXn, pXn, eXn;
STensorOperation::decomposeDevTr(q0->_backsig,devXn,trXn); // needed for chaboche model
+ pXn = 1/3*trXn; // 1st invariant of Xn
+ eXn = sqrt(1.5*devXn.dotprod()); // J2 invariant of Xn
// Initialise Phipr
static STensor3 PhiPr;
@@ -1105,7 +1115,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static double ptildepr,ptilde;
STensorOperation::decomposeDevTr(PhiPr,devPhipr,ptildepr);
ptildepr/= 3.;
-
+
// Initialise Normal
static STensor3 devN; // dev part of yield normal
double trN = 0.; // trace part of yield normal
@@ -1123,15 +1133,19 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
double dDgammaDGamma = 0.;
double Dgamma = 0.; // eqplastic strain
static fullVector<double> m(2);
- double Px(0.);
getChabocheCoeffs(m,0,q1);
double DfDGamma = 0.;
double dfdDgamma = 0.;
+ double dAdDgamma = 0.;
+ double dAdGamma = 0.;
+ double dfdGamma = 0.;
double u = 1.;
double v = 1.;
// Initialise Cx, Gt, Kt
+ double Px(0.);
+ getPlasticPotential(pXn, eXn, Px); // CHANGE -> this is wrong, change to backstress instead of phi
double Cx = m(0)*Dgamma + m(1)*Px + 1;
if (Hb>0.){Cx -= 1/Hb * dHbdT * (T-T0);}
double Gt = Ge + pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
@@ -1144,13 +1158,11 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
double PhiEqpr2 = 1.5*devPhi.dotprod();
double PhiEqpr = sqrt(PhiEqpr2); // -> second invariant
double PhiEq = PhiEqpr; // current effective deviator
- getPlasticPotential(ptildepr, PhiEqpr, Px); // CHANGE -> this is wrong, change to backstress instead of phi
// Get f and A
double f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
double A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
-
// NR Loop
if (q1->dissipationIsBlocked()){
q1->getRefToDissipationActive() = false;
@@ -1158,7 +1170,6 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
else{
if (f>_tol){
q1->getRefToDissipationActive() = true;
-
// plasticity
int ite = 0;
int maxite = 100; // maximal number of iters
@@ -1168,14 +1179,15 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
while (fabs(f) >_tol or ite <1){
// Viscosity
- double eta(0.), Deta(0.), DetaDT(0.);
+ double eta(0.),Deta(0.),DetaDT(0.);
if (_viscosity != NULL)
_viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT);
double etaOverDt = eta/this->getTimeStep();
- double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u + 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A);
+ // dAdGamma and dDgammaDGamma
+ dAdGamma = -(72.*Gt*PhiEq*PhiEq/u + 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A);
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);
@@ -1185,66 +1197,76 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// 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 //CHECK
// double Hp = ptilde*(2*_b*Gamma*kk*kk*Hb - trXn)/(3*Cx*Cx*v) * dCxdDgamma; // CHECK
-
- static STensor3 DdevPhidDgamma, term2;
+
+ static STensor3 DdevPhidDgamma, term2;
STensorOperation::zero(DdevPhidDgamma);
// 1st term of He
DdevPhidDgamma = 3*Gamma*pow(kk,1)*Hb*devPhipr; // pow(kk,2) CHANGE
DdevPhidDgamma -= devXn*(1 + 6*Gamma*(Ge + pow(kk,1)*Hb/2)); // pow(kk,2) CHANGE
DdevPhidDgamma *= 1/(Cx*Cx*u*u) * dCxdDgamma;
-
+
// 2nd term of He
- term2 = 3*Gamma*pow(kk,1)*Cx*devPhi; // pow(kk,2) CHANGE
- term2 *= 1/(Cx*Cx*u) * dHb;
-
+ term2 = (3*Gamma*pow(kk,1)*Cx*devPhi);
+ term2 *= 1/(Cx*Cx*u)*dHb; // pow(kk,2) CHANGE
+
DdevPhidDgamma -= term2;
double Stemp = STensorOperation::doubledot(devPhi,DdevPhidDgamma);
double He = 1.5*Stemp/PhiEq;
-
+
double Hp = ( 2*_b*Gamma*pow(kk,1)*Hb*ptildepr - trXn*(1 + 2*_b*Gamma*(Ke + pow(kk,1)*Hb/3)) )*1/(3*Cx*Cx*v*v)*dCxdDgamma;
Hp -= (2*_b*Gamma*pow(kk,1)*Cx*ptilde)*1/(3*Cx*Cx*v) * dHb;
-
+
dfdDgamma += a(2)*_n*pow(PhiEq,(_n-1))*He - a(1)*Hp;
-
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.));
+ DfDGamma -=pow(etaOverDt,_p)*_p*pow(Gamma,_p-1.);
+
+ // OLD NR
+ /*
+ double dGamma = -f/DfDGamma; // NR
- // double dGamma = -f/DfDGamma; // NR
+ if (Gamma + dGamma <=0.){
+ Gamma /= 2.; // NR
+ }
+ else
+ Gamma += dGamma;*/
-
- double dfdGamma = - (_n*a(2)*6.*Gt)*pow(PhiEq,_n)/u + a(1)*ptilde*2.*_b*Kt/v;
+ // NEW NR
+ dfdGamma = - (_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 dAdDgamma = (12*PhiEq*He + 8/3*_b*_b*ptilde*Hp)/(2*A);
- double dgdDgamma = 1 - kk*Gamma*dAdDgamma;
+
+ dAdDgamma = (12*PhiEq*He + 8/3*_b*_b*ptilde*Hp)/(2*A);
+ double dgdDgamma = 1 - kk*Gamma*dAdDgamma;
double dgdGamma = - dDgammaDGamma;
double g = Dgamma - kk*Gamma*A;
-
+
double dGamma = (-f + dfdDgamma*g/dgdDgamma)/(-dfdDgamma*dgdGamma/dgdDgamma + dfdGamma); // NR
double dDgamma = -(g+dgdGamma*dGamma)/dgdDgamma; // NR
-
+
if (Gamma + dGamma <=0.){
Gamma /= 2.; // NR
}
else
Gamma += dGamma;
-
+
if (Dgamma + dDgamma <=0.){
Dgamma /= 2.;
- }
+ }
else
- Dgamma += dDgamma;
-
+ Dgamma += dDgamma; // Debugging
+
//Msg::Error("Gamma = %e",Gamma);
//Msg::Error("Dgamma = %e",Dgamma);
+
+ // UPDATE FOR NEXT ITERATION
+
// Update gamma and all Hardening moduli
updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
hardening(q0,q1,T);
@@ -1256,31 +1278,37 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
ddHbddT = Hb * q1->_ipKinematic->getDDRDDT();
}
//a.print("a update");
-
+
+ // Update Normal and Backstress to get Px
+ // Can't use this because we Px needs backstress (n+1 step) which comes from normal -> which needs Phi -> which needs Cx -> which needs Px
+ // Have to put another rate equation in the NR for backstress to get it independently -> dX = kk*Hb*dalpha or something like that
+ // TBD -> also changes Phi definition and derivatives
+ // instead Px -> only as a function of backstress (nth step) is used.
+
// Update Cx, u, v
Cx = m(0)*Dgamma + m(1)*Px + 1;
if (Hb>0.){Cx -= 1/Hb * dHbdT * (T-T0);}
Gt = Ge + pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
Kt = Ke + pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
u = 1.+6.*Gt*Gamma;
- v = 1.+2.*_b*Kt*Gamma;
-
+ v = 1.+2.*_b*Kt*Gamma;
+
// Update Phi
devPhi = devPhipr + devXn*(1-1/Cx);
- ptilde = ptildepr + 1/3*trXn*(1-1/Cx);
- devPhi *= 1./u;
+ ptilde = ptildepr + 1/3*trXn*(1-1/Cx);
+ devPhi *= 1./u;
PhiEq = sqrt(1.5*devPhi.dotprod());
- ptilde *= 1./v;
-
+ ptilde /= v;
+
// Update A
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);
-
+
// Update f
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
+ double viscoTerm = etaOverDt*Gamma;
if (Gamma>0 and etaOverDt>0) f-= pow(viscoTerm,_p);
ite++;
@@ -1300,10 +1328,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// Correct Phi (effective stress tensor)
devPhi = devPhipr + devXn*(1-1/Cx);
- ptilde = ptildepr + 1/3*trXn*(1-1/Cx);
- devPhi *= 1./u;
+ ptilde = ptildepr + 1/3*trXn*(1-1/Cx);
+ devPhi *= 1./u;
PhiEq = sqrt(1.5*devPhi.dotprod());
- ptilde *= 1./v;
+ ptilde /= v;
// Correct normal
devN = devPhi;
@@ -1324,9 +1352,9 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// Correct plastic deformation tensor
STensorOperation::multSTensor3(expGN,Fp0,Fp1);
// Correct IP gamma
- updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
+ updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
- // Correct elastic deformation tensor (Ee), corotational stress
+ // Correct elastic deformation tensor, corotational stress
STensorOperation::inverseSTensor3(Fp1,Fpinv);
STensorOperation::multSTensor3(F,Fpinv,Fe);
STensorOperation::multSTensor3FirstTranspose(Fe,Fe,Ce);
@@ -1340,8 +1368,8 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// Correct backstress
static STensor3 DB; // increment
- DB = (N);
- DB *= (pow(kk,1)*Hb*Gamma); // pow(kk,2) CHANGE
+ DB = (N); // increment
+ DB *= (kk*Hb*Gamma); // pow(kk,2) CHANGE
DB += q0->_backsig;
DB *= 1/Cx;
DB -= q0->_backsig;
@@ -1361,22 +1389,18 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
}
}
+
const STensor3& KS = q1->_kirchhoff;
- // second Piola Kirchhoff stress
- // static STensor3 S;
- // STensorOperation::multSTensor3STensor43(KS,DlnDCe,S);
- // NEW
- getModifiedMandel(Ce, q0, q1); // update Mandel
+ getModifiedMandel(Ce, q0, q1); // update Mandel
const STensor3& MS = q1->_ModMandel;
// second Piola Kirchhoff stress
static STensor3 S, Ceinv;
STensorOperation::inverseSTensor3(Ce,Ceinv);
STensorOperation::multSTensor3(Ceinv,q1->_ModMandel,S);
- // NEW
-
- // first PK
+
+ // first PK
for(int i=0; i<3; i++)
for(int j=0; j<3; j++){
P(i,j) = 0.;
@@ -1391,12 +1415,8 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
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
}
@@ -1421,7 +1441,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// 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);
+ // STensorOperation::zero(dFpdT);
// I didnt make this
if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
@@ -1575,13 +1595,15 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
//7. get DGDCepr, DgammaDCepr
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
- DGDCepr(i,j) = (-dfDCepr(i,j)-dfdDgamma*kk*Gamma*dAdCepr(i,j))/DfDGamma;
+ // DGDCepr(i,j) = (-dfDCepr(i,j)-dfdDgamma*kk*Gamma*dAdCepr(i,j))/DfDGamma;
+ DGDCepr(i,j) = (-dfDCepr(i,j)-dfdDgamma*kk*Gamma*dAdCepr(i,j)/(1-kk*Gamma*dAdDgamma)) / (dfdDgamma*kk*(A+Gamma*dAdGamma)/(1.-kk*Gamma*dAdDgamma) + dfdGamma);
}
}
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
- DgammaDCepr(i,j) = kk*Gamma*dAdCepr(i,j)+ kk*dDgammaDGamma*DGDCepr(i,j);
+ // DgammaDCepr(i,j) = kk*Gamma*dAdCepr(i,j)+ kk*dDgammaDGamma*DGDCepr(i,j);
+ DgammaDCepr(i,j) = (kk*Gamma*dAdCepr(i,j)+ (1.*kk*A+kk*Gamma*dAdDgamma)*DGDCepr(i,j))/(1.-kk*Gamma*dAdDgamma);
}
}
@@ -1603,16 +1625,15 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
temp7 = dCxdCepr;
temp7 *= (1/3*1/(pow(Cx,2))*trXn/v);
- temp9 = dCxdCepr;
- temp9 *= (2/3*_b*Gamma*pow(kk,1)*Hb/pow(Cx,2)); // pow(kk,2) CHANGE
+ temp9 = 2/3*_b*Gamma*pow(kk,1)*(-Hb/pow(Cx,2)*dCxdCepr + 1/Cx*dHb*DgammaDCepr); // pow(kk,2) CHANGE
temp8 = DGDCepr;
temp8 *= (2*_b*(Ke+pow(kk,1)*Hb/(3*Cx))); // pow(kk,2) CHANGE
- temp8 -= temp9;
+ temp8 += temp9;
temp8 *= (ptildepr + 1/3*trXn*(1-1/Cx));
temp8 *= 1/pow(v,2);
- DphiPDCepr += temp7; // CHECK sign - corrected
+ DphiPDCepr += temp7;
DphiPDCepr -= temp8;
// DdevphiDCepr
@@ -1624,8 +1645,8 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
STensorOperation::zero(temp13);
STensorOperation::zero(temp14);
- temp12 = devXn; // This must be temp12?
- temp12 *= 1/pow(Cx,2); // This must be temp12?
+ temp12 = devXn;
+ temp12 *= 1/pow(Cx,2);
STensorOperation::prod(temp12, dCxdCepr, 1., temp10);
temp10 *= 1/u;
@@ -1633,14 +1654,14 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
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);
- 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
+ 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) + 1/Cx*dHb*DgammaDCepr(i,j)); // pow(kk,2) CHANGE
}
}
STensorOperation::prod(temp13, temp14, 1., temp11);
temp11 *= 1/pow(u,2);
- DdevphiDCepr += temp10; // CHECK sign - corrected
+ DdevphiDCepr += temp10;
DdevphiDCepr -= temp11;
//9. get DtrNDCepr DdevNdCepr
@@ -1720,7 +1741,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
dTrXdCepr(i,j) = pow(kk,1)*Hb*temp15_tr(i,j)/Cx - trX*dCxdCepr(i,j)/Cx; // pow(kk,2) CHANGE
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
- dXdCepr(i,j,k,l) = pow(kk,1)*Hb*temp15(i,j,k,l)/Cx - q1->_backsig(i,j)*dCxdCepr(k,l)/Cx; // CHECK, MIGHT HAVE TO DIFFERENTIATE Hb wrt gamma
+ dXdCepr(i,j,k,l) = pow(kk,1)*(Hb*temp15(i,j,k,l) + dHb*Gamma*DgammaDCepr(i,j)*N(k,l))/Cx - q1->_backsig(i,j)*dCxdCepr(k,l)/Cx;
// pow(kk,2) CHANGE
}
@@ -1795,7 +1816,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// Done below -> need dCedF
STensor43& dXdF = q1->getRefToDbackStressdF();
- static STensor3 dTrXdF, DphiPDF;
+ static STensor3 dTrXdF;
STensorOperation::multSTensor43(dXdCepr,CeprToF,dXdF);
// DphiPDCepr - use this for the pressure term in R
@@ -1805,7 +1826,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
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);
}
}
@@ -1941,7 +1962,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// 1. get dMeprDT from dKeprDT -> DcorKirprDT from TVE
static STensor3 dMeprDT, dDevMeprDT, DcorKirDT;
- double dpMeprDT(0.);
+ static double dpMeprDT;
const STensor3& dKeprDT = q1->getConstRefToDcorKirDT();
STensorOperation::zero(dMeprDT);
@@ -1966,8 +1987,8 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// 3. get dCxdT, dGxdT, dKxdT, dXdT
static double dCxdT, dGxdT, dKxdT;
-
- dCxdT = 1/(pow(Hb,2))*pow(dHbdT,2)*(T-T0) - 1/Hb*ddHbddT*(T-T0) - 1/Hb*dHbdT;
+ dCxdT = 0.;
+ if (Hb>0.){dCxdT = 1/(pow(Hb,2))*pow(dHbdT,2)*(T-T0) - 1/Hb*ddHbddT*(T-T0) - 1/Hb*dHbdT;}
dGxdT = DGDTsum + pow(kk,1)/(2*Cx)*dHbdT - (pow(kk,1)*Hb)/(2*pow(Cx,2))*dCxdT; // pow(kk,2) CHANGE
dKxdT = DKDTsum + pow(kk,1)/(3*Cx)*dHbdT - (pow(kk,1)*Hb)/(3*pow(Cx,2))*dCxdT; // pow(kk,2) CHANGE
@@ -1985,7 +2006,7 @@ 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 // CHECK the sign of the 2nd term -corrected
+ 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;
@@ -1993,7 +2014,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
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); // CHECK the sign of the 2nd term -corrected
+ dDevPhiDT(i,j) = (dDevPhiprDT(i,j) + devXn(i,j)/pow(Cx,2)*dCxdT)/u - temp19(i,j);
}
}
@@ -2022,14 +2043,23 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
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));
}
// 7. get dGammaDT, DgammaDT
- dGammaDT = (-dfdT-dfdDgamma*kk*Gamma*dAdT)/DfDGamma;
- DgammaDT = kk*Gamma*dAdT + kk*dDgammaDGamma*dGammaDT;
+ // dGammaDT = (-dfdT-dfdDgamma*kk*Gamma*dAdT)/DfDGamma;
+ dGammaDT = (-dfdT-dfdDgamma*kk*Gamma*dAdT/(1.-kk*Gamma*dAdDgamma))/(dfdDgamma*(1.*kk*A+kk*Gamma*dAdGamma)/(1.-kk*Gamma*dAdDgamma) + dfdGamma);
+ // DgammaDT = kk*Gamma*dAdT + kk*dDgammaDGamma*dGammaDT;
+ DgammaDT = (kk*Gamma*dAdT + (1*kk*A+kk*Gamma*dAdDgamma)*dGammaDT)/(1.-kk*Gamma*dAdDgamma);
+
+ // 8.0 update dCxdT, dKxdT, dGxdT
+ static double dHxdT;
+ dHxdT = dHbdT + dHb*DgammaDT;
+ if (Hb>0.){dCxdT = 1/(pow(Hb,2))*pow(dHxdT,2)*(T-T0) - 1/Hb*ddHbddT*(T-T0) - 1/Hb*dHxdT;} // -> assuming this is OK
+ dGxdT = DGDTsum + pow(kk,1)/(2*Cx)*dHxdT - (pow(kk,1)*Hb)/(2*pow(Cx,2))*dCxdT; // pow(kk,2) CHANGE
+ dKxdT = DKDTsum + pow(kk,1)/(3*Cx)*dHxdT - (pow(kk,1)*Hb)/(3*pow(Cx,2))*dCxdT; // pow(kk,2) CHANGE
// 8. update DdevphiDT and DphiPDT
static double dPxdT, temp18_2;
@@ -2060,7 +2090,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
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++)
@@ -2104,16 +2134,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
@@ -2188,14 +2214,13 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
}
}
-
- }
+ }
if(stiff){
-
+
// TVP - flux, thermalEnergy, thermalSource
if (!_thermalEstimationPreviousConfig){
-
+
static STensor3 DJDF;
static STensor43 DCinvDF;
for (int i=0; i<3; i++){
@@ -2230,9 +2255,9 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
}
// thermSrc and MechSrc
-
- // thermalEnergy
- double CpT;
+
+ // thermalEnergy
+ double CpT;
if (stiff){
double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
@@ -2242,9 +2267,9 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
else{
getCp(CpT,T);
}
-
- q1->_thermalEnergy = CpT*T;
-
+
+ q1->_thermalEnergy = CpT*T;
+
// thermalSource
if (this->getTimeStep() > 0.){
thermalSource = -CpT*(T-T0)/this->getTimeStep();
@@ -2255,22 +2280,22 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// 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; // DEBUGGING
+ mechanicalSource += Wm_TVE; // DEBUGGING
getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,&Wm);
- // mechanicalSource += Wm; // DEBUGGING
-
+ mechanicalSource += Wm; // DEBUGGING
+
// 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);
@@ -2282,7 +2307,7 @@ if (stiff){
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();
@@ -2296,10 +2321,10 @@ if (stiff){
dthermalSourcedT = 0.;
STensorOperation::zero(dthermalSourcedF);
}
-
-
+
+
// mechSourceTVP derivatives
-
+
// conversion tensor for mechSourceTVE Terms
static STensor43 DCeDFe, DEeDFe;
for (int i=0; i<3; i++)
@@ -2311,26 +2336,26 @@ if (stiff){
DCeDFe(i,j,k,l) += 2*F(p,i)*dFedF(p,j,k,l);
}
}
-
+
STensorOperation::multSTensor43(DlnDCe,DCeDFe,DEeDFe);
DEeDFe *= 0.5;
-
+
double& dWmdT_TVE = q1->getRefTodMechSrcTVEdT();
STensor3& dWmdF_TVE = q1->getRefTodMechSrcTVEdF();
mlawNonLinearTVE::getMechSourceTVE(q0,q1,T0,T,DKDTsum,DGDTsum,DEeDFe,NULL,&dWmdF_TVE,&dWmdT_TVE);
-
+
dmechanicalSourcedT = 0.;
STensorOperation::zero(dmechanicalSourceF);
- // dmechanicalSourceF += dWmdF_TVE; // DEBUGGING
- // dmechanicalSourcedT += dWmdT_TVE; // DEBUGGING
-
+ dmechanicalSourceF += dWmdF_TVE; // DEBUGGING
+ dmechanicalSourcedT += dWmdT_TVE; // DEBUGGING
+
double& dWmdT = q1->getRefTodMechSrcTVPdT();
STensor3& dWmdF = q1->getRefTodMechSrcTVPdF();
- getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,NULL,&dWmdF,&dWmdT);
-
- // dmechanicalSourceF += dWmdF; // DEBUGGING
- // dmechanicalSourcedT += dWmdT; // DEBUGGING
+ getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,NULL,&dWmdF,&dWmdT);
+
+ dmechanicalSourceF += dWmdF; // DEBUGGING
+ dmechanicalSourcedT += dWmdT; // DEBUGGING
}
};
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index 0042ecfe845da139f2bed7cba6936183d8c0a827..016efa98814b6fc3917c75e598fb7a560cd2ddc0 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -5180,4 +5180,4 @@ class NonLinearTVPDG3DIPVariable : public ThermoMechanicsDG3DIPVariableBase{ /
// NonLinearTVP Law Interface =================================================================== END
-#endif //DG3DIPVARIABLE
+#endif
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 7ddfc5fff55a96e7349214256723bce6e42889d0..12c4f2edd10cd17efc58e81981dfdbc3ab8cbb1d 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -245,8 +245,11 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
int numCurlDof = ipvcur->getNumConstitutiveCurlVariable();
_stressLaw->stress(ipvcur,ipvprev,false,checkfrac,dTangent);
+
#if 1
_stressLaw->stress(ipvcur,ipvprev,true,checkfrac,dTangent);
+ ipvcur->getConstRefToDeformationGradient().print("F Analytique"); // FLE
+ ipvcur->getConstRefToFirstPiolaKirchhoffStress().print("P Analytique"); // FLE
ipvcur->getConstRefToTangentModuli().print("dPdE Analytique");
for (int k=0; k< numNonLocalVars; k++)
{
@@ -287,6 +290,10 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
}
for (int i=0; i< numExtraDof; i++)
{
+ Msg::Info("Field Analytique: %e",ipvcur->getConstRefToField(i)); // FLE
+ ipvcur->getConstRefToFlux()[i].print("Flux Analytique"); // FLE
+ Msg::Info("ThermSource Analytique: %e",ipvcur->getConstRefToFieldSource()(i)); // FLE
+ Msg::Info("MechSource Analytique: %e",ipvcur->getConstRefToMechanicalSource()(i)); // FLE
ipvcur->getConstRefTodFluxdF()[i].print("d Flux dF Analytique");
ipvcur->getConstRefTodFieldSourcedF()[i].print("d Field source dF Analytique");
ipvcur->getConstRefTodMechanicalSourcedF()[i].print("d Mechanical source dF Analytique");
@@ -608,6 +615,8 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
}
#if 0
+ ipvcur->getConstRefToDeformationGradient().print("F Numerique"); // FLE
+ ipvcur->getConstRefToFirstPiolaKirchhoffStress().print("P Numerique"); // FLE
ipvcur->getConstRefToTangentModuli().print("dPdE Numerique");
for (int k=0; k< numNonLocalVars; k++)
{
@@ -646,6 +655,10 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
}
for (int i=0; i< numExtraDof; i++)
{
+ Msg::Info("Field Numerique: %e",ipvcur->getConstRefToField(i)); // FLE
+ ipvcur->getConstRefToFlux()[i].print("Flux Numerique"); // FLE
+ Msg::Info("ThermSource Numerique: %e",ipvcur->getConstRefToFieldSource()(i)); // FLE
+ Msg::Info("MechSource Numerique: %e",ipvcur->getConstRefToMechanicalSource()(i)); // FLE
ipvcur->getConstRefTodFluxdF()[i].print("d Flux dF Numerique");
ipvcur->getConstRefTodFieldSourcedF()[i].print("d Field source dF Numerique");
ipvcur->getConstRefTodMechanicalSourcedF()[i].print("d Mechanical source dF Numerique");