diff --git a/NonLinearSolver/internalPoints/ipField.cpp b/NonLinearSolver/internalPoints/ipField.cpp
index 770e10689a0cb36ff91e9d09356e27547ee99bf5..1e232cc6b7caae0655245948a59f4e9a174b5e2e 100644
--- a/NonLinearSolver/internalPoints/ipField.cpp
+++ b/NonLinearSolver/internalPoints/ipField.cpp
@@ -941,6 +941,7 @@ std::string IPField::ToString(const int i){
else if (i == FveAM_ZX) return "FveAM_ZX";
else if (i == epsAM) return "epsAM";
else if (i == epsR) return "epsR";
+ else if (i == epsG) return "epsG";
else if (i == Tcrys) return "Tcrys";
else if (i == Tmelt) return "Tmelt";
diff --git a/NonLinearSolver/internalPoints/ipField.h b/NonLinearSolver/internalPoints/ipField.h
index 4f0ddd782a90db25883f41682b3ee5b97ab71b96..a9e276e0391293854ae7ff4d7ce0fd7d056a993e 100644
--- a/NonLinearSolver/internalPoints/ipField.h
+++ b/NonLinearSolver/internalPoints/ipField.h
@@ -302,7 +302,7 @@ class IPField : public elementsField {
FP_XX, FP_YY, FP_ZZ, FP_XY, FP_YX, FP_YZ, FP_ZY, FP_ZX, FP_XZ,
FpR_XX, FpR_YY, FpR_ZZ, FpR_XY, FpR_YZ, FpR_ZX,
FveR_XX, FveR_YY, FveR_ZZ, FveR_XY, FveR_YZ, FveR_ZX,
- ZG, PDF, TT, WT,
+ epsG, ZG, PDF, TT, WT,
FthAM_XX, FthAM_YY, FthAM_ZZ, FthAM_XY, FthAM_YZ, FthAM_ZX,
FfAM_XX, FfAM_YY, FfAM_ZZ, FfAM_XY, FfAM_YZ, FfAM_ZX,
FpAM_XX, FpAM_YY, FpAM_ZZ, FpAM_XY, FpAM_YZ, FpAM_ZX,
diff --git a/NonLinearSolver/internalPoints/ipHyperelastic.cpp b/NonLinearSolver/internalPoints/ipHyperelastic.cpp
index 93a02553cb5cdaca09f9b49a77bd786cfe44ec5f..605fe07dace5f0e46b2fe3e3fc60c3583aa53578 100644
--- a/NonLinearSolver/internalPoints/ipHyperelastic.cpp
+++ b/NonLinearSolver/internalPoints/ipHyperelastic.cpp
@@ -13,8 +13,8 @@
#include "restartManager.h"
#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.)
+IPHyperViscoElastic::IPHyperViscoElastic(const int N):IPVariableMechanics(),_N(N),_elasticEnergy(0.),_Ee(0.),_kirchhoff(0.),_kirchhoff_vol(0.),_kirchhoff_dev(0.),
+ _irreversibleEnergy(0.),_DirreversibleEnergyDF(0.), _viscousEnergyPart(0.), _dElasticEnergyPartdF(0.), _dViscousEnergyPartdF(0.), _elasticBulkPropertyScaleFactor(1.), _elasticShearPropertyScaleFactor(1.)
{
_A.clear();
_B.clear();
@@ -25,10 +25,11 @@ IPHyperViscoElastic::IPHyperViscoElastic(const int N):IPVariableMechanics(),_N(N
}
};
IPHyperViscoElastic::IPHyperViscoElastic(const IPHyperViscoElastic& src): IPVariableMechanics(src), _Ee(src._Ee),
- _N(src._N),_A(src._A),_B(src._B),_elasticEnergy(src._elasticEnergy), _kirchhoff(src._kirchhoff),
+ _N(src._N),_A(src._A),_B(src._B),_elasticEnergy(src._elasticEnergy), _kirchhoff(src._kirchhoff), _kirchhoff_vol(src._kirchhoff_vol), _kirchhoff_dev(src._kirchhoff_dev),
_irreversibleEnergy(src._irreversibleEnergy),_DirreversibleEnergyDF(src._DirreversibleEnergyDF),
- _viscousEnergyPart(src._viscousEnergyPart), _dElasticEnergyPartdF(src._dElasticEnergyPartdF),
- _dViscousEnergyPartdF(src._dViscousEnergyPartdF)
+ _viscousEnergyPart(src._viscousEnergyPart), _dElasticEnergyPartdF(src._dElasticEnergyPartdF),
+ _dViscousEnergyPartdF(src._dViscousEnergyPartdF), _elasticBulkPropertyScaleFactor(src._elasticBulkPropertyScaleFactor), _elasticShearPropertyScaleFactor(src._elasticShearPropertyScaleFactor)
+
{
};
@@ -39,6 +40,8 @@ IPHyperViscoElastic& IPHyperViscoElastic::operator =(const IPVariable& src){
if (psrc != NULL){
_Ee = psrc->_Ee;
_kirchhoff = psrc->_kirchhoff;
+ _kirchhoff_vol = psrc->_kirchhoff_vol;
+ _kirchhoff_dev = psrc->_kirchhoff_dev;
_N = psrc->_N;
_A = psrc->_A;
_B = psrc->_B;
@@ -48,6 +51,8 @@ IPHyperViscoElastic& IPHyperViscoElastic::operator =(const IPVariable& src){
_viscousEnergyPart = psrc->_viscousEnergyPart;
_dElasticEnergyPartdF = psrc->_dElasticEnergyPartdF;
_dViscousEnergyPartdF = psrc->_dViscousEnergyPartdF;
+ _elasticBulkPropertyScaleFactor= psrc->_elasticBulkPropertyScaleFactor;
+ _elasticShearPropertyScaleFactor= psrc->_elasticShearPropertyScaleFactor;
}
return *this;
};
@@ -60,11 +65,15 @@ void IPHyperViscoElastic::restart() {
restartManager::restart(_B);
restartManager::restart(_Ee);
restartManager::restart(_kirchhoff);
+ restartManager::restart(_kirchhoff_vol);
+ restartManager::restart(_kirchhoff_dev);
restartManager::restart(_irreversibleEnergy);
restartManager::restart(_DirreversibleEnergyDF);
restartManager::restart(_viscousEnergyPart);
restartManager::restart(_dElasticEnergyPartdF);
restartManager::restart(_dViscousEnergyPartdF);
+ restartManager::restart(_elasticBulkPropertyScaleFactor);
+ restartManager::restart(_elasticShearPropertyScaleFactor);
};
void IPHyperViscoElastic::getViscoElasticStrain(int i, STensor3& Ev) const
diff --git a/NonLinearSolver/internalPoints/ipHyperelastic.h b/NonLinearSolver/internalPoints/ipHyperelastic.h
index f6c3dc4e19cc3a48c2304e4eaa00634341644d40..49302c54f245254d4ae05a127f234bb5e289563a 100644
--- a/NonLinearSolver/internalPoints/ipHyperelastic.h
+++ b/NonLinearSolver/internalPoints/ipHyperelastic.h
@@ -25,14 +25,20 @@ class IPHyperViscoElastic : public IPVariableMechanics{
double _elasticEnergy; // elastic energy stored
STensor3 _Ee; // elastic strain
STensor3 _kirchhoff; // corotational Kirchhoff stress
+ STensor3 _kirchhoff_vol; // Volumetric corotational Kirchhoff stress
+ STensor3 _kirchhoff_dev; // Deviatoric corotational Kirchhoff stress
double _irreversibleEnergy;
STensor3 _DirreversibleEnergyDF;
-
+
+ double _elasticBulkPropertyScaleFactor;
+ double _elasticShearPropertyScaleFactor;
+
protected:
//For energy sources
- double _viscousEnergyPart;
- STensor3 _dElasticEnergyPartdF, _dViscousEnergyPartdF;
+ double _viscousEnergyPart, _p_vol, _Dppr_vol_DCepr;
+ STensor3 _dElasticEnergyPartdF, _dViscousEnergyPartdF, _K_dev_pr;
+ STensor43 _DKpr_dev_DCepr;
public:
IPHyperViscoElastic(const int N);
@@ -50,8 +56,18 @@ class IPHyperViscoElastic : public IPVariableMechanics{
virtual STensor3& getRefToCorotationalKirchhoffStress() {return _kirchhoff;};
virtual const STensor3& getConstRefToCorotationalKirchhoffStress() const {return _kirchhoff;};
+ virtual STensor3& getRefToCorotationalKirchhoffVolStress() {return _kirchhoff_vol;};
+ virtual const STensor3& getConstRefToCorotationalKirchhoffVolStress() const {return _kirchhoff_vol;};
+
+ virtual STensor3& getRefToCorotationalKirchhoffDevStress() {return _kirchhoff_dev;};
+ virtual const STensor3& getConstRefToCorotationalKirchhoffDevStress() const {return _kirchhoff_dev;};
+
virtual double defoEnergy() const{return _elasticEnergy;}
virtual double& getRefToElasticEnergy() {return _elasticEnergy;};
+ virtual double& getRefTo_p_vol() {return _p_vol;};
+ virtual STensor3& getRefTo_K_dev_pr() {return _K_dev_pr;};
+ virtual double& getRefTo_Dppr_vol_DCepr() {return _Dppr_vol_DCepr;};
+ virtual STensor43& getRefTo_DKpr_dev_DCepr() {return _DKpr_dev_DCepr;};
virtual double plasticEnergy() const{return 0.;}
virtual double irreversibleEnergy() const {return _irreversibleEnergy;};
@@ -77,7 +93,11 @@ class IPHyperViscoElastic : public IPVariableMechanics{
virtual const STensor3 &getConstRefToDElasticEnergyPartdF() const {return _dElasticEnergyPartdF;}
virtual STensor3 &getRefToDViscousEnergyPartdF() {return _dViscousEnergyPartdF;}
virtual const STensor3 &getConstRefToDViscousEnergyPartdF() const {return _dViscousEnergyPartdF;}
-
+
+ virtual double &getRefToElasticBulkPropertyScaleFactor() {return _elasticBulkPropertyScaleFactor;}
+ virtual const double &getConstRefToElasticBulkPropertyScaleFactor() const {return _elasticBulkPropertyScaleFactor;}
+ virtual double &getRefToElasticShearPropertyScaleFactor() {return _elasticShearPropertyScaleFactor;}
+ virtual const double &getConstRefToElasticShearPropertyScaleFactor() const {return _elasticShearPropertyScaleFactor;}
};
@@ -115,8 +135,9 @@ class IPHyperViscoElastoPlastic : public IPHyperViscoElastic{
bool _dissipationBlocked;
bool _dissipationActive;
protected:
- STensor3 _dPlasticEnergyPartdF;
-
+ STensor3 _dPlasticEnergyPartdF, _dptilde_VoldCepr, _K_dev, _dK_devdGamma;
+ STensor43 _dK_devdCepr;
+ double _ptilde_Vol, _X_vol, _dptilde_VoldGamma, _D_dev;
public:
IPHyperViscoElastoPlastic(const J2IsotropicHardening* comp,
@@ -136,6 +157,18 @@ class IPHyperViscoElastoPlastic : public IPHyperViscoElastic{
virtual double& getRefToPlasticEnergy() {return _plasticEnergy;};
virtual double plasticEnergy() const{return _plasticEnergy;}
+
+
+ virtual double& getRefTo_ptilde_Vol() {return _ptilde_Vol;};
+ virtual double& getRefTo_X_vol() {return _X_vol;};
+ virtual double& getRefTo_dptilde_VoldGamma() {return _dptilde_VoldGamma;};
+ virtual STensor3& getRefTo_dptilde_VoldCepr() {return _dptilde_VoldCepr;};
+ virtual double& getRefTo_D_dev() {return _D_dev;};
+ virtual STensor3& getRefTo_K_dev() {return _K_dev;};
+ virtual STensor3& getRefTo_dK_devdGamma() {return _dK_devdGamma;};
+ virtual STensor43& getRefTo_dK_devdCepr() {return _dK_devdCepr;};
+
+
virtual void restart();
diff --git a/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.cpp b/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.cpp
index 3ed569e504008510754a301c0e437fd6169c7917..33e6b11e7a4926904a859d4d8dd5d5c1092de81a 100644
--- a/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.cpp
+++ b/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.cpp
@@ -75,6 +75,7 @@ IPPhenomenologicalSMP::IPPhenomenologicalSMP(const double t0, const double Tc0,
STensorOperation::zero(PR);
STensorOperation::zero(PAM);
STensorOperation::zero(PAML2);
+ STensorOperation::zero(PRL2);
};
@@ -149,6 +150,7 @@ IPPhenomenologicalSMP::IPPhenomenologicalSMP(const IPPhenomenologicalSMP &source
PR =source.PR;
PAM =source.PAM;
PAML2 =source.PAML2;
+ PRL2 =source.PRL2;
}
@@ -205,6 +207,7 @@ IPPhenomenologicalSMP& IPPhenomenologicalSMP::operator=(const IPVariable &source
PR =src->PR;
PAM =src->PAM;
PAML2 =src->PAML2;
+ PRL2 =src->PRL2;
//energy
appliedEnergy =src->appliedEnergy;
@@ -314,6 +317,7 @@ void IPPhenomenologicalSMP::restart()
restartManager::restart(PR);
restartManager::restart(PAM);
restartManager::restart(PAML2);
+ restartManager::restart(PRL2);
//energy
restartManager::restart(appliedEnergy);
diff --git a/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.h b/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.h
index 512ca6e353e8650530a28cf3a55be980297295da..26aa3297bf4614d209b722949a3ee36aa24fcd98 100644
--- a/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.h
+++ b/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.h
@@ -19,7 +19,7 @@ class IPPhenomenologicalSMP : public IPVariableMechanics//: public IPLinearTherm
IPHyperViscoElastoPlastic *qAM;
STensor3 FthG, FfG, FthR, FfR, FthI, FthAD, FthAM, FfAM;
- STensor3 PG, PR, PAM, PAML2;
+ STensor3 PG, PR, PAM, PAML2, PRL2;
#if 0
STensor3 FthAI;
@@ -80,6 +80,8 @@ public:
virtual const STensor3 &getConstRefToPAM() const {return PAM;};
virtual STensor3 &getRefToPAML2() {return PAML2;};
virtual const STensor3 &getConstRefToPAML2() const {return PAML2;};
+ virtual STensor3 &getRefToPRL2() {return PRL2;};
+ virtual const STensor3 &getConstRefToPRL2() const {return PRL2;};
virtual double &getRefTodefoEnergy() {return _SMPEnergy;};
virtual double defoEnergy() const {return _SMPEnergy;}
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
index 4096527b46682f9bd4c21986b1e5a1a1fb5536ef..5ca1ad71c2aa90e291f8961952e5f5cd8cb2e063 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
@@ -40,6 +40,7 @@ void mlawHyperViscoElastic::setViscoElasticData(const int i, const double Ei, co
if (i> _N or i<1)
Msg::Error("This setting is invalid %d > %d",i,_N);
else{
+ double _nu = (3.*_K-2.*_mu)/2./(3.*_K+_mu);
double KK = Ei/(3.*(1.-2.*_nu));
double GG = Ei/(2.*(1.+_nu));
_Ki[i-1] = KK;
@@ -97,6 +98,34 @@ void mlawHyperViscoElastic::setViscoElasticData(const std::string filename){
fclose(file);
};
+void mlawHyperViscoElastic::evaluatePhiPCorrection(const IPHyperViscoElastic *q1, double trEe, STensor3 devEpr, double &p_vol, double &Dppr_vol_DCepr, STensor3 &K_dev_pr, STensor43 &DKpr_dev_DCepr) const
+{
+ p_vol= getUpdatedBulkModulus(q1)*getVolumeCorrection()*pow(exp(getXiVolumeCorrection()*trEe)-1.,getZetaVolumeCorrection())*trEe;
+ Dppr_vol_DCepr= 0.5*(getUpdatedBulkModulus(q1)*getVolumeCorrection()*pow(exp(getXiVolumeCorrection()*trEe)-1.,getZetaVolumeCorrection()-1.)*(getXiVolumeCorrection()*getZetaVolumeCorrection()*trEe*exp(getXiVolumeCorrection()*trEe)+exp(getXiVolumeCorrection()*trEe)-1.));
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ K_dev_pr(i,j)=2.*getUpdatedShearModulus(q1)*getDevCorrection()*pow(exp(getThetaDevCorrection()*devEpr.dotprod())-1.,getPiDevCorrection())*devEpr(i,j);
+ }
+ }
+
+
+ static STensor3 temp_DKpr_dev_DCepr;
+ STensorOperation::multSTensor3STensor43(devEpr,_Idev,temp_DKpr_dev_DCepr);
+ 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++){
+ DKpr_dev_DCepr(i,j,k,l) = 2.*getUpdatedShearModulus(q1)*getDevCorrection()*(2.*getThetaDevCorrection()*getPiDevCorrection()*exp(getThetaDevCorrection()*devEpr.dotprod())* pow(exp(getThetaDevCorrection()*devEpr.dotprod())-1.,getPiDevCorrection()-1.)*(temp_DKpr_dev_DCepr(i,j))*devEpr(k,l)+pow(exp(getThetaDevCorrection()*devEpr.dotprod())-1.,getPiDevCorrection())*_Idev(i,j,k,l));
+ }
+ }
+ }
+ }
+
+
+}
+
+
double mlawHyperViscoElastic::deformationEnergy(const IPHyperViscoElastic& q) const
{
double Psy = 0.;
@@ -105,7 +134,10 @@ double mlawHyperViscoElastic::deformationEnergy(const IPHyperViscoElastic& q) co
double trEe;
static STensor3 devEe;
STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
- Psy = _K*0.5*trEe*trEe+_mu*STensorOperation::doubledot(devEe,devEe);
+
+ Psy = getUpdatedBulkModulus(&q)*(0.5*trEe*trEe+_volCorrection*pow(exp(_xivolCorrection*trEe)-1.,_zetavolCorrection))+getUpdatedShearModulus(&q)*STensorOperation::doubledot(devEe,devEe);
+
+ //Psy = getUpdatedBulkModulus(&q)*(0.5*trEe*trEe+_volCorrection*(exp(trEe)-1.)*(exp(trEe)-1.))+getUpdatedShearModulus(&q)*STensorOperation::doubledot(devEe,devEe); //this is not correct
for (int i=0; i<_N; i++){
static STensor3 devEebranch;
Psy += _Ki[i]*0.5*(q._B[i])*(q._B[i])+_Gi[i]*STensorOperation::doubledot(q._A[i],q._A[i]);
@@ -230,8 +262,8 @@ void mlawHyperViscoElastic::dViscousEnergydF(const IPHyperViscoElastic& q0, cons
mlawHyperViscoElastic::mlawHyperViscoElastic(const int num,const double E,const double nu, const double rho,
const bool matrixbyPerturbation, const double pert):
- materialLaw(num,true), _E(E),_nu(nu),_rho(rho),_tangentByPerturbation(matrixbyPerturbation),_perturbationfactor(pert),
- _viscoMethod(Maxwell),_N(0.),_order(1),_Ki(0),_ki(0),_Gi(0),_gi(0), _I4(1.,1.), _I(1.)
+ materialLaw(num,true), _rho(rho),_tangentByPerturbation(matrixbyPerturbation),_perturbationfactor(pert),
+ _viscoMethod(Maxwell),_N(0.),_order(1),_Ki(0),_ki(0),_Gi(0),_gi(0), _I4(1.,1.), _I(1.), _volCorrection(0.)
{
_Idev = _I4;
STensor3 mIon3(-1./3);
@@ -239,7 +271,8 @@ mlawHyperViscoElastic::mlawHyperViscoElastic(const int num,const double E,const
tensprod(_I,mIon3, mIIon3);
_Idev += mIIon3;
- _lambda = (E*nu)/(1.+nu)/(1.-2.*nu);
+ double _lambda = (E*nu)/(1.+nu)/(1.-2.*nu);
+
_mu = 0.5*E/(1.+nu);
_K = E/3./(1.-2.*nu);
@@ -269,12 +302,12 @@ mlawHyperViscoElastic::mlawHyperViscoElastic(const int num,const double E,const
Cel(2,1,1,2) = _mu;
};
mlawHyperViscoElastic::mlawHyperViscoElastic(const mlawHyperViscoElastic& src): materialLaw(src),
- _E(src._E),_nu(src._nu),_rho(src._rho),
- _lambda(src._lambda),_mu(src._mu),_K(src._K),Cel(src.Cel),
+ _rho(src._rho),
+ _mu(src._mu),_K(src._K),Cel(src.Cel),
_order(src._order),
_tangentByPerturbation(src._tangentByPerturbation),_perturbationfactor(src._perturbationfactor),
_N(src._N),_viscoMethod(src._viscoMethod),_Ki(src._Ki),_ki(src._ki),_Gi(src._Gi),_gi(src._gi),
- _I4(src._I4), _I(src._I), _Idev(src._Idev)
+ _I4(src._I4), _I(src._I), _Idev(src._Idev), _volCorrection(src._volCorrection)
{
};
@@ -285,8 +318,8 @@ mlawHyperViscoElastic& mlawHyperViscoElastic::operator=(const materialLaw &sourc
const mlawHyperViscoElastic* src =dynamic_cast<const mlawHyperViscoElastic*>(&source);
if(src != NULL)
{
- _E = src->_E; _nu = src->_nu; _rho = src->_rho,
- _lambda = src->_lambda; _mu = src->_mu; _K = src->_K; Cel = src->Cel;
+ _rho = src->_rho,
+ _mu = src->_mu; _K = src->_K; Cel = src->Cel;
_order = src->_order;
_tangentByPerturbation = src->_tangentByPerturbation; _perturbationfactor = src->_perturbationfactor;
_N = src->_N; _viscoMethod = src->_viscoMethod; _Ki = src->_Ki; _ki = src->_ki; _Gi = src->_Gi; _gi = src->_gi;
@@ -296,6 +329,8 @@ mlawHyperViscoElastic& mlawHyperViscoElastic::operator=(const materialLaw &sourc
double mlawHyperViscoElastic::soundSpeed() const
{
+ double _E=9*_K*_mu/(3*_K+_mu);
+ double _nu = (3.*_K-2.*_mu)/2./(3.*_K+_mu);
double factornu = (1.-_nu)/((1.+_nu)*(1.-2.*_nu));
return sqrt(_E*factornu/_rho);
}
@@ -313,8 +348,8 @@ void mlawHyperViscoElastic::updateViscoElasticFlow(const IPHyperViscoElastic *q0
STensorOperation::decomposeDevTr(DE,devDE,trDE);
// maxwell
- Ge = _mu;
- Ke = _K;
+ Ge = getUpdatedShearModulus(q1);
+ Ke = getUpdatedBulkModulus(q1);
for (int i=0; i<_Gi.size(); i++){
double dtg = dt/_gi[i];
@@ -379,8 +414,8 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
double dt = this->getTimeStep();
if (_viscoMethod == Maxwell){
- Ge = _mu;
- Ke = _K;
+ Ge = getUpdatedShearModulus(q1);
+ Ke = getUpdatedBulkModulus(q1);
for (int i=0; i<_Gi.size(); i++){
double dtg = dt/_gi[i];
double expmdtg = exp(-dtg);
@@ -400,11 +435,20 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
Ke += _Ki[i]*ztak;
}
- static STensor3 devK;
- static double p;
- STensorOperation::decomposeDevTr(Ee,devK,p);
- devK *= (2.*_mu); // deviatoric part
- p *= _K; // pressure
+ static STensor3 devK, K_dev_pr;
+ double p, trEe;
+ STensorOperation::decomposeDevTr(Ee,devK,trEe);
+ evaluatePhiPCorrection(q1,trEe, devK, q1->getRefTo_p_vol(), q1->getRefTo_Dppr_vol_DCepr(), q1->getRefTo_K_dev_pr(), q1->getRefTo_DKpr_dev_DCepr());
+ devK *= (2.*getUpdatedShearModulus(q1)); // deviatoric part
+ p = trEe;
+ p *= getUpdatedBulkModulus(q1); // pressure
+
+
+
+
+ K_dev_pr = q1->getRefTo_K_dev_pr();
+ p+=q1->getRefTo_p_vol();
+
for (int i=0; i<_Gi.size(); i++){
devK.daxpy(q1->_A[i],2.*_Gi[i]);
@@ -412,14 +456,27 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
for (int i=0; i<_Ki.size(); i++){
p += q1->_B[i]*_Ki[i];
}
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ devK(i,j) = devK(i,j)+K_dev_pr(i,j);
+ }
+ }
q1->_kirchhoff = devK;
q1->_kirchhoff(0,0) += p;
q1->_kirchhoff(1,1) += p;
q1->_kirchhoff(2,2) += p;
+ //std::cout <<" VolCorrection = "<<_volCorrection<<" Xi VolCorrection = "<<_xivolCorrection<<" Zeta VolCorrection = "<<_zetavolCorrection <<std::endl;
+
+ q1->_kirchhoff_dev = K_dev_pr;
+ q1->_kirchhoff_vol(0,0) += q1->getRefTo_p_vol();
+ q1->_kirchhoff_vol(1,1) += q1->getRefTo_p_vol();
+ q1->_kirchhoff_vol(2,2) += q1->getRefTo_p_vol();
}
+
else if (_viscoMethod == KelvinVoight){
- double invGe = 1./_mu;
+ double invGe = 1./getUpdatedShearModulus(q1);
STensor3 D(0.);
for (int i=0; i<_Gi.size(); i++){
double dtg = dt/(_gi[i]);
@@ -479,14 +536,39 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
}
}
else{
- static STensor3 devK;
- static double p;
- STensorOperation::decomposeDevTr(Ee,devK,p);
- devK *= (2.*_mu); // deviatoric part
- p *= _K; // pressure
+ static STensor3 devK, K_dev_pr;
+ double p, trEe, p_vol;
+ STensorOperation::decomposeDevTr(Ee,devK,trEe);
+ evaluatePhiPCorrection(q1,trEe, devK, q1->getRefTo_p_vol(), q1->getRefTo_Dppr_vol_DCepr(), q1->getRefTo_K_dev_pr(), q1->getRefTo_DKpr_dev_DCepr());
+ devK *= (2.*getUpdatedShearModulus(q1)); // deviatoric part
+ p = trEe;
+
+ p *= getUpdatedBulkModulus(q1); // pressure
+
+
+
+
+ K_dev_pr = q1->getRefTo_K_dev_pr();
+
+
+ p+=q1->getRefTo_p_vol();
+
+
+ q1->_kirchhoff_dev = K_dev_pr;
+ q1->_kirchhoff_vol(0,0) += q1->getRefTo_p_vol();
+ q1->_kirchhoff_vol(1,1) += q1->getRefTo_p_vol();
+ q1->_kirchhoff_vol(2,2) += q1->getRefTo_p_vol();
+
+ Ke = getUpdatedBulkModulus(q1);
+ Ge = getUpdatedShearModulus(q1);
+
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ devK(i,j) = devK(i,j)+K_dev_pr(i,j);
- Ke = _K;
- Ge = _mu;
+ }
+ }
q1->_kirchhoff = devK;
q1->_kirchhoff(0,0) += p;
@@ -517,7 +599,6 @@ void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0,
static STensor63 ddlnCddC;
STensor3& E = q1->getRefToElasticStrain();
-
STensorOperation::multSTensor3FirstTranspose(F,F,C);
if (_order == 1){
STensorOperation::logSTensor3(C,_order,E,&dlnCdC); // as ddlogCddC = 0
@@ -526,7 +607,10 @@ void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0,
STensorOperation::logSTensor3(C,_order,E,&dlnCdC,&ddlnCddC);
}
E *= 0.5; // strain
-
+ static STensor3 devE;
+ double trE;
+ STensorOperation::decomposeDevTr(E,devE,trE);
+
double Ke, Ge;
viscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge);
@@ -542,7 +626,8 @@ void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0,
static STensor43 DsecondPKdC;
static STensor43 DcorKDE;
isotropicHookTensor(Ke,Ge,DcorKDE);
-
+ double Dppr_vol_DCepr=q1->getRefTo_Dppr_vol_DCepr();
+ static STensor43 DKpr_dev_DCepr=q1->getRefTo_DKpr_dev_DCepr();
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
@@ -555,7 +640,10 @@ void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0,
}
for (int r=0; r<3; r++){
for (int s=0; s<3; s++){
- DsecondPKdC(i,j,k,l) += 0.5*DcorKDE(p,q,r,s)*dlnCdC(p,q,i,j)*dlnCdC(r,s,k,l);
+ DsecondPKdC(i,j,k,l) += 0.5*DcorKDE(p,q,r,s)*dlnCdC(p,q,i,j)*dlnCdC(r,s,k,l)+0.5*DKpr_dev_DCepr(p,q,r,s)*dlnCdC(p,q,i,j)*dlnCdC(r,s,k,l);
+
+ if(p==q and r==s)
+ DsecondPKdC(i,j,k,l) += Dppr_vol_DCepr*dlnCdC(p,q,i,j)*dlnCdC(r,s,k,l);
}
}
}
@@ -726,6 +814,291 @@ void mlawPowerYieldHyper::constitutive(
};
+void mlawPowerYieldHyper::evaluatePhiPCorrection(const IPHyperViscoElastic *q1, double _b, double v, double trEepr, double Gamma, double ptildepr_Vol, double ptildepr,double Ke, double u, STensor3 devEepr, STensor3 K_dev_pr, STensor3 devPhipr, double Ge, double &ptilde_Vol, double &X_vol, double &dptilde_VoldGamma, STensor3 &dptilde_VoldCepr, STensor3 &K_dev, double &D_dev, STensor3 &dK_devdGamma, STensor43 &dK_devdCepr) const
+{
+
+ STensorOperation::zero(dptilde_VoldCepr);
+ ptilde_Vol = ptildepr_Vol;
+ X_vol=0.;
+
+ double Dptilde_Vol=ptilde_Vol-ptildepr_Vol;
+ double trN = 2.*_b*(ptildepr+Dptilde_Vol)/(v);
+ //static STensor3 DpprDCepr;
+ int ite = 0;
+ int maxite = 100; // maximal number of iters
+
+ double J=getUpdatedBulkModulus(q1)*getVolumeCorrection()*pow(exp(getXiVolumeCorrection()*trEepr-getXiVolumeCorrection()*Gamma*trN)-1.,getZetaVolumeCorrection())*(trEepr-Gamma*trN)-ptilde_Vol;
+
+
+ while (fabs(J) >_tol or ite <1){
+
+ double dJdptilde_Vol = -getXiVolumeCorrection()*getUpdatedBulkModulus(q1)*getVolumeCorrection()*getZetaVolumeCorrection()*Gamma*(2.*_b/v)*(trEepr-Gamma*trN)*exp(getXiVolumeCorrection()*trEepr-getXiVolumeCorrection()*Gamma*trN)*pow(exp(getXiVolumeCorrection()*trEepr-getXiVolumeCorrection()*Gamma*trN)-1.,getZetaVolumeCorrection()-1.)-Gamma*getUpdatedBulkModulus(q1)*getVolumeCorrection()*(2.*_b/v)*pow(exp(getXiVolumeCorrection()*trEepr-getXiVolumeCorrection()*Gamma*trN)-1.,getZetaVolumeCorrection())-1.;
+
+
+ double dptilde_Vol = -J/dJdptilde_Vol;
+
+ ptilde_Vol += dptilde_Vol;
+
+ Dptilde_Vol = ptilde_Vol-ptildepr_Vol;
+ trN = 2.*_b*(ptildepr+Dptilde_Vol)/(v);
+ J=getUpdatedBulkModulus(q1)*getVolumeCorrection()*pow(exp(getXiVolumeCorrection()*trEepr-getXiVolumeCorrection()*Gamma*trN)-1.,getZetaVolumeCorrection())*(trEepr-Gamma*trN)-ptilde_Vol;
+
+ ite++;
+
+
+ if (fabs(J) <_tol) break;
+ if(ite > maxite){
+ // Msg::Error("No convergence for phi_p_vol in mlawPowerYieldHyper nonAssociatedFlow Maxwell iter = %d, J = %e!!",ite,J);
+ break;
+ }
+
+ }
+ X_vol = exp(getXiVolumeCorrection()*(trEepr-2.*_b*Gamma/v*(ptildepr+Dptilde_Vol)));
+
+ //dptilde_VoldGamma=-((2.*_b*getXiVolumeCorrection()*getXiVolumeCorrection()*getZetaVolumeCorrection()*getUpdatedBulkModulus(q1)*getVolumeCorrection())/(v*v)*pow(X_vol-1.,getZetaVolumeCorrection()-2.)*(getZetaVolumeCorrection()*X_vol-1.)*X_vol*(ptildepr+Dptilde_Vol))/(1.+(2.*_b*getXiVolumeCorrection()*getXiVolumeCorrection()*getZetaVolumeCorrection()*Gamma*getUpdatedBulkModulus(q1)*getVolumeCorrection())/(v)*pow(X_vol-1.,getZetaVolumeCorrection()-2.)*(getZetaVolumeCorrection()*X_vol-1.)*X_vol);
+
+ dptilde_VoldGamma= -( (2.*_b)/(v*v)*getUpdatedBulkModulus(q1)*getVolumeCorrection()*(ptildepr+Dptilde_Vol)*(getZetaVolumeCorrection()*pow(X_vol-1.,getZetaVolumeCorrection()-1. )*(trEepr-2.*_b*Gamma/v*(ptildepr+Dptilde_Vol))*X_vol*getXiVolumeCorrection()+pow(X_vol-1.,getZetaVolumeCorrection())) )/(1.+2.*_b*Gamma/v*getUpdatedBulkModulus(q1)*getVolumeCorrection()*(getZetaVolumeCorrection()*pow(X_vol-1.,getZetaVolumeCorrection()-1.)*(trEepr-2.*_b*Gamma/v*(ptildepr+Dptilde_Vol))*X_vol*getXiVolumeCorrection()+pow(X_vol-1.,getZetaVolumeCorrection()) )) ;
+
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+
+ dptilde_VoldCepr(i,j) = ((getUpdatedBulkModulus(q1)*getVolumeCorrection()*(X_vol*getXiVolumeCorrection()*getZetaVolumeCorrection()*pow(X_vol-1.,getZetaVolumeCorrection()-1.)* (trEepr-2.*_b*Gamma/v*(ptildepr+Dptilde_Vol))+pow(X_vol-1.,getZetaVolumeCorrection())))*(_I(i,j)-(2.*_b*Gamma)/(v)*Ke*_I(i,j))) / (1.+(2.*_b*Gamma)/(v)*getUpdatedBulkModulus(q1)*getVolumeCorrection()*(X_vol*getXiVolumeCorrection()*getZetaVolumeCorrection()*pow(X_vol-1.,getZetaVolumeCorrection()-1.) *(trEepr-2.*_b*Gamma/v*(ptildepr+Dptilde_Vol))+pow(X_vol-1.,getZetaVolumeCorrection())));
+ }
+ }
+
+
+
+ //################################### DEVIATORIC TERMS
+ static STensor3 DK_dev, devE, J_dev;
+ STensorOperation::zero(dK_devdGamma);
+ STensorOperation::zero(dK_devdCepr);
+ K_dev = K_dev_pr;
+ D_dev=0.;
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ DK_dev(i,j)=K_dev(i,j)-K_dev_pr(i,j);
+ }
+ }
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ devE(i,j) = devEepr(i,j)-3.*Gamma*(devPhipr(i,j)+DK_dev(i,j))/(u);
+ }
+ }
+
+ double Temp_Jdev = devE.dotprod();
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ J_dev(i,j)=2.*getUpdatedShearModulus(q1)*getDevCorrection()*pow(exp(getThetaDevCorrection()*Temp_Jdev)-1.,getPiDevCorrection())*devE(i,j)-K_dev(i,j) ;
+ }
+ }
+
+ double J_dev_tol = 0.;
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ J_dev_tol+=abs(J_dev(i,j)) ;
+ }
+ }
+ //std::cout << "###############CHECK Loop 1 J_dev_tol = "<<J_dev_tol<<std::endl;
+
+ int ite_dev = 0;
+ int maxite_dev = 100; // maximal number of iters
+ static STensor43 ddevE, dJ_dev, dJ_dev_inv;
+ static STensor3 dTemp, NK_dev;
+
+ STensorOperation::zero(dJ_dev);
+ double _tol_dev = _tol; //1.0e-15;
+ while (fabs(J_dev_tol) >_tol_dev or ite_dev <1){
+
+ ddevE = _I4;
+ ddevE*=-3.*Gamma/u;
+ STensorOperation::multSTensor3STensor43( devE, _I4, dTemp);
+ dTemp*=-6.*Gamma/u;
+
+ 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++){
+ dJ_dev(i,j,k,l) = 2.*getUpdatedShearModulus(q1)*getDevCorrection()*getThetaDevCorrection()*getPiDevCorrection()*exp(getThetaDevCorrection()*Temp_Jdev)* pow(exp(getThetaDevCorrection()*Temp_Jdev)-1.,getPiDevCorrection()-1.)*dTemp(i,j)*devE(k,l)+2.*getUpdatedShearModulus(q1)*getDevCorrection()*pow(exp(getThetaDevCorrection()*Temp_Jdev)-1.,getPiDevCorrection())*ddevE(i,j,k,l)-_I4(i,j,k,l);
+ }
+ }
+ }
+ }
+
+ STensorOperation::inverseSTensor43(dJ_dev, dJ_dev_inv);
+ STensorOperation::zero(NK_dev);
+ 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++){
+ NK_dev(i,j)+= dJ_dev_inv(i,j,k,l)*J_dev(l,k);
+ }
+ }
+ }
+ }
+
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ K_dev(i,j) = K_dev(i,j)-NK_dev(i,j);
+ }
+ }
+
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ DK_dev(i,j)=K_dev(i,j)-K_dev_pr(i,j);
+ }
+ }
+
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ devE(i,j) = devEepr(i,j)-3.*Gamma*(devPhipr(i,j)+DK_dev(i,j))/(u);
+ }
+ }
+
+ Temp_Jdev = devE.dotprod();
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ J_dev(i,j)=2.*getUpdatedShearModulus(q1)*getDevCorrection()*pow(exp(getThetaDevCorrection()*Temp_Jdev)-1.,getPiDevCorrection())*devE(i,j)-K_dev(i,j) ;
+ }
+ }
+
+ J_dev_tol=0.;
+
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ J_dev_tol+=abs(J_dev(i,j));
+ }
+ }
+
+ //std::cout << "###############CHECK DESPUES LOOP 2 J_dev_tol = "<<J_dev_tol << " ite_dev = "<< ite_dev<<std::endl;
+
+ if (fabs(J_dev_tol)<_tol_dev and ite_dev >1) break;
+ ite_dev++;
+ if(ite_dev > maxite_dev){
+ Msg::Error("No convergence for J_dev in mlawPowerYieldHyper nonAssociatedFlow Maxwell iter = %d, J_dev = %e!!",ite_dev,J_dev);
+ break;
+ }
+
+ }
+
+
+
+ //Variables
+ double M = 0., V = 0., P = 0.;
+ static STensor3 R, O, T, M_temp;
+ //Final matrices
+ static STensor3 U;
+ static STensor43 S, S_inv;
+
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ M_temp(i,j) = 2.*devE(i,j);
+ T(i,j) = (-3.)/(u*u)*(devPhipr(i,j)+DK_dev(i,j));
+ }
+ }
+
+ D_dev = exp(getThetaDevCorrection()*Temp_Jdev);
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ M+=M_temp(i,j)*T(j,i);
+ R(i,j) = -6.*Gamma/u*devE(i,j);
+ O(i,j) = 2.*getUpdatedShearModulus(q1)*getDevCorrection()*getPiDevCorrection()*getThetaDevCorrection()*D_dev*pow(D_dev-1.,getPiDevCorrection()-1.)*devE(i,j);
+ }
+ }
+
+ V = 2.*getUpdatedShearModulus(q1)*getDevCorrection()*pow(D_dev-1.,getPiDevCorrection());
+ P=-3.*Gamma/u;
+
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ U(i,j)=M*O(i,j)+V*T(i,j);
+ for (int k=0; k<3; k++){
+ for (int l=0; l<3; l++){
+ S(i,j,k,l) = -R(l,k)*O(i,j)+(1.-V*P)*_I4(i,j,k,l);
+ }
+ }
+ }
+ }
+
+ STensorOperation::inverseSTensor43(S, S_inv);
+ STensorOperation::multSTensor3STensor43(U,S_inv,dK_devdGamma);
+
+
+ //########################################dK_devdCepr
+ static STensor43 G,H;
+ static STensor43 S_dev,U_dev, S_dev_inv;
+ STensorOperation::zero(H);
+ STensorOperation::zero(G);
+ STensorOperation::zero(S_dev);
+ STensorOperation::zero(U_dev);
+ STensorOperation::zero(S_dev_inv);
+
+ 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++){
+ G(i,j,k,l) = 2.*getUpdatedShearModulus(q1)*getDevCorrection()*getPiDevCorrection()*pow(D_dev-1.,getPiDevCorrection()-1.)*devE(i,j)*devE(k,l);
+ }
+ }
+ }
+ }
+
+ 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++){
+ for (int m=0; m<3; m++){
+ for (int n=0; n<3; n++){
+ H(i,j,k,l) += (1.-6.*Ge*Gamma/u)*G(i,j,m,n)*_Idev(n,m,k,l);
+ }
+ }
+ }
+ }
+ }
+ }
+
+
+ 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++){
+ S_dev(i,j,k,l) = (1.+6.*getUpdatedShearModulus(q1)*getDevCorrection()*pow(D_dev-1.,getPiDevCorrection())*Gamma/u)*_I4(i,j,k,l)+6.*Gamma*D_dev*getThetaDevCorrection()/u*G(i,j,k,l);
+ U_dev(i,j,k,l) = 2.*D_dev*getThetaDevCorrection()*H(i,j,k,l)+2.*getUpdatedShearModulus(q1)*getDevCorrection()*pow(D_dev-1.,getPiDevCorrection())*(1.-6.*Ge*Gamma/u)*_Idev(i,j,k,l);
+ }
+ }
+ }
+ }
+
+ STensorOperation::inverseSTensor43(S_dev, S_dev_inv);
+
+
+ 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++){
+ for (int m=0; m<3; m++){
+ for (int n=0; n<3; n++){
+ dK_devdCepr(i,j,k,l) += S_dev_inv(i,j,m,n)*U_dev(n,m,k,l);
+ }
+ }
+ }
+ }
+ }
+ }
+
+}
+
+
void mlawPowerYieldHyper::updateEqPlasticDeformation(IPHyperViscoElastoPlastic *q1, const IPHyperViscoElastoPlastic *q0,
const double& nup, const double& Dgamma) const{
q1->_epspbarre = q0->_epspbarre+ Dgamma;
@@ -733,7 +1106,6 @@ void mlawPowerYieldHyper::updateEqPlasticDeformation(IPHyperViscoElastoPlastic *
q1->_epspTraction = q0->_epspTraction+ Dgamma;
double k = 1./(sqrt(1.+2.*nup*nup));
q1->_epspShear = q0->_epspShear+ Dgamma/(k*sqrt(2.));
-
};
@@ -808,20 +1180,32 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
DlnDCe = DlnDCepr;
// update A, B
- double Ge, Ke;
+ double Ge, Ke, trEepr;
+ static STensor3 devEepr;
viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge);
-
- static STensor3 PhiPr;
+
+ STensorOperation::decomposeDevTr(Ee,devEepr,trEepr);
+
+ static STensor3 PhiPr, PhiPr_Vol, PhiPr_Dev;
PhiPr = q1->_kirchhoff;
PhiPr -= q1->_backsig;
-
- static STensor3 devPhipr,devPhi; // effective dev stress predictor
- static double ptildepr,ptilde;
+
+ PhiPr_Dev = q1->_kirchhoff_dev;
+ PhiPr_Vol = q1->_kirchhoff_vol;
+
+ static STensor3 devPhipr,devPhi, devPhipr_Vol, K_dev_pr, K_dev; // effective dev stress predictor
+ static double ptildepr,ptilde, ptildepr_Vol, ptilde_Vol, ptildepr_Dev;
STensorOperation::decomposeDevTr(PhiPr,devPhipr,ptildepr);
ptildepr/= 3.;
-
+ STensorOperation::decomposeDevTr(PhiPr_Vol,devPhipr_Vol,ptildepr_Vol);
+ ptildepr_Vol/= 3.;
+ STensorOperation::decomposeDevTr(PhiPr_Dev,K_dev_pr,ptildepr_Dev);
+ ptildepr_Dev/= 3.;
+
ptilde = ptildepr; // current effective pressure
+ ptilde_Vol = ptildepr_Vol; // current effective volumetric pressure
devPhi =devPhipr;
+ K_dev = K_dev_pr;
double PhiEqpr2 = 1.5*devPhipr.dotprod();
double PhiEqpr = sqrt(PhiEqpr2);
@@ -833,7 +1217,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
double Gamma = 0.; // // flow rule parameter
double PhiEq = PhiEqpr;
-
+
// hardening
this->hardening(q0,q1);
static fullVector<double> a(3), Da(3); // yield coefficients and derivatives in respect to plastic deformation
@@ -855,6 +1239,9 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
double f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
+
+
+
double DfDGamma = 0.;
double dfdDgamma = 0.;
double u = 1.;
@@ -862,10 +1249,26 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
double A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
-
double dDgammaDGamma = 0.;
double Dgamma = 0.; // eqplastic strain
-
+ double Dptilde_vol = 0.0;
+ double dptilde_VoldGamma= 0.0;
+ double Dptilde_volDGamma=0.0;
+ double PhiEq_cor = 0.;
+ static STensor3 DK_dev, devPhipr_Dev_cor, dK_devdGamma = q1->getRefTo_dK_devdGamma();
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ DK_dev(i,j)=K_dev(i,j)-K_dev_pr(i,j);
+ }
+ }
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ devPhipr_Dev_cor(i,j) = devPhipr(i,j)+DK_dev(i,j);
+ }
+ }
+ PhiEq_cor = sqrt(1.5*devPhipr_Dev_cor.dotprod());
+
if (q1->dissipationIsBlocked()){
q1->getRefToDissipationActive() = false;
}
@@ -886,20 +1289,35 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
if (_viscosity != NULL)
_viscosity->get(q1->_epspbarre,eta,Deta);
double etaOverDt = eta/this->getTimeStep();
- double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A);
+
+
+ dptilde_VoldGamma=q1->getRefTo_dptilde_VoldGamma();
+
+ double dAdGamma = 1./(2.*A)*(-(72.*Gt*PhiEq*PhiEq)/(u)-(16.*_b*_b*_b*Kt*(ptildepr+Dptilde_vol)*(ptildepr+Dptilde_vol)/(3*v*v*v))+(8.*_b*_b*(ptildepr+Dptilde_vol))/(3*v*v)*dptilde_VoldGamma);
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ dAdGamma+= 18./(u*u)*(devPhipr_Dev_cor(i,j)*dK_devdGamma(j,i));
+ }
+ }
+ //std::cout << "###############CHECK DESPUES LOOP 2 dAdGamma = "<<dAdGamma <<std::endl;
dDgammaDGamma = kk*(A+Gamma*dAdGamma);
this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
- dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*ptilde -Da(0);
+ dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*(ptildepr+Dptilde_vol)/v -Da(0); //OK
if (Gamma>0 and etaOverDt>0)
dfdDgamma -= _p*pow(etaOverDt,_p-1.)*Deta/this->getTimeStep()*pow(Gamma,_p);
- DfDGamma = dfdDgamma*dDgammaDGamma - (_n*a(2)*6.*Gt)*pow(PhiEq,_n)/u + a(1)*ptilde*2.*_b*Kt/v;
+ DfDGamma = dfdDgamma*dDgammaDGamma - (_n*a(2)*6.*Gt)*pow(PhiEq,_n)/u + a(1)*(ptildepr+Dptilde_vol)*2.*_b*Kt/(v*v)-(a(1)*dptilde_VoldGamma)/(v);
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ DfDGamma+= 3.*_n*a(2)/(2.*u*u)*pow(PhiEq,_n-2.) *(devPhipr_Dev_cor(i,j)*dK_devdGamma(j,i));
+ }
+ }
+
if (Gamma>0 and etaOverDt>0)
DfDGamma -=pow(etaOverDt,_p)*_p*pow(Gamma,_p-1.);
double dGamma = -f/DfDGamma;
-
if (Gamma + dGamma <=0.){
Gamma /= 2.;
}
@@ -910,9 +1328,34 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
u = 1.+6.*Gt*Gamma;
v = 1.+2.*_b*Kt*Gamma;
- PhiEq = PhiEqpr/u;
- ptilde = ptildepr/v;
- A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
+
+
+ evaluatePhiPCorrection(q1, _b, v, trEepr, Gamma, ptildepr_Vol, ptildepr, Ke, u, devEepr, K_dev_pr, devPhipr, Ge, q1->getRefTo_ptilde_Vol(), q1->getRefTo_X_vol(), q1->getRefTo_dptilde_VoldGamma(), q1->getRefTo_dptilde_VoldCepr(), q1->getRefTo_K_dev(), q1->getRefTo_D_dev(), q1->getRefTo_dK_devdGamma(), q1->getRefTo_dK_devdCepr());
+
+
+
+ ptilde_Vol = q1->getRefTo_ptilde_Vol();
+ Dptilde_vol = ptilde_Vol-ptildepr_Vol;
+
+ dK_devdGamma = q1->getRefTo_dK_devdGamma();
+ K_dev = q1->getRefTo_K_dev();
+
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ DK_dev(i,j)=K_dev(i,j)-K_dev_pr(i,j);
+ }
+ }
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ devPhipr_Dev_cor(i,j) = devPhipr(i,j)+DK_dev(i,j);
+ }
+ }
+ PhiEq_cor = sqrt(1.5*devPhipr_Dev_cor.dotprod());
+ PhiEq = PhiEq_cor/u;
+
+ //PhiEq = PhiEqpr/u;
+ A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*((ptildepr+Dptilde_vol)*(ptildepr+Dptilde_vol))/(v*v));
Dgamma = kk*Gamma*A;
//Msg::Error("it = %d, u=%e, v=%e, Dgamma=%e",ite,u,v,Dgamma);
@@ -920,16 +1363,16 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
hardening(q0,q1);
getYieldCoefficients(q1,a);
- //a.print("a update");
- f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
+ f = a(2)*pow(PhiEq,_n) - (a(1)*(ptildepr+Dptilde_vol)/(v)+a(0));
double viscoTerm = etaOverDt*Gamma;
if (Gamma>0 and etaOverDt>0) f-= pow(viscoTerm,_p);
ite++;
//if (ite> maxite-5)
//Msg::Error("it = %d, DfDGamma = %e error = %e dGamma = %e, Gamma = %e",ite,DfDGamma,f,dGamma,Gamma);
-
+
+ std::cout << "NR F: f = "<<f <<" Gamma = "<<Gamma<<" dfdDgamma = "<<dfdDgamma<< " dDgammaDGamma = "<< dDgammaDGamma << " dptilde_VoldGamma = " << dptilde_VoldGamma<< " Dptilde_vol"<< Dptilde_vol<< " dptilde_VoldGamma = "<< dptilde_VoldGamma << " dGamma = "<<dGamma<< " DfDGamma = "<< DfDGamma <<" Iteration = "<<ite<<" ## Vc = " << getVolumeCorrection()<<" Xi = " << getXiVolumeCorrection()<<" Zeta = " << getZetaVolumeCorrection()<<" Dc = " << getDevCorrection()<<" Theta = " << getThetaDevCorrection()<<" Pi = " << getPiDevCorrection()<<std::endl;
if (fabs(f) <_tol) break;
if(ite > maxite){
@@ -938,12 +1381,14 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
return;
}
};
-
+ std::cout << "###################################################################################################################################################################################################"<<std::endl;
q1->_DgammaDt = Dgamma/this->getTimeStep();
// update effective stress tensor
+
+ devPhi =devPhipr_Dev_cor;
devPhi *= (1./u);
- ptilde = ptildepr/v;
+ ptilde = (ptildepr+Dptilde_vol)/v;
// update normal
devN = devPhi;
@@ -998,7 +1443,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
const STensor3& KS = q1->_kirchhoff;
// second Piola Kirchhoff stress
static STensor3 S;
- STensorOperation::multSTensor3STensor43(KS,DlnDCe,S);
+ STensorOperation::multSTensor3STensor43(KS,DlnDCe,S);
for(int i=0; i<3; i++)
for(int j=0; j<3; j++){
@@ -1034,10 +1479,14 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
static STensor3 DpprDCepr;
static STensor43 DdevKprDCepr;
STensorOperation::multSTensor3STensor43(_I,DlnDCepr,DpprDCepr);
- DpprDCepr*= (0.5*Ke);
+ static STensor3 devEe;
+ double trEe;
+
+ STensorOperation::decomposeDevTr(Ee,devEe,trEe);
+
+ DpprDCepr*= (0.5*Ke)+ q1->getRefTo_Dppr_vol_DCepr();
STensorOperation::multSTensor43(_Idev,DlnDCepr,DdevKprDCepr);
DdevKprDCepr*= Ge;
-
static STensor3 DpDCepr;
static STensor43 DdevKDCepr;
DpDCepr = DpprDCepr;
@@ -1051,23 +1500,34 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
if (Gamma >0){
// plastic
- static STensor3 dAdCepr, dfDCepr;
-
- double fact = 1.5*a(2)*_n*pow(PhiEq,_n-2.)/(u*u);
+ static STensor3 dAdCepr, dfDCepr, dptilde_VoldCepr, dgamadCepr, DGDCepr, dptildeOldpr_VoldCepr;
+ static STensor43 dK_devdCepr;
+
+ DpDCepr *=0.;
+
+ dptilde_VoldCepr = q1->getRefTo_dptilde_VoldCepr();
+ dK_devdCepr = q1->getRefTo_dK_devdCepr();
+ STensorOperation::multSTensor43(0.5*dK_devdCepr, DlnDCepr, dK_devdCepr);
+ STensorOperation::multSTensor3STensor43(0.5*dptilde_VoldCepr,DlnDCepr,dptilde_VoldCepr);
+ STensorOperation::multSTensor3STensor43(_I,DlnDCepr,dptildeOldpr_VoldCepr);
+ dptildeOldpr_VoldCepr*=0.5*Ke;
+
+ double fact = 1.5*a(2)*_n*pow(PhiEq,_n-2.)/(u*u);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
- dAdCepr(i,j) = (4.*_b*_b*ptildepr/(A*3.*v*v))*DpprDCepr(i,j);
- dfDCepr(i,j) = -(a(1)/v)*DpprDCepr(i,j);
+ dAdCepr(i,j) = (4.*_b*_b*(ptildepr+Dptilde_vol)/(A*3.*v*v))*(dptildeOldpr_VoldCepr(i,j)+dptilde_VoldCepr(i,j));
+ dfDCepr(i,j) = -(a(1)/v)*(dptildeOldpr_VoldCepr(i,j)+dptilde_VoldCepr(i,j));
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
- dAdCepr(i,j) += (9./(A*u*u))*devPhipr(k,l)*DdevKprDCepr(k,l,i,j);
- dfDCepr(i,j) += fact*devPhipr(k,l)*DdevKprDCepr(k,l,i,j);
+ dAdCepr(i,j) += (9./(A*u*u))*devPhipr_Dev_cor(k,l)*(DdevKprDCepr(k,l,i,j)+dK_devdCepr(k,l,i,j));
+ dfDCepr(i,j) += fact*devPhipr_Dev_cor(k,l)*(DdevKprDCepr(k,l,i,j)+dK_devdCepr(k,l,i,j));
}
}
}
}
+
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
DGDCepr(i,j) = (-dfDCepr(i,j)-dfdDgamma*kk*Gamma*dAdCepr(i,j))/DfDGamma;
@@ -1083,18 +1543,16 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
- DtrNDCepr(i,j) = 2.*_b/v*DpprDCepr(i,j) - 2.*_b*ptildepr*(2.*_b*Kt)/(v*v)*DGDCepr(i,j);
+ DtrNDCepr(i,j) = 2.*_b/v*(dptildeOldpr_VoldCepr(i,j)+dptilde_VoldGamma*DGDCepr(i,j)+dptilde_VoldCepr(i,j)) - 2.*_b*ptilde*(2.*_b*Kt)/(v)*DGDCepr(i,j);
+
- DdevNDCepr = (DdevKprDCepr);
- DdevNDCepr *= (3./u);
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++){
- DdevNDCepr(i,j,k,l) -= 18.*Gt/(u*u)*devPhipr(i,j)*DGDCepr(k,l);
+ DdevNDCepr(i,j,k,l) = (3./u)*(DdevKprDCepr(i,j,k,l)+dK_devdCepr(i,j,k,l))+(3./u*dK_devdGamma(i,j)-18.*Gt/(u*u)*devPhipr_Dev_cor(i,j))*DGDCepr(k,l);
}
-
static STensor43 temp1;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
@@ -1113,14 +1571,16 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
}
}
+ DpDCepr=dptildeOldpr_VoldCepr;
+
STensorOperation::multSTensor43(EprFp0,temp1,dFpDCepr);
// update
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
- DpDCepr(i,j) -= Ke*(DGDCepr(i,j)*trN+Gamma*DtrNDCepr(i,j));
+ DpDCepr(i,j) -= Ke*(DGDCepr(i,j)*trN+Gamma*DtrNDCepr(i,j))-dptilde_VoldGamma*DGDCepr(i,j)-dptilde_VoldCepr(i,j);
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
- DdevKDCepr(i,j,k,l) -= 2.*Ge*(DGDCepr(k,l)*devN(i,j)+Gamma*DdevNDCepr(i,j,k,l));
+ DdevKDCepr(i,j,k,l) -= 2.*Ge*(DGDCepr(k,l)*devN(i,j)+Gamma*DdevNDCepr(i,j,k,l))-dK_devdCepr(i,j,k,l)-DGDCepr(k,l)*dK_devdGamma(i,j);
}
}
}
@@ -1250,7 +1710,6 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
dViscousEnergydF(*q0, *q1, DlnDCe, Fe, q1->getRefToDViscousEnergyPartdF(), dFedF);
dPlasticEnergydF(N, Gamma, Dgamma, dKcorDcepr, DGDCepr, DdevNDCepr, DtrNDCepr, KS, CeprToF,q1->getRefToDPlasticEnergyPartdF());
-
STensor3& DirrEnergDF = q1->_DirreversibleEnergyDF;
if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
@@ -1283,8 +1742,6 @@ void mlawPowerYieldHyper::dDeformationEnergydF(const IPHyperViscoElastoPlastic&
}
}
-
-
void mlawPowerYieldHyper::dViscousEnergydF(const IPHyperViscoElastoPlastic& q0, const IPHyperViscoElastoPlastic& q, const STensor43 &dlnCdC, const STensor3 &Fe, STensor3 &dViscousEnergydF, const STensor43 &dFedF) const
{
STensorOperation::zero(dViscousEnergydF);
@@ -1865,7 +2322,6 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c
void mlawPowerYieldHyper::predictorCorrector(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1,
STensor3&P, const bool stiff, STensor43& Tangent,STensor43& dFedF, STensor43& dFpdF,
STensor43* elasticTangent) const{
-
if (_tangentByPerturbation){
if (_nonAssociatedFlow){
this->predictorCorrector_nonAssociatedFlow(F,q0,q1,P,false,Tangent,dFedF,dFpdF);
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.h b/NonLinearSolver/materialLaw/mlawHyperelastic.h
index 89c1dd540423296e39289f82c2ed5e38e6fdef70..de3347fe7ecc34f6911d97cc4fa2869786d5ae42 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.h
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.h
@@ -27,11 +27,15 @@ class mlawHyperViscoElastic : public materialLaw{
STensor3 _I; // second order tensor
STensor43 _Idev; // dev part of _I4
double _rho; // density
- double _E; // young modulus
- double _nu; // Poisson ratio
- double _lambda; // 1st lame parameter
double _mu; // 2nd lame parameter (=G)
double _K; // bulk modulus
+ double _volCorrection; //Vc correction on volume potential
+ double _xivolCorrection; //Xi correction on volume potential
+ double _zetavolCorrection; //Zeta correction on volume potential
+
+ double _devCorrection; //Dc correction on volume potential
+ double _thetadevCorrection; //Theta correction on volume potential
+ double _pidevCorrection; //pi correction on volume potential
STensor43 Cel; // elastic hook tensor
int _order ; // to approximated log and exp of a tensor
@@ -62,6 +66,7 @@ class mlawHyperViscoElastic : public materialLaw{
void predictorCorrector_ViscoElastic(const STensor3& F0, const STensor3& F, STensor3&P, const IPHyperViscoElastic *q0_, IPHyperViscoElastic *q1,
const bool stiff, STensor43& Tangent) const;
+ void evaluatePhiPCorrection(const IPHyperViscoElastic *q1, double trEe, STensor3 devEpr, double &p_vol, double &Dppr_vol_DCepr, STensor3 &K_dev_pr, STensor43 &DKpr_dev_DCepr) const;
#endif // SWIG
public:
@@ -75,6 +80,14 @@ class mlawHyperViscoElastic : public materialLaw{
virtual void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
virtual void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
virtual void setViscoElasticData(const std::string filename);
+ virtual void setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc) {_volCorrection=_vc,_xivolCorrection=_xivc,_zetavolCorrection=_zetavc, _devCorrection=_dc,_thetadevCorrection=_thetadc,_pidevCorrection=_pidc;};
+ virtual double getVolumeCorrection() const {if (_volCorrection<1.e-5) return 0.; else return _volCorrection;};
+ virtual double getXiVolumeCorrection() const {if (_xivolCorrection<1.e-5) return 0.5; else return _xivolCorrection;};
+ virtual double getZetaVolumeCorrection() const {if (_zetavolCorrection<1.e-5) return 2.; else return _zetavolCorrection;};
+
+ virtual double getDevCorrection() const {if (_devCorrection<1.e-5) return 0.; else return _devCorrection;};
+ virtual double getThetaDevCorrection() const {if (_thetadevCorrection<1.e-5) return 0.5; else return _thetadevCorrection;};
+ virtual double getPiDevCorrection() const {if (_pidevCorrection<1.e-5) return 2.; else return _pidevCorrection;};
#ifndef SWIG
mlawHyperViscoElastic(const mlawHyperViscoElastic& src);
@@ -94,7 +107,10 @@ class mlawHyperViscoElastic : public materialLaw{
virtual double density()const{return _rho;}
virtual const double bulkModulus() const{return _K;}
virtual const double shearModulus() const{return _mu;}
- virtual const double poissonRatio() const{return _nu;}
+ virtual const double poissonRatio() const{return (3.*_K-2.*_mu)/2./(3.*_K+_mu);}
+
+ virtual double getUpdatedBulkModulus(const IPHyperViscoElastic *q1) const { return _K*q1->getConstRefToElasticBulkPropertyScaleFactor(); }
+ virtual double getUpdatedShearModulus(const IPHyperViscoElastic *q1) const { return _mu*q1->getConstRefToElasticShearPropertyScaleFactor(); }
virtual double scaleFactor() const { return _mu;};
@@ -161,6 +177,8 @@ class mlawPowerYieldHyper : public mlawHyperViscoElastic{
virtual void predictorCorrector_associatedFlow(const STensor3& F, const IPHyperViscoElastoPlastic *q0_, IPHyperViscoElastoPlastic *q1,
STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF) const;
+
+ void evaluatePhiPCorrection(const IPHyperViscoElastic *q1, double _b, double v, double trEepr, double Gamma, double ptildepr_Vol, double ptildepr,double Ke, double u, STensor3 devEepr, STensor3 K_dev_pr, STensor3 devPhipr, double Ge, double &ptilde_Vol, double &X_vol, double &dptilde_VoldGamma, STensor3 &dptilde_VoldCepr, STensor3 &K_dev, double &D_dev, STensor3 &dK_devdGamma, STensor43 &dK_devdCepr) const;
protected:
@@ -176,6 +194,8 @@ class mlawPowerYieldHyper : public mlawHyperViscoElastic{
const IPHyperViscoElastoPlastic* q0, // previous internal variables
IPHyperViscoElastoPlastic* q1 // current internal variables
) const;
+
+
#endif // SWIG
public:
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
index 480e9295e495a01d55ea9e034bf5cfe9abf41003..57494458544e8c84d888b07a55c673d728193ec2 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
@@ -226,6 +226,7 @@ void mlawNonLinearTVE::setViscoElasticData(const int i, const double Ei, const d
if (i> _N or i<1)
Msg::Error("This setting is invalid %d > %d",i,_N);
else{
+ double _nu = (3.*_K-2.*_mu)/2./(3.*_K+_mu);
double KK = Ei/(3.*(1.-2.*_nu));
double GG = Ei/(2.*(1.+_nu));
@@ -247,6 +248,7 @@ void mlawNonLinearTVE::setViscoElasticData(const int i, const double Ei, const d
if (i> _N or i<1)
Msg::Error("This setting is invalid %d > %d",i,_N);
else{
+ double _nu = (3.*_K-2.*_mu)/2./(3.*_K+_mu);
double KK = Ei/(3.*(1.-2.*_nu));
double GG = Ei/(2.*(1.+_nu));
double AA = Alphai;
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
index 8a8198b4f641f8765cb3671542dde63d7ba29635..bfcf853ea2f78a462f3cf9e95f262ab462d914ab 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
@@ -293,6 +293,7 @@ void mlawNonLinearTVM::setViscoElasticData(const int i, const double Ei, const d
if (i> _N or i<1)
Msg::Error("This setting is invalid %d > %d",i,_N);
else{
+ double _nu = (3.*_K-2.*_mu)/2./(3.*_K+_mu);
double KK = Ei/(3.*(1.-2.*_nu));
double GG = Ei/(2.*(1.+_nu));
@@ -314,6 +315,7 @@ void mlawNonLinearTVM::setViscoElasticData(const int i, const double Ei, const d
if (i> _N or i<1)
Msg::Error("This setting is invalid %d > %d",i,_N);
else{
+ double _nu = (3.*_K-2.*_mu)/2./(3.*_K+_mu);
double KK = Ei/(3.*(1.-2.*_nu));
double GG = Ei/(2.*(1.+_nu));
double AA = Alphai;
@@ -2202,4 +2204,4 @@ void mlawNonLinearTVM::constitutive(const STensor3& F0,
// q1->_elasticEnergy=deformationEnergy(Fn,gradT); //comment it to get H1 norm in the interface e.g. LinearThermoMech
};
-*/
\ No newline at end of file
+*/
diff --git a/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.cpp b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.cpp
index 66410d06a8acd8519e3b0cd76752d129113939c5..df97bdf70f9810799b68c586b3b17fd52d073b1b 100644
--- a/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.cpp
+++ b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.cpp
@@ -50,8 +50,14 @@ mlawPhenomenologicalSMP::mlawPhenomenologicalSMP(const int num, const double rho
_temFunc_kqAM = new constantScalarFunction(1.);
_temFunc_alphaAM= new constantScalarFunction(1.);
_temFunc_cAM = new constantScalarFunction(1.);
+
+ _yieldCrystallinityG = new constantScalarFunction(1.);
+ _yieldCrystallinityR = new constantScalarFunction(1.);
+ _yieldCrystallinityAM = new constantScalarFunction(1.);
+
- _EPFunc_AM = new NeoHookeanElasticPotential(num, EAM, nuAM);
+ _EPFunc_AM = NULL;
+ _EPFunc_R = NULL;
for (int i=0; i<3; i++){
@@ -226,6 +232,21 @@ mlawPhenomenologicalSMP::mlawPhenomenologicalSMP(const mlawPhenomenologicalSMP &
if (source._temFunc_alphaAM!=NULL)
{_temFunc_alphaAM = source._temFunc_alphaAM->clone();}
+ _yieldCrystallinityG = NULL;
+ if (source._yieldCrystallinityG!=NULL)
+ {_yieldCrystallinityG=source._yieldCrystallinityG->clone();}
+
+ _yieldCrystallinityR = NULL;
+ if (source._yieldCrystallinityR!=NULL)
+ {_yieldCrystallinityR=source._yieldCrystallinityR->clone();}
+
+ _yieldCrystallinityAM = NULL;
+ if (source._yieldCrystallinityAM!=NULL)
+ {_yieldCrystallinityAM=source._yieldCrystallinityAM->clone();}
+
+
+
+
lawG = NULL;
if (source.lawG != NULL)
lawG = dynamic_cast<mlawPowerYieldHyper*>(source.lawG->clone());
@@ -240,6 +261,11 @@ mlawPhenomenologicalSMP::mlawPhenomenologicalSMP(const mlawPhenomenologicalSMP &
if (source._EPFunc_AM!=NULL)
{_EPFunc_AM = source._EPFunc_AM->clone();}
+ _EPFunc_R = NULL;
+ if (source._EPFunc_R!=NULL)
+ {_EPFunc_R = source._EPFunc_R->clone();}
+
+
}
@@ -341,10 +367,23 @@ mlawPhenomenologicalSMP& mlawPhenomenologicalSMP::operator=(const materialLaw &
if(_temFunc_cAM!=NULL) delete _temFunc_cAM; _temFunc_cAM=NULL;
if(src->_temFunc_cAM!=NULL){_temFunc_cAM=src->_temFunc_cAM->clone();}
+ if(_yieldCrystallinityG != NULL) delete _yieldCrystallinityG; _yieldCrystallinityG=NULL;
+ if (src->_yieldCrystallinityG!=NULL) {_yieldCrystallinityG=src->_yieldCrystallinityG->clone();}
+
+ if(_yieldCrystallinityR != NULL) delete _yieldCrystallinityR; _yieldCrystallinityR=NULL;
+ if (src->_yieldCrystallinityR!=NULL) {_yieldCrystallinityR=src->_yieldCrystallinityR->clone();}
+
+ if(_yieldCrystallinityAM != NULL) delete _yieldCrystallinityAM; _yieldCrystallinityAM=NULL;
+ if (src->_yieldCrystallinityAM!=NULL) {_yieldCrystallinityAM=src->_yieldCrystallinityAM->clone();}
+
if(_EPFunc_AM!=NULL) delete _EPFunc_AM; _EPFunc_AM=NULL;
if(src->_EPFunc_AM!=NULL){_EPFunc_AM=src->_EPFunc_AM->clone();}
+ if(_EPFunc_R!=NULL) delete _EPFunc_R; _EPFunc_R=NULL;
+ if(src->_EPFunc_R!=NULL){_EPFunc_R=src->_EPFunc_R->clone();}
+
}
+
return *this;
}
@@ -595,26 +634,11 @@ void mlawPhenomenologicalSMP::PhaseTransitionThermalParameters
////////
}
-void mlawPhenomenologicalSMP::getZgFunction_Sy0G(const double zg, double &wp, double &dwpdz) const
+void mlawPhenomenologicalSMP::getZgFunction_Sy0(const scalarFunction& Tfunc, const double zg, double &wp, double &dwpdz) const
{
- double tol=0.9;
- if(zg<tol)
- {
- wp=1.e-6;
- dwpdz=0.;
-
- }
- else
- {
- double n=1.;
- wp=pow((zg-tol)/(1.-tol),n)+1.e-6;
- dwpdz=pow((zg-tol)/(1.-tol),n-1.)/(1.-tol);
- if(wp>1.)
- {
- wp=1.;
- dwpdz=0.;
- }
- }
+ wp = 1.-Tfunc.getVal(zg);
+ dwpdz = -Tfunc.getDiff(zg);
+
}
@@ -695,7 +719,19 @@ void mlawPhenomenologicalSMP::funFthAD
}
}
+void mlawPhenomenologicalSMP::funEnThAI(const double T, const double TRef, const IPHyperViscoElastic& q0, IPHyperViscoElastic &q1,
+ double anisotropicThermalCoefAmplitudeBulk, double anisotropicThermalCoefAmplitudeShear,
+ double anisotropicThermalAlphaCoefTanhTempBulk, double anisotropicThermalAlphaCoefTanhTempShear) const
+{
+ double alphaTK=1.+anisotropicThermalCoefAmplitudeBulk*tanh(anisotropicThermalAlphaCoefTanhTempBulk*(T-TRef));
+ q1.getRefToElasticBulkPropertyScaleFactor()=alphaTK;
+ double alphaTG=1.+anisotropicThermalCoefAmplitudeShear*tanh(anisotropicThermalAlphaCoefTanhTempShear*(T-TRef));
+ q1.getRefToElasticShearPropertyScaleFactor()=alphaTG;
+
+}
+
+#if 0
void mlawPhenomenologicalSMP::funEnThAI
(const double T, const double TRef,
const STensor3& Ene, STensor3& strainEffect, STensor3& dStrainEffectdT, STensor43& dStrainEffectdEne,
@@ -888,7 +924,7 @@ void mlawPhenomenologicalSMP::funEnThAI
}
}
}
-
+#endif
void mlawPhenomenologicalSMP::funFthI
(const IPPhenomenologicalSMP *q0, IPPhenomenologicalSMP *q1,
@@ -1068,23 +1104,28 @@ void mlawPhenomenologicalSMP::constitutive(
STensor3& PG =q1->getRefToPG();
STensor3& PR =q1->getRefToPR();
STensor3& PAM =q1->getRefToPAM();
+ STensor3& PRL2 =q1->getRefToPRL2();
STensor3& PAML2 =q1->getRefToPAML2();
STensorOperation::zero(PG);
STensorOperation::zero(PR);
STensorOperation::zero(PAM);
+ STensorOperation::zero(PRL2);
STensorOperation::zero(PAML2);
- static STensor3 dPGdT, dPRdT, dPAMdT, dPAML2dT;
+ static STensor3 dPGdT, dPRdT, dPAMdT, dPAML2dT, dPRL2dT;
STensorOperation::zero(dPGdT);
STensorOperation::zero(dPRdT);
STensorOperation::zero(dPAMdT);
+ STensorOperation::zero(dPRL2dT);
STensorOperation::zero(dPAML2dT);
- static STensor43 TangentG, TangentR, TangentAM, TangentAML2;
+ static STensor43 TangentG, TangentR, TangentAM, TangentAML2, TangentRL2;
STensorOperation::zero(TangentG);
STensorOperation::zero(TangentR);
STensorOperation::zero(TangentAM);
STensorOperation::zero(TangentAML2);
+ STensorOperation::zero(TangentRL2);
+
double &zg =q1->getRefTozg();
double dzgdT=0.;
@@ -1145,6 +1186,7 @@ void mlawPhenomenologicalSMP::constitutive(
double PsiG=0.;
double PsiR=0.;
double PsiAM=0.;
+ double PsiRL2=0.;
double PsiAML2=0.;
double VFG = (1-_zAM)*(zg);
@@ -1202,6 +1244,7 @@ void mlawPhenomenologicalSMP::constitutive(
PR, dPRdT, TangentR,
PsiAM, dPlasAM, dPlasAMdT, dPlasAMdF, dViscAM, dViscAMdT, dViscAMdF,
PAM, dPAMdT, TangentAM,
+ PsiRL2, PRL2, dPRL2dT, TangentRL2,
PsiAML2, PAML2, dPAML2dT, TangentAML2);}
if ((zg>toleranceOnZg)and(zg<1.-toleranceOnZg))
@@ -1216,6 +1259,7 @@ void mlawPhenomenologicalSMP::constitutive(
PR, dPRdT, TangentR,
PsiAM, dPlasAM, dPlasAMdT, dPlasAMdF, dViscAM, dViscAMdT, dViscAMdF,
PAM, dPAMdT, TangentAM,
+ PsiRL2, PRL2, dPRL2dT, TangentRL2,
PsiAML2, PAML2, dPAML2dT, TangentAML2);}
if (zg>=1.-toleranceOnZg)
@@ -1230,10 +1274,20 @@ void mlawPhenomenologicalSMP::constitutive(
PR, dPRdT, TangentR,
PsiAM, dPlasAM, dPlasAMdT, dPlasAMdF, dViscAM, dViscAMdT, dViscAMdF,
PAM, dPAMdT, TangentAM,
+ PsiRL2, PRL2, dPRL2dT, TangentRL2,
PsiAML2, PAML2, dPAML2dT, TangentAML2);}
- P =VFG*PG+ VFR*PR + VFAM*PAM + VFAM*PAML2;
+ double VFAM2=VFAM;
+ double dVFAM2dT=dVFAMdT;
+ STensor3 dVFAM2dF(dVFAMdF);
+
+ double VFR2=VFR;
+ double dVFR2dT=dVFRdT;
+ STensor3 dVFR2dF(dVFRdF);
+
+
+ P =VFG*PG+ VFR*PR + VFAM*PAM + VFR2*PRL2 + VFAM2*PAML2;
//Wth & mecSource
dPlas = VFG*dPlasG + VFR*dPlasR + VFAM*dPlasAM;
@@ -1244,7 +1298,7 @@ void mlawPhenomenologicalSMP::constitutive(
dViscdT = dVFGdT*dViscG + VFG*dViscGdT + dVFRdT*dViscR + VFR*dViscRdT + dVFAMdT*dViscAM + VFAM*dViscAMdT;
dViscdF = dVFGdF*dViscG + VFG*dViscGdF + dVFRdF*dViscR + VFR*dViscRdF + dVFAMdF*dViscAM + VFAM*dViscAMdF;
- Psi = VFG*PsiG + VFR*PsiR + VFAM*(PsiAM + PsiAML2);
+ Psi = VFG*PsiG + VFR*PsiR + VFAM*PsiAM + VFAM2*PsiAML2 + VFR2*PsiRL2;
STensor3 Fninv, F0inv;
@@ -1257,17 +1311,20 @@ void mlawPhenomenologicalSMP::constitutive(
const STensor3& PR0 =q0->getConstRefToPR();
const STensor3& PAM0 =q0->getConstRefToPAM();
const STensor3& PAML20 =q0->getConstRefToPAML2();
+ const STensor3& PRL20 =q0->getConstRefToPRL2();
- STensor3 PmG, PmR, PmAM, PmAML2;
+ STensor3 PmG, PmR, PmAM, PmRL2, PmAML2;
STensorOperation::multSTensor3SecondTranspose(PG,FthG,PmG);
STensorOperation::multSTensor3SecondTranspose(PR,FthR,PmR);
STensorOperation::multSTensor3SecondTranspose(PAM,FthAM,PmAM);
STensorOperation::multSTensor3SecondTranspose(PAML2,FthAM,PmAML2);
- STensor3 PmG0, PmR0, PmAM0, PmAML20;
+ STensorOperation::multSTensor3SecondTranspose(PRL2,FthR,PmRL2);
+ STensor3 PmG0, PmR0, PmAM0, PmAML20, PmRL20;
STensorOperation::multSTensor3SecondTranspose(PG0,FthG0,PmG0);
STensorOperation::multSTensor3SecondTranspose(PR0,FthG0,PmR0);
STensorOperation::multSTensor3SecondTranspose(PAM0,FthAM0,PmAM0);
STensorOperation::multSTensor3SecondTranspose(PAML20,FthAM0,PmAML20);
+ STensorOperation::multSTensor3SecondTranspose(PRL20,FthR0,PmRL20);
double dAppl=0.;
@@ -1276,7 +1333,8 @@ void mlawPhenomenologicalSMP::constitutive(
dAppl+=VFG*(PmG(i,j) + PmG0(i,j))/2.*(FmG(i,j) - FmG0(i,j))
+VFR*(PmR(i,j) + PmR0(i,j))/2.*(FmR(i,j) - FmR0(i,j))
+VFAM*(PmAM(i,j) + PmAM0(i,j))/2.*(FmAM(i,j) - FmAM0(i,j))
- +VFAM*(PmAML2(i,j) + PmAML20(i,j))/2.*(FmAM(i,j) - FmAM0(i,j));
+ +VFAM2*(PmAML2(i,j) + PmAML20(i,j))/2.*(FmAM(i,j) - FmAM0(i,j));
+ +VFR2*(PmRL2(i,j) + PmRL20(i,j))/2.*(FmR(i,j) - FmR0(i,j));
}
}
@@ -1361,13 +1419,14 @@ void mlawPhenomenologicalSMP::constitutive(
STensorOperation::zero(dkqdT);
if(stiff)
{
- dPdT +=dVFGdT*PG + VFG*dPGdT + dVFRdT*PR + VFR*dPRdT + dVFAMdT*(PAM + PAML2) + VFAM*(dPAMdT + dPAML2dT);
+ dPdT +=dVFGdT*PG + VFG*dPGdT + dVFRdT*PR + VFR*dPRdT + dVFAMdT*PAM + dVFAM2dT*PAML2 + dVFR2dT*PRL2 + VFAM*dPAMdT + VFAM2*dPAML2dT + VFR2*dPRL2dT;
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++){
Tangent(i,j,k,l)+=dVFGdF(k,l)*PG(i,j) + VFG*TangentG(i,j,k,l) + dVFRdF(k,l)*PR(i,j) + VFR*TangentR(i,j,k,l)
- +dVFAMdF(k,l)*(PAM(i,j) + PAML2(i,j)) + VFAM*(TangentAM(i,j,k,l) + TangentAML2(i,j,k,l));
+ +dVFAMdF(k,l)*PAM(i,j) + dVFAM2dF(k,l)*PAML2(i,j)+ dVFR2dF(k,l)*PRL2(i,j) + VFAM*TangentAM(i,j,k,l) + VFAM2*TangentAML2(i,j,k,l)
+ + VFR2*TangentRL2(i,j,k,l);
}
}
}
@@ -1661,10 +1720,25 @@ void mlawPhenomenologicalSMP::constitutive1
STensor3& PR, STensor3& dPRdT, STensor43& TangentR,
double& PsiAM, double &dPlasAM, double& dPlasAMdT, STensor3& dPlasAMdF, double &dViscAM, double& dViscAMdT, STensor3& dViscAMdF,
STensor3 &PAM, STensor3 &dPAMdT, STensor43 &TangentAM,
+ double &PsiRL2, STensor3 &PRL2, STensor3 &dPRL2dT, STensor43 &TangentRL2,
double &PsiAML2, STensor3 &PAML2, STensor3 &dPAML2dT, STensor43 &TangentAML2) const
{
-
+ double anisotropicThermalCoefAmplitudeBulkG=(_alphaParam)[0];
+ double anisotropicThermalCoefAmplitudeShearG=(_alphaParam)[4];
+ double anisotropicThermalAlphaCoefTanhTempBulkG=(_alphaParam)[8];
+ double anisotropicThermalAlphaCoefTanhTempShearG=(_alphaParam)[12];
+
+ double anisotropicThermalCoefAmplitudeBulkR=(_alphaParam)[1];
+ double anisotropicThermalCoefAmplitudeShearR=(_alphaParam)[5];
+ double anisotropicThermalAlphaCoefTanhTempBulkR=(_alphaParam)[9];
+ double anisotropicThermalAlphaCoefTanhTempShearR=(_alphaParam)[13];
+
+ double anisotropicThermalCoefAmplitudeBulkAM=(_alphaParam)[3];
+ double anisotropicThermalCoefAmplitudeShearAM=(_alphaParam)[7];
+ double anisotropicThermalAlphaCoefTanhTempBulkAM=(_alphaParam)[11];
+ double anisotropicThermalAlphaCoefTanhTempShearAM=(_alphaParam)[15];
+
// Initialisation of internal variables and material properties
double &Tdot =q1->getRefToTdot();
@@ -1706,6 +1780,15 @@ void mlawPhenomenologicalSMP::constitutive1
evaluateFParti(Fn, FfR, FthR, FRn);
evaluateFParti(F0, FfR0, FthR0, FR0);
+
+ double wp=0., dwpdz=0.;
+ getZgFunction_Sy0(*_yieldCrystallinityR,q0->getzg(), wp, dwpdz);
+ q1->getRefToIPHyperViscoElastoPlasticR()._ipCompression->setScaleFactor(wp,0.);
+ q1->getRefToIPHyperViscoElastoPlasticR()._ipTraction->setScaleFactor(wp,0.);
+ funEnThAI(T0, q1->getConstRefToTRefR(), q0->getConstRefToIPHyperViscoElastoPlasticR(),q1->getRefToIPHyperViscoElastoPlasticR(),
+ anisotropicThermalCoefAmplitudeBulkR, anisotropicThermalCoefAmplitudeShearR,
+ anisotropicThermalAlphaCoefTanhTempBulkR, anisotropicThermalAlphaCoefTanhTempShearR);
+
lawR->constitutive(FR0,FRn, PVEVPR, &(q0->getConstRefToIPHyperViscoElastoPlasticR()), &(q1->getRefToIPHyperViscoElastoPlasticR()),dPVEVPRdFR,stiff);
@@ -1744,6 +1827,13 @@ void mlawPhenomenologicalSMP::constitutive1
evaluateFParti(Fn, FfAM, FthAM, FAMn);
evaluateFParti(F0, FfAM0, FthAM0, FAM0);
+ getZgFunction_Sy0(*_yieldCrystallinityAM,q0->getzg(), wp, dwpdz);
+ q1->getRefToIPHyperViscoElastoPlasticAM()._ipCompression->setScaleFactor(wp,0.);
+ q1->getRefToIPHyperViscoElastoPlasticAM()._ipTraction->setScaleFactor(wp,0.);
+ funEnThAI(T0, q1->getConstRefToTRefAM(), q0->getConstRefToIPHyperViscoElastoPlasticAM(),q1->getRefToIPHyperViscoElastoPlasticAM(),
+ anisotropicThermalCoefAmplitudeBulkAM, anisotropicThermalCoefAmplitudeShearAM,
+ anisotropicThermalAlphaCoefTanhTempBulkAM, anisotropicThermalAlphaCoefTanhTempShearAM);
+
lawAM->constitutive(FAM0,FAMn, PVEVPAM, &(q0->getConstRefToIPHyperViscoElastoPlasticAM()), &(q1->getRefToIPHyperViscoElastoPlasticAM()),dPVEVPAMdFAM,stiff);
assemblePParti(PVEVPAM, FfAM, FthAM, Fn, PAM);
@@ -1765,6 +1855,24 @@ void mlawPhenomenologicalSMP::constitutive1
assembleDDissParti(dViscVEVPAMdFVEVP, dViscVEVPAMdT, FfAM, FthAM, dFthAMdF, dFthAMdT, Fn, dViscAMdF, dViscAMdT);
+ //p-R: polynomial 2nd degree
+ static STensor3 FRL20, FRL2n, PVEVPRL2, dPVEVPRL2dT;
+ static STensor43 dPVEVPRL2dFRL2;
+ STensorOperation::zero(PVEVPRL2);
+ STensorOperation::zero(dPVEVPRL2dT);
+ STensorOperation::zero(dPVEVPRL2dFRL2);
+ STensorOperation::zero(FRL20);
+ STensorOperation::zero(FRL2n);
+
+ evaluateFParti(Fn, FfR, FthR, FRL2n);
+ evaluateFParti(F0, FfR0, FthR0, FRL20);
+ if(_EPFunc_R!=NULL)
+ {
+ _EPFunc_R->constitutive(FRL2n, PVEVPRL2, stiff, &dPVEVPRL2dFRL2);
+ PsiRL2= _EPFunc_R->get(FRL2n);
+ assemblePParti(PVEVPRL2, FfR, FthR, Fn, PRL2);
+ assembleTangentParti(PRL2, dPVEVPRL2dFRL2, dPVEVPRL2dT, FfR, FthR, dFthRdF, dFthRdT, Fn, TangentRL2, dPRL2dT);
+ }
//p-AM: polynomial 2nd degree
static STensor3 FAML20, FAML2n, PVEVPAML2, dPVEVPAML2dT;
@@ -1775,14 +1883,16 @@ void mlawPhenomenologicalSMP::constitutive1
STensorOperation::zero(FAML20);
STensorOperation::zero(FAML2n);
- evaluateFParti(Fn, _I2, FthAM, FAML2n);
- evaluateFParti(F0, _I2, FthAM0, FAML20);
+ evaluateFParti(Fn, FfAM, FthAM, FAML2n);
+ evaluateFParti(F0, FfAM0, FthAM0, FAML20);
- _EPFunc_AM->constitutive(FAML2n, PVEVPAML2, stiff, &dPVEVPAML2dFAML2);
- PsiAML2= _EPFunc_AM->get(FAML2n);
- assemblePParti(PVEVPAML2, _I2, FthAM, Fn, PAML2);
- assembleTangentParti(PAML2, dPVEVPAML2dFAML2, dPVEVPAML2dT, _I2, FthAM, dFthAMdF, dFthAMdT, Fn, TangentAML2, dPAML2dT);
-
+ if(_EPFunc_AM!=NULL)
+ {
+ _EPFunc_AM->constitutive(FAML2n, PVEVPAML2, stiff, &dPVEVPAML2dFAML2);
+ PsiAML2= _EPFunc_AM->get(FAML2n);
+ assemblePParti(PVEVPAML2, FfAM, FthAM, Fn, PAML2);
+ assembleTangentParti(PAML2, dPVEVPAML2dFAML2, dPVEVPAML2dT, FfAM, FthAM, dFthAMdF, dFthAMdT, Fn, TangentAML2, dPAML2dT);
+ }
//if(stiff)
//{
// tangentamorphousEP(Fn, F0, FmAM, stiff, FthAM, dFthAMdT, dFthAMdF, PsiAML2, PAML2, dPAML2dT, TangentAML2);
@@ -1840,9 +1950,25 @@ void mlawPhenomenologicalSMP::constitutive2
STensor3& PR, STensor3& dPRdT, STensor43& TangentR,
double& PsiAM, double &dPlasAM, double& dPlasAMdT, STensor3& dPlasAMdF, double &dViscAM, double& dViscAMdT, STensor3& dViscAMdF,
STensor3 &PAM, STensor3 &dPAMdT, STensor43 &TangentAM,
+ double &PsiRL2, STensor3 &PRL2, STensor3 &dPRL2dT, STensor43 &TangentRL2,
double &PsiAML2, STensor3 &PAML2, STensor3 &dPAML2dT, STensor43 &TangentAML2) const
{
// Initialisation of internal variables and material properties
+ double anisotropicThermalCoefAmplitudeBulkG=(_alphaParam)[0];
+ double anisotropicThermalCoefAmplitudeShearG=(_alphaParam)[4];
+ double anisotropicThermalAlphaCoefTanhTempBulkG=(_alphaParam)[8];
+ double anisotropicThermalAlphaCoefTanhTempShearG=(_alphaParam)[12];
+
+ double anisotropicThermalCoefAmplitudeBulkR=(_alphaParam)[1];
+ double anisotropicThermalCoefAmplitudeShearR=(_alphaParam)[5];
+ double anisotropicThermalAlphaCoefTanhTempBulkR=(_alphaParam)[9];
+ double anisotropicThermalAlphaCoefTanhTempShearR=(_alphaParam)[13];
+
+ double anisotropicThermalCoefAmplitudeBulkAM=(_alphaParam)[3];
+ double anisotropicThermalCoefAmplitudeShearAM=(_alphaParam)[7];
+ double anisotropicThermalAlphaCoefTanhTempBulkAM=(_alphaParam)[11];
+ double anisotropicThermalAlphaCoefTanhTempShearAM=(_alphaParam)[15];
+
double &Tdot =q1->getRefToTdot();
const STensor3& FthAD0 =q0->getConstRefToFthAD();
STensor3& FthAD =q1->getRefToFthAD();
@@ -1866,6 +1992,16 @@ void mlawPhenomenologicalSMP::constitutive2
evaluateFParti(Fn, FfG, FthG, FGn);
evaluateFParti(F0, FfG0, FthG0, FG0);
+ double wp=0., dwpdz=0.;
+ getZgFunction_Sy0(*_yieldCrystallinityG,q0->getzg(), wp, dwpdz);
+ q1->getRefToIPHyperViscoElastoPlasticG()._ipCompression->setScaleFactor(wp,0.);
+ q1->getRefToIPHyperViscoElastoPlasticG()._ipTraction->setScaleFactor(wp,0.);
+
+ funEnThAI(T0, q1->getConstRefToTRefG(), q0->getConstRefToIPHyperViscoElastoPlasticG(),q1->getRefToIPHyperViscoElastoPlasticG(),
+ anisotropicThermalCoefAmplitudeBulkG, anisotropicThermalCoefAmplitudeShearG,
+ anisotropicThermalAlphaCoefTanhTempBulkG, anisotropicThermalAlphaCoefTanhTempShearG);
+
+
lawG->constitutive(FG0,FGn, PVEVPG, &(q0->getConstRefToIPHyperViscoElastoPlasticG()), &(q1->getRefToIPHyperViscoElastoPlasticG()),dPVEVPGdFG,stiff);
assemblePParti(PVEVPG, FfG, FthG, Fn, PG);
@@ -1905,6 +2041,15 @@ void mlawPhenomenologicalSMP::constitutive2
evaluateFParti(Fn, FfR, FthR, FRn);
evaluateFParti(F0, FfR0, FthR0, FR0);
+
+ getZgFunction_Sy0(*_yieldCrystallinityR,q0->getzg(), wp, dwpdz);
+ q1->getRefToIPHyperViscoElastoPlasticR()._ipCompression->setScaleFactor(wp,0.);
+ q1->getRefToIPHyperViscoElastoPlasticR()._ipTraction->setScaleFactor(wp,0.);
+
+ funEnThAI(T0, q1->getConstRefToTRefR(), q0->getConstRefToIPHyperViscoElastoPlasticR(),q1->getRefToIPHyperViscoElastoPlasticR(),
+ anisotropicThermalCoefAmplitudeBulkR, anisotropicThermalCoefAmplitudeShearR,
+ anisotropicThermalAlphaCoefTanhTempBulkR, anisotropicThermalAlphaCoefTanhTempShearR);
+
lawR->constitutive(FR0,FRn, PVEVPR, &(q0->getConstRefToIPHyperViscoElastoPlasticR()), &(q1->getRefToIPHyperViscoElastoPlasticR()),dPVEVPRdFR,stiff);
@@ -1941,6 +2086,14 @@ void mlawPhenomenologicalSMP::constitutive2
evaluateFParti(Fn, FfAM, FthAM, FAMn);
evaluateFParti(F0, FfAM0, FthAM0, FAM0);
+ getZgFunction_Sy0(*_yieldCrystallinityAM,q0->getzg(), wp, dwpdz);
+ q1->getRefToIPHyperViscoElastoPlasticAM()._ipCompression->setScaleFactor(wp,0.);
+ q1->getRefToIPHyperViscoElastoPlasticAM()._ipTraction->setScaleFactor(wp,0.);
+
+ funEnThAI(T0, q1->getConstRefToTRefAM(), q0->getConstRefToIPHyperViscoElastoPlasticAM(),q1->getRefToIPHyperViscoElastoPlasticAM(),
+ anisotropicThermalCoefAmplitudeBulkAM, anisotropicThermalCoefAmplitudeShearAM,
+ anisotropicThermalAlphaCoefTanhTempBulkAM, anisotropicThermalAlphaCoefTanhTempShearAM);
+
lawAM->constitutive(FAM0,FAMn, PVEVPAM, &(q0->getConstRefToIPHyperViscoElastoPlasticAM()), &(q1->getRefToIPHyperViscoElastoPlasticAM()),dPVEVPAMdFAM,stiff);
assemblePParti(PVEVPAM, FfAM, FthAM, Fn, PAM);
@@ -1960,6 +2113,26 @@ void mlawPhenomenologicalSMP::constitutive2
double dViscVEVPAMdT=0.;
assembleDDissParti(dViscVEVPAMdFVEVP, dViscVEVPAMdT, FfAM, FthAM, dFthAMdF, dFthAMdT, Fn, dViscAMdF, dViscAMdT);
+ //p-R: polynomial 2nd degree
+ static STensor3 FRL20, FRL2n, PVEVPRL2, dPVEVPRL2dT;
+ static STensor43 dPVEVPRL2dFRL2;
+ STensorOperation::zero(PVEVPRL2);
+ STensorOperation::zero(dPVEVPRL2dT);
+ STensorOperation::zero(dPVEVPRL2dFRL2);
+ STensorOperation::zero(FRL20);
+ STensorOperation::zero(FRL2n);
+
+ evaluateFParti(Fn, FfR, FthR, FRL2n);
+ evaluateFParti(F0, FfR0, FthR0, FRL20);
+
+ if(_EPFunc_R!=NULL)
+ {
+ _EPFunc_R->constitutive(FRL2n, PVEVPRL2, stiff, &dPVEVPRL2dFRL2);
+ PsiRL2= _EPFunc_R->get(FRL2n);
+ assemblePParti(PVEVPRL2, FfR, FthR, Fn, PRL2);
+ assembleTangentParti(PRL2, dPVEVPRL2dFRL2, dPVEVPRL2dT, FfR, FthR, dFthRdF, dFthRdT, Fn, TangentRL2, dPRL2dT);
+ }
+
//p-AM: polynomial 2nd degree
static STensor3 FAML20, FAML2n, PVEVPAML2, dPVEVPAML2dT;
static STensor43 dPVEVPAML2dFAML2;
@@ -1969,13 +2142,17 @@ void mlawPhenomenologicalSMP::constitutive2
STensorOperation::zero(FAML20);
STensorOperation::zero(FAML2n);
- evaluateFParti(Fn, _I2, FthAM, FAML2n);
- evaluateFParti(F0, _I2, FthAM0, FAML20);
+ evaluateFParti(Fn, FfAM, FthAM, FAML2n);
+ evaluateFParti(F0, FfAM0, FthAM0, FAML20);
- _EPFunc_AM->constitutive(FAML2n, PVEVPAML2, stiff, &dPVEVPAML2dFAML2);
- PsiAML2= _EPFunc_AM->get(FAML2n);
- assemblePParti(PVEVPAML2, _I2, FthAM, Fn, PAML2);
- assembleTangentParti(PAML2, dPVEVPAML2dFAML2, dPVEVPAML2dT, _I2, FthAM, dFthAMdF, dFthAMdT, Fn, TangentAML2, dPAML2dT);
+ if(_EPFunc_AM!=NULL)
+ {
+ _EPFunc_AM->constitutive(FAML2n, PVEVPAML2, stiff, &dPVEVPAML2dFAML2);
+ PsiAML2= _EPFunc_AM->get(FAML2n);
+ assemblePParti(PVEVPAML2, FfAM, FthAM, Fn, PAML2);
+ assembleTangentParti(PAML2, dPVEVPAML2dFAML2, dPVEVPAML2dT, FfAM, FthAM, dFthAMdF, dFthAMdT, Fn, TangentAML2, dPAML2dT);
+ }
+
}
@@ -1994,9 +2171,25 @@ void mlawPhenomenologicalSMP::constitutive3
STensor3& PR, STensor3& dPRdT, STensor43& TangentR,
double& PsiAM, double &dPlasAM, double& dPlasAMdT, STensor3& dPlasAMdF, double &dViscAM, double& dViscAMdT, STensor3& dViscAMdF,
STensor3 &PAM, STensor3 &dPAMdT, STensor43 &TangentAM,
+ double &PsiRL2, STensor3 &PRL2, STensor3 &dPRL2dT, STensor43 &TangentRL2,
double &PsiAML2, STensor3 &PAML2, STensor3 &dPAML2dT, STensor43 &TangentAML2) const
{
// Initialisation of internal variables and material properties
+ double anisotropicThermalCoefAmplitudeBulkG=(_alphaParam)[0];
+ double anisotropicThermalCoefAmplitudeShearG=(_alphaParam)[4];
+ double anisotropicThermalAlphaCoefTanhTempBulkG=(_alphaParam)[8];
+ double anisotropicThermalAlphaCoefTanhTempShearG=(_alphaParam)[12];
+
+ double anisotropicThermalCoefAmplitudeBulkR=(_alphaParam)[1];
+ double anisotropicThermalCoefAmplitudeShearR=(_alphaParam)[5];
+ double anisotropicThermalAlphaCoefTanhTempBulkR=(_alphaParam)[9];
+ double anisotropicThermalAlphaCoefTanhTempShearR=(_alphaParam)[13];
+
+ double anisotropicThermalCoefAmplitudeBulkAM=(_alphaParam)[3];
+ double anisotropicThermalCoefAmplitudeShearAM=(_alphaParam)[7];
+ double anisotropicThermalAlphaCoefTanhTempBulkAM=(_alphaParam)[11];
+ double anisotropicThermalAlphaCoefTanhTempShearAM=(_alphaParam)[15];
+
double &Tdot =q1->getRefToTdot();
const STensor3& FthAD0 =q0->getConstRefToFthAD();
STensor3& FthAD =q1->getRefToFthAD();
@@ -2026,6 +2219,15 @@ void mlawPhenomenologicalSMP::constitutive3
evaluateFParti(Fn, FfG, FthG, FGn);
evaluateFParti(F0, FfG0, FthG0, FG0);
+ double wp=0., dwpdz=0.;
+ getZgFunction_Sy0(*_yieldCrystallinityG,q0->getzg(), wp, dwpdz);
+ q1->getRefToIPHyperViscoElastoPlasticG()._ipCompression->setScaleFactor(wp,0.);
+ q1->getRefToIPHyperViscoElastoPlasticG()._ipTraction->setScaleFactor(wp,0.);
+
+ funEnThAI(T0, q1->getConstRefToTRefG(), q0->getConstRefToIPHyperViscoElastoPlasticG(),q1->getRefToIPHyperViscoElastoPlasticG(),
+ anisotropicThermalCoefAmplitudeBulkG, anisotropicThermalCoefAmplitudeShearG,
+ anisotropicThermalAlphaCoefTanhTempBulkG, anisotropicThermalAlphaCoefTanhTempShearG);
+
lawG->constitutive(FG0,FGn, PVEVPG, &(q0->getConstRefToIPHyperViscoElastoPlasticG()), &(q1->getRefToIPHyperViscoElastoPlasticG()),dPVEVPGdFG,stiff);
assemblePParti(PVEVPG, FfG, FthG, Fn, PG);
@@ -2097,6 +2299,14 @@ void mlawPhenomenologicalSMP::constitutive3
evaluateFParti(Fn, FfAM, FthAM, FAMn);
evaluateFParti(F0, FfAM0, FthAM0, FAM0);
+ getZgFunction_Sy0(*_yieldCrystallinityAM,q0->getzg(), wp, dwpdz);
+ q1->getRefToIPHyperViscoElastoPlasticAM()._ipCompression->setScaleFactor(wp,0.);
+ q1->getRefToIPHyperViscoElastoPlasticAM()._ipTraction->setScaleFactor(wp,0.);
+
+ funEnThAI(T0, q1->getConstRefToTRefAM(), q0->getConstRefToIPHyperViscoElastoPlasticAM(),q1->getRefToIPHyperViscoElastoPlasticAM(),
+ anisotropicThermalCoefAmplitudeBulkAM, anisotropicThermalCoefAmplitudeShearAM,
+ anisotropicThermalAlphaCoefTanhTempBulkAM, anisotropicThermalAlphaCoefTanhTempShearAM);
+
lawAM->constitutive(FAM0,FAMn, PVEVPAM, &(q0->getConstRefToIPHyperViscoElastoPlasticAM()), &(q1->getRefToIPHyperViscoElastoPlasticAM()),dPVEVPAMdFAM,stiff);
assemblePParti(PVEVPAM, FfAM, FthAM, Fn, PAM);
@@ -2116,6 +2326,26 @@ void mlawPhenomenologicalSMP::constitutive3
double dViscVEVPAMdT=0.;
assembleDDissParti(dViscVEVPAMdFVEVP, dViscVEVPAMdT, FfAM, FthAM, dFthAMdF, dFthAMdT, Fn, dViscAMdF, dViscAMdT);
+ //p-R: polynomial 2nd degree
+ static STensor3 FRL20, FRL2n, PVEVPRL2, dPVEVPRL2dT;
+ static STensor43 dPVEVPRL2dFRL2;
+ STensorOperation::zero(PVEVPRL2);
+ STensorOperation::zero(dPVEVPRL2dT);
+ STensorOperation::zero(dPVEVPRL2dFRL2);
+ STensorOperation::zero(FRL20);
+ STensorOperation::zero(FRL2n);
+
+ evaluateFParti(Fn, FfR, FthR, FRL2n);
+ evaluateFParti(F0, FfR0, FthR0, FRL20);
+
+ if(_EPFunc_R!=NULL)
+ {
+ _EPFunc_R->constitutive(FRL2n, PVEVPRL2, stiff, &dPVEVPRL2dFRL2);
+ PsiRL2= _EPFunc_R->get(FRL2n);
+ assemblePParti(PVEVPRL2, FfR, FthR, Fn, PRL2);
+ assembleTangentParti(PRL2, dPVEVPRL2dFRL2, dPVEVPRL2dT, FfR, FthR, dFthRdF, dFthRdT, Fn, TangentRL2, dPRL2dT);
+ }
+
//p-AM: polynomial 2nd degree
static STensor3 FAML20, FAML2n, PVEVPAML2, dPVEVPAML2dT;
static STensor43 dPVEVPAML2dFAML2;
@@ -2125,13 +2355,16 @@ void mlawPhenomenologicalSMP::constitutive3
STensorOperation::zero(FAML20);
STensorOperation::zero(FAML2n);
- evaluateFParti(Fn, _I2, FthAM, FAML2n);
- evaluateFParti(F0, _I2, FthAM0, FAML20);
+ evaluateFParti(Fn, FfAM, FthAM, FAML2n);
+ evaluateFParti(F0, FfAM0, FthAM0, FAML20);
- _EPFunc_AM->constitutive(FAML2n, PVEVPAML2, stiff, &dPVEVPAML2dFAML2);
- PsiAML2= _EPFunc_AM->get(FAML2n);
- assemblePParti(PVEVPAML2, _I2, FthAM, Fn, PAML2);
- assembleTangentParti(PAML2, dPVEVPAML2dFAML2, dPVEVPAML2dT, _I2, FthAM, dFthAMdF, dFthAMdT, Fn, TangentAML2, dPAML2dT);
+ if(_EPFunc_AM!=NULL)
+ {
+ _EPFunc_AM->constitutive(FAML2n, PVEVPAML2, stiff, &dPVEVPAML2dFAML2);
+ PsiAML2= _EPFunc_AM->get(FAML2n);
+ assemblePParti(PVEVPAML2, FfAM, FthAM, Fn, PAML2);
+ assembleTangentParti(PAML2, dPVEVPAML2dFAML2, dPVEVPAML2dT, FfAM, FthAM, dFthAMdF, dFthAMdT, Fn, TangentAML2, dPAML2dT);
+ }
// TangentG.print("con1 TanG");
// TangentAM.print("con1 TanAM");
@@ -3220,14 +3453,14 @@ void mlawPhenomenologicalSMP::computePi
const STensor3 &A1i0, const STensor3 &A2i0, const double B1i0, const double B2i0,
STensor3 &A1i, STensor3 &A2i, double &B1i, double &B2i,
STensor3& devq1i, double& pq1i, STensor3& devq2i, double& pq2i,
- double anisotropicThermalCoefAmplitude, double anisotropicThermalCoefUnused,
- double anisotropicThermalAlphaCoefTanhEn, double anisotropicThermalAlphaCoefTanhTemp, double TRef,
+ double anisotropicThermalCoefAmplitudeBulk, double anisotropicThermalCoefAmplitudeShear,
+ double anisotropicThermalAlphaCoefTanhTempBulk, double anisotropicThermalAlphaCoefTanhTempShear , double TRef,
STensor3 &strainEffect, STensor3 &dStrainEffectDt, STensor43 &dStrainEffectdEne) const
{
/* double MU1i =(_ViscParami)[0];
double KAPPA1i =(_ViscParami)[1];
double TAU1i =(_ViscParami)[2];
- double MU2i =(_ViscParami)[3];
+ double MU2i =(_ViscParami)[3];double
double KAPPA2i =(_ViscParami)[4];
double TAU2i =(_ViscParami)[5];
diff --git a/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.h b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.h
index 7d0c027d555a23f98fc44d87032088d354a1148b..34ec2af98a6d030d4ea0374c5c5419d729e736ae 100644
--- a/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.h
+++ b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.h
@@ -53,7 +53,12 @@ protected:
scalarFunction *_temFunc_cG;
scalarFunction *_temFunc_cR;
scalarFunction *_temFunc_cAM;
+
+ scalarFunction *_yieldCrystallinityG;
+ scalarFunction *_yieldCrystallinityR;
+ scalarFunction *_yieldCrystallinityAM;
+ elasticPotential *_EPFunc_R;
elasticPotential *_EPFunc_AM;
public:
@@ -196,62 +201,62 @@ public:
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
return lawAM->getViscoElasticNumberOfElement();
};
- virtual void setPowerFactorG(const double n) const
+ virtual void setPowerFactorG(const double n)
{
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
lawG->setPowerFactor(n);
};
- virtual void setPowerFactorR(const double n) const
+ virtual void setPowerFactorR(const double n)
{
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
lawR->setPowerFactor(n);
};
- virtual void setPowerFactorAM(const double n) const
+ virtual void setPowerFactorAM(const double n)
{
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
lawAM->setPowerFactor(n);
};
- virtual void nonAssociatedFlowRuleFactorG(const double n) const
+ virtual void nonAssociatedFlowRuleFactorG(const double n)
{
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
lawG->nonAssociatedFlowRuleFactor(n);
};
- virtual void nonAssociatedFlowRuleFactorR(const double n) const
+ virtual void nonAssociatedFlowRuleFactorR(const double n)
{
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
lawR->nonAssociatedFlowRuleFactor(n);
};
- virtual void nonAssociatedFlowRuleFactorAM(const double n) const
+ virtual void nonAssociatedFlowRuleFactorAM(const double n)
{
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
lawAM->nonAssociatedFlowRuleFactor(n);
};
- virtual void setPlasticPoissonRatioG(const double n) const
+ virtual void setPlasticPoissonRatioG(const double n)
{
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
lawG->setPlasticPoissonRatio(n);
};
- virtual void setPlasticPoissonRatioR(const double n) const
+ virtual void setPlasticPoissonRatioR(const double n)
{
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
lawR->setPlasticPoissonRatio(n);
};
- virtual void setPlasticPoissonRatioAM(const double n) const
+ virtual void setPlasticPoissonRatioAM(const double n)
{
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
lawAM->setPlasticPoissonRatio(n);
};
- virtual void setNonAssociatedFlowG(const bool n) const
+ virtual void setNonAssociatedFlowG(const bool n)
{
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
lawG->setNonAssociatedFlow(n);
};
- virtual void setNonAssociatedFlowR(const bool n) const
+ virtual void setNonAssociatedFlowR(const bool n)
{
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
lawR->setNonAssociatedFlow(n);
};
- virtual void setNonAssociatedFlowAM(const bool n) const
+ virtual void setNonAssociatedFlowAM(const bool n)
{
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
lawAM->setNonAssociatedFlow(n);
@@ -261,7 +266,7 @@ public:
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
return lawG->getConstRefToCompressionHardening();
};
- virtual void setJ2IsotropicHardeningCompressionG(const J2IsotropicHardening& isoHard) const {
+ virtual void setJ2IsotropicHardeningCompressionG(const J2IsotropicHardening& isoHard) {
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
lawG->setCompressionHardening(isoHard);
};
@@ -269,7 +274,7 @@ public:
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
return lawG->getConstRefToTractionHardening();
};
- virtual void setJ2IsotropicHardeningTractionG(const J2IsotropicHardening& isoHard) const {
+ virtual void setJ2IsotropicHardeningTractionG(const J2IsotropicHardening& isoHard) {
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
lawG->setTractionHardening(isoHard);
};
@@ -277,7 +282,7 @@ public:
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
return lawG->getConstRefToKinematicHardening();
};
- virtual void setKinematicHardeningG(const kinematicHardening& isoHard) const {
+ virtual void setKinematicHardeningG(const kinematicHardening& isoHard) {
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
lawG->setKinematicHardening(isoHard);
};
@@ -286,7 +291,7 @@ public:
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
return lawR->getConstRefToCompressionHardening();
};
- virtual void setJ2IsotropicHardeningCompressionR(const J2IsotropicHardening& isoHard) const {
+ virtual void setJ2IsotropicHardeningCompressionR(const J2IsotropicHardening& isoHard) {
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
lawR->setCompressionHardening(isoHard);
};
@@ -294,7 +299,7 @@ public:
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
return lawR->getConstRefToTractionHardening();
};
- virtual void setJ2IsotropicHardeningTractionR(const J2IsotropicHardening& isoHard) const {
+ virtual void setJ2IsotropicHardeningTractionR(const J2IsotropicHardening& isoHard) {
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
lawR->setTractionHardening(isoHard);
};
@@ -302,7 +307,7 @@ public:
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
return lawR->getConstRefToKinematicHardening();
};
- virtual void setKinematicHardeningR(const kinematicHardening& isoHard) const {
+ virtual void setKinematicHardeningR(const kinematicHardening& isoHard) {
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
lawR->setKinematicHardening(isoHard);
};
@@ -311,7 +316,7 @@ public:
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
return lawAM->getConstRefToCompressionHardening();
};
- virtual void setJ2IsotropicHardeningCompressionAM(const J2IsotropicHardening& isoHard) const {
+ virtual void setJ2IsotropicHardeningCompressionAM(const J2IsotropicHardening& isoHard) {
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
lawAM->setCompressionHardening(isoHard);
};
@@ -319,7 +324,7 @@ public:
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
return lawAM->getConstRefToTractionHardening();
};
- virtual void setJ2IsotropicHardeningTractionAM(const J2IsotropicHardening& isoHard) const {
+ virtual void setJ2IsotropicHardeningTractionAM(const J2IsotropicHardening& isoHard) {
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
lawAM->setTractionHardening(isoHard);
};
@@ -327,23 +332,35 @@ public:
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
return lawAM->getConstRefToKinematicHardening();
};
- virtual void setKinematicHardeningAM(const kinematicHardening& isoHard) const {
+ virtual void setKinematicHardeningAM(const kinematicHardening& isoHard) {
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
lawAM->setKinematicHardening(isoHard);
};
- virtual void setViscosityEffectG(const viscosityLaw& isoHard, const double p) const {
+ virtual void setViscosityEffectG(const viscosityLaw& isoHard, const double p) {
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
lawG->setViscosityEffect(isoHard,p);
};
- virtual void setViscosityEffectR(const viscosityLaw& isoHard, const double p) const {
+ virtual void setViscosityEffectR(const viscosityLaw& isoHard, const double p) {
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
lawR->setViscosityEffect(isoHard,p);
};
- virtual void setViscosityEffectAM(const viscosityLaw& isoHard, const double p) const {
+ virtual void setViscosityEffectAM(const viscosityLaw& isoHard, const double p) {
if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
lawAM->setViscosityEffect(isoHard,p);
};
+ virtual void setVolumeCorrectionG(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc) {
+ if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
+ lawG->setVolumeCorrection(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
+ };
+ virtual void setVolumeCorrectionR(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc) {
+ if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
+ lawR->setVolumeCorrection(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
+ };
+ virtual void setVolumeCorrectionAM(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc) {
+ if(lawAM==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper AM not defined");
+ lawAM->setVolumeCorrection(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
+ };
virtual bool withEnergyDissipation() const {
if(lawG==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper G not defined");
if(lawR==NULL) Msg::Error("mlawPhenomenologicalSMP: mlawPowerYieldHyper R not defined");
@@ -377,8 +394,14 @@ public:
if(_temFunc_cG!=NULL) delete _temFunc_cG;
if(_temFunc_cR!=NULL) delete _temFunc_cR;
if(_temFunc_cAM!=NULL) delete _temFunc_cAM;
+
+ if(_yieldCrystallinityG!=NULL) delete _yieldCrystallinityG;
+ if(_yieldCrystallinityR!=NULL) delete _yieldCrystallinityR;
+ if(_yieldCrystallinityAM!=NULL) delete _yieldCrystallinityAM;
+
if(_EPFunc_AM!=NULL) delete _EPFunc_AM;
+ if(_EPFunc_R!=NULL) delete _EPFunc_R;
};
virtual materialLaw* clone() const {return new mlawPhenomenologicalSMP(*this);};
@@ -556,7 +579,25 @@ public:
_temFunc_cAM = Tfunc.clone();
}
- virtual void getZgFunction_Sy0G(const double zg, double &wp, double &dwpdz) const;
+ virtual void setFunctionYieldCrystallinityG(const scalarFunction& Tfunc)
+ {
+ if(_yieldCrystallinityG!=NULL) delete _yieldCrystallinityG;
+ _yieldCrystallinityG=Tfunc.clone();
+ }
+
+ virtual void setFunctionYieldCrystallinityR(const scalarFunction& Tfunc)
+ {
+ if(_yieldCrystallinityR!=NULL) delete _yieldCrystallinityR;
+ _yieldCrystallinityR=Tfunc.clone();
+ }
+
+ virtual void setFunctionYieldCrystallinityAM(const scalarFunction& Tfunc)
+ {
+ if(_yieldCrystallinityAM!=NULL) delete _yieldCrystallinityAM;
+ _yieldCrystallinityAM=Tfunc.clone();
+ }
+
+ virtual void getZgFunction_Sy0(const scalarFunction& Tfunc, const double zg, double &wp, double &dwpdz) const;
virtual void get_zg (const double T, const double T0, const IPPhenomenologicalSMP *q0, IPPhenomenologicalSMP *q1,
@@ -684,11 +725,16 @@ public:
_temFunc_cAM = Tfunc.clone();
}
- virtual void setElasticPotentialFunction(const elasticPotential& EPfunc)
+ virtual void setElasticPotentialFunctionAM(const elasticPotential& EPfunc)
{
if (_EPFunc_AM != NULL) delete _EPFunc_AM;
_EPFunc_AM = EPfunc.clone();
}
+ virtual void setElasticPotentialFunctionR(const elasticPotential& EPfunc)
+ {
+ if (_EPFunc_R != NULL) delete _EPFunc_R;
+ _EPFunc_R = EPfunc.clone();
+ }
void setTcTmWcWm(const double Atc, const double Btm, const double Cwc, const double Dwm,const double alphatc, const double alphatm, const double alphawc, const double alphawm)
@@ -777,6 +823,7 @@ public:
STensor3& PR, STensor3& dPRdT, STensor43& TangentR,
double& PsiAM, double &dPlasAM, double& dPlasAMdT, STensor3& dPlasAMdF, double &dViscAM, double& dViscAMdT, STensor3& dViscAMdF,
STensor3 &PAM, STensor3 &dPAMdT, STensor43 &TangentAM,
+ double &PsiRL2, STensor3 &PRL2, STensor3 &dPRL2dT, STensor43 &TangentRL2,
double &PsiAML2, STensor3 &PAML2, STensor3 &dPAML2dT, STensor43 &TangentAML2) const;
void constitutive2(
@@ -791,6 +838,7 @@ public:
STensor3& PR, STensor3& dPRdT, STensor43& TangentR,
double& PsiAM, double &dPlasAM, double& dPlasAMdT, STensor3& dPlasAMdF, double &dViscAM, double& dViscAMdT, STensor3& dViscAMdF,
STensor3 &PAM, STensor3 &dPAMdT, STensor43 &TangentAM,
+ double &PsiRL2, STensor3 &PRL2, STensor3 &dPRL2dT, STensor43 &TangentRL2,
double &PsiAML2, STensor3 &PAML2, STensor3 &dPAML2dT, STensor43 &TangentAML2) const;
void constitutive3(
@@ -805,6 +853,7 @@ public:
STensor3& PR, STensor3& dPRdT, STensor43& TangentR,
double& PsiAM, double &dPlasAM, double& dPlasAMdT, STensor3& dPlasAMdF, double &dViscAM, double& dViscAMdT, STensor3& dViscAMdF,
STensor3 &PAM, STensor3 &dPAMdT, STensor43 &TangentAM,
+ double &PsiRL2, STensor3 &PRL2, STensor3 &dPRL2dT, STensor43 &TangentRL2,
double &PsiAML2, STensor3 &PAML2, STensor3 &dPAML2dT, STensor43 &TangentAML2) const;
const STensor3& getInitialConductivityTensor()const {return _k0;};
@@ -849,11 +898,17 @@ protected:
const double T0, const double T, const double Tcrys, const double Tmelt, const double zg, const double dzgdT,
STensor3& FthI, STensor3& dFthIdT, STensor43& dFthIdF) const;
- void funEnThAI
- (const double T, const double TRefR,
- const STensor3& Ene, STensor3& strainEffect, STensor3& dStrainEffectdT, STensor43& dStrainEffectdEne,
- double anisotropicThermalCoefAmplitude, double anisotropicThermalCoefUnused,
- double anisotropicThermalAlphaCoefTanhEn, double anisotropicThermalAlphaCoefTanhTemp) const;
+ //void funEnThAI
+ //(const double T, const double TRefR,
+ // const STensor3& Ene, STensor3& strainEffect, STensor3& dStrainEffectdT, STensor43& dStrainEffectdEne,
+ // double anisotropicThermalCoefAmplitude, double anisotropicThermalCoefUnused,
+ // double anisotropicThermalAlphaCoefTanhEn, double anisotropicThermalAlphaCoefTanhTemp) const;
+
+
+ void funEnThAI(const double T, const double TRef, const IPHyperViscoElastic& q0, IPHyperViscoElastic &q1,
+ double anisotropicThermalCoefAmplitudeBulkAM, double anisotropicThermalCoefAmplitudeShearAM,
+ double anisotropicThermalAlphaCoefTanhTempBulkAM, double anisotropicThermalAlphaCoefTanhTempShearAM) const;
+
void predictorCorrectori(
const STensor3& Ffi0, STensor3& Ffi, const STensor3& Fpi0, STensor3& Fpi, const double &gammai0, double &gammai1, const STensor3& A1i0, const STensor3& A2i0,
diff --git a/dG3D/benchmarks/CMakeLists.txt b/dG3D/benchmarks/CMakeLists.txt
index 07f97d664a40d4d5fdaff2e06e249a17c0be753c..8c0e7e92076c732883f631da391f53f310e84e4f 100644
--- a/dG3D/benchmarks/CMakeLists.txt
+++ b/dG3D/benchmarks/CMakeLists.txt
@@ -162,6 +162,7 @@ add_subdirectory(restartMixedmodeDelaminationExpoLaw)
add_subdirectory(hyperViscoElastic)
add_subdirectory(powerYieldViscoElastoPlastic)
add_subdirectory(powerYieldViscoElastic)
+add_subdirectory(powerYieldViscoElastWithCorrection)
add_subdirectory(powerYieldViscoElastoPlasticSaturationDamage)
add_subdirectory(powerYieldViscoElastoPlasticFullFailure)
add_subdirectory(J2plasticExtractCohesiveLaw)
diff --git a/dG3D/benchmarks/powerYieldViscoElastWithCorrection/CMakeLists.txt b/dG3D/benchmarks/powerYieldViscoElastWithCorrection/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3667d9747b7c516dd7833cc4e59989e55b267c01
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastWithCorrection/CMakeLists.txt
@@ -0,0 +1,11 @@
+# test file
+
+set(PYFILE cylinder.py)
+
+set(FILES2DELETE
+ disp*.msh
+ stress*.msh
+ *.csv
+)
+
+add_cm3python_test(${PYFILE} "${FILES2DELETE}")
diff --git a/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.geo b/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.geo
new file mode 100644
index 0000000000000000000000000000000000000000..b7914ad13f66e2bf45173c552530ce3b3c7ba8a7
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.geo
@@ -0,0 +1,27 @@
+mm = 1;
+D = 12.*mm;
+H = 12.*mm;
+
+R = D/2.;
+lsca = 0.15*D;
+
+Point(1) = {0,0,0,lsca};
+Point(2) = {R,0,0,lsca};
+Point(3) = {0,R,0,lsca};
+
+Circle(1) = {3, 1, 2};
+Line(2) = {3, 1};
+Line(3) = {1, 2};
+Line Loop(4) = {1, -3, -2};
+Plane Surface(5) = {4};
+Extrude {0, 0, H/2} {
+ Surface{5}; Layers{5}; //Recombine;
+}
+Physical Volume(11) = {1};
+Physical Surface(1) = {5};
+Physical Surface(2) = {22};
+Physical Surface(3) = {17};
+Physical Surface(4) = {21};
+Physical Point(5) = {5};
+
+//Recombine Surface {5};
diff --git a/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.msh b/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.msh
new file mode 100644
index 0000000000000000000000000000000000000000..0f52bfc5a294857b8f304b1da6e258fb7505d374
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.msh
@@ -0,0 +1,698 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$Nodes
+132
+1 0 0 0
+2 6 0 0
+3 0 6 0
+4 0 6 6
+5 0 0 6
+6 6 0 6
+7 1.552914270611066 5.795554957735497 0
+8 2.999999999993018 5.196152422710663 0
+9 4.242640687111061 4.242640687127508 0
+10 5.196152422702976 3.000000000006332 0
+11 5.795554957733712 1.552914270617726 0
+12 0 4.500000000005166 0
+13 0 3.000000000008336 0
+14 0 1.500000000004214 0
+15 1.499999999993832 0 0
+16 2.999999999988567 0 0
+17 4.499999999994237 0 0
+18 1.552914270611066 5.795554957735497 6
+19 2.999999999993018 5.196152422710663 6
+20 4.242640687111061 4.242640687127508 6
+21 5.196152422702976 3.000000000006332 6
+22 5.795554957733712 1.552914270617726 6
+23 4.499999999994237 0 6
+24 2.999999999988567 0 6
+25 1.499999999993832 0 6
+26 0 1.500000000004214 6
+27 0 3.000000000008336 6
+28 0 4.500000000005166 6
+29 0 6 1.2
+30 0 6 2.4
+31 0 6 3.6
+32 0 6 4.800000000000001
+33 6 0 1.2
+34 6 0 2.4
+35 6 0 3.6
+36 6 0 4.800000000000001
+37 0 0 1.2
+38 0 0 2.4
+39 0 0 3.6
+40 0 0 4.800000000000001
+41 2.527448697695314 2.575668272009054 0
+42 1.627841394027458 3.980250003805291 0
+43 3.914833052601809 1.621430334871464 0
+44 1.441200958338174 1.441200958338581 0
+45 3.760905731444732 3.080330269885938 0
+46 2.682328465952599 1.220004215439306 0
+47 1.220004215434672 2.682328465952725 0
+48 3.000820586356653 3.856980823767017 0
+49 1.552914270611066 5.795554957735497 1.2
+50 1.552914270611066 5.795554957735497 2.4
+51 1.552914270611066 5.795554957735497 3.6
+52 1.552914270611066 5.795554957735497 4.800000000000001
+53 2.999999999993018 5.196152422710663 1.2
+54 2.999999999993018 5.196152422710663 2.4
+55 2.999999999993018 5.196152422710663 3.6
+56 2.999999999993018 5.196152422710663 4.800000000000001
+57 4.242640687111061 4.242640687127508 1.2
+58 4.242640687111061 4.242640687127508 2.4
+59 4.242640687111061 4.242640687127508 3.6
+60 4.242640687111061 4.242640687127508 4.800000000000001
+61 5.196152422702976 3.000000000006332 1.2
+62 5.196152422702976 3.000000000006332 2.4
+63 5.196152422702976 3.000000000006332 3.6
+64 5.196152422702976 3.000000000006332 4.800000000000001
+65 5.795554957733712 1.552914270617726 1.2
+66 5.795554957733712 1.552914270617726 2.4
+67 5.795554957733712 1.552914270617726 3.6
+68 5.795554957733712 1.552914270617726 4.800000000000001
+69 1.499999999993832 0 1.2
+70 1.499999999993832 0 2.4
+71 1.499999999993832 0 3.6
+72 1.499999999993832 0 4.800000000000001
+73 2.999999999988567 0 1.2
+74 2.999999999988567 0 2.4
+75 2.999999999988567 0 3.6
+76 2.999999999988567 0 4.800000000000001
+77 4.499999999994237 0 1.2
+78 4.499999999994237 0 2.4
+79 4.499999999994237 0 3.6
+80 4.499999999994237 0 4.800000000000001
+81 0 4.500000000005166 1.2
+82 0 4.500000000005166 2.4
+83 0 4.500000000005166 3.6
+84 0 4.500000000005166 4.800000000000001
+85 0 3.000000000008336 1.2
+86 0 3.000000000008336 2.4
+87 0 3.000000000008336 3.6
+88 0 3.000000000008336 4.800000000000001
+89 0 1.500000000004214 1.2
+90 0 1.500000000004214 2.4
+91 0 1.500000000004214 3.6
+92 0 1.500000000004214 4.800000000000001
+93 2.527448697695314 2.575668272009054 6
+94 1.627841394027458 3.980250003805291 6
+95 3.914833052601809 1.621430334871464 6
+96 1.441200958338174 1.441200958338581 6
+97 3.760905731444732 3.080330269885938 6
+98 2.682328465952599 1.220004215439306 6
+99 1.220004215434672 2.682328465952725 6
+100 3.000820586356653 3.856980823767017 6
+101 2.527448697695314 2.575668272009054 1.2
+102 2.527448697695314 2.575668272009054 2.4
+103 2.527448697695314 2.575668272009054 3.6
+104 2.527448697695314 2.575668272009054 4.800000000000001
+105 1.627841394027458 3.980250003805291 1.2
+106 1.627841394027458 3.980250003805291 2.4
+107 1.627841394027458 3.980250003805291 3.6
+108 1.627841394027458 3.980250003805291 4.800000000000001
+109 3.914833052601809 1.621430334871464 1.2
+110 3.914833052601809 1.621430334871464 2.4
+111 3.914833052601809 1.621430334871464 3.6
+112 3.914833052601809 1.621430334871464 4.800000000000001
+113 1.441200958338174 1.441200958338581 1.2
+114 1.441200958338174 1.441200958338581 2.4
+115 1.441200958338174 1.441200958338581 3.6
+116 1.441200958338174 1.441200958338581 4.800000000000001
+117 3.760905731444732 3.080330269885938 1.2
+118 3.760905731444732 3.080330269885938 2.4
+119 3.760905731444732 3.080330269885938 3.6
+120 3.760905731444732 3.080330269885938 4.800000000000001
+121 2.682328465952599 1.220004215439306 1.2
+122 2.682328465952599 1.220004215439306 2.4
+123 2.682328465952599 1.220004215439306 3.6
+124 2.682328465952599 1.220004215439306 4.800000000000001
+125 1.220004215434672 2.682328465952725 1.2
+126 1.220004215434672 2.682328465952725 2.4
+127 1.220004215434672 2.682328465952725 3.6
+128 1.220004215434672 2.682328465952725 4.800000000000001
+129 3.000820586356653 3.856980823767017 1.2
+130 3.000820586356653 3.856980823767017 2.4
+131 3.000820586356653 3.856980823767017 3.6
+132 3.000820586356653 3.856980823767017 4.800000000000001
+$EndNodes
+$Elements
+557
+1 15 2 5 5 5
+2 2 2 1 5 7 42 12
+3 2 2 1 5 11 17 43
+4 2 2 1 5 16 43 17
+5 2 2 1 5 12 42 13
+6 2 2 1 5 1 44 15
+7 2 2 1 5 1 14 44
+8 2 2 1 5 3 7 12
+9 2 2 1 5 11 2 17
+10 2 2 1 5 16 46 43
+11 2 2 1 5 13 42 47
+12 2 2 1 5 7 8 42
+13 2 2 1 5 10 11 43
+14 2 2 1 5 41 42 48
+15 2 2 1 5 10 43 45
+16 2 2 1 5 15 46 16
+17 2 2 1 5 13 47 14
+18 2 2 1 5 41 47 42
+19 2 2 1 5 41 43 46
+20 2 2 1 5 41 45 43
+21 2 2 1 5 15 44 46
+22 2 2 1 5 14 47 44
+23 2 2 1 5 8 48 42
+24 2 2 1 5 9 10 45
+25 2 2 1 5 41 48 45
+26 2 2 1 5 8 9 48
+27 2 2 1 5 41 46 44
+28 2 2 1 5 41 44 47
+29 2 2 1 5 9 45 48
+30 2 2 3 17 37 69 1
+31 2 2 3 17 1 69 15
+32 2 2 3 17 38 70 37
+33 2 2 3 17 37 70 69
+34 2 2 3 17 39 71 38
+35 2 2 3 17 38 71 70
+36 2 2 3 17 40 72 39
+37 2 2 3 17 39 72 71
+38 2 2 3 17 5 25 40
+39 2 2 3 17 40 25 72
+40 2 2 3 17 69 73 15
+41 2 2 3 17 15 73 16
+42 2 2 3 17 70 74 69
+43 2 2 3 17 69 74 73
+44 2 2 3 17 71 75 70
+45 2 2 3 17 70 75 74
+46 2 2 3 17 72 76 71
+47 2 2 3 17 71 76 75
+48 2 2 3 17 25 24 72
+49 2 2 3 17 72 24 76
+50 2 2 3 17 73 77 16
+51 2 2 3 17 16 77 17
+52 2 2 3 17 74 78 73
+53 2 2 3 17 73 78 77
+54 2 2 3 17 75 79 74
+55 2 2 3 17 74 79 78
+56 2 2 3 17 76 80 75
+57 2 2 3 17 75 80 79
+58 2 2 3 17 24 23 76
+59 2 2 3 17 76 23 80
+60 2 2 3 17 77 2 17
+61 2 2 3 17 77 33 2
+62 2 2 3 17 78 33 77
+63 2 2 3 17 78 34 33
+64 2 2 3 17 79 34 78
+65 2 2 3 17 79 35 34
+66 2 2 3 17 80 35 79
+67 2 2 3 17 80 36 35
+68 2 2 3 17 23 36 80
+69 2 2 3 17 23 6 36
+70 2 2 4 21 29 81 3
+71 2 2 4 21 3 81 12
+72 2 2 4 21 30 82 29
+73 2 2 4 21 29 82 81
+74 2 2 4 21 31 83 30
+75 2 2 4 21 30 83 82
+76 2 2 4 21 32 84 31
+77 2 2 4 21 31 84 83
+78 2 2 4 21 4 28 32
+79 2 2 4 21 32 28 84
+80 2 2 4 21 81 85 12
+81 2 2 4 21 12 85 13
+82 2 2 4 21 82 86 81
+83 2 2 4 21 81 86 85
+84 2 2 4 21 83 86 82
+85 2 2 4 21 83 87 86
+86 2 2 4 21 84 88 83
+87 2 2 4 21 83 88 87
+88 2 2 4 21 28 27 84
+89 2 2 4 21 84 27 88
+90 2 2 4 21 85 89 13
+91 2 2 4 21 13 89 14
+92 2 2 4 21 86 90 85
+93 2 2 4 21 85 90 89
+94 2 2 4 21 87 91 86
+95 2 2 4 21 86 91 90
+96 2 2 4 21 88 92 87
+97 2 2 4 21 87 92 91
+98 2 2 4 21 27 26 88
+99 2 2 4 21 88 26 92
+100 2 2 4 21 89 1 14
+101 2 2 4 21 89 37 1
+102 2 2 4 21 90 37 89
+103 2 2 4 21 90 38 37
+104 2 2 4 21 91 38 90
+105 2 2 4 21 91 39 38
+106 2 2 4 21 92 39 91
+107 2 2 4 21 92 40 39
+108 2 2 4 21 26 40 92
+109 2 2 4 21 26 5 40
+110 2 2 2 22 18 94 28
+111 2 2 2 22 22 23 95
+112 2 2 2 22 24 95 23
+113 2 2 2 22 28 94 27
+114 2 2 2 22 5 96 25
+115 2 2 2 22 5 26 96
+116 2 2 2 22 4 18 28
+117 2 2 2 22 22 6 23
+118 2 2 2 22 24 98 95
+119 2 2 2 22 27 94 99
+120 2 2 2 22 18 19 94
+121 2 2 2 22 21 22 95
+122 2 2 2 22 93 94 100
+123 2 2 2 22 21 95 97
+124 2 2 2 22 25 98 24
+125 2 2 2 22 27 99 26
+126 2 2 2 22 93 99 94
+127 2 2 2 22 93 95 98
+128 2 2 2 22 93 97 95
+129 2 2 2 22 25 96 98
+130 2 2 2 22 26 99 96
+131 2 2 2 22 19 100 94
+132 2 2 2 22 20 21 97
+133 2 2 2 22 93 100 97
+134 2 2 2 22 19 20 100
+135 2 2 2 22 93 98 96
+136 2 2 2 22 93 96 99
+137 2 2 2 22 20 97 100
+138 4 2 11 1 7 12 42 105
+139 4 2 11 1 81 49 105 12
+140 4 2 11 1 12 49 105 7
+141 4 2 11 1 49 81 105 106
+142 4 2 11 1 82 50 106 49
+143 4 2 11 1 81 49 82 106
+144 4 2 11 1 50 82 106 83
+145 4 2 11 1 83 51 107 50
+146 4 2 11 1 106 50 83 107
+147 4 2 11 1 51 83 107 84
+148 4 2 11 1 84 52 108 51
+149 4 2 11 1 107 51 84 108
+150 4 2 11 1 52 84 108 28
+151 4 2 11 1 28 18 94 52
+152 4 2 11 1 108 52 28 94
+153 4 2 11 1 11 43 17 109
+154 4 2 11 1 109 65 77 11
+155 4 2 11 1 17 11 109 77
+156 4 2 11 1 65 109 77 78
+157 4 2 11 1 110 66 78 65
+158 4 2 11 1 109 65 110 78
+159 4 2 11 1 66 110 78 79
+160 4 2 11 1 111 67 79 66
+161 4 2 11 1 110 66 111 79
+162 4 2 11 1 67 111 79 80
+163 4 2 11 1 112 68 80 67
+164 4 2 11 1 111 67 112 80
+165 4 2 11 1 68 112 80 23
+166 4 2 11 1 95 22 23 68
+167 4 2 11 1 112 68 95 23
+168 4 2 11 1 16 17 43 109
+169 4 2 11 1 77 73 109 16
+170 4 2 11 1 17 16 77 109
+171 4 2 11 1 73 77 109 78
+172 4 2 11 1 78 74 110 73
+173 4 2 11 1 109 73 78 110
+174 4 2 11 1 74 78 110 79
+175 4 2 11 1 79 75 111 74
+176 4 2 11 1 110 74 79 111
+177 4 2 11 1 75 79 111 80
+178 4 2 11 1 80 76 112 75
+179 4 2 11 1 111 75 80 112
+180 4 2 11 1 76 80 112 23
+181 4 2 11 1 23 24 95 76
+182 4 2 11 1 112 76 23 95
+183 4 2 11 1 12 13 42 105
+184 4 2 11 1 85 81 105 12
+185 4 2 11 1 13 12 85 105
+186 4 2 11 1 81 85 105 106
+187 4 2 11 1 86 82 106 81
+188 4 2 11 1 85 81 86 106
+189 4 2 11 1 82 86 106 83
+190 4 2 11 1 87 83 107 106
+191 4 2 11 1 86 83 87 106
+192 4 2 11 1 83 87 107 88
+193 4 2 11 1 88 84 108 107
+194 4 2 11 1 83 84 88 107
+195 4 2 11 1 84 88 108 27
+196 4 2 11 1 27 28 94 108
+197 4 2 11 1 84 28 27 108
+198 4 2 11 1 1 15 44 113
+199 4 2 11 1 69 37 113 1
+200 4 2 11 1 15 1 69 113
+201 4 2 11 1 37 69 113 114
+202 4 2 11 1 70 38 114 37
+203 4 2 11 1 69 37 70 114
+204 4 2 11 1 38 70 114 115
+205 4 2 11 1 71 39 115 38
+206 4 2 11 1 70 38 71 115
+207 4 2 11 1 39 71 115 116
+208 4 2 11 1 72 40 116 39
+209 4 2 11 1 71 39 72 116
+210 4 2 11 1 40 72 116 96
+211 4 2 11 1 25 5 96 40
+212 4 2 11 1 72 40 25 96
+213 4 2 11 1 1 44 14 113
+214 4 2 11 1 113 37 89 1
+215 4 2 11 1 14 1 113 89
+216 4 2 11 1 37 113 89 90
+217 4 2 11 1 114 38 90 37
+218 4 2 11 1 113 37 114 90
+219 4 2 11 1 38 114 90 91
+220 4 2 11 1 115 39 91 38
+221 4 2 11 1 114 38 115 91
+222 4 2 11 1 39 115 91 92
+223 4 2 11 1 116 40 92 39
+224 4 2 11 1 115 39 116 92
+225 4 2 11 1 40 116 92 26
+226 4 2 11 1 96 5 26 40
+227 4 2 11 1 116 40 96 26
+228 4 2 11 1 3 12 7 49
+229 4 2 11 1 81 29 49 3
+230 4 2 11 1 12 3 81 49
+231 4 2 11 1 29 81 49 82
+232 4 2 11 1 82 30 50 29
+233 4 2 11 1 49 29 82 50
+234 4 2 11 1 30 82 50 83
+235 4 2 11 1 83 31 51 30
+236 4 2 11 1 50 30 83 51
+237 4 2 11 1 31 83 51 84
+238 4 2 11 1 84 32 52 31
+239 4 2 11 1 51 31 84 52
+240 4 2 11 1 32 84 52 28
+241 4 2 11 1 28 4 18 32
+242 4 2 11 1 52 32 28 18
+243 4 2 11 1 11 17 2 77
+244 4 2 11 1 77 65 33 2
+245 4 2 11 1 11 65 77 2
+246 4 2 11 1 65 77 33 78
+247 4 2 11 1 78 66 34 33
+248 4 2 11 1 65 66 78 33
+249 4 2 11 1 66 78 34 79
+250 4 2 11 1 79 67 35 66
+251 4 2 11 1 34 66 79 35
+252 4 2 11 1 67 79 35 80
+253 4 2 11 1 80 68 36 35
+254 4 2 11 1 67 68 80 35
+255 4 2 11 1 68 80 36 23
+256 4 2 11 1 23 22 6 36
+257 4 2 11 1 68 22 23 36
+258 4 2 11 1 16 43 46 121
+259 4 2 11 1 109 73 121 16
+260 4 2 11 1 43 16 109 121
+261 4 2 11 1 73 109 121 122
+262 4 2 11 1 110 74 122 73
+263 4 2 11 1 109 73 110 122
+264 4 2 11 1 74 110 122 123
+265 4 2 11 1 111 75 123 74
+266 4 2 11 1 110 74 111 123
+267 4 2 11 1 75 111 123 124
+268 4 2 11 1 112 76 124 75
+269 4 2 11 1 111 75 112 124
+270 4 2 11 1 76 112 124 98
+271 4 2 11 1 95 24 98 76
+272 4 2 11 1 112 76 95 98
+273 4 2 11 1 13 47 42 105
+274 4 2 11 1 125 85 105 13
+275 4 2 11 1 47 13 125 105
+276 4 2 11 1 85 125 105 106
+277 4 2 11 1 126 86 106 125
+278 4 2 11 1 125 86 106 85
+279 4 2 11 1 86 126 106 87
+280 4 2 11 1 127 87 107 126
+281 4 2 11 1 126 87 107 106
+282 4 2 11 1 87 127 107 128
+283 4 2 11 1 128 88 108 107
+284 4 2 11 1 87 88 128 107
+285 4 2 11 1 88 128 108 99
+286 4 2 11 1 99 27 94 108
+287 4 2 11 1 88 27 99 108
+288 4 2 11 1 7 42 8 105
+289 4 2 11 1 105 49 53 8
+290 4 2 11 1 7 49 105 8
+291 4 2 11 1 49 105 53 106
+292 4 2 11 1 106 50 54 49
+293 4 2 11 1 53 49 106 54
+294 4 2 11 1 50 106 54 107
+295 4 2 11 1 107 51 55 50
+296 4 2 11 1 54 50 107 55
+297 4 2 11 1 51 107 55 108
+298 4 2 11 1 108 52 56 51
+299 4 2 11 1 55 51 108 56
+300 4 2 11 1 52 108 56 94
+301 4 2 11 1 94 18 19 52
+302 4 2 11 1 56 52 94 19
+303 4 2 11 1 10 43 11 109
+304 4 2 11 1 109 61 65 11
+305 4 2 11 1 10 61 109 11
+306 4 2 11 1 61 109 65 110
+307 4 2 11 1 110 62 66 65
+308 4 2 11 1 61 62 110 65
+309 4 2 11 1 62 110 66 111
+310 4 2 11 1 111 63 67 66
+311 4 2 11 1 62 63 111 66
+312 4 2 11 1 63 111 67 112
+313 4 2 11 1 112 64 68 63
+314 4 2 11 1 67 63 112 68
+315 4 2 11 1 64 112 68 95
+316 4 2 11 1 95 21 22 64
+317 4 2 11 1 68 64 95 22
+318 4 2 11 1 41 48 42 105
+319 4 2 11 1 129 101 105 48
+320 4 2 11 1 48 101 105 41
+321 4 2 11 1 101 129 105 102
+322 4 2 11 1 130 102 106 105
+323 4 2 11 1 129 102 130 105
+324 4 2 11 1 102 130 106 103
+325 4 2 11 1 131 103 107 130
+326 4 2 11 1 130 103 107 106
+327 4 2 11 1 103 131 107 104
+328 4 2 11 1 132 104 108 107
+329 4 2 11 1 131 104 132 107
+330 4 2 11 1 104 132 108 93
+331 4 2 11 1 100 93 94 108
+332 4 2 11 1 132 93 100 108
+333 4 2 11 1 10 45 43 117
+334 4 2 11 1 117 61 109 10
+335 4 2 11 1 43 10 117 109
+336 4 2 11 1 61 117 109 118
+337 4 2 11 1 118 62 110 61
+338 4 2 11 1 109 61 118 110
+339 4 2 11 1 62 118 110 119
+340 4 2 11 1 119 63 111 62
+341 4 2 11 1 110 62 119 111
+342 4 2 11 1 63 119 111 120
+343 4 2 11 1 120 64 112 63
+344 4 2 11 1 111 63 120 112
+345 4 2 11 1 64 120 112 97
+346 4 2 11 1 97 21 95 64
+347 4 2 11 1 112 64 97 95
+348 4 2 11 1 15 16 46 121
+349 4 2 11 1 73 69 121 15
+350 4 2 11 1 16 15 73 121
+351 4 2 11 1 69 73 121 122
+352 4 2 11 1 74 70 122 69
+353 4 2 11 1 73 69 74 122
+354 4 2 11 1 70 74 122 123
+355 4 2 11 1 75 71 123 70
+356 4 2 11 1 74 70 75 123
+357 4 2 11 1 71 75 123 124
+358 4 2 11 1 76 72 124 71
+359 4 2 11 1 75 71 76 124
+360 4 2 11 1 72 76 124 98
+361 4 2 11 1 24 25 98 72
+362 4 2 11 1 76 72 24 98
+363 4 2 11 1 13 14 47 125
+364 4 2 11 1 89 85 125 13
+365 4 2 11 1 14 13 89 125
+366 4 2 11 1 85 89 125 90
+367 4 2 11 1 90 86 126 125
+368 4 2 11 1 85 86 90 125
+369 4 2 11 1 86 90 126 91
+370 4 2 11 1 91 87 127 126
+371 4 2 11 1 86 87 91 126
+372 4 2 11 1 87 91 127 128
+373 4 2 11 1 92 88 128 87
+374 4 2 11 1 91 87 92 128
+375 4 2 11 1 88 92 128 99
+376 4 2 11 1 26 27 99 88
+377 4 2 11 1 92 88 26 99
+378 4 2 11 1 41 42 47 105
+379 4 2 11 1 105 101 125 47
+380 4 2 11 1 41 101 105 47
+381 4 2 11 1 101 105 125 102
+382 4 2 11 1 106 102 126 125
+383 4 2 11 1 105 102 106 125
+384 4 2 11 1 102 106 126 103
+385 4 2 11 1 107 103 127 126
+386 4 2 11 1 106 103 107 126
+387 4 2 11 1 103 107 127 128
+388 4 2 11 1 108 104 128 107
+389 4 2 11 1 107 104 128 103
+390 4 2 11 1 104 108 128 99
+391 4 2 11 1 94 93 99 108
+392 4 2 11 1 108 93 99 104
+393 4 2 11 1 41 46 43 101
+394 4 2 11 1 121 101 109 43
+395 4 2 11 1 46 101 121 43
+396 4 2 11 1 101 121 109 122
+397 4 2 11 1 122 102 110 109
+398 4 2 11 1 101 102 122 109
+399 4 2 11 1 102 122 110 103
+400 4 2 11 1 123 103 111 110
+401 4 2 11 1 122 103 123 110
+402 4 2 11 1 103 123 111 124
+403 4 2 11 1 124 104 112 111
+404 4 2 11 1 103 104 124 111
+405 4 2 11 1 104 124 112 93
+406 4 2 11 1 98 93 95 112
+407 4 2 11 1 124 93 98 112
+408 4 2 11 1 41 43 45 101
+409 4 2 11 1 109 101 117 43
+410 4 2 11 1 43 101 117 45
+411 4 2 11 1 101 109 117 102
+412 4 2 11 1 110 102 118 109
+413 4 2 11 1 109 102 118 117
+414 4 2 11 1 102 110 118 103
+415 4 2 11 1 111 103 119 110
+416 4 2 11 1 110 103 119 118
+417 4 2 11 1 103 111 119 120
+418 4 2 11 1 112 104 120 111
+419 4 2 11 1 111 104 120 103
+420 4 2 11 1 104 112 120 93
+421 4 2 11 1 95 93 97 112
+422 4 2 11 1 112 93 97 120
+423 4 2 11 1 15 46 44 113
+424 4 2 11 1 121 69 113 15
+425 4 2 11 1 46 15 121 113
+426 4 2 11 1 69 121 113 122
+427 4 2 11 1 122 70 114 69
+428 4 2 11 1 113 69 122 114
+429 4 2 11 1 70 122 114 123
+430 4 2 11 1 123 71 115 70
+431 4 2 11 1 114 70 123 115
+432 4 2 11 1 71 123 115 124
+433 4 2 11 1 124 72 116 71
+434 4 2 11 1 115 71 124 116
+435 4 2 11 1 72 124 116 98
+436 4 2 11 1 98 25 96 72
+437 4 2 11 1 116 72 98 96
+438 4 2 11 1 14 44 47 113
+439 4 2 11 1 113 89 125 14
+440 4 2 11 1 47 14 113 125
+441 4 2 11 1 89 113 125 90
+442 4 2 11 1 114 90 126 113
+443 4 2 11 1 113 90 126 125
+444 4 2 11 1 90 114 126 91
+445 4 2 11 1 115 91 127 114
+446 4 2 11 1 114 91 127 126
+447 4 2 11 1 91 115 127 128
+448 4 2 11 1 116 92 128 115
+449 4 2 11 1 115 92 128 91
+450 4 2 11 1 92 116 128 99
+451 4 2 11 1 96 26 99 116
+452 4 2 11 1 116 26 99 92
+453 4 2 11 1 8 42 48 105
+454 4 2 11 1 105 53 129 8
+455 4 2 11 1 48 8 105 129
+456 4 2 11 1 53 105 129 130
+457 4 2 11 1 106 54 130 53
+458 4 2 11 1 105 53 106 130
+459 4 2 11 1 54 106 130 107
+460 4 2 11 1 107 55 131 54
+461 4 2 11 1 130 54 107 131
+462 4 2 11 1 55 107 131 132
+463 4 2 11 1 108 56 132 55
+464 4 2 11 1 107 55 108 132
+465 4 2 11 1 56 108 132 100
+466 4 2 11 1 94 19 100 56
+467 4 2 11 1 108 56 94 100
+468 4 2 11 1 9 45 10 117
+469 4 2 11 1 117 57 61 9
+470 4 2 11 1 10 9 117 61
+471 4 2 11 1 57 117 61 118
+472 4 2 11 1 118 58 62 57
+473 4 2 11 1 61 57 118 62
+474 4 2 11 1 58 118 62 119
+475 4 2 11 1 119 59 63 58
+476 4 2 11 1 62 58 119 63
+477 4 2 11 1 59 119 63 120
+478 4 2 11 1 120 60 64 63
+479 4 2 11 1 59 60 120 63
+480 4 2 11 1 60 120 64 97
+481 4 2 11 1 97 20 21 64
+482 4 2 11 1 60 20 97 64
+483 4 2 11 1 41 45 48 101
+484 4 2 11 1 117 101 129 48
+485 4 2 11 1 45 101 117 48
+486 4 2 11 1 101 117 129 102
+487 4 2 11 1 118 102 130 117
+488 4 2 11 1 117 102 130 129
+489 4 2 11 1 102 118 130 103
+490 4 2 11 1 119 103 131 130
+491 4 2 11 1 118 103 119 130
+492 4 2 11 1 103 119 131 120
+493 4 2 11 1 120 104 132 131
+494 4 2 11 1 103 104 120 131
+495 4 2 11 1 104 120 132 93
+496 4 2 11 1 97 93 100 120
+497 4 2 11 1 120 93 100 132
+498 4 2 11 1 8 48 9 129
+499 4 2 11 1 129 53 57 8
+500 4 2 11 1 9 8 129 57
+501 4 2 11 1 53 129 57 130
+502 4 2 11 1 130 54 58 53
+503 4 2 11 1 57 53 130 58
+504 4 2 11 1 54 130 58 131
+505 4 2 11 1 131 55 59 54
+506 4 2 11 1 58 54 131 59
+507 4 2 11 1 55 131 59 132
+508 4 2 11 1 132 56 60 55
+509 4 2 11 1 59 55 132 60
+510 4 2 11 1 56 132 60 100
+511 4 2 11 1 100 19 20 56
+512 4 2 11 1 60 56 100 20
+513 4 2 11 1 41 44 46 113
+514 4 2 11 1 113 101 121 46
+515 4 2 11 1 41 101 113 46
+516 4 2 11 1 101 113 121 122
+517 4 2 11 1 114 102 122 113
+518 4 2 11 1 113 102 122 101
+519 4 2 11 1 102 114 122 103
+520 4 2 11 1 115 103 123 114
+521 4 2 11 1 114 103 123 122
+522 4 2 11 1 103 115 123 124
+523 4 2 11 1 116 104 124 115
+524 4 2 11 1 115 104 124 103
+525 4 2 11 1 104 116 124 93
+526 4 2 11 1 96 93 98 116
+527 4 2 11 1 116 93 98 124
+528 4 2 11 1 41 47 44 113
+529 4 2 11 1 125 101 113 47
+530 4 2 11 1 47 101 113 41
+531 4 2 11 1 101 125 113 102
+532 4 2 11 1 126 102 114 113
+533 4 2 11 1 125 102 126 113
+534 4 2 11 1 102 126 114 103
+535 4 2 11 1 127 103 115 114
+536 4 2 11 1 126 103 127 114
+537 4 2 11 1 103 127 115 128
+538 4 2 11 1 128 104 116 115
+539 4 2 11 1 103 104 128 115
+540 4 2 11 1 104 128 116 99
+541 4 2 11 1 99 93 96 116
+542 4 2 11 1 104 93 99 116
+543 4 2 11 1 9 48 45 117
+544 4 2 11 1 129 57 117 9
+545 4 2 11 1 48 9 129 117
+546 4 2 11 1 57 129 117 130
+547 4 2 11 1 130 58 118 57
+548 4 2 11 1 117 57 130 118
+549 4 2 11 1 58 130 118 119
+550 4 2 11 1 131 59 119 58
+551 4 2 11 1 130 58 131 119
+552 4 2 11 1 59 131 119 120
+553 4 2 11 1 132 60 120 59
+554 4 2 11 1 131 59 132 120
+555 4 2 11 1 60 132 120 100
+556 4 2 11 1 100 20 97 60
+557 4 2 11 1 120 60 100 97
+$EndElements
diff --git a/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.py b/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.py
new file mode 100644
index 0000000000000000000000000000000000000000..045b0ea765b24f29b3a0237f20fa0b89edb6cfb9
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.py
@@ -0,0 +1,116 @@
+#coding-Utf-8-*-
+
+from gmshpy import *
+from dG3Dpy import*
+
+from math import*
+
+#script to launch PBC problem with a python script
+
+# material law
+lawnum = 11 # unique number of law
+
+E = 1410. #MPa
+nu = 0.4
+K = E/3./(1.-2.*nu) # Bulk mudulus
+mu =E/2./(1.+nu) # Shear mudulus
+rho = 7850e-9 # Bulk mass
+
+
+law1 = HyperViscoElasticDG3DMaterialLaw(lawnum,rho,E,nu,False,1e-8)
+law1.setStrainOrder(11)
+
+law1.setViscoelasticMethod(0)
+N = 3;
+law1.setViscoElasticNumberOfElement(N)
+law1.setViscoElasticData(1,700.,250)
+law1.setViscoElasticData(2,240,40)
+law1.setViscoElasticData(3,130,2)
+law1.setVolumeCorrection(0.6,1.,2.,0.6,1.,2.)
+
+# geometry
+meshfile="cylinder.msh" # name of mesh file
+
+fullDg = 0 #O = CG, 1 = DG
+space1 = 0 # function space (Lagrange=0)
+beta1 = 100
+# creation of part Domain
+nfield = 11 # number of the field (physical number of surface)
+
+myfield1 = dG3DDomain(1000,nfield,space1,lawnum,fullDg,3)
+#myfield1.matrixByPerturbation(1,1,1,1e-8)
+myfield1.stabilityParameters(beta1)
+
+# solver
+sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
+soltype =1 # StaticLinear=0 (default) StaticNonLinear=1
+nstep = 200 # number of step (used only if soltype=1)
+ftime =1.e3 # Final time (used only if soltype=1)
+tol=1.e-8 # relative tolerance for NR scheme (used only if soltype=1)
+nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
+
+
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+mysolver.loadModel(meshfile)
+mysolver.addDomain(myfield1)
+mysolver.addMaterialLaw(law1)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol)
+
+mysolver.displacementBC('Face',1,2,0.)
+mysolver.displacementBC('Face',3,1,0.)
+mysolver.displacementBC('Face',4,0,0.)
+
+mysolver.displacementBC('Face',1,0,0)
+mysolver.displacementBC('Face',1,1,0)
+
+fct = PiecewiseLinearFunction()
+fct.put(0.,0.)
+fct.put(0.25*ftime,1.)
+fct.put(0.5*ftime,0.)
+fct.put(0.75*ftime,1.)
+fct.put(ftime,0.)
+
+mysolver.displacementBC('Face',2,2,fct)
+
+# build view
+mysolver.internalPointBuildView("Strain_xx",IPField.STRAIN_XX, 1, 1);
+mysolver.internalPointBuildView("Strain_yy",IPField.STRAIN_YY, 1, 1);
+mysolver.internalPointBuildView("Strain_zz",IPField.STRAIN_ZZ, 1, 1);
+mysolver.internalPointBuildView("Strain_xy",IPField.STRAIN_XY, 1, 1);
+mysolver.internalPointBuildView("Strain_yz",IPField.STRAIN_YZ, 1, 1);
+mysolver.internalPointBuildView("Strain_xz",IPField.STRAIN_XZ, 1, 1);
+mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1);
+mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1);
+mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1);
+mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1);
+mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1);
+mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1);
+mysolver.internalPointBuildView("sig_VM",IPField.SVM, 1, 1);
+mysolver.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN, 1, 1);
+mysolver.internalPointBuildView("pression",IPField.PRESSION,1,1)
+mysolver.internalPointBuildView("plastic possion ratio",IPField.PLASTIC_POISSON_RATIO,1,1)
+
+mysolver.internalPointBuildView("Ee_ZZ",IPField.Ee_ZZ,1,1)
+mysolver.internalPointBuildView("Ee_YY",IPField.Ee_YY,1,1)
+mysolver.internalPointBuildView("Ee_XX",IPField.Ee_XX,1,1)
+
+mysolver.internalPointBuildView("Eve1_ZZ",IPField.Eve1_ZZ,1,1)
+mysolver.internalPointBuildView("Eve1_YY",IPField.Eve1_YY,1,1)
+mysolver.internalPointBuildView("Eve1_XX",IPField.Eve1_XX,1,1)
+
+#archiving
+mysolver.archivingForceOnPhysicalGroup('Face',1,0)
+mysolver.archivingForceOnPhysicalGroup('Face',1,1)
+mysolver.archivingForceOnPhysicalGroup('Face',1,2)
+mysolver.archivingNodeDisplacement(5,0,1)
+mysolver.archivingNodeDisplacement(5,1,1)
+mysolver.archivingNodeDisplacement(5,2,1)
+# solve
+mysolver.solve()
+
+check = TestCheck()
+check.equal(1.802052e+03,mysolver.getArchivedForceOnPhysicalGroup("Face", 1, 2),1.e-3)
diff --git a/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp b/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp
index 5b00e77a0ae576fb05547a649aa976c93112ad59..0b23fe65707b279860a9f34fb349b0d707ab0e54 100644
--- a/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp
+++ b/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp
@@ -93,7 +93,9 @@ void NonLocalDamagePowerYieldHyperDG3DMaterialLaw::setViscoElasticData_Shear(con
void NonLocalDamagePowerYieldHyperDG3DMaterialLaw::setViscoElasticData(const std::string filename){
_nldPowerYieldHyperlaw->setViscoElasticData(filename);
}
-
+void NonLocalDamagePowerYieldHyperDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc){
+ _nldPowerYieldHyperlaw->setVolumeCorrection(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
+};
NonLocalDamagePowerYieldHyperDG3DMaterialLaw::~NonLocalDamagePowerYieldHyperDG3DMaterialLaw(){
if (_nldPowerYieldHyperlaw) delete _nldPowerYieldHyperlaw;
_nldPowerYieldHyperlaw = NULL;
@@ -306,7 +308,9 @@ void LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setViscoElasticData_S
void LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setViscoElasticData(const std::string filename){
_localPowerYieldHyperlaw->setViscoElasticData(filename);
}
-
+void LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc){
+ _localPowerYieldHyperlaw->setVolumeCorrection(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
+};
void LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setFailureCriterion(const mlawHyperelasticFailureCriterionBase& fCr){
_localPowerYieldHyperlaw->setFailureCriterion(fCr);
};
@@ -517,7 +521,9 @@ void NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setViscoElasticDat
void NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setViscoElasticData(const std::string filename){
_nldPowerYieldHyperlaw->setViscoElasticData(filename);
}
-
+void NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc){
+ _nldPowerYieldHyperlaw->setVolumeCorrection(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
+};
void NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setFailureCriterion(const mlawHyperelasticFailureCriterionBase& fCr){
_nldPowerYieldHyperlaw->setFailureCriterion(fCr);
};
diff --git a/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.h b/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.h
index 83fbed40e66ef6f9a8051c88a0c79a18e70fd35e..6c1885a179882aa0c17221ad87628ca229b2bfed 100644
--- a/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.h
+++ b/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.h
@@ -44,6 +44,7 @@ class NonLocalDamagePowerYieldHyperDG3DMaterialLaw : public dG3DMaterialLaw
void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
void setViscoElasticData(const std::string filename);
+ void setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc);
#ifndef SWIG
virtual ~NonLocalDamagePowerYieldHyperDG3DMaterialLaw();
@@ -109,6 +110,7 @@ class LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure: public dG3DMaterialL
void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
void setViscoElasticData(const std::string filename);
+ void setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc);
void setFailureCriterion(const mlawHyperelasticFailureCriterionBase& fCr);
@@ -176,6 +178,7 @@ class NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure: public dG3DMateri
void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
void setViscoElasticData(const std::string filename);
+ void setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc);
void setFailureCriterion(const mlawHyperelasticFailureCriterionBase& fCr);
void setCritialDamage(const double D);
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 0eb550cf7ffa2c71c697315d8dd003095ce25680..2cad40102a6a1c0efddf19adc50e88132a963684 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -3894,6 +3894,7 @@ double PhenomenologicalSMPDG3DIPVariable::get(const int i) const
else if(i==IPField::epsAM) {return _TMIPSMP->getRefToEquivalentPlasticDefoAM();}
else if(i==IPField::epsR) {return _TMIPSMP->getRefToEquivalentPlasticDefoR();}
+ else if(i==IPField::epsG) {return _TMIPSMP->getRefToEquivalentPlasticDefoG();}
else if(i==IPField::Tcrys) {return _TMIPSMP->getConstRefToTcrys();}
else if(i==IPField::Tmelt) {return _TMIPSMP->getConstRefToTmelt();}
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 630001205bd71a311f51535f053d04280ca3a7f3..455a815151cb3087d1c95f274e5022c52d777f7e 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -2895,6 +2895,9 @@ void HyperViscoElasticDG3DMaterialLaw::setViscoElasticData_Shear(const int i, co
void HyperViscoElasticDG3DMaterialLaw::setViscoElasticData(const std::string filename){
_viscoLaw.setViscoElasticData(filename);
};
+void HyperViscoElasticDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc ,double _dc, double _thetadc, double _pidc){
+ _viscoLaw.setVolumeCorrection(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
+};
void HyperViscoElasticDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
@@ -3006,7 +3009,9 @@ void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setViscoElasticData_Shear
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setViscoElasticData(const std::string filename){
_viscoLaw.setViscoElasticData(filename);
};
-
+void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc ,double _dc, double _thetadc, double _pidc){
+ _viscoLaw.setVolumeCorrection(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
+};
void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setCompressionHardening(const J2IsotropicHardening& comp){
_viscoLaw.setCompressionHardening(comp);
};
@@ -6477,8 +6482,20 @@ void PhenomenologicalSMPDG3DMaterialLaw::setTemperatureFunction_Cp_rubbery(const
void PhenomenologicalSMPDG3DMaterialLaw::setTemperatureFunction_Cp_AM(const scalarFunction& Tfunc)
{_lawTMSMP->setTemperatureFunction_Cp_AM(Tfunc);};
-void PhenomenologicalSMPDG3DMaterialLaw::setElasticPotentialFunction(const elasticPotential& EPfunc)
-{_lawTMSMP->setElasticPotentialFunction(EPfunc);};
+
+void PhenomenologicalSMPDG3DMaterialLaw::setFunctionYieldCrystallinityG(const scalarFunction& Tfunc)
+{_lawTMSMP->setFunctionYieldCrystallinityG(Tfunc);};
+
+void PhenomenologicalSMPDG3DMaterialLaw::setFunctionYieldCrystallinityR(const scalarFunction& Tfunc)
+{_lawTMSMP->setFunctionYieldCrystallinityR(Tfunc);};
+
+void PhenomenologicalSMPDG3DMaterialLaw::setFunctionYieldCrystallinityAM(const scalarFunction& Tfunc)
+{_lawTMSMP->setFunctionYieldCrystallinityAM(Tfunc);};
+
+void PhenomenologicalSMPDG3DMaterialLaw::setElasticPotentialFunctionR(const elasticPotential& EPfuncR)
+{_lawTMSMP->setElasticPotentialFunctionR(EPfuncR);};
+void PhenomenologicalSMPDG3DMaterialLaw::setElasticPotentialFunctionAM(const elasticPotential& EPfuncAM)
+{_lawTMSMP->setElasticPotentialFunctionAM(EPfuncAM);};
void PhenomenologicalSMPDG3DMaterialLaw::setTcTmWcWm(const double Atc, const double Btm, const double Cwc, const double Dwm,
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 234321640135c11b24053a1687bf9c35e0bb1a8a..b8d6ce18ac9232428d61188f5565dbadbe62e0e5 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -737,7 +737,7 @@ class HyperViscoElasticDG3DMaterialLaw : public dG3DMaterialLaw{
void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
void setViscoElasticData(const std::string filename);
-
+ void setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc);
#ifndef SWIG
HyperViscoElasticDG3DMaterialLaw(const HyperViscoElasticDG3DMaterialLaw &source);
virtual ~HyperViscoElasticDG3DMaterialLaw(){}
@@ -788,7 +788,8 @@ class HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw : public dG3DMaterialLaw
void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
void setViscoElasticData(const std::string filename);
-
+ void setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc);
+
void setCompressionHardening(const J2IsotropicHardening& comp);
void setTractionHardening(const J2IsotropicHardening& trac);
void setKinematicHardening(const kinematicHardening& kin);
@@ -2189,6 +2190,18 @@ public:
virtual void setViscoElasticDataAM(const std::string filename)
{
_lawTMSMP->setViscoElasticDataAM( filename);
+ }
+ virtual void setVolumeCorrectionG(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc)
+ {
+ _lawTMSMP->setVolumeCorrectionG(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
+ }
+ virtual void setVolumeCorrectionR(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc)
+ {
+ _lawTMSMP->setVolumeCorrectionR(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
+ }
+ virtual void setVolumeCorrectionAM(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc)
+ {
+ _lawTMSMP->setVolumeCorrectionAM(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
}
virtual void setJ2IsotropicHardeningCompressionG(const J2IsotropicHardening& isoHard) const {
@@ -2290,8 +2303,12 @@ public:
virtual void setTemperatureFunction_Cp_glassy(const scalarFunction& Tfunc);
virtual void setTemperatureFunction_Cp_rubbery(const scalarFunction& Tfunc);
virtual void setTemperatureFunction_Cp_AM(const scalarFunction& Tfunc);
+ virtual void setFunctionYieldCrystallinityG(const scalarFunction& Tfunc);
+ virtual void setFunctionYieldCrystallinityR(const scalarFunction& Tfunc);
+ virtual void setFunctionYieldCrystallinityAM(const scalarFunction& Tfunc);
- virtual void setElasticPotentialFunction(const elasticPotential& EPfunc);
+ virtual void setElasticPotentialFunctionR(const elasticPotential& EPfuncR);
+ virtual void setElasticPotentialFunctionAM(const elasticPotential& EPfuncAM);
virtual void setTcTmWcWm(const double Atc, const double Btm, const double Cwc, const double Dwm,
const double alphaAtc, const double alphaBtm, const double alphaCwc, const double alphaDwm);