diff --git a/NonLinearSolver/internalPoints/ipTVEMullins.cpp b/NonLinearSolver/internalPoints/ipTVEMullins.cpp
index ade8c119bb09e2105e7826ca2639407a1b6da206..921b46c650f7146b5600f229c670738b6bd79be7 100644
--- a/NonLinearSolver/internalPoints/ipTVEMullins.cpp
+++ b/NonLinearSolver/internalPoints/ipTVEMullins.cpp
@@ -17,9 +17,35 @@ IPTVEMullins::IPTVEMullins(const J2IsotropicHardening* comp,
_T(298.15),_devDE(0.),_trDE(0.),_DE(0.),_DcorKirDT(0.),_DDcorKirDTT(0.),
_mechSrcTVE(0.),_DmechSrcTVEdT(0.),_DmechSrcTVEdF(0.),
_thermalEnergy(0.),_dtShift(0.),_DdtShiftDT(0.),_DDdtShiftDDT(0.){
+
+ _devOGammai.clear();
+ _trOGammai.clear();
+ _devDOGammaiDT.clear();
+ _trDOGammaiDT.clear();
+ _devOi.clear();
+ _trOi.clear();
+ _devDOiDT.clear();
+ _trDOiDT.clear();
+ _devDDOiDTT.clear();
+ _trDDOiDTT.clear();
+ for (int i=0; i < N; i++){
+ STensor3 el(0.);
+ _devOGammai.push_back(el);
+ _trOGammai.push_back(0.);
+ _devDOGammaiDT.push_back(el);
+ _trDOGammaiDT.push_back(0.);
+ _devOi.push_back(el);
+ _trOi.push_back(0.);
+ _devDOiDT.push_back(el);
+ _trDOiDT.push_back(0.);
+ _devDDOiDTT.push_back(el);
+ _trDDOiDTT.push_back(0.);
};
IPTVEMullins::IPTVEMullins(const IPTVEMullins& src): IPHyperViscoElastoPlastic(src),
+ _devOGammai(src._devOGammai), _trOGammai(src._trOGammai),
+ _devDOGammaiDT(src._devDOGammaiDT), _trDOGammaiDT(src._trDOGammaiDT),
+ _devOi(src._devOi), _trOi(src._trOi), _devDOiDT(src._devDOiDT), _trDOiDT(src._trDOiDT), _devDDOiDTT(src._devDDOiDTT), _trDDOiDTT(src._trDDOiDTT),
_T(src._T),_devDE(src._devDE), _trDE(src._trDE), _DE(src._DE), _DcorKirDT(src._DcorKirDT), _DDcorKirDTT(src._DDcorKirDTT),
_mechSrcTVE(src._mechSrcTVE),_DmechSrcTVEdT(src._DmechSrcTVEdT),_DmechSrcTVEdF(src._DmechSrcTVEdF),
_thermalEnergy(src._thermalEnergy),
@@ -44,6 +70,16 @@ IPTVEMullins& IPTVEMullins::operator=(const IPVariable& src)
_dtShift = psrc->_dtShift;
_DdtShiftDT = psrc->_DdtShiftDT;
_DDdtShiftDDT = psrc->_DDdtShiftDDT;
+ _devOGammai = psrc->_devOGammai;
+ _trOGammai = psrc->_trOGammai;
+ _devDOGammaiDT = psrc->_devDOGammaiDT;
+ _trDOGammaiDT = psrc->_trDOGammaiDT;
+ _devOi = psrc->_devOi;
+ _trOi = psrc->_trOi;
+ _devDOiDT = psrc->_devDOiDT;
+ _trDOiDT = psrc->_trDOiDT;
+ _devDDOiDTT = psrc->_devDDOiDTT;
+ _trDDOiDTT = psrc->_trDDOiDTT;
}
return *this;
}
@@ -68,4 +104,14 @@ void IPTVEMullins::restart(){
restartManager::restart(_dtShift);
restartManager::restart(_DdtShiftDT);
restartManager::restart(_DDdtShiftDDT);
+ restartManager::restart(_devOGammai);
+ restartManager::restart(_trOGammai);
+ restartManager::restart(_devDOGammaiDT);
+ restartManager::restart(_trDOGammaiDT);
+ restartManager::restart(_devOi);
+ restartManager::restart(_trOi);
+ restartManager::restart(_devDOiDT);
+ restartManager::restart(_trDOiDT);
+ restartManager::restart(_devDDOiDTT);
+ restartManager::restart(_trDDOiDTT);
}
diff --git a/NonLinearSolver/internalPoints/ipTVEMullins.h b/NonLinearSolver/internalPoints/ipTVEMullins.h
index f8c8dc019b09ad19ed0f55bc9615482c46b0471d..e5d0f6b1bb874df17a9a2134d23ef3e1184f6dd8 100644
--- a/NonLinearSolver/internalPoints/ipTVEMullins.h
+++ b/NonLinearSolver/internalPoints/ipTVEMullins.h
@@ -30,6 +30,20 @@ class IPTVEMullins : public IPHyperViscoElastoPlastic{
double _mechSrcTVE;
double _DmechSrcTVEdT;
STensor3 _DmechSrcTVEdF;
+
+ // Viscoelastic Strain
+ std::vector<STensor3> _devOGammai; // dev viscoelastic strain for each branch
+ std::vector<double> _trOGammai; // tr viscoelastic strain for each branch
+ std::vector<STensor3> _devDOGammaiDT;
+ std::vector<double> _trDOGammaiDT;
+
+ // Viscoelastic Stress
+ std::vector<STensor3> _devOi; // dev viscoelastic stress for each branch
+ std::vector<double> _trOi; // tr viscoelastic stress for each branch
+ std::vector<STensor3> _devDOiDT;
+ std::vector<double> _trDOiDT;
+ std::vector<STensor3> _devDDOiDTT;
+ std::vector<double> _trDDOiDTT;
// DcorKirDT
STensor3 _DcorKirDT;
@@ -50,6 +64,36 @@ class IPTVEMullins : public IPHyperViscoElastoPlastic{
virtual IPTVEMullins& operator=(const IPVariable& src);
virtual ~IPTVEMullins(){};
+ virtual std::vector<STensor3>& getRefToDevViscoElasticOverStrain() {return _devOGammai;};
+ virtual const std::vector<STensor3>& getConstRefToDevViscoElasticOverStrain() const {return _devOGammai;};
+
+ virtual std::vector<double>& getRefToTrViscoElasticOverStrain() {return _trOGammai;};
+ virtual const std::vector<double>& getConstRefToTrViscoElasticOverStrain() const {return _trOGammai;};
+
+ virtual std::vector<STensor3>& getRefToDevDOGammaiDT() {return _devDOGammaiDT;};
+ virtual const std::vector<STensor3>& getConstRefToDevDOGammaiDT() const {return _devDOGammaiDT;};
+
+ virtual std::vector<double>& getRefToTrDOGammaiDT() {return _trDOGammaiDT;};
+ virtual const std::vector<double>& getConstRefToTrDOGammaiDT() const {return _trDOGammaiDT;};
+
+ virtual std::vector<STensor3>& getRefToDevViscoElasticOverStress() {return _devOi;};
+ virtual const std::vector<STensor3>& getConstRefToDevViscoElasticOverStress() const {return _devOi;};
+
+ virtual std::vector<double>& getRefToTrViscoElasticOverStress() {return _trOi;};
+ virtual const std::vector<double>& getConstRefToTrViscoElasticOverStress() const {return _trOi;};
+
+ virtual std::vector<STensor3>& getRefToDevDOiDT() {return _devDOiDT;};
+ virtual const std::vector<STensor3>& getConstRefToDevDOiDT() const {return _devDOiDT;};
+
+ virtual std::vector<double>& getRefToTrDOiDT() {return _trDOiDT;};
+ virtual const std::vector<double>& getConstRefToTrDOiDT() const {return _trDOiDT;};
+
+ virtual std::vector<STensor3>& getRefToDevDDOiDTT() {return _devDDOiDTT;};
+ virtual const std::vector<STensor3>& getConstRefToDevDDOiDTT() const {return _devDDOiDTT;};
+
+ virtual std::vector<double>& getRefToTrDDOiDTT() {return _trDDOiDTT;};
+ virtual const std::vector<double>& getConstRefToTrDDOiDTT() const {return _trDDOiDTT;};
+
virtual STensor3& getRefToDcorKirDT() {return _DcorKirDT;};
virtual const STensor3& getConstRefToDcorKirDT() const {return _DcorKirDT;};
diff --git a/NonLinearSolver/materialLaw/mlawTVEMullins.cpp b/NonLinearSolver/materialLaw/mlawTVEMullins.cpp
index 80902fcaaed42063be57fa78d2ee7fb2066d66be..6e85a546c5c43c0639e6244b82464345410eafb4 100644
--- a/NonLinearSolver/materialLaw/mlawTVEMullins.cpp
+++ b/NonLinearSolver/materialLaw/mlawTVEMullins.cpp
@@ -121,11 +121,15 @@ void mlawTVEMullins::setViscoElasticNumberOfElement(const int N){
Msg::Info("Numer of Spring/Dashpot for viscoelastic model: %d",_N);
_Ki.clear(); _ki.clear();
_Gi.clear(); _gi.clear();
+
+ _nKi.clear(); _nGi.clear();
_Ki.resize(_N,0.);
_ki.resize(_N,0.);
_Gi.resize(_N,0.);
_gi.resize(_N,0.);
+ _nKi.resize(_N,0.);
+ _nGi.resize(_N,0.);
};
void mlawTVEMullins::setViscoElasticData(const int i, const double Ei, const double taui){
@@ -139,6 +143,9 @@ void mlawTVEMullins::setViscoElasticData(const int i, const double Ei, const dou
_ki[i-1] = taui; // /(3.*(1.-2.*_nu));
_Gi[i-1] = GG;
_gi[i-1] = taui; // /(2.*(1.+_nu));
+
+ _nKi[i-1] = KK/_K;
+ _nGi[i-1] = GG/_G;
Msg::Info("setting: element=%d Ki= %e ki= %e, Gi= %e, gi= %e",i-1,KK,taui,GG,taui);
}
@@ -150,6 +157,7 @@ void mlawTVEMullins::setViscoElasticData_Bulk(const int i, const double Ki, cons
else{
_Ki[i-1] = Ki;
_ki[i-1] = ki;
+ _nKi[i-1] = KK/_K;
// Msg::Info("setting: element=%d Ki = %e ki = %e, ",i-1,Ki,ki);
}
};
@@ -160,6 +168,7 @@ void mlawTVEMullins::setViscoElasticData_Shear(const int i, const double Gi, con
else{
_Gi[i-1] = Gi;
_gi[i-1] = gi;
+ _nGi[i-1] = GG/_G;
// Msg::Info("setting: element=%d Gi = %e gi = %e, ",i-1,Gi,gi);
}
};
@@ -169,6 +178,7 @@ void mlawTVEMullins::setViscoElasticData(const std::string filename){
if (file == NULL) Msg::Error("file: %s does not exist",filename.c_str());
_Ki.clear(); _ki.clear();
_Gi.clear(); _gi.clear();
+ _nKi.clear(); _nGi.clear();
_N = 0;
Msg::Info("reading viscoelastic input data");
@@ -176,13 +186,18 @@ void mlawTVEMullins::setViscoElasticData(const std::string filename){
int ok = fscanf(file,"%d",&_N);
Msg::Info("Numer of Maxwell elements: %d",_N);
double KK, k, GG, g;
+ double tempK, tempG = 0.;
for (int i=0; i<_N; i++){
ok = fscanf(file,"%lf %lf %lf %lf",&KK,&k,&GG,&g);
+ tempK = KK/_K;
+ tempG = GG/_G;
_Ki.push_back(KK);
_ki.push_back(k);
_Gi.push_back(GG);
_gi.push_back(g);
+ _nKi.push_back(tempK);
+ _nGi.push_back(tempG);
Msg::Info("Maxwell element %d: K[%d] = %e, k[%d] = %e, G[%d] = %e, g[%d] = %e",
i,i,KK,i,k,i,GG,i,g);
}
@@ -203,10 +218,157 @@ double mlawTVEMullins::ShiftFactorWLF(const double T, double *DaDT, double *DDaD
return a_WLF;
}
-double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double Tc) const{
+void mlawTVEMullins::shiftedTime(const double t, const double T0, const double T1, double *t_shift, double *Ddt_shiftDT, double *DDdt_shiftDT,
+ const int opt_time, const int opt_order) const{
+
+// Shift t -> t is either delta_t or t (final time). This function returns either t^* or delta_t^*.
+
+// opt_time -> 0 or 1 (default), depending on full-step or mid-step integration approximation for shifted time step
+// opt_order -> 1 or 2 (default) or 4, order of Guassian quadrature (1 doesnt work, 4 probably works)
+
+ double dt = this->getTimeStep();
+ double t1 = this->getTime();
+ double t0 = t1 - dt;
+
+ // Linearly Interpolate Temperture -> This means dT/dt is constant between tn and tn+1
+ double ts, t_x0, t_x1, t_x2 = 0.;
+ double t_x6, t_x7, t_x8, t_x9 = 0.;
+ if (opt_time == 0){
+ ts = 2;
+
+ t_x0 = (t1+t0)/2;
+
+ t_x1 = (t1+t0)/2 - dt/2 * (1/sqrt(3));
+ t_x2 = (t1+t0)/2 + dt/2 * (1/sqrt(3));
+
+ t_x6 = (t1+t0)/2 - dt/2 * sqrt(3/7 - 2/7*sqrt(6/5));
+ t_x7 = (t1+t0)/2 + dt/2 * sqrt(3/7 - 2/7*sqrt(6/5));
+ t_x8 = (t1+t0)/2 - dt/2 * sqrt(3/7 + 2/7*sqrt(6/5));
+ t_x9 = (t1+t0)/2 + dt/2 * sqrt(3/7 + 2/7*sqrt(6/5));
+ } else {
+ ts = 4;
+
+ t_x0 = (3*t1+t0)/4;
+
+ t_x1 = (3*t1+t0)/4 - dt/4 * (1/sqrt(3));
+ t_x2 = (3*t1+t0)/4 + dt/4 * (1/sqrt(3));
+
+ t_x6 = (3*t1+t0)/4 - dt/4 * sqrt(3/7 - 2/7*sqrt(6/5));
+ t_x7 = (3*t1+t0)/4 + dt/4 * sqrt(3/7 - 2/7*sqrt(6/5));
+ t_x8 = (3*t1+t0)/4 - dt/4 * sqrt(3/7 + 2/7*sqrt(6/5));
+ t_x9 = (3*t1+t0)/4 + dt/4 * sqrt(3/7 + 2/7*sqrt(6/5));
+ }
+
+ double T_x0 = linearInterpTemp(t_x0, t0, t1, T0, T1);
+ double T_x1 = linearInterpTemp(t_x1, t0, t1, T0, T1);
+ double T_x2 = linearInterpTemp(t_x2, t0, t1, T0, T1);
+ double T_x6 = linearInterpTemp(t_x6, t0, t1, T0, T1);
+ double T_x7 = linearInterpTemp(t_x7, t0, t1, T0, T1);
+ double T_x8 = linearInterpTemp(t_x8, t0, t1, T0, T1);
+ double T_x9 = linearInterpTemp(t_x9, t0, t1, T0, T1);
+ // double T_x1 = 1/dt*( (t_1 - t_x1) )*T_0 + 1/dt*( (t_x1 - t_0) )*T_1 ;
+ // double T_x2 = 1/dt*( (t_1 - t_x2) )*T_0 + 1/dt*( (t_x2 - t_0) )*T_1 ;
+
+ if (t_shift!=NULL){
+
+ // double a_T1 = ShiftFactorWLF(T1);
+
+ if (opt_order == 1){
+ // 1 point Gauss rule
+ double a_x0 = ShiftFactorWLF(T_x0);
+ *t_shift = t/ts * (1/a_x0);
+ } else if (opt_order == 2){
+ // 2 point Gauss rule
+ double a_x1 = ShiftFactorWLF(T_x1);
+ double a_x2 = ShiftFactorWLF(T_x2);
+ *t_shift = t/ts * (1/a_x1 + 1/a_x2);
+ } else {
+ // 4 point Gauss rule
+ double a_x6 = ShiftFactorWLF(T_x6);
+ double a_x7 = ShiftFactorWLF(T_x7);
+ double a_x8 = ShiftFactorWLF(T_x8);
+ double a_x9 = ShiftFactorWLF(T_x9);
+ *t_shift = t/ts * (((18+sqrt(30))/36)*(1/a_x6 + 1/a_x7) + ((18-sqrt(30))/36)*(1/a_x8 + 1/a_x9));
+ }
+ }
+
+ if (Ddt_shiftDT!=NULL){
+
+ // double DaDT_T1 = 0.; double a_T1 = ShiftFactorWLF(T1,&DaDT_T1);
+
+ if (dt > 0.){
+ if (opt_order == 1){
+ // 1 point Gauss rule
+ double DaDT_x0 = 0.;
+ double a_x0 = ShiftFactorWLF(T_x0,&DaDT_x0);
+ *Ddt_shiftDT = t/ts * (DaDT_x0 * (t_x0-t0)/dt );
+ } else if (opt_order == 2){
+ // 2 point Gauss rule
+ double DaDT_x1, DaDT_x2 = 0.;
+ double a_x1 = ShiftFactorWLF(T_x1,&DaDT_x1);
+ double a_x2 = ShiftFactorWLF(T_x2,&DaDT_x2);
+ *Ddt_shiftDT = t/ts * (DaDT_x1 * (t_x1-t0)/dt + DaDT_x2 * (t_x2-t0)/dt);
+ } else {
+ // 4 point Gauss rule
+ double DaDT_x6, DaDT_x7, DaDT_x8, DaDT_x9 = 0.;
+ double a_x6 = ShiftFactorWLF(T_x6,&DaDT_x6);
+ double a_x7 = ShiftFactorWLF(T_x7,&DaDT_x7);
+ double a_x8 = ShiftFactorWLF(T_x8,&DaDT_x8);
+ double a_x9 = ShiftFactorWLF(T_x9,&DaDT_x9);
+ *Ddt_shiftDT = t/ts * ( ((18+sqrt(30))/36)*(DaDT_x6 * (t_x6-t0)/dt + DaDT_x7 * (t_x7-t0)/dt ) +
+ ((18-sqrt(30))/36)*(DaDT_x8 * (t_x8-t0)/dt + DaDT_x9 * (t_x9-t0)/dt ) );
+ }
+ } else { *Ddt_shiftDT = 0.;}
+ }
+
+ if (DDdt_shiftDT!=NULL){
+
+ // double DDaDT_T1 = 0.; double a_T1 = ShiftFactorWLF(T1,NULL,&DDaDT_T1);
+
+ if (dt > 0.){
+ if (opt_order == 1){
+ // 1 point Gauss rule
+ double DDaDT_x0 = 0.;
+ double a_x0 = ShiftFactorWLF(T_x0,NULL,&DDaDT_x0);
+ *DDdt_shiftDT = t/ts * ( DDaDT_x0 * pow( ((t_x0-t0)/dt) ,2 ) );
+ } else if (opt_order == 2){
+ // 2 point Gauss rule
+ double DDaDT_x1, DDaDT_x2 = 0.;
+ double a_x1 = ShiftFactorWLF(T_x1,NULL,&DDaDT_x1);
+ double a_x2 = ShiftFactorWLF(T_x2,NULL,&DDaDT_x2);
+ *DDdt_shiftDT = t/ts * ( DDaDT_x1 * pow( ((t_x1-t0)/dt) ,2 ) + DDaDT_x2 * pow( ((t_x2-t0)/dt) ,2 ) );
+ } else {
+ // 4 point Gauss rule
+ double DDaDT_x6, DDaDT_x7, DDaDT_x8, DDaDT_x9 = 0.;
+ double a_x6 = ShiftFactorWLF(T_x6,NULL,&DDaDT_x6);
+ double a_x7 = ShiftFactorWLF(T_x7,NULL,&DDaDT_x7);
+ double a_x8 = ShiftFactorWLF(T_x8,NULL,&DDaDT_x8);
+ double a_x9 = ShiftFactorWLF(T_x9,NULL,&DDaDT_x9);
+ *Ddt_shiftDT = t/ts * ( ((18+sqrt(30))/36)*( DDaDT_x6 * pow( (t_x6-t0)/dt, 2 ) + DDaDT_x7 * pow( (t_x7-t0)/dt, 2 )) +
+ ((18-sqrt(30))/36)*( DDaDT_x8 * pow( (t_x8-t0)/dt, 2 ) + DDaDT_x9 * pow( (t_x9-t0)/dt, 2 )) );
+ }
+ } else { *DDdt_shiftDT = 0.;}
+ }
+}
+
+
+double mlawTVEMullins::linearInterpTemp(const double t_interp, const double t0, const double t1, const double T0, const double T1) const{
- // Update the Properties to the current temperature (see below which ones are being used)
- double KT, GT, AlphaT; getK(KT,Tc); getG(GT,Tc); getAlpha(AlphaT,Tc);
+ // This function linearly interpolates temperature for the given intermediate time -> t_interp is tn + fraction of step (i.e., midstep = tn + 0.5*dt)
+ // Also, linear interp -> constant gradients
+ double dt = t1 - t0;
+ double T_interp;
+ if (dt > 0.){
+ T_interp = 1/dt*( (t1 - t_interp) )*T0 + 1/dt*( (t_interp - t0) )*T1;
+ }
+ else { T_interp = T0;}
+ return T_interp;
+}
+
+double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double T) const{
+
+ // Update the Properties to the current temperature
+ double KT, GT, AlphaT; getK(KT,T); getG(GT,T); getAlpha(AlphaT,T);
double Psy_mech = 0.;
@@ -216,7 +378,7 @@ double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double
STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
// Get thermal strain
- double Eth = 3.*AlphaT*(Tc-_Tinitial);
+ double Eth = 3.*AlphaT*(T-_Tinitial);
// Get Equilibrium freeEnergy_mech
Psy_mech = KT*0.5*(trEe-Eth)*(trEe-Eth) + GT*STensorOperation::doubledot(devEe,devEe);
@@ -247,7 +409,10 @@ double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double
}
void mlawTVEMullins::corKirInfinity(const STensor3& devEe, const double& trEe, const double T0, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const{
-
+ double KT, GT, AlphaT;
+ getK(KT,T); getG(GT,T); getAlpha(AlphaT,T);
+ CorKirDevInf = 2*GT*devEe;
+ CorKirTrInf = KT*(trEe - 3*AlphaT*(T-_Tinitial));
}
void mlawTVEMullins::getTVEdCorKirDT(const IPTVEMullins *q0, IPTVEMullins *q1, const double& T0, const double& T) const{
@@ -263,7 +428,7 @@ void mlawTVEMullins::getMechSourceTVE(const IPTVEMullins *q0, IPTVEMullins *q1,
void mlawTVEMullins::ThermoViscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
const IPTVEMullins *q0, IPTVEMullins *q1,
double& Ke, double& Ge,
- const double T0, const double T1, const bool calculateStress) const{
+ const double T0, const double T, const bool calculateStress) const{
// Get LogStrain
@@ -277,7 +442,7 @@ void mlawTVEMullins::ThermoViscoElasticPredictor(const STensor3& Ee, const STens
// Get corKirInfinity
static double CorKirTrInf;
static STensor3 CorKirDevInf;
- corKirInfinity(devEe,trEe,T0,T1,CorKirDevInf,CorKirTrInf);
+ corKirInfinity(devEe,trEe,T0,T,CorKirDevInf,CorKirTrInf);
// Initialise CorKirStress - Dev
static STensor3 devK;
@@ -295,9 +460,9 @@ void mlawTVEMullins::ThermoViscoElasticPredictor(const STensor3& Ee, const STens
double dt_shift_1, dt_shift_2 = 0.;
double Ddt_shiftDT, DDdt_shiftDT = 0.;
// dt_shift = dt;
- shiftedTime(dt, T0, T1, &dt_shift_1,NULL,NULL,0);
- shiftedTime(dt, T0, T1, &dt_shift_2,&Ddt_shiftDT,&DDdt_shiftDT);
- // double aT_1 = ShiftFactorWLF(T1);
+ shiftedTime(dt, T0, T, &dt_shift_1,NULL,NULL,0);
+ shiftedTime(dt, T0, T, &dt_shift_2,&Ddt_shiftDT,&DDdt_shiftDT);
+ // double aT_1 = ShiftFactorWLF(T);
std::vector<double> KiT, GiT, DKiDT, DGiDT, DDKiDT, DDGiDT;
KiT.resize(_N,0.); GiT.resize(_N,0.);
diff --git a/NonLinearSolver/materialLaw/mlawTVEMullins.h b/NonLinearSolver/materialLaw/mlawTVEMullins.h
index 195b9c7dd20b96ce21d7161ef234e775bd1b3682..65ada9227b82d5512381cfc48d1f23629737aed1 100644
--- a/NonLinearSolver/materialLaw/mlawTVEMullins.h
+++ b/NonLinearSolver/materialLaw/mlawTVEMullins.h
@@ -49,8 +49,8 @@ class mlawTVEMullins : public mlawPowerYieldHyper{
scalarFunction* _temFunc_Cp; // specific heat temperature function
// Scalars - Additional constants
- std::vector<double> nKi; // branch free energy scalars
- std::vector<double> nGi; // branch free energy scalars
+ std::vector<double> _nKi; // normalized branch Bulk Moduli
+ std::vector<double> _nGi; // normalized branch Shear Moduli
// Shift factor - Everything related to it
double _C1; // WLF Universal Constant 1 - set in .py file
@@ -61,7 +61,11 @@ class mlawTVEMullins : public mlawPowerYieldHyper{
protected:
// Const Functions for constitutive function - To update any variable, pass it as an argument by reference.
virtual double ShiftFactorWLF(const double T, double *DaDT = NULL, double *DDaDDT = NULL) const ;
- virtual double freeEnergyMechanical(const IPTVEMullins& q, const double Tc) const ;
+ virtual void shiftedTime(const double t, const double T0, const double T1,
+ double *t_shift = NULL, double *Ddt_shiftDT = NULL, double *DDdt_shiftDT = NULL,
+ const int opt_time = 1, const int opt_order = 2) const;
+ virtual double linearInterpTemp(const double t_interp, const double t0, const double t1, const double T0, const double T1) const;
+ virtual double freeEnergyMechanical(const IPTVEMullins& q, const double T) const ;
void corKirInfinity(const STensor3& devEe, const double& trEe, const double T0, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const;