diff --git a/NonLinearSolver/internalPoints/CMakeLists.txt b/NonLinearSolver/internalPoints/CMakeLists.txt
index 174256314067b55dc46aedc208c330c14e29b92f..9eb0f884c769a84b1c75fdf1357540bb201ec4f6 100644
--- a/NonLinearSolver/internalPoints/CMakeLists.txt
+++ b/NonLinearSolver/internalPoints/CMakeLists.txt
@@ -60,6 +60,7 @@ set(SRC
ipNonLinearTVM.cpp
ipNonLinearTVE.cpp
ipNonLinearTVP.cpp
+ ipTVEMullins.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
diff --git a/NonLinearSolver/internalPoints/ipTVEMullins.cpp b/NonLinearSolver/internalPoints/ipTVEMullins.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ade8c119bb09e2105e7826ca2639407a1b6da206
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipTVEMullins.cpp
@@ -0,0 +1,71 @@
+//
+// C++ Interface: Material Law
+//
+// Description: IPNonLinearTVM (IP Thermo-ViscoElasto-ViscoPlasto Law) with Non-Local Damage Interface (soon.....)
+//
+// Author: <Ujwal Kishore J - FLE_Knight>, (C) 2022
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "ipTVEMullins.h"
+
+IPTVEMullins::IPTVEMullins(const J2IsotropicHardening* comp,
+ const J2IsotropicHardening* trac,
+ const kinematicHardening* kin, const int N): IPHyperViscoElastoPlastic(comp,trac,kin,N),
+ _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.){
+};
+
+IPTVEMullins::IPTVEMullins(const IPTVEMullins& src): IPHyperViscoElastoPlastic(src),
+ _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),
+ _dtShift(src._dtShift),_DdtShiftDT(src._DdtShiftDT),_DDdtShiftDDT(src._DDdtShiftDDT){};
+
+IPTVEMullins& IPTVEMullins::operator=(const IPVariable& src)
+{
+ IPHyperViscoElastoPlastic::operator=(src);
+ const IPTVEMullins* psrc = dynamic_cast<const IPTVEMullins*>(&src);
+ if(psrc != NULL)
+ {
+ _T = psrc->_T;
+ _mechSrcTVE = psrc->_mechSrcTVE;
+ _DmechSrcTVEdT = psrc->_DmechSrcTVEdT;
+ _DmechSrcTVEdF = psrc->_DmechSrcTVEdF;
+ _devDE = psrc->_devDE;
+ _trDE = psrc->_trDE;
+ _DE = psrc->_DE;
+ _DcorKirDT = psrc->_DcorKirDT;
+ _DDcorKirDTT = psrc->_DDcorKirDTT;
+ _thermalEnergy = psrc->_thermalEnergy;
+ _dtShift = psrc->_dtShift;
+ _DdtShiftDT = psrc->_DdtShiftDT;
+ _DDdtShiftDDT = psrc->_DDdtShiftDDT;
+ }
+ return *this;
+}
+
+double IPTVEMullins::defoEnergy() const
+{
+ return IPHyperViscoElastic::defoEnergy();
+}
+
+void IPTVEMullins::restart(){
+ IPHyperViscoElastoPlastic::restart();
+ restartManager::restart(_devDE);
+ restartManager::restart(_trDE);
+ restartManager::restart(_DE);
+ restartManager::restart(_DcorKirDT);
+ restartManager::restart(_DDcorKirDTT);
+ restartManager::restart(_mechSrcTVE);
+ restartManager::restart(_DmechSrcTVEdT);
+ restartManager::restart(_DmechSrcTVEdF);
+ restartManager::restart(_T);
+ restartManager::restart(_thermalEnergy);
+ restartManager::restart(_dtShift);
+ restartManager::restart(_DdtShiftDT);
+ restartManager::restart(_DDdtShiftDDT);
+}
diff --git a/NonLinearSolver/internalPoints/ipTVEMullins.h b/NonLinearSolver/internalPoints/ipTVEMullins.h
new file mode 100644
index 0000000000000000000000000000000000000000..f8c8dc019b09ad19ed0f55bc9615482c46b0471d
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipTVEMullins.h
@@ -0,0 +1,88 @@
+//
+// C++ Interface: Material Law
+//
+// Description: IPNonLinearTVE (IP Thermo-ViscoElasto-ViscoPlasto Law) with Non-Local Damage Interface (soon.....)
+//
+// Author: <Ujwal Kishore J - FLE_Knight>, (C) 2022
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef IPTVEMULLINS_H_
+#define IPTVEMULLINS_H_
+#include "ipHyperelastic.h"
+#include "j2IsotropicHardening.h"
+#include "kinematicHardening.h"
+#include "STensor3.h"
+#include "STensor43.h"
+
+
+class IPTVEMullins : public IPHyperViscoElastoPlastic{
+ public:
+ // Temperature history
+ double _thermalEnergy;
+ double _T;
+ // double _fracEnergy;
+ STensor3 _devDE;
+ double _trDE;
+ STensor3 _DE;
+ double _mechSrcTVE;
+ double _DmechSrcTVEdT;
+ STensor3 _DmechSrcTVEdF;
+
+ // DcorKirDT
+ STensor3 _DcorKirDT;
+ STensor3 _DDcorKirDTT;
+
+ // Additional Internal Variables
+
+ // Additional variables
+ double _dtShift;
+ double _DdtShiftDT;
+ double _DDdtShiftDDT;
+
+ public:
+ IPTVEMullins(const J2IsotropicHardening* comp,
+ const J2IsotropicHardening* trac,
+ const kinematicHardening* kin, const int N);
+ IPTVEMullins(const IPTVEMullins& src);
+ virtual IPTVEMullins& operator=(const IPVariable& src);
+ virtual ~IPTVEMullins(){};
+
+ virtual STensor3& getRefToDcorKirDT() {return _DcorKirDT;};
+ virtual const STensor3& getConstRefToDcorKirDT() const {return _DcorKirDT;};
+
+ virtual STensor3& getRefToDDcorKirDTT() {return _DDcorKirDTT;};
+ virtual const STensor3& getConstRefToDDcorKirDTT() const {return _DDcorKirDTT;};
+
+ virtual double& getRefToTrDE() {return _trDE;};
+ virtual const double& getConstRefToTrDE() const {return _trDE;};
+ virtual STensor3& getRefToDevDE() {return _devDE;};
+ virtual const STensor3& getConstRefToDevDE() const {return _devDE;};
+
+ virtual STensor3& getRefToDE() {return _DE;};
+ virtual const STensor3& getConstRefToDE() const {return _DE;};
+
+ virtual double &getRefToMechSrcTVE() {return _mechSrcTVE;}
+ virtual const double &getConstRefToMechSrcTVE() const {return _mechSrcTVE;}
+ virtual double &getRefTodMechSrcTVEdT() {return _DmechSrcTVEdT;}
+ virtual const double &getConstRefTodMechSrcTVEdT() const {return _DmechSrcTVEdT;}
+ virtual STensor3 &getRefTodMechSrcTVEdF() {return _DmechSrcTVEdF;}
+ virtual const STensor3 &getConstRefTodMechSrcTVEdF() const {return _DmechSrcTVEdF;}
+
+ virtual double& getShiftedTimeStep() {return _dtShift;};
+ virtual const double& getConstShiftedTimeStep() const {return _dtShift;};
+ virtual double& getShiftedTimeStepTempDerivative() {return _DdtShiftDT;};
+ virtual const double& getConstShiftedTimeStepTempDerivative() const {return _DdtShiftDT;};
+ virtual double& getShiftedTimeStepTempDoubleDerivative() {return _DDdtShiftDDT;};
+ virtual const double& getConstShiftedTimeStepTempDoubleDerivative() const {return _DDdtShiftDDT;};
+
+ virtual double defoEnergy() const;
+ virtual double getThermalEnergy() const {return _thermalEnergy;}; // Added _thermalEnergy
+ // virtual double getConstRefToFractureEnergy() const {return _fracEnergy;}; // Added - Unused (edit this because its neither "const" nor "ref" here)
+
+ virtual IPVariable* clone() const{return new IPTVEMullins(*this);};
+ virtual void restart();
+};
+#endif // IPTVEMullins.h
diff --git a/NonLinearSolver/materialLaw/CMakeLists.txt b/NonLinearSolver/materialLaw/CMakeLists.txt
index 9da69791fd17edb91f649091ca268300272dc4f0..fdaa8c6db8e82fe50053ebf374eb23b22f5a9add 100644
--- a/NonLinearSolver/materialLaw/CMakeLists.txt
+++ b/NonLinearSolver/materialLaw/CMakeLists.txt
@@ -99,6 +99,7 @@ set(SRC
mlawNonLinearTVM.cpp
mlawNonLinearTVE.cpp
mlawNonLinearTVP.cpp
+ mlawTVEMullins.cpp
# src/op_eshelby.cpp
# Headers without cpp change this ??
ID.h
diff --git a/NonLinearSolver/materialLaw/mlawTVEMullins.cpp b/NonLinearSolver/materialLaw/mlawTVEMullins.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..80902fcaaed42063be57fa78d2ee7fb2066d66be
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mlawTVEMullins.cpp
@@ -0,0 +1,1088 @@
+//
+// C++ Interface: Material Law
+//
+// Description: Non-Linear Thermo-Visco-Mechanics (Thermo-ViscoElasto-ViscoPlasto Law) with Non-Local Damage Interface (soon.....)
+//
+// Author: <Ujwal Kishore J - FLE_Knight>, (C) 2022; <Van Dung Nguyen>, (C) 2014
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "mlawHyperelastic.h"
+#include "STensorOperations.h"
+#include "nonLinearMechSolver.h"
+#include "mlawTVEMullins.h"
+
+// Constructor
+mlawTVEMullins::mlawTVEMullins(const int num,const double E,const double nu, const double rho, const double tol,
+ const double Tinitial, const double Alpha, const double KThCon, const double Cp,
+ const bool matrixbyPerturbation, const double pert, const bool thermalEstimationPreviousConfig):
+ mlawPowerYieldHyper(num, E, nu, rho, tol, matrixbyPerturbation, pert),
+ _Tinitial(Tinitial),_scalarAlpha(Alpha),_scalarK(KThCon),_Cp(Cp),
+ _thermalEstimationPreviousConfig(thermalEstimationPreviousConfig){
+
+ _G = _mu;
+
+ _Tref = 298.;
+ _C1 = 0.;
+ _C2 = 0.;
+
+ // _scalarAlpha = 0.; _scalarK = 0.; _Cp = 0.;
+ _tensorAlpha = 0.;
+ _tensorAlpha(0,0)= _scalarAlpha; _tensorAlpha(1,1)= _scalarAlpha; _tensorAlpha(2,2)= _scalarAlpha;
+
+ // to be unifom in DG3D q = k' gradT with k'=-k instead of q=-kgradT
+ _tensorK = 0.;
+ _tensorK(0,0)= - _scalarK; _tensorK(1,1)= - _scalarK; _tensorK(2,2)= - _scalarK;
+
+ _temFunc_K = new constantScalarFunction(1.);
+ _temFunc_G = new constantScalarFunction(1.);
+ _temFunc_Alpha = new constantScalarFunction(1.);
+ _temFunc_KThCon = new constantScalarFunction(1.);
+ _temFunc_Cp = new constantScalarFunction(1.);
+};
+
+// Copy Constructor
+mlawTVEMullins::mlawTVEMullins(const mlawTVEMullins& src): mlawPowerYieldHyper(src),
+ _Tinitial(src._Tinitial),_scalarAlpha(src._scalarAlpha),_scalarK(src._scalarK),_Cp(src._Cp),
+ _thermalEstimationPreviousConfig(src._thermalEstimationPreviousConfig){
+
+ _G = src._G;
+ _Tref = src._Tref;
+ _testFlag = src._testFlag;
+ _C1 = src._C1;
+ _C2 = src._C2;
+ _tensorAlpha = src._tensorAlpha;
+ _tensorK = src._tensorK;
+
+ _temFunc_K = NULL; // bulk modulus
+ if (src._temFunc_K != NULL){ _temFunc_K = src._temFunc_K->clone();}
+
+ _temFunc_G = NULL; // shear modulus
+ if (src._temFunc_G!=NULL){ _temFunc_G = src._temFunc_G->clone();}
+
+ _temFunc_Alpha = NULL; // thermal dilatation coeffecient
+ if (src._temFunc_Alpha!=NULL){ _temFunc_Alpha = src._temFunc_Alpha->clone();}
+
+ _temFunc_KThCon = NULL; // thermal conductivity
+ if (src._temFunc_KThCon!=NULL){ _temFunc_KThCon = src._temFunc_KThCon->clone();}
+
+ _temFunc_Cp = NULL; // specific heat
+ if (src._temFunc_Cp!=NULL){ _temFunc_Cp = src._temFunc_Cp->clone();}
+};
+
+
+mlawTVEMullins& mlawTVEMullins::operator=(const materialLaw& source){
+
+ mlawPowerYieldHyper::operator=(source);
+ const mlawTVEMullins* src =dynamic_cast<const mlawTVEMullins*>(&source);
+ if(src != NULL){
+
+ _Tinitial = src->_Tinitial; _scalarAlpha = src->_scalarAlpha; _scalarK = src->_scalarK; _Cp = src->_Cp;
+ _thermalEstimationPreviousConfig = src->_thermalEstimationPreviousConfig;
+
+ _G = src->_G;
+ _Tref = src->_Tref;
+ _testFlag = src->_testFlag;
+ _C1 = src->_C1;
+ _C2 = src->_C2;
+ _tensorAlpha = src->_tensorAlpha;
+ _tensorK = src->_tensorK;
+
+ if(_temFunc_K != NULL) delete _temFunc_K; // bulk modulus
+ if (src->_temFunc_K != NULL){ _temFunc_K = src->_temFunc_K->clone();}
+
+ if(_temFunc_G != NULL) delete _temFunc_G; // shear modulus
+ if (src->_temFunc_G!=NULL){ _temFunc_G = src->_temFunc_G->clone();}
+
+ if(_temFunc_Alpha != NULL) delete _temFunc_Alpha; // thermal dilatation coeffecient
+ if (src->_temFunc_Alpha!=NULL){ _temFunc_Alpha = src->_temFunc_Alpha->clone();}
+
+ if(_temFunc_KThCon != NULL) delete _temFunc_KThCon; // thermal conductivity
+ if (src->_temFunc_KThCon!=NULL){ _temFunc_KThCon = src->_temFunc_KThCon->clone();}
+
+ if(_temFunc_Cp != NULL) delete _temFunc_Cp; // specific heat
+ if (src->_temFunc_Cp!=NULL){ _temFunc_Cp = src->_temFunc_Cp->clone();}
+ }
+ return *this;
+}
+
+mlawTVEMullins::~mlawTVEMullins(){
+ if(_temFunc_K != NULL) delete _temFunc_K; _temFunc_K = NULL; // bulk modulus;
+ if(_temFunc_G != NULL) delete _temFunc_G; _temFunc_G = NULL; // shear modulus;
+ if(_temFunc_Alpha != NULL) delete _temFunc_Alpha; _temFunc_Alpha = NULL; // thermal dilatation coeff;
+ if(_temFunc_KThCon != NULL) delete _temFunc_KThCon; _temFunc_KThCon = NULL; // thermal conductivity;
+ if(_temFunc_Cp != NULL) delete _temFunc_Cp; _temFunc_Cp = NULL; // specific heat;
+};
+
+void mlawTVEMullins::setViscoElasticNumberOfElement(const int N){
+ _N = N;
+ Msg::Info("Numer of Spring/Dashpot for viscoelastic model: %d",_N);
+ _Ki.clear(); _ki.clear();
+ _Gi.clear(); _gi.clear();
+
+ _Ki.resize(_N,0.);
+ _ki.resize(_N,0.);
+ _Gi.resize(_N,0.);
+ _gi.resize(_N,0.);
+};
+
+void mlawTVEMullins::setViscoElasticData(const int i, const double Ei, const double taui){
+ if (i> _N or i<1)
+ Msg::Error("This setting is invalid %d > %d",i,_N);
+ else{
+ double KK = Ei/(3.*(1.-2.*_nu));
+ double GG = Ei/(2.*(1.+_nu));
+
+ _Ki[i-1] = KK;
+ _ki[i-1] = taui; // /(3.*(1.-2.*_nu));
+ _Gi[i-1] = GG;
+ _gi[i-1] = taui; // /(2.*(1.+_nu));
+
+ Msg::Info("setting: element=%d Ki= %e ki= %e, Gi= %e, gi= %e",i-1,KK,taui,GG,taui);
+ }
+};
+
+void mlawTVEMullins::setViscoElasticData_Bulk(const int i, const double Ki, const double ki){
+ if (i> _N or i<1)
+ Msg::Error("This setting is invalid %d > %d",i,_N);
+ else{
+ _Ki[i-1] = Ki;
+ _ki[i-1] = ki;
+ // Msg::Info("setting: element=%d Ki = %e ki = %e, ",i-1,Ki,ki);
+ }
+};
+
+void mlawTVEMullins::setViscoElasticData_Shear(const int i, const double Gi, const double gi){
+ if (i> _N or i<1)
+ Msg::Error("This setting is invalid %d > %d",i,_N);
+ else{
+ _Gi[i-1] = Gi;
+ _gi[i-1] = gi;
+ // Msg::Info("setting: element=%d Gi = %e gi = %e, ",i-1,Gi,gi);
+ }
+};
+
+void mlawTVEMullins::setViscoElasticData(const std::string filename){
+ FILE* file = fopen(filename.c_str(),"r");
+ if (file == NULL) Msg::Error("file: %s does not exist",filename.c_str());
+ _Ki.clear(); _ki.clear();
+ _Gi.clear(); _gi.clear();
+ _N = 0;
+ Msg::Info("reading viscoelastic input data");
+
+ while (!feof(file)){
+ int ok = fscanf(file,"%d",&_N);
+ Msg::Info("Numer of Maxwell elements: %d",_N);
+ double KK, k, GG, g;
+
+ for (int i=0; i<_N; i++){
+ ok = fscanf(file,"%lf %lf %lf %lf",&KK,&k,&GG,&g);
+ _Ki.push_back(KK);
+ _ki.push_back(k);
+ _Gi.push_back(GG);
+ _gi.push_back(g);
+ 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);
+ }
+ }
+ fclose(file);
+};
+
+double mlawTVEMullins::ShiftFactorWLF(const double T, double *DaDT, double *DDaDDT) const{
+
+ double a_WLF = exp(- (_C1*(T - _Tref)) / (_C2 + T - _Tref) );
+
+ if(DaDT!=NULL){
+ *DaDT = 1/a_WLF * (_C1*_C2) / pow( (_C2 + T-_Tref) ,2);
+ }
+ if(DDaDDT!=NULL){
+ *DDaDDT = 1/a_WLF * (_C1*_C2) * (_C1*_C2 - 2*(_C2+ T-_Tref)) / pow( (_C2 + T-_Tref) ,4);
+ }
+ return a_WLF;
+}
+
+double mlawTVEMullins::freeEnergyMechanical(const IPTVEMullins& q, const double Tc) 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);
+
+ double Psy_mech = 0.;
+
+ // Get trEe and devEe
+ double trEe;
+ static STensor3 devEe;
+ STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
+
+ // Get thermal strain
+ double Eth = 3.*AlphaT*(Tc-_Tinitial);
+
+ // Get Equilibrium freeEnergy_mech
+ Psy_mech = KT*0.5*(trEe-Eth)*(trEe-Eth) + GT*STensorOperation::doubledot(devEe,devEe);
+
+ // Add Branch freeEnergy_mech
+ if ((_Ki.size() > 0) or (_Gi.size() > 0)){
+
+ std::vector<double> KiT,GiT; KiT.resize(_N,0.); GiT.resize(_N,0.);
+ KiT = _Ki;
+ GiT = _Gi;
+
+ for (int i=0; i<_N; i++){
+ static STensor3 devEebranch;
+ devEebranch = devEe;
+
+ // devEebranch -= q._A[i];
+ // Psy_mech += KiT[i]*0.5*(trEe-q._B[i])*(trEe-q._B[i])+GiT[i]*STensorOperation::doubledot(devEebranch,devEebranch);
+
+ double trGamma = q._trOGammai[i];
+ static STensor3 devGamma = q._devOGammai[i];
+ devEebranch -= devGamma;
+
+ Psy_mech += KiT[i]*0.5*(trEe-trGamma)*(trEe-trGamma) + GiT[i]*STensorOperation::doubledot(devEebranch,devEebranch);
+ }
+ }
+
+ return Psy_mech;
+}
+
+void mlawTVEMullins::corKirInfinity(const STensor3& devEe, const double& trEe, const double T0, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const{
+
+}
+
+void mlawTVEMullins::getTVEdCorKirDT(const IPTVEMullins *q0, IPTVEMullins *q1, const double& T0, const double& T) const{
+
+}
+
+void mlawTVEMullins::getMechSourceTVE(const IPTVEMullins *q0, IPTVEMullins *q1, const double& T0, const double& T,
+ const double& DKDTsum, const double& DGDTsum, const STensor43& DEeDFe,
+ double *Wm, STensor3 *dWmdF, double *dWmdT) const{
+
+};
+
+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{
+
+
+ // Get LogStrain
+ static STensor3 DE, devDE, devEe;
+ static double trDE, trEe;
+ DE = Ee;
+ DE -= Ee0;
+ STensorOperation::decomposeDevTr(DE,devDE,trDE);
+ STensorOperation::decomposeDevTr(Ee,devEe,trEe);
+
+ // Get corKirInfinity
+ static double CorKirTrInf;
+ static STensor3 CorKirDevInf;
+ corKirInfinity(devEe,trEe,T0,T1,CorKirDevInf,CorKirTrInf);
+
+ // Initialise CorKirStress - Dev
+ static STensor3 devK;
+ devK = CorKirDevInf;
+
+ // Initialise CorKirStress - Tr
+ static double p;
+ p = CorKirTrInf;
+
+ // Main Loop for Maxwell Branches
+ if ((_Ki.size() > 0) or (_Gi.size() > 0)){
+
+ // Get ShiftFactor and shifted time increment
+ double dt = this->getTimeStep();
+ 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);
+
+ std::vector<double> KiT, GiT, DKiDT, DGiDT, DDKiDT, DDGiDT;
+ KiT.resize(_N,0.); GiT.resize(_N,0.);
+ DKiDT.resize(_N,0.); DGiDT.resize(_N,0.); DDKiDT.resize(_N,0.); DDGiDT.resize(_N,0.);
+
+
+ // Dev Q_i
+ std::vector<STensor3>& devOi = q1->getRefToDevViscoElasticOverStress(); // devOi.resize(_N,0.);
+ std::vector<STensor3>& devOGammai = q1->getRefToDevViscoElasticOverStrain();
+
+ for (int i=0; i<_Gi.size(); i++){
+
+ double dtg = dt_shift_1/(_gi[i]);
+ double dtg2 = dt_shift_2/(_gi[i]);
+ double expdtg = exp(-dtg);
+ // double expdtgby2 = exp(-dtg/2.);
+ double expdtgby2 = exp(-dtg2);
+
+ getGi(GiT,T_set); // Removed 2*_gi[i]
+ DGiDT[i] = -_Gi[i]/(_gi[i]) * Ddt_shiftDT*expdtgby2;
+ DDGiDT[i] = -_Gi[i]/(_gi[i]) * (DDdt_shiftDT - 1/(_gi[i])*pow(Ddt_shiftDT,2))*expdtgby2;
+
+ // For Single Convolution
+ // Ge for DcorKirDE, DGe for DDcorKirDEDT (sort of)
+ Ge += GiT[i]*expdtgby2; // 2* (..) - Dont multiply by 2 if using lambda and mu in Hooks tensor function
+ DGe += DGiDT[i]*expdtgby2;
+
+ GTsum += GiT[i];
+ DGDTsum += DGiDT[i];
+
+ for (int k=0; k<3; k++){
+ for (int l=0; l<3; l++){
+
+ // Single Convolution
+ // q1->_A[i](k,l) = expdtg*q0->_A[i](k,l) + 2*GiT[i]*expdtgby2*devDE(k,l);
+ q1->_A[i](k,l) = expdtg*q0->_A[i](k,l) + expdtgby2*devDE(k,l);
+
+
+ // devOGammai[i](k,l) = devEe(k,l) - (q1->_A[i](k,l))/GiT[i];
+ // devOGammai[i](k,l) = (q1->_A[i](k,l))/GiT[i];
+ devOGammai[i](k,l) = (q1->_A[i](k,l));
+ }
+ }
+
+ devOi[i] = 2*GiT[i]*q1->_A[i]; // + q1->_C[i];
+ devK += devOi[i];
+ }
+
+
+ // Tr Q_i
+ std::vector<double>& trOi = q1 -> getRefToTrViscoElasticOverStress(); // trOi[i].resize(_N,0.);
+ std::vector<double>& trOGammai = q1 -> getRefToTrViscoElasticOverStrain();
+
+ for (int i=0; i<_Ki.size(); i++){
+
+ double dtk = dt_shift_1/(_ki[i]);
+ double dtk2 = dt_shift_2/(_ki[i]);
+ double expdtk = exp(-dtk);
+ double expdtkby2 = exp(-dtk2);
+ // double expdtkby2 = exp(-dtk/2.);
+
+ getKi(KiT,T_set); getAlphai(AlphaiT,T_set); // Removed 2*_ki[i]
+ DKiDT[i] = -_Ki[i]/(_ki[i]) * Ddt_shiftDT *expdtkby2;
+ DDKiDT[i] = -_Ki[i]/(_ki[i]) * (DDdt_shiftDT - 1/(_ki[i])*pow(Ddt_shiftDT,2)) *expdtkby2;
+
+ // Single Convolution
+ Ke += KiT[i]*expdtkby2;
+ DKe += DKiDT[i]*expdtkby2;
+
+ KTsum += KiT[i]; DKDTsum += DKiDT[i];
+
+ // q1->_B[i] = expdtk*q0->_B[i] + KiT[i]*expdtkby2*trDE;
+ q1->_B[i] = expdtk*q0->_B[i] + expdtkby2*trDE;
+
+
+ trOi[i] = KiT[i]*q1->_B[i]; //+ q1->_D[i];
+
+ // trOGammai
+ // trOGammai[i] = trEe - (q1->_B[i])/KiT[i];
+ // trOGammai[i] = (q1->_B[i])/KiT[i];
+ trOGammai[i] = (q1->_B[i]);
+
+ p += trOi[i];
+ }
+ }
+
+ if (calculateStress){
+ q1->_kirchhoff = devK;
+ q1->_kirchhoff(0,0) += p;
+ q1->_kirchhoff(1,1) += p;
+ q1->_kirchhoff(2,2) += p;
+ }
+
+}
+
+void mlawTVEMullins::isotropicHookTensor(const double K, const double G, STensor43& L) const{
+ double lambda = K - 2.*G/3.;
+ static STensor3 I(1.);
+ 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++){
+ L(i,j,k,l) = lambda*I(i,j)*I(k,l)+ G*(I(i,k)*I(j,l)+I(i,l)*I(j,k));
+ }
+ }
+ }
+ }
+};
+
+
+void mlawTVEMullins::constitutive(const STensor3& F0,
+ const STensor3& F,
+ STensor3& P,
+ const IPVariable *q0i,
+ IPVariable *q1i,
+ STensor43& Tangent,
+ const bool stiff,
+ STensor43* elasticTangent,
+ const bool dTangent,
+ STensor63* dCalgdeps) const{
+
+ const IPTVEMullins *q0=dynamic_cast<const IPTVEMullins *>(q0i);
+ IPTVEMullins *q1 = dynamic_cast<IPTVEMullins *>(q1i);
+
+ double Tc = _Tref;
+ static SVector3 gradT, temp2;
+ static STensor3 temp3;
+ static STensor33 temp33;
+ static double tempVal;
+ static STensor43 dFpdF, dFedF;
+ static STensor3 dFpdT, dFedT, dGammadF;
+ double dGammadT;
+
+ if(elasticTangent!=NULL) Msg::Error("mlawTVEMullins elasticTangent not defined");
+ predictorCorrector_ThermoViscoElastic(F0,F,P,q0,q1,Tangent,dFedF,dFedT,_Tref,_Tref,gradT,gradT,temp2,temp3,temp3,temp2,temp33,tempVal,
+ tempVal,temp3,tempVal,tempVal,temp3,stiff);
+};
+
+void mlawTVEMullins::constitutive(
+ const STensor3& F0, // initial deformation gradient (input @ time n)
+ const STensor3& F, // updated deformation gradient (input @ time n+1)
+ STensor3& P, // updated 1st Piola-Kirchhoff stress tensor (output)
+ const IPVariable *q0i, // array of initial internal variable
+ IPVariable *q1i, // updated array of internal variable (in ipvcur on output),
+ STensor43& Tangent, // mechanical tangents (output)
+ const double& T0, // previous temperature
+ const double& T, // temperature
+ const SVector3& gradT0, // previous temperature gradient
+ const SVector3& gradT, // temperature gradient
+ SVector3& fluxT, // temperature flux
+ STensor3& dPdT, // mechanical-thermal coupling
+ STensor3& dfluxTdgradT, // thermal tengent
+ SVector3& dfluxTdT,
+ STensor33& dfluxTdF, // thermal-mechanical coupling
+ double& thermalSource, // - Cp*DTdt
+ double& dthermalSourcedT, // thermal source
+ STensor3& dthermalSourcedF,
+ double& mechanicalSource, // mechanical source --> convert to heat
+ double& dmechanicalSourcedT,
+ STensor3& dmechanicalSourceF,
+ const bool stiff,
+ STensor43* elasticTangent) const{
+
+ const IPTVEMullins *q0 = dynamic_cast<const IPTVEMullins *>(q0i);
+ IPTVEMullins *q1 = dynamic_cast<IPTVEMullins *>(q1i);
+ if(elasticTangent!=NULL) Msg::Error("mlawTVECaruthers elasticTangent not defined");
+
+ if (!_tangentByPerturbation){
+ static STensor43 dFedF;
+ static STensor3 dFedT;
+ STensorOperation::unity(dFedF);
+ STensorOperation::zero(dFedT);
+ predictorCorrector_ThermoViscoElastic(F0,F,P,q0,q1,Tangent,dFedF,dFedT,T0,T,gradT0,gradT,fluxT,dPdT,dfluxTdgradT,dfluxTdT,dfluxTdF,
+ thermalSource,dthermalSourcedT,dthermalSourcedF,
+ mechanicalSource,dmechanicalSourcedT,dmechanicalSourceF,stiff);
+ }
+ else{
+ predictorCorrector_ThermoViscoElastic(F0,F,P,q0,q1,T0,T,gradT0,gradT,fluxT,thermalSource,mechanicalSource);
+
+ static STensor3 Fplus, Pplus;
+ static SVector3 fluxTPlus, gradTplus;
+ static double thermalSourcePlus;
+ static double mechanicalSourcePlus;
+ static IPTVEMullins qPlus(*q0);
+
+ // perturb F
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ Fplus = (F);
+ Fplus(i,j) += _perturbationfactor;
+ predictorCorrector_ThermoViscoElastic(F0,Fplus,Pplus,q0,&qPlus,T0,T,gradT0,gradT,fluxTPlus,thermalSourcePlus,mechanicalSourcePlus);
+ for (int k=0; k<3; k++){
+ for (int l=0; l<3; l++){
+ Tangent(k,l,i,j) = (Pplus(k,l) - P(k,l))/(_perturbationfactor);
+ }
+ dfluxTdF(k,i,j) = (fluxTPlus(k) - fluxT(k))/_perturbationfactor;
+ }
+ dthermalSourcedF(i,j) = (thermalSourcePlus - thermalSource)/_perturbationfactor;
+ dmechanicalSourceF(i,j) = (mechanicalSourcePlus - mechanicalSource)/_perturbationfactor;
+ }
+ }
+
+ // perturb gradT
+ double gradTpert = _perturbationfactor*T0/1e-3;
+ for (int i=0; i<3; i++){
+ gradTplus = (gradT);
+ gradTplus(i) += (gradTpert);
+ predictorCorrector_ThermoViscoElastic(F0,F,Pplus,q0,&qPlus,T0,T,gradT0,gradTplus,fluxTPlus,thermalSourcePlus,mechanicalSourcePlus);
+ for (int k=0; k<3; k++){
+ dfluxTdgradT(k,i) = (fluxTPlus(k) - fluxT(k))/gradTpert;
+ }
+ }
+
+ // perturb on T
+ double Tplus = T+ _perturbationfactor*T0;
+ predictorCorrector_ThermoViscoElastic(F0,F,Pplus,q0,&qPlus,T0,Tplus,gradT0,gradT,fluxTPlus,thermalSourcePlus,mechanicalSourcePlus);
+ for (int k=0; k<3; k++){
+ for (int l=0; l<3; l++){
+ dPdT(k,l) = (Pplus(k,l) - P(k,l))/(_perturbationfactor*T0);
+ }
+ dfluxTdT(k) = (fluxTPlus(k) - fluxT(k))/(_perturbationfactor*T0);
+ }
+
+ dthermalSourcedT = (thermalSourcePlus - thermalSource)/(_perturbationfactor*T0);
+ dmechanicalSourcedT = (mechanicalSourcePlus - mechanicalSource)/(_perturbationfactor*T0);
+
+ }
+};
+
+void mlawTVEMullins::predictorCorrector_ThermoViscoElastic(
+ const STensor3& F0, // initial deformation gradient (input @ time n)
+ const STensor3& F, // updated deformation gradient (input @ time n+1)
+ STensor3& P, // updated 1st Piola-Kirchhoff stress tensor (output)
+ const IPTVEMullins *q0, // array of initial internal variables
+ IPTVEMullins *q1, // updated array of internal variables (in ipvcur on output),
+ const double& T0, // previous temperature
+ const double& T, // temperature
+ const SVector3& gradT0, // previous temperature gradient
+ const SVector3& gradT, // temperature gradient
+ SVector3& fluxT, // temperature flux
+ double& thermalSource,
+ double& mechanicalSource) const{
+ // temp variables
+ static STensor43 Tangent, dFedF;
+ static STensor3 dPdT, dfluxTdgradT, dthermalSourcedF, dmechanicalSourceF, dFedT;
+ static STensor33 dfluxTdF;
+ static SVector3 dfluxTdT;
+ static double dthermalSourcedT, dmechanicalSourcedT;
+ predictorCorrector_ThermoViscoElastic(F0,F,P,q0,q1,Tangent,dFedF,dFedT,T0,T,gradT0,gradT,fluxT,dPdT,dfluxTdgradT,dfluxTdT,dfluxTdF,
+ thermalSource,dthermalSourcedT,dthermalSourcedF,
+ mechanicalSource,dmechanicalSourcedT,dmechanicalSourceF,false);
+}
+
+void mlawTVEMullins::predictorCorrector_ThermoViscoElastic(
+ const STensor3& F0, // initial deformation gradient (input @ time n)
+ const STensor3& F, // updated deformation gradient (input @ time n+1)
+ STensor3& P, // updated 1st Piola-Kirchhoff stress tensor (output)
+ const IPTVEMullins *q0, // array of initial internal variables
+ IPTVEMullins *q1, // updated array of internal variableS (in ipvcur on output),
+ STensor43& Tangent, // mechanical tangents (output)
+ STensor43& dFedF, // elastic tangent
+ STensor3& dFedT, // elastic tangent
+ const double& T0, // previous temperature
+ const double& T, // temperature
+ const SVector3& gradT0, // previous temperature gradient
+ const SVector3& gradT, // temperature gradient
+ SVector3& fluxT, // temperature flux
+ STensor3& dPdT, // mechanical-thermal coupling
+ STensor3& dfluxTdgradT, // thermal tengent
+ SVector3& dfluxTdT,
+ STensor33& dfluxTdF, // thermal-mechanical coupling
+ double& thermalSource, // - Cp*dTdt
+ double& dthermalSourcedT, // thermal source
+ STensor3& dthermalSourcedF,
+ double& mechanicalSource, // mechanical source --> convert to heat
+ double& dmechanicalSourcedT,
+ STensor3& dmechanicalSourceF,
+ const bool stiff) const{
+
+ // const bool _thermalEstimationPreviousConfig(false);
+
+ // Get time step
+ double dt = this->getTimeStep();
+
+ // Make Identity Tensors
+ static const STensor3 I3(1.);
+ static const STensor43 I4(1.,1.);
+
+ // Calculate/update Log Strain
+ static STensor3 C;
+ static STensor43 dlnCdC;
+ static STensor63 ddlnCddC;
+
+ STensor3& E = q1->getRefToElasticStrain();
+
+ STensorOperation::multSTensor3FirstTranspose(F,F,C);
+ if (_order == 1){
+ STensorOperation::logSTensor3(C,_order,E,&dlnCdC); // as ddlogCddC = 0
+ }
+ else{
+ STensorOperation::logSTensor3(C,_order,E,&dlnCdC,&ddlnCddC);
+ }
+ E *= 0.5; // strain
+
+ // Stresses and Effective Moduli
+ double Ke, Ge = 0.; double Kde, Gde = 0.; double KTsum, GTsum = 0.;
+ double DKe, DGe = 0.; double DKde, DGde = 0.; double DKDTsum, DGDTsum = 0.; // double Wm = 0.;
+
+ ThermoViscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T); // Updates the values of moduli to the current temperatures and calculate stresses.
+ Kde = 0.; Gde = 0.;
+
+ const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
+ static STensor3 secondPK;
+ STensorOperation::multSTensor3STensor43(corKir,dlnCdC,secondPK);
+ STensorOperation::multSTensor3(F,secondPK,P); // 1st PK
+
+ // elastic energy
+ q1->_elasticEnergy = freeEnergyMechanical(*q1,T); // deformationEnergy(*q1,T);
+
+ // thermal energy
+ q1->_thermalEnergy = CpT*T;
+
+
+ // fluxT
+ 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);
+
+
+ // thermalSource
+ if (this->getTimeStep() > 0.){
+ thermalSource = -CpT*(T-T0)/this->getTimeStep();
+ }
+ else
+ thermalSource = 0.;
+
+
+ // mechanicalSource
+ mechanicalSource = 0.;
+ static STensor3 dmechanicalSourceE;
+
+ double T_temp = 0.;
+ if (_testFlag == 1 || _testFlag == 3) {T_temp = T0;} else if(_testFlag == 2) {T_temp = setTemp(T0,T,2);} // CHANGE the latter - case 2 of double convolution
+
+ // CorKir derivatives wrt T
+ static STensor3 DcorKirDT;
+ static STensor3 DDcorKirDT; // for mech source derivatives
+ STensorOperation::zero(DcorKirDT);
+ STensorOperation::zero(DDcorKirDT);
+
+
+ // Add corKirinf terms to DcorKirDT
+ static STensor3 corKirDevInf, devEe;
+ static double corKirTrInf, trEe;
+ STensorOperation::decomposeDevTr(E,devEe,trEe);
+
+ static STensor3 DEe;
+ DEe = E;
+ DEe -= q0->_Ee;
+ q1->_DE = DEe;
+ STensor3& devDE = q1->getRefToDevDE();
+ double& trDE = q1->getRefToTrDE();
+ STensorOperation::decomposeDevTr(DEe,devDE,trDE);
+
+ corKirInfinity( devEe, trEe, T, corKirDevInf, corKirTrInf );
+
+ for (int i=0; i<3; i++){
+ // DcorKirDT(i,i) += corKirTr/3/(Ke+Kde*(T-T_temp))*(DKe + DKde*(T-T_temp)) - 3*KT*(DAlphaDT*(T-_Tinitial) + AlphaT); // 1st formulation
+ DcorKirDT(i,i) += corKirTrInf*DKDT/KT - 3*KT*(DAlphaDT*(T-_Tinitial) + AlphaT);
+ DDcorKirDT(i,i) += corKirTrInf*DDKDT/KT - 6*DKDT*(DAlphaDT*(T-_Tinitial) + AlphaT)- 3*KT*(DDAlphaDT*(T-_Tinitial) + 2*DAlphaDT);
+ for (int j=0; j<3; j++){
+ // DcorKirDT(i,j) += corKirDev(i,j)/(Ge+Gde*(T-T_temp))*(DGe + DGde*(T-T_temp)) ; // 1st formulation
+ DcorKirDT(i,j) += corKirDevInf(i,j)*DGDT/GT;
+ DDcorKirDT(i,j) += corKirDevInf(i,j)*DDGDT/GT;
+ }
+ }
+
+ // Add visco terms
+ if ((_Ki.size() > 0) or (_Gi.size() > 0)){
+
+ // Update the Branch Moduli to the current temperature
+ std::vector<double> KiT, GiT, DKiDT, DGiDT, DDKiDT, DDGiDT;
+ KiT.resize(_N,0.);GiT.resize(_N,0.);DKiDT.resize(_N,0.);DGiDT.resize(_N,0.);DDKiDT.resize(_N,0.);DDGiDT.resize(_N,0.);
+
+ getKi(KiT,T_set,&DKiDT,&DDKiDT); getGi(GiT,T_set,&DGiDT,&DDGiDT);
+
+ static STensor3 DOiDT;
+ static STensor3 DDOiDT;
+ static STensor3 DOGammaiDT;
+
+ static STensor3 DQiDT;
+ static STensor3 DDQiDT;
+ static STensor3 DDGammaiDT;
+
+ STensor3 DDOGammaiDT;
+ STensor3 DOGammai, Oi, temp_WmDQDE, temp_WmDGammaDE;
+ STensor3 DGammai, Qi, temp_WmDqDE, temp_WmDgammaDE;
+
+ const std::vector<STensor3>& devOi = q1->getConstRefToDevViscoElasticOverStress();
+ const std::vector<double>& trOi = q1->getConstRefToTrViscoElasticOverStress();
+ const std::vector<STensor3>& devOi_0 = q0->getConstRefToDevViscoElasticOverStress();
+ const std::vector<double>& trOi_0 = q0->getConstRefToTrViscoElasticOverStress();
+
+ const std::vector<STensor3>& devOGammai = q1->getConstRefToDevViscoElasticOverStrain();
+ const std::vector<double>& trOGammai = q1->getConstRefToTrViscoElasticOverStrain();
+ const std::vector<STensor3>& devOGammai_0 = q0->getConstRefToDevViscoElasticOverStrain();
+ const std::vector<double>& trOGammai_0 = q0->getConstRefToTrViscoElasticOverStrain();
+
+ std::vector<STensor3>& devDOGammaiDT = q1->getRefToDevDOGammaiDT();
+ std::vector<double>& trDOGammaiDT = q1->getRefToTrDOGammaiDT();
+ std::vector<STensor3>& devDOiDT = q1->getRefToDevDOiDT();
+ std::vector<double>& trDOiDT = q1->getRefToTrDOiDT();
+ std::vector<STensor3>& devDDOiDTT = q1->getRefToDevDDOiDTT();
+ std::vector<double>& trDDOiDTT = q1->getRefToTrDDOiDTT();
+
+ double dt_shift, Ddt_shiftDT, DDdt_shiftDT = 0.;
+ shiftedTime(dt,T0,T,&dt_shift,&Ddt_shiftDT,&DDdt_shiftDT); // the convolution term with mid step approximation exponential
+ q1->_dtShift = dt_shift; q1->_DdtShiftDT = Ddt_shiftDT; q1->_DDdtShiftDDT = DDdt_shiftDT;
+
+ double dt_shift_rec, Ddt_shiftDT_rec, DDdt_shiftDT_rec = 0.;
+ shiftedTime(dt,T0,T,&dt_shift_rec,&Ddt_shiftDT_rec,&DDdt_shiftDT_rec,0); // the recursive term exponential
+
+ double dtg,dtk,dtg2,dtk2,expdtk,expdtg,expdtkby2,expdtgby2 = 0.;
+
+ if(stiff){
+ dmechanicalSourcedT = 0;
+ STensorOperation::zero(dmechanicalSourceE);
+ }
+
+ // Construct Oi, DOiDT and DDOiDT tensors
+ for (int i=0; i<_Gi.size(); i++){
+
+ dtk = dt_shift_rec/(_ki[i]); dtg = dt_shift_rec/(_gi[i]); expdtk = exp(-dtk); expdtg = exp(-dtg);
+
+ dtk2 = dt_shift/(_ki[i]); dtg2 = dt_shift/(_gi[i]);
+ expdtkby2 = exp(-dtk2); // exp(-dtk/2.);
+ expdtgby2 = exp(-dtg2); // exp(-dtk/2.);
+
+ STensorOperation::zero(Oi); STensorOperation::zero(DOiDT); STensorOperation::zero(DDOiDT);
+ STensorOperation::zero(DOGammai); STensorOperation::zero(DOGammaiDT); STensorOperation::zero(DDOGammaiDT);
+
+ getKi(KiT,T_set); getGi(GiT,T_set); getAlphai(AlphaiT,T_set);
+ DKiDT[i] = -_Ki[i]/(_ki[i]) * Ddt_shiftDT*expdtkby2; // Removed _ki/2 and _gi/2
+ DDKiDT[i] = -_Ki[i]/(_ki[i]) * (DDdt_shiftDT - 1/(_ki[i])*pow(Ddt_shiftDT,2))*expdtkby2;
+ DGiDT[i] = -_Gi[i]/(_gi[i]) * Ddt_shiftDT*expdtgby2;
+ DDGiDT[i] = -_Gi[i]/(_gi[i]) * (DDdt_shiftDT - 1/(_gi[i])*pow(Ddt_shiftDT,2))*expdtgby2;
+
+ for (int j=0; j<3; j++){
+
+ DOiDT(j,j) += (trOi[i]-trOi_0[i]*expdtk)*(-1/(_ki[i])) * Ddt_shiftDT ;
+ DOiDT(j,j) += trOi_0[i]*expdtk*(-1/(_ki[i])) * Ddt_shiftDT_rec ;
+
+ Oi(j,j) += trOi[i];
+ // Oi(j,j) -= T*DOiDT(j,j);
+
+ DDOiDT(j,j) += (trOi[i]-trOi_0[i]*expdtk)*( (-1/(_ki[i])) * DDdt_shiftDT + (-1/(_ki[i]))*Ddt_shiftDT*(-1/(_ki[i]))*Ddt_shiftDT);
+ DDOiDT(j,j) += trOi_0[i]*expdtk*( (-1/(_ki[i])) * DDdt_shiftDT_rec + (-1/(_ki[i]))*Ddt_shiftDT_rec*(-1/(_ki[i]))*Ddt_shiftDT_rec);
+
+ DOGammai(j,j) += 1/3*trOGammai[i];
+ DOGammai(j,j) -= 1/3*trOGammai_0[i];
+
+ DOGammaiDT(j,j) += 1/3*(trOGammai[i]-trOGammai_0[i]*expdtk)*(-1/(_ki[i])* Ddt_shiftDT);
+ DOGammaiDT(j,j) += 1/3*trOGammai_0[i]*expdtk*(-1/(_ki[i])* Ddt_shiftDT_rec);
+
+ trDOGammaiDT[i] = DOGammaiDT(j,j);
+ DDOGammaiDT(j,j) = DOGammaiDT(j,j) - q0->_trDOGammaiDT[i];
+
+ trDOiDT[i] = DOiDT(j,j);
+ trDDOiDTT[i] = DDOiDT(j,j);
+
+ for (int k=0; k<3; k++){
+
+ DOiDT(j,k) += (devOi[i](j,k)-devOi_0[i](j,k)*expdtg)*(-1/(_gi[i])) * Ddt_shiftDT ;
+ DOiDT(j,k) += devOi_0[i](j,k)*expdtg*(-1/(_gi[i])) * Ddt_shiftDT_rec ;
+
+ Oi(j,k) += devOi[i](j,k);
+ // Oi(j,k) -= T*DOiDT(j,k);
+
+ DDOiDT(j,k) += (devOi[i](j,k)-devOi_0[i](j,k)*expdtg)*( (-1/(_gi[i])) * DDdt_shiftDT + (-1/(_gi[i]))*Ddt_shiftDT*(-1/(_gi[i]))*Ddt_shiftDT );
+ DDOiDT(j,k) += devOi_0[i](j,k)*expdtg*( (-1/(_gi[i])) * DDdt_shiftDT_rec + (-1/(_gi[i]))*Ddt_shiftDT_rec*(-1/(_gi[i]))*Ddt_shiftDT_rec );
+
+ DOGammai(j,k) += devOGammai[i](j,k);
+ DOGammai(j,k) -= devOGammai_0[i](j,k);
+
+ DOGammaiDT(j,k) += (devOGammai[i](j,k)-devOGammai_0[i](j,k)*expdtg)*(-1/(_gi[i])* Ddt_shiftDT);
+ DOGammaiDT(j,k) += devOGammai_0[i](j,k)*expdtg*(-1/(_gi[i])* Ddt_shiftDT_rec);
+
+ devDOGammaiDT[i](j,k) = DOGammaiDT(j,k);
+ DDOGammaiDT(j,k) = DOGammaiDT(j,k) - q0->_devDOGammaiDT[i](j,k);
+
+ devDOiDT[i](j,k) = DOiDT(j,k);
+ devDDOiDTT[i](j,k) = DDOiDT(j,k);
+ }
+ }
+
+ // Corrected Viscous Tensors
+ STensorOperation::zero(Qi); STensorOperation::zero(DQiDT); STensorOperation::zero(DDQiDT);
+ STensorOperation::zero(DGammai); STensorOperation::zero(DDGammaiDT);
+ DQiDT -= DOiDT; DDQiDT -= DDOiDT; DDGammaiDT -= DDOGammaiDT;
+ for (int j=0; j<3; j++){
+ Qi(j,j) += (_Ki[i]*trEe - trOi[i]);
+ for (int k=0; k<3; k++){
+ Qi(j,k) += (2*_Gi[i]*devEe(j,k) - devOi[i](j,k));
+ }
+ }
+ DGammai += DEe;
+ DGammai -= DOGammai;
+ Qi -= (T*DQiDT);
+
+ // corKir Derivatives
+ DcorKirDT += DOiDT;
+ DDcorKirDT += DDOiDT;
+
+ // mechanicalSource - TVE contributions
+
+ static STensor43 temp_DDQiDTDE;
+ static STensor43 temp_DDgammiDtDE;
+
+ double Ki_term = (_Ki[i]*(1-expdtkby2)+T*DKiDT[i]);
+ double Gi_term = (_Gi[i]*(1-expdtgby2)+T*DGiDT[i]);
+ isotropicHookTensor( Ki_term, Gi_term , temp_DDQiDTDE );
+
+ double Ki_term_2 = 1/3*expdtkby2;
+ double Gi_term_2 = 0.5*expdtgby2;
+ isotropicHookTensor( Ki_term_2, Gi_term_2, temp_DDgammiDtDE );
+
+ // For convenience -> Add the TVE Term here
+ STensorOperation::zero(temp_WmDqDE); STensorOperation::zero(temp_WmDgammaDE);
+
+ if (this->getTimeStep() > 0){
+
+ mechanicalSource += (STensorOperation::doubledot(Qi,DGammai)/this->getTimeStep()); // viscoelastic heat (added to mechSource)
+
+ if (stiff){
+
+ // dmechanicalSourceT
+ dmechanicalSourcedT -= (STensorOperation::doubledot(DDQiDT,DGammai))*T/this->getTimeStep(); // 1st Term
+ dmechanicalSourcedT -= (STensorOperation::doubledot(Qi,DOGammaiDT))/this->getTimeStep(); // 2nd Term
+
+ // dmechanicalSourceE
+ for (int l=0; l<3; l++){
+ for (int m=0; m<3; m++){
+ temp_WmDqDE(l,m) = 0.;
+ for (int n=0; n<3; n++){
+ for (int p=0; p<3; p++){
+ temp_WmDqDE(l,m) += temp_DDQiDTDE(l,m,n,p)*DGammai(n,p)/this->getTimeStep();
+ }
+ }
+ }
+ }
+
+ for (int l=0; l<3; l++){
+ for (int m=0; m<3; m++){
+ temp_WmDgammaDE(l,m) = 0.;
+ for (int n=0; n<3; n++){
+ for (int p=0; p<3; p++){
+ temp_WmDgammaDE(l,m) += Qi(n,p)*I4(n,p,l,m)/this->getTimeStep();
+ temp_WmDgammaDE(l,m) -= Qi(n,p)*temp_DDgammiDtDE(n,p,l,m)/this->getTimeStep();
+ }
+ }
+ }
+ }
+
+ dmechanicalSourceE += temp_WmDqDE; // 1st Term
+ dmechanicalSourceE += temp_WmDgammaDE; // 2nd Term
+
+ }
+ }
+ }
+ }
+
+ // put DcorKirDT in the IP
+ STensor3& _DcorKirDT = q1->getRefToDcorKirDT();
+ _DcorKirDT = DcorKirDT;
+
+ // ThermoElastic Heat
+ if (this->getTimeStep() > 0){
+ mechanicalSource += (STensorOperation::doubledot(DcorKirDT,DEe)*T/this->getTimeStep()); // thermoelastic heat (added to mechSource)
+ }
+
+
+ if (stiff){
+
+ // Compute Mechanical Tangents - dSdC and dPdF
+
+ static STensor43 DcorKDE;
+ // isotropicHookTensor(Ke,Ge,DcorKDE); // Hookes tensor
+ isotropicHookTensor(Ke+Kde*(T-T_temp),Ge+Gde*(T-T_temp),DcorKDE); // Reformulate
+
+ static STensor43 DCeDFe; //(0.);
+ static STensor43 DEeDFe(0.);
+ STensorOperation::unity(dFedF);
+ 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++){
+ 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);
+ }
+ }
+ }
+ }
+ }
+
+ STensorOperation::multSTensor43(dlnCdC,DCeDFe,DEeDFe);
+ DEeDFe *= 0.5;
+
+ static STensor43 DsecondPKdC;
+ 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++){
+ DsecondPKdC(i,j,k,l) = 0.;
+ for (int p=0; p<3; p++){
+ for (int q=0; q<3; q++){
+ if (_order != 1){
+ DsecondPKdC(i,j,k,l) += corKir(p,q)*ddlnCddC(p,q,i,j,k,l);
+ }
+ for (int r=0; r<3; r++){
+ for (int s=0; s<3; s++){
+ DsecondPKdC(i,j,k,l) += 0.5*DcorKDE(p,q,r,s)*dlnCdC(p,q,i,j)*dlnCdC(r,s,k,l);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ STensorOperation::zero(Tangent);
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ for (int p=0; p<3; p++){
+ Tangent(i,j,i,p) += secondPK(p,j);
+ for (int q=0; q<3; q++){
+ for (int s=0; s<3; s++){
+ for (int k=0; k<3; k++){
+ Tangent(i,j,p,q) += 2.*F(i,k)*DsecondPKdC(k,j,q,s)*F(p,s);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Compute Mechano-thermal Tangents - dSdT and dPdT
+ static STensor3 DSDT;
+ STensorOperation::zero(DSDT);
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ for (int r=0; r<3; r++){
+ for (int s=0; s<3; s++){
+ DSDT(i,j) += DcorKirDT(r,s)*dlnCdC(r,s,i,j);
+ }
+ }
+ }
+ }
+
+ STensorOperation::zero(dPdT);
+ STensorOperation::multSTensor3(F,DSDT,dPdT);
+
+ // fluxT derivatives
+ dfluxTdT = fluxT;
+ dfluxTdT *= (DKThConDT/KThConT);
+ dfluxTdgradT = Cinv;
+ dfluxTdgradT *= (-KThConT*J);
+ STensorOperation::zero(dfluxTdF);
+
+ if (!_thermalEstimationPreviousConfig){
+
+ static STensor3 DJDF;
+ static STensor43 DCinvDF;
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ DJDF(i,j) = J*Finv(j,i);
+ for (int k=0; k<3; k++){
+ for (int l=0; l<3; l++){
+ DCinvDF(i,j,k,l) = 0.;
+ for (int p=0; p<3; p++){
+ for (int a=0; a<3; a++){
+ for (int b=0; b<3; b++){
+ DCinvDF(i,j,k,l) -= 2.*F(k,p)*Cinv(i,a)*Cinv(j,b)*I4(a,b,p,l);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ for (int k=0; k<3; k++){
+ 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);
+ }
+ }
+ }
+ }
+ }
+
+
+ // thermal source derivatives
+ static STensor3 DCpDF;
+ STensorOperation::zero(DCpDF); // CHANGE for DCpDF
+ if (this->getTimeStep() > 0){
+ dthermalSourcedT = -CpT/this->getTimeStep() - DCpDT*(T-T0)/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();
+ }
+ }
+ }
+ else{
+ dthermalSourcedT = 0.;
+ STensorOperation::zero(dthermalSourcedF);
+ }
+
+
+
+ // mechanical source derivatives --- CHANGED for TVE terms
+ if (_N<1){
+ STensorOperation::zero(dmechanicalSourceE);
+ dmechanicalSourcedT = 0;
+ }
+ STensorOperation::zero(dmechanicalSourceF);
+
+
+ if (this->getTimeStep() > 0){
+
+ // dmechanicalSourceE
+
+ // 1st term
+ static STensor3 CorKir_WmDE; STensorOperation::zero(CorKir_WmDE);
+ static STensor43 DDcorKirDTDE;
+ // isotropicHookTensor(DKe,DGe,DDcorKirDTDE);
+ isotropicHookTensor(DKDTsum,DGDTsum,DDcorKirDTDE);
+ 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++){
+ CorKir_WmDE(i,j) += DDcorKirDTDE(i,j,k,l)*DEe(k,l);
+ }
+ }
+ }
+ }
+
+ // 2nd term
+ static STensor3 DcorKirDT_I;
+ STensorOperation::multSTensor3STensor43(DcorKirDT,I4,DcorKirDT_I);
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ dmechanicalSourceE(i,j) += T*( CorKir_WmDE(i,j) + DcorKirDT_I(i,j) )/this->getTimeStep();
+ }
+ }
+
+ STensorOperation::multSTensor3STensor43(dmechanicalSourceE,DEeDFe,dmechanicalSourceF);
+
+
+ // dmechanicalSourcedT
+ dmechanicalSourcedT += (STensorOperation::doubledot(DcorKirDT,DEe)/this->getTimeStep()); // 1st term
+ dmechanicalSourcedT += T*(STensorOperation::doubledot(DDcorKirDT,DEe)/this->getTimeStep()); // 2nd term
+
+ }
+ }
+};
+
+
+// ========================================================================================================================================
+// Unused Member Functions
+
+double mlawTVEMullins::soundSpeed() const { // Correct this for current temperature
+ return sqrt((_K+4.*_G/3.)/_rho);
+}
diff --git a/NonLinearSolver/materialLaw/mlawTVEMullins.h b/NonLinearSolver/materialLaw/mlawTVEMullins.h
new file mode 100644
index 0000000000000000000000000000000000000000..195b9c7dd20b96ce21d7161ef234e775bd1b3682
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mlawTVEMullins.h
@@ -0,0 +1,291 @@
+//
+// C++ Interface: Material Law
+//
+// Description: Non-Linear Thermo-Visco-Mechanics (Thermo-ViscoElasto-ViscoPlasto Law) with Non-Local Damage Interface (soon.....)
+//
+// Author: <Ujwal Kishore J - FLE_Knight>, (C) 2023;
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#ifndef MLAWTVEMULLINS_H_
+#define MLAWTVEMULLINS_H_
+
+#include "mlaw.h"
+#include "STensor3.h"
+#include "STensor43.h"
+#include "STensor63.h"
+#include "ipTVEMullins.h"
+#include "viscosityLaw.h"
+#include "mlawHyperelasticFailureCriterion.h"
+#include "scalarFunction.h"
+#include "mlawHyperelastic.h"
+
+class mlawTVEMullins : public mlawPowerYieldHyper{
+ protected:
+
+ // Material Properties at Room Temp. - Mechanical
+ double _G; // 2nd Lame parameter (=mu) Shear Modulus
+
+ // Material Properties at Room Temp. - Thermal
+ double _scalarAlpha; // (scalar) Thermal Dilatation Coefficient (=alpha) - Tensor exists for Anisotropic Dilatation
+ double _scalarK; // (scalar) Thermal Conductivity (K) - Tensor exists for Anisotropic Conductivity - Generic Fourier's Law
+ double _Cp; // (scalar) Specific Heat at const pressure or volume or deformation???? - WHAT?
+
+ STensor3 _tensorAlpha; // (tensor3) Thermal Dilatation Coefficient (=alpha)
+ STensor3 _tensorK; // (tensor3) Thermal Conductivity (K)
+
+ // Misc.
+ bool _thermalEstimationPreviousConfig; // flap for thermalEstimation using PreviousConfig
+ int _testFlag; // for various solutions to convolution
+
+ // Temperature Settings and Functions
+ double _Tinitial; // Initial Temperature
+ double _Tref; // Reference Tempereture - Set this regardless of thermo problem, can check the code at this constant temperature.
+ scalarFunction* _temFunc_K; // bulk modulus temperature function
+ scalarFunction* _temFunc_G; // shear modulus temperature function
+ scalarFunction* _temFunc_Alpha; // alpha temperature function
+ scalarFunction* _temFunc_KThCon; // thermal conductivity temperature function
+ 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
+
+ // Shift factor - Everything related to it
+ double _C1; // WLF Universal Constant 1 - set in .py file
+ double _C2; // WLF Universal Constant 2 - set in .py file
+
+
+ #ifndef SWIG
+ 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 ;
+
+ void corKirInfinity(const STensor3& devEe, const double& trEe, const double T0, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const;
+
+ void ThermoViscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
+ const IPTVEMullins *q0, IPTVEMullins *q1,
+ double& Ke, double& Ge, double& Kde, double& Gde,
+ double& DKe, double& DGe, double& DKde, double& DGde,
+ double& KTsum, double& GTsum, double& DKDTsum, double& DGDTsum,
+ const double T0, const double T1, const bool calculateStress = true) const;
+
+ void isotropicHookTensor(const double K, const double G, STensor43& L) const;
+
+ void getMechSourceTVE(const IPTVEMullins *q0, IPTVEMullins *q1, const double& T0, const double& T,
+ const double& DKDTsum, const double& DGDTsum, const STensor43& DEeDFe,
+ double *Wm = NULL, STensor3 *dWmdF = NULL, double *dWmdT = NULL) const;
+
+ void getTVEdCorKirDT(const IPTVEMullins *q0, IPTVEMullins *q1, const double& T0, const double& T) const;
+
+ void predictorCorrector_ThermoViscoElastic( // For thermomech - with _tengentPerturbation
+ const STensor3& F0, // initial deformation gradient (input @ time n)
+ const STensor3& F, // updated deformation gradient (input @ time n+1)
+ STensor3 &P, // updated 1st Piola-Kirchhoff stress tensor (output)
+ const IPTVEMullins *q0, // array of initial internal variables
+ IPTVEMullins *q1, // updated array of internal variables (in ipvcur on output),
+ const double& T0, // previous temperature
+ const double& T, // temperature
+ const SVector3& gradT0, // previous temperature gradient
+ const SVector3& gradT, // temperature gradient
+ SVector3& fluxT, // temperature flux
+ double& thermalSource,
+ double& mechanicalSource) const;
+
+ void predictorCorrector_ThermoViscoElastic( // For thermomech
+ const STensor3& F0, // initial deformation gradient (input @ time n)
+ const STensor3& F, // updated deformation gradient (input @ time n+1)
+ STensor3& P, // updated 1st Piola-Kirchhoff stress tensor (output)
+ const IPTVEMullins *q0, // array of initial internal variables
+ IPTVEMullins *q1, // updated array of internal variableS (in ipvcur on output),
+ STensor43& Tangent, // mechanical tangents (output)
+ STensor43& dFedF, // elastic tangent
+ STensor3& dFedT, // elastic tangent
+ const double& T0, // previous temperature
+ const double& T, // temperature
+ const SVector3 &gradT0, // previous temperature gradient
+ const SVector3 &gradT, // temperature gradient
+ SVector3 &fluxT, // temperature flux
+ STensor3 &dPdT, // mechanical-thermal coupling
+ STensor3 &dfluxTdgradT, // thermal tengent
+ SVector3 &dfluxTdT,
+ STensor33 &dfluxTdF, // thermal-mechanical coupling
+ double& thermalSource, // - Cp*dTdt
+ double& dthermalSourcedT, // thermal source
+ STensor3 &dthermalSourcedF,
+ double& mechanicalSource, // mechanical source --> convert to heat
+ double& dmechanicalSourcedT,
+ STensor3& dmechanicalSourceF,
+ const bool stiff) const;
+ #endif // SWIG
+
+ public:
+ mlawTVEMullins(const int num,const double E,const double nu, const double rho, const double tol, const double Tinitial, const double Alpha, const double KThCon, const double Cp,
+ const bool matrixbyPerturbation = false, const double pert = 1e-8, const bool thermalEstimationPreviousConfig = true);
+
+ // Copy Contructor
+ mlawTVEMullins(const mlawTVEMullins& src);
+
+ // Operator Assignment
+ mlawTVEMullins& operator=(const materialLaw& source);
+
+ // Virtual Destructor
+ virtual ~mlawTVEMullins();
+
+ // Non-Const Functions for strain_order and VE data that are called from .py file. ---- These WILL change the values of the variables defined in the class object.
+ virtual void setViscoElasticNumberOfElement(const int N);
+ virtual void setViscoElasticData(const int i, const double Ei, const double taui);
+ virtual void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
+ virtual void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
+ virtual void setViscoElasticData(const std::string filename);
+
+ // Non-Const Functions for temp and temp_functions that are called from .py file. ---- These WILL change the values of the variables defined in the class object.
+ virtual void setReferenceTemperature(const double Tr){_Tref = Tr;};
+ void setShiftFactorConstantsWLF(const double C1, const double C2){_C1 = C1; _C2 = C2;};
+ void setTemperatureFunction_BulkModulus(const scalarFunction& Tfunc){
+ if (_temFunc_K != NULL) delete _temFunc_K;
+ _temFunc_K = Tfunc.clone();
+ }
+ void setTemperatureFunction_ShearModulus(const scalarFunction& Tfunc){
+ if (_temFunc_G != NULL) delete _temFunc_G;
+ _temFunc_G = Tfunc.clone();
+ }
+ void setTemperatureFunction_ThermalDilationCoefficient(const scalarFunction& Tfunc){
+ if (_temFunc_Alpha != NULL) delete _temFunc_Alpha;
+ _temFunc_Alpha = Tfunc.clone();
+ }
+ void setTemperatureFunction_ThermalConductivity(const scalarFunction& Tfunc){
+ if (_temFunc_KThCon != NULL) delete _temFunc_KThCon;
+ _temFunc_KThCon = Tfunc.clone();
+ }
+ void setTemperatureFunction_SpecificHeat(const scalarFunction& Tfunc){
+ if (_temFunc_Cp != NULL) delete _temFunc_Cp;
+ _temFunc_Cp = Tfunc.clone();
+ }
+
+#ifndef SWIG
+
+ // Const Functions to get temp-updated values of material properties - Args by reference, or return the variable.
+
+ virtual std::vector<double> get_nKi(const double _K, const std::vector<double> _Ki) const{
+ std::vector<double> nKi;
+ nKi.resize(_N,0.);
+ for (int i=0; i<_Ki.size();i++){
+ nKi[i] = _Ki[i]/_K;
+ }
+ return nKi;
+ }
+ virtual void getK(double& KT, const double Tc, double *dK = NULL, double *ddK = NULL) const{
+ KT = _K*_temFunc_K->getVal(Tc);
+ if(dK!=NULL){
+ *dK=_K*_temFunc_K->getDiff(Tc); // Here, *dK can be imagined as *(dK+0) or dK[0] if dK were a vector. See dKi below.
+ }
+ if(ddK!=NULL){
+ *ddK = _K*_temFunc_K->getDoubleDiff(Tc);
+ }
+ }
+ virtual void getG(double& GT, const double Tc, double *dG = NULL, double *ddG = NULL) const{
+ GT = _G*_temFunc_G->getVal(Tc);
+ if(dG!=NULL){
+ *dG = _G*_temFunc_G->getDiff(Tc);
+ }
+ if(ddG!=NULL){
+ *ddG = _G*_temFunc_G->getDoubleDiff(Tc);
+ }
+ }
+ virtual void getAlpha(double& alphaT, const double Tc, double *dAlpha = NULL, double *ddAlpha = NULL) const{
+ alphaT = _scalarAlpha*_temFunc_Alpha->getVal(Tc);
+ if(dAlpha!=NULL){
+ *dAlpha = _scalarAlpha*_temFunc_Alpha->getDiff(Tc);
+ }
+ if(ddAlpha!=NULL){
+ *ddAlpha = _scalarAlpha*_temFunc_Alpha->getDoubleDiff(Tc);
+ }
+ }
+ virtual void getKTHCon(double& KThConT,const double Tc, double *dKThCon = NULL) const{
+ KThConT = _scalarK*_temFunc_KThCon->getVal(Tc);
+ if(dKThCon!=NULL){
+ *dKThCon = _scalarK*_temFunc_KThCon->getDiff(Tc);
+ }
+ }
+ virtual void getCp(double& CpT, const double Tc, double *dCp = NULL, double *ddCp = NULL) const{
+ CpT = _Cp*_temFunc_Cp->getVal(Tc);
+ if(dCp!=NULL){
+ *dCp = _Cp*_temFunc_Cp->getDiff(Tc);
+ }
+ if(ddCp!=NULL){
+ *ddCp = _Cp*_temFunc_Cp->getDoubleDiff(Tc);
+ }
+ }
+
+ // Some more Member Functions
+
+ virtual matname getType() const{return materialLaw::nonlinearTVE;} // function of materialLaw
+ virtual void createIPState(IPStateBase* &ips,const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const{
+ Msg::Error("mlawTVEMullins::createIPState has not been defined");
+ };
+ virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}; // this law is initialized so nothing to do
+
+ // virtual bool withEnergyDissipation() const {if (_N == 0) return false; else return true;};
+ // virtual int getViscoElasticNumberOfElement() const{return _N;};
+
+ virtual double soundSpeed() const; // default but you can redefine it for your case
+
+ // Added from LinearThermomech
+ virtual double getInitialTemperature() const {return _Tinitial;};
+ virtual double getInitialExtraDofStoredEnergyPerUnitField() const {return getExtraDofStoredEnergyPerUnitField(_Tinitial);}
+ virtual double getExtraDofStoredEnergyPerUnitField(double T) const {
+ double Cp; getCp(Cp,T);
+ return Cp;
+ };
+
+ // Const functions to directly export the thermal_cond and alpha when called for ipvcur in dG3DMaterialLaw.cpp
+ virtual const STensor3& getConductivityTensor() const {return _tensorK;}; // isotropic only for now
+ virtual const STensor3& getStiff_alphaDilatation() const {return _tensorAlpha;}; // isotropic only for now
+
+ // Default Constitutive function - placeholder
+ virtual void constitutive(
+ const STensor3& F0, // initial deformation gradient (input @ time n)
+ const STensor3& Fn, // updated deformation gradient (input @ time n+1)
+ STensor3& P, // updated 1st Piola-Kirchhoff stress tensor (output)
+ const IPVariable *q0, // array of initial internal variables
+ IPVariable *q1, // updated array of internal variables (in ipvcur on output),
+ STensor43& Tangent, // constitutive tangents (output)
+ const bool stiff, // if true compute the tangents
+ STensor43* elasticTangent = NULL,
+ const bool dTangent =false,
+ STensor63* dCalgdeps = NULL) const;
+
+ // ThermoMech Constitutive function - J2Full
+ virtual void constitutive(
+ const STensor3& F0, // initial deformation gradient (input @ time n)
+ const STensor3& F1, // updated deformation gradient (input @ time n+1)
+ STensor3& P, // updated 1st Piola-Kirchhoff stress tensor (output)
+ const IPVariable *q0, // array of initial internal variable
+ IPVariable *q1, // updated array of internal variable (in ipvcur on output),
+ STensor43 &Tangent, // constitutive mechanical tangents (output)
+ const double& T0, // previous temperature
+ const double& T, // temperature
+ const SVector3 &gradT0, // previous temperature gradient
+ const SVector3 &gradT, // temperature gradient
+ SVector3 &fluxT, // temperature flux
+ STensor3 &dPdT, // mechanical-thermal coupling
+ STensor3 &dfluxTdgradT, // thermal tengent
+ SVector3 &dfluxTdT,
+ STensor33 &dfluxTdF, // thermal-mechanical coupling
+ double &thermalSource, // - Cp*dTdt
+ double &dthermalSourcedT, // thermal source
+ STensor3 &dthermalSourcedF,
+ double &mechanicalSource, // mechanical source--> convert to heat
+ double &dmechanicalSourcedT,
+ STensor3 & dmechanicalSourceF,
+ const bool stiff,
+ STensor43* elasticTangent = NULL) const;
+
+ virtual materialLaw* clone() const {return new mlawTVEMullins(*this);};
+ virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{}; // do nothing
+ #endif // SWIG
+};
+#endif //mlawTVECaruthers
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 0eb550cf7ffa2c71c697315d8dd003095ce25680..8be63c8954866dd0799ea42d9642fbf4398b759d 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -5435,3 +5435,152 @@ void NonLinearTVPDG3DIPVariable::restart(){
}
// NonLinearTVP Law Interface =================================================================== END
+
+// TVEMullins Law Interface =================================================================== BEGIN
+
+// 1) Constructor Changed
+TVEMullinsDG3DIPVariable::TVEMullinsDG3DIPVariable(const mlawTVEMullins& viscoLaw, double Tinitial, const bool createBodyForceHO, const bool oninter) :
+ ThermoMechanicsDG3DIPVariableBase(createBodyForceHO, oninter){
+ _ipViscoElastic = new IPTVEMullins(viscoLaw.getConstRefToCompressionHardening(),
+ viscoLaw.getConstRefToTractionHardening(),
+ viscoLaw.getConstRefToKinematicHardening(),
+ viscoLaw.getViscoElasticNumberOfElement());
+ this->getRefToTemperature()=Tinitial;
+};
+
+// 2) Copy Constructor Changed
+TVEMullinsDG3DIPVariable::TVEMullinsDG3DIPVariable(const TVEMullinsDG3DIPVariable& src) :
+ ThermoMechanicsDG3DIPVariableBase(src), _ipViscoElastic(src._ipViscoElastic){
+ _ipViscoElastic = NULL;
+ if (src._ipViscoElastic != NULL){
+ _ipViscoElastic = dynamic_cast<IPTVECaruthers*>(src._ipViscoElastic->clone());
+ }
+};
+
+// 3) Operator Assignment Changed
+TVEMullinsDG3DIPVariable& TVEMullinsDG3DIPVariable::operator = (const IPVariable& src){
+ ThermoMechanicsDG3DIPVariableBase::operator =(src);
+ const TVEMullinsDG3DIPVariable* psrc = dynamic_cast<const TVEMullinsDG3DIPVariable*>(&src);
+ if (psrc != NULL){
+ if (psrc->_ipViscoElastic != NULL){
+ if (_ipViscoElastic!= NULL){
+ _ipViscoElastic->operator=(*dynamic_cast<const IPVariable*>(psrc->_ipViscoElastic));
+ }
+ else{
+ _ipViscoElastic = dynamic_cast<IPNonLinearTVE*>(psrc->_ipViscoElastic->clone());
+ }
+ }
+ }
+ return *this;
+};
+
+// 4) Virtual Destructor - Unchanged
+TVEMullinsDG3DIPVariable::~TVEMullinsDG3DIPVariable(){
+ if (_ipViscoElastic){delete _ipViscoElastic; _ipViscoElastic=NULL;};
+};
+
+// 5) Member Function - get(i) - Changed - Refer to HyperVisco...DG3DIP for the Original
+double TVEMullinsDG3DIPVariable::get(const int i) const{
+
+ if (i == IPField::Ee_XX){return _ipViscoElastic->_Ee(0,0);}
+ else if (i == IPField::Ee_XY){return _ipViscoElastic->_Ee(0,1);}
+ else if (i == IPField::Ee_XZ){return _ipViscoElastic->_Ee(0,2);}
+ else if (i == IPField::Ee_YY){return _ipViscoElastic->_Ee(1,1);}
+ else if (i == IPField::Ee_YZ){return _ipViscoElastic->_Ee(1,2);}
+ else if (i == IPField::Ee_ZZ){return _ipViscoElastic->_Ee(2,2);}
+ else if (i == IPField::Eve1_XX){
+ if (_ipViscoElastic->_A.size() > 0)
+ return _ipViscoElastic->_A[0](0,0) + _ipViscoElastic->_B[0]/3.;
+ else
+ return 0.;
+ }
+ else if (i == IPField::Eve1_XY){
+ if (_ipViscoElastic->_A.size() > 0)
+ return _ipViscoElastic->_A[0](0,1);
+ else return 0.;
+ }
+ else if (i == IPField::Eve1_XZ){
+ if (_ipViscoElastic->_A.size() > 0)
+ return _ipViscoElastic->_A[0](0,2);
+ return 0.;
+ }
+ else if (i == IPField::Eve1_YY){
+ if (_ipViscoElastic->_A.size() > 0)
+ return _ipViscoElastic->_A[0](1,1) + _ipViscoElastic->_B[0]/3.;
+ else return 0.;
+ }
+ else if (i == IPField::Eve1_YZ){
+ if (_ipViscoElastic->_A.size() > 0)
+ return _ipViscoElastic->_A[0](1,2);
+ else return 0.;
+ }
+ else if (i == IPField::Eve1_ZZ){
+ if (_ipViscoElastic->_A.size() > 0)
+ return _ipViscoElastic->_A[0](2,2) + _ipViscoElastic->_B[0]/3.;
+ else return 0.;
+ }
+ else if (i == IPField::Eve2_XX){ // Eve2 - WHAT is it?
+ if (_ipViscoElastic->_A.size() > 1)
+ return _ipViscoElastic->_A[1](0,0) + _ipViscoElastic->_B[1]/3.;
+ else return 0.;
+ }
+ else if (i == IPField::Eve2_XY){
+ if (_ipViscoElastic->_A.size() > 1)
+ return _ipViscoElastic->_A[1](0,1);
+ else return 0.;
+ }
+ else if (i == IPField::Eve2_XZ){
+ if (_ipViscoElastic->_A.size() > 1)
+ return _ipViscoElastic->_A[1](0,2);
+ else return 0.;
+ }
+ else if (i == IPField::Eve2_YY){
+ if (_ipViscoElastic->_A.size() > 1)
+ return _ipViscoElastic->_A[1](1,1) + _ipViscoElastic->_B[1]/3.;
+ else return 0.;
+ }
+ else if (i == IPField::Eve2_YZ){
+ if (_ipViscoElastic->_A.size() > 1)
+ return _ipViscoElastic->_A[1](1,2);
+ else return 0.;
+ }
+ else if (i == IPField::Eve2_ZZ){
+ if (_ipViscoElastic->_A.size() > 1)
+ return _ipViscoElastic->_A[1](2,2) + _ipViscoElastic->_B[1]/3.;
+ else return 0.;
+ }
+
+ // Additional IPFields from ThermoMechanicsDG3DIPVariableBase - Adopted here directly
+ else if (i == IPField::TEMPERATURE){return getConstRefToField(0);}
+ else if (i == IPField::THERMALFLUX_X){return getConstRefToFlux()[0](0);}
+ else if (i == IPField::THERMALFLUX_Y){return getConstRefToFlux()[0](1);}
+ else if (i == IPField::THERMALFLUX_Z){return getConstRefToFlux()[0](2);}
+ else if (i == IPField::FIELD_SOURCE){return getConstRefToFieldSource()(0);}
+ else if (i == IPField::MECHANICAL_SOURCE){return getConstRefToMechanicalSource()(0);}
+
+ else
+ {
+ return dG3DIPVariable::get(i);
+ }
+};
+
+// 5) Member Function - defoEnergy() - Unchanged
+double TVEMullinsDG3DIPVariable::defoEnergy() const{
+ return _ipViscoElastic->defoEnergy();
+}
+
+double TVEMullinsDG3DIPVariable::plasticEnergy() const{
+ return getIPTVECaruthers()->plasticEnergy();
+}
+
+// 7) Other Member Functions - damageEnergy(), etc, may have to be added later.
+
+// 8) Member Function - restart() - Changed only 2nd line
+void TVEMullinsDG3DIPVariable::restart(){
+// dG3DIPVariable::restart();
+ ThermoMechanicsDG3DIPVariableBase::restart();
+ restartManager::restart(_ipViscoElastic);
+ return;
+}
+
+// TVEMullins Law Interface =================================================================== END
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index 016efa98814b6fc3917c75e598fb7a560cd2ddc0..b16f5aadea83b08f8184d2145112071018b301e9 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -67,6 +67,8 @@
#include "mlawNonLinearTVE.h"
#include "ipNonLinearTVP.h"
#include "mlawNonLinearTVP.h"
+#include "ipTVEMullins.h"
+#include "mlawTVEMullins.h"
#if defined(HAVE_TORCH)
#include <torch/torch.h>
@@ -5180,4 +5182,56 @@ class NonLinearTVPDG3DIPVariable : public ThermoMechanicsDG3DIPVariableBase{ /
// NonLinearTVP Law Interface =================================================================== END
-#endif
+// TVECaruthers Law Interface =================================================================== BEGIN
+
+class TVEMullinsDG3DIPVariable : public ThermoMechanicsDG3DIPVariableBase{ // Changed
+
+ protected:
+ IPTVEMullins * _ipViscoElastic;
+ TVEMullinsDG3DIPVariable(const bool createBodyForceHO, const bool oninter):ThermoMechanicsDG3DIPVariableBase(createBodyForceHO, oninter) {
+ Msg::Error("TVEMullinssDG3DIPVariable: Cannot use the default constructor of ThermoMechanicsDG3DIPVariableBase");}
+
+ public:
+ TVEMullinsDG3DIPVariable(const mlawTVEMullins& viscoLaw, double Tinitial, const bool createBodyForceHO, const bool oninter); // Changed
+ TVEMullinsDG3DIPVariable(const TVEMullinsDG3DIPVariable& src);
+ virtual TVEMullinsDG3DIPVariable& operator = (const IPVariable& src);
+ virtual ~TVEMullinsDG3DIPVariable();
+
+ virtual void setLocation(const IPVariable::LOCATION loc) {
+ dG3DIPVariable::setLocation(loc);
+ if (_ipViscoElastic)
+ _ipViscoElastic->setLocation(loc);
+ };
+
+ virtual IPVariable* getInternalData() {return _ipViscoElastic;};
+ virtual const IPVariable* getInternalData() const {return _ipViscoElastic;};
+
+ // Specific functions - getIP
+ IPTVEMullins* getIPTVEMullins() {return _ipViscoElastic;};
+ const IPTVEMullins* getIPTVEMullins() const {return _ipViscoElastic;};
+
+ virtual double get(const int i) const;
+ virtual double defoEnergy() const;
+
+// Added from LinearThermoMechanicsDG3DIPVariable
+ virtual double getInternalEnergyExtraDofDiffusion() const {return _ipViscoElastic->getThermalEnergy();}; // Added
+ virtual double plasticEnergy() const; // Added
+ virtual double damageEnergy() const {return _ipViscoElastic->damageEnergy();}; // Added
+ // virtual int fractureEnergy(double* arrayEnergy) const{arrayEnergy[0]=_linearTMIP.getConstRefToFractureEnergy(); return 1;} // need it to get H1 norm in the interface e.g. LinearThermoMech
+ // Fracture energy for WHAT?? Added but unused?
+
+// for Path-Following based on irreversible energy
+ virtual double irreversibleEnergy() const {return _ipViscoElastic->irreversibleEnergy();};
+ virtual double& getRefToIrreversibleEnergy() {return _ipViscoElastic->getRefToIrreversibleEnergy();};
+
+ virtual const STensor3& getConstRefToDIrreversibleEnergyDDeformationGradient() const {return _ipViscoElastic->getConstRefToDIrreversibleEnergyDF();};
+ virtual STensor3& getRefToDIrreversibleEnergyDDeformationGradient() {return _ipViscoElastic->getRefToDIrreversibleEnergyDF();};
+
+ virtual IPVariable* clone() const {return new TVEMullinsDG3DIPVariable(*this);};
+ virtual void restart();
+};
+
+// TVEMullins Law Interface =================================================================== END
+
+#endif //DG3DIPVARIABLE
+
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index f6e6feb02ad0ac2500bc13b7e3dedd4bb4bab07c..8a080ad7c0571a74c0b7db8d285d2c5864f4b0cb 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -8720,3 +8720,175 @@ void NonLinearTVPDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp
}
// NonLinearTVP Law Interface =================================================================== END
+
+
+// TVEMullins Law Interface =================================================================== BEGIN
+
+TVEMullinsDG3DMaterialLaw::TVEMullinsDG3DMaterialLaw(const int num, const double rho, const double E, const double nu,
+ const double Tinitial, const double Alpha, const double KThCon, const double Cp,
+ const bool matrixbyPerturbation, const double pert, const bool thermalEstimationPreviousConfig)
+ : dG3DMaterialLaw(num,rho), _viscoLaw(num,E,nu,rho,0.,Tinitial,Alpha, KThCon, Cp, matrixbyPerturbation,pert, thermalEstimationPreviousConfig)
+{
+ fillElasticStiffness(E, nu, elasticStiffness);
+};
+
+TVEMullinsDG3DMaterialLaw::TVEMullinsDG3DMaterialLaw(const TVEMullinsDG3DMaterialLaw& src): dG3DMaterialLaw(src), _viscoLaw(src._viscoLaw){};
+
+void TVEMullinsDG3DMaterialLaw::setStrainOrder(const int order){
+ _viscoLaw.setStrainOrder(order);
+};
+void TVEMullinsDG3DMaterialLaw::setViscoelasticMethod(const int method){
+ _viscoLaw.setViscoelasticMethod(method);
+};
+void TVEMullinsDG3DMaterialLaw::setViscoElasticNumberOfElement(const int N){
+ _viscoLaw.setViscoElasticNumberOfElement(N);
+};
+void TVEMullinsDG3DMaterialLaw::setViscoElasticData(const int i, const double Ei, const double taui){
+ _viscoLaw.setViscoElasticData(i,Ei,taui);
+};
+void TVEMullinsDG3DMaterialLaw::setViscoElasticData_Bulk(const int i, const double Ki, const double ki){
+ _viscoLaw.setViscoElasticData_Bulk(i,Ki,ki);
+};
+void TVEMullinsDG3DMaterialLaw::setViscoElasticData_Shear(const int i, const double Gi, const double gi){
+ _viscoLaw.setViscoElasticData_Shear(i,Gi,gi);
+};
+void TVEMullinsDG3DMaterialLaw::setViscoElasticData(const std::string filename){
+ _viscoLaw.setViscoElasticData(filename);
+};
+
+void TVEMullinsDG3DMaterialLaw::setReferenceTemperature(const double Tref){
+ _viscoLaw.setReferenceTemperature(Tref);
+};
+void TVEMullinsDG3DMaterialLaw::setShiftFactorConstantsWLF(const double C1, const double C2){
+ _viscoLaw.setShiftFactorConstantsWLF(C1,C2);
+};
+
+// Added extra argument in the NonLinearTVEDG3DIPVariable contructor - Tinitial (added from LinearThermoMech)
+void TVEMullinsDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+ // check interface element
+ bool inter=true;
+ const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+ if(iele==NULL) inter=false;
+ IPVariable* ipvi = new TVEMullinsDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),hasBodyForce,inter);
+ IPVariable* ipv1 = new TVEMullinsDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),hasBodyForce,inter);
+ IPVariable* ipv2 = new TVEMullinsDG3DIPVariable(_viscoLaw,_viscoLaw.getInitialTemperature(),hasBodyForce,inter);
+ if(ips != NULL) delete ips;
+ ips = new IP3State(state_,ipvi,ipv1,ipv2);
+}
+
+void TVEMullinsDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+ if(ipv !=NULL) delete ipv;
+ bool inter=true;
+ const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+ if(iele == NULL){
+ inter=false;
+ }
+ ipv = new TVEMullinsDG3DMaterialLaw(_viscoLaw,_viscoLaw.getInitialTemperature(),hasBodyForce,inter);
+}
+
+double TVEMullinsDG3DMaterialLaw::getExtraDofStoredEnergyPerUnitField(double T) const{
+ return _viscoLaw.getExtraDofStoredEnergyPerUnitField(T);
+}
+
+double TVEMullinsDG3DMaterialLaw::getInitialExtraDofStoredEnergyPerUnitField() const{
+ return _viscoLaw.getInitialExtraDofStoredEnergyPerUnitField();
+};
+
+void TVEMullinsDG3DMaterialLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const{
+ // get ipvariable
+ TVEMullinsDG3DIPVariable* ipvcur;
+ const TVEMullinsDG3DIPVariable* ipvprev;
+ FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
+ if(ipvtmp !=NULL)
+ {
+ ipvcur = static_cast<TVEMullinsDG3DIPVariable*>(ipvtmp->getIPvBulk());
+ const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
+ ipvprev = static_cast<const TVEMullinsDG3DIPVariable*>(ipvtmp2->getIPvBulk());
+ }
+ else
+ {
+ ipvcur = static_cast<TVEMullinsDG3DIPVariable*>(ipv);
+ ipvprev = static_cast<const TVEMullinsDG3DIPVariable*>(ipvp);
+ }
+ _viscoLaw.checkInternalState(ipvcur->getInternalData(), ipvprev->getInternalData());
+};
+
+void TVEMullinsDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac, const bool dTangent){
+
+ // get ipvariable
+ TVEMullinsDG3DIPVariable* ipvcur;
+ const TVEMullinsDG3DIPVariable* ipvprev;
+ FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
+ if(ipvtmp !=NULL)
+ {
+ ipvcur = static_cast<TVEMullinsDG3DIPVariable*>(ipvtmp->getIPvBulk());
+ const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
+ ipvprev = static_cast<const TVEMullinsDG3DIPVariable*>(ipvtmp2->getIPvBulk());
+ }
+ else
+ {
+ ipvcur = static_cast<TVEMullinsDG3DIPVariable*>(ipv);
+ ipvprev = static_cast<const TVEMullinsDG3DIPVariable*>(ipvp);
+ }
+
+// Added a bunch of ip variables and derivatives from LinearThermoMechanicsDG3DMaterialLaw in dG3DMaterialLaw.cpp
+ // Current and Previous Deformation Gradient (F) and Temperature (T), Current Temp Gradient (gradT) and other current derivatives.
+ // THe new variables need to be added in the constitutive function of material law.
+ const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
+ const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
+ const double T = ipvcur->getConstRefToTemperature();
+ const double T0 = ipvprev->getConstRefToTemperature();
+ const SVector3& gradT0 = ipvprev->getConstRefToGradT();
+ const SVector3& gradT = ipvcur->getConstRefToGradT();
+ STensor3& dPdT = ipvcur->getRefTodPdT();
+ SVector3& fluxT = ipvcur->getRefToThermalFlux();
+ STensor3& dqdgradT = ipvcur->getRefTodThermalFluxdGradT();
+ SVector3& dqdT = ipvcur->getRefTodThermalFluxdT();
+ STensor33& dqdF = ipvcur->getRefTodThermalFluxdF();
+//const STensor3 *linearK()=ipvcur-> getConstRefTolinearK();
+
+ STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress();
+ STensor43& Tangent = ipvcur->getRefToTangentModuli();
+ const double& w0 = ipvprev->getConstRefToThermalSource();
+ double& w = ipvcur->getRefToThermalSource();
+ double& dwdt = ipvcur->getRefTodThermalSourcedField();
+ STensor3& dwdf = ipvcur->getRefTodThermalSourcedF();
+ const SVector3& N = ipvcur->getRefToReferenceOutwardNormal();
+ const SVector3& ujump = ipvcur->getRefToJump()(0);
+ const double& Tjump = ipvcur->getRefToFieldJump()(0);
+ double& Cp = ipvcur->getRefToExtraDofFieldCapacityPerUnitField()(0);
+ const double& mechSource0 = ipvprev->getConstRefToMechanicalSource()(0);
+ double& mechSource = ipvcur->getRefToMechanicalSource()(0);
+ double& dmechSourcedt = ipvcur->getRefTodMechanicalSourcedField()(0,0);
+ STensor3& dmechSourcedf = ipvcur->getRefTodMechanicalSourcedF()[0];
+
+ // data for J2 law --------------------------------------------- WHAT is J2law here?
+ IPTVEMullins* q1 = ipvcur->getTVEMullins();
+ const IPTVEMullins* q0 = ipvprev->getTVEMullins();
+ // static STensor63 dCalgdeps;
+ // STensor43& elasticL = ipvcur->getRefToElasticTangentModuli();
+
+ // Compute Stress
+ _viscoLaw.constitutive(F0,F1,P,q0,q1,Tangent,T0,T,gradT0,gradT,fluxT,dPdT,dqdgradT,dqdT,dqdF,w,dwdt,dwdf,mechSource,dmechSourcedt,dmechSourcedf,stiff);
+
+ ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
+ ipvcur->setRefToLinearConductivity(_viscoLaw.getConductivityTensor());
+ ipvcur->setRefToStiff_alphaDilatation(_viscoLaw.getStiff_alphaDilatation());
+
+// Add the following lines here from LinearThermoMechanicsDG3DMaterialLaw ??
+ // Fracture hasn't been defined as of yet.
+/*
+ if (checkfrac==true) // to get the H1 norm in the interface
+ q1->_fracEnergy=defoEnergy(ujump,Tjump,N);
+ else
+ q1->_elasticEnergy=defoEnergy(F1, gradT); */
+ Cp = _viscoLaw.getExtraDofStoredEnergyPerUnitField(T);
+
+}
+
+// defoEnergy functions from LinearThermoMechanicsDG3DMaterialLaw - add them?
+
+// TVEMullins Law Interface =================================================================== END
+
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 4e0e4a9ec9d229b933c3dc6f1c81607356c9ee75..6b2fd2919534ab24b01f4a25e5c457e53528895f 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -52,6 +52,7 @@
#include "mlawNonLinearTVM.h"
#include "mlawNonLinearTVE.h"
#include "mlawNonLinearTVP.h"
+#include "mlawTVEMullins.h"
#endif
#if defined(HAVE_TORCH)
#include <torch/torch.h>
@@ -3048,4 +3049,79 @@ class NonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{ // public Material
// NonLinearTVP Law Interface =================================================================== END
+// TVEMullins Law Interface =================================================================== BEGIN
+
+class TVEMullinsDG3DMaterialLaw : public dG3DMaterialLaw{ // public MaterialLaw
+ #ifndef SWIG
+ protected:
+ mlawTVEMullins _viscoLaw;
+
+ #endif //SWIG
+ public:
+ TVEMullinsDG3DMaterialLaw(const int num, const double rho, const double E, const double nu,
+ const double Tinitial, const double Alpha, const double KThCon, const double Cp,
+ const bool matrixbyPerturbation = false, const double pert = 1e-8, const bool thermalEstimationPreviousConfig = false);
+
+ void setStrainOrder(const int order);
+ void setViscoelasticMethod(const int method);
+ void setViscoElasticNumberOfElement(const int N);
+ void setViscoElasticData(const int i, const double Ei, const double taui);
+ void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
+ void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
+ void setViscoElasticData(const std::string filename);
+
+ void setReferenceTemperature(const double Tref);
+ void setTestOption(const int testFlag);
+ void setShiftFactorConstantsWLF(const double C1, const double C2);
+ void setTemperatureFunction_BulkModulus(const scalarFunction& Tfunc);
+ void setTemperatureFunction_ShearModulus(const scalarFunction& Tfunc);
+ void setTemperatureFunction_ThermalDilationCoefficient(const scalarFunction& Tfunc);
+ void setTemperatureFunction_ThermalConductivity(const scalarFunction& Tfunc);
+ void setTemperatureFunction_SpecificHeat(const scalarFunction& Tfunc);
+
+ #ifndef SWIG
+ TVEMullinsDG3DMaterialLaw(const TVEMullinsDG3DMaterialLaw &source);
+ virtual ~TVEMullinsDG3DMaterialLaw(){}
+ // set the time of _j2law
+ virtual void setTime(const double t,const double dtime){
+ dG3DMaterialLaw::setTime(t,dtime);
+ _viscoLaw.setTime(t,dtime);
+ }
+ virtual materialLaw::matname getType() const {return _viscoLaw.getType();}
+ virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const;
+ virtual void createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const;
+ virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}
+ virtual double soundSpeed() const{return _viscoLaw.soundSpeed();} // or change value ??
+ virtual void stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff=true, const bool checkfrac=true, const bool dTangent=false);
+ virtual const double bulkModulus() const {return _viscoLaw.bulkModulus();}
+ virtual double scaleFactor() const{return _viscoLaw.scaleFactor();}
+
+// Added from LinearThermoMechanicsDG3DMaterialLaw
+ virtual int getNumberOfExtraDofsDiffusion() const {return 1;}
+ virtual double getExtraDofStoredEnergyPerUnitField(double T) const;
+ virtual double getInitialExtraDofStoredEnergyPerUnitField() const;
+ // virtual double defoEnergy(const STensor3& Fn,const SVector3& gradT) const; // NO definition in .cpp yet
+ // virtual double defoEnergy( const SVector3&ujump,const double &Tjump, const SVector3& N) const; // NO definition in .cpp yet
+//
+ virtual materialLaw* clone() const{return new TVEMullinsDG3DMaterialLaw(*this);};
+ virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const;
+ virtual bool withEnergyDissipation() const {return _viscoLaw.withEnergyDissipation();};
+ virtual void setMacroSolver(const nonLinearMechSolver* sv){
+ dG3DMaterialLaw::setMacroSolver(sv);
+ _viscoLaw.setMacroSolver(sv);
+ };
+
+ virtual const materialLaw *getConstNonLinearSolverMaterialLaw() const
+ {
+ return &_viscoLaw;
+ }
+ virtual materialLaw *getNonLinearSolverMaterialLaw()
+ {
+ return &_viscoLaw;
+ }
+ #endif //SWIG
+};
+
+// TVEMullins Law Interface =================================================================== END
+
#endif