Select Git revision
mlawNonLocalDamageHyperelastic.cpp
mlawNonLocalDamageHyperelastic.cpp 12.96 KiB
//
// C++ Interface: material law
//
// Description: mlawNonLocalDamageQuadYieldHyper
//
// Author: V.D. Nguyen, (C) 2014
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "mlawNonLocalDamageHyperelastic.h"
mlawNonLocalDamagePowerYieldHyper::mlawNonLocalDamagePowerYieldHyper(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),cLLaw(NULL),damLaw(NULL){};
mlawNonLocalDamagePowerYieldHyper::mlawNonLocalDamagePowerYieldHyper(const mlawNonLocalDamagePowerYieldHyper &source):mlawPowerYieldHyper(source){
cLLaw = NULL;
damLaw = NULL;
if(source.cLLaw != NULL)
{
cLLaw=source.cLLaw->clone();
}
if(source.damLaw != NULL)
{
damLaw=source.damLaw->clone();
};
}
void mlawNonLocalDamagePowerYieldHyper::setCLengthLaw(const CLengthLaw &_cLLaw){
if (cLLaw) delete cLLaw;
cLLaw = _cLLaw.clone();
};
void mlawNonLocalDamagePowerYieldHyper::setDamageLaw(const DamageLaw &_damLaw){
if (damLaw) delete damLaw;
damLaw = _damLaw.clone();
};
mlawNonLocalDamagePowerYieldHyper::~mlawNonLocalDamagePowerYieldHyper(){
if (cLLaw!= NULL) delete cLLaw;
if (damLaw!= NULL) delete damLaw;
};
void mlawNonLocalDamagePowerYieldHyper::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
IPVariable* ipvi = new IPHyperViscoElastoPlasticNonLocalDamage(_compression,_traction,_shear,_kinematic,_N,cLLaw,damLaw);
IPVariable* ipv1 = new IPHyperViscoElastoPlasticNonLocalDamage(_compression,_traction,_shear,_kinematic,_N,cLLaw,damLaw);
IPVariable* ipv2 = new IPHyperViscoElastoPlasticNonLocalDamage(_compression,_traction,_shear,_kinematic,_N,cLLaw,damLaw);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void mlawNonLocalDamagePowerYieldHyper::createIPState(IPHyperViscoElastoPlasticNonLocalDamage *ivi, IPHyperViscoElastoPlasticNonLocalDamage *iv1, IPHyperViscoElastoPlasticNonLocalDamage *iv2) const
{
}
void mlawNonLocalDamagePowerYieldHyper::createIPVariable(IPHyperViscoElastoPlasticNonLocalDamage *&ipv,const MElement *ele,const int nbFF,const IntPt *GP, const int gpt) const
{
}
double mlawNonLocalDamagePowerYieldHyper::soundSpeed() const
{
return mlawPowerYieldHyper::soundSpeed();
}
void mlawNonLocalDamagePowerYieldHyper::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P,const IPHyperViscoElastoPlasticNonLocalDamage *q0, IPHyperViscoElastoPlasticNonLocalDamage *q1,STensor43 &Tangent,
const bool stiff) const
{
mlawPowerYieldHyper::constitutive(F0,Fn,P,q0,q1,Tangent,stiff);
}
void mlawNonLocalDamagePowerYieldHyper::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)
// contains the initial values on input
const IPHyperViscoElastoPlasticNonLocalDamage *ipvprev, // array of initial internal variable
IPHyperViscoElastoPlasticNonLocalDamage *ipvcur, // updated array of internal variable (in ipvcur on output),
STensor43 &Tangent, // constitutive tangents (output)
STensor3 &dLocalPlasticStrainDStrain,
STensor3 &dStressDNonLocalPlasticStrain,
double &dLocalPlasticStrainDNonLocalPlasticStrain,
const bool stiff // if true compute the tangents
) const
{
double p0 = ipvprev->getCurrentPlasticStrain();
cLLaw->computeCL(p0, ipvcur->getRefToIPCLength());
static STensor43 dFedF, dFpdF;
static STensor3 Peff;
mlawPowerYieldHyper::predictorCorrector(Fn,ipvprev,ipvcur,Peff,stiff,Tangent,dFedF,dFpdF);
const STensor3& Fe = ipvcur->_Fe;
const STensor3& dgammadF = ipvcur->_DgammaDF;
const STensor3 &Fp = ipvcur->getConstRefToFp();
double ene = ipvcur->defoEnergy();
//Msg::Info("enery = %e",ene);
damLaw->computeDamage(ipvcur->getEffectivePlasticStrain(),
ipvprev->getEffectivePlasticStrain(),
ene, Fe, Fp, Peff, Cel,
ipvprev->getConstRefToIPDamage(),ipvcur->getRefToIPDamage());
double D = ipvcur->getDamage();
P = Peff;
P*=(1.-D);
if(stiff)
{
// we need to correct partial P/partial F: (1-D) partial P/partial F - Peff otimes partial D partial F
Tangent*=(1.-D);
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++)
{
Tangent(i,j,k,l)-=Peff(i,j)*ipvcur->getConstRefToDDamageDFe()(m,n)*dFedF(m,n,k,l);
}
}
}
}
}
}
// partial p/partial F
dLocalPlasticStrainDStrain = dgammadF;
// -hat{P} partial D/partial tilde p
dStressDNonLocalPlasticStrain = Peff*(-1*ipvcur->getDDamageDp());
// partial p partial tilde p (0 if no MFH)
dLocalPlasticStrainDNonLocalPlasticStrain = 0.;
}
}
mlawNonLocalDamagePowerYieldHyperWithFailure::mlawNonLocalDamagePowerYieldHyperWithFailure(const int num,const double E,const double nu, const double rho,
const double tol, const bool matrixbyPerturbation , const double pert):
mlawPowerYieldHyperWithFailure(num,E,nu,rho,tol,matrixbyPerturbation,pert), _saturated(false){};
void mlawNonLocalDamagePowerYieldHyperWithFailure::clearAllCLengthLaw(){
for (int i=0; i< cLLaw.size(); i++){
if (cLLaw[i]!= NULL) delete cLLaw[i];
}
cLLaw.clear();
};
void mlawNonLocalDamagePowerYieldHyperWithFailure::clearAllDamageLaw(){
for (int i=0; i< damLaw.size(); i++){
if (damLaw[i]!= NULL) delete damLaw[i];
}
damLaw.clear();
};
void mlawNonLocalDamagePowerYieldHyperWithFailure::setCLengthLaw(const CLengthLaw &_cLLaw){
cLLaw.push_back(_cLLaw.clone());
};
void mlawNonLocalDamagePowerYieldHyperWithFailure::setDamageLaw(const DamageLaw &_damLaw){
damLaw.push_back(_damLaw.clone());
};
mlawNonLocalDamagePowerYieldHyperWithFailure::mlawNonLocalDamagePowerYieldHyperWithFailure(const mlawNonLocalDamagePowerYieldHyperWithFailure &source):
mlawPowerYieldHyperWithFailure(source),_saturated(source._saturated){
cLLaw.clear();
damLaw.clear();
for (int i=0; i< source.cLLaw.size(); i++){
if(source.cLLaw[i] != NULL)
{
cLLaw.push_back(source.cLLaw[i]->clone());
}
if(source.damLaw[i] != NULL)
{
damLaw.push_back(source.damLaw[i]->clone());
};
}
}
mlawNonLocalDamagePowerYieldHyperWithFailure::~mlawNonLocalDamagePowerYieldHyperWithFailure(){
for (int i=0; i< cLLaw.size(); i++){
if (cLLaw[i]!= NULL) delete cLLaw[i];
if (damLaw[i]!= NULL) delete damLaw[i];
}
cLLaw.clear();
damLaw.clear();
};
void mlawNonLocalDamagePowerYieldHyperWithFailure::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
IPVariable* ipvi = new IPHyperViscoElastoPlasticMultipleNonLocalDamage(_compression,_traction,_shear,_kinematic,_N,cLLaw,damLaw);
IPVariable* ipv1 = new IPHyperViscoElastoPlasticMultipleNonLocalDamage(_compression,_traction,_shear,_kinematic,_N,cLLaw,damLaw); IPVariable* ipv2 = new IPHyperViscoElastoPlasticMultipleNonLocalDamage(_compression,_traction,_shear,_kinematic,_N,cLLaw,damLaw);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
void mlawNonLocalDamagePowerYieldHyperWithFailure::createIPState(IPHyperViscoElastoPlasticMultipleNonLocalDamage *ivi, IPHyperViscoElastoPlasticMultipleNonLocalDamage *iv1, IPHyperViscoElastoPlasticMultipleNonLocalDamage *iv2) const
{
}
void mlawNonLocalDamagePowerYieldHyperWithFailure::createIPVariable(IPHyperViscoElastoPlasticMultipleNonLocalDamage *&ipv,const MElement *ele,const int nbFF,const IntPt *GP, const int gpt) const
{
}
double mlawNonLocalDamagePowerYieldHyperWithFailure::soundSpeed() const
{
return mlawPowerYieldHyperWithFailure::soundSpeed();
}
void mlawNonLocalDamagePowerYieldHyperWithFailure::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P,const IPHyperViscoElastoPlasticMultipleNonLocalDamage *q0,
IPHyperViscoElastoPlasticMultipleNonLocalDamage *q1,STensor43 &Tangent,
const bool stiff) const
{
mlawPowerYieldHyperWithFailure::constitutive(F0,Fn,P,q0,q1,Tangent,stiff);
}
void mlawNonLocalDamagePowerYieldHyperWithFailure::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 IPHyperViscoElastoPlasticMultipleNonLocalDamage *ipvprev, // array of initial internal variable
IPHyperViscoElastoPlasticMultipleNonLocalDamage *ipvcur, // updated array of internal variable (in ipvcur on output),
STensor43 &Tangent, // constitutive tangents (output)
std::vector<STensor3> &dLocalVariableDStrain,
std::vector<STensor3> &dStressDNonLocalVariable,
fullMatrix<double> &dLocalVariableDNonLocalVariable,
const bool stiff // if true compute the tangents
) const
{
double p0 = ipvprev->getCurrentPlasticStrain();
// compute non-local length scales for first and second damage variable
cLLaw[0]->computeCL(p0, ipvcur->getRefToIPCLength(0));
cLLaw[1]->computeCL(p0, ipvcur->getRefToIPCLength(1));
if (isSaturatedHardening()){
// saturate when defined critical for damage law
if (ipvprev->getDamage(1) >= damLaw[1]->getCriticalDamage()){
ipvcur->saturateHardening(ipvprev);
}
}
static STensor43 dFedF, dFpdF;
static STensor3 Peff;
mlawPowerYieldHyperWithFailure::predictorCorrector(Fn,ipvprev,ipvcur,Peff,stiff,Tangent,dFedF,dFpdF);
// get result from effective law
const STensor3& Fe = ipvcur->_Fe;
const STensor3& dgammadF = ipvcur->_DgammaDF;
const STensor3 &Fp = ipvcur->getConstRefToFp();
const STensor3& dgFdF = ipvcur->getConstRefToDFailurePlasticityDF();
double ene = ipvcur->defoEnergy();
//Msg::Info("enery = %e",ene);
// saturation damage always develops
damLaw[0]->computeDamage(ipvcur->getEffectivePlasticStrain(),
ipvprev->getEffectivePlasticStrain(),
ene, Fe, Fp, Peff, Cel,
ipvprev->getConstRefToIPDamage(0),ipvcur->getRefToIPDamage(0));
if (ipvprev->getDamage(1) < ipvcur->getCriticalDamage()){
// set for current critical damage
ipvcur->setCriticalDamage(ipvcur->getCriticalDamage());
// damage evolution
damLaw[1]->computeDamage(ipvcur->getNonLocalFailurePlasticity(),
ipvprev->getNonLocalFailurePlasticity(),
ene, Fe, Fp, Peff, Cel,
ipvprev->getConstRefToIPDamage(1),ipvcur->getRefToIPDamage(1));
}
else{
// damage stop increasing
STensor3 dDDFe(0.);
IPDamage& curDama1 = ipvcur->getRefToIPDamage(1);
curDama1.setValues(ipvprev->getDamage(1),ipvcur->getMaximalP(1),0.,0.,dDDFe);
}
// computue total softening
double D1 = ipvcur->getDamage(0);
double D2 = ipvcur->getDamage(1);
double D = 1. - (1.- D1)*(1.-D2);
P = Peff;
P*= (1.- D);
if(stiff)
{
// we need to correct partial P/partial F: (1-D) partial P/partial F - Peff otimes partial D partial F
Tangent*=(1.-D);
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++)
{
Tangent(i,j,k,l)-=Peff(i,j)*((1.-D2)*ipvcur->getConstRefToDDamageDFe(0)(m,n) +(1.-D1)*ipvcur->getConstRefToDDamageDFe(1)(m,n)) *dFedF(m,n,k,l);
}
}
}
}
}
}
// -hat{P} partial D/partial tilde p
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
// partial p/partial F
dLocalVariableDStrain[0](i,j) = dgammadF(i,j);
dLocalVariableDStrain[1](i,j) = dgFdF(i,j);
dStressDNonLocalVariable[0](i,j) = -1.*Peff(i,j)*ipvcur->getDDamageDp(0)*(1.- D2);
dStressDNonLocalVariable[1](i,j) = -1.*Peff(i,j)*(1.-D1)*ipvcur->getDDamageDp(1);
}
}
// partial p partial tilde p (0 if no MFH)
dLocalVariableDNonLocalVariable.setAll(0.);
}
}