diff --git a/NonLinearSolver/internalPoints/ipHyperelastic.cpp b/NonLinearSolver/internalPoints/ipHyperelastic.cpp
index 605fe07dace5f0e46b2fe3e3fc60c3583aa53578..da6d2b9d101c42035c31fb1c9a0ec8191aed5147 100644
--- a/NonLinearSolver/internalPoints/ipHyperelastic.cpp
+++ b/NonLinearSolver/internalPoints/ipHyperelastic.cpp
@@ -13,7 +13,7 @@
#include "restartManager.h"
#include "STensorOperations.h"
-IPHyperViscoElastic::IPHyperViscoElastic(const int N):IPVariableMechanics(),_N(N),_elasticEnergy(0.),_Ee(0.),_kirchhoff(0.),_kirchhoff_vol(0.),_kirchhoff_dev(0.),
+IPHyperViscoElastic::IPHyperViscoElastic(const int N):IPVariableMechanics(),_N(N),_elasticEnergy(0.),_Ee(0.),_kirchhoff(0.),
_irreversibleEnergy(0.),_DirreversibleEnergyDF(0.), _viscousEnergyPart(0.), _dElasticEnergyPartdF(0.), _dViscousEnergyPartdF(0.), _elasticBulkPropertyScaleFactor(1.), _elasticShearPropertyScaleFactor(1.)
{
_A.clear();
@@ -25,7 +25,7 @@ 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), _kirchhoff_vol(src._kirchhoff_vol), _kirchhoff_dev(src._kirchhoff_dev),
+ _N(src._N),_A(src._A),_B(src._B),_elasticEnergy(src._elasticEnergy), _kirchhoff(src._kirchhoff),
_irreversibleEnergy(src._irreversibleEnergy),_DirreversibleEnergyDF(src._DirreversibleEnergyDF),
_viscousEnergyPart(src._viscousEnergyPart), _dElasticEnergyPartdF(src._dElasticEnergyPartdF),
_dViscousEnergyPartdF(src._dViscousEnergyPartdF), _elasticBulkPropertyScaleFactor(src._elasticBulkPropertyScaleFactor), _elasticShearPropertyScaleFactor(src._elasticShearPropertyScaleFactor)
@@ -40,8 +40,6 @@ 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;
@@ -65,8 +63,6 @@ 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);
diff --git a/NonLinearSolver/internalPoints/ipHyperelastic.h b/NonLinearSolver/internalPoints/ipHyperelastic.h
index 49302c54f245254d4ae05a127f234bb5e289563a..e15a2742f16b6b43abeafc0ac32f111fc67f4f2e 100644
--- a/NonLinearSolver/internalPoints/ipHyperelastic.h
+++ b/NonLinearSolver/internalPoints/ipHyperelastic.h
@@ -25,8 +25,6 @@ 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;
@@ -36,9 +34,8 @@ class IPHyperViscoElastic : public IPVariableMechanics{
protected:
//For energy sources
- double _viscousEnergyPart, _p_vol, _Dppr_vol_DCepr;
- STensor3 _dElasticEnergyPartdF, _dViscousEnergyPartdF, _K_dev_pr;
- STensor43 _DKpr_dev_DCepr;
+ double _viscousEnergyPart;
+ STensor3 _dElasticEnergyPartdF, _dViscousEnergyPartdF;
public:
IPHyperViscoElastic(const int N);
@@ -56,18 +53,11 @@ 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;};
@@ -135,9 +125,7 @@ class IPHyperViscoElastoPlastic : public IPHyperViscoElastic{
bool _dissipationBlocked;
bool _dissipationActive;
protected:
- STensor3 _dPlasticEnergyPartdF, _dptilde_VoldCepr, _K_dev, _dK_devdGamma;
- STensor43 _dK_devdCepr;
- double _ptilde_Vol, _X_vol, _dptilde_VoldGamma, _D_dev;
+ STensor3 _dPlasticEnergyPartdF;
public:
IPHyperViscoElastoPlastic(const J2IsotropicHardening* comp,
@@ -159,17 +147,6 @@ class IPHyperViscoElastoPlastic : public IPHyperViscoElastic{
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();
virtual bool inPostFailureStage() const {return _inPostFailureStage;};
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
index fb0d7c1f1a8f88b6467573ce167fbd6b8736bf66..7a5c788b19900f5c6ae52013e86a2c2057ca23e9 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
@@ -98,31 +98,19 @@ 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
+void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &B_d, STensor3 &dB_vddev) 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.));
+ A_v = getVolumeCorrection()*pow(exp(getXiVolumeCorrection()/3.*tr*tr)-1.,getZetaVolumeCorrection());
- 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));
- }
- }
- }
- }
+ dA_vdE = getVolumeCorrection()*getZetaVolumeCorrection()*pow(exp(getXiVolumeCorrection()/3.*tr*tr)-1.,getZetaVolumeCorrection()-1.)*exp(getXiVolumeCorrection()/3.*tr*tr) *getXiVolumeCorrection()*2./3.*tr;
+
+ B_d = getDevCorrection()*pow(exp(getThetaDevCorrection()*dev.dotprod())-1.,getPiDevCorrection());
+ STensorOperation::zero(dB_vddev);
+
+ dB_vddev=dev;
+ dB_vddev*=2.*getPiDevCorrection()*getDevCorrection()*pow(exp(getThetaDevCorrection()*dev.dotprod())-1.,getPiDevCorrection()-1.)* exp(getThetaDevCorrection()*dev.dotprod())*getThetaDevCorrection();
}
@@ -135,9 +123,9 @@ double mlawHyperViscoElastic::deformationEnergy(const IPHyperViscoElastic& q) co
static STensor3 devEe;
STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
- 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*pow(exp(_xivolCorrection*trEe)-1.,_zetavolCorrection))+getUpdatedShearModulus(&q)*STensorOperation::doubledot(devEe,devEe); //this is not correct
- //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]);
@@ -404,7 +392,7 @@ void mlawHyperViscoElastic::updateViscoElasticFlow(const IPHyperViscoElastic *q0
};
void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1,
- double& Ke, double& Ge) const{
+ double& Ke, double& Ge, double& A_v, double& B_d, double& dApr, STensor3& dB_vddev) const{
if ((_Ki.size() > 0) or (_Gi.size() > 0)){
static STensor3 DE, devDE;
static double trDE;
@@ -435,20 +423,21 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
Ke += _Ki[i]*ztak;
}
- 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
+ static STensor3 devK, devE;
+ STensorOperation::zero(devK);
+ STensorOperation::zero(devE);
+ double p, trEe;
+ STensorOperation::decomposeDevTr(Ee,devE,trEe);
+ evaluatePhiPCorrection(trEe, devE, A_v, dApr, B_d, dB_vddev);
+ devK = devE;
+ devK *= (2.*getUpdatedShearModulus(q1)*(1.+B_d)); // deviatoric part
- K_dev_pr = q1->getRefTo_K_dev_pr();
- p+=q1->getRefTo_p_vol();
+ p = trEe;
+ p *= getUpdatedBulkModulus(q1)*(1.+A_v); // pressure
for (int i=0; i<_Gi.size(); i++){
devK.daxpy(q1->_A[i],2.*_Gi[i]);
@@ -457,22 +446,11 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
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){
@@ -536,39 +514,26 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
}
}
else{
- 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();
-
+ static STensor3 devK, devE;
+ double p, trEe;
+ STensorOperation::zero(devK);
+ STensorOperation::zero(devE);
+
+ STensorOperation::decomposeDevTr(Ee,devE,trEe);
- 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();
+ STensorOperation::zero(dB_vddev);
+ evaluatePhiPCorrection(trEe, devE, A_v, dApr, B_d, dB_vddev);
+
+ devK = devE;
+ devK *= (2.*getUpdatedShearModulus(q1)*(1.+B_d)); // deviatoric part
+
+
+ p = trEe;
+ p *= getUpdatedBulkModulus(q1)*(1.+A_v); // pressure
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);
- }
- }
q1->_kirchhoff = devK;
q1->_kirchhoff(0,0) += p;
@@ -612,22 +577,26 @@ void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0,
STensorOperation::decomposeDevTr(E,devE,trE);
double Ke, Ge;
- viscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge);
+ double A_v_pr=0., B_d_pr = 0., dApr=0.;
+ static STensor3 dB_devpr;
+ STensorOperation::zero(dB_devpr);
+ viscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
static STensor3 secondPK;
STensorOperation::multSTensor3STensor43(corKir,dlnCdC,secondPK);
STensorOperation::multSTensor3(F,secondPK,P);
// first PK
-
+
+
+
q1->_elasticEnergy=deformationEnergy(*q1);
q1->getRefToViscousEnergyPart()=viscousEnergy(*q0,*q1)+q0->getConstRefToViscousEnergyPart();
if (stiff){
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();
+ isotropicHookTensor(Ke+getUpdatedBulkModulus(q1)*A_v_pr,Ge+getUpdatedShearModulus(q1)*B_d_pr,DcorKDE);
+
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
@@ -640,10 +609,7 @@ 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)+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);
+ DsecondPKdC(i,j,k,l) += 0.5*(DcorKDE(p,q,r,s)+2.*getUpdatedShearModulus(q1)*devE(p,q)*dB_devpr(r,s)+ getUpdatedBulkModulus(q1)*dApr*trE*_I(p,q)*_I(r,s))*dlnCdC(p,q,i,j)*dlnCdC(r,s,k,l);
}
}
}
@@ -652,7 +618,6 @@ void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0,
}
}
}
-
STensorOperation::zero(Tangent);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
@@ -813,291 +778,17 @@ 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
+void mlawPowerYieldHyper::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &B_d, STensor3 &dB_vddev) 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);
- }
- }
- }
- }
+ A_v = getVolumeCorrection()*pow(exp(1./3.*getXiVolumeCorrection()*tr*tr)-1.,getZetaVolumeCorrection());
+ dA_vdE =2./3.*getVolumeCorrection()*getXiVolumeCorrection()*getZetaVolumeCorrection()*pow(exp(1./3.*getXiVolumeCorrection()*tr*tr)-1.,getZetaVolumeCorrection()-1.)* exp(1./3.*getXiVolumeCorrection()*tr*tr)*tr;
- 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);
- }
- }
- }
- }
- }
- }
-
+ B_d = getDevCorrection()*pow(exp(getThetaDevCorrection()*dev.dotprod())-1.,getPiDevCorrection());
+ dB_vddev =2.*getDevCorrection()*getPiDevCorrection()*getThetaDevCorrection()*pow(exp(getThetaDevCorrection()*dev.dotprod())-1.,getPiDevCorrection()-1.)*exp(getThetaDevCorrection()*dev.dotprod())*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++){
- 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{
@@ -1182,25 +873,45 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
// update A, B
double Ge, Ke, trEepr;
static STensor3 devEepr;
- viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge);
+ STensorOperation::zero(devEepr);
+
+
+
+
- STensorOperation::decomposeDevTr(Ee,devEepr,trEepr);
+
+ double A_v_pr, B_d_pr, dApr;
+ static STensor3 dB_devpr;
+ STensorOperation::zero(dB_devpr);
+ viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
+
+ STensorOperation::decomposeDevTr(Ee,devEepr,trEepr);
+ //std::cout << "trEepr = "<<trEepr <<std::endl;
static STensor3 PhiPr, PhiPr_Vol, PhiPr_Dev;
PhiPr = q1->_kirchhoff;
PhiPr -= q1->_backsig;
- PhiPr_Dev = q1->_kirchhoff_dev;
- PhiPr_Vol = q1->_kirchhoff_vol;
+ STensorOperation::zero(PhiPr_Dev);
+ STensorOperation::zero(PhiPr_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;
+ static STensor3 devPhipr,devPhi, K_dev_pr, K_dev; // effective dev stress predictor
+ STensorOperation::zero(K_dev_pr);
+ STensorOperation::zero(K_dev);
+
+ double ptildepr,ptilde, ptildepr_Vol, ptilde_Vol;
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.;
+
+
+
+
+
+ ptildepr_Vol = getUpdatedBulkModulus(q1)*A_v_pr*trEepr;
+ K_dev_pr = devEepr;
+ K_dev_pr *= 2.*getUpdatedShearModulus(q1)*B_d_pr;
+
+
ptilde = ptildepr; // current effective pressure
ptilde_Vol = ptildepr_Vol; // current effective volumetric pressure
@@ -1239,8 +950,6 @@ 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.;
@@ -1251,16 +960,38 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
double dDgammaDGamma = 0.;
double Dgamma = 0.; // eqplastic strain
- double Dptilde_vol = 0.0;
- double dptilde_VoldGamma= 0.0;
- double Dptilde_volDGamma=0.0;
+ double Dptilde_vol = 0.;
+ double dptilde_VoldGamma= 0.;
double PhiEq_cor = 0.;
- static STensor3 DK_dev, devPhipr_Dev_cor, dK_devdGamma = q1->getRefTo_dK_devdGamma();
+
+
+ double A_v=0.;
+ double B_d = 0.;
+ double dAdtrEe=0.;
+ static STensor3 dBddevEe;
+ STensorOperation::zero(dBddevEe);
+
+ double trEe = 0.;
+ double J = 0.;
+ trEe = trEepr-2.*_b*Gamma*((ptildepr+Dptilde_vol)/(v));
+
+ static STensor3 DK_dev, devPhipr_Dev_cor, dK_devdGamma;
+ STensorOperation::zero(DK_dev);
+ static STensor3 J_dev, devEe;
+ STensorOperation::zero(devEe);
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 j=0; j<3; j++){
+ devEe(i,j) = devEepr(i,j)-3.*Gamma*(devPhipr(i,j)+DK_dev(i,j))/(u);
+ }
}
+ STensorOperation::zero(J_dev);
+
+
+
+
+ STensorOperation::zero(devPhipr_Dev_cor);
+ STensorOperation::zero(dK_devdGamma);
+
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
@@ -1290,24 +1021,20 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
_viscosity->get(q1->_epspbarre,eta,Deta);
double etaOverDt = eta/this->getTimeStep();
-
- 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);
+ double dAdGamma = 1./(2.*A)*(-(72.*Gt*PhiEq*PhiEq)/(u)-(16.*_b*_b*_b*Kt*ptilde*ptilde/(3*v))+((8.*_b*_b*ptilde)/(3*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));
+ dAdGamma+= 1./(2.*A)*(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)*(ptildepr+Dptilde_vol)/v -Da(0); //OK
+ dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*ptilde -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)*(ptildepr+Dptilde_vol)*2.*_b*Kt/(v*v)-(a(1)*dptilde_VoldGamma)/(v);
+ DfDGamma = dfdDgamma*dDgammaDGamma - (_n*a(2)*6.*Gt)*pow(PhiEq,_n)/u + a(1)*ptilde*2.*_b*Kt/(v)-((a(1)/v)*dptilde_VoldGamma);
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));
@@ -1328,34 +1055,132 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
u = 1.+6.*Gt*Gamma;
v = 1.+2.*_b*Kt*Gamma;
-
+
+ double dtrEedptilde = -2.*_b*Gamma/v;
+ J=getUpdatedBulkModulus(q1);
+ int ite_vol = 0;
+ int maxite_vol = 100; // maximal number of iters
+ double dJdptilde_Vol =0.;
+ double Np_vol;
+
+ while (1){
+
+ evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, B_d, dBddevEe);
+ J=getUpdatedBulkModulus(q1)*A_v*trEe-ptilde_Vol;
+ dJdptilde_Vol = getUpdatedBulkModulus(q1)*( dAdtrEe*dtrEedptilde*trEe+A_v*dtrEedptilde)-1.;
- 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 = ptilde_Vol-J/dJdptilde_Vol;
+
+ Dptilde_vol = ptilde_Vol-ptildepr_Vol;
+ trEe = trEepr-2.*_b*Gamma*(ptildepr+Dptilde_vol)/(v);
+ J=getUpdatedBulkModulus(q1)*A_v*trEe-ptilde_Vol;
+
- 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);
- }
+ if ((fabs(J)/getUpdatedBulkModulus(q1) <_tol and ite_vol >2) or ite_vol > maxite_vol)
+ {
+ dptilde_VoldGamma = -(2.*_b*getUpdatedBulkModulus(q1)/(v*v)*(ptildepr+Dptilde_vol)*(dAdtrEe*trEe+A_v))/(1.+2.*_b*Gamma*getUpdatedBulkModulus(q1)/(v)*(dAdtrEe*trEe+A_v));
+ if (ite_vol > maxite_vol)
+ Msg::Error("No convergence for phi_p_vol in mlawPowerYieldHyper nonAssociatedFlow Maxwell iter = %d, J = %e!!",ite,J);
+ break;
+ }
+ ite_vol++;
}
- 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);
+
+
+ //################################### DEVIATORIC TERMS
+
+ static STensor43 ddevEedK_dev;
+ STensorOperation::zero(ddevEedK_dev);
+
+
+
+ double J_dev_tol=getUpdatedShearModulus(q1);
+ int ite_dev = 0;
+ int maxite_dev = 100; // maximal number of iters
+
+ static STensor43 dJ_dev, dJ_dev_inv;
+ static STensor3 dTemp, NK_dev;
+ STensorOperation::zero(dJ_dev);
+ STensorOperation::zero(dJ_dev_inv);
+ STensorOperation::zero(dTemp);
+ STensorOperation::zero(NK_dev);
+
+ while (1){
+
+ evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, B_d, dBddevEe);
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ J_dev(i,j)=2.*getUpdatedShearModulus(q1)*B_d*devEe(i,j)-K_dev(i,j) ;
+ }
+ }
+
+ ddevEedK_dev = _I4;
+ ddevEedK_dev *= -3.*Gamma/u;
+ STensorOperation::zero(dTemp);
+ STensorOperation::multSTensor3STensor43(dBddevEe,ddevEedK_dev,dTemp);
+ 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)*devEe(i,j)*dTemp(k,l)+ 2.*getUpdatedShearModulus(q1)*B_d*ddevEedK_dev(i,j,k,l)-_I4(i,j,k,l);
+ }
+ }
+ }
+ }
+
+ STensorOperation::inverseSTensor43(dJ_dev, dJ_dev_inv);
+ STensorOperation::zero(NK_dev);
+ STensorOperation::multSTensor43STensor3(dJ_dev_inv,J_dev,NK_dev);
+
+
+ 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); // tau_c
+ DK_dev(i,j)=K_dev(i,j)-K_dev_pr(i,j); //Delta Tau_c =tau_c - tau_c_pr
+ devEe(i,j) = devEepr(i,j)-3.*Gamma*(devPhipr(i,j)+DK_dev(i,j))/(u);
+ devPhipr_Dev_cor(i,j) = devPhipr(i,j)+DK_dev(i,j); //Phipr+Delta tau_c
+ J_dev(i,j)=2.*getUpdatedShearModulus(q1)*B_d*devEe(i,j)-K_dev(i,j) ;
+ }
+ }
+ J_dev_tol=J_dev.norm0();
+
+ if ((fabs(J_dev_tol)/getUpdatedShearModulus(q1)<_tol and ite_dev >2) or ite_dev > maxite_dev)
+ {
+ //#############dK_devdGamma
+ double M_d = 0.;
+ static STensor3 T_d, U_d;
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ T_d(i,j) = -3./(u*u)*(devPhipr_Dev_cor(i,j));
+ M_d+=dBddevEe(i,j)*T_d(j,i);
+ }
+ }
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ U_d(i,j)=2.*getUpdatedShearModulus(q1)*devEe(i,j)*M_d+2.*getUpdatedShearModulus(q1)*B_d*T_d(i,j);
+ }
+ }
+ STensorOperation::multSTensor43STensor3(dJ_dev_inv,U_d,dK_devdGamma, -1.);
+
+ //#############dK_devdGamma
+ 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;
+
}
+ ite_dev++;
+
}
+
+
+
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));
+ ptilde = (ptildepr+Dptilde_vol)/v;
+ A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
Dgamma = kk*Gamma*A;
//Msg::Error("it = %d, u=%e, v=%e, Dgamma=%e",ite,u,v,Dgamma);
@@ -1364,7 +1189,8 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
hardening(q0,q1);
getYieldCoefficients(q1,a);
- f = a(2)*pow(PhiEq,_n) - (a(1)*(ptildepr+Dptilde_vol)/(v)+a(0));
+ f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
+
double viscoTerm = etaOverDt*Gamma;
if (Gamma>0 and etaOverDt>0) f-= pow(viscoTerm,_p);
@@ -1372,7 +1198,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
//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;
+ //std::cout << "NR F: f = "<<f <<" Gamma = "<<Gamma<<" dfdDgamma = "<<dfdDgamma<< " dDgammaDGamma = "<< dDgammaDGamma << " dptilde_VoldGamma = " << dptilde_VoldGamma<< " Dptilde_vol"<< Dptilde_vol<< " ptilde = "<<ptilde << " 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){
@@ -1381,16 +1207,13 @@ 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+Dptilde_vol)/v;
-
- // update normal
devN = devPhi;
devN *= 3.;
trN = 2.*_b*ptilde;
@@ -1444,7 +1267,6 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
// second Piola Kirchhoff stress
static STensor3 S;
STensorOperation::multSTensor3STensor43(KS,DlnDCe,S);
-
for(int i=0; i<3; i++)
for(int j=0; j<3; j++){
P(i,j) = 0.;
@@ -1478,15 +1300,47 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
if (stiff){
static STensor3 DpprDCepr;
static STensor43 DdevKprDCepr;
- STensorOperation::multSTensor3STensor43(_I,DlnDCepr,DpprDCepr);
- static STensor3 devEe;
- double trEe;
-
- STensorOperation::decomposeDevTr(Ee,devEe,trEe);
-
- DpprDCepr*= (0.5*Ke)+ q1->getRefTo_Dppr_vol_DCepr();
+
+
+
+ STensorOperation::multSTensor3STensor43(_I,DlnDCepr,DpprDCepr);
+ DpprDCepr*= (0.5*Ke)+ 0.5*getUpdatedBulkModulus(q1)*(dApr*trEepr+A_v_pr);
+
+
+
+ static STensor43 DKpr_dev_DCepr, DKpr_dev_DCepr_Temp, Temp1, Temp2, Temp4;
+ STensorOperation::zero(DKpr_dev_DCepr);
+ STensorOperation::zero(Temp1);
+ STensorOperation::zero(Temp2);
+ STensorOperation::zero(Temp4);
+
STensorOperation::multSTensor43(_Idev,DlnDCepr,DdevKprDCepr);
- DdevKprDCepr*= Ge;
+
+ DdevKprDCepr*= Ge;//+getUpdatedShearModulus(q1)*B_d_pr;
+
+ 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++){
+ Temp1(i,j,k,l) += getUpdatedShearModulus(q1)*(devEepr(i,j)*dB_devpr(k,l)+B_d_pr*_Idev(i,j,k,l));
+ }
+ }
+ }
+ }
+
+ STensorOperation::multSTensor43(Temp1,_Idev,Temp2);
+ STensorOperation::multSTensor43(Temp2,DlnDCepr,Temp4);
+
+ 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++){
+ DdevKprDCepr(i,j,k,l) += Temp4(i,j,k,l);
+ }
+ }
+ }
+ }
+
static STensor3 DpDCepr;
static STensor43 DdevKDCepr;
DpDCepr = DpprDCepr;
@@ -1498,36 +1352,83 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
static STensor43 DdevNDCepr;
static STensor3 DGDCepr;
+
if (Gamma >0){
// plastic
static STensor3 dAdCepr, dfDCepr, dptilde_VoldCepr, dgamadCepr, DGDCepr, dptildeOldpr_VoldCepr;
- static STensor43 dK_devdCepr;
+ static STensor43 dK_devdCepr, DdevKOldprDCepr;
+ STensorOperation::zero(dAdCepr);
+ STensorOperation::zero(dfDCepr);
+ STensorOperation::zero(dptilde_VoldCepr);
+ STensorOperation::zero(dgamadCepr);
+ STensorOperation::zero(DGDCepr);
+ STensorOperation::zero(dptildeOldpr_VoldCepr);
+ STensorOperation::zero(dK_devdCepr);
+ STensorOperation::zero(DdevKOldprDCepr);
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ dptilde_VoldCepr(i,j) = ( getUpdatedBulkModulus(q1)*(trEe*dAdtrEe+A_v)*(_I(i,j)-2.*_b*Gamma/v*Ke*_I(i,j)) )/(1.+2.*_b*Gamma*getUpdatedBulkModulus(q1)/v*(trEe*dAdtrEe+A_v));
+ }
+ }
- 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);
+ //########################################dK_devdEepr
+ static STensor43 H_d;
+ static STensor43 S_dev, S_dev_inv;
+ STensorOperation::zero(H_d);
+ STensorOperation::zero(S_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++){
+ for (int m=0; m<3; m++){
+ for (int n=0; n<3; n++){
+ H_d(i,j,k,l) += 2.*getUpdatedShearModulus(q1)*(1.-6.*Ge*Gamma/u)*(devEe(i,j)*dBddevEe(m,n)*_Idev(n,m,k,l));
+ S_dev(i,j,k,l) += 6.*getUpdatedShearModulus(q1)*Gamma/u*(devEe(i,j)*dBddevEe(m,n)*_Idev(n,m,k,l));
+ }
+ }
+ H_d(i,j,k,l) += 2.*getUpdatedShearModulus(q1)*(1.-6.*Ge*Gamma/u)*B_d*_Idev(i,j,k,l);
+ S_dev(i,j,k,l) += (1.+6.*getUpdatedShearModulus(q1)*Gamma/u*B_d)*_I4(i,j,k,l);
+ }
+ }
+ }
+ }
+
+ STensorOperation::inverseSTensor43(S_dev, S_dev_inv);
+ STensorOperation::zero(dK_devdCepr);
+ STensorOperation::multSTensor43(S_dev_inv,H_d,dK_devdCepr);
+
+ //########################################dK_devdEepr
+
+ // d/dEpr-> d/dCpr
+ STensorOperation::multSTensor43(dK_devdCepr, DlnDCepr, dK_devdCepr);
+ dK_devdCepr*=0.5;
+ STensorOperation::multSTensor3STensor43(dptilde_VoldCepr,DlnDCepr,dptilde_VoldCepr);
+ dptilde_VoldCepr*=0.5;
+
STensorOperation::multSTensor3STensor43(_I,DlnDCepr,dptildeOldpr_VoldCepr);
dptildeOldpr_VoldCepr*=0.5*Ke;
+ STensorOperation::multSTensor43(_Idev,DlnDCepr,DdevKOldprDCepr);
+ DdevKOldprDCepr*= Ge;
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+Dptilde_vol)/(A*3.*v*v))*(dptildeOldpr_VoldCepr(i,j)+dptilde_VoldCepr(i,j));
+ dAdCepr(i,j) = (4.*_b*_b*ptilde/(A*3.*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_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));
+ dAdCepr(i,j) += (9./(A*u*u))*devPhipr_Dev_cor(k,l)*(DdevKOldprDCepr(k,l,i,j)+dK_devdCepr(k,l,i,j));
+ dfDCepr(i,j) += fact*devPhipr_Dev_cor(k,l)*(DdevKOldprDCepr(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;
@@ -1536,7 +1437,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
- DgamaDCepr(i,j) = kk*Gamma*dAdCepr(i,j)+ kk*dDgammaDGamma*DGDCepr(i,j);
+ DgamaDCepr(i,j) = kk*Gamma*dAdCepr(i,j)+ dDgammaDGamma*DGDCepr(i,j); // remove the extra kk
}
}
@@ -1550,7 +1451,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
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) = (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);
+ DdevNDCepr(i,j,k,l) = (3./u)*(DdevKOldprDCepr(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;
@@ -1570,17 +1471,17 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
EprFp0(i,j,k,l) += dexpAdA(i,s,k,l)*Fp0(s,j);
}
}
-
- DpDCepr=dptildeOldpr_VoldCepr;
STensorOperation::multSTensor43(EprFp0,temp1,dFpDCepr);
+
+ DpDCepr = dptildeOldpr_VoldCepr;
// 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))-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))-dK_devdCepr(i,j,k,l)-DGDCepr(k,l)*dK_devdGamma(i,j);
+ DdevKDCepr(i,j,k,l) = DdevKOldprDCepr(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);
}
}
}
@@ -1700,7 +1601,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
Tangent(i,j,k,l) = 0.;
for (int m=0; m<3; m++){
for (int n=0; n<3; n++){
- Tangent(i,j,k,l) += dFedF(i,m,k,l)*S(m,n)*Fpinv(j,n);
+ Tangent(i,j,k,l) +=dFedF(i,m,k,l)*S(m,n)*Fpinv(j,n);
Tangent(i,j,k,l) += Fe(i,m)*dSdF(m,n,k,l)*Fpinv(j,n);
Tangent(i,j,k,l) += Fe(i,m)*S(m,n)*DinvFpDF(j,n,k,l);
}
@@ -1818,7 +1719,13 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c
// update A, B
double Ge, Ke;
- viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge);
+
+
+
+ double A_v_pr=0., B_d_pr = 0., dApr=0.;
+ static STensor3 dB_devpr;
+ STensorOperation::zero(dB_devpr);
+ viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
static STensor3 devKpr; // dev corotational kirchoff stress predictor
static double ppr; // pressure predictor
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.h b/NonLinearSolver/materialLaw/mlawHyperelastic.h
index de3347fe7ecc34f6911d97cc4fa2869786d5ae42..d56585037ce48fd50bab08d436c1b211c7c8dc42 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.h
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.h
@@ -59,14 +59,16 @@ class mlawHyperViscoElastic : public materialLaw{
void updateViscoElasticFlow(const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1, double& Ke, double& Ge) const;
void viscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1,
- double& Ke, double& Ge) const;
+ double& Ke, double& Ge, double& A_v, double& B_d, double& dAprd, STensor3& dB_vddev) const;
void isotropicHookTensor(const double K, const double G, STensor43& L) const;
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;
+ void evaluatePhiPCorrection(double trEe, const STensor3 &devEpr, double &A_v, double &dAprdEepr, double &B_d, STensor3 &dB_vddevEe) const;
+ //void evaluateA_dA(double trEe, double &A_v, double &dAprdEepr) const;
+
#endif // SWIG
public:
@@ -82,12 +84,12 @@ class mlawHyperViscoElastic : public materialLaw{
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 getXiVolumeCorrection() const {if (_xivolCorrection<1.e-5) return 1.; else return _xivolCorrection;};
+ virtual double getZetaVolumeCorrection() const {if (_zetavolCorrection<1.e-5) return 1.; 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;};
+ virtual double getThetaDevCorrection() const {if (_thetadevCorrection<1.e-5) return 1.; else return _thetadevCorrection;};
+ virtual double getPiDevCorrection() const {if (_pidevCorrection<1.e-5) return 1.; else return _pidevCorrection;};
#ifndef SWIG
mlawHyperViscoElastic(const mlawHyperViscoElastic& src);
@@ -178,7 +180,8 @@ class mlawPowerYieldHyper : public mlawHyperViscoElastic{
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;
+ void evaluatePhiPCorrection(double trEe, const STensor3 &devEpr, double &A_v, double &dAprdEepr, double &B_d, STensor3 &dB_vddevEe) const;
+
protected:
diff --git a/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.py b/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.py
index 045b0ea765b24f29b3a0237f20fa0b89edb6cfb9..0d14fb42b8d56b525e98b443e18bbc1f01bed586 100644
--- a/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.py
+++ b/dG3D/benchmarks/powerYieldViscoElastWithCorrection/cylinder.py
@@ -26,7 +26,7 @@ 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.)
+law1.setVolumeCorrection(8.,6.,2.,8.,6.,2.)
# geometry
meshfile="cylinder.msh" # name of mesh file
diff --git a/dG3D/benchmarks/smpPhenomenological_Mech/SMP_mech.py b/dG3D/benchmarks/smpPhenomenological_Mech/SMP_mech.py
index 2af42de66eba302f9dfdbd41aa33293c750efe51..5432043aa301c2162b3a4fc29461a5e3206e878b 100644
--- a/dG3D/benchmarks/smpPhenomenological_Mech/SMP_mech.py
+++ b/dG3D/benchmarks/smpPhenomenological_Mech/SMP_mech.py
@@ -54,8 +54,8 @@ cfAM = 0.0
zAM = 0.05 #0.685 #0.685
tau_y0G = 1000.e6
-tau_y0R = 6.e6 #0.6e6
-tau_y0AM= 6.e6
+tau_y0R = 0.6e6 #0.6e6
+tau_y0AM= 0.6e6
HR1 = 2.e5
HR2 = 1.
HR3 = 1
@@ -68,11 +68,11 @@ HG3 = 1.
###
##extra branch
-corJAM=corJR=corJG=8. #8.
+corJAM=corJR=corJG=8.
xicorJAM=xicorJR=xicorJG=6. #6.
-zetacorJAM=zetacorJR=zetacorJG=250. #250. #Has to be > 2.
-dcorJAM=dcorJR=dcorJG=10. #10.
-thetacorJAM=thetacorJG=thetacorJR=0.225 #0.225 #0.225
+zetacorJAM=zetacorJR=zetacorJG=4. #250. #Has to be > 2.
+dcorJAM=dcorJR=dcorJG=10.
+thetacorJAM=thetacorJG=thetacorJR=6. #0.225 #0.225 #0.225
picorJAM=picorJG=picorJR=2. #2.
#1.e6
########
@@ -182,8 +182,8 @@ law1.setYieldPowerFactorAM(3.5);
#law1.setElasticPotentialFunctionAM(EPFunc_AM)
#law1.setElasticPotentialFunctionR(EPFunc_R)
-#law1.setVolumeCorrectionAM(corJAM, xicorJAM, zetacorJAM, dcorJAM, thetacorJAM, picorJAM)
-#law1.setVolumeCorrectionR(corJR, xicorJR, zetacorJR, dcorJR, thetacorJR, picorJR)
+law1.setVolumeCorrectionAM(corJAM, xicorJAM, zetacorJAM, dcorJAM, thetacorJAM, picorJAM)
+law1.setVolumeCorrectionR(corJR, xicorJR, zetacorJR, dcorJR, thetacorJR, picorJR)
#law1.setVolumeCorrectionG(corJG, xicorJG, zetacorJG, dcorJG, thetacorJG, picorJG)
law1.setTcTmWcWm(Atc, Btm, Cwc, Dwm, alphatc,alphatm, alphawc, alphawm)
@@ -229,7 +229,7 @@ mysolver.addDomain(myfield1)
mysolver.addMaterialLaw(law1)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
-mysolver.snlData(nstep,ftime/1.,tol,tol)
+mysolver.snlData(nstep,ftime/1.,tol,tol/100.)
mysolver.snlManageTimeStep(50, 3, 2, 10)
mysolver.stepBetweenArchiving(nstepArch)
@@ -435,7 +435,7 @@ mysolver.solve()
# ===Test==========================================================================================
check= TestCheck()
-check.equal(-2.859195e-05,mysolver.getArchivedNodalValue(7,0,mysolver.displacement),1e-5)
+check.equal(-2.590210e-5,mysolver.getArchivedNodalValue(7,0,mysolver.displacement),1e-5)
@@ -610,7 +610,7 @@ mysolver.solve()
# ===Test==========================================================================================
check2 = TestCheck()
-check2.equal(0.,mysolver.getArchivedNodalValue(7,0,mysolver.displacement),1e-1)
-check2.equal(0.,mysolver.getArchivedNodalValue(7,2,mysolver.displacement),1e-1)
+check2.equal(-1.433048e-05,mysolver.getArchivedNodalValue(7,0,mysolver.displacement),1e-1)
+check2.equal(6.559040e-05,mysolver.getArchivedNodalValue(7,2,mysolver.displacement),1e-1)
diff --git a/dG3D/benchmarks/smpPhenomenological_Thermo/SMP_thermo.py b/dG3D/benchmarks/smpPhenomenological_Thermo/SMP_thermo.py
index 746e1513956c0594573bacf5e4f71a3b2df32b42..753f9110ba9aa456de8e29d64b9a2ea9b974cc22 100644
--- a/dG3D/benchmarks/smpPhenomenological_Thermo/SMP_thermo.py
+++ b/dG3D/benchmarks/smpPhenomenological_Thermo/SMP_thermo.py
@@ -54,8 +54,8 @@ cfAM = 0.0
zAM = 0.05 #0.685 #0.685
tau_y0G = 1000.e6
-tau_y0R = 6.e6 #0.6e6
-tau_y0AM= 6.e6
+tau_y0R = 0.6e6
+tau_y0AM= 0.6e6
HR1 = 2.e5
HR2 = 1.
HR3 = 1
@@ -68,11 +68,11 @@ HG3 = 1.
###
##extra branch
-corJAM=corJR=corJG=8. #8.
+corJAM=corJR=corJG=8.
xicorJAM=xicorJR=xicorJG=6. #6.
-zetacorJAM=zetacorJR=zetacorJG=250. #250. #Has to be > 2.
-dcorJAM=dcorJR=dcorJG=10. #10.
-thetacorJAM=thetacorJG=thetacorJR=0.225 #0.225 #0.225
+zetacorJAM=zetacorJR=zetacorJG=4. #250. #Has to be > 2.
+dcorJAM=dcorJR=dcorJG=10.
+thetacorJAM=thetacorJG=thetacorJR=6. #0.225 #0.225 #0.225
picorJAM=picorJG=picorJR=2. #2.
#1.e6
########
@@ -432,7 +432,7 @@ mysolver.solve()
check = TestCheck()
-check.equal(-1.563127e+02,mysolver.getArchivedForceOnPhysicalGroup("Face", 5678, 3),1.e-6)
+check.equal(-1.50000e+02,mysolver.getArchivedForceOnPhysicalGroup("Face", 5678, 3),1.e-6)
check.equal(3.452332e2,mysolver.getArchivedNodalValue(1,3,mysolver.displacement),1.e-6)
check.equal(0.,mysolver.getArchivedNodalValue(7,0,mysolver.displacement),1.e-6)