diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp b/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp
index 3b9e0ae670bb61dc0824d24e4a26401cff828c15..ea98eac368ec0303d9001656c24239b4033daad3 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp
@@ -15,10 +15,14 @@
IPNonLinearTVE::IPNonLinearTVE(const J2IsotropicHardening* comp,
const J2IsotropicHardening* trac,
- const kinematicHardening* kin, const int N): IPHyperViscoElastoPlastic(comp,trac,kin,N), // IPHyperViscoElastic(N),
+ const kinematicHardening* kin, const int N, const mullinsEffect* mullins): IPHyperViscoElastoPlastic(comp,trac,kin,N), // IPHyperViscoElastic(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.){
+
+ _ipMullinsEffect = NULL;
+ if (mullins!= NULL) mullins->createIPVariable(_ipMullinsEffect);
+
_devOGammai.clear();
_trOGammai.clear();
_devDOGammaiDT.clear();
@@ -63,7 +67,15 @@ IPNonLinearTVE::IPNonLinearTVE(const IPNonLinearTVE& src): IPHyperViscoElastoPla
_mechSrcTVE(src._mechSrcTVE),_DmechSrcTVEdT(src._DmechSrcTVEdT),_DmechSrcTVEdF(src._DmechSrcTVEdF),
_C(src._C),_D(src._D), _E(src._E), _F(src._F), _G(src._G),
_thermalEnergy(src._thermalEnergy),
- _dtShift(src._dtShift),_DdtShiftDT(src._DdtShiftDT),_DDdtShiftDDT(src._DDdtShiftDDT){};
+ _dtShift(src._dtShift),_DdtShiftDT(src._DdtShiftDT),_DDdtShiftDDT(src._DDdtShiftDDT){
+
+
+ if (src._ipMullinsEffect != NULL)
+ _ipMullinsEffect = dynamic_cast<IPMullinsEffect*>(src._ipMullinsEffect->clone());
+ else
+ _ipMullinsEffect = NULL;
+
+};
IPNonLinearTVE& IPNonLinearTVE::operator=(const IPVariable& src)
{
@@ -100,10 +112,25 @@ IPNonLinearTVE& IPNonLinearTVE::operator=(const IPVariable& src)
_dtShift = psrc->_dtShift;
_DdtShiftDT = psrc->_DdtShiftDT;
_DDdtShiftDDT = psrc->_DDdtShiftDDT;
+
+ if ( psrc->_ipMullinsEffect != NULL) {
+ if (_ipMullinsEffect == NULL){
+ _ipMullinsEffect = dynamic_cast<IPMullinsEffect*>(psrc->_ipMullinsEffect->clone());
+ }
+ else{
+ _ipMullinsEffect->operator=(*dynamic_cast<const IPMullinsEffect*>(psrc->_ipMullinsEffect));
+ }
+ }
+
}
return *this;
}
+IPNonLinearTVE::~IPNonLinearTVE(){
+ if (_ipMullinsEffect != NULL) delete _ipMullinsEffect;
+ _ipMullinsEffect = NULL;
+};
+
double IPNonLinearTVE::defoEnergy() const
{
return IPHyperViscoElastic::defoEnergy();
@@ -148,4 +175,18 @@ void IPNonLinearTVE::restart(){
restartManager::restart(_dtShift);
restartManager::restart(_DdtShiftDT);
restartManager::restart(_DDdtShiftDDT);
-}
\ No newline at end of file
+
+ if (_ipMullinsEffect != NULL)
+ restartManager::restart(_ipMullinsEffect);
+}
+
+const IPMullinsEffect& IPNonLinearTVE::getConstRefToIPMullinsEffect() const{
+ if(_ipMullinsEffect==NULL)
+ Msg::Error("IPHyperViscoElastoPlastic: _ipMullinsEffect not initialized");
+ return *_ipMullinsEffect;
+};
+IPMullinsEffect& IPNonLinearTVE::getRefToIPMullinsEffect(){
+ if(_ipMullinsEffect==NULL)
+ Msg::Error("IPHyperViscoElastoPlastic: _ipMullinsEffect not initialized");
+ return *_ipMullinsEffect;
+};
\ No newline at end of file
diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVE.h b/NonLinearSolver/internalPoints/ipNonLinearTVE.h
index c6e6d68223f0ba862d0b726fac748d0acc7af394..c2e331df6e6fbc4b5cf6c817b70464d757bb20c5 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVE.h
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVE.h
@@ -16,7 +16,7 @@
#include "kinematicHardening.h"
#include "STensor3.h"
#include "STensor43.h"
-
+#include "mullinsEffect.h"
// class IPNonLinearTVE : public IPHyperViscoElastic{
class IPNonLinearTVE : public IPHyperViscoElastoPlastic{
@@ -60,13 +60,16 @@ class IPNonLinearTVE : public IPHyperViscoElastoPlastic{
double _DdtShiftDT;
double _DDdtShiftDDT;
+ // Mullins Effect IP
+ IPMullinsEffect* _ipMullinsEffect;
+
public:
IPNonLinearTVE(const J2IsotropicHardening* comp,
const J2IsotropicHardening* trac,
- const kinematicHardening* kin, const int N);
+ const kinematicHardening* kin, const int N, const mullinsEffect* mullins);
IPNonLinearTVE(const IPNonLinearTVE& src);
virtual IPNonLinearTVE& operator=(const IPVariable& src);
- virtual ~IPNonLinearTVE(){};
+ virtual ~IPNonLinearTVE();
virtual std::vector<STensor3>& getRefToDevViscoElasticOverStrain() {return _devOGammai;};
virtual const std::vector<STensor3>& getConstRefToDevViscoElasticOverStrain() const {return _devOGammai;};
@@ -126,6 +129,9 @@ class IPNonLinearTVE : public IPHyperViscoElastoPlastic{
virtual double& getShiftedTimeStepTempDoubleDerivative() {return _DDdtShiftDDT;};
virtual const double& getConstShiftedTimeStepTempDoubleDerivative() const {return _DDdtShiftDDT;};
+ virtual const IPMullinsEffect& getConstRefToIPMullinsEffect() const;
+ virtual IPMullinsEffect& getRefToIPMullinsEffect();
+
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)
diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVP.cpp b/NonLinearSolver/internalPoints/ipNonLinearTVP.cpp
index 3806a37177e51636a5d13e5a81f0ffc04a67efb8..43c98ae52bac099b08ec2ec20d4f5bd17ce0ec0f 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVP.cpp
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVP.cpp
@@ -12,7 +12,7 @@
#include "ipNonLinearTVP.h"
IPNonLinearTVP::IPNonLinearTVP(const J2IsotropicHardening* comp, const J2IsotropicHardening* trac,
- const kinematicHardening* kin, const int N):IPNonLinearTVE(comp,trac,kin,N),
+ const kinematicHardening* kin, const int N, const mullinsEffect* mullins):IPNonLinearTVE(comp,trac,kin,N,mullins),
_DgammaDT(0.),_DGammaDT(0.), _DGammaDF(0.), _DirreversibleEnergyDT(0.),
_psiTVM(0.), _DpsiTVMdT(0.), _DDpsiTVMdTT(0.),
_ModMandel(0.),_DModMandelDT(0.), _DbackSigDT(0.), _DModMandelDF(0.), _DbackSigDF(0.), _DphiPDF(0.), _GammaN(0.), _dGammaNdT(0.), _dGammaNdF(0.),
diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVP.h b/NonLinearSolver/internalPoints/ipNonLinearTVP.h
index 2389f05c3b77497ced36d597cdf561a3ce186942..56e806bb0bc739facd0cfa1bd2e86f77ef8b3b16 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVP.h
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVP.h
@@ -57,7 +57,7 @@ class IPNonLinearTVP : public IPNonLinearTVE{
public:
IPNonLinearTVP(const J2IsotropicHardening* comp,
const J2IsotropicHardening* trac,
- const kinematicHardening* kin, const int N);
+ const kinematicHardening* kin, const int N, const mullinsEffect* mullins);
IPNonLinearTVP(const IPNonLinearTVP& src);
virtual IPNonLinearTVP& operator =(const IPVariable &source);
virtual ~IPNonLinearTVP();
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
index ca9ba845f2927e5a77f3762ad2110c97756ac202..0de0bd8dae1b751e3ffebf8d76352273b86649cf 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
@@ -1024,10 +1024,11 @@ void mlawNonLinearTVE::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
DKe = DKDT; DGe = DGDT; DKDTsum = DKDT; DGDTsum = DGDT;
// Kelvin-Voight
+ /*
double invGe = 1./GT;
double invKe = 1./KT;
STensor3 D(0.);
- double V= 0.;
+ double V= 0.;*/
// Main Loop for Maxwell Branches
if ((_Ki.size() > 0) or (_Gi.size() > 0)){
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.h b/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
index c8eed1b2bd409aed251beddcdc1018920e736bc3..c009c283110e24201207cd411fcd75bd9513aaaa 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
@@ -81,6 +81,9 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
// 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
+
+ // Mullins Effect
+ mullinsEffect* _mullinsEffect;
#ifndef SWIG
protected:
@@ -93,6 +96,7 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
virtual double freeEnergyMechanical(const IPNonLinearTVE& q, const double T) const ;
virtual double setTemp(const double T0, const double T1, const int para) const ;
+
void corKirInfinity(const IPNonLinearTVE *q1, const STensor3& devEe, const double& trEe, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const;
void ThermoViscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
const IPNonLinearTVE *q0, IPNonLinearTVE *q1,
@@ -164,6 +168,11 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
// Virtual Destructor
virtual ~mlawNonLinearTVE();
+ void setMullinsEffect(const mullinsEffect& mullins){
+ if (_mullinsEffect) delete _mullinsEffect;
+ _mullinsEffect = mullins.clone();
+ };
+
// 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);
@@ -222,6 +231,8 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
// Const Functions to get temp-updated values of material properties - Args by reference, or return the variable.
+ virtual mullinsEffect* getConstRefToMullinsEffect() const {return _mullinsEffect;};
+
virtual std::vector<double> getBetai(const double _K, const std::vector<double> _Ki) const{
std::vector<double> Beta_i;
Beta_i.resize(_N,0.);
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 8ce50726a688ff37116e673383591de6f96ec008..88095add61d8a9d8d6859662ac0db7bf1bac5dc8 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -5112,7 +5112,8 @@ NonLinearTVEDG3DIPVariable::NonLinearTVEDG3DIPVariable(const mlawNonLinearTVE& v
_ipViscoElastic = new IPNonLinearTVE(viscoLaw.getConstRefToCompressionHardening(),
viscoLaw.getConstRefToTractionHardening(),
viscoLaw.getConstRefToKinematicHardening(),
- viscoLaw.getViscoElasticNumberOfElement());
+ viscoLaw.getViscoElasticNumberOfElement(),
+ viscoLaw.getConstRefToMullinsEffect());
this->getRefToTemperature()=Tinitial;
};
@@ -5286,7 +5287,8 @@ NonLinearTVPDG3DIPVariable::NonLinearTVPDG3DIPVariable(const mlawNonLinearTVP& v
_ipViscoElastoPlastic = new IPNonLinearTVP(viscoEPLaw.getConstRefToCompressionHardening(),
viscoEPLaw.getConstRefToTractionHardening(),
viscoEPLaw.getConstRefToKinematicHardening(),
- viscoEPLaw.getViscoElasticNumberOfElement());
+ viscoEPLaw.getViscoElasticNumberOfElement(),
+ viscoEPLaw.getConstRefToMullinsEffect());
this->getRefToTemperature()=Tinitial;
};
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 83a36302c2311349833ab4f8a8f54d45b6a1b92b..e169513e4920b4a204a98a59b76581e1aa733a36 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -246,7 +246,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
_stressLaw->stress(ipvcur,ipvprev,false,checkfrac,dTangent);
-#if 0
+#if 1
_stressLaw->stress(ipvcur,ipvprev,true,checkfrac,dTangent);
ipvcur->getConstRefToDeformationGradient().print("F Analytique"); // FLE
ipvcur->getConstRefToFirstPiolaKirchhoffStress().print("P Analytique"); // FLE
@@ -8324,6 +8324,10 @@ void NonLinearTVEDG3DMaterialLaw::setViscoElasticData(const std::string filename
_viscoLaw.setViscoElasticData(filename);
};
+void NonLinearTVEDG3DMaterialLaw::setMullinsEffect(const mullinsEffect& mullins){
+ _viscoLaw.setMullinsEffect(mullins);
+};
+
void NonLinearTVEDG3DMaterialLaw::setReferenceTemperature(const double Tref){
_viscoLaw.setReferenceTemperature(Tref);
};
@@ -8598,6 +8602,9 @@ void NonLinearTVPDG3DMaterialLaw::setKinematicHardening(const kinematicHardening
void NonLinearTVPDG3DMaterialLaw::setViscosityEffect(const viscosityLaw& v, const double p){
_viscoLaw.setViscosityEffect(v,p);
};
+void NonLinearTVPDG3DMaterialLaw::setMullinsEffect(const mullinsEffect& mullins){
+ _viscoLaw.setMullinsEffect(mullins);
+};
void NonLinearTVPDG3DMaterialLaw::setYieldPowerFactor(const double n){
_viscoLaw.setPowerFactor(n);
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 0f7b100ab0319f6772afdfc776949952f5c4c202..bfbee81b40323c9cf913860568bdeede4c4a13f8 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -2917,6 +2917,8 @@ class NonLinearTVEDG3DMaterialLaw : public dG3DMaterialLaw{ // public Material
void setTemperatureFunction_BranchBulkModuli(const scalarFunction& Tfunc,int i);
void setTemperatureFunction_BranchShearModuli(const scalarFunction& Tfunc,int i);
void setTemperatureFunction_BranchThermalDilationCoefficient(const scalarFunction& Tfunc,int i);
+
+ void setMullinsEffect(const mullinsEffect& mullins);
// For now use these Scalars - Fix by directly adding to constructor
/* void setThermalExpansionCoefficient(const double Alpha); // Fix this - for both scalar and tensor
@@ -3012,7 +3014,8 @@ class NonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{ // public Material
void setTractionHardening(const J2IsotropicHardening& trac);
void setKinematicHardening(const kinematicHardening& kin);
void setViscosityEffect(const viscosityLaw& v, const double p);
-
+ void setMullinsEffect(const mullinsEffect& mullins);
+
void setYieldPowerFactor(const double n);
void nonAssociatedFlowRuleFactor(const double b);
void setPlasticPoissonRatio(const double nup);