From 2669e3eecfe5866a8f30efd07be5f718288d6acf Mon Sep 17 00:00:00 2001
From: Van Dung Nguyen <vdg.nguyen@gmail.com>
Date: Fri, 25 Aug 2017 18:23:40 +0200
Subject: [PATCH] add dissipation
---
.../internalPoints/ipHyperelastic.cpp | 61 +++++--
.../internalPoints/ipHyperelastic.h | 12 +-
.../ipNonLocalDamageHyperelastic.cpp | 15 +-
.../ipNonLocalDamageHyperelastic.h | 8 +
.../materialLaw/mlawHyperelastic.cpp | 158 +++++++++++-----
.../materialLaw/mlawHyperelastic.h | 4 +-
.../mlawNonLocalDamageHyperelastic.cpp | 169 ++++++++----------
.../mlawNonLocalDamageHyperelastic.h | 13 --
.../materialLaw/mlawNonLocalDamageJ2Hyper.cpp | 14 +-
.../nlsolver/nonLinearMechSolver.cpp | 5 +-
dG3D/src/dG3DIPVariable.h | 16 +-
...nonLocalDamageHyperelasticDG3DIPVariable.h | 46 +++++
dG3D/src/pathFollowingTerms.cpp | 1 -
13 files changed, 341 insertions(+), 181 deletions(-)
diff --git a/NonLinearSolver/internalPoints/ipHyperelastic.cpp b/NonLinearSolver/internalPoints/ipHyperelastic.cpp
index ca34ee507..87549c92c 100644
--- a/NonLinearSolver/internalPoints/ipHyperelastic.cpp
+++ b/NonLinearSolver/internalPoints/ipHyperelastic.cpp
@@ -12,7 +12,8 @@
#include "ipHyperelastic.h"
#include "restartManager.h"
-IPHyperViscoElastic::IPHyperViscoElastic(const int N):IPVariableMechanics(),_N(N),_elasticEnergy(0.),_Ee(0.),_kirchhoff(0.){
+IPHyperViscoElastic::IPHyperViscoElastic(const int N):IPVariableMechanics(),_N(N),_elasticEnergy(0.),_Ee(0.),_kirchhoff(0.),
+ _irreversibleEnergy(0.),_DirreversibleEnergyDF(0.){
_A.clear();
_B.clear();
for (int i=0; i< N; i++){
@@ -22,7 +23,8 @@ IPHyperViscoElastic::IPHyperViscoElastic(const int N):IPVariableMechanics(),_N(N
}
};
IPHyperViscoElastic::IPHyperViscoElastic(const IPHyperViscoElastic& src): IPVariableMechanics(src), _Ee(src._Ee),
- _N(src._N),_A(src._A),_B(src._B),_elasticEnergy(src._elasticEnergy), _kirchhoff(src._kirchhoff){};
+ _N(src._N),_A(src._A),_B(src._B),_elasticEnergy(src._elasticEnergy), _kirchhoff(src._kirchhoff),
+ _irreversibleEnergy(src._irreversibleEnergy),_DirreversibleEnergyDF(src._DirreversibleEnergyDF){};
IPHyperViscoElastic& IPHyperViscoElastic::operator =(const IPVariable& src){
IPVariableMechanics::operator =(src);
@@ -34,6 +36,8 @@ IPHyperViscoElastic& IPHyperViscoElastic::operator =(const IPVariable& src){
_A = psrc->_A;
_B = psrc->_B;
_elasticEnergy = psrc->_elasticEnergy;
+ _irreversibleEnergy = psrc->_irreversibleEnergy;
+ _DirreversibleEnergyDF = psrc->_DirreversibleEnergyDF;
}
return *this;
};
@@ -46,6 +50,8 @@ void IPHyperViscoElastic::restart() {
restartManager::restart(_B);
restartManager::restart(_Ee);
restartManager::restart(_kirchhoff);
+ restartManager::restart(_irreversibleEnergy);
+ restartManager::restart(_DirreversibleEnergyDF);
}
@@ -113,21 +119,42 @@ IPHyperViscoElastoPlastic& IPHyperViscoElastoPlastic::operator =(const IPVariabl
_gF = ps->_gF;
_dgFdF = ps->_dgFdF;
- if (_ipCompression != NULL) delete _ipCompression;
- _ipCompression = NULL;
- if ( ps->_ipCompression != NULL) _ipCompression = dynamic_cast<IPJ2IsotropicHardening*>(ps->_ipCompression->clone());
-
- if (_ipTraction != NULL) delete _ipTraction;
- _ipTraction = NULL;
- if (ps->_ipTraction != NULL) _ipTraction = dynamic_cast<IPJ2IsotropicHardening*>(ps->_ipTraction->clone());
-
- if (_ipShear != NULL) delete _ipShear;
- _ipShear = NULL;
- if (ps->_ipShear != NULL) _ipShear = dynamic_cast<IPJ2IsotropicHardening*>(ps->_ipShear->clone());
-
- if (_ipKinematic!= NULL) delete _ipKinematic;
- _ipKinematic = NULL;
- if (ps->_ipKinematic != NULL) _ipKinematic = dynamic_cast<IPKinematicHardening*>(ps->_ipKinematic->clone());
+ if ( ps->_ipCompression != NULL) {
+ if (_ipCompression == NULL){
+ _ipCompression = dynamic_cast<IPJ2IsotropicHardening*>(ps->_ipCompression->clone());
+ }
+ else{
+ _ipCompression->operator=(*dynamic_cast<const IPVariable*>(ps->_ipCompression));
+ }
+ }
+
+ if ( ps->_ipTraction != NULL) {
+ if (_ipTraction == NULL){
+ _ipTraction = dynamic_cast<IPJ2IsotropicHardening*>(ps->_ipTraction->clone());
+ }
+ else{
+ _ipTraction->operator=(*dynamic_cast<const IPVariable*>(ps->_ipTraction));
+ }
+ }
+
+ if ( ps->_ipShear != NULL) {
+ if (_ipShear == NULL){
+ _ipShear = dynamic_cast<IPJ2IsotropicHardening*>(ps->_ipShear->clone());
+ }
+ else{
+ _ipShear->operator=(*dynamic_cast<const IPVariable*>(ps->_ipShear));
+ }
+ }
+
+ if ( ps->_ipKinematic != NULL) {
+ if (_ipKinematic == NULL){
+ _ipKinematic = dynamic_cast<IPKinematicHardening*>(ps->_ipKinematic->clone());
+ }
+ else{
+ _ipKinematic->operator=(*dynamic_cast<const IPVariable*>(ps->_ipKinematic));
+ }
+ }
+
}
return *this;
};
diff --git a/NonLinearSolver/internalPoints/ipHyperelastic.h b/NonLinearSolver/internalPoints/ipHyperelastic.h
index 849442746..0b665e40f 100644
--- a/NonLinearSolver/internalPoints/ipHyperelastic.h
+++ b/NonLinearSolver/internalPoints/ipHyperelastic.h
@@ -26,6 +26,9 @@ class IPHyperViscoElastic : public IPVariableMechanics{
STensor3 _Ee; // elastic strain
STensor3 _kirchhoff; // corotational Kirchhoff stress
+ double _irreversibleEnergy;
+ STensor3 _DirreversibleEnergyDF;
+
public:
IPHyperViscoElastic(const int N);
IPHyperViscoElastic(const IPHyperViscoElastic& src);
@@ -43,6 +46,13 @@ class IPHyperViscoElastic : public IPVariableMechanics{
virtual double defoEnergy() const{return _elasticEnergy;}
virtual double& getRefToElasticEnergy() {return _elasticEnergy;};
virtual double plasticEnergy() const{return 0.;}
+
+ virtual double irreversibleEnergy() const {return _irreversibleEnergy;};
+ virtual double & getRefToIrreversibleEnergy() {return _irreversibleEnergy;};
+
+ virtual STensor3& getRefToDIrreversibleEnergyDF() {return _DirreversibleEnergyDF;};
+ virtual const STensor3& getConstRefToDIrreversibleEnergyDF() const{return _DirreversibleEnergyDF;};
+
virtual void restart();
};
@@ -73,7 +83,7 @@ class IPHyperViscoElastoPlastic : public IPHyperViscoElastic{
double _r; // failure onset
double _gF; //
STensor3 _dgFdF;
-
+
public:
IPHyperViscoElastoPlastic(const J2IsotropicHardening* comp,
diff --git a/NonLinearSolver/internalPoints/ipNonLocalDamageHyperelastic.cpp b/NonLinearSolver/internalPoints/ipNonLocalDamageHyperelastic.cpp
index 4ba897818..ae5d47515 100644
--- a/NonLinearSolver/internalPoints/ipNonLocalDamageHyperelastic.cpp
+++ b/NonLinearSolver/internalPoints/ipNonLocalDamageHyperelastic.cpp
@@ -76,7 +76,7 @@ void IPHyperViscoElastoPlasticLocalDamage::restart()
IPHyperViscoElastoPlasticNonLocalDamage::IPHyperViscoElastoPlasticNonLocalDamage(const J2IsotropicHardening* comp,
const J2IsotropicHardening* trac,const J2IsotropicHardening* shear, const kinematicHardening* kin, const int N,
const CLengthLaw *cll, const DamageLaw *daml): IPHyperViscoElastoPlasticLocalDamage(comp,trac,shear,kin,N,daml),
- _nonlocalPlasticStrain (0)
+ _nonlocalPlasticStrain(0.), _DirreversibleEnergyDNonLocalVariable(0.)
{
ipvCL=NULL;
if(cll ==NULL) Msg::Error("IPHyperViscoElastoPlasticNonLocalDamage::IPHyperViscoElastoPlasticNonLocalDamage has no cll");
@@ -91,7 +91,8 @@ IPHyperViscoElastoPlasticNonLocalDamage::~IPHyperViscoElastoPlasticNonLocalDamag
};
IPHyperViscoElastoPlasticNonLocalDamage::IPHyperViscoElastoPlasticNonLocalDamage(const IPHyperViscoElastoPlasticNonLocalDamage &source):
- IPHyperViscoElastoPlasticLocalDamage(source), _nonlocalPlasticStrain(source._nonlocalPlasticStrain){
+ IPHyperViscoElastoPlasticLocalDamage(source), _nonlocalPlasticStrain(source._nonlocalPlasticStrain),
+ _DirreversibleEnergyDNonLocalVariable(source._DirreversibleEnergyDNonLocalVariable){
ipvCL = NULL;
if(source.ipvCL != NULL)
{
@@ -104,6 +105,7 @@ IPHyperViscoElastoPlasticNonLocalDamage& IPHyperViscoElastoPlasticNonLocalDamage
if(src != NULL)
{
_nonlocalPlasticStrain=src->_nonlocalPlasticStrain;
+ _DirreversibleEnergyDNonLocalVariable = src->_DirreversibleEnergyDNonLocalVariable;
if(src->ipvCL != NULL)
{
if (ipvCL != NULL){
@@ -125,6 +127,7 @@ void IPHyperViscoElastoPlasticNonLocalDamage::restart()
IPHyperViscoElastoPlasticLocalDamage::restart();
restartManager::restart(ipvCL);
restartManager::restart(_nonlocalPlasticStrain);
+ restartManager::restart(_DirreversibleEnergyDNonLocalVariable);
};
IPHyperViscoElastoPlasticMultipleLocalDamage::IPHyperViscoElastoPlasticMultipleLocalDamage(const J2IsotropicHardening* comp,
@@ -192,10 +195,12 @@ IPHyperViscoElastoPlasticMultipleNonLocalDamage::IPHyperViscoElastoPlasticMultip
IPHyperViscoElastoPlasticMultipleLocalDamage(comp,trac,shear,kin,N,dl),
_nonlocalEqPlasticStrain(0.),_nonlocalFailurePlasticity(0.){
ipvCL.clear();
+ _DirreversibleEnergyDNonLocalVariable.clear();
for (int i=0; i< cll.size(); i++){
IPCLength* ipvLength=NULL;
cll[i]->createIPVariable(ipvLength);
ipvCL.push_back(ipvLength);
+ _DirreversibleEnergyDNonLocalVariable.push_back(0.);
}
};
IPHyperViscoElastoPlasticMultipleNonLocalDamage::~IPHyperViscoElastoPlasticMultipleNonLocalDamage(){
@@ -208,7 +213,7 @@ IPHyperViscoElastoPlasticMultipleNonLocalDamage::IPHyperViscoElastoPlasticMultip
IPHyperViscoElastoPlasticMultipleLocalDamage(source){
_nonlocalEqPlasticStrain= source._nonlocalEqPlasticStrain;
_nonlocalFailurePlasticity = source._nonlocalFailurePlasticity;
-
+ _DirreversibleEnergyDNonLocalVariable = source._DirreversibleEnergyDNonLocalVariable;
ipvCL.clear();
for (int i=0; i< numNonLocalVariable; i++){
if(source.ipvCL[i] != NULL){
@@ -222,6 +227,7 @@ IPHyperViscoElastoPlasticMultipleNonLocalDamage& IPHyperViscoElastoPlasticMultip
if (ipsource != NULL){
_nonlocalEqPlasticStrain = ipsource->_nonlocalEqPlasticStrain;
_nonlocalFailurePlasticity = ipsource->_nonlocalFailurePlasticity;
+ _DirreversibleEnergyDNonLocalVariable = ipsource->_DirreversibleEnergyDNonLocalVariable;
for (int i=0; i< ipsource->numNonLocalVariable; i++){
if(ipsource->ipvCL[i] != NULL){
@@ -239,6 +245,9 @@ void IPHyperViscoElastoPlasticMultipleNonLocalDamage::restart(){
IPHyperViscoElastoPlasticMultipleLocalDamage::restart();
restartManager::restart(_nonlocalEqPlasticStrain);
restartManager::restart(_nonlocalFailurePlasticity);
+ for (int i=0; i< numNonLocalVariable; i++){
+ restartManager::restart(_DirreversibleEnergyDNonLocalVariable[i]);
+ }
for (int i=0; i< numNonLocalVariable; i++){
ipvCL[i]->restart();
}
diff --git a/NonLinearSolver/internalPoints/ipNonLocalDamageHyperelastic.h b/NonLinearSolver/internalPoints/ipNonLocalDamageHyperelastic.h
index 06cacfd60..6ec8f8f6d 100644
--- a/NonLinearSolver/internalPoints/ipNonLocalDamageHyperelastic.h
+++ b/NonLinearSolver/internalPoints/ipNonLocalDamageHyperelastic.h
@@ -97,6 +97,7 @@ class IPHyperViscoElastoPlasticNonLocalDamage : public IPHyperViscoElastoPlastic
protected:
IPCLength* ipvCL;
double _nonlocalPlasticStrain;
+ double _DirreversibleEnergyDNonLocalVariable;
public:
IPHyperViscoElastoPlasticNonLocalDamage(const J2IsotropicHardening* comp,
@@ -109,6 +110,9 @@ class IPHyperViscoElastoPlasticNonLocalDamage : public IPHyperViscoElastoPlastic
virtual IPHyperViscoElastoPlasticNonLocalDamage& operator=(const IPVariable &source);
virtual IPVariable* clone() const{return new IPHyperViscoElastoPlasticNonLocalDamage(*this);}
+
+ virtual const double& getDIrreversibleEnergyDNonLocalVariable() const {return _DirreversibleEnergyDNonLocalVariable;};
+ virtual double& getRefToDIrreversibleEnergyDNonLocalVariable() {return _DirreversibleEnergyDNonLocalVariable;};
virtual void restart();
virtual double getCurrentPlasticStrain() const { return _epspbarre;}
@@ -218,6 +222,7 @@ class IPHyperViscoElastoPlasticMultipleNonLocalDamage : public IPHyperViscoElast
// bareps
double _nonlocalEqPlasticStrain;
double _nonlocalFailurePlasticity;
+ std::vector<double> _DirreversibleEnergyDNonLocalVariable;
public:
IPHyperViscoElastoPlasticMultipleNonLocalDamage(const J2IsotropicHardening* comp,
@@ -228,6 +233,9 @@ class IPHyperViscoElastoPlasticMultipleNonLocalDamage : public IPHyperViscoElast
virtual IPHyperViscoElastoPlasticMultipleNonLocalDamage& operator=(const IPVariable& source);
virtual IPVariable* clone() const{return new IPHyperViscoElastoPlasticMultipleNonLocalDamage(*this);}
+
+ virtual const double& getDIrreversibleEnergyDNonLocalVariable(const int i) const {return _DirreversibleEnergyDNonLocalVariable[i];};
+ virtual double& getRefToDIrreversibleEnergyDNonLocalVariable(const int i) {return _DirreversibleEnergyDNonLocalVariable[i];};
virtual void restart();
virtual double getCurrentPlasticStrain() const { return _epspbarre;}
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
index ad7edaf4f..05d19df95 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
@@ -96,6 +96,23 @@ void mlawHyperViscoElastic::setViscoElasticData(const std::string filename){
fclose(file);
};
+double mlawHyperViscoElastic::deformationEnergy(const STensor3 &C) const
+{
+ static STensor3 logCdev;
+ STensorOperation::logSTensor3(C,_order,logCdev,NULL);
+ double trace = logCdev.trace();
+ double lnJ = 0.5*trace;
+ logCdev(0,0)-=trace/3.;
+ logCdev(1,1)-=trace/3.;
+ logCdev(2,2)-=trace/3.;
+
+ double Psy = _K*0.5*lnJ*lnJ+_mu*0.25*dot(logCdev,logCdev);
+ for (int i=0; i<_N; i++){
+ Psy += _Ki[i]*0.5*lnJ*lnJ+_Gi[i]*0.25*dot(logCdev,logCdev);
+ }
+ return Psy;;
+}
+
mlawHyperViscoElastic::mlawHyperViscoElastic(const int num,const double E,const double nu, const double rho,
const bool matrixbyPerturbation, const double pert):
@@ -380,6 +397,8 @@ void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0,
STensorOperation::multSTensor3STensor43(corKir,dlnCdC,secondPK);
STensorOperation::multSTensor3(F,secondPK,P);
// first PK
+
+ q1->_elasticEnergy=deformationEnergy(C);
if (stiff){
static STensor43 DsecondPKdC;
@@ -479,27 +498,6 @@ void mlawPowerYieldHyper::setViscosityEffect(const viscosityLaw& v, const double
_viscosity = v.clone();
};
-double mlawPowerYieldHyper::deformationEnergy(const STensor3 &C) const
-{
- double Jac= sqrt((C(0,0) * (C(1,1) * C(2,2) - C(1,2) * C(2,1)) -
- C(0,1) * (C(1,0) * C(2,2) - C(1,2) * C(2,0)) +
- C(0,2) * (C(1,0) * C(2,1) - C(1,1) * C(2,0))));
- double lnJ = log(Jac);
- STensor3 logCdev;
- STensorOperation::logSTensor3(C,_order,logCdev,NULL);
- double trace = logCdev.trace();
- logCdev(0,0)-=trace/3.;
- logCdev(1,1)-=trace/3.;
- logCdev(2,2)-=trace/3.;
-
- double Psy = _K*0.5*lnJ*lnJ+_mu*0.25*dot(logCdev,logCdev);
- for (int i=0; i<_N; i++){
- Psy += _Ki[i]*0.5*lnJ*lnJ+_Gi[i]*0.25*dot(logCdev,logCdev);
- }
- return Psy;;
-}
-
-
void mlawPowerYieldHyper::hardening(IPHyperViscoElastoPlastic* q) const{
//Msg::Error("epspCompression = %e, epspTRaction = %e, epspShear = %e",q->_epspCompression,q->_epspTraction,q->_epspShear);
if (_compression != NULL && q->_ipCompression != NULL){
@@ -531,7 +529,7 @@ void mlawPowerYieldHyper::tangent_full_perturbation(
static STensor43 tmpSTensor43;
static STensor3 Fplus, Pplus;
- static IPHyperViscoElastoPlastic q11(*q1);
+ static IPHyperViscoElastoPlastic q11(*q0);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
@@ -539,6 +537,7 @@ void mlawPowerYieldHyper::tangent_full_perturbation(
Fplus(i,j)+=_perturbationfactor;
this->predictorCorrector(Fplus,q0,&q11,Pplus,false,tmpSTensor43,tmpSTensor43,tmpSTensor43);
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;
@@ -680,7 +679,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
q1->_nup = (9.-2.*_b)/(18.+2.*_b);
double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
- double Gamma = 0.; // q1->_Gamma; // flow rule parameter
+ double Gamma = 0.; // // flow rule parameter
double PhiEq = PhiEqpr;
@@ -712,12 +711,13 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
double dDgammaDGamma = 0.;
+ double Dgamma = 0.; // eqplastic strain
if (f>0){
// plasticity
int ite = 0;
int maxite = 100; // maximal number of iters
- double Dgamma = 0.; // eqplastic strain
+
//Msg::Error("plasticity occurs f = %e",f);
@@ -829,7 +829,10 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
q1->_kirchhoff(0,0) += (ptilde);
q1->_kirchhoff(1,1) += (ptilde);
q1->_kirchhoff(2,2) += (ptilde);
- };
+ }
+ else{
+ N *= 0.;
+ }
const STensor3& KS = q1->_kirchhoff;
@@ -845,10 +848,17 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
- // update current internal variables
+ // defo energy
q1->_elasticEnergy=deformationEnergy(Ce);
-
-
+
+ // dissipation
+ double& irrEnerg = q1->getRefToIrreversibleEnergy();
+ irrEnerg = q0->irreversibleEnergy();
+ if (Dgamma > 0){
+ double dotKSN = dot(KS,N);
+ irrEnerg += Gamma*dotKSN;
+ }
+
if (stiff){
static STensor3 DpprDCepr;
static STensor43 DdevKprDCepr;
@@ -864,6 +874,9 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
static STensor43 dFpDCepr;
static STensor3 DgamaDCepr;
+ static STensor3 DtrNDCepr;
+ static STensor43 DdevNDCepr;
+ static STensor3 DGDCepr;
if (Gamma >0){
// plastic
@@ -883,7 +896,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
}
}
- static STensor3 DGDCepr;
+
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;
@@ -896,14 +909,11 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
}
}
-
- static STensor3 DtrNDCepr;
+
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
DtrNDCepr(i,j) = 2.*_b/v*DpprDCepr(i,j) - 2.*_b*ptildepr*(2.*_b*Kt)/(v*v)*DGDCepr(i,j);
-
- static STensor43 DdevNDCepr;
DdevNDCepr = (DdevKprDCepr);
DdevNDCepr *= (3./u);
for (int i=0; i<3; i++)
@@ -949,6 +959,9 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
// elastic
DgamaDCepr *= 0.;
dFpDCepr *= 0.;
+ DtrNDCepr *= 0.;
+ DdevNDCepr *= 0.;
+ DGDCepr *= 0.;
}
static STensor43 dKcorDcepr;
@@ -977,9 +990,15 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
STensor3& DgammaDF = q1->_DgammaDF;
static STensor43 DKcorDF;
- STensorOperation::multSTensor3STensor43(DgamaDCepr,CeprToF,DgammaDF);
STensorOperation::multSTensor43(dKcorDcepr,CeprToF,DKcorDF);
- STensorOperation::multSTensor43(dFpDCepr,CeprToF,dFpdF);
+ if (Gamma > 0){
+ STensorOperation::multSTensor3STensor43(DgamaDCepr,CeprToF,DgammaDF);
+ STensorOperation::multSTensor43(dFpDCepr,CeprToF,dFpdF);
+ }
+ else{
+ DgammaDF *= 0.;
+ dFpdF *= 0.;
+ }
static STensor43 DinvFpDF; //
for (int i=0; i<3; i++)
@@ -1055,7 +1074,29 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
}
}
}
-
+
+
+ // irreversible energy
+ STensor3& DirrEnergDF = q1->_DirreversibleEnergyDF;
+ 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{
+ DirrEnergDF *= 0.;
+ }
};
};
@@ -1288,7 +1329,16 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c
for(int l=0; l<3; l++)
P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
+ // defo energy
q1->_elasticEnergy=deformationEnergy(Ce);
+
+ // dissipation
+ double& irrEnerg = q1->getRefToIrreversibleEnergy();
+ irrEnerg = q0->irreversibleEnergy();
+ if (Dgamma > 0){
+ double dotKSN = dot(KS,N);
+ irrEnerg += Gamma*dotKSN;
+ }
if (stiff){
static STensor3 DpprDCepr;
@@ -1305,8 +1355,10 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c
static STensor43 dFpDCepr;
static STensor3 DgamaDCepr;
-
- if (Dgamma > 0.){
+ 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++){
@@ -1335,7 +1387,7 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c
DsigVMDCepr(i,j) = -invJ(1,0)*DPhigammaDCepr(i,j) - invJ(1,1)*DPhiSigVMDCepr(i,j);
}
- static STensor3 DGDCepr;
+
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)) +
@@ -1343,11 +1395,11 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c
}
- static STensor3 DtrNDCepr;
+
DtrNDCepr = (DgamaDCepr);
DtrNDCepr *= -Da(1);
- static STensor43 DdevNDCepr;
+
DdevNDCepr = (DdevKprDCepr);
double B = a(2)*pow(sigVM,_n-2.);
double u = 1.+3.*Ge*Gamma*_n*B;
@@ -1406,6 +1458,9 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c
else{
DgamaDCepr *= 0.;
dFpDCepr *= 0.;
+ DGDCepr *= 0.;
+ DdevNDCepr *= 0.;
+ DtrNDCepr *= 0.;
}
@@ -1518,6 +1573,27 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c
}
}
}
+ // irreversible energy
+ STensor3& DirrEnergDF = q1->_DirreversibleEnergyDF;
+ 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{
+ DirrEnergDF *= 0.;
+ }
}
};
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.h b/NonLinearSolver/materialLaw/mlawHyperelastic.h
index ebd6603f0..1974c61bc 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.h
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.h
@@ -43,6 +43,7 @@ class mlawHyperViscoElastic : public materialLaw{
std::vector<double> _gi;
protected:
+ virtual double deformationEnergy(const STensor3 &C) const ;
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,
@@ -129,9 +130,6 @@ class mlawPowerYieldHyper : public mlawHyperViscoElastic{
protected:
-
- virtual double deformationEnergy(const STensor3 &C) const ;
-
virtual void hardening(IPHyperViscoElastoPlastic* q) const;
virtual void updateEqPlasticDeformation(IPHyperViscoElastoPlastic *q1, const IPHyperViscoElastoPlastic *q0,
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.cpp
index d3bb8b050..100ab7cab 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.cpp
@@ -152,40 +152,25 @@ void mlawNonLocalDamagePowerYieldHyper::constitutive(
}
-mlawNonLocalDamagePowerYieldHyperWithFailure::mlawNonLocalDamagePowerYieldHyperWithFailure(const int num,const double E,const double nu, const double rho,
+mlawLocalDamagePowerYieldHyperWithFailure::mlawLocalDamagePowerYieldHyperWithFailure(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(){
+void mlawLocalDamagePowerYieldHyperWithFailure::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){
+
+void mlawLocalDamagePowerYieldHyperWithFailure::setDamageLaw(const DamageLaw &_damLaw){
damLaw.push_back(_damLaw.clone());
};
-mlawNonLocalDamagePowerYieldHyperWithFailure::mlawNonLocalDamagePowerYieldHyperWithFailure(const mlawNonLocalDamagePowerYieldHyperWithFailure &source):
+mlawLocalDamagePowerYieldHyperWithFailure::mlawLocalDamagePowerYieldHyperWithFailure(const mlawLocalDamagePowerYieldHyperWithFailure &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());
- }
+ for (int i=0; i< source.damLaw.size(); i++){
if(source.damLaw[i] != NULL)
{
damLaw.push_back(source.damLaw[i]->clone());
@@ -193,62 +178,45 @@ mlawNonLocalDamagePowerYieldHyperWithFailure::mlawNonLocalDamagePowerYieldHyperW
}
}
-mlawNonLocalDamagePowerYieldHyperWithFailure::~mlawNonLocalDamagePowerYieldHyperWithFailure(){
- for (int i=0; i< cLLaw.size(); i++){
- if (cLLaw[i]!= NULL) delete cLLaw[i];
+mlawLocalDamagePowerYieldHyperWithFailure::~mlawLocalDamagePowerYieldHyperWithFailure(){
+ for (int i=0; i< damLaw.size(); 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
+void mlawLocalDamagePowerYieldHyperWithFailure::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);
+ IPVariable* ipvi = new IPHyperViscoElastoPlasticMultipleLocalDamage(_compression,_traction,_shear,_kinematic,_N,damLaw);
+ IPVariable* ipv1 = new IPHyperViscoElastoPlasticMultipleLocalDamage(_compression,_traction,_shear,_kinematic,_N,damLaw);
+ IPVariable* ipv2 = new IPHyperViscoElastoPlasticMultipleLocalDamage(_compression,_traction,_shear,_kinematic,_N,damLaw);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
-void mlawNonLocalDamagePowerYieldHyperWithFailure::createIPState(IPHyperViscoElastoPlasticMultipleNonLocalDamage *ivi, IPHyperViscoElastoPlasticMultipleNonLocalDamage *iv1, IPHyperViscoElastoPlasticMultipleNonLocalDamage *iv2) const
+void mlawLocalDamagePowerYieldHyperWithFailure::createIPState(IPHyperViscoElastoPlasticMultipleLocalDamage *ivi, IPHyperViscoElastoPlasticMultipleLocalDamage *iv1, IPHyperViscoElastoPlasticMultipleLocalDamage *iv2) const
{
}
-void mlawNonLocalDamagePowerYieldHyperWithFailure::createIPVariable(IPHyperViscoElastoPlasticMultipleNonLocalDamage *&ipv,const MElement *ele,const int nbFF,const IntPt *GP, const int gpt) const
+void mlawLocalDamagePowerYieldHyperWithFailure::createIPVariable(IPHyperViscoElastoPlasticMultipleLocalDamage *&ipv,const MElement *ele,const int nbFF,const IntPt *GP, const int gpt) const
{
}
-double mlawNonLocalDamagePowerYieldHyperWithFailure::soundSpeed() const
+double mlawLocalDamagePowerYieldHyperWithFailure::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(
+void mlawLocalDamagePowerYieldHyperWithFailure::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),
+ const IPHyperViscoElastoPlasticMultipleLocalDamage *ipvprev, // array of initial internal variable
+ IPHyperViscoElastoPlasticMultipleLocalDamage *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
@@ -270,8 +238,8 @@ void mlawNonLocalDamagePowerYieldHyperWithFailure::constitutive(
//Msg::Info("enery = %e",ene);
// saturation damage always develops
- damLaw[0]->computeDamage(ipvcur->getEffectivePlasticStrain(),
- ipvprev->getEffectivePlasticStrain(),
+ damLaw[0]->computeDamage(ipvcur->getConstRefToEqPlasticStrain(),
+ ipvprev->getConstRefToEqPlasticStrain(),
ene, Fe, Fp, Peff, Cel,
ipvprev->getConstRefToIPDamage(0),ipvcur->getRefToIPDamage(0));
@@ -280,8 +248,8 @@ void mlawNonLocalDamagePowerYieldHyperWithFailure::constitutive(
// set for current critical damage
ipvcur->setCriticalDamage(ipvcur->getCriticalDamage());
// damage evolution
- damLaw[1]->computeDamage(ipvcur->getNonLocalFailurePlasticity(),
- ipvprev->getNonLocalFailurePlasticity(),
+ damLaw[1]->computeDamage(ipvcur->getFailurePlasticity(),
+ ipvprev->getFailurePlasticity(),
ene, Fe, Fp, Peff, Cel,
ipvprev->getConstRefToIPDamage(1),ipvcur->getRefToIPDamage(1));
}
@@ -306,18 +274,13 @@ void mlawNonLocalDamagePowerYieldHyperWithFailure::constitutive(
// 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++)
- {
+ 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)-=Peff(i,j)*(ipvcur->getDDamageDp(0)*dgammadF(k,l)*(1.- D2)+(1.-D1)*ipvcur->getDDamageDp(1)*dgFdF(k,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);
}
}
@@ -325,25 +288,7 @@ void mlawNonLocalDamagePowerYieldHyperWithFailure::constitutive(
}
}
}
-
-
- // -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.);
-
}
-
}
@@ -497,6 +442,20 @@ void mlawNonLocalDamagePowerYieldHyperWithFailure::constitutive(
double D = 1. - (1.- D1)*(1.-D2);
P = Peff;
P*= (1.- D);
+
+ ipvcur->_elasticEnergy = (1-D)*ene;
+
+
+
+ double dIrrPlatic = ipvcur->irreversibleEnergy();
+ dIrrPlatic -= ipvprev->irreversibleEnergy();
+ double Dprev = 1- (1-ipvprev->getDamage(0))*(1-ipvprev->getDamage(1));
+
+ // dissipation
+ double& irrEnerg = ipvcur->getRefToIrreversibleEnergy();
+ irrEnerg = ipvprev->irreversibleEnergy();
+ irrEnerg += ((1-D)*dIrrPlatic + ene*(D - Dprev));
+
if(stiff)
@@ -504,18 +463,12 @@ void mlawNonLocalDamagePowerYieldHyperWithFailure::constitutive(
// 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++)
- {
+ 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);
}
}
@@ -539,8 +492,28 @@ void mlawNonLocalDamagePowerYieldHyperWithFailure::constitutive(
// partial p partial tilde p (0 if no MFH)
dLocalVariableDNonLocalVariable.setAll(0.);
-
+
+
+ static STensor3 DIrrPlaticDF;
+ DIrrPlaticDF = ipvcur->getConstRefToDIrreversibleEnergyDF();
+
+ STensor3& DirrEnergDF = ipvcur->getRefToDIrreversibleEnergyDF();
+ for(int i=0;i<3;i++){
+ for(int j=0;j<3;j++){
+ DirrEnergDF(i,j) = (1-D)*DIrrPlaticDF(i,j);
+ for(int k=0;k<3;k++){
+ for(int l=0;l<3;l++){
+ for(int m=0; m<3; m++){
+ DirrEnergDF(i,j) += (D - Dprev)*Peff(k,m)*Fp(l,m)*dFedF(k,l,i,j);
+ };
+ }
+ }
+ }
+ }
+
+ ipvcur->getRefToDIrreversibleEnergyDNonLocalVariable(0) = (ene- dIrrPlatic)*(ipvcur->getDDamageDp(0)*(1.- D2));
+ ipvcur->getRefToDIrreversibleEnergyDNonLocalVariable(1) = (ene- dIrrPlatic)*(ipvcur->getDDamageDp(1)*(1.- D1));
+
}
-
}
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.h b/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.h
index 0e40899c5..bac77a429 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.h
+++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageHyperelastic.h
@@ -108,19 +108,9 @@ class mlawLocalDamagePowerYieldHyperWithFailure : public mlawPowerYieldHyperWith
virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}; // this law is initialized so nothing to do
virtual double soundSpeed() const; // default but you can redefine it for your case
- virtual const std::vector<CLengthLaw*>& getCLengthLaw() const {return cLLaw; };
virtual const std::vector<DamageLaw*>& getDamageLaw() const {return damLaw; };
// specific function
public:
- virtual void 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 IPHyperViscoElastoPlasticMultipleLocalDamage *q0, // array of initial internal variable
- IPHyperViscoElastoPlasticMultipleLocalDamage *q1, // updated array of internal variable (in ipvcur on output),
- STensor43 &Tangent, // constitutive tangents (output)
- const bool stiff // if true compute the tangents
- ) const;
virtual void constitutive(
const STensor3& F0, // initial deformation gradient (input @ time n)
@@ -129,9 +119,6 @@ class mlawLocalDamagePowerYieldHyperWithFailure : public mlawPowerYieldHyperWith
const IPHyperViscoElastoPlasticMultipleLocalDamage *q0, // array of initial internal variable
IPHyperViscoElastoPlasticMultipleLocalDamage *q1, // 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;
#endif // SWIG
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2Hyper.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2Hyper.cpp
index 267d29500..346c5bd3d 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2Hyper.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageJ2Hyper.cpp
@@ -282,7 +282,19 @@ void mlawNonLocalDamageJ2Hyper::constitutive(
}
else if (this->_pfApproxMethod == mlawJ2linear::Implicit){
double H = ipvcur->getConstRefToIPJ2IsotropicHardening().getDR();
- DdissEnergDF = Peff;
+ for(int i=0;i<3;i++){
+ for(int j=0;j<3;j++){
+ DdissEnergDF(i,j) = 0.;
+ for(int k=0;k<3;k++){
+ for(int l=0;l<3;l++){
+ for (int m=0; m<3; m++){
+ DdissEnergDF(i,j) += Peff(k,m)*Fp(l,m)*dFedF(k,l,i,j);
+ }
+ }
+ }
+ }
+ }
+
DdissEnergDF*= (D-D0);
double DdissDD = ene - Sy*dp;
DdissEnergDF.daxpy(DDamageDF,DdissDD);
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index 07ca904bc..32e693885 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -1058,7 +1058,10 @@ void nonLinearMechSolver::oneStepPreSolvePathFollowing(const double curtime, con
currentPathFollowingIncr = _loadStep;
}
else if (pathsys->getControlType() == pathFollowingSystemBase::LOCAL_CONTROL){
- double ff = double(_numNROptimalLocal)/double(_numNROptimal);
+ double ff = 1.;
+ if (_timeStepAdaptation){
+ ff = double(_numNROptimalLocal)/double(_numNROptimal);
+ }
_localStep = _localStepPrev*_timeStepAdaptationFractor*ff;
currentPathFollowingIncr = _localStep;
}
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index 1bfa746e1..f31dad5de 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -362,6 +362,13 @@ public:
if (_ipViscoElastic)
_ipViscoElastic->setLocation(loc);
};
+
+ // for path following based on irreversibe energt
+ virtual double irreversibleEnergy() const {return _ipViscoElastic->irreversibleEnergy();};
+ virtual double& getRefToIrreversibleEnergy() {return _ipViscoElastic->getRefToIrreversibleEnergy();};
+
+ virtual const STensor3& getConstRefToDIrreversibleEnergyDDeformationGradient() const {return _ipViscoElastic->getConstRefToDIrreversibleEnergyDF();};
+ virtual STensor3& getRefToDIrreversibleEnergyDDeformationGradient() {return _ipViscoElastic->getRefToDIrreversibleEnergyDF();};
virtual double defoEnergy() const;
@@ -390,11 +397,16 @@ class HyperViscoElastoPlasticdG3DIPVariable : public dG3DIPVariable{
virtual double get(const int i) const;
+ // for path following based on irreversibe energt
+ virtual double irreversibleEnergy() const {return _ipViscoElastoPlastic->irreversibleEnergy();};
+ virtual double& getRefToIrreversibleEnergy() {return _ipViscoElastoPlastic->getRefToIrreversibleEnergy();};
+
+ virtual const STensor3& getConstRefToDIrreversibleEnergyDDeformationGradient() const {return _ipViscoElastoPlastic->getConstRefToDIrreversibleEnergyDF();};
+ virtual STensor3& getRefToDIrreversibleEnergyDDeformationGradient() {return _ipViscoElastoPlastic->getRefToDIrreversibleEnergyDF();};
+
virtual double defoEnergy() const {return _ipViscoElastoPlastic->defoEnergy();};
virtual double plasticEnergy() const {return _ipViscoElastoPlastic->plasticEnergy();};
- // for path following based on irreversibe energt
- virtual double irreversibleEnergy() const {return _ipViscoElastoPlastic->irreversibleEnergy();};
virtual IPVariable* clone() const {return new HyperViscoElastoPlasticdG3DIPVariable(*this);};
virtual void restart();
diff --git a/dG3D/src/nonLocalDamageHyperelasticDG3DIPVariable.h b/dG3D/src/nonLocalDamageHyperelasticDG3DIPVariable.h
index 162afa477..c3aae4e35 100644
--- a/dG3D/src/nonLocalDamageHyperelasticDG3DIPVariable.h
+++ b/dG3D/src/nonLocalDamageHyperelasticDG3DIPVariable.h
@@ -54,6 +54,32 @@ class nonLocalDamageHyperelasticDG3DIPVariable: public nonLocalDamageDG3DIPVaria
else
Msg::Fatal("the non-local damge id = %d is not defined",idex);
}
+
+ // for path following based on irreversibe energt
+ virtual double irreversibleEnergy() const {return _nldHyperipv->irreversibleEnergy();};
+ virtual double& getRefToIrreversibleEnergy() {return _nldHyperipv->getRefToIrreversibleEnergy();};
+
+ virtual const STensor3& getConstRefToDIrreversibleEnergyDDeformationGradient() const
+ {return _nldHyperipv->getConstRefToDIrreversibleEnergyDF();};
+ virtual STensor3& getRefToDIrreversibleEnergyDDeformationGradient()
+ {return _nldHyperipv->getRefToDIrreversibleEnergyDF();};
+
+ virtual const double& getConstRefToDIrreversibleEnergyDNonLocalVariable(const int index) const
+ {
+ if (index == 0)
+ {
+ return _nldHyperipv->getDIrreversibleEnergyDNonLocalVariable();
+ }
+ else{Msg::Fatal("the non-local variable %d is not defined",index);}
+ };
+ virtual double& getRefToDIrreversibleEnergyDDNonLocalVariable(const int index)
+ {
+ if (index == 0)
+ {
+ return _nldHyperipv->getRefToDIrreversibleEnergyDNonLocalVariable();
+ }
+ else{Msg::Fatal("the non-local variable %d is not defined",index);}
+ };
virtual IPVariable* clone() const {return new nonLocalDamageHyperelasticDG3DIPVariable(*this);};
virtual void restart();
@@ -97,6 +123,26 @@ class nonLocalDamageHyperelasticDG3DIPVariableWithFailure: public nonLocalDamage
{
return _nldHyperipv->getConstRefToCharacteristicLength(idex);
}
+
+ // for path following based on irreversibe energt
+ virtual double irreversibleEnergy() const {return _nldHyperipv->irreversibleEnergy();};
+ virtual double& getRefToIrreversibleEnergy() {return _nldHyperipv->getRefToIrreversibleEnergy();};
+
+ virtual const STensor3& getConstRefToDIrreversibleEnergyDDeformationGradient() const {
+ return _nldHyperipv->getConstRefToDIrreversibleEnergyDF();
+ };
+ virtual STensor3& getRefToDIrreversibleEnergyDDeformationGradient(){
+ return _nldHyperipv->getRefToDIrreversibleEnergyDF();
+ };
+
+ virtual const double& getConstRefToDIrreversibleEnergyDNonLocalVariable(const int index) const
+ {
+ return _nldHyperipv->getDIrreversibleEnergyDNonLocalVariable(index);
+ };
+ virtual double& getRefToDIrreversibleEnergyDDNonLocalVariable(const int index)
+ {
+ return _nldHyperipv->getRefToDIrreversibleEnergyDNonLocalVariable(index);
+ };
virtual IPVariable* clone() const {return new nonLocalDamageHyperelasticDG3DIPVariableWithFailure(*this);};
virtual void restart();
diff --git a/dG3D/src/pathFollowingTerms.cpp b/dG3D/src/pathFollowingTerms.cpp
index deeee1467..17ee48375 100644
--- a/dG3D/src/pathFollowingTerms.cpp
+++ b/dG3D/src/pathFollowingTerms.cpp
@@ -441,7 +441,6 @@ void dG3DNonLocalDissipationPathFollowingBulkLinearTerm::get(MElement *ele,int n
int numNonLocalVar = ipv->getNumberNonLocalVariable();
for (int inl =0; inl < numNonLocalVar; inl++){
const double & DdamEnergyDNonLocalVar = ipv->getConstRefToDIrreversibleEnergyDNonLocalVariable(inl);
-
for (int k=0; k< nbFF; k++){
m(k+inl*nbFF+3*nbFF) += DdamEnergyDNonLocalVar*Vals[k+ inl*nbFF+3*nbFF]*detJ*weight;
}
--
GitLab