Select Git revision
nodalBasis.h
-
Christophe Geuzaine authored
This reverts commit 211fa521, reversing changes made to 8ba56d73.
Christophe Geuzaine authoredThis reverts commit 211fa521, reversing changes made to 8ba56d73.
mlawHyperelastic.cpp 78.77 KiB
//
// C++ Interface: material law
//
// Author: <Van Dung Nguyen>, (C) 2014
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "mlawHyperelastic.h"
#include "STensorOperations.h"
#include "nonLinearMechSolver.h"
void mlawHyperViscoElastic::setStrainOrder(const int order){
_order = order;
Msg::Info("order = %d is used to approximate log and exp");
};
void mlawHyperViscoElastic::setViscoelasticMethod(const int method){
_viscoMethod = (viscoelasticType)method;
if (_viscoMethod == Maxwell)
Msg::Info("generalized maxwell model is used for viscoelasticity");
else if (_viscoMethod == KelvinVoight)
Msg::Info("generalized Voigt-Kelvin model is used for viscoelasticity");
else
Msg::Error("this method has not been implemented");
};
void mlawHyperViscoElastic::setViscoElasticNumberOfElement(const int N){
_N = N;
Msg::Info("Numer of Spring/Dashpot for viscoelastic model: %d",_N);
_Ki.clear(); _ki.clear();
_Gi.clear(); _gi.clear();
_Ki.resize(_N,0.);
_ki.resize(_N,0.);
_Gi.resize(_N,0.);
_gi.resize(_N,0.);
};
void mlawHyperViscoElastic::setViscoElasticData(const int i, const double Ei, const double taui){
if (i> _N or i<1)
Msg::Error("This setting is invalid %d > %d",i,_N);
else{
double _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;
_ki[i-1] = taui;
_Gi[i-1] = GG;
_gi[i-1] = taui;
Msg::Info("setting: element=%d Ki= %e ki = %e, Gi=%e, gi=%e",i-1,KK,taui,GG,taui);
}
};
void mlawHyperViscoElastic::setViscoElasticData_Bulk(const int i, const double Ki, const double ki){
if (i> _N or i<1)
Msg::Error("This setting is invalid %d > %d",i,_N);
else{
_Ki[i-1] = Ki;
_ki[i-1] = ki;
// Msg::Info("setting: element=%d Ki= %e ki = %e, ",i-1,Ki,ki);
}
};
void mlawHyperViscoElastic::setViscoElasticData_Shear(const int i, const double Gi, const double gi){
if (i> _N or i<1)
Msg::Error("This setting is invalid %d > %d",i,_N);
else{
_Gi[i-1] = Gi;
_gi[i-1] = gi;
// Msg::Info("setting: element=%d: Gi=%e, gi=%e",i-1,Gi,gi);
}
};
void mlawHyperViscoElastic::setViscoElasticData(const std::string filename){
FILE* file = fopen(filename.c_str(),"r");
if (file == NULL) Msg::Error("file: %s does not exist",filename.c_str());
_Ki.clear(); _ki.clear();
_Gi.clear(); _gi.clear();
_N = 0;
Msg::Info("reading viscoelastic input data");
while (!feof(file)){
int ok = fscanf(file,"%d",&_N);
Msg::Info("Numer of Maxwell elements: %d",_N);
double KK, k, GG, g;
for (int i=0; i<_N; i++){
ok = fscanf(file,"%lf %lf %lf %lf",&KK,&k,&GG,&g);
_Ki.push_back(KK);
_ki.push_back(k);
_Gi.push_back(GG);
_gi.push_back(g);
Msg::Info("Maxwell element %d: K[%d] = %e, k[%d] = %e, G[%d] = %e, g[%d] = %e",
i,i,KK,i,k,i,GG,i,g);
}
}
fclose(file);
};
void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &B_d, STensor3 &dB_vddev) const
{
int method =1;
if(method ==0)
{
A_v = getVolumeCorrection()*pow(exp(getXiVolumeCorrection()/3.*tr*tr)-1.,getZetaVolumeCorrection());
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();
}
else
{
A_v = getVolumeCorrection()*(tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection())+tanh(getZetaVolumeCorrection()));
dA_vdE = getVolumeCorrection()*getXiVolumeCorrection()*2./3.*(1.-tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection())*tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection()));
B_d = getDevCorrection()*(tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection())+tanh(getPiDevCorrection()));
STensorOperation::zero(dB_vddev);
dB_vddev=dev;
dB_vddev*=2.*getPiDevCorrection()*getDevCorrection()*(1.-tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection())*tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection()) );
}
}
double mlawHyperViscoElastic::deformationEnergy(const IPHyperViscoElastic& q) const
{
double Psy = 0.;
if (_viscoMethod == Maxwell)
{
double trEe;
static STensor3 devEe;
STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
double Av, dAv,B_d;
static STensor3 dB_d;
evaluatePhiPCorrection(trEe, devEe, Av, dAv, B_d, dB_d);
Psy = getUpdatedBulkModulus(&q)*Av*(0.5*trEe*trEe)+getUpdatedShearModulus(&q)*B_d*STensorOperation::doubledot(devEe,devEe); //this is not correct we should do an integral
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]);
}
}
else if (_viscoMethod == KelvinVoight)
{
// need to be complete
Psy = 0.5*STensorOperation::doubledot(q._Ee,q._kirchhoff);
for (int i=0; i<_N; i++){
static STensor3 Ev;
q.getViscoElasticStrain(i+1,Ev);
Psy += 0.5*STensorOperation::doubledot(Ev,q._kirchhoff);
}
}
else{
Msg::Error("visco elastic method %d has not been defined");
}
return Psy;;
}
void mlawHyperViscoElastic::dDeformationEnergydF(const IPHyperViscoElastic& q, const STensor3 &P, STensor3 &dElasticEnergydF) const
{
STensorOperation::zero(dElasticEnergydF);
dElasticEnergydF=P;
}
double mlawHyperViscoElastic::viscousEnergy(const IPHyperViscoElastic& q0,const IPHyperViscoElastic& q) const
{
double DViscous = 0.;
if (_viscoMethod == Maxwell)
{
double trEe, trEe0;
static STensor3 devEe, devEe0;
STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
STensorOperation::decomposeDevTr(q0._Ee,devEe0,trEe0);
for (int i=0; i<_N; i++){
static STensor3 dqdev;
dqdev = devEe;
dqdev -= q._A[i];
dqdev -= devEe0;
dqdev += q0._A[i];
dqdev*=2.*_Gi[i];
double dtrq= trEe-q._B[i]-trEe0+q0._B[i];
dtrq*=_Ki[i]*3.;
DViscous+= STensorOperation::doubledot(q._A[i],dqdev);
DViscous+= q._B[i]*dtrq/3.;
}
}
//else
//{
// Msg::Error("visco elastic method %d has not been defined");
//}
return DViscous;
}
void mlawHyperViscoElastic::dViscousEnergydF(const IPHyperViscoElastic& q0, const IPHyperViscoElastic& q, const STensor43 &dlnCdC, const STensor3 &Fe, STensor3 &dViscousEnergydF) const
{
static STensor3 dDeltaViscousdEve;
STensorOperation::zero(dDeltaViscousdEve);
STensorOperation::zero(dViscousEnergydF);
if (_viscoMethod == Maxwell)
{
double trEe, trEe0;
static STensor3 devEe, devEe0;
STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
STensorOperation::decomposeDevTr(q0._Ee,devEe0,trEe0);
double dt=this->getTimeStep();
for (int i=0; i<_N; i++)
{
double dtg = dt/(_gi[i]);
double ratiog=exp(-dtg/2.);
double dtk = dt/(_ki[i]);
double ratiok=exp(-dtk/2.);
static STensor3 dqdev, GA2;
dqdev = devEe;
dqdev -= q._A[i];
dqdev -= devEe0;
dqdev += q0._A[i];
dqdev*=2.*_Gi[i];
double dtrq= trEe-q._B[i]-trEe0+q0._B[i];
dtrq*=_Ki[i]*3.;
GA2=q._A[i];
GA2*=2.*_Gi[i];
dDeltaViscousdEve=dqdev;
dDeltaViscousdEve-=GA2;
dDeltaViscousdEve*=ratiog;
dDeltaViscousdEve+=GA2;
double ppart=(dtrq/3.-q._B[i]*_Ki[i])*ratiok+q._B[i]*_Ki[i];
for (int k=0; k<3; k++){
dDeltaViscousdEve(k,k)+=ppart;
}
}
//HERE WE NEED TO GET WITH RESPECT TO DF -> dlog
}
//else
//{
// Msg::Error("visco elastic method %d has not been defined");
//}
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 N=0.; N<3.;N++)
{
dViscousEnergydF(i,J)+=0.5*dDeltaViscousdEve(K,L)*(dlnCdC(K,L,J,N)*Fe(i,N)+dlnCdC(K,L,N,J)*Fe(i,N));
}
}
}
}
}
}
mlawHyperViscoElastic::mlawHyperViscoElastic(const int num,const double E,const double nu, const double rho,
const bool matrixbyPerturbation, const double pert):
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.), _xivolCorrection(1.), _zetavolCorrection(1.), _devCorrection(0.), _thetadevCorrection(1.), _pidevCorrection(1.)
{
_Idev = _I4;
STensor3 mIon3(-1./3);
STensor43 mIIon3;
tensprod(_I,mIon3, mIIon3);
_Idev += mIIon3;
double _lambda = (E*nu)/(1.+nu)/(1.-2.*nu);
_mu = 0.5*E/(1.+nu);
_K = E/3./(1.-2.*nu);
Cel=0.;
Cel(0,0,0,0) = _lambda + 2.*_mu;
Cel(1,1,0,0) = _lambda;
Cel(2,2,0,0) = _lambda;
Cel(0,0,1,1) = _lambda;
Cel(1,1,1,1) = _lambda + 2.*_mu;
Cel(2,2,1,1) = _lambda;
Cel(0,0,2,2) = _lambda;
Cel(1,1,2,2) = _lambda;
Cel(2,2,2,2) = _lambda + 2.*_mu;
Cel(1,0,1,0) = _mu;
Cel(2,0,2,0) = _mu;
Cel(0,1,0,1) = _mu;
Cel(2,1,2,1) = _mu;
Cel(0,2,0,2) = _mu;
Cel(1,2,1,2) = _mu;
Cel(0,1,1,0) = _mu;
Cel(0,2,2,0) = _mu;
Cel(1,0,0,1) = _mu;
Cel(1,2,2,1) = _mu;
Cel(2,0,0,2) = _mu;
Cel(2,1,1,2) = _mu;
};
mlawHyperViscoElastic::mlawHyperViscoElastic(const mlawHyperViscoElastic& src): materialLaw(src),
_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), _volCorrection(src._volCorrection), _xivolCorrection(src._xivolCorrection),
_zetavolCorrection(src._zetavolCorrection), _devCorrection(src._devCorrection), _thetadevCorrection(src._thetadevCorrection),
_pidevCorrection(src._pidevCorrection)
{
};
mlawHyperViscoElastic& mlawHyperViscoElastic::operator=(const materialLaw &source)
{
materialLaw::operator=(source);
const mlawHyperViscoElastic* src =dynamic_cast<const mlawHyperViscoElastic*>(&source);
if(src != NULL)
{
_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;
_zetavolCorrection=src->_zetavolCorrection; _devCorrection=src->_devCorrection;
_thetadevCorrection=src->_thetadevCorrection; _pidevCorrection=src->_pidevCorrection;
}
return *this;
}
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);
}
void mlawHyperViscoElastic::updateViscoElasticFlow(const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1, double& Ke, double & Ge) const{
if ((_Ki.size() > 0) or (_Gi.size() > 0)){
double dt = this->getTimeStep();
if (_viscoMethod == Maxwell){
static STensor3 DE, devDE;
static double trDE;
DE = q1->_Ee;
DE -= q0->_Ee;
STensorOperation::decomposeDevTr(DE,devDE,trDE);
// maxwell
Ge = getUpdatedShearModulus(q1);
Ke = getUpdatedBulkModulus(q1);
for (int i=0; i<_Gi.size(); i++){
double dtg = dt/_gi[i];
double expmdtg = exp(-dtg);
double ztag = exp(-dtg/2.);
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
q1->_A[i](k,l) = expmdtg*q0->_A[i](k,l) + ztag*devDE(k,l);
}
}
Ge += _Gi[i]*ztag;
}
for (int i=0; i<_Ki.size(); i++){
double dtk = dt/_ki[i];
double expmdtk = exp(-dtk);
double ztak = exp(-dtk/2.);
q1->_B[i] = q0->_B[i]*expmdtk +ztak*trDE;
Ke += _Ki[i]*ztak;
}
}
else if (_viscoMethod == KelvinVoight){
static STensor3 DK, devDK;
static double trDK;
DK = q1->_kirchhoff;
DK -= q0->_kirchhoff;
STensorOperation::decomposeDevTr(DK,devDK,trDK);
/*update internal variable from stress increment*/
for (int i=0; i< _Gi.size(); i++)
{
STensorOperation::zero(q1->_A[i]);
double dtg = dt/(_gi[i]);
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
q1->_A[i](k,l) += exp(-dtg)*q0->_A[i](k,l) + exp(-dtg/2.)*devDK(k,l)/(2.*_Gi[i]);
}
}
}
for (int i=0; i< _Ki.size(); i++){
q1->_B[i] = 0.;
double dtk = dt/(_ki[i]);
q1->_B[i] += exp(-dtk)*q0->_B[i] + exp(-dtk/2.)*trDK/(3.*_Ki[i]);
}
}
else{
Msg::Error("visco elastic method %d has not been defined");
}
}
};
void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1,
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;
DE = Ee;
DE -= Ee0;
STensorOperation::decomposeDevTr(DE,devDE,trDE);
double dt = this->getTimeStep();
if (_viscoMethod == Maxwell){
Ge = getUpdatedShearModulus(q1);
Ke = getUpdatedBulkModulus(q1);
for (int i=0; i<_Gi.size(); i++){
double dtg = dt/_gi[i];
double expmdtg = exp(-dtg);
double ztag = exp(-dtg/2.);
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
q1->_A[i](k,l) = expmdtg*q0->_A[i](k,l) + ztag*devDE(k,l);
}
}
Ge += _Gi[i]*ztag;
}
for (int i=0; i<_Ki.size(); i++){
double dtk = dt/_ki[i];
double expmdtk = exp(-dtk);
double ztak = exp(-dtk/2.);
q1->_B[i] = q0->_B[i]*expmdtk +ztak*trDE;
Ke += _Ki[i]*ztak;
}
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
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]);
}
for (int i=0; i<_Ki.size(); i++){
p += q1->_B[i]*_Ki[i];
}
q1->_kirchhoff = devK;
q1->_kirchhoff(0,0) += p;
q1->_kirchhoff(1,1) += p;
q1->_kirchhoff(2,2) += p;
}
else if (_viscoMethod == KelvinVoight){
double invGe = 1./getUpdatedShearModulus(q1);
STensor3 D(0.);
for (int i=0; i<_Gi.size(); i++){
double dtg = dt/(_gi[i]);
invGe += (1.-exp(-dtg/2.))/_Gi[i];
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
D(k,l) += q0->_A[i](k,l)*(exp(-dtg)-1.);
}
}
}
Ge = 1./invGe;
double invKe = 1./_K;
double V= 0.;
for (int i=0; i<_Ki.size(); i++){
double dtk = dt/(_ki[i]);
invKe += (1.-exp(-dtk/2))/_Ki[i];
V += q0->_B[i]*(exp(-dtk)-1.);
}
Ke = 1./invKe;
// stress increment
static STensor3 DdevK;
DdevK = devDE; // dev corotational kirchoff stress predictor
DdevK += D;
DdevK *= (2.*Ge);
double Dp = Ke*(trDE+ V); // pressure predictor
/*update internal variable from stress increment*/
for (int i=0; i< _Gi.size(); i++){
STensorOperation::zero(q1->_A[i]);
double dtg = dt/(_gi[i]);
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
q1->_A[i](k,l) += exp(-dtg)*q0->_A[i](k,l) + exp(-dtg/2.)*DdevK(k,l)/(2.*_Gi[i]);
}
}
}
for (int i=0; i< _Ki.size(); i++){
q1->_B[i] = 0.;
double dtk = dt/(_ki[i]);
q1->_B[i] += exp(-dtk)*q0->_B[i] + exp(-dtk/2.)*Dp/(_Ki[i]);
}
// corotational Kirchhoff stress tenor from previous and increment
STensor3& corK = q1->_kirchhoff;
corK = q0->_kirchhoff;
corK += DdevK;
corK(0,0) += Dp;
corK(1,1) += Dp;
corK(2,2) += Dp;
}
else{
Msg::Error("visco elastic method %d has not been defined");
}
}
else{
static STensor3 devK, devE;
double p, trEe;
STensorOperation::zero(devK);
STensorOperation::zero(devE);
STensorOperation::decomposeDevTr(Ee,devE,trEe);
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);
q1->_kirchhoff = devK;
q1->_kirchhoff(0,0) += p;
q1->_kirchhoff(1,1) += p;
q1->_kirchhoff(2,2) += p;
}
};
void mlawHyperViscoElastic::isotropicHookTensor(const double K, const double G, STensor43& L) const{
double lambda = K - 2.*G/3.;
static STensor3 I(1.);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
L(i,j,k,l) = lambda*I(i,j)*I(k,l)+ G*(I(i,k)*I(j,l)+I(i,l)*I(j,k));
}
}
}
}
};
void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0, const STensor3& F, STensor3&P, const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1,
const bool stiff, STensor43& Tangent) const{
static STensor3 C;
static STensor43 dlnCdC;
static STensor63 ddlnCddC;
STensor3& E = q1->getRefToElasticStrain();
STensorOperation::multSTensor3FirstTranspose(F,F,C);
if (_order == 1){
STensorOperation::logSTensor3(C,_order,E,&dlnCdC); // as ddlogCddC = 0
}
else{
STensorOperation::logSTensor3(C,_order,E,&dlnCdC,&ddlnCddC);
}
E *= 0.5; // strain
static STensor3 devE;
double trE;
STensorOperation::decomposeDevTr(E,devE,trE);
double 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+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++){
for (int l=0; l<3; l++){
DsecondPKdC(i,j,k,l) = 0.;
for (int p=0; p<3; p++){
for (int q=0; q<3; q++){
if (_order != 1){
DsecondPKdC(i,j,k,l) += corKir(p,q)*ddlnCddC(p,q,i,j,k,l);
}
for (int r=0; r<3; r++){
for (int s=0; s<3; s++){
DsecondPKdC(i,j,k,l) += 0.5*(DcorKDE(p,q,r,s)+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);
}
}
}
}
}
}
}
}
STensorOperation::zero(Tangent);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int p=0; p<3; p++){
Tangent(i,j,i,p) += secondPK(p,j);
for (int q=0; q<3; q++){
for (int s=0; s<3; s++){
for (int k=0; k<3; k++){
Tangent(i,j,p,q) += 2.*F(i,k)*DsecondPKdC(k,j,q,s)*F(p,s);
}
}
}
}
}
}
dDeformationEnergydF(*q1, P, q1->getRefToDElasticEnergyPartdF());
dViscousEnergydF(*q0,*q1,dlnCdC, F, q1->getRefToDViscousEnergyPartdF());
};
};
void mlawHyperViscoElastic::constitutive(const STensor3& F0,
const STensor3& F, STensor3&P,
const IPVariable *q0i,
IPVariable *q1i,
STensor43& Tangent,
const bool stiff,
STensor43* elasticTangent,
const bool dTangent,
STensor63* dCalgdeps) const{
const IPHyperViscoElastic *q0=dynamic_cast<const IPHyperViscoElastic *>(q0i);
IPHyperViscoElastic *q1 = dynamic_cast<IPHyperViscoElastic *>(q1i);
if(elasticTangent!=NULL) Msg::Error("mlawHyperViscoElastic elasticTangent not defined");
if (!_tangentByPerturbation){
this->predictorCorrector_ViscoElastic(F0,F,P,q0,q1,stiff,Tangent);
}
else{
this->predictorCorrector_ViscoElastic(F0,F,P,q0,q1,false,Tangent);
if (stiff){
static STensor3 Fplus, Pplus;
static STensor43 Ttemp;
static IPHyperViscoElastic qtemp(*q1);
for (int i= 0; i<3; i++){
for (int j=0; j<3; j++){
Fplus = F;
Fplus(i,j) += _perturbationfactor;
this->predictorCorrector_ViscoElastic(F0,Fplus,Pplus,q0,&qtemp,false,Ttemp);
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
Tangent(k,l,i,j) = (Pplus(k,l)-P(k,l))/_perturbationfactor;
}
}
q1->getRefToDElasticEnergyPartdF()(i,j)=(qtemp.getConstRefToElasticEnergyPart()-q1->getConstRefToElasticEnergyPart())/_perturbationfactor;
q1->getRefToDViscousEnergyPartdF()(i,j)=(qtemp.getConstRefToViscousEnergyPart()-q1->getConstRefToViscousEnergyPart())/_perturbationfactor;
}
}
}
}
};
void mlawPowerYieldHyper::setCompressionHardening(const J2IsotropicHardening& comp){
if (_compression) delete _compression;
_compression = comp.clone();
};
void mlawPowerYieldHyper::setTractionHardening(const J2IsotropicHardening& trac){
if (_traction) delete _traction;
_traction = trac.clone();
};
void mlawPowerYieldHyper::setKinematicHardening(const kinematicHardening& kin){
if (_kinematic) delete _kinematic;
_kinematic = kin.clone();
};
void mlawPowerYieldHyper::setViscosityEffect(const viscosityLaw& v, const double p){
_p = p;
if (_viscosity) delete _viscosity;
_viscosity = v.clone();
};
void mlawPowerYieldHyper::hardening(const IPHyperViscoElastoPlastic* q0, IPHyperViscoElastoPlastic* q) const{
//Msg::Error("epspCompression = %e, epspTRaction = %e, epspShear = %e",q->_epspCompression,q->_epspTraction,q->_epspShear);
if (_compression != NULL && q->_ipCompression != NULL){
_compression->hardening(q0->_epspCompression, *q0->_ipCompression, q->_epspCompression,*q->_ipCompression);
}
if (_traction!= NULL && q->_ipTraction != NULL){
_traction->hardening(q0->_epspTraction,*q0->_ipTraction, q->_epspTraction,*q->_ipTraction);
}
if (_kinematic!= NULL && q->_ipKinematic != NULL)
_kinematic->hardening(q0->_epspbarre,*q0->_ipKinematic,q->_epspbarre,*q->_ipKinematic);
};
void mlawPowerYieldHyper::tangent_full_perturbation(
STensor43 &T_,
STensor43& dFedF,
STensor43& dFpdF,
const STensor3 &P,
const STensor3 &F,
const IPHyperViscoElastoPlastic* q0,
IPHyperViscoElastoPlastic* q1
) const{
static STensor43 tmpSTensor43;
static STensor3 Fplus, Pplus;
static IPHyperViscoElastoPlastic q11(*q0);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
Fplus = F;
Fplus(i,j)+=_perturbationfactor;
this->predictorCorrector(Fplus,q0,&q11,Pplus,false,tmpSTensor43,tmpSTensor43,tmpSTensor43,NULL);
q1->_DgammaDF(i,j) = (q11._epspbarre - q1->_epspbarre)/_perturbationfactor;
q1->_DirreversibleEnergyDF(i,j) = (q11._irreversibleEnergy - q1->_irreversibleEnergy)/_perturbationfactor;
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
T_(k,l,i,j) = (Pplus(k,l)-P(k,l))/_perturbationfactor;
dFpdF(k,l,i,j) = (q11._Fp(k,l)-q1->_Fp(k,l))/_perturbationfactor;
dFedF(k,l,i,j) = (q11._Fe(k,l)-q1->_Fe(k,l))/_perturbationfactor;
}
}
q1->getRefToDElasticEnergyPartdF()(i,j)=(q11.getConstRefToElasticEnergyPart()-q1->getConstRefToElasticEnergyPart())/_perturbationfactor;
q1->getRefToDViscousEnergyPartdF()(i,j)=(q11.getConstRefToViscousEnergyPart()-q1->getConstRefToViscousEnergyPart())/_perturbationfactor;
q1->getRefToDPlasticEnergyPartdF()(i,j)=(q11.getConstRefToPlasticEnergyPart()-q1->getConstRefToPlasticEnergyPart())/_perturbationfactor;
}
}
};
void mlawPowerYieldHyper::constitutive(
const STensor3& F0, // initial deformation gradient (input @ time n)
const STensor3& Fn, // updated deformation gradient (input @ time n+1)
STensor3 &P, // updated 1st Piola-Kirchhoff stress tensor (output)
const IPVariable *q0i, // array of initial internal variable
IPVariable *q1i, // updated array of internal variable (in ipvcur on output),
STensor43 &Tangent, // constitutive tangents (output)
const bool stiff, // if true compute the tangents
STensor43* elasticTangent,
const bool dTangent,
STensor63* dCalgdeps) const{
const IPHyperViscoElastoPlastic *q0=dynamic_cast<const IPHyperViscoElastoPlastic *>(q0i);
IPHyperViscoElastoPlastic *q1 = dynamic_cast<IPHyperViscoElastoPlastic *>(q1i);
static STensor43 dFedF, dFpdF;
if (_tangentByPerturbation){
this->predictorCorrector(Fn,q0,q1,P,false,Tangent,dFedF,dFpdF,elasticTangent);
if (stiff)
this->tangent_full_perturbation(Tangent,dFedF,dFpdF,P,Fn,q0,q1);
}
else{
this->predictorCorrector(Fn,q0,q1,P,stiff,Tangent,dFedF,dFpdF,elasticTangent);
}
};
void mlawPowerYieldHyper::updateEqPlasticDeformation(IPHyperViscoElastoPlastic *q1, const IPHyperViscoElastoPlastic *q0,
const double& nup, const double& Dgamma) const{
q1->_epspbarre = q0->_epspbarre+ Dgamma;
q1->_epspCompression = q0->_epspCompression+ Dgamma;
q1->_epspTraction = q0->_epspTraction+ Dgamma;
double k = 1./(sqrt(1.+2.*nup*nup));
q1->_epspShear = q0->_epspShear+ Dgamma/(k*sqrt(2.));
};
void mlawPowerYieldHyper::getYieldCoefficients(const IPHyperViscoElastoPlastic *q, fullVector<double>& coeffs) const{
double sigc = q->_ipCompression->getR();
double sigt = q->_ipTraction->getR();
double m = sigt/sigc;
//Msg::Error("m = %e, _n = %e",m,_n);
//Msg::Error("sigc = %e, sigt = %e",sigc,sigt);
coeffs.resize(3);
coeffs(2) = pow(sigc,-_n);
coeffs(1) = 3.*(pow(m,_n)-1.)/(m+1.)/sigc;
coeffs(0) = (pow(m,_n)+m)/(m+1);
};
void mlawPowerYieldHyper::getYieldCoefficientDerivatives(const IPHyperViscoElastoPlastic *q, const double& nup, fullVector<double>& Dcoeffs) const{
double sigc(0.), Hc(0.);
sigc = q->_ipCompression->getR();
Hc = q->_ipCompression->getDR();
double sigt(0.), Ht(0.);
sigt = q->_ipTraction->getR();
Ht = q->_ipTraction->getDR();
Dcoeffs.resize(3);
double m = sigt/sigc;
double Dm = (Ht*sigc- sigt*Hc)/(sigc*sigc);
double Da1Dm = 3./sigc*(_n*pow(m,_n-1.)/(m+1.) - (pow(m,_n)-1.)/(m+1.)/(m+1.));
Dcoeffs(2) = -_n*pow(sigc,-_n-1.)*Hc;
Dcoeffs(1) = Da1Dm*Dm -3.*(pow(m,_n)-1.)/(m+1.)/(sigc*sigc)*Hc;
Dcoeffs(0) = ((_n*pow(m,_n-1)+1.)/(m+1) - (pow(m,_n)+m)/(m+1.)/(m+1.))*Dm;
};
void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1,
STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF) const{
/* compute elastic predictor */
STensor3& Fp1 = q1->_Fp;
const STensor3& Fp0 = q0->_Fp;
Fp1 = Fp0; // plastic deformation tensor
q1->_epspbarre = q0->_epspbarre; // plastic equivalent strain
q1->_epspCompression = q0->_epspCompression;
q1->_epspTraction = q0->_epspTraction;
q1->_epspShear = q0->_epspShear;
q1->_backsig = q0->_backsig; // backstress tensor
q1->_DgammaDt = 0.;
static STensor3 Fpinv, Ce, Fepr;
STensorOperation::inverseSTensor3(Fp1,Fpinv);
STensorOperation::multSTensor3(F,Fpinv,Fepr);
STensorOperation::multSTensor3FirstTranspose(Fepr,Fepr,Ce);
static STensor3 invFp0; // plastic predictor
invFp0= Fpinv;
STensor3& Fe = q1->_Fe;
Fe = Fepr;
static STensor43 DlnDCepr, DlnDCe;
static STensor63 DDlnDDCe;
static STensor43 dexpAdA; // estimation of dexpA/dA
STensor3& Ee = q1->_Ee;
STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCepr,&DDlnDDCe);
Ee *= 0.5;
DlnDCe = DlnDCepr;
// update A, B
double Ge, Ke, trEepr;
static STensor3 devEepr;
STensorOperation::zero(devEepr);
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;
STensorOperation::zero(PhiPr_Dev);
STensorOperation::zero(PhiPr_Vol);
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.;
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
devPhi =devPhipr;
K_dev = K_dev_pr;
double PhiEqpr2 = 1.5*devPhipr.dotprod();
double PhiEqpr = sqrt(PhiEqpr2);
// plastic poisson ratio
q1->_nup = (9.-2.*_b)/(18.+2.*_b);
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
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
this->getYieldCoefficients(q1,a);
double Hb =0.;
if (q1->_ipKinematic != NULL)
Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
double Gt= Ge + kk*Hb/2.;
double Kt = Ke + kk*Hb/3.;
//a.print("a init");
static STensor3 devN; // dev part of yield normal
double trN = 0.; // trace part of yield normal
static STensor3 N; // yield normal
STensorOperation::zero(devN);
STensorOperation::zero(N);
double f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
double DfDGamma = 0.;
double dfdDgamma = 0.;
double u = 1.;
double v = 1.;
double A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
double dDgammaDGamma = 0.;
double Dgamma = 0.; // eqplastic strain
double Dptilde_vol = 0.;
double dptilde_VoldGamma= 0.;
double PhiEq_cor = 0.;
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++){
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++){
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;
}
else{
if (f>_tol){
q1->getRefToDissipationActive() = true;
// plasticity
int ite = 0;
int maxite = 100; // maximal number of iters
//Msg::Error("plasticity occurs f = %e",f);
//double f0 = fabs(f);
while (fabs(f) >_tol or ite <1){
double eta(0.),Deta(0.);
if (_viscosity != NULL)
_viscosity->get(q1->_epspbarre,eta,Deta);
double etaOverDt = eta/this->getTimeStep();
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+= 1./(2.*A)*(18./(u*u)*(devPhipr_Dev_cor(i,j)*dK_devdGamma(j,i)));
}
}
dDgammaDGamma = kk*(A+Gamma*dAdGamma);
this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
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)*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));
}
}
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.;
}
else
Gamma += dGamma;
//Msg::Error("Gamma = %e",Gamma);
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.;
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;
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++;
}
//################################### 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;
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);
updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
hardening(q0,q1);
getYieldCoefficients(q1,a);
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);
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<< " 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){
Msg::Error("No convergence for plastic correction in mlawPowerYieldHyper nonAssociatedFlow Maxwell iter = %d, f = %e!!",ite,f);
P(0,0) = P(1,1) = P(2,2) = sqrt(-1.);
return;
}
};
q1->_DgammaDt = Dgamma/this->getTimeStep();
// update effective stress tensor
devPhi =devPhipr_Dev_cor;
devPhi *= (1./u);
devN = devPhi;
devN *= 3.;
trN = 2.*_b*ptilde;
N = devN;
N(0,0) += trN/3.;
N(1,1) += trN/3.;
N(2,2) += trN/3.;
// estimate exp(GammaN)
static STensor3 expGN;
static STensor3 GammaN;
GammaN = N;
GammaN *= Gamma;
STensorOperation::expSTensor3(GammaN,_order,expGN,&dexpAdA);
// update plastic deformation tensor
STensorOperation::multSTensor3(expGN,Fp0,Fp1);
// update IP
updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
// update elastic deformation tensor, corotational stress
STensorOperation::inverseSTensor3(Fp1,Fpinv);
STensorOperation::multSTensor3(F,Fpinv,Fe);
STensorOperation::multSTensor3FirstTranspose(Fe,Fe,Ce);
STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCe,&DDlnDDCe);
Ee *= 0.5;
// update A, B
updateViscoElasticFlow(q0,q1,Ke,Ge);
// backstress
static STensor3 DB; // increment
DB = (N); // increment
DB *= (kk*Hb*Gamma);
q1->_backsig += DB; // update
// corotationaal Kirchhoff stress
q1->_kirchhoff = devPhi;
q1->_kirchhoff += q1->_backsig;
q1->_kirchhoff(0,0) += (ptilde);
q1->_kirchhoff(1,1) += (ptilde);
q1->_kirchhoff(2,2) += (ptilde);
}
else{
q1->getRefToDissipationActive() = false;
}
}
const STensor3& KS = q1->_kirchhoff;
// 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.;
for(int k=0; k<3; k++)
for(int l=0; l<3; l++)
P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
}
// defo energy
q1->getRefToElasticEnergy()=deformationEnergy(*q1);
q1->getRefToViscousEnergyPart()=viscousEnergy(*q0,*q1)+q0->getConstRefToViscousEnergyPart();
q1->getRefToPlasticEnergy() = q0->plasticEnergy();
if (Gamma > 0){
double dotKSN = dot(KS,N);
q1->getRefToPlasticEnergy() += Gamma*dotKSN;
}
if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
q1->getRefToIrreversibleEnergy() = q1->defoEnergy();
}
else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
(this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
q1->getRefToIrreversibleEnergy() = q1->plasticEnergy();
}
else{
q1->getRefToIrreversibleEnergy() = 0.;
}
if (stiff){
static STensor3 DpprDCepr;
static STensor43 DdevKprDCepr;
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;//+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;
DdevKDCepr = DdevKprDCepr;
static STensor43 dFpDCepr;
static STensor3 DgamaDCepr;
static STensor3 DtrNDCepr;
static STensor43 DdevNDCepr;
static STensor3 DGDCepr;
if (Gamma >0){
// plastic
static STensor3 dAdCepr, dfDCepr, dptilde_VoldCepr, dgamadCepr, DGDCepr, dptildeOldpr_VoldCepr;
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));
}
}
//########################################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*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)*(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;
}
}
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
DgamaDCepr(i,j) = kk*Gamma*dAdCepr(i,j)+ dDgammaDGamma*DGDCepr(i,j); // remove the extra kk
}
}
for (int i=0; i<3; i++)
for (int j=0; j<3; 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);
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) = (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;
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) = N(i,j)*DGDCepr(k,l)+ Gamma*DdevNDCepr(i,j,k,l)+ Gamma/3.*_I(i,j)*DtrNDCepr(k,l);
static STensor43 EprFp0;
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++){
EprFp0(i,j,k,l) = 0.;
for (int s=0; s<3; s++){
EprFp0(i,j,k,l) += dexpAdA(i,s,k,l)*Fp0(s,j);
}
}
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) = 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);
}
}
}
}
}
else{
// elastic
STensorOperation::zero(DgamaDCepr);
STensorOperation::zero(dFpDCepr);
STensorOperation::zero(DtrNDCepr);
STensorOperation::zero(DdevNDCepr);
STensorOperation::zero(DGDCepr);
}
static STensor43 dKcorDcepr;
dKcorDcepr = DdevKDCepr;
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++){
dKcorDcepr(i,j,k,l) += _I(i,j)*DpDCepr(k,l);
}
}
}
}
static STensor43 CeprToF;
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++){
CeprToF(i,j,k,l) = 2.*Fepr(k,i)*invFp0(j,l);
}
}
}
}
STensor3& DgammaDF = q1->_DgammaDF;
static STensor43 DKcorDF;
STensorOperation::multSTensor43(dKcorDcepr,CeprToF,DKcorDF);
if (Gamma > 0){
STensorOperation::multSTensor3STensor43(DgamaDCepr,CeprToF,DgammaDF);
STensorOperation::multSTensor43(dFpDCepr,CeprToF,dFpdF);
}
else{
STensorOperation::zero(DgammaDF);
STensorOperation::zero(dFpdF);
}
static STensor43 DinvFpDF; //
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DinvFpDF(i,s,k,l) = 0.;
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
DinvFpDF(i,s,k,l) -= Fpinv(i,m)*dFpdF(m,j,k,l)*Fpinv(j,s);
}
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dFedF(m,j,k,l) = _I(m,k)*Fpinv(l,j);
for (int s=0; s<3; s++)
dFedF(m,j,k,l) += F(m,s)*DinvFpDF(s,j,k,l);
}
static STensor63 DlnDF;
if (_order != 1){
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
for (int p=0; p<3; p++){
for (int q=0; q<3; q++){
DlnDF(i,j,k,l,p,q) = 0.;
for (int r=0; r<3; r++){
for (int s=0; s<3; s++){
for (int a=0; a<3; a++){
DlnDF(i,j,k,l,p,q) += DDlnDDCe(i,j,k,l,r,s)*2.*Fe(a,r)*dFedF(a,s,p,q);
}
}
}
}
}
}
}
}
}
}
else{
STensorOperation::zero(DlnDF);
}
static STensor43 dSdF;
STensorOperation::zero(dSdF);
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++){
dSdF(i,j,k,l) += DKcorDF(m,n,k,l)*DlnDCe(m,n,i,j);
dSdF(i,j,k,l) += KS(m,n)*DlnDF(m,n,i,j,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++){
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) += 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);
}
}
}
dDeformationEnergydF(*q1, P, q1->getRefToDElasticEnergyPartdF(), Fp1, dFedF);
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){
DirrEnergDF=q1->getConstRefToDElasticEnergyPartdF();
}
else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
(this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
DirrEnergDF=q1->getConstRefToDPlasticEnergyPartdF();
}
else{
STensorOperation::zero(DirrEnergDF);
}
};
};
void mlawPowerYieldHyper::dDeformationEnergydF(const IPHyperViscoElastoPlastic& q, const STensor3 &P, STensor3 &dElasticEnergydF, const STensor3 &Fp, const STensor43 &dFedF) const
{
STensorOperation::zero(dElasticEnergydF);
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++){
dElasticEnergydF(i,j) += P(k,m)*Fp(l,m)*dFedF(k,l,i,j);
}
}
}
}
}
}
void mlawPowerYieldHyper::dViscousEnergydF(const IPHyperViscoElastoPlastic& q0, const IPHyperViscoElastoPlastic& q, const STensor43 &dlnCdC, const STensor3 &Fe, STensor3 &dViscousEnergydF, const STensor43 &dFedF) const
{
STensorOperation::zero(dViscousEnergydF);
static STensor3 dViscousEnergyPartdFe;
STensorOperation::zero(dViscousEnergyPartdFe);
mlawHyperViscoElastic::dViscousEnergydF(dynamic_cast<const IPHyperViscoElastic &>(q0),dynamic_cast<const IPHyperViscoElastic &>(q),dlnCdC, Fe, dViscousEnergyPartdFe);
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++)
{
dViscousEnergydF(i,J)+=dViscousEnergyPartdFe(k,L)*dFedF(k,L,i,J);
}
}
}
}
}
void mlawPowerYieldHyper::dPlasticEnergydF(const STensor3 &N, double Gamma, double Dgamma, const STensor43 &dKcorDcepr, const STensor3 &DGDCepr, const STensor43 &DdevNDCepr, const STensor3 &DtrNDCepr, const STensor3 &KS, const STensor43 &CeprToF, STensor3 &dPlasticEnergydF) const
{
STensorOperation::zero(dPlasticEnergydF);
if (Dgamma > 0){
static STensor3 DirrEnergDCepr;
double dotKSN = dot(KS,N);
DirrEnergDCepr = DGDCepr;
DirrEnergDCepr *= dotKSN;
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++){
DirrEnergDCepr(i,j) += (Gamma*dKcorDcepr(k,l,i,j)*N(k,l)+ Gamma*KS(k,l)*(DdevNDCepr(k,l,i,j)+_I(k,l)*DtrNDCepr(i,j)/3.));
}
}
}
}
STensorOperation::multSTensor3STensor43(DirrEnergDCepr,CeprToF,dPlasticEnergydF);
}
}
void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1,
STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF) const{
/* compute elastic predictor */
STensor3& Fp1 = q1->_Fp;
const STensor3& Fp0 = q0->_Fp;
Fp1 = Fp0; // plastic deformation tensor
q1->_epspbarre = q0->_epspbarre; // plastic equivalent strain
q1->_epspCompression = q0->_epspCompression;
q1->_epspTraction = q0->_epspTraction;
q1->_epspShear = q0->_epspShear;
q1->_backsig = q0->_backsig; // backstress tensor
q1->_DgammaDt = 0.; // plastic rate --> failure
static STensor3 Fpinv, Ce, Fepr;
STensorOperation::inverseSTensor3(Fp1,Fpinv);
STensorOperation::multSTensor3(F,Fpinv,Fepr);
STensorOperation::multSTensor3FirstTranspose(Fepr,Fepr,Ce);
static STensor3 invFp0; // plastic predictor
invFp0= Fpinv;
STensor3& Fe = q1->_Fe;
Fe = Fepr;
static STensor43 DlnDCepr, DlnDCe;
static STensor63 DDlnDDCe;
static STensor43 dexpAdA; // estimation of dexpA/dA
STensor3& Ee = q1->_Ee;
STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCepr,&DDlnDDCe);
Ee *= 0.5;
DlnDCe = DlnDCepr;
// update A, B
double Ge, Ke;
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
STensorOperation::decomposeDevTr(q1->_kirchhoff,devKpr,ppr);
ppr /= 3.;
double keqpr = sqrt(1.5*devKpr.dotprod());
static STensor3 devK;
devK= devKpr; // dev corotational kirchoff stress
double p = ppr; // pressure
// hardening
this->hardening(q0,q1);
fullVector<double> a(3), Da(3); // yield coefficients and derivatives in respect to plastic deformation
this->getYieldCoefficients(q1,a);
static STensor3 devN; // dev part of yield normal
double trN = 0.; // trace part of yield normal
static STensor3 N; // yield normal
STensorOperation::zero(devN);
STensorOperation::zero(N);
double sigVM = keqpr;
double Dgamma = 0.;
double Gamma = 0.;
double g0 = 0.;
double z = a(1)/a(2)*pow(sigVM,1.-_n);
double A = sqrt(1.5+1.*z*z/(3.*_n*_n));
static fullMatrix<double> invJ(2,2); // inversed jacobian
invJ.setAll(0.);
double dAdz =0.;
double dzdDgamma = 0.;
double dzdsigVM = 0.;
double dg0dsigVM = -1./(3.*Ge);
if (q1->dissipationIsBlocked()){
q1->getRefToDissipationActive() = false;
}
else{
double f = a(2)*pow(sigVM,_n) - (a(1)*ppr+a(0));
if (f>_tol){
q1->getRefToDissipationActive() = true;
// plasticity
int ite = 0;
int maxite = 1000; // maximal number of iters
double val = 1.5*a(2)*_n*pow(3.*fabs(p),_n-2.);
q1->_nup = (val*p+ a(1)/3.)/(val*2.*p - a(1)/3.);
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
A*=kk;
fullVector<double> R(2); // residual of eqs to estimate Dgamma, Dpression and Gamma
R(0) = Dgamma - g0*A;
R(1) = a(2)*pow(sigVM,_n) - a(1)*(ppr+Ke*z*g0/_n) - a(0);
double res = R.norm();
while(res>_tol or ite<1){
this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
static fullMatrix<double> J(2,2);
J.setAll(0.);
dAdz = kk*z/(3.*_n*_n*A);
dzdDgamma = (Da(1)*a(2) - a(1)*Da(2))/(a(2)*a(2))*pow(sigVM,1.-_n);
dzdsigVM = a(1)/a(2)*(1.-_n)*pow(sigVM,-_n);
J(0,0) = 1. - g0*dAdz*dzdDgamma;
J(0,1) = -dg0dsigVM*A-g0*dAdz*dzdsigVM;
J(1,0) = Da(2)*pow(sigVM,_n) - Da(1)*(ppr+Ke*z*g0/_n) - Da(0) - a(1)*Ke*dzdDgamma*g0/_n;
J(1,1) = a(2)*_n*pow(sigVM,_n-1) - a(1)*Ke/_n*(dzdsigVM*g0+z*dg0dsigVM);
double detJ = J(0,0)*J(1,1) - J(1,0)*J(0,1);
if (detJ == 0.) Msg::Error("the corrected system can not be solved: mlawPowerYieldHyper::predictorCorrector_associatedFlow");
invJ(0,0) = J(1,1)/detJ;
invJ(0,1) = -J(0,1)/detJ;
invJ(1,0) = -J(1,0)/detJ;
invJ(1,1) = J(0,0)/detJ;
double solDgamma = -invJ(0,0)*R(0)-invJ(0,1)*R(1);
double solsigVM = -invJ(1,0)*R(0)-invJ(1,1)*R(1);
// bool isSolved = J.luSolve(residual,sol);
// if (!isSolved) Msg::Error("the corrected system can not be solved: mlawPowerYieldHyper::predictorCorrector");
//sol.print("sol");
if (solDgamma >0.1)
Dgamma+= 0.1;
else if (Dgamma+ solDgamma< 0.)
Dgamma /= 2.;
else
Dgamma += solDgamma;
if (sigVM + solsigVM< 0.)
sigVM /= 2.;
else
sigVM += solsigVM;
updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
this->hardening(q0,q1);
this->getYieldCoefficients(q1,a);
g0 = (keqpr- sigVM)/(3.*Ge);
z = a(1)/a(2)*pow(sigVM,1.-_n);
A = kk*sqrt(1.5+1.*z*z/(3.*_n*_n));
R(0) = Dgamma - g0*A;
R(1) = a(2)*pow(sigVM,_n) - a(1)*(ppr+Ke*z*g0/_n) - a(0);
res = R.norm();
ite++;
//Msg::Info("iter = %d, res = %e, Dgamma = %e, sigVM/keqpr = %e",ite,res, Dgamma, sigVM/keqpr);
if(ite > maxite){
Msg::Error("No convergence for plastic correction in mlawPowerYieldHyper iter = %d, res = %e!!",ite,res);
P(0,0) = P(1,1) = P(2,2) = sqrt(-1.);
break;
}
}
Gamma = z*g0/(_n*a(1));
// update
p += Ke*Gamma*a(1);
double ff =a(2)*_n*pow(sigVM,_n-2.);
devK*= 1./(1.+3.*Ge*Gamma*ff);
// estimate yield normal
devN = devK;
devN *= 1.5*ff;
trN = -a(1);
N = devN;
N(0,0) += trN/3.;
N(1,1) += trN/3.;
N(2,2) += trN/3.;
// estimate exp(GammaN)
static STensor3 expGN;
static STensor3 GammaN;
GammaN = N;
GammaN *= Gamma;
STensorOperation::expSTensor3(GammaN,_order,expGN,&dexpAdA);
// update plastic deformation tensor
STensorOperation::multSTensor3(expGN,Fp0,Fp1);
// update IP
updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
// update elastic deformation tensor, corotational stress
STensorOperation::inverseSTensor3(Fp1,Fpinv);
STensorOperation::multSTensor3(F,Fpinv,Fe);
STensorOperation::multSTensor3FirstTranspose(Fe,Fe,Ce);
STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCe,&DDlnDDCe);
Ee *= 0.5;
updateViscoElasticFlow(q0,q1,Ke,Ge);
}
else{
q1->getRefToDissipationActive() = false;
}
}
// corotational Kirchhoff stress tenor
STensor3& KS = q1->_kirchhoff;
KS = devK;
KS(0,0) += p;
KS(1,1) += p;
KS(2,2) += p;
// first Piola Kirchhoff stress
static STensor3 S;
S*=(0.);
for(int i=0; i<3; i++)
for(int j=0; j<3; j++)
for(int k=0; k<3; k++)
for(int l=0; l<3; l++)
S(i,j)+= KS(k,l)*DlnDCe(k,l,i,j);
STensorOperation::zero(P);
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++)
P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
// defo energy
q1->_elasticEnergy=deformationEnergy(*q1);
q1->getRefToPlasticEnergy() = q0->plasticEnergy();
if (Gamma > 0){
double dotKSN = dot(KS,N);
q1->getRefToPlasticEnergy() += Gamma*dotKSN;
}
if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
q1->getRefToIrreversibleEnergy() = q1->defoEnergy();
}
else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
(this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
q1->getRefToIrreversibleEnergy() = q1->plasticEnergy();
}
else{
q1->getRefToIrreversibleEnergy() = 0.;
}
if (stiff){
static STensor3 DpprDCepr;
static STensor43 DdevKprDCepr;
STensorOperation::multSTensor3STensor43(_I,DlnDCepr,DpprDCepr);
DpprDCepr*= (0.5*Ke);
STensorOperation::multSTensor43(_Idev,DlnDCepr,DdevKprDCepr);
DdevKprDCepr*= Ge;
static STensor3 DpDCepr;
static STensor43 DdevKDCepr;
DpDCepr = DpprDCepr;
DdevKDCepr = DdevKprDCepr;
static STensor43 dFpDCepr;
static STensor3 DgamaDCepr;
static STensor3 DGDCepr;
static STensor3 DtrNDCepr;
static STensor43 DdevNDCepr;
if (Dgamma > 0.){
static STensor3 dg0dCepr;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
dg0dCepr(i,j) = 0.;
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
dg0dCepr(i,j) += devKpr(k,l)*DdevKprDCepr(k,l,i,j)/(2.*Ge*keqpr);
}
static STensor3 DPhigammaDCepr;
DPhigammaDCepr = dg0dCepr;
DPhigammaDCepr *= -A;
static STensor3 DPhiSigVMDCepr;
DPhiSigVMDCepr = (dg0dCepr);
DPhiSigVMDCepr *= Ke*z/_n;
DPhiSigVMDCepr += DpprDCepr;
DPhiSigVMDCepr *= -a(1);
static STensor3 DsigVMDCepr;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
DgamaDCepr(i,j) = -invJ(0,0)*DPhigammaDCepr(i,j) - invJ(0,1)*DPhiSigVMDCepr(i,j);
DsigVMDCepr(i,j) = -invJ(1,0)*DPhigammaDCepr(i,j) - invJ(1,1)*DPhiSigVMDCepr(i,j);
}
for (int i=0; i<3; i++)
for (int j=0; j<3; j++){
DGDCepr(i,j) = (g0/(_n*a(1)))*(dzdDgamma*DgamaDCepr(i,j)+dzdsigVM*DsigVMDCepr(i,j)) +
z/(_n*a(1))*(dg0dsigVM*DsigVMDCepr(i,j)+ dg0dCepr(i,j)) - z*g0/(_n*a(1)*a(1))*Da(1)*DgamaDCepr(i,j);
}
DtrNDCepr = (DgamaDCepr);
DtrNDCepr *= -Da(1);
DdevNDCepr = (DdevKprDCepr);
double B = a(2)*pow(sigVM,_n-2.);
double u = 1.+3.*Ge*Gamma*_n*B;
double u2 = u*u;
static STensor3 DdevNDB;
DdevNDB = (devKpr);
DdevNDB *= (1.5*_n/u2);
static STensor3 DdevNDGamma;
DdevNDGamma = (devKpr);
DdevNDGamma *= (-1.5*_n*B*3.*Ge*_n*B/u2);
DdevNDCepr *= (1.5*_n*B/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) +=DdevNDB(i,j)*(Da(2)*pow(sigVM,_n-2.)*DgamaDCepr(k,l)+
a(2)*(_n-2.)*pow(sigVM,_n-3.)*DsigVMDCepr(k,l));
DdevNDCepr(i,j,k,l) += DdevNDGamma(i,j)*DGDCepr(k,l);
}
// compute dFpdF
static STensor43 temp1;
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) = N(i,j)*DGDCepr(k,l)+ Gamma*DdevNDCepr(i,j,k,l)+ Gamma/3.*_I(i,j)*DtrNDCepr(k,l);
static STensor43 EprFp0;
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++){
EprFp0(i,j,k,l) = 0.;
for (int s=0; s<3; s++){
EprFp0(i,j,k,l) += dexpAdA(i,s,k,l)*Fp0(s,j);
}
}
STensorOperation::multSTensor43(EprFp0,temp1,dFpDCepr);
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));
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));
}
}
}
}
}
else{
STensorOperation::zero(DgamaDCepr);
STensorOperation::zero(dFpDCepr);
STensorOperation::zero(DGDCepr);
STensorOperation::zero(DdevNDCepr);
STensorOperation::zero(DtrNDCepr);
}
static STensor43 dKcorDcepr;
dKcorDcepr = DdevKDCepr;
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++){
dKcorDcepr(i,j,k,l) += _I(i,j)*DpDCepr(k,l);
}
}
}
}
static STensor43 CeprToF;
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++){
CeprToF(i,j,k,l) = 2.*Fepr(k,i)*invFp0(j,l);
}
}
}
}
STensor3& DgammaDF = q1->_DgammaDF;
static STensor43 DKcorDF;
STensorOperation::multSTensor3STensor43(DgamaDCepr,CeprToF,DgammaDF);
STensorOperation::multSTensor43(dKcorDcepr,CeprToF,DKcorDF);
STensorOperation::multSTensor43(dFpDCepr,CeprToF,dFpdF);
static STensor43 DinvFpDF; //
for (int i=0; i<3; i++)
for (int s=0; s<3; s++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
DinvFpDF(i,s,k,l) = 0.;
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
DinvFpDF(i,s,k,l) -= Fpinv(i,m)*dFpdF(m,j,k,l)*Fpinv(j,s);
}
for (int m=0; m<3; m++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++){
dFedF(m,j,k,l) = _I(m,k)*Fpinv(l,j);
for (int s=0; s<3; s++)
dFedF(m,j,k,l) += F(m,s)*DinvFpDF(s,j,k,l);
}
static STensor63 DlnDF;
if (_order != 1){
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
for (int p=0; p<3; p++){
for (int q=0; q<3; q++){
DlnDF(i,j,k,l,p,q) = 0.;
for (int r=0; r<3; r++){
for (int s=0; s<3; s++){
for (int a=0; a<3; a++){
DlnDF(i,j,k,l,p,q) += DDlnDDCe(i,j,k,l,r,s)*2.*Fe(a,r)*dFedF(a,s,p,q);
}
}
}
}
}
}
}
}
}
}
else{
STensorOperation::zero(DlnDF);
}
static STensor43 dSdF;
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++){
dSdF(i,j,k,l) = 0.;
for (int m=0; m<3; m++)
for (int n=0; n<3; n++){
dSdF(i,j,k,l) += DKcorDF(m,n,k,l)*DlnDCe(m,n,i,j);
if (_order != 1){
dSdF(i,j,k,l) += KS(m,n)*DlnDF(m,n,i,j,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++){
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) += 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);
}
}
}
STensor3& DirrEnergDF = q1->_DirreversibleEnergyDF;
if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
for(int i=0;i<3;i++){
for(int j=0;j<3;j++){
DirrEnergDF(i,j) = 0.;
for(int k=0;k<3;k++){
for(int l=0;l<3;l++){
for (int m=0; m<3; m++){
DirrEnergDF(i,j) += P(k,m)*Fp1(l,m)*dFedF(k,l,i,j);
}
}
}
}
}
}
else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
(this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
if (Dgamma > 0){
static STensor3 DirrEnergDCepr;
double dotKSN = dot(KS,N);
DirrEnergDCepr = DGDCepr;
DirrEnergDCepr *= dotKSN;
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++){
DirrEnergDCepr(i,j) += (Gamma*dKcorDcepr(k,l,i,j)*N(k,l)+ Gamma*KS(k,l)*(DdevNDCepr(k,l,i,j)+_I(k,l)*DtrNDCepr(i,j)/3.));
}
}
}
}
STensorOperation::multSTensor3STensor43(DirrEnergDCepr,CeprToF,DirrEnergDF);
}
else{
STensorOperation::zero(DirrEnergDF);
}
}
else{
STensorOperation::zero(DirrEnergDF);
}
}
};
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);
}
else{
this->predictorCorrector_associatedFlow(F,q0,q1,P,false,Tangent,dFedF,dFpdF);
}
if (stiff){
tangent_full_perturbation(Tangent,dFedF,dFpdF,P,F,q0,q1);
}
}
else{
if (_nonAssociatedFlow){
this->predictorCorrector_nonAssociatedFlow(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF);
}
else{
this->predictorCorrector_associatedFlow(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF);
}
}
// compute mechanical tengent
if (elasticTangent!= NULL)
{
STensor3 Ptmp(0.);
static IPHyperViscoElastoPlastic q1tmp(*q1);
q1tmp.operator=(*(dynamic_cast<const IPVariable*> (q1)));
mlawHyperViscoElastic::predictorCorrector_ViscoElastic(q0->getConstRefToFe(), q1->getConstRefToFe(), Ptmp, q0, &q1tmp,
true, *elasticTangent);
}
};
mlawPowerYieldHyper::mlawPowerYieldHyper(const int num,const double E,const double nu, const double rho,
const double tol,
const bool matrixbyPerturbation, const double pert):
mlawHyperViscoElastic(num,E,nu,rho,matrixbyPerturbation,pert),_tol(tol),
_n(1.5),_nonAssociatedFlow(true),_b(0.3),
_viscosity(NULL),_p(1.),_compression(NULL),_traction(NULL),_kinematic(NULL){
};
mlawPowerYieldHyper::mlawPowerYieldHyper(const mlawPowerYieldHyper& src):mlawHyperViscoElastic(src),
_tol(src._tol),_n(src._n),_nonAssociatedFlow(src._nonAssociatedFlow),_b(src._b), _p(src._p){
_viscosity = NULL;
if (src._viscosity) _viscosity = src._viscosity->clone();
_compression = NULL;
if (src._compression) _compression = src._compression->clone();
_traction = NULL;
if (src._traction) _traction = src._traction->clone();
_kinematic= NULL;
if (src._kinematic) _kinematic = src._kinematic->clone();
};
mlawPowerYieldHyper& mlawPowerYieldHyper::operator=(const materialLaw& source){
mlawHyperViscoElastic::operator=(source);
const mlawPowerYieldHyper* src =dynamic_cast<const mlawPowerYieldHyper*>(&source);
if(src != NULL){
_tol = src->_tol;
_n = src->_n;
_nonAssociatedFlow = src->_nonAssociatedFlow;
_b = src->_b;
_p = src->_p;
if(_viscosity != NULL) delete _viscosity;
if (src->_viscosity != NULL){ _viscosity = src->_viscosity->clone();}
if(_compression != NULL) delete _compression;
if (src->_compression!=NULL){ _compression = src->_compression->clone();}
if(_traction != NULL) delete _traction;
if (src->_traction!=NULL){ _traction = src->_traction->clone();}
if(_kinematic != NULL) delete _kinematic;
if (src->_kinematic!=NULL){ _kinematic = src->_kinematic->clone();}
}
return *this;
}
mlawPowerYieldHyper::~mlawPowerYieldHyper(){
if (_compression) delete _compression; _compression = NULL;
if (_traction) delete _traction; _traction = NULL;
if (_kinematic) delete _kinematic; _kinematic = NULL;
if (_viscosity) delete _viscosity; _viscosity = NULL;
};
void mlawPowerYieldHyper::setPowerFactor(const double n) {
_n = n;
};
void mlawPowerYieldHyper::nonAssociatedFlowRuleFactor(const double b){
_b = b;
};
void mlawPowerYieldHyper::setPlasticPoissonRatio(const double nup){
_b = 4.5*(1.-2.*nup)/(1.+nup);
};
void mlawPowerYieldHyper::setNonAssociatedFlow(const bool flag) {
_nonAssociatedFlow = flag;
if (_nonAssociatedFlow){
Msg::Info("non associated flow is used");
}
else{
Msg::Info("associated flow is used, kinematic hardening is not considered");
}
};
void mlawPowerYieldHyper::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
IPVariable* ipvi = new IPHyperViscoElastoPlastic(_compression,_traction,_kinematic,_N);
IPVariable* ipv1 = new IPHyperViscoElastoPlastic(_compression,_traction,_kinematic,_N);
IPVariable* ipv2 = new IPHyperViscoElastoPlastic(_compression,_traction,_kinematic,_N);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
};
mlawPowerYieldHyperWithFailure::mlawPowerYieldHyperWithFailure(const int num,const double E,const double nu, const double rho,
const double tol, const bool matrixbyPerturbation, const double pert):
mlawPowerYieldHyper(num,E,nu,rho,tol,matrixbyPerturbation,pert){
_failureCriterion = new mlawHyperelasticFailureTrivialCriterion(num);
};
mlawPowerYieldHyperWithFailure::mlawPowerYieldHyperWithFailure(const mlawPowerYieldHyperWithFailure& src):
mlawPowerYieldHyper(src){
_failureCriterion = NULL;
if (src._failureCriterion){
_failureCriterion = src._failureCriterion->clone();
}
};
mlawPowerYieldHyperWithFailure::~mlawPowerYieldHyperWithFailure(){
if (_failureCriterion != NULL) {
delete _failureCriterion;
_failureCriterion = NULL;
}
};
void mlawPowerYieldHyperWithFailure::checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{
IPHyperViscoElastoPlastic *q1 = dynamic_cast<IPHyperViscoElastoPlastic*>(ipv);
const IPHyperViscoElastoPlastic *q0 = dynamic_cast<const IPHyperViscoElastoPlastic*>(ipvprev);
if (q1 != NULL){
double& r = q1->getRefToFailureOnset();
r = q0->getFailureOnset();
double failCr = _failureCriterion->getFailureCriterion(q0,q1);
if (failCr > r){
r = failCr;
}
// check failure
if (q0->inPostFailureStage()){
q1->getRefToInPostFailureStage() = true;
}
else{
if (r > 0.){
q1->getRefToInPostFailureStage() = true;
}
else{
q1->getRefToInPostFailureStage() = false;
}
}
}
};
void mlawPowerYieldHyperWithFailure::predictorCorrector(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1,
STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF,
STensor43* elasticTangent) const{
mlawPowerYieldHyper::predictorCorrector(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF, elasticTangent);
// compute failure r and DrDF
double& gF = q1->getRefToFailurePlasticity();
gF = q0->getFailurePlasticity();
if (q1->dissipationIsBlocked()){
if (stiff){
STensorOperation::zero(q1->getRefToDFailurePlasticityDF());
}
}
else{
if (q1->inPostFailureStage()){
gF += q1->_epspbarre - q0->_epspbarre;
if (stiff){
q1->getRefToDFailurePlasticityDF() = q1->_DgammaDF;
}
}
else{
if (stiff){
STensorOperation::zero(q1->getRefToDFailurePlasticityDF());
}
}
}
};
void mlawPowerYieldHyperWithFailure::setFailureCriterion(const mlawHyperelasticFailureCriterionBase& fCr){
if (_failureCriterion) delete _failureCriterion;
_failureCriterion = fCr.clone();
};