diff --git a/NonLinearSolver/Domain/partDomain.h b/NonLinearSolver/Domain/partDomain.h
index 49b6877986e51ef226e9c084f37f8997008e6e3e..75762d54e671e66c432c9f20d2f2c536348987e6 100644
--- a/NonLinearSolver/Domain/partDomain.h
+++ b/NonLinearSolver/Domain/partDomain.h
@@ -229,12 +229,12 @@ class partDomain
virtual void createIPState(AllIPState::ipstateContainer& map, const bool *state) const;
- virtual void computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff)=0;
+ virtual bool computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff)=0;
virtual void setDeformationGradientGradient(AllIPState *aips, const STensor33 &GM, const IPStateBase::whichState ws)=0;
virtual void setdFmdFM(AllIPState *aips, std::map<Dof, STensor3> &dUdF,const IPStateBase::whichState ws)=0;
virtual void setdFmdGM(AllIPState *aips, std::map<Dof, STensor33> &dUdG,const IPStateBase::whichState ws)=0;
virtual void setValuesForBodyForce(AllIPState *aips, const IPStateBase::whichState ws)=0;
- virtual void computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
+ virtual bool computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
materialLaw *mlaw,fullVector<double> &disp, bool stiff)=0;
virtual void setGaussIntegrationRule()=0;
virtual bool getFormulation() const{return _fullDg;}
@@ -382,13 +382,13 @@ class dgPartDomain : public partDomain
ScalarTermBase<double>* getScalarTermPFBound() const {return scalarTermPFBound;};
LinearTermBase<double>* getLinearTermPFBound() const {return linearTermPFBound;};
- virtual void computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff)=0;
- virtual void computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
+ virtual bool computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff)=0;
+ virtual bool computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
partDomain* efMinus, partDomain *efPlus,materialLaw *mlawminus,
materialLaw *mlawplus,fullVector<double> &dispm,
fullVector<double> &dispp,
const bool virt, bool stiff,const bool checkfrac=true)=0;
- virtual void computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
+ virtual bool computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
materialLaw *mlaw,fullVector<double> &disp, bool stiff)=0;
virtual void initiallyBreakAllInterfaceIP(AllIPState *aips) const;
diff --git a/NonLinearSolver/internalPoints/ipField.cpp b/NonLinearSolver/internalPoints/ipField.cpp
index 34b6c99cbf7a35bf9aea3dd550ce749fa0fea620..cad06a1bdb708f3a9e133d133d7ee0a5286b3dd7 100644
--- a/NonLinearSolver/internalPoints/ipField.cpp
+++ b/NonLinearSolver/internalPoints/ipField.cpp
@@ -1093,13 +1093,15 @@ std::string IPField::ToString(const int i){
}
};
-void IPField::compute1state(IPStateBase::whichState ws, bool stiff){
+bool IPField::compute1state(IPStateBase::whichState ws, bool stiff){
std::vector<partDomain*>* domainVector = _solver->getDomainVector();
unknownField* ufield = _solver->getUnknownField();
for(std::vector<partDomain*>::iterator itdom=domainVector->begin(); itdom!=domainVector->end(); ++itdom){
partDomain *dom = *itdom;
- dom->computeIPVariable(_AIPS,ufield,ws,stiff);
+ bool ok = dom->computeIPVariable(_AIPS,ufield,ws,stiff);
+ if (!ok) return false;
}
+ return true;
}
diff --git a/NonLinearSolver/internalPoints/ipField.h b/NonLinearSolver/internalPoints/ipField.h
index 9ef8b3b46042298a7c260003c73cd44a98419a4c..9804de242a5b0f85fc8ddbf05bd849c3b8512d69 100644
--- a/NonLinearSolver/internalPoints/ipField.h
+++ b/NonLinearSolver/internalPoints/ipField.h
@@ -582,7 +582,7 @@ class IPField : public elementsField {
const std::vector<partDomain*>& getDomainVector() const;
- void compute1state(IPStateBase::whichState ws, bool stiff);
+ bool compute1state(IPStateBase::whichState ws, bool stiff);
void setDeformationGradientGradient(const STensor33 &GM,const IPStateBase::whichState ws);
void setdFmdFM(std::map<Dof, STensor3> &dUdF, const IPStateBase::whichState ws);
void setdFmdGM(std::map<Dof, STensor33> &dUdG, const IPStateBase::whichState ws);
diff --git a/NonLinearSolver/internalPoints/ipHyperelastic.cpp b/NonLinearSolver/internalPoints/ipHyperelastic.cpp
index 75bbe43292600553a3e779d55ba4ea0820c99b4d..80e933732a4f4bbb42238953287a40a25a6707fd 100644
--- a/NonLinearSolver/internalPoints/ipHyperelastic.cpp
+++ b/NonLinearSolver/internalPoints/ipHyperelastic.cpp
@@ -14,7 +14,8 @@
#include "STensorOperations.h"
IPHyperViscoElastic::IPHyperViscoElastic(const int N):IPVariableMechanics(),_N(N),_elasticEnergy(0.),_Ee(0.),_kirchhoff(0.),
- _irreversibleEnergy(0.),_DirreversibleEnergyDF(0.), _viscousEnergyPart(0.), _dElasticEnergyPartdF(0.), _dViscousEnergyPartdF(0.), _elasticBulkPropertyScaleFactor(1.), _elasticShearPropertyScaleFactor(1.), _intA(0.), _intB(0.)
+ _irreversibleEnergy(0.),_DirreversibleEnergyDF(0.), _viscousEnergyPart(0.), _dElasticEnergyPartdF(0.), _dViscousEnergyPartdF(0.), _elasticBulkPropertyScaleFactor(1.), _elasticShearPropertyScaleFactor(1.),
+ _intA(0.), _intB(0.), _DintA(0.), _DintB(0.)
{
_A.clear();
_B.clear();
@@ -30,7 +31,8 @@ IPHyperViscoElastic::IPHyperViscoElastic(const IPHyperViscoElastic& src): IPVari
_N(src._N),_A(src._A),_B(src._B),_elasticEnergy(src._elasticEnergy), _kirchhoff(src._kirchhoff),
_irreversibleEnergy(src._irreversibleEnergy),_DirreversibleEnergyDF(src._DirreversibleEnergyDF),
_viscousEnergyPart(src._viscousEnergyPart), _dElasticEnergyPartdF(src._dElasticEnergyPartdF),
- _dViscousEnergyPartdF(src._dViscousEnergyPartdF), _elasticBulkPropertyScaleFactor(src._elasticBulkPropertyScaleFactor), _elasticShearPropertyScaleFactor(src._elasticShearPropertyScaleFactor),_intA(src._intA),_intB(src._intB), _psi_branch(src._psi_branch)
+ _dViscousEnergyPartdF(src._dViscousEnergyPartdF), _elasticBulkPropertyScaleFactor(src._elasticBulkPropertyScaleFactor), _elasticShearPropertyScaleFactor(src._elasticShearPropertyScaleFactor),
+ _intA(src._intA),_intB(src._intB), _psi_branch(src._psi_branch), _DintA(src._DintA),_DintB(src._DintB)
{
@@ -55,6 +57,8 @@ IPHyperViscoElastic& IPHyperViscoElastic::operator =(const IPVariable& src){
_elasticShearPropertyScaleFactor= psrc->_elasticShearPropertyScaleFactor;
_intA= psrc->_intA;
_intB= psrc->_intB;
+ _DintA= psrc->_DintA;
+ _DintB= psrc->_DintB;
_psi_branch= psrc->_psi_branch;
}
return *this;
@@ -77,6 +81,8 @@ void IPHyperViscoElastic::restart() {
restartManager::restart(_elasticShearPropertyScaleFactor);
restartManager::restart(_intA);
restartManager::restart(_intB);
+ restartManager::restart(_DintA);
+ restartManager::restart(_DintB);
restartManager::restart(_psi_branch);
};
diff --git a/NonLinearSolver/internalPoints/ipHyperelastic.h b/NonLinearSolver/internalPoints/ipHyperelastic.h
index 40a02848b706e865ec745005191b7f29be4b685b..e1ccba37b8d58702163381fc7aaf4aa286c91424 100644
--- a/NonLinearSolver/internalPoints/ipHyperelastic.h
+++ b/NonLinearSolver/internalPoints/ipHyperelastic.h
@@ -33,6 +33,8 @@ class IPHyperViscoElastic : public IPVariableMechanics{
double _elasticShearPropertyScaleFactor;
double _intA; // intA = psi_ExtraBranch/Kinf
double _intB; // intB = psi_ExtraBranch/Ginf
+ double _DintA; // derivative of _intA
+ double _DintB; // derivative of _intB
std::vector<double> _psi_branch; // N elements
diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVM.cpp b/NonLinearSolver/internalPoints/ipNonLinearTVM.cpp
index cd70e403b1808d18a74c9a01e4e2b51df52b7586..be3d6fa710da72276bb524fffa4f80dee755e31a 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVM.cpp
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVM.cpp
@@ -31,6 +31,7 @@ IPNonLinearTVM::IPNonLinearTVM(const J2IsotropicHardening* comp,
_Av_TVE_vector.clear(); _dAv_dTrEe_TVE_vector.clear();
_Bd_TVE_vector.clear(); _dBd_dTrEe_TVE_vector.clear(); _dBd_dDevEe_TVE_vector.clear();
_intAv_TVE_vector.clear(); _intBd_TVE_vector.clear();
+ _DintAv_TVE_vector.clear(); _DintBd_TVE_vector.clear();
for (int i=0; i < N; i++){
STensor3 el(0.);
@@ -50,6 +51,8 @@ IPNonLinearTVM::IPNonLinearTVM(const J2IsotropicHardening* comp,
_dBd_dTrEe_TVE_vector.push_back(0.);
_intAv_TVE_vector.push_back(0.);
_intBd_TVE_vector.push_back(0.);
+ _DintAv_TVE_vector.push_back(0.);
+ _DintBd_TVE_vector.push_back(0.);
_DA_DdevE.push_back(el43);
_DDA_DTDdevE.push_back(el43);
}
@@ -75,7 +78,8 @@ IPNonLinearTVM::IPNonLinearTVM(const IPNonLinearTVM& src): IPHyperViscoElastoPla
_intAv_1_TVE(src._intAv_1_TVE), _intBd_1_TVE(src._intBd_1_TVE), _intAv_2_TVE(src._intAv_2_TVE), _intBd_2_TVE(src._intBd_2_TVE),
_dBd_dDevEe(src._dBd_dDevEe), _dBd_dDevEe_TVE(src._dBd_dDevEe_TVE),
_Av_TVE_vector(src._Av_TVE_vector), _Bd_TVE_vector(src._Bd_TVE_vector), _dAv_dTrEe_TVE_vector(src._dAv_dTrEe_TVE_vector), _dBd_dTrEe_TVE_vector(src._dBd_dTrEe_TVE_vector),
- _intAv_TVE_vector(src._intAv_TVE_vector), _intBd_TVE_vector(src._intBd_TVE_vector), _dBd_dDevEe_TVE_vector(src._dBd_dDevEe_TVE_vector),
+ _intAv_TVE_vector(src._intAv_TVE_vector), _intBd_TVE_vector(src._intBd_TVE_vector),_DintAv_TVE_vector(src._DintAv_TVE_vector), _DintBd_TVE_vector(src._DintBd_TVE_vector),
+ _dBd_dDevEe_TVE_vector(src._dBd_dDevEe_TVE_vector),
_viscousDissipatedEnergy(src._viscousDissipatedEnergy), _mullinsDissipatedEnergy(src._mullinsDissipatedEnergy), _psiInfCorrector(src._psiInfCorrector){
if (src._ipMullinsEffect != NULL)
@@ -143,6 +147,8 @@ IPNonLinearTVM& IPNonLinearTVM::operator=(const IPVariable& src)
_dBd_dTrEe_TVE_vector = psrc->_dBd_dTrEe_TVE_vector;
_intAv_TVE_vector = psrc->_intAv_TVE_vector;
_intBd_TVE_vector = psrc->_intBd_TVE_vector;
+ _DintAv_TVE_vector = psrc->_DintAv_TVE_vector;
+ _DintBd_TVE_vector = psrc->_DintBd_TVE_vector;
_dBd_dDevEe_TVE_vector = psrc->_dBd_dDevEe_TVE_vector;
_viscousDissipatedEnergy = psrc->_viscousDissipatedEnergy;
_mullinsDissipatedEnergy = psrc->_mullinsDissipatedEnergy;
@@ -227,6 +233,8 @@ void IPNonLinearTVM::restart(){
restartManager::restart(_dBd_dTrEe_TVE_vector);
restartManager::restart(_intAv_TVE_vector);
restartManager::restart(_intBd_TVE_vector);
+ restartManager::restart(_DintAv_TVE_vector);
+ restartManager::restart(_DintBd_TVE_vector);
restartManager::restart(_dBd_dDevEe_TVE_vector);
restartManager::restart(_viscousDissipatedEnergy);
restartManager::restart(_mullinsDissipatedEnergy);
@@ -245,4 +253,4 @@ IPMullinsEffect& IPNonLinearTVM::getRefToIPMullinsEffect(){
if(_ipMullinsEffect==NULL)
Msg::Error("IPNonLinearTVM: _ipMullinsEffect not initialized");
return *_ipMullinsEffect;
-};
\ No newline at end of file
+};
diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVM.h b/NonLinearSolver/internalPoints/ipNonLinearTVM.h
index 9b7ada89f44c0d0f634509549357e0531b1899e5..755f0cdd36d5d89c68c2e58fa3701882acad9aa9 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVM.h
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVM.h
@@ -42,7 +42,7 @@ class IPNonLinearTVM : public IPHyperViscoElastoPlastic{
double _intAv_1_TVE, _intAv_2_TVE;
STensor3 _dBd_dDevEe;
STensor3 _dBd_dDevEe_TVE, _intBd_1_TVE, _intBd_2_TVE;
- std::vector<double> _Av_TVE_vector, _Bd_TVE_vector, _dAv_dTrEe_TVE_vector, _dBd_dTrEe_TVE_vector, _intAv_TVE_vector, _intBd_TVE_vector;
+ std::vector<double> _Av_TVE_vector, _Bd_TVE_vector, _dAv_dTrEe_TVE_vector, _dBd_dTrEe_TVE_vector, _intAv_TVE_vector, _intBd_TVE_vector, _DintAv_TVE_vector, _DintBd_TVE_vector;
std::vector<STensor3> _dBd_dDevEe_TVE_vector;
// DcorKirDT
@@ -119,4 +119,4 @@ class IPNonLinearTVM : public IPHyperViscoElastoPlastic{
virtual IPVariable* clone() const{return new IPNonLinearTVM(*this);};
virtual void restart();
};
-#endif // IPNonLinearTVM.h
\ No newline at end of file
+#endif // IPNonLinearTVM.h
diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVP.cpp b/NonLinearSolver/internalPoints/ipNonLinearTVP.cpp
index 891478f1a8298feb168f23613399c65b88c1c1e6..aec94de8f42fb64858adc4c1ded5468450f67901 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVP.cpp
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVP.cpp
@@ -18,7 +18,7 @@ IPNonLinearTVP::IPNonLinearTVP(const J2IsotropicHardening* comp, const J2Isotrop
_ModMandel(0.),_DModMandelDT(0.), _DbackSigDT(0.), _DModMandelDF(0.), _DbackSigDF(0.), _DphiPDF(0.), _GammaN(0.), _dGammaNdT(0.), _dGammaNdF(0.),
_mechSrcTVP(0.),_DmechSrcTVPdT(0.),_DmechSrcTVPdF(0.),
_dPlasticEnergyPartdT(0.),_dElasticEnergyPartdT(0.),_dViscousEnergyPartdT(0.),
- _IsoHardForce_simple(0.), _dIsoHardForcedT_simple(0.), _dIsoHardForcedF_simple(0.)
+ _IsoHardForce_simple(0.), _dIsoHardForcedT_simple(0.), _dIsoHardForcedF_simple(0.), _DphiDF(0.),_DphiDT(0.)
{
_IsoHardForce.clear(); _dIsoHardForcedT.clear(); _dIsoHardForcedF.clear();
@@ -44,7 +44,7 @@ IPNonLinearTVP::IPNonLinearTVP(const IPNonLinearTVP& src):IPNonLinearTVM(src),
_dPlasticEnergyPartdT(src._dPlasticEnergyPartdT),_dElasticEnergyPartdT(src._dElasticEnergyPartdT),
_dViscousEnergyPartdT(src._dViscousEnergyPartdT),
_IsoHardForce_simple(src._IsoHardForce_simple),
- _dIsoHardForcedT_simple(src._dIsoHardForcedT_simple), _dIsoHardForcedF_simple(src._dIsoHardForcedF_simple)
+ _dIsoHardForcedT_simple(src._dIsoHardForcedT_simple), _dIsoHardForcedF_simple(src._dIsoHardForcedF_simple), _DphiDF(src._DphiDF),_DphiDT(src._DphiDT)
{
};
@@ -84,7 +84,8 @@ IPNonLinearTVP& IPNonLinearTVP::operator =(const IPVariable &source){
_IsoHardForce_simple = ps->_IsoHardForce_simple;
_dIsoHardForcedT_simple = ps->_dIsoHardForcedT_simple;
_dIsoHardForcedF_simple = ps->_dIsoHardForcedF_simple;
-
+ _DphiDF = ps->_DphiDF;
+ _DphiDT = ps->_DphiDT;
}
return *this;
};
@@ -125,5 +126,7 @@ void IPNonLinearTVP::restart() {
restartManager::restart(_IsoHardForce_simple);
restartManager::restart(_dIsoHardForcedT_simple);
restartManager::restart(_dIsoHardForcedF_simple);
+ restartManager::restart(_DphiDF);
+ restartManager::restart(_DphiDT);
return;
}
diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVP.h b/NonLinearSolver/internalPoints/ipNonLinearTVP.h
index eeef7ca85d602a773c54d36f8b0f36fdc5361bd8..8cc1e9aed6195f5003e2f9345e44321baec0b579 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVP.h
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVP.h
@@ -42,12 +42,14 @@ class IPNonLinearTVP : public IPNonLinearTVM{
STensor3 _ModMandel;
STensor3 _DModMandelDT;
STensor3 _DbackSigDT;
+ STensor3 _DphiDT;
STensor3 _DmechSrcTVPdF;
STensor3 _dIsoHardForcedF_simple;
std::vector<STensor3> _dIsoHardForcedF;
STensor43 _DbackSigDF;
STensor43 _DModMandelDF;
+ STensor43 _DphiDF;
STensor3 _DphiPDF; // need this for mechSrc
STensor3 _GammaN;
@@ -98,6 +100,9 @@ class IPNonLinearTVP : public IPNonLinearTVM{
virtual STensor3& getRefToDphiPDF() {return _DphiPDF;}
virtual const STensor3& getConstRefToDphiPDF() const {return _DphiPDF;}
+ virtual STensor3& getRefToGammaN() {return _GammaN;}
+ virtual const STensor3& getConstRefToGammaN() const {return _GammaN;}
+
virtual double &getRefToDElasticEnergyPartdT() {return _dElasticEnergyPartdT;}
virtual const double &getConstRefToDElasticEnergyPartdT() const {return _dElasticEnergyPartdT;}
virtual double &getRefToDViscousEnergyPartdT() {return _dViscousEnergyPartdT;}
@@ -114,4 +119,4 @@ class IPNonLinearTVP : public IPNonLinearTVM{
};
-#endif // IPNonLinearTVP.h
\ No newline at end of file
+#endif // IPNonLinearTVP.h
diff --git a/NonLinearSolver/internalPoints/ipNonLocalPorosity.cpp b/NonLinearSolver/internalPoints/ipNonLocalPorosity.cpp
index bcb9c6f034d462da122dfb2253667a926597a4c0..612b94ca8306e3c6ecda88945f492f0645846a98 100644
--- a/NonLinearSolver/internalPoints/ipNonLocalPorosity.cpp
+++ b/NonLinearSolver/internalPoints/ipNonLocalPorosity.cpp
@@ -21,7 +21,8 @@ IPNonLocalPorosity::IPNonLocalPorosity() : IPVariableMechanics(), _elasticEnerg
_currentOutwardNormal(NULL),_referenceOutwardNormal(NULL),_corKir(0.), _yieldStress(0.),
_failed(false), _hatq(0.),_nonlocalEplmatrix(0.),_nonlocalHatQ(0.), _fVinitial(0.),
_eplMacro(0.), _eplMacroNonLocal(0.), _lnfV(0.),_lnfVtilde(0.), _T(288.), _dissipationActive(false),
- _thermalEnergy(0.),_ThermoElasticCoupling(0.), _DplasticEnergyDEplmatrix(0.)
+ _thermalEnergy(0.),_ThermoElasticCoupling(0.), _DplasticEnergyDEplmatrix(0.),
+ _successFlag(true)
{
// Default constructor without arguments is forbiden as it is incorrect and incomplete
ipvJ2IsotropicHardening=NULL;
@@ -48,7 +49,8 @@ IPNonLocalPorosity::IPNonLocalPorosity(double fVinitial, const J2IsotropicHarden
_currentOutwardNormal(NULL), _referenceOutwardNormal(NULL),_deformationGradient(NULL),
_corKir(0.),_failed(false),_hatq(0.),_nonlocalEplmatrix(0.),_nonlocalHatQ(0.), _fVinitial(fVinitial),
_eplMacro(0.),_eplMacroNonLocal(0.), _lnfV(0.),_lnfVtilde(0.),_T(288.), _dissipationActive(false),
- _thermalEnergy(0.),_ThermoElasticCoupling(0.), _DplasticEnergyDEplmatrix(0.)
+ _thermalEnergy(0.),_ThermoElasticCoupling(0.), _DplasticEnergyDEplmatrix(0.),
+ _successFlag(true)
{
// Real constructor
// Initialisation of pointers for internal laws Ipvs and other related variables
@@ -158,7 +160,8 @@ IPNonLocalPorosity::IPNonLocalPorosity(const IPNonLocalPorosity &source) :
_eplMacro(source._eplMacro), _eplMacroNonLocal(source._eplMacroNonLocal),
_lnfV(source._lnfV),_lnfVtilde(source._lnfVtilde),_T(source._T), _dissipationActive(source._dissipationActive),
_thermalEnergy(source._thermalEnergy),_ThermoElasticCoupling(source._ThermoElasticCoupling),
- _DplasticEnergyDEplmatrix(source._DplasticEnergyDEplmatrix)
+ _DplasticEnergyDEplmatrix(source._DplasticEnergyDEplmatrix),
+ _successFlag(source._successFlag)
{
// Copy all pointers for internal laws Ipvs
@@ -342,6 +345,7 @@ IPNonLocalPorosity& IPNonLocalPorosity::operator=(const IPVariable &source)
_thermalEnergy = src->_thermalEnergy;
_ThermoElasticCoupling = src->_ThermoElasticCoupling;
_DplasticEnergyDEplmatrix = src->_DplasticEnergyDEplmatrix;
+ _successFlag = src->_successFlag;
}
return *this;
};
@@ -417,6 +421,7 @@ void IPNonLocalPorosity::restart()
restartManager::restart(_thermalEnergy);
restartManager::restart(_ThermoElasticCoupling);
restartManager::restart(_DplasticEnergyDEplmatrix);
+ restartManager::restart(_successFlag);
// Pointer to variables store somewhere else
_currentOutwardNormal = NULL;
_referenceOutwardNormal = NULL;
diff --git a/NonLinearSolver/internalPoints/ipNonLocalPorosity.h b/NonLinearSolver/internalPoints/ipNonLocalPorosity.h
index 2406081b20ac3b098fc9faa9fc3e316743fdd4d8..acbde87d7fa796337c2b96f18a02045bfb3f34c8 100644
--- a/NonLinearSolver/internalPoints/ipNonLocalPorosity.h
+++ b/NonLinearSolver/internalPoints/ipNonLocalPorosity.h
@@ -95,6 +95,8 @@ protected:
double _T; // temperature
double _thermalEnergy; // thermal energy
double _ThermoElasticCoupling; //
+
+ bool _successFlag; // true if contitutive law is successfully passed in GP
public:
// Constructor & destructor
@@ -341,7 +343,8 @@ protected:
};
-
+ virtual bool constitutiveSuccessFlag() const {return _successFlag;};
+ virtual bool& getRefToConstitutiveSuccessFlag() {return _successFlag;};
// - Coalesncence and void state ipv
virtual const IPCoalescence& getConstRefToIPCoalescence() const{
diff --git a/NonLinearSolver/internalPoints/ipvariable.h b/NonLinearSolver/internalPoints/ipvariable.h
index 744401d56ef7093f8992816865ea46db9f2044d2..adb7cf7ec4bfa149b83f21dc122197ebaf09b768 100644
--- a/NonLinearSolver/internalPoints/ipvariable.h
+++ b/NonLinearSolver/internalPoints/ipvariable.h
@@ -67,6 +67,8 @@ class IPVariable
virtual void setValueForBodyForce(const STensor43& K, const int ModuliType)=0;
virtual void computeBodyForce(const STensor33& G, SVector3& B) const=0;
virtual const STensor43 &getConstRefToTangentModuli() const=0;
+
+ virtual bool getConstitutiveSuccessFlag() const {return true;}; // true if contitutive law is successfully passed in GP
#if defined(HAVE_MPI)
// using in multiscale analysis with MPI
@@ -120,6 +122,8 @@ class IPVariableMechanics : public IPVariable{
virtual void computeBodyForce(const STensor33& G, SVector3& B) const { Msg::Error("ComputeBodyForce is not defined in IPVariableMechanics"); }
virtual const STensor43 &getConstRefToTangentModuli() const {Msg::Error("getConstRefToTangentModuli IPVariableMechanics"); static STensor43 dummy; return dummy;};
+
+
virtual void restart(){
IPVariable::restart();
}
@@ -231,6 +235,26 @@ template<class Tbulk,class Tfrac> class IPVariable2Enhanced : public Tbulk{
}
return;
}
+
+ virtual bool getConstitutiveSuccessFlag() const
+ {
+ if (_ipvbulk!=NULL && _ipvfrac!=NULL)
+ {
+ return _ipvbulk->getConstitutiveSuccessFlag() and _ipvfrac->getConstitutiveSuccessFlag();
+ }
+ else if (_ipvbulk!=NULL)
+ {
+ return _ipvbulk->getConstitutiveSuccessFlag();
+ }
+ else if (_ipvfrac!= NULL)
+ {
+ return _ipvfrac->getConstitutiveSuccessFlag();
+ }
+ else
+ {
+ return true;
+ }
+ }
virtual IPVariable* clone() const =0;
#if defined(HAVE_MPI)
diff --git a/NonLinearSolver/materialLaw/kinematicHardening.cpp b/NonLinearSolver/materialLaw/kinematicHardening.cpp
index ec2a7b7c2121eeba524895ff8dc89bf377a5f707..5cc01f665d5ad65c0f9b52c919e19afbe33fdf7e 100644
--- a/NonLinearSolver/materialLaw/kinematicHardening.cpp
+++ b/NonLinearSolver/materialLaw/kinematicHardening.cpp
@@ -189,12 +189,13 @@ void PolynomialKinematicHardening::hardening(double p0, const IPKinematicHardeni
Msg::Error("order %d is not implemented PolynomialJ2IsotropicHaderning::hardening",_order);
}
- R *= _temFunc_K->getVal(T);
- dR *= _temFunc_K->getVal(T);
- ddR *= _temFunc_K->getVal(T);
- dRdT = R*_temFunc_K->getDiff(T);
- ddRdTT = R*_temFunc_K->getDoubleDiff(T);
-
+ R *= _temFunc_K->getVal(T); // Hb
+ dR *= _temFunc_K->getVal(T); // dHb
+ ddR *= _temFunc_K->getVal(T); // ddHb
+ dRdT = R*_temFunc_K->getDiff(T)/_temFunc_K->getVal(T); // dHbdT
+ ddRdTT = R*_temFunc_K->getDoubleDiff(T)/_temFunc_K->getVal(T); // ddHbddT
+ ddRdT = dR*_temFunc_K->getDiff(T)/_temFunc_K->getVal(T); // ddHbdT = ddHbdTdgamma
+
ipv.set(R,dR,ddR,dRdT,ddRdTT,ddRdT);
}
@@ -405,7 +406,6 @@ void exponentialKinematicHardeningLaw::hardening(double p0, const IPKinematicHar
ipv.set(R,dR,ddR,dRdT,ddRdTT,ddRdT);
}
-// added
void exponentialKinematicHardeningLaw::hardening(double p0, const IPKinematicHardening &ipvprev, double p, double T, IPKinematicHardening &ipv) const
{
double R(0.), dR(0.), ddR(0.), dRdT(0.), ddRdTT(0.), ddRdT(0.);
@@ -430,7 +430,7 @@ void exponentialKinematicHardeningLaw::hardening(double p0, const IPKinematicHar
}
ipv.set(R,dR,ddR,dRdT,ddRdTT,ddRdT); //added dRdT
-} //
+}
kinematicHardening * exponentialKinematicHardeningLaw::clone() const
{
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
index 048f868c610d773c3a88ce77c850fd155069a4a9..e15627afd345e1775f0abda078bd5edca516a235 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
@@ -193,7 +193,7 @@ void mlawHyperViscoElastic::setViscoElasticData(const std::string filename){
};
void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &intA, double &B_d, STensor3 &dB_vddev, double &intB,
- double* dB_dTrEe, double* psiInfCorrector) const{
+ double* dB_dTrEe, double* psiInfCorrector, double* DintA, double* DintB) const{
if(_extraBranchNLType == expType)
{
A_v = getVolumeCorrection()*pow(exp(getXiVolumeCorrection()/3.*tr*tr)-1.,getZetaVolumeCorrection());
@@ -327,47 +327,78 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
}
else if(_extraBranchNLType == TensionCompressionRegularisedType){
+
+ // To make it simpler for the sake of it
+ double V0 = getVolumeCorrection(); // This code doesnt use it
+ double V1 = getXiVolumeCorrection();
+ double V2 = getZetaVolumeCorrection();
+ double V3 = _volCorrection2;
+ double V4 = _xivolCorrection2;
+ double V5 = _zetavolCorrection2;
+ double D0 = getDevCorrection(); // This code doesnt use it
+ double D1 = getThetaDevCorrection();
+ double D2 = getPiDevCorrection();
+ double D3 = _devCorrection2;
+ double D4 = _thetadevCorrection2;
+ double D5 = _pidevCorrection2;
+ double Ci = _compCorrection;
// Regularising function
double m = _tensionCompressionRegularisation;
double expmtr = exp(-m*tr);
// if (exp(-m*tr)<1.e+10){ expmtr = exp(-m*tr);}
double sigmoid = 1/(1.+expmtr);
-
- double x = getXiVolumeCorrection()*tr*tr + getZetaVolumeCorrection();
- A_v = -1. + getVolumeCorrection()/sqrt(x)*(sigmoid + _compCorrection*(1.-sigmoid));
+
+ // A_v
+ double x1 = V1*tr*tr + V2;
+ double x2 = V4*tr*tr + V5;
+ A_v = -1. + sigmoid*(1./sqrt(x1)+V0) + (1.-sigmoid)*Ci*(1/sqrt(x2)+V3);
+
+ dA_vdE = - sigmoid*V1*tr/pow(x1,1.5) - (1.-sigmoid)*(V4*Ci*tr/pow(x2,1.5))
+ + (m*expmtr/pow((1.+expmtr),2.))*(1./sqrt(x1)+V0) - (Ci*m*expmtr/pow((1.+expmtr),2.))*(1/sqrt(x2)+V3);
- dA_vdE = - getXiVolumeCorrection()*tr/pow(x,1.5) *(sigmoid + _compCorrection*(1.-sigmoid))
- + 1./sqrt(x)* (m*expmtr/pow((1.+expmtr),2.)) *(1.-_compCorrection);
- dA_vdE *= getVolumeCorrection();
-
- if(getXiVolumeCorrection()>0.){
- double integrand = 1./getXiVolumeCorrection()*sqrt(x); // integral of A_v * trEe
- integrand -= ( 1./getXiVolumeCorrection()*sqrt(getZetaVolumeCorrection()) ); // value at trEe = 0.
- intA = - tr*tr/2. + getVolumeCorrection()*integrand*(sigmoid + _compCorrection*(1.-sigmoid));
+ if(V1>0. && V4>0.){
+ double integrand1 = 1./V1*sqrt(x1) + V0*tr*tr/2.; // integral of A_v * trEe
+ double integrand2 = 1./V4*sqrt(x2) + V3*tr*tr/2.;
+ integrand1 -= ( 1./V1*sqrt(V2) ); // value at trEe = 0.
+ integrand2 -= ( 1./V4*sqrt(V5) ); // value at trEe = 0.
+ intA = - tr*tr/2. + sigmoid*integrand1 + Ci*(1.-sigmoid)*integrand2;
+
+ if(DintA!=NULL){
+ *DintA = (m*expmtr/pow((1.+expmtr),2.))*integrand1 - (Ci*m*expmtr/pow((1.+expmtr),2.))*integrand2;
+ }
}
else{
intA = - tr*tr/2.;
+ // Msg::Error(" Inside extraBranch_Hyper Type = TensionCompressionRegularisedType, V1 = 0. or V4 = 0., incompatibility with Mullin's effect");
}
-
- double y = getThetaDevCorrection()*dev.dotprod() + getPiDevCorrection();
- B_d = -1. + getDevCorrection()/sqrt(y)*(sigmoid + _compCorrection*(1.-sigmoid));
+
+ // B_d
+ double y1 = D1*dev.dotprod() + D2;
+ double y2 = D4*dev.dotprod() + D5;
+ B_d = -1. + sigmoid*(1./sqrt(y1)+D0) + (1.-sigmoid)*Ci*(1/sqrt(y2)+D3);
STensorOperation::zero(dB_vddev);
dB_vddev = dev;
- dB_vddev *= -getThetaDevCorrection()/pow(y,1.5)*(sigmoid + _compCorrection*(1.-sigmoid));
- dB_vddev *= getDevCorrection();
-
+ dB_vddev *= (-sigmoid*D1/pow(y1,1.5) - (1.-sigmoid)*(D4*Ci/pow(y2,1.5)));
+
if(dB_dTrEe!=NULL){
- *dB_dTrEe = getDevCorrection()/sqrt(y)* (m*expmtr/pow((1.+expmtr),2.)) *( 1.-_compCorrection);
+ *dB_dTrEe = (m*expmtr/pow((1.+expmtr),2.))*(1./sqrt(y1)+D0) - (Ci*m*expmtr/pow((1.+expmtr),2.))*(1/sqrt(y2)+D3);
}
- if(getThetaDevCorrection()>0.){
- double integrand = 1./getThetaDevCorrection()*sqrt(y); // integral of B_d * devEe
- integrand -= (1./getThetaDevCorrection()*sqrt(getPiDevCorrection()) ); // value at devEe = 0.
- intB = -dev.dotprod()/2. + getDevCorrection()*integrand*(sigmoid + _compCorrection*(1.-sigmoid));
+ if(D1>0. && D4>0.){
+ double integrand1 = 1./D1*sqrt(y1) + D0*dev.dotprod()/2.; // integral of B_d * devEe
+ double integrand2 = 1./D4*sqrt(y2) + D3*dev.dotprod()/2.;
+ integrand1 -= ( 1./D1*sqrt(D2) ); // value at devEe = 0.
+ integrand2 -= ( 1./D4*sqrt(D5) ); // value at devEe = 0.
+ intB = -dev.dotprod()/2. + sigmoid*integrand1 + Ci*(1.-sigmoid)*integrand2;
+
+ if(DintB!=NULL){
+ *DintB = (m*expmtr/pow((1.+expmtr),2.))*integrand1 - (Ci*m*expmtr/pow((1.+expmtr),2.))*integrand2;
+ }
}
else{
intB = -dev.dotprod()/2.;
+ // Msg::Error(" Inside extraBranch_Hyper Type = TensionCompressionRegularisedType, D1 = 0. or D4 = 0., incompatibility with Mullin's effect");
}
}
/*else if(_extraBranchNLType == hyper_exp_TCasymm_Type){
@@ -445,6 +476,7 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
double D3 = _devCorrection2;
double D4 = _thetadevCorrection2;
double D5 = _pidevCorrection2;
+ double Ci = _compCorrection;
// Regularising function
double m = _tensionCompressionRegularisation;
@@ -468,6 +500,10 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
integrand1 -= ( 1./V1*sqrt(V2) ); // value at trEe = 0.
intA = - tr*tr/2. + sigmoid*(integrand1 + integrand2) + _compCorrection*(1.-sigmoid)*integrand1;
// intA *= getVolumeCorrection();
+
+ if(DintA!=NULL){
+ *DintA = (m*expmtr/pow((1.+expmtr),2.))*integrand1 - (Ci*m*expmtr/pow((1.+expmtr),2.))*integrand2;
+ }
}
else{
intA = - tr*tr/2.;
@@ -496,6 +532,10 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
integrand1 -= (1./D1*sqrt(D2) ); // value at devEe = 0.
intB = -dev.dotprod()/2. + sigmoid*(integrand1 + integrand2) + _compCorrection*(1.-sigmoid)*integrand1;
// intB *= getDevCorrection();
+
+ if(DintB!=NULL){
+ *DintB = (m*expmtr/pow((1.+expmtr),2.))*integrand1 - (Ci*m*expmtr/pow((1.+expmtr),2.))*integrand2;
+ }
}
else{
intB = -dev.dotprod()/2.;
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.h b/NonLinearSolver/materialLaw/mlawHyperelastic.h
index 853e17338dbd0c728adcd270f1b3cee4351159c5..78a2656834578f4ed2c167a2a79d7556b4c81bdc 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.h
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.h
@@ -88,7 +88,7 @@ class mlawHyperViscoElastic : public materialLaw{
const bool stiff, STensor43& Tangent) const;
void evaluatePhiPCorrection(double trEe, const STensor3 &devEpr, double &A_v, double &dAprdEepr, double &intA, double &B_d, STensor3 &dB_vddevEe, double &intB,
- double* dB_dTrEe = NULL, double* psiInfCorrector = NULL) const;
+ double* dB_dTrEe = NULL, double* psiInfCorrector = NULL, double* DintA = NULL, double* DintB = NULL) const;
//void evaluateA_dA(double trEe, double &A_v, double &dAprdEepr) const;
#endif // SWIG
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
index d76a3995f44178fab8202dc165efbe061610dd29..7db46ba61b9be195ea518d152d4eb08fc27de70f 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
@@ -15,11 +15,11 @@
mlawNonLinearTVENonLinearTVP::mlawNonLinearTVENonLinearTVP(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):
- mlawNonLinearTVP(num, E, nu, rho, tol, Tinitial, Alpha, KThCon, Cp, matrixbyPerturbation, pert, thermalEstimationPreviousConfig){
+ const bool matrixbyPerturbation, const double pert, const double stressIteratorTol, const bool thermalEstimationPreviousConfig):
+ mlawNonLinearTVP(num, E, nu, rho, tol, Tinitial, Alpha, KThCon, Cp, matrixbyPerturbation, pert, thermalEstimationPreviousConfig), _stressIteratorTol(stressIteratorTol) {
};
-mlawNonLinearTVENonLinearTVP::mlawNonLinearTVENonLinearTVP(const mlawNonLinearTVENonLinearTVP& src): mlawNonLinearTVP(src) {
+mlawNonLinearTVENonLinearTVP::mlawNonLinearTVENonLinearTVP(const mlawNonLinearTVENonLinearTVP& src): mlawNonLinearTVP(src), _stressIteratorTol(src._stressIteratorTol) {
};
mlawNonLinearTVENonLinearTVP& mlawNonLinearTVENonLinearTVP::operator=(const materialLaw& source){
@@ -27,6 +27,7 @@ mlawNonLinearTVENonLinearTVP& mlawNonLinearTVENonLinearTVP::operator=(const mate
mlawNonLinearTVP::operator=(source);
const mlawNonLinearTVENonLinearTVP* src =dynamic_cast<const mlawNonLinearTVENonLinearTVP*>(&source);
if(src != NULL){
+ _stressIteratorTol = src->_stressIteratorTol;
}
return *this;
};
@@ -34,21 +35,23 @@ mlawNonLinearTVENonLinearTVP& mlawNonLinearTVENonLinearTVP::operator=(const mate
mlawNonLinearTVENonLinearTVP::~mlawNonLinearTVENonLinearTVP(){
};
-double mlawNonLinearTVENonLinearTVP::freeEnergyPlasticity(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double T0, const double T) const{
+double mlawNonLinearTVENonLinearTVP::freeEnergyPlasticity(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double& T0, const double& T) const{
double Psy_TVP(0.);
// 1. Kinematic Hardening
double Hb(0.);
if (q1->_ipKinematic != NULL){
- Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
+ Hb = q1->_ipKinematic->getR(); // kinematic hardening parameter
}
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
static STensor3 alpha;
- alpha = q1->_backsig;
- alpha *= 1/(pow(kk,2)*Hb);
-
+ STensorOperation::zero(alpha);
+ if(Hb>0.){
+ alpha = q1->_backsig;
+ alpha *= 1/(pow(kk,1)*Hb); //pow(kk,2) debug
+ }
STensorOperation::doubleContractionSTensor3(q1->_backsig,alpha,Psy_TVP);
// 2. Isotropic Hardening
@@ -58,7 +61,7 @@ double mlawNonLinearTVENonLinearTVP::freeEnergyPlasticity(const IPNonLinearTVP *
Psy_TVP += intR;
- // If mullins, then calculate psi_VE derivatives
+ // If mullins, then calculate psi_VP derivatives
// if (_mullinsEffect != NULL && q._ipMullinsEffect != NULL){
// *DpsiDT = 0.;
// STensorOperation::zero(*DpsiDE);
@@ -66,6 +69,57 @@ double mlawNonLinearTVENonLinearTVP::freeEnergyPlasticity(const IPNonLinearTVP *
return Psy_TVP;
}
+void mlawNonLinearTVENonLinearTVP::freeEnergyPlasticityDerivatives(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double& T0, const double& T, STensor3& dPsy_TVPdF, double& dPsy_TVPdT) const{
+
+ STensorOperation::zero(dPsy_TVPdF);
+ dPsy_TVPdT = 0.;
+
+ // 1. Kinematic Hardening
+ double Hb(0.), dHb(0.), dHbdT(0.);
+ if (q1->_ipKinematic != NULL){
+ Hb = q1->_ipKinematic->getR(); // kinematic hardening parameter
+ dHb = q1->_ipKinematic->getDR(); // kinematic hardening parameter derivative (dHb/dgamma)
+ dHbdT = q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file
+ }
+
+ double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
+ if(Hb>0.){
+
+ static STensor3 term1, term2;
+ STensorOperation::multSTensor3STensor43(q1->_backsig,q1->_DbackSigDF,term1);
+ term2 = q1->_DgammaDF;
+ term2 *= (-1./(2*Hb)*dHb*q1->_backsig.dotprod());
+
+ dPsy_TVPdF = term1;
+ dPsy_TVPdF += term2;
+ dPsy_TVPdF *= 1./(pow(kk,1)*Hb); //pow(kk,2) debug
+
+ double term3, term4;
+ STensorOperation::doubleContractionSTensor3(q1->_backsig,q1->_DbackSigDT,term3);
+ term4 = -(dHbdT + dHb*q1->_DgammaDT)/(2*Hb)*q1->_backsig.dotprod();
+ dPsy_TVPdT = (term3+term4)/(pow(kk,1)*Hb); //pow(kk,2) debug
+
+ }
+
+
+ // 2. Isotropic Hardening
+
+ // Get Hc, etc.
+ double intR = q1->_ipCompression->getIntegR();
+ double sigc0 = q1->_ipCompression->getR0(); // Initial yield stress
+ intR -= sigc0*q1->_epspbarre;
+
+ double R(0.), dsigc0dT(0.), sigc(0.);
+ dsigc0dT = sigc0*_temFunc_Sy0->getDiff(T);
+ sigc = q1->_ipCompression->getR();
+
+ // Get R
+ R = sigc - sigc0;
+
+ dPsy_TVPdF += R*q1->_DgammaDF;
+ dPsy_TVPdT += (R*q1->_DgammaDT + _temFunc_Sy0->getDiff(T)/ _temFunc_Sy0->getVal(T)*intR - dsigc0dT*q1->_epspbarre);
+}
+
void mlawNonLinearTVENonLinearTVP::constitutive( // This function is just a placeholder, defined in the pure virtual class -> does nothing, its never called.
const STensor3& F0, // previous deformation gradient (input @ time n)
const STensor3& F1, // current deformation gradient (input @ time n+1)
@@ -372,10 +426,11 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
double Hb(0.), dHb(0.), ddHb(0.), dHbdT(0.), ddHbdgammadT(0.), ddHbddT(0.);
if (q1->_ipKinematic != NULL){
- Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
- dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter derivative (dHb/dgamma)
- dHbdT = Hb * q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file
- ddHbddT = Hb * q1->_ipKinematic->getDDRDTT();
+ Hb = q1->_ipKinematic->getR(); // kinematic hardening parameter
+ dHb = q1->_ipKinematic->getDR(); // kinematic hardening parameter derivative (dHb/dgamma)
+ dHbdT = q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file
+ ddHbddT = q1->_ipKinematic->getDDRDTT();
+ ddHbdgammadT = q1->_ipKinematic->getDDRDT();
}
double eta(0.),Deta(0.),DetaDT(0.);
@@ -648,10 +703,11 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
mlawNonLinearTVP::hardening(q0,q1,T);
mlawNonLinearTVP::getYieldCoefficients(q1,a);
if (q1->_ipKinematic != NULL){
- Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
- dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter derivative (dHb/dgamma ??)
- dHbdT = Hb * q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
- ddHbddT = Hb * q1->_ipKinematic->getDDRDTT();
+ Hb = q1->_ipKinematic->getR(); // kinematic hardening parameter
+ dHb = q1->_ipKinematic->getDR(); // kinematic hardening parameter derivative (dHb/dgamma ??)
+ dHbdT = q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
+ ddHbddT = q1->_ipKinematic->getDDRDTT();
+ ddHbdgammadT = q1->_ipKinematic->getDDRDT();
}
//a.print("a update");
@@ -732,11 +788,12 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
// Correct Normal, H = expGN
STensorOperation::decomposeDevTr(N,devN,trN);
- // Get GammaN for mechSrcTVP
+ // Get GammaN for mechSrcTVP -> done within the loop already
+ /*
STensorOperation::zero(GammaN);
GammaN = N;
GammaN *= Gamma;
- q1->_GammaN = GammaN;
+ q1->_GammaN = GammaN;*/
// Correct plastic deformation tensor
STensorOperation::multSTensor3(expGN,Fp0,Fp1);
@@ -765,7 +822,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
Ge_Tensor = _I4;
Ge_Tensor *= Ge*2; // *2 because the function doesnt do it
Ge_Tensor += Bd_stiffnessTerm;
-
+
// Correct backstress;
devX = (pow(kk,1)*Hb*Gamma*devN + devXn); // pow(kk,2) DEBUG
devX *= 1./Cxdev;
@@ -774,14 +831,14 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
q1->_backsig(0,0) += trX/3.;
q1->_backsig(1,1) += trX/3.;
q1->_backsig(2,2) += trX/3.;
-
+
// Correct Mandel
q1->_ModMandel = devPhi;
q1->_ModMandel += q1->_backsig;
q1->_ModMandel(0,0) += (ptilde);
q1->_ModMandel(1,1) += (ptilde);
q1->_ModMandel(2,2) += (ptilde); */
-
+
} // 2nd if
else{
q1->getRefToDissipationActive() = false;
@@ -808,32 +865,6 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
}
- // Elastic energy and Mullin's effect
- double& Dpsi_DT = q1->getRefTo_Dpsi_DT();
- STensor3& Dpsi_DE = q1->getRefTo_Dpsi_DE();
- q1->_elasticEnergy = mlawNonLinearTVM::freeEnergyMechanical(*q0,*q1,T0,T,&Dpsi_DT,&Dpsi_DE);
-
- // TVP energy
- double psi_TVP = freeEnergyPlasticity(q0,q1,T0,T);
- // q1->_elasticEnergy += psi_TVP;
-
- // Mullins Effect
- if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){
- mlawNonLinearTVM::calculateMullinsEffectScaler(q0, q1, T, &Dpsi_DT);
- }
- double eta_mullins = 1.;
-
- if (q1->_ipMullinsEffect != NULL){
-
- // Put Stress in IP
- q1->_P_cap = P;
-
- // Modify Stress for mullins
- eta_mullins = q1->_ipMullinsEffect->getEta();
- P *= eta_mullins;
- q1->_elasticEnergy *= eta_mullins;
- }
-
// Msg::Error(" Inside TVP, after correction, after updating q1->_elasticEnergy = %e, eta = %e !!",q1->_elasticEnergy, eta_mullins);
// q1->getRefToElasticEnergy()=deformationEnergy(*q1);
@@ -865,8 +896,13 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
// ThermSrc and MechSrc after dPdF and dPdT - But need the following tensors for the mechSrcTVP function call
static STensor3 DphiPDF;
+ STensorOperation::zero(DphiPDF);
STensorOperation::zero(dFpdF); //dFpdT
- static STensor43 DEeDCepr, DEeDF;
+ static STensor43 CeprToF, DEeDCepr, DEeDF, dGammaNdCepr; // conversion tensors
+ STensorOperation::zero(CeprToF);
+ STensorOperation::zero(DEeDCepr);
+ STensorOperation::zero(DEeDF);
+ STensorOperation::zero(dGammaNdCepr);
// I didnt make this
if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
@@ -974,7 +1010,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
//4. get DdevphiDCepr and DphiPprDCepr; and DphiEprDCepr
static STensor3 DphiPDCepr, DphiEDCepr, DphiEDdevPhi;
- static STensor43 DdevphiDCepr;
+ static STensor43 DdevphiDCepr, DphiDCepr;
DphiPDCepr = Dp_pr_DCepr;
DdevphiDCepr = DdevMeprDCepr;
@@ -986,9 +1022,16 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
STensorOperation::multSTensor3STensor43(DphiEDdevPhi,DdevphiDCepr,DphiEDCepr);
+ DphiDCepr = DdevphiDCepr;
+ for (int i=0; i<3; i++)
+ for (int s=0; s<3; s++)
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++)
+ DphiDCepr(i,s,k,l) += _I(i,s)*DphiPDCepr(k,l);
+
// 6. to 11. (inside if loop-> Gamma > 0.)
static STensor3 dAdCepr, dfDCepr, dgDCepr, DGDCepr, DgammaDCepr, DtrNDCepr, dTrXdCepr;
- static STensor43 DdevNDCepr, dFpDCepr, dDevXdCepr, dXdCepr, dGammaNdCepr;
+ static STensor43 DdevNDCepr, dFpDCepr, dDevXdCepr, dXdCepr;
if (Gamma >0.){
@@ -1117,6 +1160,13 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
STensorOperation::multSTensor43(Dho2_u_inv,DdevphiDCepr_RHS,DdevphiDCepr);
STensorOperation::multSTensor3STensor43(DphiEDdevPhi,DdevphiDCepr,DphiEDCepr);
+ DphiDCepr = DdevphiDCepr;
+ for (int i=0; i<3; i++)
+ for (int s=0; s<3; s++)
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++)
+ DphiDCepr(i,s,k,l) += _I(i,s)*DphiPDCepr(k,l);
+
//9. get DtrNDCepr DdevNdCepr
DtrNDCepr = DphiPDCepr;
DtrNDCepr *= (2.*_b);
@@ -1218,7 +1268,6 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
}
// 12. get CeprToF conversion tensor
- static STensor43 CeprToF;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
@@ -1249,13 +1298,16 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
// 15. get DgammaDF
STensor3& DgammaDF = q1->_DgammaDF;
STensor3& DGammaDF = q1->_DGammaDF;
+ STensor43& DphiDF = q1->_DphiDF;
STensorOperation::zero(DgammaDF);
STensorOperation::zero(DGammaDF);
+ STensorOperation::zero(DphiDF);
if (Gamma > 0){
STensorOperation::multSTensor3STensor43(DgammaDCepr,CeprToF,DgammaDF);
STensorOperation::multSTensor3STensor43(DGDCepr,CeprToF,DGammaDF);
STensorOperation::multSTensor43(dFpDCepr,CeprToF,dFpdF);
}
+ STensorOperation::multSTensor43(DphiDCepr,CeprToF,DphiDF);
// 16. Everything else - Conversion Tensors
static STensor43 DinvFpDF, dCedF, dCeinvDF;
@@ -1366,16 +1418,26 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
// 3. get dPhipprDT and dDevPhiprDT
double dPhiPdT(0.);
static STensor3 dDevPhiDT;
+ STensor3& DphiDT = q1->_DphiDT;
dPhiPdT = dpMeprDT;
dDevPhiDT = dDevMeprDT;
+ DphiDT = dDevPhiDT;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ DphiDT(i,j) += dPhiPdT*_I(i,j);
+
// 4. get dPhiEdT
double dPhiEdT(0.);
STensorOperation::doubleContractionSTensor3(DphiEDdevPhi,dDevPhiDT,dPhiEdT);
// 5. to 11. (inside if loop-> Gamma > 0.)
- double dAdT(0.), dfdT(0.), dgdT(0.), dGammaDT(0.), DgammaDT(0.), DtrNdT(0.), DtrXdT(0.);
+ double& dGammaDT = q1->_DGammaDT;
+ double& DgammaDT = q1->_DgammaDT;
+ dGammaDT = 0.; DgammaDT =0.;
+
+ double dAdT(0.), dfdT(0.), dgdT(0.), DtrNdT(0.), DtrXdT(0.);
static STensor3 DdevNdT, dFpDT, DdevXdT;
double eta(0.), Deta(0.), DetaDT(0.);
@@ -1489,7 +1551,12 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
}
STensorOperation::multSTensor3STensor43(DdevphiDT_RHS,Dho2_u_inv,dDevPhiDT);
-
+
+ DphiDT = dDevPhiDT;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ DphiDT(i,j) += dPhiPdT*_I(i,j);
+
// =============================== ROUND 2 - END ========================================================================= //
// 12. get DtrNdT, DdevNdT
@@ -1616,53 +1683,6 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
}
}
-
- // Important conversion tensors
- STensorOperation::zero(DEeDCepr);
- STensorOperation::zero(DEeDF);
-
- // DEeDCepr to DEeDF
- // STensorOperation::multSTensor43(_I4, DlnDCepr, DEeDCepr);
- DEeDCepr = DlnDCepr;
- DEeDCepr *= 0.5;
- DEeDCepr -= dGammaNdCepr;
- STensorOperation::multSTensor43(DEeDCepr, CeprToF, DEeDF);
-
- // Note: DEeDT = - dGammaNdT
-
- // Check DEeDF
-
- // Compute Tangents for Mullin's Effect
- if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){
-
-
- // 1st Term
- Tangent *= q1->_ipMullinsEffect->getEta();
- dPdT *= q1->_ipMullinsEffect->getEta();
-
- // 2nd Term
- static double Deta_mullins_DT(0.);
- static STensor3 Deta_mullins_DF, Deta_mullins_DEe;
-
- // Deta_mullins_DF
- STensorOperation::multSTensor3STensor43(q1->_DpsiDE,DEeDF,Deta_mullins_DF);
- Deta_mullins_DF *= q1->_DmullinsDamage_Dpsi_cap;
-
- // Deta_mullins_DT
- STensorOperation::doubleContractionSTensor3(q1->_DpsiDE,dGammaNdT,Deta_mullins_DT);
- Deta_mullins_DT *= -q1->_DmullinsDamage_Dpsi_cap;
- Deta_mullins_DT += q1->_DmullinsDamage_DT;
-
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- dPdT(i,j) += q1->_P_cap(i,j) * Deta_mullins_DT;
- for (int k=0; k<3; k++)
- for (int l=0; l<3; l++)
- Tangent(i,j,k,l) += q1->_P_cap(i,j) * Deta_mullins_DF(k,l);
- }
-
- }
-
// TVP - flux derivatives
// fluxT
dfluxTdT = fluxT;
@@ -1733,14 +1753,14 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
double& Wm_TVE = q1->getRefToMechSrcTVE(); // TVE
mechanicalSource = 0.;
- // mlawNonLinearTVM::getTVEdCorKirDT(q0,q1,T0,T); // update dCorKirdT
+ double mechanicalSource_cap = 0.; // in case of mullins, undamaged mechanicalSource
- // 1) ThermoElastic Heat -> _DcorKirDT needs to corrected for plasticity ( DcorKirDT = dcorKir_dT + dCorKirDEe * dEeDT
+ // 1) ThermoElastic Heat -> _DcorKirDT needs to corrected for plasticity ( DcorKirDT = dcorKir_dT + dCorKirDEe * dEeDT )
static STensor3 DEe;
DEe = q1->_Ee;
DEe -= q0->_Ee;
if (this->getTimeStep() > 0){
- mechanicalSource += (STensorOperation::doubledot(q1->_DcorKirDT,DEe)*T/this->getTimeStep());
+ // mechanicalSource += (STensorOperation::doubledot(q1->_DcorKirDT,DEe)*T/this->getTimeStep());
}
// 2) Viscoelastic Contribution to mechSrc
@@ -1748,27 +1768,70 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
STensor3& dWmdF_TVE = q1->getRefTodMechSrcTVEdF();
static STensor3 dWmdE_TVE;
mlawNonLinearTVM::getMechSource_nonLinearTVE_term(q0,q1,T0,T,Wm_TVE,stiff,dWmdE_TVE,dWmdT_TVE);
- mechanicalSource += Wm_TVE;
+ // mechanicalSource += Wm_TVE;
// 3) Viscoplastic Contribution to mechSrc
- getMechSourceTVP(F0,F,q0,q1,T0,T,Fepr,Cepr,DphiPDF,&Wm_TVP);
+ double& dWmdT_TVP = q1->getRefTodMechSrcTVPdT();
+ STensor3& dWmdF_TVP = q1->getRefTodMechSrcTVPdF();
+ mlawNonLinearTVP::getMechSourceTVP(F0,F,q0,q1,T0,T,Fepr,Cepr,DphiPDF,Wm_TVP,dWmdF_TVP,dWmdT_TVP);
mechanicalSource += Wm_TVP;
- // Compute mechanicalSource for Mullin's Effect
+ // --------------------------------------------------------------------------------------------
+ // TVE energy and derivatives
+ double& Dpsi_DT = q1->getRefTo_Dpsi_DT();
+ STensor3& Dpsi_TVEdE = q1->getRefTo_Dpsi_DE();
+ q1->_elasticEnergy = mlawNonLinearTVM::freeEnergyMechanical(*q0,*q1,T0,T,&Dpsi_DT,&Dpsi_TVEdE);
+
+ // TVP energy and derivatives
+ double psi_TVP = freeEnergyPlasticity(q0,q1,T0,T);
+ double Dpsi_TVPdT(0.);
+ static STensor3 Dpsi_TVPdF;
+ freeEnergyPlasticityDerivatives(q0,q1,T0,T,Dpsi_TVPdF,Dpsi_TVPdT);
+ // q1->_elasticEnergy += psi_TVP; // Mullins TVP
+ // Dpsi_DT += Dpsi_TVPdT; // add plastic part // Mullins TVP
+ Dpsi_DT -= dot(Dpsi_TVEdE,q1->_dGammaNdT); // add plastic dependency of elastic strain // Mullins TVP
+
+ // Mullins Effect
if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){
+ mlawNonLinearTVM::calculateMullinsEffectScaler(q0, q1, T, &Dpsi_DT);
+ }
+ double eta_mullins = 1.;
+ if (q1->_ipMullinsEffect != NULL){
+
+ // Put Stress in IP
+ q1->_P_cap = P;
+
+ // Modify Stress for mullins
+ eta_mullins = q1->_ipMullinsEffect->getEta();
+ P *= eta_mullins;
+ // q1->_elasticEnergy *= eta_mullins;
+
+ mechanicalSource_cap = mechanicalSource;
mechanicalSource *= q1->_ipMullinsEffect->getEta();
if (this->getTimeStep() > 0){
mechanicalSource -= q1->_ipMullinsEffect->getDpsiNew_DpsiMax()*(q1->_ipMullinsEffect->getpsiMax() - q0->_ipMullinsEffect->getpsiMax())/this->getTimeStep();
}
}
- // freeEnergy and elastic energy
- double& psiTVM = q1->getRefToFreeEnergyTVM();
- // getFreeEnergyTVM(q0,q1,T0,T,&psiTVM,NULL);
if(stiff){
+ // --------------------------------------------------------------------------------------------
+ // Important conversion tensors
+ STensorOperation::zero(DEeDCepr);
+ STensorOperation::zero(DEeDF);
+
+ // DEeDCepr to DEeDF
+ // STensorOperation::multSTensor43(_I4, DlnDCepr, DEeDCepr);
+ DEeDCepr = DlnDCepr;
+ DEeDCepr *= 0.5;
+ DEeDCepr -= dGammaNdCepr;
+ STensorOperation::multSTensor43(DEeDCepr, CeprToF, DEeDF);
+
+ // Note: DEeDT = - dGammaNdT
+ // --------------------------------------------------------------------------------------------
+
// thermSrc and MechSrc Derivatives
// thermal source derivatives
@@ -1814,27 +1877,88 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
}
}
}
- STensorOperation::multSTensor3STensor43(dmechanicalSourceE,DEeDF,dmechanicalSourceF);
+ // STensorOperation::multSTensor3STensor43(dmechanicalSourceE,DEeDF,dmechanicalSourceF);
- dmechanicalSourcedT += (STensorOperation::doubledot(q1->_DcorKirDT,DEe)/this->getTimeStep());
- dmechanicalSourcedT += T*(STensorOperation::doubledot(DDcorKirDTT,DEe)/this->getTimeStep());
+ // dmechanicalSourcedT += (STensorOperation::doubledot(q1->_DcorKirDT,DEe)/this->getTimeStep());
+ // dmechanicalSourcedT += T*(STensorOperation::doubledot(DDcorKirDTT,DEe)/this->getTimeStep());
+ // dmechanicalSourcedT -= T*(STensorOperation::doubledot(q1->_DcorKirDT,q1->_dGammaNdT)/this->getTimeStep());
}
// TVE Term
STensorOperation::multSTensor3STensor43(dWmdE_TVE, DEeDF, dWmdF_TVE);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
- dWmdT_TVE += - dWmdE_TVE(i,j) * q1->_dGammaNdT(i,j); // the first term in dWmdT_TVE from pure TVE temperature dependency is obtained above.
+ dWmdT_TVE += - dWmdE_TVE(i,j) * q1->_dGammaNdT(i,j); // the second term in dWmdT_TVE from dependency of Ee on plasticity
+
+ // dmechanicalSourceF += dWmdF_TVE;
+ // dmechanicalSourcedT += dWmdT_TVE;
+
// TVP Term
- double& dWmdT_TVP = q1->getRefTodMechSrcTVPdT();
- STensor3& dWmdF_TVP = q1->getRefTodMechSrcTVPdF();
- getMechSourceTVP(F0,F,q0,q1,T0,T,Fepr,Cepr,DphiPDF,NULL,&dWmdF_TVP,&dWmdT_TVP);
-
- dmechanicalSourceF += dWmdF_TVE;
- dmechanicalSourcedT += dWmdT_TVE;
- dmechanicalSourceF += dWmdF_TVP;
- dmechanicalSourcedT += dWmdT_TVP;
+ dmechanicalSourceF += dWmdF_TVP;
+ dmechanicalSourcedT += dWmdT_TVP;
+
+ // Mullins Effect Derivatives
+
+ // Compute Tangents for Mullin's Effect
+ if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){
+
+ // 1st Term
+ Tangent *= q1->_ipMullinsEffect->getEta();
+ dPdT *= q1->_ipMullinsEffect->getEta();
+ dmechanicalSourceF *= q1->_ipMullinsEffect->getEta();
+ dmechanicalSourcedT *= q1->_ipMullinsEffect->getEta();
+
+
+ // 2nd Term
+ double Deta_mullins_DT(0.);
+ static STensor3 Deta_mullins_DF, Deta_mullins_DEe;
+
+ // Deta_mullins_DF
+ STensorOperation::multSTensor3STensor43(q1->_DpsiDE,DEeDF,Deta_mullins_DF);
+ Deta_mullins_DF *= q1->_DmullinsDamage_Dpsi_cap; // elastic part
+ // Deta_mullins_DF += Dpsi_TVPdF*q1->_DmullinsDamage_Dpsi_cap; // plastic part // Mullins TVP
+
+ // Deta_mullins_DT
+ // STensorOperation::doubleContractionSTensor3(q1->_DpsiDE,q1->_dGammaNdT,Deta_mullins_DT); // Remove - Mullins TVP
+ // Deta_mullins_DT *= -q1->_DmullinsDamage_Dpsi_cap; // Remove - Mullins TVP
+ // Deta_mullins_DT += q1->_DmullinsDamage_DT; // Remove - Mullins TVP
+ Deta_mullins_DT = q1->_DmullinsDamage_DT; // + q1->_DmullinsDamage_Dpsi_cap*q1->_DpsiDT; // Mullins TVP
+
+ // 2nd Term - 1st PK
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dPdT(i,j) += q1->_P_cap(i,j) * Deta_mullins_DT;
+ for (int k=0; k<3; k++)
+ for (int l=0; l<3; l++)
+ Tangent(i,j,k,l) += q1->_P_cap(i,j) * Deta_mullins_DF(k,l);
+ }
+
+
+ // 2nd Term - MechSrc
+ dmechanicalSourcedT += mechanicalSource_cap*Deta_mullins_DT;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dmechanicalSourceF(i,j) += mechanicalSource_cap*Deta_mullins_DF(i,j);
+ }
+
+ static STensor3 DDpsiNew_DpsiMaxDF, DpsiMax_DF;
+ STensorOperation::multSTensor3STensor43(q1->_DpsiDE,DEeDF,DDpsiNew_DpsiMaxDF);
+ STensorOperation::multSTensor3STensor43(q1->_DpsiDE,DEeDF,DpsiMax_DF);
+ DDpsiNew_DpsiMaxDF *= q1->_ipMullinsEffect->getDDpsiNew_DpsiMaxDpsi();
+ DpsiMax_DF *= q1->_ipMullinsEffect->getDpsiMax_Dpsi();
+
+ if (this->getTimeStep() > 0){
+ dmechanicalSourcedT -= (q1->_ipMullinsEffect->getDDpsiNew_DpsiMaxDT()*(q1->_ipMullinsEffect->getpsiMax() - q0->_ipMullinsEffect->getpsiMax())/this->getTimeStep()
+ + q1->_ipMullinsEffect->getDpsiNew_DpsiMax()*q1->_ipMullinsEffect->getDpsiMax_Dpsi()*q1->_DpsiDT/this->getTimeStep());
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dmechanicalSourceF(i,j) -= DDpsiNew_DpsiMaxDF(i,j)*(q1->_ipMullinsEffect->getpsiMax() - q0->_ipMullinsEffect->getpsiMax())/this->getTimeStep();
+ dmechanicalSourceF(i,j) -= q1->_ipMullinsEffect->getDpsiNew_DpsiMax()*DpsiMax_DF(i,j)/this->getTimeStep();
+ }
+ }
+
+ } // mullin's
}
}
@@ -2068,9 +2192,11 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
// At every iteration it will calculate Ee, corKir and N
+ // Initialise GammaN in IP -> for mechSrc later
+ STensor3& GammaN = q1->getRefToGammaN();
+
// Initialise Hinv
- static STensor3 Hinv, GammaN;
- static STensor3 Ceinvpr;
+ static STensor3 Hinv, Ceinvpr;
STensorOperation::inverseSTensor3(expGN,Hinv);
STensorOperation::inverseSTensor3(Cepr,Ceinvpr);
@@ -2126,7 +2252,7 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
// Get Hb
double Hb(0.);
if (q1->_ipKinematic != NULL)
- Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
+ Hb = q1->_ipKinematic->getR(); // kinematic hardening parameter
double kk = (1./sqrt(1.+2.*q1->_nup*q1->_nup));
// Initialise backStress
@@ -2148,6 +2274,9 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
for (int j=0; j<3; j++)
J(i,j) += Phi(i,j) - MS(i,j) + q1->_backsig(i,j);
+ static STensor3 J0; //debug
+ J0 = J; //debug
+
// Initialise J_tol
double J_tol = 0.;
J_tol = J.norm0();
@@ -2177,7 +2306,8 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
int ite = 0;
int maxite = 100;
- double tol = 1e-9;
+ double tol = _stressIteratorTol; // 1e-6;
+
while (fabs(J_tol) > tol or ite <1){
// numerical Dho4inv (Jacobian dJ/dPhi)
@@ -2318,7 +2448,7 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
// update mandel
getModifiedMandel(Ce, q0, q1);
MS = q1->_ModMandel;
-
+
// update backStress
STensorOperation::decomposeDevTr(N,devN,trN);
devX = pow(kk,1)*Hb*Gamma*devN + devXn; // pow(kk,2) DEBUG
@@ -2352,6 +2482,25 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
}
}*/
+ // debug
+ /*
+ 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++)
+ dJdPhi_plus(k,l,i,j) = (J(k,l) - J0(k,l))/(_perturbationfactor);
+
+ STensorOperation::inverseSTensor43(dJdPhi_plus,dJdPhi_plus_inv);
+ J0 = J;*/
+
+ 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++){
+ // Msg::Error(" Dho2inv(%d,%d,%d,%d)=%d",i,j,k,l,Dho2inv(i,j,k,l));
+ // Msg::Error(" dJdPhi_plus_inv(%d,%d,%d,%d)=%d",i,j,k,l,dJdPhi_plus_inv(i,j,k,l));
+ }
+
J_tol = J.norm0();
J_tol/=(_K+_G);
// J_tol/=(1.e+10);
@@ -2364,8 +2513,9 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
if (fabs(J_tol) <_tol) break;
if(ite > maxite){
- Msg::Error("No convergence for iterated Phi mlawNonLinearTVENonLinearTVP nonAssociatedFlow iter = %d, J_tol = %e!!",ite,J_tol);
- Phi(0,0) = Phi(1,1) = Phi(2,2) = sqrt(-1.);
+ Msg::Error("No convergence for iterated Phi mlawNonLinearTVENonLinearTVP nonAssociatedFlow iter = %d, J_tol = %e!!",ite,J_tol);
+ // Msg::Error("Exiting the function Iterated Phi!!");
+ // Phi(0,0) = Phi(1,1) = Phi(2,2) = sqrt(-1.);
return;
}
} // while
@@ -2398,10 +2548,11 @@ void mlawNonLinearTVENonLinearTVP::TEST_predictorCorrector_TVP_nonAssociatedFlow
double Hb(0.), dHb(0.), ddHb(0.), dHbdT(0.), ddHbdgammadT(0.), ddHbddT(0.);
if (q1->_ipKinematic != NULL){
- Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
- dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter derivative (dHb/dgamma)
- dHbdT = Hb * q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
- ddHbddT = Hb * q1->_ipKinematic->getDDRDTT();
+ Hb = q1->_ipKinematic->getR(); // kinematic hardening parameter
+ dHb = q1->_ipKinematic->getDR(); // kinematic hardening parameter derivative (dHb/dgamma)
+ dHbdT = q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
+ ddHbddT = q1->_ipKinematic->getDDRDTT();
+ ddHbdgammadT = q1->_ipKinematic->getDDRDT();
}
double eta(0.),Deta(0.),DetaDT(0.);
@@ -2639,10 +2790,11 @@ void mlawNonLinearTVENonLinearTVP::TEST_predictorCorrector_TVP_nonAssociatedFlow
mlawNonLinearTVP::hardening(q0,q1,T);
mlawNonLinearTVP::getYieldCoefficients(q1,a);
if (q1->_ipKinematic != NULL){
- Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
- dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter derivative (dHb/dgamma ??)
- dHbdT = Hb * q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
- ddHbddT = Hb * q1->_ipKinematic->getDDRDTT();
+ Hb = q1->_ipKinematic->getR(); // kinematic hardening parameter
+ dHb = q1->_ipKinematic->getDR(); // kinematic hardening parameter derivative (dHb/dgamma ??)
+ dHbdT = q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
+ ddHbddT = q1->_ipKinematic->getDDRDTT();
+ ddHbdgammadT = q1->_ipKinematic->getDDRDT();
}
//a.print("a update");
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.h b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.h
index 981f958f0050801dd03f7826850831f351619b60..35e26d4bedf3070726ad9676bba1cae7b2b2f600 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.h
@@ -17,9 +17,12 @@
class mlawNonLinearTVENonLinearTVP : public mlawNonLinearTVP{
protected:
+ double _stressIteratorTol;
protected:
- virtual double freeEnergyPlasticity(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double T0, const double T) const;
+ virtual double freeEnergyPlasticity(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double& T0, const double& T) const;
+
+ virtual void freeEnergyPlasticityDerivatives(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double& T0, const double& T, STensor3& dPsy_TVPdF, double& dPsy_TVPdT) const;
virtual void getIterated_DPhi(const double& T0, const double& T, const IPNonLinearTVP *q0, IPNonLinearTVP *q1,
const double& Gamma, const double& Cxtr, const double& Cxdev,
@@ -137,7 +140,7 @@ class mlawNonLinearTVENonLinearTVP : public mlawNonLinearTVP{
public:
mlawNonLinearTVENonLinearTVP(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);
+ const bool matrixbyPerturbation = false, const double pert = 1e-8, const double stressIteratorTol = 1.e-9, const bool thermalEstimationPreviousConfig = true);
#ifndef SWIG
mlawNonLinearTVENonLinearTVP(const mlawNonLinearTVENonLinearTVP& src);
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
index e235103549c8c13c2e9b3a05870bfc8ab9d4b9e4..412ae54f97fa561f9c56149936fd8b353f0a2a95 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
@@ -190,13 +190,13 @@ void mlawNonLinearTVM::setViscoElasticData(const int i, const double Ei, const d
}
};
-void mlawNonLinearTVM::setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2){
+void mlawNonLinearTVM::setCorrectionsAllBranchesTVE(const int i, const double V0, const double V1, const double V2, const double D0, const double D1, const double D2){
if (i> _N or i<1)
Msg::Error("This setting is invalid %d > %d",i,_N);
else{
- _V1[i-1] = V1; _V2[i-1] = V2; _D1[i-1] = D1; _D2[i-1] = D2;
+ _V0[i-1] = V0; _V1[i-1] = V1; _V2[i-1] = V2; _D0[i-1] = D0; _D1[i-1] = D1; _D2[i-1] = D2;
- Msg::Info("setting: element=%d, V1 = %e, V2 = %e, D1 = %e, D2 = %e",i-1,V1,V2,D1,D2);
+ Msg::Info("setting: element=%d, V0 = %e, V1 = %e, V2 = %e, D0 = %e, D1 = %e, D2 = %e",i-1,V0,V1,V2,D0,D1,D2);
}
};
@@ -489,22 +489,23 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
// _intA, _intB have CTasymm
// Regularising function
- double m = _tensionCompressionRegularisation;
- double expmtr = exp(-m*(trEe-Eth));
- double sigmoid = 1/(1.+expmtr);
- double CTasymm = sigmoid + _compCorrection*(1.-sigmoid);
- double dCTasymmDtrE = m*expmtr/pow((1.+expmtr),2.)*(1.-_compCorrection);
- double psy_inf = (KT*0.5*(trEe-Eth)*(trEe-Eth) + GT*STensorOperation::doubledot(devEe,devEe) + KT*q._intA + 2.*GT*q._intB)/CTasymm;
+ // double m = _tensionCompressionRegularisation;
+ // double expmtr = exp(-m*(trEe-Eth));
+ // double sigmoid = 1/(1.+expmtr);
+ // double CTasymm = sigmoid + _compCorrection*(1.-sigmoid);
+ // double dCTasymmDtrE = m*expmtr/pow((1.+expmtr),2.)*(1.-_compCorrection);
+ // double psy_inf = (KT*0.5*(trEe-Eth)*(trEe-Eth) + GT*STensorOperation::doubledot(devEe,devEe) + KT*q._intA + 2.*GT*q._intB)/CTasymm;
+ // retains only the A and B integrands assuming integrands1 = integrands2 -> cant divide by CTasymm if integrands1 != integrands2
(*DpsiDT) += dKdT*q._intA + (KT*(q._elasticBulkPropertyScaleFactor) *(trEe-Eth) * (-3*dAlphaTdT*(T-_Tinitial)-3*AlphaT)) +
2.*dGdT*q._intB;
- (*DpsiDT) += psy_inf*dCTasymmDtrE*(-3*dAlphaTdT*(T-_Tinitial)-3*AlphaT);
+ (*DpsiDT) += (KT*q._DintA + 2*GT*q._DintB)*(-3*dAlphaTdT*(T-_Tinitial)-3*AlphaT); // psy_inf*dCTasymmDtrE*(-3*dAlphaTdT*(T-_Tinitial)-3*AlphaT);
(*DpsiDE) += q._corKirExtra; // already has CTasymm
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
- (*DpsiDE)(j,k) += psy_inf*dCTasymmDtrE*_I(j,k);
+ (*DpsiDE)(j,k) += (KT*q._DintA + 2*GT*q._DintB)*_I(j,k); //psy_inf*dCTasymmDtrE*_I(j,k);
}
}
@@ -610,12 +611,14 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
(*DpsiDT) += (_Ki[i]*q._Av_TVE_vector[i]*q._B[i]*dTrEei_DT +
2.*_Gi[i]*q._Bd_TVE_vector[i]*STensorOperation::doubledot(q._A[i],dDevEei_DT)); // * CTasymm;
- (*DpsiDT) += (_Ki[i]*q._intAv_TVE_vector[i] + 2*_Gi[i]*q._intBd_TVE_vector[i])/CTasymm * dCTasymmDtrE * dTrEei_DT;
+ // (*DpsiDT) += (_Ki[i]*q._intAv_TVE_vector[i] + 2*_Gi[i]*q._intBd_TVE_vector[i])/CTasymm * dCTasymmDtrE * dTrEei_DT;
+ (*DpsiDT) += (_Ki[i]*q._DintAv_TVE_vector[i] + 2*_Gi[i]*q._DintBd_TVE_vector[i]) * dTrEei_DT;
for (int j=0; j<3; j++)
for (int k=0; k<3; k++){
(*DpsiDE)(j,k) += _Ki[i]*q._Av_TVE_vector[i]*q._B[i]*exp_mid_k*_I(j,k); // * CTasymm;
- (*DpsiDE)(j,k) += (_Ki[i]*q._intAv_TVE_vector[i] + 2.*_Gi[i]*q._intBd_TVE_vector[i])/CTasymm * dCTasymmDtrE * exp_mid_k*_I(j,k);
+ // (*DpsiDE)(j,k) += (_Ki[i]*q._intAv_TVE_vector[i] + 2.*_Gi[i]*q._intBd_TVE_vector[i])/CTasymm * dCTasymmDtrE * exp_mid_k*_I(j,k);
+ (*DpsiDE)(j,k) += (_Ki[i]*q._DintAv_TVE_vector[i] + 2.*_Gi[i]*q._DintBd_TVE_vector[i])* exp_mid_k*_I(j,k);
for (int p=0; p<3; p++)
for (int r=0; r<3; r++)
(*DpsiDE)(j,k) += 2.*_Gi[i]*q._Bd_TVE_vector[i]*q._A[i](p,r)*exp_mid_g*_Idev(p,r,j,k); // * CTasymm;
@@ -676,13 +679,13 @@ void mlawNonLinearTVM::evaluateElasticCorrection(const double trE, const STensor
double &A_v, double &dA_vdE, double &intA, double &dAdT,
double &B_d, STensor3 &dB_vddev, double &intB, double &dBdT,
double *ddAdTT, double *ddBdTT, double *ddA_vdTdE, STensor3 *ddB_vdTdevE,
- double *dB_dTrEe, double *dB_dTdTrEe, double* psiInfCorrector) const{
+ double *dB_dTrEe, double *dB_dTdTrEe, double* psiInfCorrector, double* DintA, double* DintB) const{
if (_extraBranchNLType != TensionCompressionRegularisedType && _extraBranchNLType != hyper_exp_TCasymm_Type){
- mlawHyperViscoElastic::evaluatePhiPCorrection(trE,devE,A_v,dA_vdE,intA,B_d,dB_vddev,intB,NULL,psiInfCorrector);
+ mlawHyperViscoElastic::evaluatePhiPCorrection(trE,devE,A_v,dA_vdE,intA,B_d,dB_vddev,intB,NULL,psiInfCorrector,DintA,DintB);
}
else{
- mlawHyperViscoElastic::evaluatePhiPCorrection(trE,devE,A_v,dA_vdE,intA,B_d,dB_vddev,intB,dB_dTrEe,psiInfCorrector);
+ mlawHyperViscoElastic::evaluatePhiPCorrection(trE,devE,A_v,dA_vdE,intA,B_d,dB_vddev,intB,dB_dTrEe,psiInfCorrector,DintA,DintB);
// *dB_dTrEe *= _temFunc_elasticCorrection_Shear->getVal(T);
// if(dB_dTdTrEe!=NULL){
// *dB_dTdTrEe = (*dB_dTrEe)*_temFunc_elasticCorrection_Bulk->getDiff(T);
@@ -745,15 +748,15 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
double A(0.), B(0.);
double dA_dTrEe(0.), dA_dTdTrEe(0.), dB_dTrEe(0.), dB_dTdTrEe(0.), psiInfCorrector(0.);
STensor3 dB_dDevEe, ddB_dTdDevEe;
- double intA(0.), intB(0.), dA_dT(0.), dB_dT(0.), ddA_dTT(0.), ddB_dTT(0.);
+ double intA(0.), intB(0.), dA_dT(0.), dB_dT(0.), ddA_dTT(0.), ddB_dTT(0.), DintA(0.), DintB(0.);
if (_extraBranchNLType != TensionCompressionRegularisedType && _extraBranchNLType != hyper_exp_TCasymm_Type){
evaluateElasticCorrection(trEe, devEe, T, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, &ddA_dTT, &ddB_dTT, &dA_dTdTrEe, &ddB_dTdDevEe,
- NULL, NULL, &psiInfCorrector);
+ NULL, NULL, &psiInfCorrector, &DintA, &DintB);
}
else {
evaluateElasticCorrection(trEe, devEe, T, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, &ddA_dTT, &ddB_dTT, &dA_dTdTrEe, &ddB_dTdDevEe,
- &dB_dTrEe, &dB_dTdTrEe, &psiInfCorrector);
+ &dB_dTrEe, &dB_dTdTrEe, &psiInfCorrector, &DintA, &DintB);
}
if (A <= -1. + 1.e-5){ // saturated
@@ -762,6 +765,7 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
intA += A*0.5 *(trEe*trEe-trEe_0*trEe_0);
}
q1->_intA = intA;
+ q1->_DintA = DintA;
if (B <= -1. + 1.e-5){ // saturated
static STensor3 devEe_0;
@@ -771,6 +775,7 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
intA += B*0.5 *(devEe.dotprod()-devEe_0.dotprod());
}
q1->_intB = intB;
+ q1->_DintB = DintB;
q1->_psiInfCorrector = psiInfCorrector; // Additional term to correct thermal expansion term in free energy
@@ -888,7 +893,7 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
void mlawNonLinearTVM::extraBranch_nonLinearTVE(const int i, const STensor3& Ee, const double& T,
double& Av, double& dA_dTrEe, double& intA, double& intA1, double& intA2, double& dA_dT, double& ddA_dTT, double& ddA_dTdTrEe,
double& Bd, STensor3& dB_dDevEe, double &intB, STensor3 &intB1, STensor3 &intB2, double &dB_dT, double &ddB_dTT, STensor3& ddB_dTdDevEe,
- double* dB_dTrEe) const{
+ double& DintA, double& DintB, double* dB_dTrEe) const{
// Similar code to previous implementations of extraBranch except we dont do -1 + Av, have to make this a unique implementation.
@@ -904,6 +909,7 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const int i, const STensor3& Ee,
// int(1/Bd*dB_dDevEe dDevEe) = Gi*intB2 // For 2nd convolution integral and stiffness --> actually it is, int(1/Gi * dGi/dDevEe dDevEe)
// ddA_dTdTrEe - mechSrc??
// ddB_dTdDevEe - mechSrc??
+ // DintA, DintB - the additional derivatives of the sigmoid in the free energy approximation associated with tension-compression assymmetry
static STensor3 dev;
double tr;
@@ -1032,6 +1038,9 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const int i, const STensor3& Ee,
integrand -= ( 1./getXiVolumeCorrection()*sqrt(getZetaVolumeCorrection()) ); // value at trEe = 0.
intA = sigmoid*integrand + _compCorrection*(1.-sigmoid)*integrand;
intA *= getVolumeCorrection();
+
+ DintA = (m*expmtr/pow((1.+expmtr),2.))*integrand - _compCorrection*(m*expmtr/pow((1.+expmtr),2.))*integrand;
+ DintA *= getVolumeCorrection();
}
else{
intA = 1.;
@@ -1056,6 +1065,9 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const int i, const STensor3& Ee,
integrand -= (1./getThetaDevCorrection()*sqrt(getPiDevCorrection()) ); // value at devEe = 0.
intB = sigmoid*integrand + _compCorrection*(1.-sigmoid)*integrand;
intB *= getDevCorrection();
+
+ DintB = (m*expmtr/pow((1.+expmtr),2.))*integrand - _compCorrection*(m*expmtr/pow((1.+expmtr),2.))*integrand;
+ DintB *= getDevCorrection();
}
else{
intB = 1.;
@@ -1073,46 +1085,51 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const int i, const STensor3& Ee,
double sigmoid = 1/(1.+expmtr);
// Av
- double x = _V1[i]*tr*tr + _V2[i];
- Av = sigmoid*1./sqrt(x) + (1.-sigmoid)*(_Ci[i]/sqrt(x));
- Av *= getVolumeCorrection();
+ double x1 = _V1[i]*tr*tr + _V2[i];
+ double x2 = _V4[i]*tr*tr + _V5[i];
+ Av = sigmoid*(1./sqrt(x1)+_V0[i]) + (1.-sigmoid)*_Ci[i]*(1/sqrt(x2)+_V3[i]);
- dA_dTrEe = - sigmoid*_V1[i]*tr/pow(x,1.5) - (1.-sigmoid)*(_V1[i]*_Ci[i]*tr/pow(x,1.5))
- + (m*expmtr/pow((1.+expmtr),2.))*1./sqrt(x) - (_Ci[i]*m*expmtr/pow((1.+expmtr),2.))*1./sqrt(x);
- dA_dTrEe *= getVolumeCorrection();
+ dA_dTrEe = - sigmoid*_V1[i]*tr/pow(x1,1.5) - (1.-sigmoid)*(_V4[i]*_Ci[i]*tr/pow(x2,1.5))
+ + (m*expmtr/pow((1.+expmtr),2.))*(1./sqrt(x1)+_V0[i]) - (_Ci[i]*m*expmtr/pow((1.+expmtr),2.))*(1/sqrt(x2)+_V3[i]);
- if(_V1[i]>0.){
- double integrand = 1./_V1[i]*sqrt(x); // integral of A_v * trEe
- integrand -= ( 1./_V1[i]*sqrt(_V2[i]) ); // value at trEe = 0.
- intA = sigmoid*integrand + _Ci[i]*(1.-sigmoid)*integrand;
- intA *= getVolumeCorrection();
+ if(_V1[i]>0. && _V4[i]>0.){
+ double integrand1 = 1./_V1[i]*sqrt(x1) + _V0[i]*tr*tr/2.; // integral of A_v * trEe
+ double integrand2 = 1./_V4[i]*sqrt(x2) + _V3[i]*tr*tr/2.;
+ integrand1 -= ( 1./_V1[i]*sqrt(_V2[i]) ); // value at trEe = 0.
+ integrand2 -= ( 1./_V4[i]*sqrt(_V5[i]) ); // value at trEe = 0.
+ intA = sigmoid*integrand1 + _Ci[i]*(1.-sigmoid)*integrand2;
+
+ DintA = (m*expmtr/pow((1.+expmtr),2.))*integrand1 - _Ci[i]*(m*expmtr/pow((1.+expmtr),2.))*integrand2;
}
else{
intA = 1.;
+ // Msg::Error(" Inside extraBranch_nonLinearTVE Type 4, _V1 = 0. or _V4 = 0., incompatibility with Mullin's effect");
}
// Bd
- double y = _D1[i]*dev.dotprod() + _D2[i];
- Bd = sigmoid*1./sqrt(y) + (1.-sigmoid)*(_Ci[i]/sqrt(y));
- Bd *= getDevCorrection();
+ double y1 = _D1[i]*dev.dotprod() + _D2[i];
+ double y2 = _D4[i]*dev.dotprod() + _D5[i];
+ Bd = sigmoid*(1./sqrt(y1)+_D0[i]) + (1.-sigmoid)*_Ci[i]*(1/sqrt(y2)+_D3[i]);
STensorOperation::zero(dB_dDevEe);
dB_dDevEe = dev;
- dB_dDevEe *= (-sigmoid*_D1[i]/pow(y,1.5) - (1.-sigmoid)*(_D1[i]*_Ci[i]/pow(y,1.5)));
- dB_dDevEe *= getDevCorrection();
+ dB_dDevEe *= (-sigmoid*_D1[i]/pow(y1,1.5) - (1.-sigmoid)*(_D4[i]*_Ci[i]/pow(y2,1.5)));
if(dB_dTrEe!=NULL){
- *dB_dTrEe = (m*expmtr/pow((1.+expmtr),2.))*1./sqrt(y) - (_Ci[i]*m*expmtr/pow((1.+expmtr),2.))*1./sqrt(y);
- *dB_dTrEe *= getDevCorrection();
+ *dB_dTrEe = ( (m*expmtr/pow((1.+expmtr),2.))*(1./sqrt(y1)+_D0[i]) - (_Ci[i]*m*expmtr/pow((1.+expmtr),2.))*(1/sqrt(y2)+_D3[i]) );
}
- if(_D1[i]>0.){
- double integrand = 1./_D1[i]*sqrt(y); // integral of B_d * devEe
- integrand -= (1./_D1[i]*sqrt(_D2[i]) ); // value at devEe = 0.
- intB = sigmoid*integrand + _Ci[i]*(1.-sigmoid)*integrand;
- intB *= getDevCorrection();
+ if(_D1[i]>0. && _D4[i]>0.){
+ double integrand1 = 1./_D1[i]*sqrt(y1) + _D0[i]*dev.dotprod()/2.; // integral of B_d * devEe
+ double integrand2 = 1./_D4[i]*sqrt(y2) + _D3[i]*dev.dotprod()/2.;
+ integrand1 -= ( 1./_D1[i]*sqrt(_D2[i]) ); // value at devEe = 0.
+ integrand2 -= ( 1./_D4[i]*sqrt(_D5[i]) ); // value at devEe = 0.
+ intB = sigmoid*integrand1 + _Ci[i]*(1.-sigmoid)*integrand2;
+
+ DintB = (m*expmtr/pow((1.+expmtr),2.))*integrand1 - _Ci[i]*(m*expmtr/pow((1.+expmtr),2.))*integrand2;
}
else{
intB = 1.;
+ // Msg::Error(" Inside extraBranch_nonLinearTVE Type 4, _D1 = 0. or _D4 = 0., incompatibility with Mullin's effect");
}
}
/*else if (_ExtraBranch_TVE_option == 5){
@@ -1186,7 +1203,6 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const int i, const STensor3& Ee,
// Regularising function
double m = _tensionCompressionRegularisation;
double expmtr = exp(-m*tr);
- // if (exp(-m*tr)<1.e+10){ expmtr = exp(-m*tr);}
double sigmoid = 1/(1.+expmtr);
// Av
@@ -1208,6 +1224,8 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const int i, const STensor3& Ee,
integrand1 -= ( 1./_V1[i]*sqrt(_V2[i]) ); // value at trEe = 0.
intA = sigmoid*(integrand1 + integrand2) + _Ci[i]*(1.-sigmoid)*integrand1;
// intA *= getVolumeCorrection();
+
+ DintA = (m*expmtr/pow((1.+expmtr),2.))*integrand1 - _Ci[i]*(m*expmtr/pow((1.+expmtr),2.))*integrand2;
}
else{
intA = 1.;
@@ -1237,6 +1255,8 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const int i, const STensor3& Ee,
integrand1 -= (1./_D1[i]*sqrt(_D2[i]) ); // value at devEe = 0.
intB = sigmoid*(integrand1 + integrand2) + _Ci[i]*(1.-sigmoid)*integrand1;
// intB *= getDevCorrection();
+
+ DintB = (m*expmtr/pow((1.+expmtr),2.))*integrand1 - _Ci[i]*(m*expmtr/pow((1.+expmtr),2.))*integrand2;
}
else{
intB = 1.;
@@ -1331,20 +1351,22 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
double A(0.), B(0.), dA_dTrEe(0.), psiInfCorrector(0.);
static STensor3 dB_dDevEe;
- double intA(0.), intB(0.), dA_dT(0.), dB_dT(0.);
+ double intA(0.), intB(0.), dA_dT(0.), dB_dT(0.), DintA(0.), DintB(0.);
// Here, A, B and related terms first calculated as (1+A) and (1+B)
if(_extraBranchNLType != TensionCompressionRegularisedType && _extraBranchNLType != hyper_exp_TCasymm_Type){
- evaluateElasticCorrection(eff_trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, NULL ,NULL, NULL, NULL, NULL, &psiInfCorrector);
+ evaluateElasticCorrection(eff_trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, NULL ,NULL, NULL, NULL, NULL, &psiInfCorrector, &DintA, &DintB);
}
else{
double dB_dTrEe(0.);
- evaluateElasticCorrection(eff_trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, NULL ,NULL, NULL, NULL, &dB_dTrEe, &psiInfCorrector);
+ evaluateElasticCorrection(eff_trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, NULL ,NULL, NULL, NULL, &dB_dTrEe, &psiInfCorrector, &DintA, &DintB);
q1->_dBd_dTrEe = dB_dTrEe;
}
q1->_intA = intA;
q1->_intB = intB;
+ q1->_DintA = DintA;
+ q1->_DintB = DintB;
q1->_psiInfCorrector = psiInfCorrector; // Additional term to correct thermal expansion term in free energy
q1->_elasticBulkPropertyScaleFactor = A;
q1->_elasticShearPropertyScaleFactor = B;
@@ -1391,7 +1413,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
// Do non-linear TVE extraBranch here
double Av(1.), Bd(1.);
- double dA_dTrEe(0.),intA(0.),intA1(0.),intA2(0.),dA_dT(0.),ddA_dTT(0.),ddA_dTdTrEe(0.),intB(0.),dB_dT(0.),ddB_dTT(0.);
+ double dA_dTrEe(0.),intA(0.),intA1(0.),intA2(0.),dA_dT(0.),ddA_dTT(0.),ddA_dTdTrEe(0.),intB(0.),dB_dT(0.),ddB_dTT(0.),DintA(0.),DintB(0.);
static STensor3 intB1,intB2,dB_dDevEe,ddB_dTdDevEe;
STensorOperation::zero(intB1);
STensorOperation::zero(intB2);
@@ -1405,7 +1427,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
// NEW
if(_ExtraBranch_TVE_option == 1){
- mlawNonLinearTVM::extraBranch_nonLinearTVE(0,Ee,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe);
+ mlawNonLinearTVM::extraBranch_nonLinearTVE(0,Ee,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,DintA,DintB);
q1->_Av_TVE = Av; q1->_Bd_TVE = Bd; q1->_dAv_dTrEe_TVE = dA_dTrEe; q1->_dBd_dDevEe_TVE = dB_dDevEe;
q1->_intAv_TVE = intA; q1->_intBd_TVE = intB; q1->_intAv_1_TVE = intA1; q1->_intBd_1_TVE = intB1; // this is zero
q1->_intAv_2_TVE = intA2; q1->_intBd_2_TVE = intB2; // this is zero
@@ -1509,12 +1531,12 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
branchElasticStrain(2,2) += 1./3.*q1->_B[i];
if (_ExtraBranch_TVE_option == 2){
mlawNonLinearTVM::extraBranch_nonLinearTVE(i,branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
- Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe);
+ Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,DintA,DintB);
}
else{
double dB_dTrEe(0.);
mlawNonLinearTVM::extraBranch_nonLinearTVE(i,branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
- Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,&dB_dTrEe);
+ Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,DintA,DintB,&dB_dTrEe);
q1->_dBd_dTrEe_TVE_vector[i] = dB_dTrEe;
}
@@ -1524,6 +1546,8 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
q1->_dBd_dDevEe_TVE_vector[i] = dB_dDevEe;
q1->_intAv_TVE_vector[i] = intA;
q1->_intBd_TVE_vector[i] = intB;
+ q1->_DintAv_TVE_vector[i] = DintA;
+ q1->_DintBd_TVE_vector[i] = DintB;
// Evaluate Bd_stiffnessTerm and/or Bd_stiffnessTerm2
for (int k=0; k<3; k++)
@@ -1694,6 +1718,11 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
for (int q=0; q<3; q++){
DDA_DTDdevE[i](k,l,r,s) += exp_mid_g * q1->_dBd_dDevEe_TVE_vector[i](p,q)*dDevEei_DT(p,q)*_I4(k,l,r,s);
}
+
+ if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
+ DDA_DTDdevE[i](k,l,r,s) += exp_mid_g*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT*_I4(k,l,r,s);
+ }
+
}
}
} // if - extraBranchTVE_option
@@ -1952,6 +1981,9 @@ void mlawNonLinearTVM::getMechSource_nonLinearTVE_term(const IPNonLinearTVM *q0,
double T_set = setTemp(T0,T,1); // 1->T1, 2->T_mid, 3->T_midshift
double AlphaT0, AlphaT, DAlphaDT;
getAlpha(AlphaT,T_set,&DAlphaDT); getAlpha(AlphaT0,T0);
+ double eff_trEe = trEe - 3.*(AlphaT*(T-_Tinitial));
+ double eff_trEe0 = trEe0 - 3.*(AlphaT0*(T0-_Tinitial));
+ double Deff_trEe_DT = -3.*(DAlphaDT*(T-_Tinitial) + AlphaT);
double eff_trDE = q1->_trDE - 3.*(AlphaT*(T-_Tinitial) - AlphaT0*(T0-_Tinitial));
double Deff_trDE_DT = -3.*(DAlphaDT*(T-_Tinitial) + AlphaT);
// NEW
@@ -2010,12 +2042,12 @@ void mlawNonLinearTVM::getMechSource_nonLinearTVE_term(const IPNonLinearTVM *q0,
for (int l=0; l<3; l++){
Qi[i](k,l) = _Ki[i]*q1->_Av_TVE_vector[i]*q1->_B[i]*_I(k,l) + 2*_Gi[i]*q1->_Bd_TVE_vector[i]*q1->_A[i](k,l);
DQiDT[i](k,l) = _Ki[i]*q1->_DB_DT[i]*_I(k,l) + 2*_Gi[i]*q1->_DA_DT[i](k,l);
- dot_Gammai[i](k,l) = 1/3*(trEe-q1->_B[i] - trEe0 + q0->_B[i])*_I(k,l) + (devEe(k,l)-q1->_A[i](k,l) - devEe0(k,l) + q0->_A[i](k,l) );
+ dot_Gammai[i](k,l) = 1/3*(eff_trEe - q1->_B[i] - eff_trEe0 + q0->_B[i])*_I(k,l) + (devEe(k,l)-q1->_A[i](k,l) - devEe0(k,l) + q0->_A[i](k,l) ); // NEW
}
}
if (this->getTimeStep() > 0){
dot_Gammai[i] *= 1/this->getTimeStep();
- Wm += STensorOperation::doubledot(Qi[i],dot_Gammai[i]); // - T*STensorOperation::doubledot(DQiDT[i],dot_Gammai[i]);
+ Wm += STensorOperation::doubledot(Qi[i],dot_Gammai[i]); // - T*STensorOperation::doubledot(DQiDT[i],dot_Gammai[i]); // OLD
}
// TVE Mechanical Source Term derivatives
@@ -2051,7 +2083,7 @@ void mlawNonLinearTVM::getMechSource_nonLinearTVE_term(const IPNonLinearTVM *q0,
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
DDQiDTT[i](k,l) = _Ki[i]*q1->_DDB_DTT[i]*_I(k,l) + 2*_Gi[i]*q1->_DDA_DTT[i](k,l);
- dot_DGammaiDT[i](k,l) = - 1/3*dTrEei_DT*_I(k,l) - dDevEei_DT(k,l);
+ dot_DGammaiDT[i](k,l) = 1/3*(Deff_trEe_DT-dTrEei_DT)*_I(k,l) - dDevEei_DT(k,l); // NEW
for (int r=0; r<3; r++){
for (int s=0; s<3; s++){
DQiDEe[i](k,l,r,s) = _Ki[i]*q1->_DB_DtrE[i]*_I(k,l)*_I(r,s);
@@ -2080,15 +2112,15 @@ void mlawNonLinearTVM::getMechSource_nonLinearTVE_term(const IPNonLinearTVM *q0,
// dWmdT += (STensorOperation::doubledot(Qi[i],dot_DGammaiDT[i]) - T*STensorOperation::doubledot(DQiDT[i],dot_DGammaiDT[i])
- // - T*STensorOperation::doubledot(DDQiDTT[i],dot_Gammai[i]));
+ // - T*STensorOperation::doubledot(DDQiDTT[i],dot_Gammai[i])); // OLD
- dWmdT += (STensorOperation::doubledot(DQiDT[i][i],dot_Gammai[i]) + STensorOperation::doubledot(Qi[i],dot_DGammaiDT[i]));
+ dWmdT += (STensorOperation::doubledot(DQiDT[i],dot_Gammai[i]) + STensorOperation::doubledot(Qi[i],dot_DGammaiDT[i]));
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
for (int r=0; r<3; r++)
for (int s=0; s<3; s++){
- // dWmdE(r,s) += ((DQiDEe[i](k,l,r,s) - T*DDQiDTDEe[i](k,l,r,s))*dot_Gammai[i](k,l) + (Qi[i](k,l) - T*DQiDT[i](k,l))*dot_DGammaiDEe[i](k,l,r,s));
+ // dWmdE(r,s) += ((DQiDEe[i](k,l,r,s) - T*DDQiDTDEe[i](k,l,r,s))*dot_Gammai[i](k,l) + (Qi[i](k,l) - T*DQiDT[i](k,l))*dot_DGammaiDEe[i](k,l,r,s)); // OLD
dWmdE(r,s) += (DQiDEe[i](k,l,r,s)*dot_Gammai[i](k,l) + Qi[i](k,l)*dot_DGammaiDEe[i](k,l,r,s));
}
}
@@ -2482,7 +2514,7 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
// Modify Stress for mullins
eta_mullins = q1->_ipMullinsEffect->getEta();
P *= eta_mullins;
- q1->_elasticEnergy *= eta_mullins;
+ // q1->_elasticEnergy *= eta_mullins;
}
// Msg::Error(" Inside TVE predCorr, after updating q1->_elasticEnergy = %e!!",q1->_elasticEnergy);
@@ -2544,6 +2576,7 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
// mechanicalSource
mechanicalSource = 0.;
+ double mechanicalSource_cap = 0.; // in case of mullins, undamaged mechanicalSource
// 1) ThermoElastic Heat
if (this->getTimeStep() > 0){
@@ -2564,11 +2597,11 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
dWmdT_TVE *= 1.;
q1->_mechSrcTVE = Wm_TVE;
}
- mechanicalSource += Wm_TVE;
+ // mechanicalSource += Wm_TVE; // TVE term
// Compute mechanicalSource for Mullin's Effect
if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){
-
+ mechanicalSource_cap = mechanicalSource;
mechanicalSource *= q1->_ipMullinsEffect->getEta();
if (this->getTimeStep() > 0){
mechanicalSource -= q1->_ipMullinsEffect->getDpsiNew_DpsiMax()*(q1->_ipMullinsEffect->getpsiMax() - q0->_ipMullinsEffect->getpsiMax())/this->getTimeStep();
@@ -2740,7 +2773,7 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
static const STensor43 I4_2(1.,1.);
static STensor3 DcorKirDT_I;
STensorOperation::multSTensor3STensor43(DcorKirDT,I4_2,DcorKirDT_I);
-
+
if (this->getTimeStep() > 0){
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
@@ -2751,12 +2784,12 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
}
}
}
- dmechanicalSourceE += dWmdE_TVE;
+ // dmechanicalSourceE += dWmdE_TVE; // TVE term
STensorOperation::multSTensor3STensor43(dmechanicalSourceE,DEeDFe,dmechanicalSourceF);
dmechanicalSourcedT += (STensorOperation::doubledot(DcorKirDT,DEe)/this->getTimeStep());
dmechanicalSourcedT += T*(STensorOperation::doubledot(DDcorKirDT,DEe)/this->getTimeStep());
- dmechanicalSourcedT += dWmdT_TVE;
+ // dmechanicalSourcedT += dWmdT_TVE; // TVE term
}
// Compute Tangents for Mullin's Effect
@@ -2772,10 +2805,14 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
static STensor3 Deta_mullins_DF;
STensorOperation::multSTensor3STensor43(q1->_DpsiDE,DEeDFe,Deta_mullins_DF);
Deta_mullins_DF *= q1->_DmullinsDamage_Dpsi_cap;
+
+ double Deta_mullins_DT(0.);
+ Deta_mullins_DT = q1->_DmullinsDamage_DT; // + q1->_DmullinsDamage_Dpsi_cap*q1->_DpsiDT;
+
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
- dPdT(i,j) += q1->_P_cap(i,j) * q1->_DmullinsDamage_DT;
- dmechanicalSourceF(i,j) += mechanicalSource*Deta_mullins_DF(i,j);
+ dPdT(i,j) += q1->_P_cap(i,j) * Deta_mullins_DT;
+ dmechanicalSourceF(i,j) += mechanicalSource_cap*Deta_mullins_DF(i,j);
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
Tangent(i,j,k,l) += q1->_P_cap(i,j) * Deta_mullins_DF(k,l);
@@ -2787,7 +2824,7 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
DDpsiNew_DpsiMaxDF *= q1->_ipMullinsEffect->getDDpsiNew_DpsiMaxDpsi();
DpsiMax_DF *= q1->_ipMullinsEffect->getDpsiMax_Dpsi();
- dmechanicalSourcedT += q1->_DmullinsDamage_DT*mechanicalSource;
+ dmechanicalSourcedT += mechanicalSource_cap*Deta_mullins_DT; // q1->_DmullinsDamage_DT;
if (this->getTimeStep() > 0){
dmechanicalSourcedT -= (q1->_ipMullinsEffect->getDDpsiNew_DpsiMaxDT()*(q1->_ipMullinsEffect->getpsiMax() - q0->_ipMullinsEffect->getpsiMax())/this->getTimeStep()
+ q1->_ipMullinsEffect->getDpsiNew_DpsiMax()*q1->_ipMullinsEffect->getDpsiMax_Dpsi()*q1->_DpsiDT/this->getTimeStep());
@@ -2853,19 +2890,21 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
double A(0.), B(0.), dA_dTrEe(0.), psiInfCorrector(0.);
static STensor3 dB_dDevEe;
- double intA(0.), intB(0.), dA_dT(0.), dB_dT(0.);
+ double intA(0.), intB(0.), dA_dT(0.), dB_dT(0.), DintA(0.), DintB(0.);
if(_extraBranchNLType != TensionCompressionRegularisedType && _extraBranchNLType != hyper_exp_TCasymm_Type){
- evaluateElasticCorrection(eff_trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, NULL ,NULL, NULL, NULL, NULL, &psiInfCorrector);
+ evaluateElasticCorrection(eff_trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, NULL ,NULL, NULL, NULL, NULL, &psiInfCorrector, &DintA, &DintB);
}
else{
double dB_dTrEe(0.);
- evaluateElasticCorrection(eff_trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, NULL ,NULL, NULL, NULL, &dB_dTrEe, &psiInfCorrector);
+ evaluateElasticCorrection(eff_trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, NULL ,NULL, NULL, NULL, &dB_dTrEe, &psiInfCorrector, &DintA, &DintB);
q1->_dBd_dTrEe = dB_dTrEe;
}
q1->_intA = intA;
q1->_intB = intB;
+ q1->_DintA = DintA;
+ q1->_DintB = DintB;
q1->_psiInfCorrector = psiInfCorrector; // Additional term to correct thermal expansion term in free energy
q1->_elasticBulkPropertyScaleFactor = A;
q1->_elasticShearPropertyScaleFactor = B;
@@ -2912,7 +2951,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
// Do non-linear TVE extraBranch here
double Av(1.), Bd(1.);
- double dA_dTrEe(0.),intA(0.),intA1(0.),intA2(0.),dA_dT(0.),ddA_dTT(0.),ddA_dTdTrEe(0.),intB(0.),dB_dT(0.),ddB_dTT(0.);
+ double dA_dTrEe(0.),intA(0.),intA1(0.),intA2(0.),dA_dT(0.),ddA_dTT(0.),ddA_dTdTrEe(0.),intB(0.),dB_dT(0.),ddB_dTT(0.),DintA(0.),DintB(0.);
static STensor3 intB1,intB2,dB_dDevEe,ddB_dTdDevEe;
STensorOperation::zero(intB1);
STensorOperation::zero(intB2);
@@ -2921,7 +2960,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
if (_useExtraBranch_TVE){
if(_ExtraBranch_TVE_option == 1){
- mlawNonLinearTVM::extraBranch_nonLinearTVE(0,Ee,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe);
+ mlawNonLinearTVM::extraBranch_nonLinearTVE(0,Ee,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,DintA,DintB);
q1->_Av_TVE = Av; q1->_Bd_TVE = Bd; q1->_dAv_dTrEe_TVE = dA_dTrEe; q1->_dBd_dDevEe_TVE = dB_dDevEe;
q1->_intAv_TVE = intA; q1->_intBd_TVE = intB;
q1->_intAv_1_TVE = intA1; q1->_intBd_1_TVE = intB1; // this is zero
@@ -3014,12 +3053,12 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
branchElasticStrain(2,2) += 1./3.*q1->_B[i];
if (_ExtraBranch_TVE_option == 2){
mlawNonLinearTVM::extraBranch_nonLinearTVE(i,branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
- Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe);
+ Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,DintA,DintB);
}
else{
double dB_dTrEe(0.);
mlawNonLinearTVM::extraBranch_nonLinearTVE(i,branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
- Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,&dB_dTrEe);
+ Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,DintA,DintB,&dB_dTrEe);
q1->_dBd_dTrEe_TVE_vector[i] = dB_dTrEe;
}
@@ -3029,7 +3068,9 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
q1->_dBd_dDevEe_TVE_vector[i] = dB_dDevEe;
q1->_intAv_TVE_vector[i] = intA;
q1->_intBd_TVE_vector[i] = intB;
-
+ q1->_DintAv_TVE_vector[i] = DintA;
+ q1->_DintBd_TVE_vector[i] = DintB;
+
// Evaluate Bd_stiffnessTerm
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
index 75ab05a0bba921e48a93c3c33b59abbf9d176b3c..245bf095e0800c64167a3399bab7340f385fc247 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
@@ -100,7 +100,7 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
double &B_d, STensor3 &dB_vddev, double &intB, double &dBdT,
double *ddAdTT = NULL, double *ddBdTT = NULL,
double *ddA_vdTdE = NULL, STensor3 *ddB_vdTdevE = NULL,
- double *dB_dTrEe = NULL, double *dB_dTdTrEe = NULL, double *psiInfCorrector = NULL) const;
+ double *dB_dTrEe = NULL, double *dB_dTdTrEe = NULL, double *psiInfCorrector = NULL, double* DintA = NULL, double* DintB = NULL) const;
virtual void extraBranchLaw(const STensor3& Ee, const double& T, const IPNonLinearTVM *q0, IPNonLinearTVM *q1,
STensor3& sig, const bool stiff = false, STensor43* DsigDEe = NULL, STensor3* DsigDT = NULL,
@@ -113,8 +113,8 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
virtual void extraBranch_nonLinearTVE(const int i, const STensor3& Ee, const double& T,
double& Av, double& dA_dTrEe, double& intA, double& intA1, double& intA2, double& dA_dT, double& ddA_dTT, double& ddA_dTdTrEe,
double& Bd, STensor3& dB_dDevEe, double &intB, STensor3 &intB1, STensor3 &intB2, double &dB_dT, double &ddB_dTT, STensor3& ddB_dTdDevEe,
- double* dB_dTrEe = NULL) const;
-
+ double& DintA, double& DintB, double* dB_dTrEe = NULL) const;
+
void corKirInfinity(const IPNonLinearTVM *q1, const STensor3& devEe, const double& trEe, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const;
void ThermoViscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
@@ -202,7 +202,7 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
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 setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2);
+ virtual void setCorrectionsAllBranchesTVE(const int i, const double V0, const double V1, const double V2, const double D0, const double D1, const double D2);
virtual void setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double D3);
virtual void setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
virtual void setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci);
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index e45d224a9fd182848c587a32018f2a3480410fba..8c38474fb3228fbdaa1fe44c6cb93371ab87fa3e 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -298,7 +298,7 @@ void mlawNonLinearTVP::getFreeEnergyTVM(const IPNonLinearTVP *q0, IPNonLinearTVP
// get Kinematic Hardening stuff
double Hb(0.), dHbdT(0.), ddHbddT(0.);
if (q1->_ipKinematic != NULL){
- Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
+ Hb = q1->_ipKinematic->getR(); // kinematic hardening parameter
dHbdT = q1->_ipKinematic->getDRDT();
ddHbddT = q1->_ipKinematic->getDDRDTT();
}
@@ -435,6 +435,7 @@ void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNo
// dRdT
dRdT = dsigcdT - dsigc0dT;
+ dRdT += Hc*q1->_DgammaDT;
// dRdF
STensorOperation::zero(dRdF);
@@ -445,12 +446,12 @@ void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNo
// ddRdTT
if(ddRdTT!=NULL){
- *ddRdTT = ddsigcdTT - ddsigc0dTT;
+ *ddRdTT = (ddsigcdTT - ddsigc0dTT) + dHcdT*DgammaDT; // DDgammaDTT = 0.
}
// ddRdTF
if(ddRdFdT!=NULL){
- STensor3 _ddRdFdT;
+ static STensor3 _ddRdFdT;
STensorOperation::zero(_ddRdFdT);
for (int i=0; i<3; i++)
@@ -463,7 +464,7 @@ void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNo
void mlawNonLinearTVP::getMechSourceTVP(const STensor3& F0, const STensor3& F,
const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double T0, const double& T, const STensor3& Fepr, const STensor3& Cepr,
- const STensor3& DphiPDF, double *Wm, STensor3 *dWmdF, double *dWmdT) const{
+ const STensor3& DphiPDF, double& mechSourceTVP, STensor3& dmechSourceTVPdF, double& dmechSourceTVPdT) const{
const double& DgammaDT = q1->_DgammaDT;
const STensor3& DgammaDF = q1->_DgammaDF;
@@ -471,105 +472,51 @@ void mlawNonLinearTVP::getMechSourceTVP(const STensor3& F0, const STensor3& F,
const double& gamma1 = q1->_epspbarre;
const STensor3& Fp1 = q1->_Fp;
const STensor3& Me = q1->_ModMandel;
- const STensor3& X1 = q1->_backsig;
- const STensor3& X0 = q0->_backsig;
- const STensor3& dXdT = q1->_DbackSigDT; // For debugging = 0. //DEBUGGING
- const STensor3& dXdT_0 = q0->_DbackSigDT;
- const STensor43& dXdF_0 = q0->_DbackSigDF;
+ const STensor3& X = q1->_backsig;
+ const STensor3& dXdT = q1->_DbackSigDT;
const STensor43& dXdF = q1->_DbackSigDF;
const STensor3& dMeDT = q1->_DModMandelDT;
const STensor43& dMedF = q1->_DModMandelDF;
const STensor43& dGammaNdF = q1->_dGammaNdF;
const STensor3& dGammaNdT = q1->_dGammaNdT;
+ const STensor43& DphiDF = q1->_DphiDF;
+ const STensor3& DphiDT = q1->_DphiDT;
// get Dgamma
double Dgamma = gamma1 - gamma0;
- // get Hb - Kinematic hardening parameter
- double Hb(0.), dHb(0), dHbdT(0.), ddHbddT(0.);
- if (q1->_ipKinematic != NULL){
- Hb = q1->_ipKinematic->getDR();
- dHb = q1->_ipKinematic->getDDR();
- dHbdT = q1->_ipKinematic->getDRDT();
- ddHbddT = q1->_ipKinematic->getDDRDTT();
- }
-
// get Dp - dont use Dp = Fp_dot*Fpinv, USE Dp = Gamma*N/this->getTimeStep()
- static STensor3 Dp;
- STensorOperation::zero(Dp);
- if (this->getTimeStep() > 0){
- Dp = q1->_GammaN;
- Dp *= 1/this->getTimeStep();
- }
-
- // get alpha
- double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
- static STensor3 alpha0, alpha1;
- alpha0 = X0; alpha1 = X1;
- if (Hb>0.){
- alpha0 *= 1/(pow(kk,2)* (q0->_ipKinematic->getDR())) ; // using previous Hb for alpha0
- alpha1 *= 1/(pow(kk,2)*Hb);
- }
- alpha1 -= alpha0;
+ const STensor3& Dp = q1->_GammaN;
// get TempX - convenient backStress tensor
- STensor3 tempX;
- tempX = X1;
- tempX -= T*dXdT;
+ static STensor3 tempX;
+ tempX = X;
+ // tempX -= T*dXdT;
// get R, dRdT, dRdF
double R(0.), dRdT(0.), ddRdTT(0.);
- STensor3 dRdF, ddRdFdT; STensorOperation::zero(dRdF); STensorOperation::zero(ddRdFdT);
-
- // Old Formulation
- /*
- std::vector<double> ddIsoHardForcedTT; ddIsoHardForcedTT.resize(3,0.);
- std::vector<STensor3> ddIsoHardForcedFdT; ddIsoHardForcedFdT.resize(3,0.);
-
- getIsotropicHardeningForce(q0,q1,T0,T,DphiPDF,&ddIsoHardForcedTT,&ddIsoHardForcedFdT); // For analytical derivatives; for numerical, put the last 2 args as NULL
- const std::vector<double>& IsoHardForce = q1->_IsoHardForce;
- const std::vector<double>& dIsoHardForcedT = q1->_dIsoHardForcedT;
- const std::vector<STensor3>& dIsoHardForcedF = q1->_dIsoHardForcedF;
-
- for (int i=0; i<3; i++){
- R += IsoHardForce[i]
- dRdT += dIsoHardForcedT[i];
- ddRdTT += ddIsoHardForcedTT[i]; // 0.;
- dRdF += dIsoHardForcedF[i];
- ddRdFdT += ddIsoHardForcedFdT[i]; // 0.;
- }*/
+ static STensor3 dRdF, ddRdFdT;
+ STensorOperation::zero(dRdF); STensorOperation::zero(ddRdFdT);
// New Formulation
-
- double ddIsoHardForcedTT(0.);
- STensor3 ddIsoHardForcedFdT;
- STensorOperation::zero(ddIsoHardForcedFdT);
-
- getIsotropicHardeningForce(q0,q1,T0,T,&ddIsoHardForcedTT,&ddIsoHardForcedFdT);
+ getIsotropicHardeningForce(q0,q1,T0,T,&ddRdTT,&ddRdFdT);
R = q1->_IsoHardForce_simple;
dRdT = q1->_dIsoHardForcedT_simple;
dRdF = q1->_dIsoHardForcedF_simple;
// for convenience
- R -= T*dRdT;
+ // R -= T*dRdT;
// mechSourceTVP
- if(Wm!=NULL){
- double mechSourceTVP = 0.;
-
- mechSourceTVP += STensorOperation::doubledot(Me,Dp);
- if (!_TaylorQuineyFlag){
- if (this->getTimeStep() > 0){
- alpha1 *= 1/this->getTimeStep();
- mechSourceTVP -= STensorOperation::doubledot(tempX,alpha1);
- mechSourceTVP -= (R*Dgamma)/this->getTimeStep();
- }
- else{
- mechSourceTVP *= _TaylorQuineyFactor;
- }
+ mechSourceTVP = 0.;
+ if (this->getTimeStep() > 0){
+ mechSourceTVP += STensorOperation::doubledot(Me,Dp);
+ mechSourceTVP -= STensorOperation::doubledot(X,Dp); // dot_alpha = Dp
+ mechSourceTVP *= 1/this->getTimeStep();
- *Wm = mechSourceTVP;
- }
+ if (!_TaylorQuineyFlag){
+ mechSourceTVP -= (R*Dgamma)/this->getTimeStep();
+ }
}
// dDpdF and dDpdT
@@ -577,98 +524,45 @@ void mlawNonLinearTVP::getMechSourceTVP(const STensor3& F0, const STensor3& F,
// dDpdT = dGammaNdT/this->getTimeStep();
- // dAlphadF - from dXdF
- static STensor43 dAlphadF;
- STensorOperation::zero(dAlphadF);
- if (Hb>0.){
- 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++){
- dAlphadF(i,j,k,l) += dXdF(i,j,k,l)*1/(pow(kk,2)*Hb) - X1(i,j)*DgammaDF(k,l)*1/(pow(kk*Hb,2))*dHb;
- }
- }
-
- // dAlphadT - from dXdT
- static STensor3 dAlphadT;
- if (Hb>0.){
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- dAlphadT(i,j) = dXdT(i,j)*1/(pow(kk,2)*Hb) - X1(i,j)*DgammaDT*1/(pow(kk*Hb,2))*dHbdT;
- }
- }
-
- // ddXdTF -> make it -- CHANGE!!
- static STensor43 ddXdTdF; // ddGammaDTF and ddQDTF are very difficult to estimate
- static STensor3 ddXdTT; // ddGammaDTT and ddQDTT are very difficult to estimate
- STensorOperation::zero(ddXdTdF);
- STensorOperation::zero(ddXdTT);
- double delT = T-T0;
- if (delT>1e-10){
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- // ddXdTT(i,j) += 2*dXdT(i,j)/(T-T0) - 2*(X1(i,j)-X0(i,j))/pow((T-T0),2); //DEBUGGING
- for (int k=0; k<3; k++)
- for (int l=0; l<3; l++){
- // ddXdTdF(i,j,k,l) += (dXdF(i,j,k,l)-dXdF_0(i,j,k,l))/(T-T0); //DEBUGGING
- }
- }
- }
-
- if(dWmdF!=NULL){
- STensor3 dmechSourceTVPdF(0.);
- if (this->getTimeStep() > 0){
+ STensorOperation::zero(dmechSourceTVPdF);
+ dmechSourceTVPdT = 0.;
+ if (this->getTimeStep() > 0){
+ 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++){
+ dmechSourceTVPdF(k,l) += DphiDF(i,j,k,l)*Dp(i,j) + (Me(i,j)-X(i,j))*dGammaNdF(i,j,k,l)/this->getTimeStep();
+ // dmechSourceTVPdF(k,l) += dMedF(i,j,k,l)*Dp(i,j) + Me(i,j)*dGammaNdF(i,j,k,l)/this->getTimeStep();
+ // dmechSourceTVPdF(k,l) += - (dXdF(i,j,k,l) - T*ddXdTdF(i,j,k,l))*Dp(i,j) - tempX(i,j)*dGammaNdF(i,j,k,l)/this->getTimeStep();
+ // dmechSourceTVPdF(k,l) += - (dXdF(i,j,k,l))*Dp(i,j) - tempX(i,j)*dGammaNdF(i,j,k,l)/this->getTimeStep();
+ }
- 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++){
- dmechSourceTVPdF(i,j) += dMedF(i,j,k,l)*Dp(i,j) + Me(i,j)*dGammaNdF(i,j,k,l)/this->getTimeStep();
- }
- }
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ dmechSourceTVPdT += DphiDT(i,j)*Dp(i,j) + (Me(i,j)-X(i,j))*dGammaNdT(i,j)/this->getTimeStep();
+ // dmechSourceTVPdT += dMeDT(i,j)*Dp(i,j) + Me(i,j)*dGammaNdT(i,j)/this->getTimeStep();
+ // dmechSourceTVPdT += T*ddXdTT(i,j)*Dp(i,j) - tempX(i,j)*dGammaNdT(i,j)/this->getTimeStep();// old
+ // dmechSourceTVPdT -= dXdT(i,j)*Dp(i,j) - tempX(i,j)*dGammaNdT(i,j)/this->getTimeStep();
- if (!_TaylorQuineyFlag){
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- dmechSourceTVPdF(i,j) = - (dRdF(i,j) - T*ddRdFdT(i,j))*Dgamma/this->getTimeStep() - R*DgammaDF(i,j)/this->getTimeStep();
- for (int k=0; k<3; k++)
- for (int l=0; l<3; l++){
- dmechSourceTVPdF(i,j) += - (dXdF(i,j,k,l) - T*ddXdTdF(i,j,k,l))*alpha1(i,j) - tempX(i,j)*dAlphadF(i,j,k,l)/this->getTimeStep();
- }
- }
- }
- else{
- dmechSourceTVPdF *= _TaylorQuineyFactor;
- }
+ }
+ if (!_TaylorQuineyFlag){
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++){
+ // dmechSourceTVPdF(i,j) += - (dRdF(i,j) - T*ddRdFdT(i,j))*Dgamma/this->getTimeStep() - R*DgammaDF(i,j)/this->getTimeStep(); // old
+ dmechSourceTVPdF(i,j) += - dRdF(i,j)*Dgamma/this->getTimeStep() - R*DgammaDF(i,j)/this->getTimeStep(); // new
}
- *dWmdF = dmechSourceTVPdF;
- }
-
- if(dWmdT!=NULL){
- double dmechSourceTVPdT(0.);
- if (this->getTimeStep() > 0){
-
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- dmechSourceTVPdT += dMeDT(i,j)*Dp(i,j) + Me(i,j)*dGammaNdT(i,j)/this->getTimeStep();
- }
-
- if (!_TaylorQuineyFlag){
-
- dmechSourceTVPdT += T*ddRdTT*Dgamma/this->getTimeStep() - R*DgammaDT/this->getTimeStep();
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++){
- dmechSourceTVPdT += T*ddXdTT(i,j)*alpha1(i,j) - tempX(i,j)*dAlphadT(i,j)/this->getTimeStep();
- }
- }
- else{
- dmechSourceTVPdT *= _TaylorQuineyFactor;
- }
- }
- *dWmdT = dmechSourceTVPdT;
+ // dmechSourceTVPdT += T*ddRdTT*Dgamma/this->getTimeStep() - R*DgammaDT/this->getTimeStep(); // old
+ dmechSourceTVPdT += - dRdT*Dgamma/this->getTimeStep() - R*DgammaDT/this->getTimeStep(); // new
+ }
+ else{
+ dmechSourceTVPdF *= _TaylorQuineyFactor;
+ dmechSourceTVPdT *= _TaylorQuineyFactor;
+ }
}
+
+ double debug_CHECKER = 0.;
};
@@ -982,10 +876,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
double Hb(0.), dHb(0.), ddHb(0.), dHbdT(0.), ddHbdgammadT(0.), ddHbddT(0.);
if (q1->_ipKinematic != NULL){
- Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
- dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter derivative (dHb/dgamma)
- dHbdT = Hb * q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
- ddHbddT = Hb * q1->_ipKinematic->getDDRDTT();
+ Hb = q1->_ipKinematic->getR(); // kinematic hardening parameter
+ dHb = q1->_ipKinematic->getDR(); // kinematic hardening parameter derivative (dHb/dgamma)
+ dHbdT = q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
+ ddHbddT = q1->_ipKinematic->getDDRDTT();
}
double eta(0.),Deta(0.),DetaDT(0.);
@@ -1266,10 +1160,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
hardening(q0,q1,T);
getYieldCoefficients(q1,a);
if (q1->_ipKinematic != NULL){
- Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
- dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter derivative (dHb/dgamma ??)
- dHbdT = Hb * q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
- ddHbddT = Hb * q1->_ipKinematic->getDDRDTT();
+ Hb = q1->_ipKinematic->getR(); // kinematic hardening parameter
+ dHb = q1->_ipKinematic->getDR(); // kinematic hardening parameter derivative (dHb/dgamma ??)
+ dHbdT = q1->_ipKinematic->getDRDT(); // make sure to normalise DRDT while defining in python file //NEW
+ ddHbddT = q1->_ipKinematic->getDDRDTT();
}
//a.print("a update");
@@ -1509,6 +1403,8 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
// mechanical Source
double& Wm_TVP = q1->getRefToMechSrcTVP(); // TVP
double& Wm_TVE = q1->getRefToMechSrcTVE(); // TVE
+ double& dWmdT_TVP = q1->getRefTodMechSrcTVPdT();
+ STensor3& dWmdF_TVP = q1->getRefTodMechSrcTVPdF();
mechanicalSource = 0.;
@@ -1517,7 +1413,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
DEe = q1->_Ee;
DEe -= q0->_Ee;
if (this->getTimeStep() > 0){
- mechanicalSource += (STensorOperation::doubledot(q1->_DcorKirDT,DEe)*T/this->getTimeStep());
+ mechanicalSource += (STensorOperation::doubledot(q1->_DcorKirDT,DEe)*T/this->getTimeStep());
}
// 2) Viscoelastic Contribution to mechSrc
@@ -1529,7 +1425,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
mechanicalSource += Wm_TVE;
// 3) Viscoplastic Contribution to mechSrc
- getMechSourceTVP(F0,F,q0,q1,T0,T,Fepr,Cepr,DphiPDF,&Wm_TVP);
+ getMechSourceTVP(F0,F,q0,q1,T0,T,Fepr,Cepr,DphiPDF,Wm_TVP,dWmdF_TVP,dWmdT_TVP);
mechanicalSource += Wm_TVP;
// freeEnergy and elastic energy
@@ -2547,10 +2443,6 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
dWmdT_TVE += - dWmdE_TVE(i,j) * dGammaNdT(i,j); // the first term in dWmdT_TVE from pure TVE temperature dependency is obtained above.
// TVP Term
- double& dWmdT_TVP = q1->getRefTodMechSrcTVPdT();
- STensor3& dWmdF_TVP = q1->getRefTodMechSrcTVPdF();
- getMechSourceTVP(F0,F,q0,q1,T0,T,Fepr,Cepr,DphiPDF,NULL,&dWmdF_TVP,&dWmdT_TVP);
-
dmechanicalSourceF += dWmdF_TVE;
dmechanicalSourcedT += dWmdT_TVE;
dmechanicalSourceF += dWmdF_TVP; // dWmdF_TVP;
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.h b/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
index 0c86d67d788c77f9f65a3903d13f655620886d11..5ec7f1f49ef1b0a3baeabe24558ed46c2b5a1455 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
@@ -48,7 +48,7 @@ class mlawNonLinearTVP : public mlawNonLinearTVM{
virtual void getMechSourceTVP(const STensor3& F0, const STensor3& F, const IPNonLinearTVP *q0, IPNonLinearTVP *q1,
const double T0, const double& T, const STensor3& Fepr, const STensor3& Cepr,
- const STensor3& DphiPDF, double *Wm = NULL, STensor3 *dWmdF = NULL, double *dWmdT = NULL) const;
+ const STensor3& DphiPDF, double& mechSourceTVP, STensor3& dmechSourceTVPdF, double& dmechSourceTVPdT) const;
virtual void getFreeEnergyTVM(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double& T0, const double& T,
double *psiTVM = NULL, double *DpsiTVMdT = NULL, double *DDpsiTVMdTT = NULL) const;
@@ -1373,4 +1373,4 @@ if (stiff){
} // stiff
};
-*/
\ No newline at end of file
+*/
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp b/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp
index deb62d95705e968f4fe4d667cd4f735b96b85d65..af333371280f9d2e3307ac501c4939d3eb4d33f1 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp
@@ -1192,7 +1192,8 @@ void mlawNonLocalPorosity::predictorCorrector_NonLocalPorosity(
static STensor3 DFpDT;
static double dpprDT;
-
+
+ ipvcur->getRefToConstitutiveSuccessFlag() = true;
// Compute kcor, Fp and fv with an elastic predictor and eventually a plastic corrector
// Elastic predictor
bool okElast = elasticPredictor(F1,*Fp0,corKirpr,kcorprDev,ppr,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,Fepr,Cepr,Epr,Lepr,dLepr,stiff,
@@ -1200,6 +1201,7 @@ void mlawNonLocalPorosity::predictorCorrector_NonLocalPorosity(
if (!okElast)
{
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
return;
}
double kCorPrEq = sqrt(1.5*kcorprDev.dotprod());
@@ -1310,6 +1312,7 @@ void mlawNonLocalPorosity::predictorCorrector_NonLocalPorosity(
if (!okElast)
{
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
return;
}
//check yield condition
@@ -1496,6 +1499,7 @@ void mlawNonLocalPorosity::predictorCorrector_NonLocalPorosity(
&Tcur,&DcorKirprDT,&DkcorprDevDT,&dpprDT);
if (!okElast)
{
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
return;
}
@@ -1586,6 +1590,7 @@ void mlawNonLocalPorosity::predictorCorrector_NonLocalPorosity(
if (plastic and !correctorOK){
// In case of failed plastic correction
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
P(0,0) = sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
return;
}
@@ -2102,11 +2107,13 @@ void mlawNonLocalPorosity::predictorCorrectorLocal(
bool correctorOK = false;
bool plastic = false;
+ ipvcur->getRefToConstitutiveSuccessFlag() = true;
// elastic predictor
bool okElast = this->elasticPredictor(F1,*Fp0,corKirpr,kcorprDev,ppr,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,Fepr,Cepr,Epr,Lepr,dLepr,stiff);
if (!okElast)
{
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
return;
}
@@ -2197,6 +2204,7 @@ void mlawNonLocalPorosity::predictorCorrectorLocal(
okElast = elasticPredictor(Fcur,*Fp0,corKirpr,kcorprDev,ppr,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,Fepr,Cepr,Epr,Lepr,dLepr,estimateStiff);
if (!okElast)
{
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
return;
}
@@ -2265,6 +2273,7 @@ void mlawNonLocalPorosity::predictorCorrectorLocal(
}
}
if (plastic and !correctorOK){
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
P(0,0) = sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
return;
}
@@ -5580,12 +5589,14 @@ void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const ST
DFpDNonLocalVar.resize(_numNonLocalVar);
DLocalVarDEpr.resize(_numNonLocalVar);
};
-
+
+ ipvcur->getRefToConstitutiveSuccessFlag() = true;
// elastic predictor
bool okElast = elasticPredictor(F1,*Fp0,corKirpr,kcorprDev,ppr,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,Fepr,Cepr,Epr,Lepr,dLepr,stiff,
&T,&DcorKirprDT,&DkcorprDevDT,&dpprDT);
if (!okElast)
{
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
return;
}
@@ -5703,6 +5714,7 @@ void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const ST
if (!okElast)
{
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
return;
}
//check yield condition
@@ -5873,6 +5885,7 @@ void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const ST
&Tcur,&DcorKirprDT,&DkcorprDevDT,&dpprDT);
if (!okElast)
{
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
return;
}
@@ -5950,6 +5963,7 @@ void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const ST
if (plastic and !correctorOK){
P(0,0) = sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
return;
}
@@ -7797,6 +7811,8 @@ void mlawNonLocalPorosity::I1J2J3_predictorCorrector_NonLocal(const STensor3& F0
DFpDNonLocalVar.resize(_numNonLocalVar);
DLocalVarDEpr.resize(_numNonLocalVar);
};
+
+ ipvcur->getRefToConstitutiveSuccessFlag() = true;
// elastic predictor
bool okElast = elasticPredictor(F1,*Fp0,corKirpr,kcorprDev,ppr,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,Fepr,Cepr,Epr,Lepr,dLepr,true,
@@ -7804,6 +7820,7 @@ void mlawNonLocalPorosity::I1J2J3_predictorCorrector_NonLocal(const STensor3& F0
if (!okElast)
{
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
return;
}
// elastic only
@@ -7896,6 +7913,7 @@ void mlawNonLocalPorosity::I1J2J3_predictorCorrector_NonLocal(const STensor3& F0
&Tcur,&DcorKirprDT,&DkcorprDevDT,&dpprDT);
if (!okElast)
{
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
return;
}
@@ -8029,6 +8047,7 @@ void mlawNonLocalPorosity::I1J2J3_predictorCorrector_NonLocal(const STensor3& F0
&Tcur,&DcorKirprDT,&DkcorprDevDT,&dpprDT);
if (!okElast)
{
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
return;
}
@@ -8088,6 +8107,7 @@ void mlawNonLocalPorosity::I1J2J3_predictorCorrector_NonLocal(const STensor3& F0
if (plastic and !correctorOK){
P(0,0) = sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
return;
}
@@ -8273,12 +8293,14 @@ void mlawNonLocalPorosity::I1J2J3_predictorCorrectorLocal(
bool correctorOK = false;
bool plastic = false;
+ ipvcur->getRefToConstitutiveSuccessFlag() = true;
// elastic predictor
bool okElast = this->elasticPredictor(F1,*Fp0,corKirpr,kcorprDev,ppr,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,Fepr,Cepr,Epr,Lepr,dLepr,true);
if (!okElast)
{
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
return;
}
// elastic only
@@ -8348,6 +8370,7 @@ void mlawNonLocalPorosity::I1J2J3_predictorCorrectorLocal(
okElast = elasticPredictor(Fcur,*Fp0,corKirpr,kcorprDev,ppr,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,Fepr,Cepr,Epr,Lepr,dLepr,estimateStiff);
if (!okElast)
{
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
P(0,0) = P(1,1) = P(2,2)= sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
return;
}
@@ -8399,6 +8422,7 @@ void mlawNonLocalPorosity::I1J2J3_predictorCorrectorLocal(
}
if (plastic and !correctorOK){
P(0,0) = sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
+ ipvcur->getRefToConstitutiveSuccessFlag() = false;
return;
}
}
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp
index e5e99c20b488206ddd17104e22d3fa224338493b..dc35d65bf5c5a9eb8975e27b8649acf567e9a634 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp
@@ -2369,4 +2369,294 @@ void mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::I1J2J3_getYieldNormal(STen
STensorOperation::multSTensor3STensor43(DNpEffDT, _aniM, *DNpDT);
}
}
-};
\ No newline at end of file
+};
+
+
+void mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const
+{
+ // check condition if check coalescence is carried out
+ if (getNumOfYieldSurfaces() > 1)
+ {
+ Msg::Error("check coalescenece has not implemented for num yield surface >1");
+ return;
+ }
+
+ // Get ipvs
+ const IPCoalescence* q0Thom = &q0->getConstRefToIPCoalescence();
+ IPCoalescence* q1Thom = &q1->getRefToIPCoalescence();
+
+ // check active
+ if (q0Thom->getCoalescenceOnsetFlag() and q1->dissipationIsActive())
+ {
+ q1Thom->getRefToCoalescenceActiveFlag() = true;
+ }
+ else
+ {
+ q1Thom->getRefToCoalescenceActiveFlag() = false;
+ }
+
+ bool willUseNormalToCheck = false;
+ if (_checkCoalescenceWithNormal)
+ {
+ if (q1->getLocation() == IPVariable::INTERFACE_MINUS or q1->getLocation() == IPVariable::INTERFACE_PLUS)
+ {
+
+ willUseNormalToCheck = true;
+ const SVector3& normal = q1->getConstRefToCurrentOutwardNormal();
+ if (normal.norm() <=0.)
+ {
+ Msg::Error("coalescence cannot check with zero normal!!!");
+ return;
+ };
+ }
+ else
+ {
+ willUseNormalToCheck = false;
+ if (!_withBothYieldSurfacesInBulk)
+ {
+ return;
+ }
+ }
+ }
+
+
+ if (q0Thom->getCoalescenceOnsetFlag())
+ {
+ // if onset already occurs
+ q1Thom->getRefToCoalescenceOnsetFlag() = q0Thom->getCoalescenceOnsetFlag();
+ q1Thom->getRefToCoalescenceMode() = q0Thom->getCoalescenceMode();
+ // copy data
+ q1Thom->setIPvAtCoalescenceOnset(*q0Thom->getIPvAtCoalescenceOnset());
+
+ q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft();
+ q1Thom->getRefToLodeParameterOffset() = q0Thom->getLodeParameterOffset();
+ q1Thom->getRefToYieldOffset() = q0Thom->getYieldOffset();
+ }
+ else if (q1->dissipationIsActive())
+ {
+ const STensor3& kCor = q1->getConstRefToCorotationalKirchhoffStress();
+ const STensor3& F = q1->getConstRefToDeformationGradient();
+ const double R = q1->getCurrentViscoplasticYieldStress();
+ double J = STensorOperation::determinantSTensor3(F);
+ static STensor3 corSig;
+ corSig = kCor;
+ if (_stressFormulation == CORO_CAUCHY)
+ {
+ corSig *= (1./J);
+ }
+
+ // check at interface
+ bool interfaceNeckActive = true;
+ bool interfaceShearActive = true;
+ if (willUseNormalToCheck)
+ {
+ double CftThomason;
+ _mlawCoales->getCoalescenceLaw()->computeConcentrationFactor(q0,q1,CftThomason);
+ double CftShear;
+ _mlawShear->getCoalescenceLaw()->computeConcentrationFactor(q0,q1,CftShear);
+
+ SVector3 normal = q1->getConstRefToCurrentOutwardNormal();
+ normal.normalize();
+
+ // Compute Sig from Kcor following: Sig = invFeT * sigCor * FeT = Fe*sigCor*invFe
+ static STensor3 invFep, FeSigCor , Fe, Sig;
+ STensorOperation::inverseSTensor3(q1->getConstRefToFp(), invFep);
+ STensorOperation::multSTensor3(F, invFep, Fe);
+ STensorOperation::inverseSTensor3(Fe, invFep);
+ STensorOperation::multSTensor3(Fe, corSig, FeSigCor);
+ STensorOperation::multSTensor3SecondTranspose(FeSigCor, invFep, Sig);
+
+ // traction vector
+ static SVector3 tractionForce;
+ STensorOperation::multSTensor3SVector3(Sig,normal,tractionForce);
+ double tauN = STensorOperation::dot(tractionForce,normal);
+ // Deduce tangent part
+ double tauTanSq = 0.;
+ for (int i=0; i<3; i++)
+ {
+ double temp = tractionForce(i)-tauN*normal(i);
+ tauTanSq += temp*temp;
+ }
+ // Compute effective stress
+ double tauEff = 0.;
+ if (tauN > 0.)
+ {
+ tauEff = sqrt(tauN*tauN + _beta*tauTanSq);
+ }
+ else
+ {
+ tauEff = sqrt(_beta*tauTanSq);
+ }
+
+ if (tauEff - CftThomason*R >0)
+ {
+ interfaceNeckActive = true;
+ }
+ else
+ {
+ interfaceNeckActive = false;
+ }
+ if (sqrt(3.*tauTanSq) - CftShear*R >0)
+ {
+ interfaceShearActive = true;
+ }
+ else
+ {
+ interfaceShearActive = false;
+ }
+ }
+
+ static STensor3 sigEff;
+ STensorOperation::multSTensor43STensor3(_aniM, corSig, sigEff);
+ double Gurson = _mlawGrowth->I1J2J3_yieldFunction(sigEff,R,q0,q1,T);
+ double Thomason = _mlawCoales->I1J2J3_yieldFunction(sigEff,R,q0,q1,T);
+ double Shear = _mlawShear->I1J2J3_yieldFunction(sigEff,R,q0,q1,T);
+
+ // active plasticity
+ // check coalescence is active
+ bool& CoalesFlag = q1Thom->getRefToCoalescenceOnsetFlag();
+ int& CoalesMode = q1Thom->getRefToCoalescenceMode();
+ CoalesFlag = false;
+ CoalesMode = 0;
+
+ if (_useTwoYieldRegularization)
+ {
+ double FG = Gurson +1.;
+ double FT = Thomason +1.;
+ double FS = Shear +1.;
+ if (FG < 0. or FT < 0. or FS < 0)
+ {
+ Msg::Error("FT,FS FG must be all positive but FT=%e, FS = %e, FG=%e",FT,FS,FG);
+ }
+ // check necking
+ if (std::max(FT,FS)> FG)
+ {
+ double ratio = FG/FT;
+ if (pow(ratio,_yieldRegExponent) < _coalesTol && interfaceNeckActive)
+ {
+ CoalesFlag = true;
+ CoalesMode = 1;
+ printf("Necking coalescence at Gurson = %e, Thomason =%e, Shear = %e\n",
+ Gurson,Thomason,Shear);
+ }
+ else
+ {
+ ratio = FG/FS;
+ if (pow(ratio,_yieldRegExponent) < _coalesTol && interfaceShearActive)
+ {
+ CoalesFlag = true;
+ CoalesMode = 2;
+ printf("Shear coalescence at Gurson = %e, Thomason =%e, Shear = %e\n",
+ Gurson,Thomason,Shear);
+ }
+ else
+ {
+ CoalesFlag = false;
+ CoalesMode = 0;
+ }
+ }
+ }
+ }
+ else
+ {
+ if (Thomason > std::max(Gurson,Shear) && interfaceNeckActive)
+ {
+ CoalesFlag = true;
+ CoalesMode = 1;
+ printf("Necking coalescence at Gurson = %e, Thomason =%e, Shear = %e\n",
+ Gurson,Thomason,Shear);
+ }
+ else if (Shear > std::max(Gurson,Thomason) && interfaceShearActive)
+ {
+ CoalesFlag = true;
+ CoalesMode = 2;
+ printf("Shear coalescence at Gurson = %e, Thomason =%e, Shear = %e\n",
+ Gurson,Thomason,Shear);
+ }
+ else
+ {
+ CoalesFlag = false;
+ CoalesMode = 0;
+ }
+ }
+
+ // if coalesnce occurs
+ if (CoalesFlag)
+ {
+ q1Thom->setIPvAtCoalescenceOnset(*q1);
+ // if offset is used to satisfy Thomason yield surface at the onset of coalescence
+ if (_useTwoYieldRegularization)
+ {
+ // never offset with yield regulation
+ q1Thom->getRefToCrackOffsetOnCft() = 1.;
+ q1Thom->getRefToLodeParameterOffset() = 1.;
+ q1Thom->getRefToYieldOffset() = 0.;
+ }
+ else
+ {
+ if (_withCftOffset)
+ {
+ if (CoalesMode==1)
+ {
+ q1Thom->getRefToCrackOffsetOnCft() = 1. + Thomason;
+ }
+ else if (CoalesMode==2)
+ {
+ q1Thom->getRefToCrackOffsetOnCft() = 1. + Shear;
+ }
+ q1Thom->getRefToLodeParameterOffset() = 1.;
+ q1Thom->getRefToYieldOffset() = 0.;
+ #ifdef _DEBUG
+ printf("coalescence ocurrs using Cft offset = %e\n",q1Thom->getCrackOffsetOnCft());
+ #endif //_DEBUG
+ }
+ else if (_withLodeOffset)
+ {
+ // Lode offset is not used for material law, because law already accounts for Lode variable
+ q1Thom->getRefToCrackOffsetOnCft() = 1.;
+ q1Thom->getRefToLodeParameterOffset()= 1.;
+ q1Thom->getRefToYieldOffset() = 0.;
+ }
+ else if (_withYieldOffset)
+ {
+ q1Thom->getRefToCrackOffsetOnCft() = 1.;
+ q1Thom->getRefToLodeParameterOffset() = 1.;
+ if (CoalesMode==1)
+ {
+ q1Thom->getRefToYieldOffset() = Thomason;
+ }
+ else if (CoalesMode==2)
+ {
+ q1Thom->getRefToYieldOffset() = Shear;
+ }
+ #ifdef _DEBUG
+ printf("coalescence ocurrs using yield offset = %e\n",q1Thom->getYieldOffset());
+ #endif //_DEBUG
+ }
+ else
+ {
+ q1Thom->getRefToCrackOffsetOnCft() = 1.;
+ q1Thom->getRefToLodeParameterOffset() = 1.;
+ q1Thom->getRefToYieldOffset() = 0.;
+ #ifdef _DEBUG
+ printf("coalescence ocurrs, not offset is used\n");
+ #endif //_DEBUG
+ }
+ }
+ }
+ else
+ {
+ q1Thom->getRefToCrackOffsetOnCft() = 1.;
+ q1Thom->getRefToLodeParameterOffset() = 1.;
+ q1Thom->getRefToYieldOffset() = 0.;
+ }
+ }
+};
+void mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::forceCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const
+{
+ Msg::Error("mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::forceCoalescence is not implemented!!!");
+};
+bool mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::checkCrackCriterion(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const STensor3& Fcur, const SVector3& normal, double delayFactor, const double* T) const
+{
+ Msg::Error("mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS::checkCrackCriterion is not implemented!!!");
+}
\ No newline at end of file
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.h b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.h
index 16ba5b8ee1a1d3359f529d938b6fd4fb1e38e31a..0cee29a52174713e9847ef78a02c9540e4fbee7b 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.h
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.h
@@ -311,6 +311,9 @@ class mlawNonLocalPorousCoupledHill48LawWithMPSAndMSS : public mlawNonLocalPorou
bool diff=false, STensor43* DNpDkcor=NULL, STensor3* DNpDR=NULL, std::vector<STensor3>* DNpDY=NULL,
bool withThermic=false, STensor3* DNpDT=NULL
) const;
+ virtual void checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const;
+ virtual void forceCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const;
+ virtual bool checkCrackCriterion(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const STensor3& Fcur, const SVector3& normal, double delayFactor, const double* T) const;
};
#endif // MLAWNONLOCALPOROUSCOALESCENCE_H_
diff --git a/NonLinearSolver/materialLaw/mullinsEffect.cpp b/NonLinearSolver/materialLaw/mullinsEffect.cpp
index 492a802249ea7e1f2653eb083b252ac57653f965..50ec2a4da39a0777d63dcd795314f19fe9937c33 100644
--- a/NonLinearSolver/materialLaw/mullinsEffect.cpp
+++ b/NonLinearSolver/materialLaw/mullinsEffect.cpp
@@ -469,9 +469,11 @@ void linearScaler::mullinsEffectScaling(double _psi0, double _psi, const IPMulli
DetaDpsi_cap = r/_psi_max + r*psi_cap/pow(_psi_max,2) * (-Dpsi_max_Dpsi_cap);
DetaDT = -drdT*(1.-psi_cap/_psi_max) + r*(dPsi_capDT*1/_psi_max - psi_cap/pow(_psi_max,2)*Dpsi_max_DT);
+
+ // DetaDT = -drdT*(1.-psi_cap/_psi_max) + r*(- psi_cap/pow(_psi_max,2)*Dpsi_max_Dpsi_cap*dPsi_capDT); // NEW - In paper
double DetaDpsi_max = -r*psi_cap/pow(_psi_max,2);
- DpsiNew_DpsiMax = -r*pow(psi_cap,2)/(2.*pow(_psi_max,2));
+ DpsiNew_DpsiMax = -r*pow(psi_cap,2)/(2.*pow(_psi_max,2)); // definition
DDpsiNew_DpsiMaxDpsi = DetaDpsi_max + r*pow(psi_cap,2)/pow(_psi_max,3)*Dpsi_max_Dpsi_cap;
DDpsiNew_DpsiMaxDT = DDpsiNew_DpsiMaxDpsi*dPsi_capDT - drdT*(pow(psi_cap,2)/(2.*pow(_psi_max,2)));
}
diff --git a/NonLinearSolver/nlsolver/dofManagerMultiSystems.cpp b/NonLinearSolver/nlsolver/dofManagerMultiSystems.cpp
index 54b9cd12fbab09cb0189bf99735d20f72841dbdc..fbe8dcb2075b4929b1a776065f15ba7f716a740c 100644
--- a/NonLinearSolver/nlsolver/dofManagerMultiSystems.cpp
+++ b/NonLinearSolver/nlsolver/dofManagerMultiSystems.cpp
@@ -412,6 +412,15 @@ void dofManagerMultiSystems::clearRHSfixed()
_vmanager[i]->clearRHSfixed();
}
}
+
+double dofManagerMultiSystems::getNorm1ReactionForce() const
+{
+ double val = 0;
+ for(int i=0;i<_vmanager.size();i++){
+ val += _vmanager[i]->getNorm1ReactionForce();
+ }
+ return val;
+}
void dofManagerMultiSystems::getReactionForceOnFixedPhysical(archiveForce& af) const
{
double forceallsys = 0;
diff --git a/NonLinearSolver/nlsolver/dofManagerMultiSystems.h b/NonLinearSolver/nlsolver/dofManagerMultiSystems.h
index 4d4dd3c4fff690b358952741d44ea15899d42e60..de71fbfc26102f4e05731372222f8ae342095655 100644
--- a/NonLinearSolver/nlsolver/dofManagerMultiSystems.h
+++ b/NonLinearSolver/nlsolver/dofManagerMultiSystems.h
@@ -158,6 +158,7 @@ class dofManagerMultiSystems : public nlsDofManager
virtual void preAllocateEntries();
// operations
virtual void clearRHSfixed();
+ virtual double getNorm1ReactionForce() const;
virtual void getReactionForceOnFixedPhysical(archiveForce& af) const;
virtual void getForces(std::vector<archiveForce> &vaf) const;
virtual double getKineticEnergy(const int systemSizeWithoutRigidContact,FilterDofSet &fildofcontact) const;
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index 1e7677b7a94ecbb0284a199b0f1e2f5a2b5681f2..08b3fb3184966c60a6ba96a2a68083ec56af32ef 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -9851,7 +9851,21 @@ double nonLinearMechSolver::computeRightHandSide()
if (_pathFollowing and _pfManager._pathFollowingMethod == pathFollowingManager::LOCAL_BASED){
this->computePathFollowingConstraint();
}
-
+
+ double ffixNorm = pAssembler->getNorm1ReactionForce();
+ int ok = (std::isnan(ffixNorm) or std::isinf(ffixNorm))?1:0;
+ #if defined(HAVE_MPI)
+ if (this->getNumRanks() > 1)
+ {
+ MPI_Allreduce(MPI_IN_PLACE,&ok, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+ }
+ #endif //HAVE_MPI
+ if (ok>0)
+ {
+ Msg::Error("reaction force is nan or inf: rank = %d, val = %e", Msg::GetCommRank(),ffixNorm);
+ return sqrt(-1.);
+ }
+
// save Fint component to archive ( write in file when Fint = Fext)
// return lsys->normInfRightHandSide();
return pAssembler->normInfRightHandSide();
@@ -10094,7 +10108,20 @@ int nonLinearMechSolver::NewtonRaphson(const int numstep){
if(computedUdF()){
_ipf->setDeformationGradientGradient(this->getMicroBC()->getSecondOrderKinematicalVariable(), IPStateBase::current);
}
- _ipf->compute1state(IPStateBase::current, true);
+ bool ok = _ipf->compute1state(IPStateBase::current, true);
+ #if defined(HAVE_MPI)
+ if (this->getNumRanks() > 1)
+ {
+ int val = (ok)? 0: 1;
+ MPI_Allreduce(MPI_IN_PLACE,&val, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+ ok = (val==0);
+ }
+ #endif //HAVE_MPI
+ if (!ok)
+ {
+ Msg::Warning("The evaluation of constutitive laws is failed at some GPs");
+ return _timeManager->getMaxNbIterations(); // for isnan case
+ }
// Compute Right Hand Side (Fext-Fint)
double normFinf = this->computeRightHandSide();
// norm0 = norm(Fext) + norm(Fint) to have a relative convergence
@@ -10130,7 +10157,20 @@ int nonLinearMechSolver::NewtonRaphson(const int numstep){
return _timeManager->getMaxNbIterations();
}
// update ipvariable
- _ipf->compute1state(IPStateBase::current,true);
+ ok = _ipf->compute1state(IPStateBase::current,true);
+ #if defined(HAVE_MPI)
+ if (this->getNumRanks() > 1)
+ {
+ int val = (ok)? 0: 1;
+ MPI_Allreduce(MPI_IN_PLACE,&val, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+ ok = (val==0);
+ }
+ #endif //HAVE_MPI
+ if (!ok)
+ {
+ Msg::Warning("The evaluation of constutitive laws is failed at some GPs");
+ return _timeManager->getMaxNbIterations(); // for isnan case
+ }
if (!_iterativeNR) break;
// new forces
@@ -10160,7 +10200,20 @@ int nonLinearMechSolver::NewtonRaphson(const int numstep){
}
if(ls == 1) break; // processus is converged or failed
// update ipvariable
- _ipf->compute1state(IPStateBase::current,true);
+ ok = _ipf->compute1state(IPStateBase::current,true);
+ #if defined(HAVE_MPI)
+ if (this->getNumRanks() > 1)
+ {
+ int val = (ok)? 0: 1;
+ MPI_Allreduce(MPI_IN_PLACE,&val, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+ ok = (val==0);
+ }
+ #endif //HAVE_MPI
+ if (!ok)
+ {
+ Msg::Warning("The evaluation of constutitive laws is failed at some GPs");
+ return _timeManager->getMaxNbIterations(); // for isnan case
+ }
// new RightHandSide
this->computeRightHandSide();
}
@@ -10565,7 +10618,20 @@ int nonLinearMechSolver::NewtonRaphsonPathFollowing(const int numstep)
break;
}
// update ipvariable
- _ipf->compute1state(IPStateBase::current,true);
+ bool ok = _ipf->compute1state(IPStateBase::current,true);
+ #if defined(HAVE_MPI)
+ if (this->getNumRanks() > 1)
+ {
+ int val = (ok)? 0: 1;
+ MPI_Allreduce(MPI_IN_PLACE,&val, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+ ok = (val==0);
+ }
+ #endif //HAVE_MPI
+ if (!ok)
+ {
+ Msg::Warning("The evaluation of constutitive laws is failed at some GPs");
+ return _timeManager->getMaxNbIterations(); // for isnan case
+ }
// update load
double loadParam = pathsys->getControlParameter();
double dloadParam = fabs(pathsys->getControlParameterStep());
diff --git a/NonLinearSolver/nlsolver/staticDofManager.cpp b/NonLinearSolver/nlsolver/staticDofManager.cpp
index db51d8a0edafa302d6261b4c43d9108b37c1bda8..e021c3d83e28768a19d41619a4171c3e8e74b7b7 100644
--- a/NonLinearSolver/nlsolver/staticDofManager.cpp
+++ b/NonLinearSolver/nlsolver/staticDofManager.cpp
@@ -1405,6 +1405,16 @@ void staticDofManager::clearRHSfixed()
itR->second = 0.;
};
+double staticDofManager::getNorm1ReactionForce() const
+{
+ double maxVal = 0;
+ for(std::map<Dof,double>::const_iterator itR = RHSfixed.begin();itR != RHSfixed.end(); ++itR)
+ {
+ maxVal += fabs(itR->second);
+ }
+ return maxVal;
+};
+
void staticDofManager::getReactionForceOnFixedPhysical(archiveForce& af) const
{
af.fval =0.;
diff --git a/NonLinearSolver/nlsolver/staticDofManager.h b/NonLinearSolver/nlsolver/staticDofManager.h
index 9089df1148c76a5f86939de401a8525671ff358f..e71fea8286f16c789229b1d00bf421f3480833dd 100644
--- a/NonLinearSolver/nlsolver/staticDofManager.h
+++ b/NonLinearSolver/nlsolver/staticDofManager.h
@@ -124,6 +124,7 @@ class nlsDofManager
virtual void preAllocateEntries() = 0;
// operations
virtual void clearRHSfixed() = 0;
+ virtual double getNorm1ReactionForce() const = 0;
virtual void getReactionForceOnFixedPhysical(archiveForce& af) const = 0;
virtual void getForces(std::vector<archiveForce> &vaf) const = 0;
virtual double getKineticEnergy(const int systemSizeWithoutRigidContact,FilterDofSet &fildofcontact) const = 0;
@@ -265,7 +266,7 @@ class staticDofManager : public nlsDofManager
// zero for RHSfixed
virtual void clearRHSfixed();
-
+ virtual double getNorm1ReactionForce() const;
virtual void getReactionForceOnFixedPhysical(archiveForce& af) const;
virtual void getForces(std::vector<archiveForce> &vaf) const;
virtual double getKineticEnergy(const int systemSizeWithoutRigidContact,FilterDofSet &fildofcontact) const;
diff --git a/README.md b/README.md
index 0a546b2d5824accf6726f7d456d7efc042047255..a1bf7cc7a3ad4a8b3619d938abb1abe8a496820b 100644
--- a/README.md
+++ b/README.md
@@ -176,8 +176,9 @@ To use torch you can follow the instructions either [2.1. Without GPU capability
```bash
cd && cd local
```
+Here below to update with the last stable version from ```https://pytorch.org/```
```bash
-wget https://download.pytorch.org/libtorch/cpu/libtorch-shared-with-deps-2.1.1%2Bcpu.zip
+wget https://download.pytorch.org/libtorch/cpu/libtorch-shared-with-deps-2.3.0%2Bcpu.zip
```
2. Unzip the archive, this will be extracted in ```$HOME/local/libtorch``` by default so make sure that directory does not exist:
```bash
diff --git a/dG3D/benchmarks/CMakeLists.txt b/dG3D/benchmarks/CMakeLists.txt
index c1aa55d376a128f263e4431660c207cdde3bcd57..586558b5567154cd8f0eca3f1e50d639ccdd954a 100644
--- a/dG3D/benchmarks/CMakeLists.txt
+++ b/dG3D/benchmarks/CMakeLists.txt
@@ -224,6 +224,7 @@ add_subdirectory(defoDefoContactSlaveEdge2DQuadPert)
add_subdirectory(defoDefoContactSlaveEdge2D2ndOrderTri)
add_subdirectory(defoDefoContactSlaveVolume3DLinear)
add_subdirectory(defoDefoContactSlaveFace3DLinear)
+add_subdirectory(impactSphere)
add_subdirectory(Gurson_Cube_I1J2J3)
add_subdirectory(GursonMultipleNonlocal_PathFollowing)
add_subdirectory(Thomason_cube_I1J2J3)
diff --git a/dG3D/benchmarks/impactSphere/CMakeLists.txt b/dG3D/benchmarks/impactSphere/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b4c79a585554fd46d4d54148676cfefad1997936
--- /dev/null
+++ b/dG3D/benchmarks/impactSphere/CMakeLists.txt
@@ -0,0 +1,10 @@
+# test file
+
+set(PYFILE impact.py)
+
+set(FILES2DELETE
+ *.msh
+ *.csv
+)
+
+add_cm3python_test(${PYFILE} "${FILES2DELETE}")
diff --git a/dG3D/benchmarks/impactSphere/VEVP.py b/dG3D/benchmarks/impactSphere/VEVP.py
new file mode 100755
index 0000000000000000000000000000000000000000..cb6af3a13bec6a3b123190e68abf8c2e874b8b4e
--- /dev/null
+++ b/dG3D/benchmarks/impactSphere/VEVP.py
@@ -0,0 +1,200 @@
+#coding-Utf-8-*-
+from gmshpy import *
+from dG3Dpy import *
+
+import pickle
+import numpy as np
+
+
+def FillMatPara(x_E,x_P,VEVP): #x = [E0,nu0,K1,G1,....,tk,tg, aP]
+ tk = np.power(10.0, x_E[int(len(x_E)-3)])
+ tg = np.power(10.0, x_E[int(len(x_E)-2)])
+ MatPara = []
+ MatPara.append([x_E[0], x_E[1]])
+ Ki = []
+ Gi = []
+ N = int((len(x_E)-5)/2)
+ for i in range(N):
+ a = 1.0/np.power(10.0, i)
+ Ki.append([np.power(10.0,x_E[2+2*i]), a*tk])
+ Gi.append([np.power(10.0,x_E[2+2*i+1]), a*tg])
+ MatPara.append(Ki)
+ MatPara.append(Gi)
+
+
+ if(VEVP):
+ Compression = []
+ Tension = []
+ ViscoPlastic = []
+
+ Compression.append(x_P[0])
+ Tension.append(x_P[4])
+ for i in range(1,3):
+ Compression.append(np.power(10.0, x_P[i]))
+ Tension.append(np.power(10.0, x_P[4+i]))
+ Compression.append(x_P[3])
+ Tension.append(x_P[7])
+
+ for i in range(4):
+ ViscoPlastic.append(x_P[8+i])
+
+ MatPara.append(Compression)
+ MatPara.append(Tension)
+ MatPara.append(ViscoPlastic)
+ return MatPara
+############################################################################
+
+
+
+
+def RunCase(LoadCase,MatPara, VEVP):
+
+ nstep = 800 # number of step (used only if soltype=1)
+ # === Problem parameters ============================================================================
+ E0 = MatPara[0][0]
+ nu0 = MatPara[0][1]
+ Nt = len(MatPara[1]) #the number of branches in the generalized Maxwell model
+ aP = 1.0
+
+ if(VEVP):
+ indxplastic = 3
+ sy0c = MatPara[indxplastic][0] #MPa, compressive yield stress
+ hc = MatPara[indxplastic][1]
+ hc2 = MatPara[indxplastic][2]
+ kc = MatPara[indxplastic][3]
+
+ sy0t = MatPara[indxplastic+1][0]*sy0c #MPa, compressive yield stress
+ ht = MatPara[indxplastic+1][1]
+ ht2 = MatPara[indxplastic+1][2]
+ kt = MatPara[indxplastic+1][3]
+
+ alpha = MatPara[indxplastic+2][0]
+ beta = MatPara[indxplastic+2][1]
+ eta = MatPara[indxplastic+2][2]
+ p = MatPara[indxplastic+2][3]
+
+ #=============================================
+ lawnum = 11
+ rho = 7850e-9
+
+ if(VEVP):
+ law1 = HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(lawnum,rho,E0,nu0)
+ else:
+ law1 = HyperViscoElasticDG3DMaterialLaw(lawnum,rho,E0,nu0)
+
+ law1.setViscoelasticMethod(0)
+ law1.setViscoElasticNumberOfElement(Nt)
+ for i in range(Nt):
+ law1.setViscoElasticData_Bulk(i+1, MatPara[1][i][0], aP*MatPara[1][i][1])
+ law1.setViscoElasticData_Shear(i+1, MatPara[2][i][0], aP*MatPara[2][i][1])
+ law1.setStrainOrder(5)
+ #-----------------------------------------------------------------------------
+ if(VEVP):
+ hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, hc, hc2, kc)
+ hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, ht2, kt)
+
+ law1.setCompressionHardening(hardenc)
+ law1.setTractionHardening(hardent)
+
+ law1.setYieldPowerFactor(alpha)
+ law1.setNonAssociatedFlow(True)
+ law1.nonAssociatedFlowRuleFactor(beta)
+ etac = constantViscosityLaw(1,eta)
+ law1.setViscosityEffect(etac,p)
+
+
+ if(LoadCase[0] == 0):
+ ft1 = 5.732
+ eps = -0.0157
+ if(VEVP):
+ ft1 = 106.5
+ eps = -0.296
+ ftime = LoadCase[2]
+ else:
+ strainrate = LoadCase[1]
+ strain_end = LoadCase[2]
+ ftime = strain_end/strainrate # Final time (used only if soltype=1)
+ # print("######################",strain_end,strainrate,ftime)
+ # ===Solver parameters=============================================================================
+ sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
+ soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
+ tol=1.e-5 # relative tolerance for NR scheme (used only if soltype=1)
+ nstepArch=10 # Number of step between 2 archiving (used only if soltype=1)
+ fullDg = 0 #O = CG, 1 = DG
+ space1 = 0 # function space (Lagrange=0)
+ beta1 = 100
+
+ # ===Domain creation===============================================================================
+ # creation of ElasticField
+ myfield1 = dG3D1DDomain(11,11,space1,lawnum, False) #for non local
+ myfield1.setState(1) #UniaxialStrain=0, UniaxialStress =1,ConstantTriaxiality=2
+ L = 5.
+
+ # ===Solver creation===============================================================================
+ # geometry & mesh
+ geofile = "line.geo"
+ meshfile = "line.msh"
+
+ # creation of Solver
+ mysolver = nonLinearMechSolver(1000)
+ #mysolver.createModel(geofile,meshfile,1,1)
+ mysolver.loadModel(meshfile)
+ mysolver.addDomain(myfield1)
+ mysolver.addMaterialLaw(law1)
+ mysolver.Scheme(soltype)
+ mysolver.Solver(sol)
+ mysolver.snlData(nstep,ftime,tol)
+ mysolver.stepBetweenArchiving(nstepArch)
+
+ mysolver.snlManageTimeStep(15,5,2.,20) # maximal 20 times
+
+# ===Boundary conditions===========================================================================
+ mysolver.displacementBC("Edge",11,1,0.)
+ mysolver.displacementBC("Edge",11,2,0.)
+ mysolver.displacementBC("Node",1,0,0.)
+
+ disp = PiecewiseLinearFunction()
+ disp.put(0.,0.)
+ if(LoadCase[0] != 0):
+ disp.put(ftime, LoadCase[0]*strain_end*L)
+ else:
+ disp.put(ft1,eps*L)
+ disp.put(ftime,eps*L)
+
+ mysolver.displacementBC("Node",2,0,disp)
+ mysolver.archivingAverageValue(IPField.F_XX)
+ mysolver.archivingAverageValue(IPField.P_XX)
+
+ # ===Solving=======================================================================================
+ mysolver.SetVerbosity(0)
+ mysolver.solve()
+ return
+
+
+#def test():
+# print("It works!!")
+
+
+#MatPara=
+#RunCase(L_Case,MatPara,VEVP)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/dG3D/benchmarks/impactSphere/impact.py b/dG3D/benchmarks/impactSphere/impact.py
new file mode 100755
index 0000000000000000000000000000000000000000..c31e3a8c47ac4e05214232cbf6be3f7ce37edf9d
--- /dev/null
+++ b/dG3D/benchmarks/impactSphere/impact.py
@@ -0,0 +1,245 @@
+#coding-Utf-8-*-
+from gmshpy import *
+from dG3Dpy import*
+# import functions from VEVP.py
+from VEVP import *
+import math
+
+
+# geometry ===================================================================
+# Input .msh
+#meshfile = "square.msh" # name of mesh file
+
+meshfile = 'impactSphere.msh'
+geofile = 'impactSphere.geo' # name of mesh file
+order = 1 #order of elements, better use one with contact
+rigidContact = False #True to use a contact below, False to clamp the plate
+
+
+#parameters
+nstep = 400 # number of step
+ftime = 1e-3 # Final time v0 = 0.1mm/s, pour la demi structure v = v0/2
+nstepArch = 4 # Number of step between 2 archiving (used only if soltype=1)
+tol = 1e-4
+
+penalty= 1.e3
+
+# material law for impactor
+EImpactor = 2e5
+nuImpactor = 0.3
+mm=1.
+RextImpactor = 60.*mm
+RintImpactor = 50.*mm
+lawnum3 = 3
+VolImpactor = math.pi*(RextImpactor**3-RintImpactor**3)*4./3./8.
+MassImpactor = 5.*(1e-3/mm)/4. #in tonnes if mm=1
+rhoImpactor = MassImpactor/VolImpactor
+velocity = -1.*(2.*9.81*1.8)**0.5*(mm/1.e-3)
+print("Density Impactor = ", rhoImpactor, "; Velocity Impactor = ", velocity)
+
+law3 = J2LinearDG3DMaterialLaw(lawnum3,rhoImpactor,EImpactor,nuImpactor,1000.*EImpactor,10.)
+
+#VEVP Law for plates as to generate the sve: line 70 000 of the MCMC
+x=[1374.2643289981168, 0.2826036623239532, 4.220196044260534, -0.06622730551411246, -0.2361729746734328, 1.0658488572731906, 2.955587774357835, 0.8846323146785388, 2.1134546059887156, 0.3900795755310005, -0.07807917229891632, 1.3340258222804844, 3.8978259552106973, 0.8147820955100668, 1.7082119337043935, 3.1818481562348735, 3.1460140132039993, 0.952562223405303, 1.0351529596236166, 0.8971311953490693, 0.04770622785629489, 29.696376534920287, -2.5325714558968597, 1.1573320881902116, 299.0940006076202, 0.7771046411287497, 2.1446624378399077, -3.183873130566032, 9648.947729066591, 3.7047493821391746, 0.2175363944799915, 5.182554847703984, 0.17938765700173523]
+
+
+numVE = 21
+number_of_parametrs = len(x)
+x_E = list(x)[0:numVE]
+x_P = list(x)[numVE:number_of_parametrs]
+VEVP = True # "Compression"
+MatPara=FillMatPara(x_E,x_P,VEVP)
+
+E0 = MatPara[0][0]
+nu0 = MatPara[0][1]
+Nt = len(MatPara[1]) #the number of branches in the generalized Maxwell model
+aP = 1.0
+
+indxplastic = 3
+sy0c = MatPara[indxplastic][0] #MPa, compressive yield stress
+hc = MatPara[indxplastic][1]
+hc2 = MatPara[indxplastic][2]
+kc = MatPara[indxplastic][3]
+
+sy0t = MatPara[indxplastic+1][0]*sy0c #MPa, compressive yield stress
+ht = MatPara[indxplastic+1][1]
+ht2 = MatPara[indxplastic+1][2]
+kt = MatPara[indxplastic+1][3]
+
+alpha = MatPara[indxplastic+2][0]
+beta = 1.5*MatPara[indxplastic+2][1]
+eta = MatPara[indxplastic+2][2]
+p = MatPara[indxplastic+2][3]
+
+lawnum1 = 11
+rho = 1200e-12
+
+law1 = HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(lawnum1,rho,E0,nu0)
+law1.setViscoelasticMethod(0)
+law1.setViscoElasticNumberOfElement(Nt)
+for i in range(Nt):
+ law1.setViscoElasticData_Bulk(i+1, MatPara[1][i][0], aP*MatPara[1][i][1])
+ law1.setViscoElasticData_Shear(i+1, MatPara[2][i][0], aP*MatPara[2][i][1])
+law1.setStrainOrder(5)
+
+hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, hc, hc2, kc)
+hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, ht2, kt)
+
+law1.setCompressionHardening(hardenc)
+law1.setTractionHardening(hardent)
+
+law1.setYieldPowerFactor(alpha)
+law1.setNonAssociatedFlow(True)
+law1.nonAssociatedFlowRuleFactor(beta)
+etac = constantViscosityLaw(1,eta)
+law1.setViscosityEffect(etac,p)
+if (order==1):
+ law1.setUseBarF(True)
+
+# lattice here to change
+lawnum2=2
+law2 = HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(lawnum2,rho*0.02,E0*0.02,nu0)
+law2.setViscoelasticMethod(0)
+law2.setViscoElasticNumberOfElement(Nt)
+for i in range(Nt):
+ law2.setViscoElasticData_Bulk(i+1, MatPara[1][i][0]*0.02, aP*MatPara[1][i][1])
+ law2.setViscoElasticData_Shear(i+1, MatPara[2][i][0]*0.02, aP*MatPara[2][i][1])
+law2.setStrainOrder(5)
+
+hardenc2 = LinearExponentialJ2IsotropicHardening(3, sy0c*0.02, hc*0.02, hc2*0.02, kc)
+hardent2 = LinearExponentialJ2IsotropicHardening(4, sy0t*0.02, ht*0.02, ht2*0.02, kt)
+
+law2.setCompressionHardening(hardenc2)
+law2.setTractionHardening(hardent2)
+
+law2.setYieldPowerFactor(alpha)
+law2.setNonAssociatedFlow(True)
+law2.nonAssociatedFlowRuleFactor(beta)
+eta2 = constantViscosityLaw(2,eta)
+law2.setViscosityEffect(eta2,p)
+
+if (order==1):
+ law2.setUseBarF(True) #maybe no need if no locking
+
+
+# data law def ============================================================
+
+# solver info ============================================================
+sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
+soltype = 4 # StaticLinear=0 (default) StaticNonLinear=1
+fullDg = 0 #O = CG, 1 = DG
+space1 = 0 # function space (Lagrange=0)
+# solver info ============================================================
+
+# creation of Domain
+#Impactor
+nfield3 = 1001 # number of the field (physical number of volume in 3D)
+myfield3 = dG3DDomain(1000,nfield3,space1,lawnum3,fullDg)
+
+#Plates
+nfield1_1 = 51 # number of the field (physical number of volume in 3D)
+myfield1_1 = dG3DDomain(1001,nfield1_1,space1,lawnum1,fullDg)
+myfield1_1.strainSubstep(2,10)
+
+nfield1_2 = 52 # number of the field (physical number of volume in 3D)
+myfield1_2 = dG3DDomain(1001,nfield1_2,space1,lawnum1,fullDg)
+myfield1_2.strainSubstep(2,10)
+
+#lattice
+nfield2 = 53 # number of the field (physical number of volume in 3D)
+myfield2 = dG3DDomain(1001,nfield2,space1,lawnum2,fullDg)
+myfield2.strainSubstep(2,10)
+
+
+
+
+# creation of ElasticField ==============================================================
+
+# Solver ==============================================================
+mysolver = nonLinearMechSolver(1002)
+mysolver.createModel(geofile,meshfile,3,order)
+#mysolver.loadModel(meshfile)
+mysolver.addDomain(myfield1_1)
+mysolver.addDomain(myfield1_2)
+mysolver.addDomain(myfield2)
+mysolver.addDomain(myfield3)
+mysolver.addMaterialLaw(law1)
+mysolver.addMaterialLaw(law2)
+mysolver.addMaterialLaw(law3)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol,tol/100.)
+mysolver.stepBetweenArchiving(nstepArch)
+mysolver.implicitSpectralRadius(0.05)
+mysolver.snlManageTimeStep(50, 3, 2, 10)
+#mysolver.options("-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps")
+
+# Solver ==============================================================
+
+# BC ==============================================================
+# blocked in z bottom if no rigid contact with plane
+if rigidContact==False:
+ mysolver.displacementBC("Face",200,2,0.)
+
+# blocked in y side
+mysolver.displacementBC("Face",101,1,0.)
+mysolver.displacementBC("Face",1001,1,0.)
+# blocked in x side
+mysolver.displacementBC("Face",100,0,0.)
+mysolver.displacementBC("Face",1000,0,0.)
+
+mysolver.initialBC("Volume","Velocity",1001,2,velocity)
+
+# BC ==============================================================
+
+# Contact
+
+
+defoDefoContact1=dG3DNodeOnSurfaceContactDomain(1002, 2, 2000, 2,201, 10.*penalty, 3.*mm)
+mysolver.defoDefoContactInteraction(defoDefoContact1)
+defoDefoContact2=dG3DNodeOnSurfaceContactDomain(1003, 2, 201, 2,2000, penalty, 3.*mm)
+mysolver.defoDefoContactInteraction(defoDefoContact2)
+
+
+
+if rigidContact==True:
+ flaw1 = CoulombFrictionLaw(1,0.,0.0,penalty,penalty/1000.)
+ contact1 = dG3DRigidPlaneContactDomain(1004,2,10001,3,nfield1_1,11001,0.,0.,1.,penalty,3.*mm,1e3)
+ contact1.setFriction(False)
+ contact1.addFrictionLaw(flaw1)
+ mysolver.contactInteraction(contact1)
+
+ #mysolver.displacementRigidContactBC(10001,2,0.)
+
+ #mysolver.archivingForceOnPhysicalGroup("Face", 10001, 2,1)
+ #mysolver.archivingRigidContact(10001,2,0, 1)
+ #mysolver.archivingRigidContactForce(11001, 2, 1)
+#Contact
+
+
+# save
+mysolver.internalPointBuildView("svm",IPField.SVM, 1, nstepArch)
+# save stress tensor
+mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, nstepArch)
+mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, nstepArch)
+mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, nstepArch)
+mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, nstepArch)
+mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, nstepArch)
+mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, nstepArch)
+# save platic strain
+mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, nstepArch)
+# save bottom force
+mysolver.archivingForceOnPhysicalGroup("Face", 100, 2,1)
+
+#save impactor displacement, velocity etc ...
+mysolver.archivingNodeDisplacement(1002, 2,1)
+mysolver.archivingNodeVelocity(1002, 2,1);
+mysolver.archivingNodeAcceleration(1002, 2, 1);
+
+
+mysolver.solve()
+
+check = TestCheck()
+check.equal(4.059531e+02,mysolver.getArchivedForceOnPhysicalGroup("Face", 100, 2),1.e-6)
+
+
diff --git a/dG3D/benchmarks/impactSphere/impactSphere.geo b/dG3D/benchmarks/impactSphere/impactSphere.geo
new file mode 100644
index 0000000000000000000000000000000000000000..1f35bd3568a18e44717a47bc6df53acd61ae086e
--- /dev/null
+++ b/dG3D/benchmarks/impactSphere/impactSphere.geo
@@ -0,0 +1,133 @@
+SetFactory("OpenCASCADE");
+
+Mesh.Optimize = 1;
+Mesh.SecondOrderLinear = 1;
+mm = 1.; // Unit
+//Geometry.AutoCoherence = 0.000000001*mm;
+//Geometry.Tolerance = 0.000000001*mm; // adjust value here for correct merge result
+Mesh.CharacteristicLengthMin = 1.*mm;
+Mesh.CharacteristicLengthMax = 25.*mm;
+
+
+
+
+
+
+RextLattice = 100*mm;
+HLattice = 4.*7.5*mm;
+LLattice = 7.5*mm; //mesh size
+
+
+LowerPlate = 2.*mm;
+UpperPlate = 2.*mm;
+LPlate = 7.5*mm; //mesh size
+
+RextImpactor = 60*mm;
+RintImpactor = 50.*mm;
+LImpactor = 10*mm; //mesh size
+
+
+// lower plate
+plate()+=newv; Cylinder(newv) = {0,0,0, 0,0,LowerPlate, RextLattice, Pi/2.};
+// upper plate
+plate()+=newv; Cylinder(newv) = {0,0,LowerPlate+HLattice, 0,0,UpperPlate, RextLattice, Pi/2.};
+
+Characteristic Length { PointsOf{ Volume{ plate[] }; } } = LPlate;
+
+// lattice part
+lattice()+=newv; Cylinder(newv) = {0,0,LowerPlate, 0,0,HLattice, RextLattice, Pi/2.};
+Characteristic Length { PointsOf{ Volume{ lattice[] }; } } = LLattice;
+
+
+
+
+//impactor
+//+
+impactortmp()+=newv; Sphere(newv) = {0.-0.0001*mm, 0.-0.0001*mm, LowerPlate+HLattice+UpperPlate+RextImpactor+0.000001*mm, RextImpactor, -Pi/2, 0., Pi/2.};
+//+
+impactortmp()+=newv; Sphere(newv) = {0.-0.0001*mm, 0.-0.0001*mm, LowerPlate+HLattice+UpperPlate+RextImpactor+0.000001*mm, RintImpactor, -Pi/2, 0., Pi/2.};
+
+impactor()+=newv; BooleanDifference(newv) = { Volume{impactortmp(0)}; Delete; }{ Volume{impactortmp(1)}; Delete; };
+Characteristic Length { PointsOf{ Volume{ impactor[] }; } } = LImpactor;
+
+
+Coherence;
+
+surflowerplate[] = Abs(Boundary{ Volume{plate(0)}; });
+Printf("surf lower plate ",surflowerplate());
+surfupperplate[] = Abs(Boundary{ Volume{plate(1)}; });
+Printf("surf upper plate ",surfupperplate());
+surflattice[] = Abs(Boundary{ Volume{lattice(0)}; });
+Printf("surf lattice ",surflattice());
+surfimpactor[] = Abs(Boundary{ Volume{impactor(0)}; });
+
+Printf("surf impactor ",surfimpactor());
+
+ptimpactor[] = PointsOf{ Surface{surfimpactor(1)}; };
+Printf("pts impactor ",ptimpactor());
+
+
+
+
+
+OXZ=101;
+Physical Surface(OXZ) = {surflowerplate(3), surflattice(3),surfupperplate(3)};
+OYZ=100;
+Physical Surface(OYZ) = {surflowerplate(4), surflattice(4),surfupperplate(4)};
+
+OXYLOW=200;
+Physical Surface(OXYLOW) = {surflowerplate(2)};
+OXYUP=201;
+Physical Surface(OXYUP) = {surfupperplate(1)};
+
+LOWERPLATE=51;
+Physical Volume(LOWERPLATE) = {plate(0)};
+UPPERPLATE=52;
+Physical Volume(UPPERPLATE) = {plate(1)};
+
+LATTICE=53;
+Physical Volume(LATTICE) = {lattice(0)};
+
+
+OXZIMPACTOR=1001;
+Physical Surface(OXZIMPACTOR) = {surfimpactor(3)};
+OYZIMPACTOR=1000;
+Physical Surface(OYZIMPACTOR) = {surfimpactor(2)};
+
+OXYIMPACTORLOW=2000;
+Physical Surface(OXYIMPACTORLOW) ={surfimpactor(0)};
+OXYIMPACTORUP=2001;
+Physical Surface(OXYIMPACTORUP) = {surfimpactor(1)};
+
+IMPACTOR=1001;
+Physical Volume(IMPACTOR) = {impactor(0)};
+
+PTIMPACTOR = 1002;
+Physical Point(PTIMPACTOR) = {ptimpactor(2)};
+
+
+//rigid contact
+X_contact=-0.1*mm;
+Y_contact=-0.1*mm;
+Z_contact=-0.0000001*mm;
+L_contact=RextLattice+2.*mm;
+Point(11001) = {X_contact,Y_contact,Z_contact};
+Point(11002) = {X_contact,Y_contact+L_contact,Z_contact};
+Point(11003) = {X_contact+L_contact,Y_contact+L_contact,Z_contact};
+Point(11004) = {X_contact+L_contact,Y_contact,Z_contact};
+Line(10001) = {11001, 11002};
+Line(10002) = {11002, 11003};
+Line(10003) = {11003, 11004};
+Line(10004) = {11004, 11001};
+Line Loop(10001) = {10004, 10001, 10002, 10003};
+Physical Point(11001) = {11001};
+Plane Surface(10001) = {10001};
+Transfinite Line {10001,10002,10003,10004} = 2 Using Progression 1;
+Transfinite Surface {10001};
+Recombine Surface {10001};
+RIGIDCONTACT=10001;
+Physical Surface(RIGIDCONTACT) = {10001};
+
+Coherence;
+
+
diff --git a/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cube.geo b/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cube.geo
index e322959204440d14a93a5e085dce5e6ab743d194..7b69582119311ad0f199829d872fc7eb4f36c967 100644
--- a/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cube.geo
+++ b/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cube.geo
@@ -1,5 +1,5 @@
// Cube.geo
-mm=1.0e-3; // Units
+mm=1.0; // Units
n=1; // Number of layers (in thickness / along z)
m = 1; // Number of element + 1 along y
p = 1; // Number of element + 1 along x
diff --git a/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cubeTVP.py b/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cubeTVP.py
index b9de9ab09a709683bfe6b7b5b5194dfe25eb2843..f6b4f958964f031502eb52dc8f688ea8392b7a68 100644
--- a/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cubeTVP.py
+++ b/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cubeTVP.py
@@ -91,7 +91,7 @@ hardenk.setCoefficients(1,1.) # 300.)
hardenk.setCoefficients(2,2.) # 200.)
hardenk.setCoefficients(3,5.)
-law1 = NonLinearTVENonLinearTVPDG3DMaterialLaw(lawnum,rho,E,nu,1e-6,Tinitial,Alpha,KThCon,Cp,False,1e-8)
+law1 = NonLinearTVENonLinearTVPDG3DMaterialLaw(lawnum,rho,E,nu,1e-6,Tinitial,Alpha,KThCon,Cp,False,1e-8,1e-8)
law1.setCompressionHardening(hardenc)
law1.setTractionHardening(hardent)
@@ -119,22 +119,29 @@ law1.setViscoElasticNumberOfElement(N)
if N > 0:
for i in range(1, N+1):
law1.setViscoElasticData(i, relSpec[i], relSpec[i+N])
- law1.setCorrectionsAllBranchesTVE(i, 2000., 1.686, 200., 1.686) # at 23°C -> TPU
+ law1.setCorrectionsAllBranchesTVE(i, 1., 2000., 1.686, 1., 200., 1.686) # at 23°C -> TPU
law1.setAdditionalCorrectionsAllBranchesTVE(i,0.22,0.,0.,0.22,0.,0.)
law1.setCompressionCorrectionsAllBranchesTVE(i, 1.) # 1.5)
law1.setTemperatureFunction_ThermalDilationCoefficient(Alpha_tempfunc)
law1.setTemperatureFunction_ThermalConductivity(KThCon_tempfunc)
+mullins = linearScaler(4, 0.55)
+a = -0.9/(90.)
+x0 = 273.15-20
+linear_tempfunc = linearScalarFunction(x0,0.9,a)
+# mullins.setTemperatureFunction_r(linear_tempfunc)
+# law1.setMullinsEffect(mullins)
+
eta = constantViscosityLaw(1,1000.) #(1,3e4)
law1.setViscosityEffect(eta,0.1)
# eta.setTemperatureFunction(negExp_tempfunc)
-law1.setIsotropicHardeningCoefficients(1.,1.,1.)
+# law1.setIsotropicHardeningCoefficients(1.,1.,1.)
# solver
sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
-nstep = 1000 # 100 # number of step (used only if soltype=1)
+nstep = 100 # 100 # number of step (used only if soltype=1)
ftime = 0.65/2.78e-3 # 4. # Final time (used only if soltype=1)
tol=1.e-5 # relative tolerance for NR scheme (used only if soltype=1)
nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
@@ -157,10 +164,11 @@ myfield1 = ThermoMechanicsDG3DDomain(1000,nfield,space1,lawnum,fullDg,ThermoMech
myfield1.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
myfield1.stabilityParameters(beta1)
myfield1.ThermoMechanicsStabilityParameters(beta1,bool(1))
+# myfield1.strainSubstep(2, 10)
# creation of Solver
mysolver = nonLinearMechSolver(1000)
-mysolver.createModel(geofile,meshfile,3,1)
+mysolver.createModel(geofile,meshfile,3,2)
mysolver.loadModel(meshfile)
mysolver.addDomain(myfield1)
mysolver.addMaterialLaw(law1)
@@ -168,9 +176,10 @@ mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstep,ftime,tol)
mysolver.stepBetweenArchiving(nstepArch)
+mysolver.snlManageTimeStep(15,5,2.,20) # maximal 20 times
# BCs
-length = 1.e-3
+length = 1.
umax = -0.65*length # 0.04*length
fu_L = LinearFunctionTime(0,0,ftime,umax); # Linear Displacement with time
@@ -218,4 +227,4 @@ mysolver.solve()
# Check
check = TestCheck()
# check.equal(-2.917493481719865e-05,mysolver.getArchivedForceOnPhysicalGroup("Face", 30, 0),1.e-4)
-check.equal(-2.911073148114768e-05,mysolver.getArchivedForceOnPhysicalGroup("Face", 30, 0),1.e-4)
+check.equal(2.911073148114768e+01,mysolver.getArchivedForceOnPhysicalGroup("Face", 30, 0),1.e-4)
diff --git a/dG3D/src/dG3D1DDomain.cpp b/dG3D/src/dG3D1DDomain.cpp
index 2bdc1ddc9e676efeee714ede54184692952ec854..386f8bbd0c96426aa0bd8855650307727cd06aeb 100644
--- a/dG3D/src/dG3D1DDomain.cpp
+++ b/dG3D/src/dG3D1DDomain.cpp
@@ -201,10 +201,10 @@ void dG3D1DDomain::createTerms(unknownField *uf,IPField*ip)
}
-void dG3D1DDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
+bool dG3D1DDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
materialLaw *mlaw__,fullVector<double> &disp, bool stiff)
{
- if (!getElementErosionFilter()(e)) return;
+ if (!getElementErosionFilter()(e)) return true;
dG3DMaterialLaw *mlaw;
if(mlaw__->getType() == materialLaw::fracture)
@@ -237,5 +237,10 @@ void dG3D1DDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichSt
else {
Msg::Error("state must be set via dG3D1DDomain::setState");
}
+ if (!(ipv->getConstitutiveSuccessFlag()))
+ {
+ return false;
+ }
}
+ return true;
};
diff --git a/dG3D/src/dG3D1DDomain.h b/dG3D/src/dG3D1DDomain.h
index ab5836607e4ac456988e20922de7043bb2715518..516a6a8595a59ec540fa0208aa2fa0c54d88b60d 100644
--- a/dG3D/src/dG3D1DDomain.h
+++ b/dG3D/src/dG3D1DDomain.h
@@ -46,13 +46,14 @@ class dG3D1DDomain : public dG3DDomain{
virtual void createTerms(unknownField *uf,IPField*ip);
virtual partDomain* clone() const {return new dG3D1DDomain(*this);};
- virtual void computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
+ virtual bool computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
materialLaw *mlaw,fullVector<double> &disp, bool stiff);
- virtual void computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
+ virtual bool computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
partDomain* efMinus, partDomain *efPlus,materialLaw *mlawminus,
materialLaw *mlawplus,fullVector<double> &dispm,
fullVector<double> &dispp, const bool virt, bool stiff,const bool checkfrac=true){
Msg::Error("dG3D1DDomain::computeIpv has not been implemented");
+ return false;
};
#endif //SWIG
diff --git a/dG3D/src/dG3DDomain.cpp b/dG3D/src/dG3DDomain.cpp
index 4e9a31c6422d043b7fd4519c602e83b2f5a6f690..2d71a3ce558b98081ce9e14036bedb31214c00f7 100644
--- a/dG3D/src/dG3DDomain.cpp
+++ b/dG3D/src/dG3DDomain.cpp
@@ -1923,7 +1923,7 @@ void dG3DDomain::setdFmdFM(AllIPState *aips, std::map<Dof, STensor3> &dUdF,const
else
mlawbulkp= static_cast<dG3DMaterialLaw*>(mlawPlus);
if (mlawbulkm->getUseBarF()!=mlawbulkp->getUseBarF()) //check with the bulk law
- Msg::Error("dG3DDomain::computeIpv All the constitutive laws need to use Fbar or not");
+ Msg::Error("dG3DDomain::setdFmdFM All the constitutive laws need to use Fbar or not");
bool useBarF=mlawbulkm->getUseBarF();
for(elementGroup::elementContainer::const_iterator it=gi->begin(); it!=gi->end();++it)
@@ -1981,7 +1981,7 @@ void dG3DDomain::setdFmdGM(AllIPState *aips, std::map<Dof, STensor33> &dUdG,cons
else
mlawbulkp= static_cast<dG3DMaterialLaw*>(mlawPlus);
if (mlawbulkm->getUseBarF()!=mlawbulkp->getUseBarF()) //check with the bulk law
- Msg::Error("dG3DDomain::computeIpv All the constitutive laws need to use Fbar or not");
+ Msg::Error("dG3DDomain::setdFmdGM All the constitutive laws need to use Fbar or not");
bool useBarF=mlawbulkm->getUseBarF();
for(elementGroup::elementContainer::const_iterator it=gi->begin(); it!=gi->end();++it)
@@ -2041,7 +2041,7 @@ void dG3DDomain::computeAllIPStrain(AllIPState *aips,const unknownField *ufield,
else
mlawbulkp= static_cast<dG3DMaterialLaw*>(mlawPlus);
if (mlawbulkm->getUseBarF()!=mlawbulkp->getUseBarF()) //check with the bulk law
- Msg::Error("dG3DDomain::computeIpv All the constitutive laws need to use Fbar or not");
+ Msg::Error("dG3DDomain::computeAllIPStrain All the constitutive laws need to use Fbar or not");
bool useBarF=mlawbulkm->getUseBarF();
//
for(elementGroup::elementContainer::const_iterator it=gi->begin(); it!=gi->end();++it)
@@ -2274,7 +2274,7 @@ void dG3DDomain::computeAllIPStrain(AllIPState *aips,const unknownField *ufield,
};
-void dG3DDomain::computeIPVariable(AllIPState *aips,const unknownField *ufield, const IPStateBase::whichState ws, bool stiff)
+bool dG3DDomain::computeIPVariable(AllIPState *aips,const unknownField *ufield, const IPStateBase::whichState ws, bool stiff)
{
_evalStiff = false; // directly access to computeIPv is for matrix by perturbation otherwise compute ipvarible thanks to this function
if (ws == IPStateBase::initial)
@@ -2286,7 +2286,8 @@ void dG3DDomain::computeIPVariable(AllIPState *aips,const unknownField *ufield,
#if defined(HAVE_MPI)
if ((this->getMaterialLawMinus()->isNumeric() or this->getMaterialLawPlus()->isNumeric()) and this->_otherRanks.size()>0)
{
- this->computeIPVariableMPI(aips,ufield,ws,stiff);
+ bool ok = this->computeIPVariableMPI(aips,ufield,ws,stiff);
+ return ok;
}
else
#endif //HAVE_MPI
@@ -2311,7 +2312,11 @@ void dG3DDomain::computeIPVariable(AllIPState *aips,const unknownField *ufield,
dispp.resize(R.size());
ufield->get(R,dispp);
int npts = integBound->getIntPoints(it->second,&GP);
- this->computeIpv(aips,ie, GP,ws,getMinusDomain(),getPlusDomain(),getMaterialLawMinus(),getMaterialLawPlus(),dispm,dispp,false,stiff);
+ bool success = this->computeIpv(aips,ie, GP,ws,getMinusDomain(),getPlusDomain(),getMaterialLawMinus(),getMaterialLawPlus(),dispm,dispp,false,stiff);
+ if (!success)
+ {
+ return false;
+ }
}
// bulk
for (elementGroup::elementContainer::const_iterator it = g->begin(); it != g->end(); ++it)
@@ -2321,13 +2326,18 @@ void dG3DDomain::computeIPVariable(AllIPState *aips,const unknownField *ufield,
getFunctionSpace()->getKeys(e,R);
dispm.resize(R.size());
ufield->get(R,dispm);
- this->computeIpv(aips,e,ws,getMaterialLaw(),dispm,stiff);
+ bool success = this->computeIpv(aips,e,ws,getMaterialLaw(),dispm,stiff);
+ if (!success)
+ {
+ return false;
+ }
}
}
}
_evalStiff = true;
-}
+ return true;
+};
void dG3DDomain::computeStrain(MElement *e, const int npts_bulk, IntPt *GP,
AllIPState::ipstateElementContainer *vips,
@@ -2588,10 +2598,10 @@ void dG3DDomain::computeStrain(MElement *e, const int npts_bulk, IntPt *GP,
}
}
}
-void dG3DDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
+bool dG3DDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
materialLaw *mlaw__,fullVector<double> &disp, bool stiff)
{
- if (!getElementErosionFilter()(e)) return;
+ if (!getElementErosionFilter()(e)) return true;
dG3DMaterialLaw *mlaw;
AllIPState::ipstateElementContainer *vips = aips->getIPstate(e->getNum());
@@ -2716,7 +2726,14 @@ void dG3DDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichStat
{
Msg::Error("This substepping method %d has not been implemented",_subSteppingMethod);
}
+
+ if (!ipv->getConstitutiveSuccessFlag())
+ {
+ Msg::Error("constitive estimation is failed at ele %d GP %d",e->getNum(), j);
+ return false;
+ }
}
+ return true;
}
void dG3DDomain::setValuesForBodyForce(AllIPState *aips, const IPStateBase::whichState ws){
@@ -3710,14 +3727,14 @@ void dG3DDomain::computeStrain(AllIPState *aips,MInterfaceElement *ie, IntPt *GP
Msg::Error("No virtual interface for 3D");
}
}
-void dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
+bool dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
partDomain* efMinus, partDomain *efPlus,materialLaw *mlawmi,
materialLaw *mlawp,fullVector<double> &dispm,
fullVector<double> &dispp,
const bool virt, bool stiff,const bool checkfrac)
{
MElement *ele = dynamic_cast<MElement*>(ie);
- if (!getElementErosionFilter()(ele)) return;
+ if (!getElementErosionFilter()(ele)) return true;
if(!virt)
{
@@ -3736,7 +3753,10 @@ void dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,co
else
mlawbulkp= static_cast<dG3DMaterialLaw*>(mlawplus);
if (mlawbulkm->getUseBarF()!=mlawbulkp->getUseBarF()) //check with the bulk law
+ {
Msg::Error("dG3DDomain::computeIpv All the constitutive laws need to use Fbar or not");
+ return false;
+ }
bool useBarF=mlawbulkm->getUseBarF();
if(_evalStiff){ this->computeStrain(aips,ie, GP,ws,efMinus,efPlus,dispm,dispp,virt, useBarF);}
@@ -3762,34 +3782,35 @@ void dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,co
if (_subSteppingMethod == 0){
// no substepting
- if (mlawminus->getNum() != mlawplus->getNum()){
- // if minus and plus material are different
- if (_planeStressState){
- mlawminus->stress3DTo2D(ipvm,ipvmprev, stiff, checkfrac,needDTangent);
- mlawplus->stress3DTo2D(ipvp,ipvpprev, stiff, checkfrac,needDTangent);
- }
- else{
- mlawminus->stress(ipvm,ipvmprev, stiff, checkfrac,needDTangent);
- mlawplus->stress(ipvp,ipvpprev, stiff, checkfrac,needDTangent);
- }
+ if (mlawminus->getNum() != mlawplus->getNum())
+ {
+ // if minus and plus material are different
+ if (_planeStressState){
+ mlawminus->stress3DTo2D(ipvm,ipvmprev, stiff, checkfrac,needDTangent);
+ mlawplus->stress3DTo2D(ipvp,ipvpprev, stiff, checkfrac,needDTangent);
+ }
+ else{
+ mlawminus->stress(ipvm,ipvmprev, stiff, checkfrac,needDTangent);
+ mlawplus->stress(ipvp,ipvpprev, stiff, checkfrac,needDTangent);
+ }
}
else{
// if both part are the same
- if (_planeStressState)
- mlawminus->stress3DTo2D(ipvm,ipvmprev,stiff,checkfrac,needDTangent);
- else
+ if (_planeStressState)
+ mlawminus->stress3DTo2D(ipvm,ipvmprev,stiff,checkfrac,needDTangent);
+ else
mlawminus->stress(ipvm,ipvmprev,stiff,checkfrac,needDTangent);
- if (_averageStrainBased){
- // if the computation is performed based on the average strain
- // positive IP is equal to minus IP
+ if (_averageStrainBased){
+ // if the computation is performed based on the average strain
+ // positive IP is equal to minus IP
ipvp->operator=(*(dynamic_cast<const IPVariable*>(ipvm)));
- }
- else{
- if (_planeStressState)
- mlawplus->stress3DTo2D(ipvp,ipvpprev,stiff,checkfrac,needDTangent);
- else
- mlawplus->stress(ipvp,ipvpprev,stiff,checkfrac,needDTangent);
+ }
+ else{
+ if (_planeStressState)
+ mlawplus->stress3DTo2D(ipvp,ipvpprev,stiff,checkfrac,needDTangent);
+ else
+ mlawplus->stress(ipvp,ipvpprev,stiff,checkfrac,needDTangent);
}
}
}
@@ -3834,17 +3855,17 @@ void dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,co
if (iter == _subStepNumber){
// in last substep
- if (_planeStressState)
- mlawminus->stress3DTo2D(ipvm,ipvTemp,stiff,checkfrac,needDTangent);
- else
- mlawminus->stress(ipvm,ipvTemp,stiff,checkfrac,needDTangent);
+ if (_planeStressState)
+ mlawminus->stress3DTo2D(ipvm,ipvTemp,stiff,checkfrac,needDTangent);
+ else
+ mlawminus->stress(ipvm,ipvTemp,stiff,checkfrac,needDTangent);
break;
}
else{
- if (_planeStressState)
+ if (_planeStressState)
mlawminus->stress3DTo2D(ipvm,ipvTemp,false,checkfrac,needDTangent);
- else
- mlawminus->stress(ipvm,ipvTemp,false,checkfrac,needDTangent);
+ else
+ mlawminus->stress(ipvm,ipvTemp,false,checkfrac,needDTangent);
dgipvTemp->operator=(*(dynamic_cast<const IPVariable*>(ipvm)));
}
// next step
@@ -3897,18 +3918,18 @@ void dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,co
if (iter == _subStepNumber){
// in last substep
- if (_planeStressState)
- mlawplus->stress3DTo2D(ipvp,ipvTemp,stiff,checkfrac,needDTangent);
- else
- mlawplus->stress(ipvp,ipvTemp,stiff,checkfrac,needDTangent);
+ if (_planeStressState)
+ mlawplus->stress3DTo2D(ipvp,ipvTemp,stiff,checkfrac,needDTangent);
+ else
+ mlawplus->stress(ipvp,ipvTemp,stiff,checkfrac,needDTangent);
break;
}
else{
- if (_planeStressState)
+ if (_planeStressState)
mlawplus->stress3DTo2D(ipvp,ipvTemp,false,checkfrac,needDTangent);
- else
- mlawplus->stress(ipvp,ipvTemp,false,checkfrac,needDTangent);
- dgipvTemp->operator=(*(dynamic_cast<const IPVariable*>(ipvp)));
+ else
+ mlawplus->stress(ipvp,ipvTemp,false,checkfrac,needDTangent);
+ dgipvTemp->operator=(*(dynamic_cast<const IPVariable*>(ipvp)));
}
// next step
iter++;
@@ -3935,6 +3956,12 @@ void dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,co
}
else{ Msg::Error("Moduli used to compute body force is not defined");}
}
+
+ if (!(ipvm->getConstitutiveSuccessFlag()) or !(ipvp->getConstitutiveSuccessFlag()))
+ {
+ Msg::Error("Constitutive evaluation is failed at interface ele %d, GP %d %d",ele->getNum(), j, j+npts);
+ return false;
+ }
}
}
else
@@ -3942,6 +3969,7 @@ void dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,co
Msg::Error("No virtual interface for 3D");
}
+ return true;
}
int dG3DDomain::getStressDimension() const{
@@ -5443,7 +5471,7 @@ MElement* dG3DDomain::createInterface(IElement *ie1, IElement *ie2) const
#if defined(HAVE_MPI)
-void dG3DDomain::computeIPVariableMPI(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff){
+bool dG3DDomain::computeIPVariableMPI(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff){
IntPt *GP;
MPI_Status status;
if (Msg::GetCommRank() == _rootRank){
@@ -5765,6 +5793,7 @@ void dG3DDomain::computeIPVariableMPI(AllIPState *aips,const unknownField *ufiel
}
}
};
+ return true;
};
#endif
diff --git a/dG3D/src/dG3DDomain.h b/dG3D/src/dG3DDomain.h
index 9a155091da846f32a3493d53bc5010a05332ae93..84a883ac31c1cff6f55da32c0a2e1a0df80280fd 100644
--- a/dG3D/src/dG3DDomain.h
+++ b/dG3D/src/dG3DDomain.h
@@ -177,13 +177,13 @@ class dG3DDomain : public dgPartDomain{
virtual void computeStrain(MElement *e, const int npts, IntPt *GP,
AllIPState::ipstateElementContainer *vips,
IPStateBase::whichState ws, fullVector<double> &disp, bool useBarF) const;
- virtual void computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
+ virtual bool computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
materialLaw *mlaw,fullVector<double> &disp, bool stiff);
virtual void computeStrain(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
partDomain* efMinus, partDomain *efPlus, fullVector<double> &dispm,
fullVector<double> &dispp, const bool virt,
bool useBarF);
- virtual void computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
+ virtual bool computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
partDomain* efMinus, partDomain *efPlus,materialLaw *mlawminus,
materialLaw *mlawplus,fullVector<double> &dispm,
fullVector<double> &dispp, const bool virt, bool stiff,const bool checkfrac=true);
@@ -200,7 +200,7 @@ class dG3DDomain : public dgPartDomain{
virtual void setdFmdGM(AllIPState *aips, std::map<Dof, STensor33> &dUdG,const IPStateBase::whichState ws);
virtual void setValuesForBodyForce(AllIPState *aips, const IPStateBase::whichState ws);
virtual void computeAllIPStrain(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, const bool virt);
- virtual void computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
+ virtual bool computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
virtual void setMaterialLaw(const std::map<int,materialLaw*> &maplaw);
virtual materialLaw* getMaterialLaw(){return _mlaw;}
virtual const materialLaw* getMaterialLaw() const{return _mlaw;}
@@ -257,7 +257,7 @@ class dG3DDomain : public dgPartDomain{
virtual QuadratureBase* getQuadratureRulesForNeumannBC(nonLinearNeumannBC &neu) const;
virtual partDomain* clone() const {return new dG3DDomain(*this);};
#if defined(HAVE_MPI)
- virtual void computeIPVariableMPI(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
+ virtual bool computeIPVariableMPI(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
#endif
virtual double getNonLocalStabilityParameter() const { return _nonLocalBeta;}
virtual bool getNonLocalContinuity() const { return _nonLocalContinuity;}
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 19d1a3d283f99a7df6d5e0fb7849dbfa208bee3a..2df4b80e90cca1781bdf87293f8707daad8b3da0 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -10573,8 +10573,8 @@ void NonLinearTVMDG3DMaterialLaw::setTensionCompressionRegularisation(const doub
void NonLinearTVMDG3DMaterialLaw::setExtraBranchNLType(const int method){
_viscoLaw.setExtraBranchNLType(method);
};
-void NonLinearTVMDG3DMaterialLaw::setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2){
- _viscoLaw.setCorrectionsAllBranchesTVE(i,V1,V2,D1,D2);
+void NonLinearTVMDG3DMaterialLaw::setCorrectionsAllBranchesTVE(const int i, const double V0,const double V1, const double V2, const double D0,const double D1, const double D2){
+ _viscoLaw.setCorrectionsAllBranchesTVE(i,V0,V1,V2,D0,D1,D2);
};
void NonLinearTVMDG3DMaterialLaw::setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
_viscoLaw.setAdditionalCorrectionsAllBranchesTVE(i,V3,V4,V5,D3,D4,D5);
@@ -11243,8 +11243,8 @@ void NonLinearTVPDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp
NonLinearTVENonLinearTVPDG3DMaterialLaw::NonLinearTVENonLinearTVPDG3DMaterialLaw(const int num, const double rho, const double E, const double nu, const double tol,
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,tol,Tinitial,Alpha, KThCon, Cp, matrixbyPerturbation,pert, thermalEstimationPreviousConfig)
+ const bool matrixbyPerturbation, const double pert, const double stressIteratorTol, const bool thermalEstimationPreviousConfig):
+ dG3DMaterialLaw(num,rho), _viscoLaw(num,E,nu,rho,tol,Tinitial,Alpha, KThCon, Cp, matrixbyPerturbation,pert, stressIteratorTol, thermalEstimationPreviousConfig)
{
fillElasticStiffness(E, nu, elasticStiffness); // Why are some of the variables defined in this constructor? --- WHAT??
};
@@ -11404,8 +11404,8 @@ void NonLinearTVENonLinearTVPDG3DMaterialLaw::setExtraBranchNLType(const int met
void NonLinearTVENonLinearTVPDG3DMaterialLaw::setTensionCompressionRegularisation(const double tensionCompressionRegularisation){
_viscoLaw.setTensionCompressionRegularisation(tensionCompressionRegularisation);
};
-void NonLinearTVENonLinearTVPDG3DMaterialLaw::setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2){
- _viscoLaw.setCorrectionsAllBranchesTVE(i,V1,V2,D1,D2);
+void NonLinearTVENonLinearTVPDG3DMaterialLaw::setCorrectionsAllBranchesTVE(const int i, const double V0,const double V1, const double V2, const double D0,const double D1, const double D2){
+ _viscoLaw.setCorrectionsAllBranchesTVE(i,V0,V1,V2,D0,D1,D2);
};
void NonLinearTVENonLinearTVPDG3DMaterialLaw::setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
_viscoLaw.setAdditionalCorrectionsAllBranchesTVE(i,V3,V4,V5,D3,D4,D5);
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 1a598d0735df5a45b4be6dc37bd910eb2c33085c..fbc79e13836c29ba1a4b7439c81e4d8b09015dc8 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -3329,7 +3329,7 @@ class NonLinearTVMDG3DMaterialLaw : public dG3DMaterialLaw{ // public Material
void setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3);
void setExtraBranchNLType(const int type);
void setTensionCompressionRegularisation(const double tensionCompressionRegularisation);
- void setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2);
+ void setCorrectionsAllBranchesTVE(const int i, const double V0,const double V1, const double V2, const double D0, const double D1, const double D2);
void setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
void setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci);
@@ -3602,7 +3602,7 @@ class NonLinearTVENonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{ // pub
public:
NonLinearTVENonLinearTVPDG3DMaterialLaw(const int num, const double rho, const double E, const double nu, 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 = false);
+ const bool matrixbyPerturbation = false, const double pert = 1e-8, const double stressIteratorTol = 1e-9, const bool thermalEstimationPreviousConfig = false);
void setStrainOrder(const int order);
void setViscoelasticMethod(const int method);
@@ -3639,7 +3639,7 @@ class NonLinearTVENonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{ // pub
void setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3);
void setExtraBranchNLType(const int type);
void setTensionCompressionRegularisation(const double tensionCompressionRegularisation);
- void setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2);
+ void setCorrectionsAllBranchesTVE(const int i, const double V0,const double V1, const double V2, const double D0, const double D1, const double D2);
void setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
void setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci);
diff --git a/dG3D/src/nonLocalDamageDG3DIPVariable.h b/dG3D/src/nonLocalDamageDG3DIPVariable.h
index a8ceb589db42d43a352c52a118df2de93768cdc4..fb78db048e1d0a6166c97caa3d2f8c16755bb870 100644
--- a/dG3D/src/nonLocalDamageDG3DIPVariable.h
+++ b/dG3D/src/nonLocalDamageDG3DIPVariable.h
@@ -1942,7 +1942,10 @@ class nonLocalPorosityDG3DIPVariable : public dG3DIPVariable
return _nldPorousipv->getRefToDIrreversibleEnergyDNonLocalVariable(index);
};
-
+ virtual bool getConstitutiveSuccessFlag() const
+ {
+ return _nldPorousipv->getConstitutiveSuccessFlag();
+ }
virtual void restart();
virtual void setRefToLinearConductivity(const STensor3 &eT) {setConstRefToLinearK(eT,0,0);}
virtual void setRefToStiff_alphaDilatation(const STensor3 &eT) {setConstRefToLinearSymmetrizationCoupling(eT,0);}
diff --git a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
index ae1a7e3e1d283653c484c185bf3bec2d8f4bbaf6..3ed76344f3eeb66231db3d496213a3c839f08df4 100644
--- a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
+++ b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
@@ -1449,6 +1449,15 @@ void NonLocalDamageGursonHill48DG3DMaterialLaw::setYieldParameters(const std::ve
{
static_cast<mlawNonLocalDamageGursonHill48*>(_nldPorous)->setYieldParameters(params);
}
+void NonLocalDamageGursonHill48DG3DMaterialLaw::setYieldParameters(const fullVector<double>& params)
+{
+ std::vector<double> params_vec(params.size());
+ for (int i=0; i< params.size(); i++)
+ {
+ params_vec[i] = params(i);
+ }
+ static_cast<mlawNonLocalDamageGursonHill48*>(_nldPorous)->setYieldParameters(params_vec);
+};
void NonLocalDamageGursonHill48DG3DMaterialLaw::setTemperatureFunction_q1(const scalarFunction& Tfunc)
{
diff --git a/dG3D/src/nonLocalDamageDG3DMaterialLaw.h b/dG3D/src/nonLocalDamageDG3DMaterialLaw.h
index 8f927536996bbcd22c12b619fb9f7485280adb41..2dd6ef67368e84db72cdcfd7ee7373bc2be8b541 100644
--- a/dG3D/src/nonLocalDamageDG3DMaterialLaw.h
+++ b/dG3D/src/nonLocalDamageDG3DMaterialLaw.h
@@ -508,6 +508,7 @@ public:
void setTemperatureFunction_q3(const scalarFunction& Tfunc);
void setYieldParameters(const std::vector<double>& params);
+ void setYieldParameters(const fullVector<double>& params);
#ifndef SWIG
NonLocalDamageGursonHill48DG3DMaterialLaw(const NonLocalDamageGursonHill48DG3DMaterialLaw &source);
diff --git a/dgshell/src/dgShellDomain.cpp b/dgshell/src/dgShellDomain.cpp
index 7155664f71642cf47fb004d1e7ee19689649951a..b9bc07a895d3b034109002e04c639d9feab02f44 100644
--- a/dgshell/src/dgShellDomain.cpp
+++ b/dgshell/src/dgShellDomain.cpp
@@ -369,7 +369,7 @@ void dgLinearShellDomain::ensureSameValuesBothSideIPvariable(IPField *ip)
}
}
-void dgLinearShellDomain::computeIPVariable(AllIPState *aips,const unknownField *ufield, const IPStateBase::whichState ws, bool stiff)
+bool dgLinearShellDomain::computeIPVariable(AllIPState *aips,const unknownField *ufield, const IPStateBase::whichState ws, bool stiff)
{
IntPt *GP;
fullVector<double> dispm;
@@ -388,7 +388,8 @@ void dgLinearShellDomain::computeIPVariable(AllIPState *aips,const unknownField
dispp.resize(R.size());
ufield->get(R,dispp);
int npts = integBound->getIntPoints(ie,&GP);
- this->computeIpv(aips,ie, GP,ws,dom,dom,_mlaw,_mlaw,dispm,dispp,false, stiff);
+ bool ok = this->computeIpv(aips,ie, GP,ws,dom,dom,_mlaw,_mlaw,dispm,dispp,false, stiff);
+ if (!ok) return false;
}
// Virtual interface element
for(elementGroup::elementContainer::const_iterator it=gib->begin(); it!=gib->end();++it){
@@ -397,7 +398,8 @@ void dgLinearShellDomain::computeIPVariable(AllIPState *aips,const unknownField
_spaceMinus->getKeys(ie->getElem(0),R);
dispm.resize(R.size());
ufield->get(R,dispm);
- this->computeIpv(aips,ie, GP,ws,dom,dom,_mlaw,_mlaw,dispm,dispm,true, stiff);
+ bool ok = this->computeIpv(aips,ie, GP,ws,dom,dom,_mlaw,_mlaw,dispm,dispm,true, stiff);
+ if (!ok) return false;
}
// bulk
@@ -407,11 +409,14 @@ void dgLinearShellDomain::computeIPVariable(AllIPState *aips,const unknownField
_space->getKeys(e,R);
dispm.resize(R.size());
ufield->get(R,dispm);
- this->computeIpv(aips,e,ws,_mlaw,dispm, stiff);
+ bool ok = this->computeIpv(aips,e,ws,_mlaw,dispm, stiff);
+ if (!ok) return false;
}
+
+ return true;
}
-void dgLinearShellDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,
+bool dgLinearShellDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,
const IPStateBase::whichState ws,partDomain* efMinus, partDomain *efPlus,
materialLaw *mlawmi, materialLaw *mlawp,
fullVector<double> &dispm,
@@ -622,9 +627,10 @@ void dgLinearShellDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, Int
mbulk->stress(ips, stiff);
}
}
+ return true;
}
-void dgLinearShellDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,materialLaw *mlaw__,
+bool dgLinearShellDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,materialLaw *mlaw__,
fullVector<double> &disp, bool stiff)
{
shellMaterialLaw *mlaw = static_cast<shellMaterialLaw*>(mlaw__);
@@ -649,6 +655,7 @@ void dgLinearShellDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::
if(ws != IPStateBase::initial)
mlaw->stress(ips, stiff);
}
+ return true;
}
FilterDof* dgLinearShellDomain::createFilterDof(const int comp) const
@@ -867,7 +874,7 @@ dgNonLinearShellDomain::dgNonLinearShellDomain(const dgNonLinearShellDomain &sou
_evalStiff(source._evalStiff){}
-void dgNonLinearShellDomain::computeIPVariable(AllIPState *aips,const unknownField *ufield, const IPStateBase::whichState ws, bool stiff)
+bool dgNonLinearShellDomain::computeIPVariable(AllIPState *aips,const unknownField *ufield, const IPStateBase::whichState ws, bool stiff)
{
_evalStiff = false; // directly access to computeIPv is for matrix by perturbation otherwise compute ipvarible thanks to this function
if(ws == IPStateBase::initial){
@@ -891,7 +898,8 @@ void dgNonLinearShellDomain::computeIPVariable(AllIPState *aips,const unknownFie
dispp.resize(R.size());
ufield->get(R,dispp);
int npts = integBound->getIntPoints(ie,&GP);
- this->computeIpv(aips,ie, GP,ws,dom,dom,_mlaw,_mlaw,dispm,dispp,false, stiff);
+ bool ok = this->computeIpv(aips,ie, GP,ws,dom,dom,_mlaw,_mlaw,dispm,dispp,false, stiff);
+ if (!ok) return false;
}
// Virtual interface element
for(elementGroup::elementContainer::const_iterator it=gib->begin(); it!=gib->end();++it){
@@ -900,7 +908,8 @@ void dgNonLinearShellDomain::computeIPVariable(AllIPState *aips,const unknownFie
_spaceMinus->getKeys(ie->getElem(0),R);
dispm.resize(R.size());
ufield->get(R,dispm);
- this->computeIpv(aips,ie, GP,ws,dom,dom,_mlaw,_mlaw,dispm,dispm,true, stiff);
+ bool ok = this->computeIpv(aips,ie, GP,ws,dom,dom,_mlaw,_mlaw,dispm,dispm,true, stiff);
+ if (!ok) return false;
}
// bulk
@@ -910,13 +919,15 @@ void dgNonLinearShellDomain::computeIPVariable(AllIPState *aips,const unknownFie
_space->getKeys(e,R);
dispm.resize(R.size());
ufield->get(R,dispm);
- this->computeIpv(aips,e,ws,_mlaw,dispm, stiff);
+ bool ok = this->computeIpv(aips,e,ws,_mlaw,dispm, stiff);
+ if (!ok) return false;
}
}
_evalStiff = true;
+ return true;
}
-void dgNonLinearShellDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
+bool dgNonLinearShellDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
materialLaw *mlaw__,fullVector<double> &disp, bool stiff)
{
dgNonLinearShellMaterialLaw *mlaw;
@@ -945,8 +956,9 @@ void dgNonLinearShellDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBas
mlaw->ensurePlaneStress(ips,_evalStiff);
mlaw->stress(ips, stiff);
}
+ return true;
}
-void dgNonLinearShellDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
+bool dgNonLinearShellDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
partDomain* efMinus, partDomain *efPlus,materialLaw *mlawmi,
materialLaw *mlawp,fullVector<double> &dispm,
fullVector<double> &dispp,
@@ -1207,6 +1219,7 @@ void dgNonLinearShellDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie,
mbulk->stress(ips,stiff, false); // no fracture at extremities for now
}
}
+ return true;
}
void dgNonLinearShellDomain::initialIPVariable(AllIPState *aips,const unknownField *ufield,
@@ -1687,7 +1700,7 @@ void interDomainBetweenShell::setGaussIntegrationRule()
integBulk = NULL; // no element
}
-void interDomainBetweenShell::computeIPVariable(AllIPState *aips,const unknownField *ufield,
+bool interDomainBetweenShell::computeIPVariable(AllIPState *aips,const unknownField *ufield,
const IPStateBase::whichState ws, bool stiff)
{
IntPt *GP;
@@ -1709,8 +1722,10 @@ void interDomainBetweenShell::computeIPVariable(AllIPState *aips,const unknownFi
ufield->get(R,dispp);
// shellMaterialLaw* smlawminus = dynamic_cast<shellMaterialLaw*>(_mlawMinus);
// shellMaterialLaw* smlawplus = dynamic_cast<shellMaterialLaw*>(_mlawPlus);
- this->computeIpv(aips,iele,GP,ws,_domMinus,_domPlus,mlawMinus,mlawPlus,dispm,dispp,false, stiff);
+ bool ok = this->computeIpv(aips,iele,GP,ws,_domMinus,_domPlus,mlawMinus,mlawPlus,dispm,dispp,false, stiff);
+ if (!ok) return false;
}
+ return true;
}
void interDomainBetweenShell::setMaterialLaw(const std::map<int,materialLaw*>& maplaw)
diff --git a/dgshell/src/dgShellDomain.h b/dgshell/src/dgShellDomain.h
index 31f14b23b5f3eb35f89384ebe2d1b4d5305841e9..2e05e568cf41dc5fd138cf90a89d6f5ce2874d0b 100644
--- a/dgshell/src/dgShellDomain.h
+++ b/dgshell/src/dgShellDomain.h
@@ -52,14 +52,14 @@ class dgLinearShellDomain : public dgPartDomain{
virtual void setdFmdGM(AllIPState *aips, std::map<Dof, STensor33> &dUdG,const IPStateBase::whichState ws){};
virtual void setValuesForBodyForce(AllIPState *aips, const IPStateBase::whichState ws){};
- virtual void computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
+ virtual bool computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
materialLaw *mlaw,fullVector<double> &disp, bool stiff);
- virtual void computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
+ virtual bool computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
partDomain* efMinus, partDomain *efPlus,materialLaw *mlawminus,
materialLaw *mlawplus,fullVector<double> &dispm,
fullVector<double> &dispp,
const bool virt, bool stiff, const bool checkfrac=true);
- virtual void computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
+ virtual bool computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
virtual void checkFailure(IPField* ipf) const{
// Msg::Warning("This function checkFailure is not defined ");
}
@@ -120,14 +120,14 @@ class dgNonLinearShellDomain : public dgLinearShellDomain{
dgNonLinearShellDomain(const dgNonLinearShellDomain &source);
~dgNonLinearShellDomain(){}
#ifndef SWIG
- virtual void computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
+ virtual bool computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
virtual void setDeformationGradientGradient(AllIPState *aips, const STensor33 &GM,const IPStateBase::whichState ws){};
virtual void setdFmdFM(AllIPState *aips, std::map<Dof, STensor3> &dUdF,const IPStateBase::whichState ws){};
virtual void setdFmdGM(AllIPState *aips, std::map<Dof, STensor33> &dUdG,const IPStateBase::whichState ws){};
- virtual void computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
+ virtual bool computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
materialLaw *mlaw,fullVector<double> &disp, bool stiff);
- virtual void computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
+ virtual bool computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
partDomain* efMinus, partDomain *efPlus,materialLaw *mlawminus,
materialLaw *mlawplus,fullVector<double> &dispm,
fullVector<double> &dispp,
@@ -155,7 +155,7 @@ class interDomainBetweenShell : public dgLinearShellDomain{
#ifndef SWIG
// void initializeTerms(FunctionSpaceBase* space1_,FunctionSpaceBase* space2_,unknownField *uf,IPField*ip);
virtual void setGaussIntegrationRule();
- virtual void computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
+ virtual bool computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
virtual void setMaterialLaw(const std::map<int,materialLaw*>& maplaw);
virtual materialLaw* getMaterialLaw(){return _mlaw;}
virtual const materialLaw* getMaterialLaw() const{return _mlaw;}
diff --git a/dgshell/src/nonLinearInterDomainShell.cpp b/dgshell/src/nonLinearInterDomainShell.cpp
index 6a0a961e6fb42bfd7f92b1fb64058c05e87dca0c..0f8b73d796711f8bba8a2b2e9340be576b9c2c35 100644
--- a/dgshell/src/nonLinearInterDomainShell.cpp
+++ b/dgshell/src/nonLinearInterDomainShell.cpp
@@ -138,7 +138,7 @@ void nonLinearInterDomainShell::setGaussIntegrationRule()
}
-void nonLinearInterDomainShell::computeIPVariable(AllIPState *aips,const unknownField *ufield,
+bool nonLinearInterDomainShell::computeIPVariable(AllIPState *aips,const unknownField *ufield,
const IPStateBase::whichState ws, bool stiff)
{
IntPt *GP;
@@ -158,7 +158,8 @@ void nonLinearInterDomainShell::computeIPVariable(AllIPState *aips,const unknown
_spacePlus->getKeys(iele->getElem(1),R);
dispp.resize(R.size());
ufield->get(R,dispp);
- this->computeIpv(aips,iele,GP,ws,_domMinus,_domPlus,mlawMinus,mlawPlus,dispm,dispp,false, stiff);
+ bool ok = this->computeIpv(aips,iele,GP,ws,_domMinus,_domPlus,mlawMinus,mlawPlus,dispm,dispp,false, stiff);
+ if (!ok) return false;
// set initial value change this
if(ws == IPStateBase::initial)
{
diff --git a/dgshell/src/nonLinearInterDomainShell.h b/dgshell/src/nonLinearInterDomainShell.h
index 91604a34fba5237b41aabe8411d06534d1c06471..06ce826bb06c02a656903275ee32fcb41d9a9d5c 100644
--- a/dgshell/src/nonLinearInterDomainShell.h
+++ b/dgshell/src/nonLinearInterDomainShell.h
@@ -31,7 +31,7 @@ class nonLinearInterDomainShell : public dgNonLinearShellDomain{
~nonLinearInterDomainShell(){}
// void initializeTerms(FunctionSpaceBase* space1_,FunctionSpaceBase* space2_,unknownField *uf,IPField*ip);
virtual void setGaussIntegrationRule();
- virtual void computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
+ virtual bool computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, bool stiff);
virtual void setMaterialLaw(const std::map<int,materialLaw*>& maplaw);
virtual materialLaw* getMaterialLaw(){return _mlaw;}
virtual const materialLaw* getMaterialLaw() const{return _mlaw;}