diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
index 222d0f41b6e3c7d4b18e36deb51131a6b618d8b9..52203c66306d22634176a89e0f90859419b45565 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
@@ -858,6 +858,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
double DfDGamma = 0.;
double dfdDgamma = 0.;
+ double dfdGamma = 0.; //FLE
double u = 1.;
double v = 1.;
@@ -887,17 +888,41 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
if (_viscosity != NULL)
_viscosity->get(q1->_epspbarre,eta,Deta);
double etaOverDt = eta/this->getTimeStep();
+ double dPhiPdGamma = -2*Kt*_b*ptilde/v;
+ double dPhiEdGamma = -6*Gt*PhiEq/u;
+ double dAdGamma_TVP = (12*PhiEq*dPhiEdGamma + 8*_b*_b*ptilde*dPhiPdGamma/3)/(2.*A);
double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A);
dDgammaDGamma = kk*(A+Gamma*dAdGamma);
this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
+
+
+ // DEBUG // FLE
+ double a0 = a(0);
+ double a1 = a(1);
+ double a2 = a(2);
+ double da0 = Da(0);
+ double da1 = Da(1);
+ double da2 = Da(2);
+
dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*ptilde -Da(0);
if (Gamma>0 and etaOverDt>0)
dfdDgamma -= _p*pow(etaOverDt,_p-1.)*Deta/this->getTimeStep()*pow(Gamma,_p);
+
+ dfdGamma = _n*a(2)*pow(PhiEq,(_n-1))*dPhiEdGamma - a(1)*dPhiPdGamma;
+ if (Gamma>0 and etaOverDt>0)
+ dfdGamma -= pow(etaOverDt,_p)*_p*pow(Gamma,(_p-1.));
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.);
+
+ Msg::Error("dDgammaDGamma = %e, iter = %d",dDgammaDGamma, ite);
+ Msg::Error("dfdDgamma = %e, iter = %d",dfdDgamma, ite);
+ Msg::Error("dfdGamma = %e, iter = %d",dfdGamma, ite);
+ Msg::Error("DfDGamma = %e, iter = %d",DfDGamma, ite);
+ Msg::Error("Gamma = %e, iter = %d",Gamma, ite);
+ Msg::Error("Dgamma = %e, iter = %d",Dgamma, ite);
double dGamma = -f/DfDGamma;
@@ -1064,6 +1089,14 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
if (Gamma >0){
// plastic
static STensor3 dAdCepr, dfDCepr;
+
+ // DEBUG FLE
+ double a0 = a(0);
+ double a1 = a(1);
+ double a2 = a(2);
+ double da0 = Da(0);
+ double da1 = Da(1);
+ double da2 = Da(2);
double fact = 1.5*a(2)*_n*pow(PhiEq,_n-2.)/(u*u);
for (int i=0; i<3; i++){
@@ -1088,7 +1121,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
- DgamaDCepr(i,j) = kk*Gamma*dAdCepr(i,j)+ kk*dDgammaDGamma*DGDCepr(i,j);
+ DgamaDCepr(i,j) = kk*Gamma*dAdCepr(i,j)+ dDgammaDGamma*DGDCepr(i,j);
}
}
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
index 4d054cc05cc5c091ef6654e967bfe44d72c96e35..136d84f42ddcff181fb71a5da0e184419fd29351 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
@@ -1853,6 +1853,7 @@ void mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
+ // Msg::Info("dfluxTdF->J: Jacobian= %e",J);
dfluxTdF(i,j,k) += (DJDF(j,k)*fluxT(i)/J);
for (int l=0; l<3; l++){
dfluxTdF(i,j,k) -= (J*DCinvDF(i,l,j,k)*gradT(l)*KThConT);
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index 6ea7b7adda34c0681bdbf6365ef26289eb74b75f..aceaccba6b03109330304b429ff7994f7c5f9908 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -219,21 +219,21 @@ void mlawNonLinearTVP::getYieldCoefficientTempDers(const IPNonLinearTVP *q1, con
DDcoeffsDTT.resize(3);
double m = sigt/sigc;
double DmDT = (dsigtdT*sigc- sigt*dsigcdT)/(sigc*sigc);
- double DDmDTT = (ddsigtdTT*sigc - sigt*ddsigcdTT)/(sigc*sigc) - 2*(dsigtdT*sigc- sigt*dsigcdT)/(sigc*sigc*sigc)*dsigcdT;
+ 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 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);
+ 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;
+ 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(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;
+ 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(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;
};
@@ -1229,6 +1229,15 @@ void mlawNonLinearTVP::plasticCorrector_TVP(const STensor3& F0, const STensor3&
dDgammaDGamma = kk*(A+Gamma*dAdGamma); // mistake in the paper (VD 2016)
this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
+
+ // DEBUG // FLE
+ double a0 = a(0);
+ double a1 = a(1);
+ double a2 = a(2);
+ double da0 = Da(0);
+ double da1 = Da(1);
+ double da2 = Da(2);
+
dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*ptilde -Da(0);
dfdDgamma += a(2)*_n*pow(PhiEq,(_n-1))*He - a(1)*Hp;
@@ -1501,7 +1510,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static STensor3 devXn, devX;
static double trXn, pXn, eXn, trX;
STensorOperation::decomposeDevTr(q0->_backsig,devXn,trXn); // needed for chaboche model
- pXn = 1/3*trXn; // 1st invariant of Xn
+ pXn = 1./3.*trXn; // 1st invariant of Xn
eXn = sqrt(1.5*devXn.dotprod()); // J2 invariant of Xn
// Initialise Phipr
@@ -1541,15 +1550,15 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
double v = 1.;
// Initialise Cx, Gt, Kt
- double Cxdev = 1. + 3*m(1)*pow(kk,1)*Hb; // pow(kk,2) DEBUG
- double Cxtr = 1. + 2*_b/3 *m(1)*pow(kk,1)*Hb; // pow(kk,2) DEBUG
+ double Cxdev = 1. + 3.*m(1)*pow(kk,1.)*Hb; // pow(kk,2) DEBUG
+ double Cxtr = 1. + 2.*_b/3. *m(1)*pow(kk,1.)*Hb; // pow(kk,2) DEBUG
double dCxdevdDgamma(0.), dCxtrdDgamma(0.), dudDgamma(0.), dvdDgamma(0.), dudGamma(0.), dvdGamma(0.);
if (Hb>0.){
- Cxdev -= m(0)*dHb*Dgamma/Hb + 1/Hb * dHbdT * (T-T0);
- Cxtr -= m(0)*dHb*Dgamma/Hb + 1/Hb * dHbdT * (T-T0);
+ Cxdev -= m(0)*dHb*Dgamma/Hb + 1./Hb * dHbdT * (T-T0);
+ Cxtr -= m(0)*dHb*Dgamma/Hb + 1./Hb * dHbdT * (T-T0);
}
- double Gt = Ge + pow(kk,1) * Hb/(2*Cxdev); // pow(kk,2) DEBUG
- double Kt = Ke + pow(kk,1) * Hb/(3*Cxtr); // pow(kk,2) DEBUG
+ double Gt = Ge + pow(kk,1.) * Hb/(2.*Cxdev); // pow(kk,2) DEBUG
+ double Kt = Ke + pow(kk,1.) * Hb/(3.*Cxtr); // pow(kk,2) DEBUG
// Initialise ptilde and devPhi (phi)
ptilde = ptildepr; // current effective pressure -> first invariant
@@ -1588,31 +1597,31 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// Get dCxdevdDgamma, dCxtrdDgamma, dudDgamma dvdDgamma, dudGamma, dvdGamma
dCxdevdDgamma = 0.; dCxtrdDgamma = 0.;
if (Hb>0.){
- dCxdevdDgamma = -m(0)*1/Hb*dHb - m(0)*Dgamma*(-1/pow(Hb,2)*pow(dHb,2) + 1/Hb*ddHb) - (T-T0)*(-1/pow(Hb,2)*dHb*dHbdT + 1/Hb*ddHbdgammadT);
+ dCxdevdDgamma = -m(0)*1./Hb*dHb - m(0)*Dgamma*(-1./pow(Hb,2.)*pow(dHb,2.) + 1./Hb*ddHb) - (T-T0)*(-1./pow(Hb,2.)*dHb*dHbdT + 1./Hb*ddHbdgammadT);
dCxtrdDgamma = dCxdevdDgamma;
}
- dCxdevdDgamma += 3*m(1)*pow(kk,1)*dHb; // pow(kk,2) DEBUG
- dCxtrdDgamma += (2*_b/3)*m(1)*pow(kk,1)*dHb; // pow(kk,2) DEBUG
+ dCxdevdDgamma += 3.*m(1)*pow(kk,1.)*dHb; // pow(kk,2.) DEBUG
+ dCxtrdDgamma += (2.*_b/3.)*m(1)*pow(kk,1.)*dHb; // pow(kk,2.) DEBUG
- dudDgamma = 6*Gamma*( pow(kk,1)/(2*Cxdev)*dHb - (pow(kk,1)*Hb)/(2*pow(Cxdev,2))*dCxdevdDgamma ); // pow(kk,2) DEBUG
- dvdDgamma = 2*_b*Gamma*( pow(kk,1)/(3*Cxtr)*dHb - (pow(kk,1)*Hb)/(3*pow(Cxtr,2))*dCxtrdDgamma ); // pow(kk,2) DEBUG
+ dudDgamma = 6.*Gamma*( pow(kk,1.)/(2.*Cxdev)*dHb - (pow(kk,1.)*Hb)/(2.*pow(Cxdev,2.))*dCxdevdDgamma ); // pow(kk,2.) DEBUG
+ dvdDgamma = 2.*_b*Gamma*( pow(kk,1.)/(3.*Cxtr)*dHb - (pow(kk,1.)*Hb)/(3.*pow(Cxtr,2.))*dCxtrdDgamma ); // pow(kk,2.) DEBUG
- dudGamma = 6*(Ge + pow(kk,1)*Hb/(2*Cxdev)); // pow(kk,2) DEBUG
- dvdGamma = 2*_b*(Ke + pow(kk,1)*Hb/(3*Cxtr)); // pow(kk,2) DEBUG
+ dudGamma = 6.*(Ge + pow(kk,1.)*Hb/(2.*Cxdev)); // pow(kk,2.) DEBUG
+ dvdGamma = 2.*_b*(Ke + pow(kk,1.)*Hb/(3.*Cxtr)); // pow(kk,2.) DEBUG
// Get Hp = dPhiPdDgamma, dPhiPdGamma
STensor3 dHp_RHS = _I;
- dHp_RHS *= (1/3*trXn/(pow(Cxtr,2))*dCxtrdDgamma - ptilde*dvdDgamma);
+ dHp_RHS *= (1./3.*trXn/(pow(Cxtr,2.))*dCxtrdDgamma - ptilde*dvdDgamma);
/*
- double Hp = STensorOperation::doubledot(Dho4_v_inv,dHp_RHS); // 1/v*(1/3*trXn/(pow(Cxtr,2))*dCxtrdDgamma - ptilde*dvdDgamma);
- Hp /= 3;*/
- double Hp = 1/v*(1/3*trXn/(pow(Cxtr,2))*dCxtrdDgamma - ptilde*dvdDgamma);
+ double Hp = STensorOperation::doubledot(Dho4_v_inv,dHp_RHS); // 1/v*(1./3.*trXn/(pow(Cxtr,2.))*dCxtrdDgamma - ptilde*dvdDgamma);
+ Hp /= 3.;*/
+ double Hp = (1./3.*trXn/(pow(Cxtr,2.))*dCxtrdDgamma - ptilde*dvdDgamma)/v;
STensor3 dPhiPdGamma_RHS, temp; STensorOperation::multSTensor3STensor43(_I,Dho3,dPhiPdGamma_RHS);
temp = _I;
temp *= ptilde*dvdGamma;
- dPhiPdGamma_RHS *= 2*_b/3*ptilde;
+ dPhiPdGamma_RHS *= 2.*_b/3.*ptilde;
dPhiPdGamma_RHS -= temp;
/*
@@ -1626,13 +1635,13 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
STensor3 DdevPhidGamma, DdevPhidGamma_RHS;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
- DdevPhidDgamma_RHS(i,j) = devXn(i,j)/(pow(Cxdev,2))*dCxdevdDgamma - dudDgamma*devPhi(i,j);
+ DdevPhidDgamma_RHS(i,j) = devXn(i,j)/(pow(Cxdev,2.))*dCxdevdDgamma - dudDgamma*devPhi(i,j);
DdevPhidGamma_RHS(i,j) = -dudGamma*devPhi(i,j);
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
- DdevPhidDgamma_RHS(i,j) += 2*_b/3*Gamma*Dho3(i,j,k,l)*_I(k,l)*Hp; //ADD
+ DdevPhidDgamma_RHS(i,j) += (2.*_b/3.)*Gamma*Dho3(i,j,k,l)*_I(k,l)*Hp; //ADD
- DdevPhidGamma_RHS(i,j) += 2*_b/3*Gamma*Dho3(i,j,k,l)*_I(k,l)*dPhiPdGamma; //ADD
+ DdevPhidGamma_RHS(i,j) += (2.*_b/3.)*Gamma*Dho3(i,j,k,l)*_I(k,l)*dPhiPdGamma; //ADD
DdevPhidGamma_RHS(i,j) += Dho3(i,j,k,l)*N(k,l); // 3*Dho3(i,j,k,l)*devPhi(k,l);
}
@@ -1736,25 +1745,44 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
}*/
// dAdGamma and dDgammaDGamma
- dAdDgamma = (12*PhiEq*He + 8/3*_b*_b*ptilde*Hp)/(2.*A);
- dAdGamma = (12*PhiEq*dPhiEdGamma + 8/3*_b*_b*ptilde*dPhiPdGamma)/(2.*A);
+
+ dAdDgamma = (12.*PhiEq*He + 8.*_b*_b*ptilde*Hp/3.)/(2.*A);
+ dAdGamma = (12.*PhiEq*dPhiEdGamma + 8.*_b*_b*ptilde*dPhiPdGamma/3.)/(2.*A);
dDgammaDGamma = kk*(A+Gamma*dAdGamma); // mistake in the paper (VD 2016)
this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
+
+ // DEBUG FLE
+ double a0 = a(0);
+ double a1 = a(1);
+ double a2 = a(2);
+ double da0 = Da(0);
+ double da1 = Da(1);
+ double da2 = Da(2);
+
dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*ptilde -Da(0);
- dfdDgamma += a(2)*_n*pow(PhiEq,(_n-1))*He - a(1)*Hp;
+ 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!! WHY??
- dfdGamma = _n*a(2)*pow(PhiEq,(_n-1))*dPhiEdGamma - a(1)*dPhiPdGamma;
+ dfdGamma = _n*a(2)*pow(PhiEq,(_n-1.))*dPhiEdGamma - a(1)*dPhiPdGamma;
if (Gamma>0 and etaOverDt>0)
dfdGamma -= pow(etaOverDt,_p)*_p*pow(Gamma,(_p-1.));
DfDGamma = dfdDgamma*dDgammaDGamma + dfdGamma;
+
+ /*
+ Msg::Error("dDgammaDGamma = %e, iter = %d",dDgammaDGamma, ite);
+ Msg::Error("dfdDgamma = %e, iter = %d",dfdDgamma, ite);
+ Msg::Error("dfdGamma = %e, iter = %d",dfdGamma, ite);
+ Msg::Error("DfDGamma = %e, iter = %d",DfDGamma, ite);
+ Msg::Error("Gamma = %e, iter = %d",Gamma, ite);
+ Msg::Error("Dgamma = %e, iter = %d",Dgamma, ite);*/
+
- double dgdDgamma = 1 - kk*Gamma*dAdDgamma;
+ double dgdDgamma = 1. - kk*Gamma*dAdDgamma;
double dgdGamma = - dDgammaDGamma;
double g = Dgamma - kk*Gamma*A;
@@ -1792,14 +1820,14 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
//a.print("a update");
// Update Cx, u, v
- Cxdev = 1. + 3*m(1)*pow(kk,1)*Hb; // pow(kk,2) DEBUG
- Cxtr = 1. + 2*_b/3 *m(1)*pow(kk,1)*Hb; // pow(kk,2) DEBUG
+ Cxdev = 1. + 3.*m(1)*pow(kk,1.)*Hb; // pow(kk,2) DEBUG
+ Cxtr = 1. + 2.*_b/3. *m(1)*pow(kk,1.)*Hb; // pow(kk,2) DEBUG
if (Hb>0.){
- Cxdev -= m(0)*dHb*Dgamma/Hb - 1/Hb * dHbdT * (T-T0);
- Cxtr -= m(0)*dHb*Dgamma/Hb - 1/Hb * dHbdT * (T-T0);
+ Cxdev -= m(0)*dHb*Dgamma/Hb - 1./Hb * dHbdT * (T-T0);
+ Cxtr -= m(0)*dHb*Dgamma/Hb - 1./Hb * dHbdT * (T-T0);
}
- Gt = Ge + pow(kk,1) * Hb/(2*Cxdev); // pow(kk,2) DEBUG
- Kt = Ke + pow(kk,1) * Hb/(3*Cxtr); // pow(kk,2) DEBUG
+ Gt = Ge + pow(kk,1.) * Hb/(2.*Cxdev); // pow(kk,2) DEBUG
+ Kt = Ke + pow(kk,1.) * Hb/(3.*Cxtr); // pow(kk,2) DEBUG
u = 1.+6.*Gt*Gamma;
v = 1.+2.*_b*Kt*Gamma;
@@ -1828,7 +1856,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// Update f
f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
double viscoTerm = etaOverDt*Gamma;
- if (Gamma>0 and etaOverDt>0) f-= pow(viscoTerm,_p);
+ if (Gamma>0. and etaOverDt>0.) f-= pow(viscoTerm,_p);
ite++;
//if (ite> maxite-5)
@@ -1846,14 +1874,14 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
q1->_DgammaDt = Dgamma/this->getTimeStep();
// Correct Cx, u, v
- Cxdev = 1. + 3*m(1)*pow(kk,1)*Hb; // pow(kk,2) DEBUG
- Cxtr = 1. + 2*_b/3 *m(1)*pow(kk,1)*Hb; // pow(kk,2) DEBUG
+ Cxdev = 1. + 3.*m(1)*pow(kk,1.)*Hb; // pow(kk,2.) DEBUG
+ Cxtr = 1. + 2.*_b/3. *m(1)*pow(kk,1.)*Hb; // pow(kk,2.) DEBUG
if (Hb>0.){
- Cxdev -= m(0)*dHb*Dgamma/Hb - 1/Hb * dHbdT * (T-T0);
- Cxtr -= m(0)*dHb*Dgamma/Hb - 1/Hb * dHbdT * (T-T0);
+ Cxdev -= m(0)*dHb*Dgamma/Hb - 1./Hb * dHbdT * (T-T0);
+ Cxtr -= m(0)*dHb*Dgamma/Hb - 1./Hb * dHbdT * (T-T0);
}
- Gt = Ge + pow(kk,1) * Hb/(2*Cxdev); // pow(kk,2) DEBUG
- Kt = Ke + pow(kk,1) * Hb/(3*Cxtr); // pow(kk,2) DEBUG
+ Gt = Ge + pow(kk,1.) * Hb/(2.*Cxdev); // pow(kk,2.) DEBUG
+ Kt = Ke + pow(kk,1.) * Hb/(3.*Cxtr); // pow(kk,2.) DEBUG
u = 1.+6.*Gt*Gamma;
v = 1.+2.*_b*Kt*Gamma;
@@ -1897,9 +1925,9 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// DKDTsum and DGDTsum = sum of bulk/shear moduli derivatives wrt T
// Correct backstress;
- devX = (pow(kk,1)*Hb*Gamma*devN + devXn); // pow(kk,2) DEBUG
- devX *= 1/Cxdev;
- trX = (pow(kk,1)*Hb*Gamma*trN + trXn)*1/Cxtr; // pow(kk,2) DEBUG
+ devX = (pow(kk,1.)*Hb*Gamma*devN + devXn); // pow(kk,2) DEBUG
+ devX *= 1./Cxdev;
+ trX = (pow(kk,1.)*Hb*Gamma*trN + trXn)*1./Cxtr; // pow(kk,2) DEBUG
q1->_backsig = devX;
q1->_backsig(0,0) += trX/3.;
q1->_backsig(1,1) += trX/3.;
@@ -2073,9 +2101,9 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
DdevphiprDCepr = DdevMeprDCepr;
DphiPDCepr = DphiPprDCepr;
- DphiPDCepr *= 1/v;
+ DphiPDCepr *= 1./v;
DdevphiDCepr = DdevphiprDCepr;
- DdevphiDCepr *= 1/u;
+ DdevphiDCepr *= 1./u;
//5. get DphiEprDCepr from DdevphiDCepr
static STensor3 DphiEDCepr, DphiEDdevPhi;
@@ -2083,55 +2111,6 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
DphiEDdevPhi *= 1.5/(PhiEq);
STensorOperation::multSTensor3STensor43(DphiEDdevPhi,DdevphiDCepr,DphiEDCepr);
-
- //0.0 Numerical Predictor Derivatives
- /*
- double pKepr_plus, pMepr_plus, pKepr, pMepr;
- static STensor3 Cepr_plus, Eepr_plus, Kepr_plus, Mepr_plus, devKepr, devMepr, devKepr_plus, devMepr_plus, dpKepr_plus_dCepr_plus, dpMepr_plus_dCepr_plus;
- static STensor43 DlnDCe_plus, dKepr_plus_dCepr_plus, dMepr_plus_dCepr_plus, ddevKepr_plus_dCepr_plus, ddevMepr_plus_dCepr_plus;
- static STensor63 DDlnDDCe_plus;
- double KTsum_plus, GTsum_plus, DKDTsum_plus, DGDTsum_plus;
-
- STensorOperation::decomposeDevTr(Kepr,devKepr,pKepr);
- STensorOperation::decomposeDevTr(Mepr,devMepr,pMepr);
- pKepr = Kepr.trace()/3;
- pMepr = Mepr.trace()/3;
-
- for (int i=0; i<3; i++){
- for (int j=0; j<3; j++){
- Cepr_plus = Cepr;
- Cepr_plus(i,j) += 0.5*_perturbationfactor;
- Cepr_plus(j,i) += 0.5*_perturbationfactor;
-
- STensorOperation::logSTensor3(Cepr_plus,_order,Eepr_plus,&DlnDCe_plus,&DDlnDDCe_plus);
- Eepr_plus *= 0.5;
-
- for (int k=0; k<3; k++){
- for (int l=0; l<3; l++){
-
- // Numerical A, B
- ThermoViscoElasticPredictor(Eepr_plus,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum_plus,GTsum_plus,DKDTsum_plus,DGDTsum_plus,T0,T);
- // getTVEdCorKirDT(q0,q1,T0,T); // update dCorKirdT
- Kepr_plus = q1->_kirchhoff;
- getModifiedMandel(Cepr_plus, q0, q1); // update Mandel
- Mepr_plus = q1->_ModMandel;
- STensorOperation::decomposeDevTr(Kepr_plus,devKepr_plus,pKepr_plus);
- STensorOperation::decomposeDevTr(Mepr_plus,devMepr_plus,pMepr_plus);
- pKepr_plus = Kepr_plus.trace()/3;
- pMepr_plus = Mepr_plus.trace()/3;
-
- dpKepr_plus_dCepr_plus(k,l) = (pKepr_plus - pKepr)/(_perturbationfactor);
- dpMepr_plus_dCepr_plus(k,l) = (pMepr_plus - pMepr)/(_perturbationfactor);
-
- dKepr_plus_dCepr_plus(k,l,i,j) = (Kepr_plus(k,l) - Kepr(k,l))/(_perturbationfactor);
- dMepr_plus_dCepr_plus(k,l,i,j) = (Mepr_plus(k,l) - Mepr(k,l))/(_perturbationfactor);
-
- ddevKepr_plus_dCepr_plus(k,l,i,j) = (devKepr_plus(k,l) - devKepr(k,l))/(_perturbationfactor);
- ddevMepr_plus_dCepr_plus(k,l,i,j) = (devMepr_plus(k,l) - devMepr(k,l))/(_perturbationfactor);
- }
- }
- }
- }*/
// 6. to 11. (inside if loop-> Gamma > 0.)
static STensor3 dAdCepr, dfDCepr;
@@ -2244,7 +2223,14 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
}
}
}*/
-
+ // DEBUG FLE
+ double a0 = a(0);
+ double a1 = a(1);
+ double a2 = a(2);
+ double da0 = Da(0);
+ double da1 = Da(1);
+ double da2 = Da(2);
+
//6. get dAdCepr from DdevphiDCepr
double fact = a(2)*_n*pow(PhiEq,_n-1.);
for (int i=0; i<3; i++){
@@ -2284,7 +2270,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// update DphiPDCepr
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
- DphiPDCepr(i,j) += ( -dvdCepr(i,j)*ptilde + 1/3*trXn/(pow(Cxtr,2))*dCxtrdCepr(i,j) )/v;
+ DphiPDCepr(i,j) += ( -dvdCepr(i,j)*ptilde + 1./3.*trXn/(pow(Cxtr,2.))*dCxtrdCepr(i,j) )/v;
}
// update DdevphiDCepr, DphiEDCepr
@@ -2295,10 +2281,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
- G3(i,j,k,l) = G2(i,j,k,l) + DdevphiprDCepr(i,j,k,l) + devXn(i,j)*dCxdevdCepr(k,l)/(pow(Cxdev,2)) - devPhi(i,j) * dudCepr(k,l);
+ G3(i,j,k,l) = G2(i,j,k,l) + DdevphiprDCepr(i,j,k,l) + devXn(i,j)*dCxdevdCepr(k,l)/(pow(Cxdev,2.)) - devPhi(i,j) * dudCepr(k,l);
for (int p=0; p<3; p++)
for (int q=0; q<3; q++)
- G3(i,j,k,l) += Dho3(i,j,p,q)* ( N(p,q)*DGDCepr(k,l) + 2*_b*Gamma/3*_I(p,q)*DphiPDCepr(k,l) );
+ G3(i,j,k,l) += Dho3(i,j,p,q)* ( N(p,q)*DGDCepr(k,l) + 2.*_b*Gamma/3.*_I(p,q)*DphiPDCepr(k,l) );
}
STensorOperation::multSTensor43(Dho4_u_inv,G3,DdevphiDCepr);
@@ -2307,10 +2293,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
//9. get DtrNDCepr DdevNdCepr
DtrNDCepr = DphiPDCepr;
- DtrNDCepr *= (2*_b);
+ DtrNDCepr *= (2.*_b);
DdevNDCepr = DdevphiDCepr;
- DdevNDCepr *= 3;
+ DdevNDCepr *= 3.;
// 10. get dFpDCepr
@@ -2341,15 +2327,15 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static STensor43 DdevXdCepr;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
- DtrXdCepr(i,j) = pow(kk,1)*Hb*( DGDCepr(i,j)*trN + Gamma*DtrNDCepr(i,j) ) - dCxtrdCepr(i,j)*trX + pow(kk,1)*Gamma*dHb*DgammaDCepr(i,j)*trN; // pow(kk,2) DEBUG
+ DtrXdCepr(i,j) = pow(kk,1.)*Hb*( DGDCepr(i,j)*trN + Gamma*DtrNDCepr(i,j) ) - dCxtrdCepr(i,j)*trX + pow(kk,1.)*Gamma*dHb*DgammaDCepr(i,j)*trN; // pow(kk,2.) DEBUG
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
- DdevXdCepr(i,j,k,l) = pow(kk,1)*Hb*( DGDCepr(i,j)*devN(k,l) + Gamma*DdevNDCepr(i,j,k,l) ) -
- dCxdevdCepr(i,j)*devX(k,l) + pow(kk,1)*Gamma*dHb*DgammaDCepr(i,j)*devN(k,l); // pow(kk,2) DEBUG
+ DdevXdCepr(i,j,k,l) = pow(kk,1.)*Hb*( DGDCepr(i,j)*devN(k,l) + Gamma*DdevNDCepr(i,j,k,l) ) -
+ dCxdevdCepr(i,j)*devX(k,l) + pow(kk,1.)*Gamma*dHb*DgammaDCepr(i,j)*devN(k,l); // pow(kk,2.) DEBUG
}
}
- DdevXdCepr *= 1/Cxdev;
- DtrXdCepr *= 1/Cxtr;
+ DdevXdCepr *= 1./Cxdev;
+ DtrXdCepr *= 1./Cxtr;
STensorOperation::zero(dXdCepr);
dXdCepr = DdevXdCepr;
@@ -2357,7 +2343,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
- dXdCepr(i,j,k,l) += 1/3* _I(i,j)*DtrXdCepr(k,l);
+ dXdCepr(i,j,k,l) += 1./3. * _I(i,j)*DtrXdCepr(k,l);
// 11. update DpDCepr, DdevKDCepr
@@ -2541,7 +2527,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
STensorOperation::decomposeDevTr(dMeprDT,dDevMeprDT,dpMeprDT);
// dpMeprDT *= 1/3; WHAT IS THE PROBLEM HERE? WHY DOES IT GIVE ZERO?
- dpMeprDT = dMeprDT.trace()/3;
+ dpMeprDT = dMeprDT.trace()/3.;
DcorKirDT = dKeprDT; // update later
@@ -2585,7 +2571,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
dPhiPdT /= v;
dDevPhiDT = dDevPhiprDT;
- dDevPhiDT *= 1/u;
+ dDevPhiDT *= 1./u;
/*
// Use Dho for DdevphiDT
@@ -2620,7 +2606,8 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
if (Gamma >0){
// plastic
-
+
+ /*
//0.0 Numerical Derivatives wrt T
double trN_plus, pKepr_plus, pMepr_plus, pKepr, pMepr, trX_plus, ptilde_plus, PhiEq_plus, A_plus, f_plus, Dgamma_plus, Gamma_plus,
@@ -2710,9 +2697,9 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
dAdT_plus = (A_plus - A)/(_perturbationfactor*T0); // ( 4./3.*_b*_b*ptilde_plus*DphiP_plusDT_plus + 6.*PhiEq_plus*DphiE_plusDT_plus )/A_plus; //
dfdT_plus = (f_plus - f)/(_perturbationfactor*T0); // DaDT_plus(2)*pow(PhiEq_plus,_n) + a_plus(2)*_n*pow(PhiEq_plus,_n-1.)*DphiE_plusDT_plus - DaDT_plus(1)*ptilde_plus -a_plus(1)*DphiP_plusDT_plus - DaDT_plus(0); //
- /*if (this->getTimeStep()>0){
- dfdT_plus -= (DetaDT_plus*Gamma_plus/this->getTimeStep())*_p*pow((eta_plus*Gamma_plus/this->getTimeStep()),(_p-1));
- }*/
+ // if (this->getTimeStep()>0){
+ // dfdT_plus -= (DetaDT_plus*Gamma_plus/this->getTimeStep())*_p*pow((eta_plus*Gamma_plus/this->getTimeStep()),(_p-1));
+ // }
DGDT_plus = (Gamma_plus - Gamma)/(_perturbationfactor*T0); // (-dfdT_plus-dfdDgamma_plus*kk*Gamma_plus*dAdT_plus)/DfDGamma_plus; //
DgammaDT_plus = (Dgamma_plus - Dgamma)/(_perturbationfactor*T0); // kk*Gamma_plus*dAdT_plus+ dDgammaDGamma_plus*DGDT_plus;
@@ -2737,16 +2724,23 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
dMe_plus_dT_plus(i,j) = (Me_plus(i,j)-MS(i,j))/(_perturbationfactor*T0);
}
- }
+ }*/
// 6. get DaDT, dAdT and dfdT
fullVector<double> DaDT(3), DDaDTT(3);
getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
+
+ double DaDT0 = DaDT(0);
+ double DaDT1 = DaDT(1);
+ double DaDT2 = DaDT(2);
+ double DDaDTT0 = DDaDTT(0);
+ double DDaDTT1 = DDaDTT(1);
+ double DDaDTT2 = DDaDTT(2);
dAdT = (4./3*_b*_b*ptilde*dPhiPdT + 6.*PhiEq*dPhiEdT)/A;
- 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));
+ dfdT -= (DetaDT*Gamma/this->getTimeStep())*_p*pow((eta*Gamma/this->getTimeStep()),(_p-1.));
}
// 7. get dGammaDT, DgammaDT
@@ -2757,48 +2751,45 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
static double dCxdevdT, dCxtrdT;
dCxdevdT = 0.; dCxtrdT = 0.;
if (Hb>0.){
- dCxdevdT = -m(0)*(-1/pow(Hb,2)*dHb*dHbdT + 1/Hb*ddHbdgammadT)*Dgamma -1/Hb*dHbdT + (1/pow(kk,1)*dHbdT*dHbdT - 1/Hb*ddHbddT)*(T-T0); // pow(kk,2) DEBUG
+ dCxdevdT = -m(0)*(-1./pow(Hb,2.)*dHb*dHbdT + 1./Hb*ddHbdgammadT)*Dgamma -1./Hb*dHbdT + (1./pow(kk,1.)*dHbdT*dHbdT - 1./Hb*ddHbddT)*(T-T0); // pow(kk,2) DEBUG
}
dCxtrdT = dCxdevdT;
- dCxdevdT += 3*m(1)*pow(kk,1)*dHbdT; // pow(kk,2) DEBUG
- dCxtrdT += 2*_b/3*m(1)*pow(kk,1)*dHbdT; // pow(kk,2) DEBUG
+ dCxdevdT += 3.*m(1)*pow(kk,1.)*dHbdT; // pow(kk,2.) DEBUG
+ dCxtrdT += 2.*_b/3.*m(1)*pow(kk,1.)*dHbdT; // pow(kk,2.) DEBUG
if (Hb>0.){
- dCxdevdT -= m(0)*1/Hb*dHb*DgammaDT;
- dCxtrdT -= m(0)*1/Hb*dHb*DgammaDT;
+ dCxdevdT -= m(0)*1./Hb*dHb*DgammaDT;
+ dCxtrdT -= m(0)*1./Hb*dHb*DgammaDT;
}
static double dudT, dvdT;
- // dudT = 6*Gamma*(DGDTsum + pow(kk,1)/(3*Cxdev)*dHbdT - (pow(kk,1)*Hb)/(3*pow(Cxdev,2))*dCxdevdT); // pow(kk,2) DEBUG
- // dvdT = 2*_b*Gamma*(DKDTsum + pow(kk,1)/(3*Cxtr)*dHbdT - (pow(kk,1)*Hb)/(2*pow(Cxdev,2))*dCxtrdT); // pow(kk,2) DEBUG*/
-
- dudT = 6*Gamma*(DGDTsum + pow(kk,1)/(3*Cxdev)*dHbdT - (pow(kk,1)*Hb)/(3*pow(Cxdev,2))*dCxdevdT); // pow(kk,2) DEBUG
- dvdT = 2*_b*Gamma*(DKDTsum + pow(kk,1)/(3*Cxtr)*dHbdT - (pow(kk,1)*Hb)/(2*pow(Cxdev,2))*dCxtrdT); // pow(kk,2) DEBUG
+ dudT = 6.*Gamma*(DGDTsum + pow(kk,1.)/(2.*Cxdev)*dHbdT - (pow(kk,1.)*Hb)/(2.*pow(Cxdev,2))*dCxdevdT); // pow(kk,2.) DEBUG
+ dvdT = 2.*_b*Gamma*(DKDTsum + pow(kk,1.)/(3.*Cxtr)*dHbdT - (pow(kk,1.)*Hb)/(3.*pow(Cxtr,2))*dCxtrdT); // pow(kk,2.) DEBUG
- dudT += 6*dGammaDT*(Ge+pow(kk,1)*Hb/(2*Cxdev)); // pow(kk,2) DEBUG
- dvdT += 2*_b*dGammaDT*(Ke+pow(kk,1)*Hb/(3*Cxtr)); // pow(kk,2) DEBUG
+ dudT += 6.*dGammaDT*(Ge+pow(kk,1)*Hb/(2.*Cxdev)); // pow(kk,2) DEBUG
+ dvdT += 2.*_b*dGammaDT*(Ke+pow(kk,1)*Hb/(3.*Cxtr)); // pow(kk,2) DEBUG
// 8. update DdevphiDT and DphiPDT
- dPhiPdT = (dPhipprDT - dvdT*ptilde + 1/3*trXn/(pow(Cxtr,2))*dCxtrdT)/v;
+ dPhiPdT = (dPhipprDT - dvdT*ptilde + 1./3.*trXn/(pow(Cxtr,2.))*dCxtrdT)/v;
static STensor3 G2T, G3T;
getG2TTensor(Cepr,expGN,dKeprDT,G2T);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
- G3T(i,j) = G2T(i,j) + dDevPhiprDT(i,j) + devXn(i,j)/pow(Cxdev,2)*dCxdevdT - dudT*devPhi(i,j);
+ G3T(i,j) = G2T(i,j) + dDevPhiprDT(i,j) + devXn(i,j)/pow(Cxdev,2.)*dCxdevdT - dudT*devPhi(i,j);
for (int p=0; p<3; p++)
for (int q=0; q<3; q++)
- G3T(i,j) += Dho3(i,j,p,q)*(dGammaDT*N(p,q) + 2*_b*Gamma/3*dPhiPdT*_I(p,q));
+ G3T(i,j) += Dho3(i,j,p,q)*(dGammaDT*N(p,q) + 2.*_b*Gamma/3.*dPhiPdT*_I(p,q));
}
STensorOperation::multSTensor3STensor43(G3T,Dho4_u_inv,dDevPhiDT);
// 9. get DtrNdT DdevNdT
DtrNdT = dPhiPdT;
- DtrNdT *= (2*_b);
+ DtrNdT *= (2.*_b);
DdevNdT = dDevPhiDT;
- DdevNdT *= 3;
+ DdevNdT *= 3.;
// 10. get dFpdT
static STensor3 dGammaNdT;
@@ -2830,19 +2821,19 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// 10.1 dXdT - backstress temperature derivative -> correction - CHECK!!!
static double DtrXdT;
static STensor3 DdevXdT;
- DtrXdT = pow(kk,1)*Hb*(dGammaDT*trN + Gamma*DtrNdT) - dCxtrdT*trX + pow(kk,1)*Gamma*(dHbdT + dHb*DgammaDT)*trN; // pow(kk,2) DEBUG
+ DtrXdT = pow(kk,1.)*Hb*(dGammaDT*trN + Gamma*DtrNdT) - dCxtrdT*trX + pow(kk,1.)*Gamma*(dHbdT + dHb*DgammaDT)*trN; // pow(kk,2.) DEBUG
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
- DdevXdT(i,j) = pow(kk,1)*Hb*( dGammaDT*devN(i,j) + Gamma*DdevNdT(i,j) ) - dCxdevdT*devX(i,j) + pow(kk,1)*Gamma*(dHbdT + dHb*DgammaDT)*devN(i,j); // pow(kk,2) DEBUG
+ DdevXdT(i,j) = pow(kk,1.)*Hb*( dGammaDT*devN(i,j) + Gamma*DdevNdT(i,j) ) - dCxdevdT*devX(i,j) + pow(kk,1.)*Gamma*(dHbdT + dHb*DgammaDT)*devN(i,j); // pow(kk,2.) DEBUG
}
- DdevXdT *= 1/Cxdev;
- DtrXdT *= 1/Cxtr;
+ DdevXdT *= 1./Cxdev;
+ DtrXdT *= 1./Cxtr;
STensorOperation::zero(dXdT);
dXdT = DdevXdT;
- dXdT(0,0) = DtrXdT/3;
- dXdT(1,1) = DtrXdT/3;
- dXdT(2,2) = DtrXdT/3;
+ dXdT(0,0) = DtrXdT/3.;
+ dXdT(1,1) = DtrXdT/3.;
+ dXdT(2,2) = DtrXdT/3.;
// 11. DcorKirDT
for (int i=0; i<3; i++){
@@ -3050,7 +3041,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
for (int l=0; l<3; l++){
DCeDFe(i,j,k,l) = 0.;
for (int p=0; p<3; p++){
- DCeDFe(i,j,k,l) += 2*F(p,i)*dFedF(p,j,k,l);
+ DCeDFe(i,j,k,l) += 2.*F(p,i)*dFedF(p,j,k,l);
}
}
@@ -3252,11 +3243,11 @@ void mlawNonLinearTVP::getDho3(const double& u, const double& v, const double& G
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
Dho4(i,j,k,l) = 0.;
- Dho4_u(i,j,k,l) += u*_I4(i,j,k,l) - 3*Gamma*Dho3(i,j,k,l);
- Dho4_v(i,j) -= 2*_b/3*Gamma*Dho3(i,j,k,l)*_I(k,l);
+ Dho4_u(i,j,k,l) += u*_I4(i,j,k,l) - 3.*Gamma*Dho3(i,j,k,l);
+ Dho4_v(i,j) -= 2.*_b/3.*Gamma*Dho3(i,j,k,l)*_I(k,l);
for (int p=0; p<3; p++)
for (int q=0; q<3; q++){
- Dho4(i,j,k,l) += (u*_I4(i,j,p,q) - 3*Gamma*Dho3(i,j,p,q))*_Idev(p,q,k,l) + (v*_I4(i,j,p,q) - 2*_b/3*Gamma*Dho3(i,j,p,q))*(_I4(p,q,k,l) -_Idev(p,q,k,l));
+ Dho4(i,j,k,l) += (u*_I4(i,j,p,q) - 3.*Gamma*Dho3(i,j,p,q))*_Idev(p,q,k,l) + (v*_I4(i,j,p,q) - 2.*_b/3.*Gamma*Dho3(i,j,p,q))*(_I4(p,q,k,l) -_Idev(p,q,k,l));
}
}
}
@@ -3307,7 +3298,7 @@ void mlawNonLinearTVP::getIterated_DPhi(const double& u, const double& v, const
STensorOperation::zero(J_constant);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
- J_constant(i,j) += Phipr(i,j) + 1/3*trXn*(1-1/Cxtr)*_I(i,j) + devXn(i,j)*(1-1/Cxdev);
+ J_constant(i,j) += Phipr(i,j) + 1./3.*trXn*(1.-1./Cxtr)*_I(i,j) + devXn(i,j)*(1.-1./Cxdev);
J(i,j) += devPhi(i,j)*u + ptilde*_I(i,j)*v - D(i,j);
}
}
@@ -3383,7 +3374,7 @@ void mlawNonLinearTVP::getIterated_DPhi(const double& u, const double& v, const
STensorOperation::zero(N);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
- N(i,j) = 3*devPhi(i,j) + 2*_b/3*ptilde*_I(i,j);
+ N(i,j) = 3.*devPhi(i,j) + 2.*_b/3.*ptilde*_I(i,j);
// update H = exp(GN) = A in dexpAdA; Hinv
GammaN = N;
@@ -3446,7 +3437,7 @@ void mlawNonLinearTVP::getIterated_DPhi(const double& u, const double& v, const
if(ite > maxite){
Msg::Error("No convergence for iterated Phi mlawNonLinearTVP nonAssociatedFlow Maxwell iter = %d, J_tol = %e!!",ite,J_tol);
- // Phi(0,0) = Phi(1,1) = Phi(2,2) = sqrt(-1.);
+ Phi(0,0) = Phi(1,1) = Phi(2,2) = sqrt(-1.);
return;
}
} // while
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 70a0e5340cbf23be25331486518cc419dfb81a7d..ccd3d00f2226dee21ea87e73f8c566f18baa44e2 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -288,19 +288,20 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
Msg::Info("dEMFieldSource dp Analytique: %e",ipvcur->getConstRefTodEMFieldSourcedNonLocalVariable()(j,i));
}
}
- /*
+
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");
+ // 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");
ipvcur->getConstRefTodPdField()[i].print("dP d Field Analytique");
- ipvcur->getConstRefTodEMFieldSourcedF()[i].print("dEMFieldSource dF Analytique");
-
+ // ipvcur->getConstRefTodEMFieldSourcedF()[i].print("dEMFieldSource dF Analytique");
+
+ /*
for (int k=0; k< numExtraDof; k++)
{
ipvcur->getConstRefTodFluxdField()[k][i].print("dFludx dField Analytique");
@@ -328,8 +329,8 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
for (unsigned int k = 0; k < numNonLocalVars; ++k)
{
Msg::Info("dp dExtraDofDiffusionField Analytique: %e",ipvcur->getRefTodLocalVariableDExtraDofDiffusionField()(k,i));
- }
- }*/
+ }*/
+ }
for (unsigned int i = 0; i < numCurlDof; ++i)
{
ipvcur->getConstRefTodPdVectorPotential()[i].print("dPdMagneticVectorPotential Analytique");
@@ -654,18 +655,19 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
Msg::Info("dEMFieldSource dp Numerique: %e",ipvcur->getConstRefTodEMFieldSourcedNonLocalVariable()(j,i));
}
}
- /*
+
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");
+ // 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");
ipvcur->getConstRefTodPdField()[i].print("dP d Field Numerique");
- ipvcur->getConstRefTodEMFieldSourcedF()[i].print("dEMFieldSource dF Numerique");
+ // ipvcur->getConstRefTodEMFieldSourcedF()[i].print("dEMFieldSource dF Numerique");
+ /*
for (int k=0; k< numExtraDof; k++)
{
ipvcur->getConstRefTodFluxdField()[k][i].print("dFludx dField Numerique");
@@ -694,8 +696,8 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
for (unsigned int k = 0; k < numNonLocalVars; ++k)
{
Msg::Info("dp dExtraDofDiffusionField Numerique: %e",ipvcur->getRefTodLocalVariableDExtraDofDiffusionField()(k,i));
- }
- }*/
+ }*/
+ }
for (unsigned int i = 0; i < numCurlDof; ++i)
{
const unsigned int extradof_T = 0; // Thermal field EM source