diff --git a/NonLinearSolver/internalPoints/ipNucleation.h b/NonLinearSolver/internalPoints/ipNucleation.h
index cc3273b5a0f70de4e4a4bc80beb531169603fe24..7382f614e945067c19d6f2446e9806eebb3bfd56 100644
--- a/NonLinearSolver/internalPoints/ipNucleation.h
+++ b/NonLinearSolver/internalPoints/ipNucleation.h
@@ -39,6 +39,12 @@ class IPNucleation{
virtual double& getRefToNucleatedPorosity() {return nucleatedfV;};
virtual double& getRefToNucleatedPorosityIncrement() {return DeltafV;};
virtual double& getRefToDNucleatedPorosityDMatrixPlasticDeformation() {return DDeltafVDHatP;};
+
+ virtual void resetNucleation(const IPNucleation& ipv) {
+ nucleatedfV = ipv.nucleatedfV;
+ DeltafV = 0.;
+ DDeltafVDHatP = 0.;
+ };
};
class IPLinearNucleationRate : public IPNucleation{
diff --git a/NonLinearSolver/internalPoints/ipVoidState.cpp b/NonLinearSolver/internalPoints/ipVoidState.cpp
index 49020c96a7f99a1b2031bb37d4895e8af0b6b58a..c3cfce6975e39a903180f8fa5b28c9938ed7d56e 100644
--- a/NonLinearSolver/internalPoints/ipVoidState.cpp
+++ b/NonLinearSolver/internalPoints/ipVoidState.cpp
@@ -15,16 +15,21 @@
#include <math.h>
IPSphericalVoidState::IPSphericalVoidState(const double fV0, const double lambda0):
- Chi(pow(1.5*fV0*lambda0,1./3.)),DChiDyieldfV(0.),DChiDhatQ(0.),DChiDhatD(0.),DChiDhatP(0.),lambda(lambda0){
+ Chi(pow(1.5*fV0*lambda0,1./3.)),DChiDyieldfV(0.),DChiDhatQ(0.),DChiDhatD(0.),DChiDhatP(0.),lambda(lambda0)
+{
+ roughChi = Chi;
}
IPSphericalVoidState::IPSphericalVoidState(const IPSphericalVoidState& src): IPVoidState(src),
- Chi(src.Chi),DChiDyieldfV(src.DChiDyieldfV),DChiDhatQ(src.DChiDhatQ),DChiDhatD(src.DChiDhatD),DChiDhatP(src.DChiDhatP),lambda(src.lambda){}
+ Chi(src.Chi),DChiDyieldfV(src.DChiDyieldfV),DChiDhatQ(src.DChiDhatQ),DChiDhatD(src.DChiDhatD),DChiDhatP(src.DChiDhatP),lambda(src.lambda),
+ roughChi(src.roughChi){}
IPSphericalVoidState& IPSphericalVoidState::operator =(const IPVoidState& src){
IPVoidState::operator =(src);
const IPSphericalVoidState* psrc = dynamic_cast<const IPSphericalVoidState*>(&src);
- if (psrc != NULL){
+ if (psrc != NULL)
+ {
+ roughChi = psrc->roughChi;
Chi = psrc->Chi;
DChiDyieldfV = psrc->DChiDyieldfV;
DChiDhatQ = psrc->DChiDhatQ;
@@ -36,6 +41,7 @@ IPSphericalVoidState& IPSphericalVoidState::operator =(const IPVoidState& src){
};
void IPSphericalVoidState::restart(){
+ restartManager::restart(roughChi);
restartManager::restart(Chi);
restartManager::restart(DChiDyieldfV);
restartManager::restart(DChiDhatD);
@@ -49,16 +55,20 @@ IPSpheroidalVoidState::IPSpheroidalVoidState(const double fV0, const double lamb
Chi(pow(3.*gamma0*fV0*lambda0/W0,1./3.)),DChiDyieldfV(0.),DChiDhatQ(0.),DChiDhatD(0.),DChiDhatP(0.),
W(W0),DWDyieldfV(0.),DWDhatQ(0.),DWDhatD(0.),DWDhatP(0.),
gamma(gamma0),DgammaDyieldfV(0.),DgammaDhatQ(0.),DgammaDhatD(0.),DgammaDhatP(0.),
- lambda(lambda0){}
+ lambda(lambda0)
+{
+ roughChi = Chi;
+}
IPSpheroidalVoidState::IPSpheroidalVoidState(const IPSpheroidalVoidState& src):IPVoidState(src),
Chi(src.Chi),DChiDyieldfV(src.DChiDyieldfV),DChiDhatQ(src.DChiDhatQ),DChiDhatD(src.DChiDhatD),DChiDhatP(src.DChiDhatP),
W(src.W),DWDyieldfV(src.DWDyieldfV),DWDhatQ(src.DWDhatQ),DWDhatD(src.DWDhatD),DWDhatP(src.DWDhatP),
gamma(src.gamma),DgammaDyieldfV(src.DgammaDyieldfV),DgammaDhatQ(src.DgammaDhatQ),DgammaDhatD(src.DgammaDhatD),DgammaDhatP(src.DgammaDhatP),
- lambda(src.lambda){};
+ lambda(src.lambda),roughChi(src.roughChi){};
IPSpheroidalVoidState& IPSpheroidalVoidState::operator =(const IPVoidState& src){
IPVoidState::operator =(src);
const IPSpheroidalVoidState* psrc = dynamic_cast<const IPSpheroidalVoidState*>(&src);
if (psrc != NULL){
+ roughChi = psrc->roughChi;
Chi = psrc->Chi;
DChiDyieldfV = psrc->DChiDyieldfV;
DChiDhatQ = psrc->DChiDhatQ;
@@ -80,7 +90,7 @@ IPSpheroidalVoidState& IPSpheroidalVoidState::operator =(const IPVoidState& src)
};
void IPSpheroidalVoidState::restart(){
-
+ restartManager::restart(roughChi);
restartManager::restart(Chi);
restartManager::restart(DChiDyieldfV);
restartManager::restart(DChiDhatD);
diff --git a/NonLinearSolver/internalPoints/ipVoidState.h b/NonLinearSolver/internalPoints/ipVoidState.h
index aed88ceaebee31ef5b6005f4ce81d904ac2e010d..bbeac3264220b54b68f6ace6b6a368cef260ab4e 100644
--- a/NonLinearSolver/internalPoints/ipVoidState.h
+++ b/NonLinearSolver/internalPoints/ipVoidState.h
@@ -27,18 +27,21 @@ class IPVoidState{
virtual void restart() = 0;
// Geometrical parameters and their derivatives
+ virtual double getRoughVoidLigamentRatio() const =0;
+ virtual double& getRefToRoughVoidLigamentRatio() =0;
+
virtual double getVoidLigamentRatio() const = 0;
virtual double getDVoidLigamentRatioDYieldPorosity() const = 0;
virtual double getDVoidLigamentRatioDVolumetricPlasticDeformation() const =0;
virtual double getDVoidLigamentRatioDDeviatoricPlasticDeformation() const =0;
virtual double getDVoidLigamentRatioDMatrixPlasticDeformation() const =0;
-
+
virtual double& getRefToVoidLigamentRatio() =0;
virtual double& getRefToDVoidLigamentRatioDYieldPorosity() =0;
virtual double& getRefToDVoidLigamentRatioDVolumetricPlasticDeformation() =0;
virtual double& getRefToDVoidLigamentRatioDDeviatoricPlasticDeformation() =0;
virtual double& getRefToDVoidLigamentRatioDMatrixPlasticDeformation() =0;
-
+
virtual double getVoidAspectRatio() const =0;
virtual double getDVoidAspectRatioDYieldPorosity() const =0;
virtual double getDVoidAspectRatioDVolumetricPlasticDeformation() const =0;
@@ -81,12 +84,16 @@ class IPNoVoidState : public IPVoidState{
virtual void restart() {};
// Geometrical parameters and their derivatives
+ virtual double getRoughVoidLigamentRatio() const {return 0.;};
+ virtual double& getRefToRoughVoidLigamentRatio() {Msg::Error("IPNoVoidState::getRefToRoughVoidLigamentRatio cannot be used");};
+
virtual double getVoidLigamentRatio() const {return 0.;};
virtual double getDVoidLigamentRatioDYieldPorosity() const {return 0.;};
virtual double getDVoidLigamentRatioDVolumetricPlasticDeformation() const {return 0.;};;
virtual double getDVoidLigamentRatioDDeviatoricPlasticDeformation() const {return 0.;};;
virtual double getDVoidLigamentRatioDMatrixPlasticDeformation() const {return 0.;};;
-
+
+
virtual double& getRefToVoidLigamentRatio() {Msg::Error("IPNoVoidState::getRefToVoidLigamentRatio cannot be used");};
virtual double& getRefToDVoidLigamentRatioDYieldPorosity() {Msg::Error("IPNoVoidState::getRefToDVoidLigamentRatioDYieldPorosity cannot be used");};
virtual double& getRefToDVoidLigamentRatioDVolumetricPlasticDeformation() {Msg::Error("IPNoVoidState::getRefToDVoidLigamentRatioDVolumetricPlasticDeformation cannot be used");};
@@ -127,6 +134,7 @@ class IPSphericalVoidState : public IPVoidState{
protected:
// Geometrical parameters and their derivatives
// Ligament ratio
+ double roughChi;
double Chi;
double DChiDyieldfV;
double DChiDhatD;
@@ -145,6 +153,9 @@ class IPSphericalVoidState : public IPVoidState{
virtual IPVoidState* clone() const {return new IPSphericalVoidState(*this);};
virtual void restart();
+ virtual double getRoughVoidLigamentRatio() const {return roughChi;};
+ virtual double& getRefToRoughVoidLigamentRatio() {return roughChi;};
+
// Geometrical parameters and their derivatives
virtual double getVoidLigamentRatio() const {return Chi;};
virtual double getDVoidLigamentRatioDYieldPorosity() const {return DChiDyieldfV;};
@@ -190,6 +201,7 @@ class IPSpheroidalVoidState : public IPVoidState{
protected:
// Geometrical parameters and their derivatives
// Ligament ratio
+ double roughChi;
double Chi;
double DChiDyieldfV;
double DChiDhatD;
@@ -221,6 +233,9 @@ class IPSpheroidalVoidState : public IPVoidState{
virtual IPVoidState* clone() const {return new IPSpheroidalVoidState(*this);};
virtual void restart();
+ virtual double getRoughVoidLigamentRatio() const {return roughChi;};
+ virtual double& getRefToRoughVoidLigamentRatio() {return roughChi;};
+
// Geometrical parameters and their derivatives
virtual double getVoidLigamentRatio() const {return Chi;};
virtual double getDVoidLigamentRatioDYieldPorosity() const {return DChiDyieldfV;};
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp
index 4c2ef3395046814be0c5fa210bb49f87555824b9..b280d93b88eade0c8915aee0ccd193d9b75da73d 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp
@@ -29,9 +29,6 @@ mlawNonLocalDamageGurson::mlawNonLocalDamageGurson(const int num,const double E,
if (_correctedRegularizedFunction != NULL) delete _correctedRegularizedFunction;
_correctedRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 1.e-8, 0.98*_fV_failure,0.999*_fV_failure);
- if (_localRegularizedFunction != NULL) delete _localRegularizedFunction;
- _localRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 1.e-8, 0.9*_localfVFailure,0.999*_localfVFailure);
-
_failedTol = 0.995;
_temFunc_q1 = new constantScalarFunction(1.);
@@ -249,15 +246,13 @@ void mlawNonLocalDamageGurson::plasticFlowDerivatives( fullVector<double>& gradN
};
-void mlawNonLocalDamageGurson::computeResidual(SVector3& res, STensor3& J,
+void mlawNonLocalDamageGurson::computeResidual_NonLocalPorosity(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
+ const double yieldfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
- const double* T
- ) const{
-
+ const double* T) const
+{
double f0 = q1->getInitialPorosity();
double R0 = _j2IH->getYield0(); // reference value
@@ -291,44 +286,41 @@ void mlawNonLocalDamageGurson::computeResidual(SVector3& res, STensor3& J,
computeYieldDerivatives(gradf,kcorEq,pcor,R,yieldfV,T);
plasticFlowDerivatives(gradNs,gradNv,kcorEq,pcor,R,yieldfV,T);
- double Dres0DfV = gradf(3)*DyieldfVDfV;
- double Dres1DfV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DyieldfVDfV;
- double Dres2DfV = 0.;
-
// Dres0DDeltaHatD
- J(0,0) = gradf(0)*DkcorEqDDeltaHatD + Dres0DfV*DDeltafVDDeltaHatD;
+ J(0,0) = gradf(0)*DkcorEqDDeltaHatD;
// Dres0DDeltaHatQ
- J(0,1) = gradf(1)*DpcorDDeltaHatQ + Dres0DfV*DDeltafVDDeltaHatQ;
+ J(0,1) = gradf(1)*DpcorDDeltaHatQ;
// Dres0DDeltaHatP
- J(0,2) = gradf(2)*H + Dres0DfV*DDeltafVDDeltaHatP;
+ J(0,2) = gradf(2)*H;
// Dres1DDeltaHatD
- J(1,0) = R0*(Nv + DeltaHatD*gradNv(0)*DkcorEqDDeltaHatD - DeltaHatQ*gradNs(0)*DkcorEqDDeltaHatD) + Dres1DfV*DDeltafVDDeltaHatD;
+ J(1,0) = R0*(Nv + DeltaHatD*gradNv(0)*DkcorEqDDeltaHatD - DeltaHatQ*gradNs(0)*DkcorEqDDeltaHatD);
// Dres1DDeltaHatQ
- J(1,1) = R0*(DeltaHatD*gradNv(1)*DpcorDDeltaHatQ- Ns - DeltaHatQ*gradNs(1)*DpcorDDeltaHatQ) + Dres1DfV*DDeltafVDDeltaHatQ;
+ J(1,1) = R0*(DeltaHatD*gradNv(1)*DpcorDDeltaHatQ- Ns - DeltaHatQ*gradNs(1)*DpcorDDeltaHatQ);
// Dres1DDeltaHatP
- J(1,2) = R0*(DeltaHatD*gradNv(2) - DeltaHatQ*gradNs(2))*H + Dres1DfV*DDeltafVDDeltaHatP;
+ J(1,2) = R0*(DeltaHatD*gradNv(2) - DeltaHatQ*gradNs(2))*H;
// Dres2DDeltaHatd
- J(2,0) = (-kcorEq - DkcorEqDDeltaHatD*DeltaHatD)/R0 + Dres2DfV*DDeltafVDDeltaHatD;
+ J(2,0) = (-kcorEq - DkcorEqDDeltaHatD*DeltaHatD)/R0;
// Dres2DDeltaHatq
- J(2,1) = (-pcor - DpcorDDeltaHatQ*DeltaHatQ)/R0 + Dres2DfV*DDeltafVDDeltaHatQ;
+ J(2,1) = (-pcor - DpcorDDeltaHatQ*DeltaHatQ)/R0;
// Dres2DDeltaHatp
- J(2,2) = (1-f0)*(R + H*DeltaHatP)/R0+ Dres2DfV*DDeltafVDDeltaHatP;
+ J(2,2) = (1-f0)*(R + H*DeltaHatP)/R0;
};
-void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
- const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const double DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
- const double* T, const double* DRDT,
- const STensor3* DdevKcorprDT , const double* DpcorprDT, const double* DDeltafVDT,
- double *dres0dT, double *dres1dT, double *dres2dT
- ) const{
+void mlawNonLocalDamageGurson::computeDResidual_NonLocalPorosity(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
+ double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
+ const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
+ const double yieldfV, const STensor3& DyieldfVDKcorpr, const double DyieldfVDtildefV,
+ const STensor43& DKcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
+ const double* T , const double* DRDT,
+ const STensor3* DdevKcorprDT , const double* DpcorprDT, const STensor3* DKcorprDT,
+ double *dres0dT, double *dres1dT, double *dres2dT
+ ) const
+{
double f0 = q1->getInitialPorosity();
double R0 = _j2IH->getYield0(); // reference value
@@ -339,6 +331,9 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dEpr, STensor3 &d
KT=bulkModulus(*T);
muT=shearModulus(*T);
}
+
+ static STensor3 DyieldfVDEpr;
+ STensorOperation::multSTensor3STensor43(DyieldfVDKcorpr,DKcorprDEpr,DyieldfVDEpr);
static STensor3 DkcorEqprDEpr;
STensorOperation::multSTensor3STensor43(devKcorpr,DdevKcorprDEpr,DkcorEqprDEpr);
@@ -359,21 +354,21 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dEpr, STensor3 &d
// compute Dres0DEpr
double Dres0DkcorEq = gradf(0);
double Dres0Dpcor = gradf(1);
- double Dres0DDeltafV = gradf(3)*DyieldfVDfV;
+ double Dres0DyieldfV = gradf(3);
dres0dEpr = DkcorEqDEpr;
dres0dEpr*= Dres0DkcorEq;
dres0dEpr.daxpy(DpcorDEpr,Dres0Dpcor);
- dres0dEpr.daxpy(DDeltafVDEpr,Dres0DDeltafV);
+ dres0dEpr.daxpy(DyieldfVDEpr,Dres0DyieldfV);
// compute Dres1DEpr
double Dres1DkcorEq = R0*(DeltaHatD*gradNv(0)- DeltaHatQ*gradNs(0));
double Dres1Dpcor = R0*(DeltaHatD*gradNv(1)- DeltaHatQ*gradNs(1));
- double Dres1DDeltafV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DyieldfVDfV;
+ double Dres1DyieldfV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3));
dres1dEpr = DkcorEqDEpr;
dres1dEpr *= Dres1DkcorEq;
dres1dEpr.daxpy(DpcorDEpr,Dres1Dpcor);
- dres1dEpr.daxpy(DDeltafVDEpr,Dres1DDeltafV);
+ dres1dEpr.daxpy(DyieldfVDEpr,Dres1DyieldfV);
// compute Dres2DEpr
double Dres2DkcorEq = (-DeltaHatD)/R0;
@@ -392,31 +387,41 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dEpr, STensor3 &d
double dKdT=bulkModulusDerivative(*T);
double dmudT=shearModulusDerivative(*T);
- static double DkcorEqDT, DpcorDT;
+ static double DkcorEqDT, DpcorDT, DyieldfVDT;
DkcorEqDT = STensorOperation::doubledot(devKcorpr,*DdevKcorprDT)*1.5/kcorEqpr - 3.*dmudT*DeltaHatD;
DpcorDT = *DpcorprDT -dKdT* DeltaHatQ;
-
-
+ DyieldfVDT = STensorOperation::doubledot(DyieldfVDKcorpr,*DKcorprDT);
+
double partDres0DT = gradf(2)*(*DRDT) + gradf(4); // yield temperature-dependent
- *dres0dT = Dres0DkcorEq*DkcorEqDT + Dres0Dpcor*DpcorDT + partDres0DT;
+ *dres0dT = Dres0DkcorEq*DkcorEqDT + Dres0Dpcor*DpcorDT + Dres0DyieldfV*DyieldfVDT + partDres0DT;
double partDres1DT = R0*(DeltaHatD*gradNv(2)- DeltaHatQ*gradNs(2))*(*DRDT) + R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4));
- *dres1dT = Dres1DkcorEq*DkcorEqDT + Dres1Dpcor*DpcorDT + partDres1DT;
+ *dres1dT = Dres1DkcorEq*DkcorEqDT + Dres1Dpcor*DpcorDT + Dres1DyieldfV*DyieldfVDT + partDres1DT;
double partDres2DT = ((1-f0)*(*DRDT)*DeltaHatP)/R0;
- *dres2dT = Dres2DkcorEq*DkcorEqDT + Dres2Dpcor*DpcorDT+ partDres2DT;
+ *dres2dT = Dres2DkcorEq*DkcorEqDT + Dres2Dpcor*DpcorDT + partDres2DT;
}
};
+void mlawNonLocalDamageGurson::computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double* T
+ ) const{
+ computeResidual_NonLocalPorosity(res,J,kcorEqpr,pcorpr,R,H,yieldfV,q0,q1,DeltaHatD,DeltaHatQ,DeltaHatP,T);
+};
+
void mlawNonLocalDamageGurson::computeDResidual_multipleNonLocalVariables(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- std::vector<double> &dres0dtildefV, std::vector<double> &dres1dtildefV, std::vector<double> &dres2dtildefV,
+ std::vector<double> &Dres0DNonlocalVars, std::vector<double> &Dres1DNonlocalVars, std::vector<double> &Dres2DNonlocalVars,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const std::vector<double>& DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDKcorpr, const std::vector<double>& DyieldfVDNonlocalVars,
+ const STensor43& DKcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T, const double* DRDT,
- const STensor3* DdevKcorprDT , const double* DpcorprDT, const double* DDeltafVDT,
+ const STensor3* DdevKcorprDT , const double* DpcorprDT, const STensor3* DKcorprVDT,
double *dres0dT, double *dres1dT, double *dres2dT
) const{
@@ -429,6 +434,9 @@ void mlawNonLocalDamageGurson::computeDResidual_multipleNonLocalVariables(STenso
KT=bulkModulus(*T);
muT=shearModulus(*T);
}
+
+ static STensor3 DyieldfVDEpr;
+ STensorOperation::multSTensor3STensor43(DyieldfVDKcorpr,DKcorprDEpr,DyieldfVDEpr);
static STensor3 DkcorEqprDEpr;
STensorOperation::multSTensor3STensor43(devKcorpr,DdevKcorprDEpr,DkcorEqprDEpr);
@@ -449,21 +457,21 @@ void mlawNonLocalDamageGurson::computeDResidual_multipleNonLocalVariables(STenso
// compute Dres0DEpr
double Dres0DkcorEq = gradf(0);
double Dres0Dpcor = gradf(1);
- double Dres0DDeltafV = gradf(3)*DyieldfVDfV;
+ double Dres0DyieldfV = gradf(3);
dres0dEpr = DkcorEqDEpr;
dres0dEpr*= Dres0DkcorEq;
dres0dEpr.daxpy(DpcorDEpr,Dres0Dpcor);
- dres0dEpr.daxpy(DDeltafVDEpr,Dres0DDeltafV);
+ dres0dEpr.daxpy(DyieldfVDEpr,Dres0DyieldfV);
// compute Dres1DEpr
double Dres1DkcorEq = R0*(DeltaHatD*gradNv(0)- DeltaHatQ*gradNs(0));
double Dres1Dpcor = R0*(DeltaHatD*gradNv(1)- DeltaHatQ*gradNs(1));
- double Dres1DDeltafV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DyieldfVDfV;
+ double Dres1DyieldfV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3));
dres1dEpr = DkcorEqDEpr;
dres1dEpr *= Dres1DkcorEq;
dres1dEpr.daxpy(DpcorDEpr,Dres1Dpcor);
- dres1dEpr.daxpy(DDeltafVDEpr,Dres1DDeltafV);
+ dres1dEpr.daxpy(DyieldfVDEpr,Dres1DyieldfV);
// compute Dres2DEpr
double Dres2DkcorEq = (-DeltaHatD)/R0;
@@ -475,9 +483,9 @@ void mlawNonLocalDamageGurson::computeDResidual_multipleNonLocalVariables(STenso
// compute nonlocal derivative
for (int i=0; i< _numNonLocalVar; i++){
- dres0dtildefV[i] = gradf(3)*DyieldfVDtildefV[i];
- dres1dtildefV[i] = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DyieldfVDtildefV[i];
- dres2dtildefV[i] = 0.;
+ Dres0DNonlocalVars[i] = gradf(3)*DyieldfVDNonlocalVars[i];
+ Dres1DNonlocalVars[i] = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DyieldfVDNonlocalVars[i];
+ Dres2DNonlocalVars[i] = 0.;
}
if (isThermomechanicallyCoupled()){
@@ -487,19 +495,88 @@ void mlawNonLocalDamageGurson::computeDResidual_multipleNonLocalVariables(STenso
static double DkcorEqDT, DpcorDT;
DkcorEqDT = STensorOperation::doubledot(devKcorpr,*DdevKcorprDT)*1.5/kcorEqpr - 3.*dmudT*DeltaHatD;
DpcorDT = *DpcorprDT -dKdT* DeltaHatQ;
-
+
+ double DyieldfVDT = STensorOperation::doubledot(DyieldfVDKcorpr,*DKcorprVDT);
double partDres0DT = gradf(2)*(*DRDT) + gradf(4); // yield temperature-dependent
- *dres0dT = Dres0DkcorEq*DkcorEqDT + Dres0Dpcor*DpcorDT + partDres0DT;
+ *dres0dT = Dres0DkcorEq*DkcorEqDT + Dres0Dpcor*DpcorDT + Dres0DyieldfV*DyieldfVDT + partDres0DT;
double partDres1DT = R0*(DeltaHatD*gradNv(2)- DeltaHatQ*gradNs(2))*(*DRDT) + R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4));
- *dres1dT = Dres1DkcorEq*DkcorEqDT + Dres1Dpcor*DpcorDT + partDres1DT;
+ *dres1dT = Dres1DkcorEq*DkcorEqDT + Dres1Dpcor*DpcorDT + Dres0DyieldfV*DyieldfVDT + partDres1DT;
double partDres2DT = ((1-f0)*(*DRDT)*DeltaHatP)/R0;
*dres2dT = Dres2DkcorEq*DkcorEqDT + Dres2Dpcor*DpcorDT+ partDres2DT;
}
};
+void mlawNonLocalDamageGurson::computeResidualLocal(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV, const double DyieldfVDfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
+ const double* T) const
+{
+ double f0 = q1->getInitialPorosity();
+ double R0 = _j2IH->getYield0(); // reference value
+
+ double KT= _K;
+ double muT = _mu;
+ if (isThermomechanicallyCoupled()){
+ KT = bulkModulus(*T);
+ muT =shearModulus(*T);
+ };
+
+ // corrected value
+ double kcorEq = kcorEqpr - 3.*muT*DeltaHatD;
+ double pcor = pcorpr - KT*DeltaHatQ;
+
+ // yield
+ res(0) = yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
+
+ double Ns, Nv; // plastic flow components
+ plasticFlow(Ns,Nv,kcorEq,pcor,R,yieldfV,T);
+ // normality
+ res(1) = R0*(DeltaHatD*Nv- DeltaHatQ*Ns);
+
+ // plastic energy balance
+ res(2) = ((1-f0)*R*DeltaHatP - kcorEq*DeltaHatD - pcor*DeltaHatQ)/R0;
+
+
+ double DkcorEqDDeltaHatD = -3.*muT;
+ double DpcorDDeltaHatQ = -KT;
+
+ static fullVector<double> gradf, gradNs, gradNv;
+ computeYieldDerivatives(gradf,kcorEq,pcor,R,yieldfV,T);
+ plasticFlowDerivatives(gradNs,gradNv,kcorEq,pcor,R,yieldfV,T);
+
+ double Dres0DfV = gradf(3)*DyieldfVDfV;
+ double Dres1DfV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DyieldfVDfV;
+ double Dres2DfV = 0.;
+
+ // Dres0DDeltaHatD
+ J(0,0) = gradf(0)*DkcorEqDDeltaHatD + Dres0DfV*DDeltafVDDeltaHatD;
+ // Dres0DDeltaHatQ
+ J(0,1) = gradf(1)*DpcorDDeltaHatQ + Dres0DfV*DDeltafVDDeltaHatQ;
+ // Dres0DDeltaHatP
+ J(0,2) = gradf(2)*H + Dres0DfV*DDeltafVDDeltaHatP;
+
+ // Dres1DDeltaHatD
+ J(1,0) = R0*(Nv + DeltaHatD*gradNv(0)*DkcorEqDDeltaHatD - DeltaHatQ*gradNs(0)*DkcorEqDDeltaHatD) + Dres1DfV*DDeltafVDDeltaHatD;
+ // Dres1DDeltaHatQ
+ J(1,1) = R0*(DeltaHatD*gradNv(1)*DpcorDDeltaHatQ- Ns - DeltaHatQ*gradNs(1)*DpcorDDeltaHatQ) + Dres1DfV*DDeltafVDDeltaHatQ;
+ // Dres1DDeltaHatP
+ J(1,2) = R0*(DeltaHatD*gradNv(2) - DeltaHatQ*gradNs(2))*H + Dres1DfV*DDeltafVDDeltaHatP;
+
+ // Dres2DDeltaHatd
+ J(2,0) = (-kcorEq - DkcorEqDDeltaHatD*DeltaHatD)/R0 + Dres2DfV*DDeltafVDDeltaHatD;
+ // Dres2DDeltaHatq
+ J(2,1) = (-pcor - DpcorDDeltaHatQ*DeltaHatQ)/R0 + Dres2DfV*DDeltafVDDeltaHatQ;
+ // Dres2DDeltaHatp
+ J(2,2) = (1-f0)*(R + H*DeltaHatP)/R0+ Dres2DfV*DDeltafVDDeltaHatP;
+
+};
+
void mlawNonLocalDamageGurson::computeDResidualLocal(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
const double yieldfV, const double DyieldfVDfV,
@@ -591,17 +668,6 @@ double mlawNonLocalDamageGurson::getOnsetCriterion(const IPNonLocalPorosity* q1)
};
-void mlawNonLocalDamageGurson::nucleation(const double p1, const double p0, const double f0, IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const{
- // Update nucleation state
- for (int i = 0; i < _gdnLawContainer.size(); i++){
- // Get IpNucl
- const IPNucleation& ipvnucprev = q0->getConstRefToIPNucleation(i);
- IPNucleation& ipvnuc = q1->getRefToIPNucleation(i);
- // Compute...
- _gdnLawContainer[i]->nucleate(p1,p0,f0,ipvnucprev,ipvnuc);
- }
-};
-
void mlawNonLocalDamageGurson::checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double *T) const{
// check coalescence
_coalescenceLaw->checkCoalescence(q0,q1,&q0->getConstRefToIPCoalescence(),&q1->getRefToIPCoalescence());
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h
index 6a8f0cf8d92aedb1b98a1133a666a12c694a5b4a..cfe467bfb5cab2e029541afcff983e25449f3611 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h
+++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h
@@ -83,28 +83,35 @@ class mlawNonLocalDamageGurson : public mlawNonLocalPorosity
virtual double yieldFunction(const double kcorEq, const double pcor, const double R, const double yieldfV,
const IPNonLocalPorosity* q0, const IPNonLocalPorosity* q1, const double* T = NULL) const;
- virtual void computeResidual(SVector3& res, STensor3& J,
+ virtual void computeResidual_NonLocalPorosity(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
+ const double yieldfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
const double* T = NULL
) const;
- virtual void computeDResidual(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
- const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const double DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
- const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
- double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
+ virtual void computeDResidual_NonLocalPorosity(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
+ double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
+ const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
+ const double yieldfV, const STensor3& DyieldfVDKcorpr, const double DyieldfVDtildefV,
+ const STensor43& DKcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
+ const double* T = NULL, const double* DRDT= NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const STensor3* DKcorprDT = NULL,
+ double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
) const;
-
+
+ virtual void computeResidualLocal(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV, const double DyieldfVDfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
+ const double* T = NULL
+ ) const;
virtual void computeDResidualLocal(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
const double yieldfV, const double DyieldfVDfV,
@@ -117,31 +124,27 @@ class mlawNonLocalDamageGurson : public mlawNonLocalPorosity
) const;
virtual void computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J,
- const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double* T = NULL
- ) const{
- computeResidual(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,0.,0.,0.,T);
- };
-
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double* T = NULL
+ ) const;
virtual void computeDResidual_multipleNonLocalVariables(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- std::vector<double> &dres0dtildefV, std::vector<double> &dres1dtildefV, std::vector<double> &dres2dtildefV,
+ std::vector<double> &Dres0DNonlocalVars, std::vector<double> &Dres1DNonlocalVars, std::vector<double> &Dres2DNonlocalVars,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const std::vector<double>& DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDKcorpr, const std::vector<double>& DyieldfVDNonlocalVars,
+ const STensor43& DkcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const STensor3* DKcorprVDT = NULL,
double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
) const;
virtual double getOnsetCriterion(const IPNonLocalPorosity* q1) const;
- virtual void nucleation(const double p1, const double p0, const double f0, IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const;
virtual void checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const;
virtual void checkFailed(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const;
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp b/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp
index 33fecc25f030f0de6257231dfd5e3bfece5cfa11..ff24f3b544e862c33a993237917c09fec67b0de4 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp
@@ -81,11 +81,9 @@ mlawNonLocalPorosity::mlawNonLocalPorosity(const int num,const double E,const do
_triFunc_kw = new constantScalarFunction(1.);
//---------------------------end
if (fVinitial > 0){
- _localRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(std::min(1.e-16,fVinitial), 0., 0.9*_localfVFailure,0.999*_localfVFailure);
_correctedRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(std::min(1.e-16,fVinitial), 0., 0.9*_fV_failure,0.999*_fV_failure);
}
else{
- _localRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(0.,-1.e-16,0.9*_localfVFailure,0.999*_localfVFailure);
_correctedRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(0.,-1.e-16, 0.9*_fV_failure,0.999*_fV_failure);
}
//
@@ -172,10 +170,6 @@ mlawNonLocalPorosity::mlawNonLocalPorosity(const mlawNonLocalPorosity &source) :
}
//-----------------------------end
- _localRegularizedFunction = NULL;
- if (source._localRegularizedFunction != NULL){
- _localRegularizedFunction = source._localRegularizedFunction->clone();
- }
_correctedRegularizedFunction = NULL;
if (source._correctedRegularizedFunction != NULL){
@@ -218,8 +212,6 @@ mlawNonLocalPorosity::~mlawNonLocalPorosity()
if (_triFunc_kw != NULL) delete _triFunc_kw; _triFunc_kw = NULL;
//-----------------------end
- if (_localRegularizedFunction != NULL) delete _localRegularizedFunction;
- _localRegularizedFunction = NULL;
if (_correctedRegularizedFunction != NULL) delete _correctedRegularizedFunction;
_correctedRegularizedFunction = NULL;
if (_temFunc_ThermalConductivity != NULL) delete _temFunc_ThermalConductivity; _temFunc_ThermalConductivity = NULL;
@@ -312,10 +304,6 @@ void mlawNonLocalPorosity::setFailureTolerance(const double NewFailureTol){
_failedTol = NewFailureTol;
};
-void mlawNonLocalPorosity::setLocalRegularizedFunction(const scalarFunction& fct){
- if (_localRegularizedFunction) delete _localRegularizedFunction;
- _localRegularizedFunction = fct.clone();
-};
void mlawNonLocalPorosity::setCorrectedRegularizedFunction(const scalarFunction& fct){
if (_correctedRegularizedFunction) delete _correctedRegularizedFunction;
@@ -586,12 +574,11 @@ void mlawNonLocalPorosity::elasticPredictor( const STensor3& F1, const STensor3&
void mlawNonLocalPorosity::updateInternalVariables(IPNonLocalPorosity *q1, const IPNonLocalPorosity* q0,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, const double DeltafV,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T) const{
double hatPn = q0->getMatrixPlasticStrain();
double hatQn = q0->getVolumetricPlasticStrain();
double hatDn = q0->getDeviatoricPlasticStrain();
- double fVn = q0->getLocalPorosity();
double& hatP = q1->getRefToMatrixPlasticStrain();
hatP = hatPn + DeltaHatP;
@@ -602,9 +589,6 @@ void mlawNonLocalPorosity::updateInternalVariables(IPNonLocalPorosity *q1, const
double& hatD = q1->getRefToDeviatoricPlasticStrain();
hatD = hatDn + DeltaHatD;
- double& fV = q1->getRefToLocalPorosity();
- fV = fVn + DeltafV;
-
// update hardening
if (isThermomechanicallyCoupled()){
_j2IH->hardening(hatPn,q0->getConstRefToIPJ2IsotropicHardening(),hatP,*T, q1->getRefToIPJ2IsotropicHardening());
@@ -618,7 +602,7 @@ void mlawNonLocalPorosity::updateInternalVariables(IPNonLocalPorosity *q1, const
void mlawNonLocalPorosity::updatePlasticState(IPNonLocalPorosity *q1, const IPNonLocalPorosity* q0, const STensor3& F1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, const double DeltafV,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double yield,
STensor3& Fe1, STensor3& Ce, STensor3& Ee, STensor43& Le, STensor63& dLe,
STensor3& DGNp, STensor3& devDGNp, double& trDGNp, STensor43& Dexp, const double* T) const{
@@ -773,20 +757,23 @@ void mlawNonLocalPorosity::predictorCorrector_NonLocalPorosity(
_cLLaw[0]->computeCL(ipvprev->getNonLocalPorosity(), ipvcur->getRefToIPCLength(0));
// Hardening law
- if (isThermomechanicallyCoupled()){
+ if (isThermomechanicallyCoupled())
+ {
_j2IH->hardening(ipvprev->getMatrixPlasticStrain(),ipvprev->getConstRefToIPJ2IsotropicHardening(), ipvcur->getMatrixPlasticStrain(),T, ipvcur->getRefToIPJ2IsotropicHardening());
}
- else{
+ else
+ {
_j2IH->hardening(ipvprev->getMatrixPlasticStrain(),ipvprev->getConstRefToIPJ2IsotropicHardening(), ipvcur->getMatrixPlasticStrain(),ipvcur->getRefToIPJ2IsotropicHardening());
}
// Yield stress
double& R = ipvcur->getRefToCurrentViscoplasticYieldStress();
R = ipvcur->getConstRefToIPJ2IsotropicHardening().getR();
- // nucleation from previous
- // cannot copy because nucleation void part must be set to zero initially
- nucleation(ipvcur->getMatrixPlasticStrain(),ipvprev->getMatrixPlasticStrain(),ipvprev->getLocalPorosity(),ipvcur,ipvprev);
-
+ // no nucleation
+ for (int i=0; i< _gdnLawContainer.size(); i++)
+ {
+ ipvcur->getRefToIPNucleation(i).resetNucleation(ipvprev->getConstRefToIPNucleation(i));
+ }
// copy void state
ipvcur->getRefToIPVoidState() = ipvprev->getConstRefToIPVoidState();
// coalescene from previous
@@ -839,11 +826,13 @@ void mlawNonLocalPorosity::predictorCorrector_NonLocalPorosity(
ipvcur->getRefToDPlasticEnergyDMatrixPlasticStrain() = 0.;
STensorOperation::zero(ipvcur->getRefToDPlasticEnergyDF());
- if (isThermomechanicallyCoupled()){
+ if (isThermomechanicallyCoupled())
+ {
(*dLocalPorosityDT)=0.;
STensorOperation::zero(ipvcur->getRefToDPlasticEnergyDT());
}
- for (int i=0; i< _numNonLocalVar; i++){
+ for (int i=0; i< _numNonLocalVar; i++)
+ {
STensorOperation::zero(ipvcur->getRefToDPlasticEnergyDNonLocalVariable(i));
}
}
@@ -859,8 +848,8 @@ void mlawNonLocalPorosity::predictorCorrector_NonLocalPorosity(
if (plastic)
{
- double DfVtildeStarDtildefVPrev = ipvprev->getConstRefToIPCoalescence().getAccelerateRate();
- correctorOK = plasticCorrectorLocal(F1,ipvprev->getCorrectedPorosity(),DfVtildeStarDtildefVPrev,corKirpr,kcorprDev,ppr,ipvprev,ipvcur,
+ double DfVtildeStarDtildefV = ipvcur->getConstRefToIPCoalescence().getAccelerateRate();
+ correctorOK = plasticCorrectorLocal(F1,ipvprev->getCorrectedPorosity(),DfVtildeStarDtildefV,corKirpr,kcorprDev,ppr,ipvprev,ipvcur,
Fe,Ce,Ee,Le,dLe,
stiff,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,
DcorKirDEpr, DFpDEpr,
@@ -917,7 +906,7 @@ void mlawNonLocalPorosity::predictorCorrector_NonLocalPorosity(
if (plastic){
//Msg::Info("plastic occurs");
- correctorOK = plasticCorrectorLocal(Fcur,ipvTemp.getCorrectedPorosity(),DfVtildeStarDtildefVPrev,corKirpr,kcorprDev,ppr,&ipvTemp,ipvcur,
+ correctorOK = plasticCorrectorLocal(Fcur,ipvTemp.getCorrectedPorosity(),DfVtildeStarDtildefV,corKirpr,kcorprDev,ppr,&ipvTemp,ipvcur,
Fe,Ce,Ee,Le,dLe,
estimateStiff,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,
DcorKirDEpr, DFpDEpr,
@@ -1245,10 +1234,34 @@ void mlawNonLocalPorosity::predictorCorrector_NonLocalPorosity(
+ void mlawNonLocalPorosity::voidCharacteristicsEvoltion_NonLocalPorosity(const STensor3& kcor,
+ const IPNonLocalPorosity *q0, IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
+ double& DfVDDeltaHatD, double& DfVDDeltaHatQ, double& DfVDDeltaHatP,
+ STensor3& DfVDkcorpr) const
+{
+ double hatPn = q0->getMatrixPlasticStrain(); // previous converge value
+ double fn = q0->getLocalPorosity();
+
+ // Update nucleation state
+ nucleation(hatPn+DeltaHatP,hatPn,fn,q1,q0);
+
+ // void evolution
+ localPorosityGrowth(q1, DeltaHatD,DeltaHatQ,DeltaHatP,
+ kcor, q0,
+ true, &DfVDDeltaHatD,&DfVDDeltaHatQ,&DfVDDeltaHatP,&DfVDkcorpr);
+
+ // update void state data
+ _voidEvolutionLaw->evolve(q0->getYieldPorosity(),q1->getYieldPorosity(),DeltaHatD,DeltaHatQ,DeltaHatP,
+ q0,q1, q0->getConstRefToIPVoidState(), q1->getRefToIPVoidState());
+ // update coalescene data
+ _coalescenceLaw->computeIPCoalescence(q0,q1,&q0->getConstRefToIPCoalescence(),&q1->getRefToIPCoalescence());
+};
-bool mlawNonLocalPorosity::plasticCorrector_NonLocalPorosity(const STensor3& F1, const double& tildefVstar,const double& DtildefVstarDtildefV,
+bool mlawNonLocalPorosity::plasticCorrector_NonLocalPorosity(const STensor3& F1,
+ const double& yieldfV,const double& DyieldfVDtildefV,
const STensor3& Kcorpr, const STensor3& devKcorpr, const double& pcorpr,
const IPNonLocalPorosity *q0, IPNonLocalPorosity *q1,
STensor3& Fe, STensor3& Ce, STensor3& Ee, STensor43& Le, STensor63& dLe,
@@ -1261,27 +1274,25 @@ bool mlawNonLocalPorosity::plasticCorrector_NonLocalPorosity(const STensor3& F1,
const STensor3* DKcorprDT, const STensor3* DdevKcorprDT, const double* DpcorprDT,
STensor3* DKcorDT, STensor3* DFpDT, double* DfVDT
) const{
-
- // previous converged state
- double hatPn = q0->getMatrixPlasticStrain();
- double hatQn = q0->getVolumetricPlasticStrain();
- double hatDn = q0->getDeviatoricPlasticStrain();
- double fVn = q0->getLocalPorosity();
-
+ //
+ // yieldfV depends only on nonlocal porosity
+ //
+ double hatPn = q0->getMatrixPlasticStrain(); // previous converge value
+ double fn = q0->getLocalPorosity();
+ // unknowns
double DeltaHatD = 1.e-8;
double DeltaHatQ = 0.;
double DeltaHatP = 1.e-8;
- double DeltafV = 0.;
- double DDeltafVDDeltaHatD = 0.;
- double DDeltafVDDeltaHatQ = 0.;
- double DDeltafVDDeltaHatP = 0.;
- static STensor3 DDeltafVDkcorpr;
-
- // to get initial derivatives
- double kcorEqpr = sqrt(1.5*devKcorpr.dotprod());
- localPorosityGrowth(DeltafV, DeltaHatD,DeltaHatQ,DeltaHatP,
- Kcorpr,q0,q1, true,&DDeltafVDDeltaHatD,&DDeltafVDDeltaHatQ,&DDeltafVDDeltaHatP,&DDeltafVDkcorpr);
+ // dependence of local porosity on unknowns
+ double DfVDDeltaHatD = 0.;
+ double DfVDDeltaHatQ = 0.;
+ double DfVDDeltaHatP = 0.;
+ // and stress tensor due to kw
+ static STensor3 DfVDkcorpr;
+
+ voidCharacteristicsEvoltion_NonLocalPorosity(Kcorpr,q0,q1,DeltaHatD,DeltaHatQ,DeltaHatP,
+ DfVDDeltaHatD,DfVDDeltaHatQ,DfVDDeltaHatP,DfVDkcorpr);
// First residual
static SVector3 res;
@@ -1291,10 +1302,10 @@ bool mlawNonLocalPorosity::plasticCorrector_NonLocalPorosity(const STensor3& F1,
double & yield = q1->getRefToCurrentViscoplasticYieldStress();
yield = q1->getConstRefToIPJ2IsotropicHardening().getR();
double H = q1->getConstRefToIPJ2IsotropicHardening().getDR();
-
- computeResidual(res,J,kcorEqpr,pcorpr,yield,H,tildefVstar,0.,q0,q1,
- DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,
- DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
+
+ double kcorEqpr = sqrt(1.5*devKcorpr.dotprod());
+ computeResidual_NonLocalPorosity(res,J,kcorEqpr,pcorpr,yield,H,yieldfV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,T);
double muT = _mu;
if (isThermomechanicallyCoupled()){
@@ -1339,24 +1350,10 @@ bool mlawNonLocalPorosity::plasticCorrector_NonLocalPorosity(const STensor3& F1,
else
DeltaHatP -= sol(2);
- // Update porosity
- // Update nucleation state
- nucleation(hatPn+DeltaHatP,hatPn,q0->getLocalPorosity(),q1,q0);
-
- // void evolution
- localPorosityGrowth(DeltafV, DeltaHatD,DeltaHatQ,DeltaHatP,
- Kcorpr, q0,q1,
- true, &DDeltafVDDeltaHatD,&DDeltafVDDeltaHatQ,&DDeltafVDDeltaHatP,&DDeltafVDkcorpr);
-
+ voidCharacteristicsEvoltion_NonLocalPorosity(Kcorpr,q0,q1,DeltaHatD,DeltaHatQ,DeltaHatP,
+ DfVDDeltaHatD,DfVDDeltaHatQ,DfVDDeltaHatP,DfVDkcorpr);
// update internal variable
- updateInternalVariables(q1,q0,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
-
- // update void state data
- _voidEvolutionLaw->evolve(q0->getYieldPorosity(),q1->getYieldPorosity(),DeltaHatD,DeltaHatQ,DeltaHatP,
- q0,q1, q0->getConstRefToIPVoidState(), q1->getRefToIPVoidState());
- // update coalescene data
- _coalescenceLaw->computeIPCoalescence(q0,q1,&q0->getConstRefToIPCoalescence(),&q1->getRefToIPCoalescence());
-
+ updateInternalVariables(q1,q0,DeltaHatD,DeltaHatQ,DeltaHatP,T);
yield = q1->getConstRefToIPJ2IsotropicHardening().getR();
H = q1->getConstRefToIPJ2IsotropicHardening().getDR();
@@ -1370,19 +1367,18 @@ bool mlawNonLocalPorosity::plasticCorrector_NonLocalPorosity(const STensor3& F1,
H += dR/delta_t;
}
- computeResidual(res,J,kcorEqpr,pcorpr,yield,H,tildefVstar,0.,q0,q1,
- DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,
- DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
+ computeResidual_NonLocalPorosity(res,J,kcorEqpr,pcorpr,yield,H,yieldfV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,T);
f = res.norm();
- //Msg::Info("ite = %d maxite = %d norm = %e ,DeltaHatD = %e,DeltaHatQ = %e,DeltaHatP = %e,DeltafV =%e !!",ite,_maxite,f,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV);
+ //Msg::Info("ite = %d maxite = %d norm = %e ,DeltaHatD = %e,DeltaHatQ = %e,DeltaHatP = %e !!",ite,_maxite,f,DeltaHatD,DeltaHatQ,DeltaHatP);
if (f < _tol) break;
ite++;
if((ite > _maxite) or STensorOperation::isnan(f))
{
- Msg::Warning("No convergence in NON-local plast. correct. ite= %d: norm= %e, DeltaHatD= %e, DeltaHatQ= %e, DeltaHatP= %e, DeltafV= %e, fv_y = %e, Cft= %e ",
- ite,f,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV, q1->getYieldPorosity(),q1->getConstRefToIPCoalescence().getConcentrationFactor());
+ Msg::Warning("No convergence in NON-local plast. correct. ite= %d: norm= %e, DeltaHatD= %e, DeltaHatQ= %e, DeltaHatP= %e, fv_y = %e, Cft= %e ",
+ ite,f,DeltaHatD,DeltaHatQ,DeltaHatP, q1->getYieldPorosity(),q1->getConstRefToIPCoalescence().getConcentrationFactor());
return false;
}
}
@@ -1392,37 +1388,33 @@ bool mlawNonLocalPorosity::plasticCorrector_NonLocalPorosity(const STensor3& F1,
double trDGNp;
static STensor43 ENp;
- updatePlasticState(q1,q0,F1,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,devKcorpr,kcorEqpr,pcorpr,
+ updatePlasticState(q1,q0,F1,DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,
yield,Fe,Ce,Ee,Le,dLe,DGNp,devDGNp,trDGNp,ENp,T);
double f0=q0->getInitialPorosity();
q1->getRefToDPlasticEnergyDMatrixPlasticStrain() =(1.-f0)*H*DeltaHatP+(1.-f0)*yield;
/* Stiffness computation*/
if(stiff)
{
- static STensor3 DDeltafVDEpr; // partial derivatives
- STensorOperation::multSTensor3STensor43(DDeltafVDkcorpr,DKcorprDEpr,DDeltafVDEpr);
double dRdT= 0.;
- double DDeltafVDT =0.; // partial derivatives
-
if (isThermomechanicallyCoupled()){
dRdT = q1->getConstRefToIPJ2IsotropicHardening().getDRDT();
- DDeltafVDT = STensorOperation::doubledot(DDeltafVDkcorpr,*DKcorprDT);
}
-
/* Compute internal variables derivatives from residual derivatives */
static STensor3 dres0dEpr, dres1dEpr, dres2dEpr;
double dres0dtildefV, dres1dtildefV, dres2dtildefV;
double dres0dT, dres1dT,dres2dT;
+
+ static STensor3 DyieldfVDKcorpr(0.); // fV
- computeDResidual(dres0dEpr,dres1dEpr,dres2dEpr,
+ computeDResidual_NonLocalPorosity(dres0dEpr,dres1dEpr,dres2dEpr,
dres0dtildefV,dres1dtildefV,dres2dtildefV,
devKcorpr,kcorEqpr,pcorpr,yield,
- tildefVstar,0.,DtildefVstarDtildefV,
- DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
+ yieldfV,DyieldfVDKcorpr,DyieldfVDtildefV,
+ DKcorprDEpr,DdevKcorprDEpr,DpcorprDEpr,
q0,q1,
DeltaHatD,DeltaHatQ,DeltaHatP,
- T,&dRdT,DdevKcorprDT,DpcorprDT,&DDeltafVDT,
+ T,&dRdT,DdevKcorprDT,DpcorprDT,DKcorprDT,
&dres0dT,&dres1dT,&dres2dT);
static STensor3 dDeltaHatDdEpr, dDeltaHatQdEpr, dDeltaHatPdEpr;
@@ -1476,20 +1468,20 @@ bool mlawNonLocalPorosity::plasticCorrector_NonLocalPorosity(const STensor3& F1,
dDeltaHatDdT = sol(0);
dDeltaHatQdT = sol(1);
dDeltaHatPdT = sol(2);
- *DfVDT = DDeltafVDT+ DDeltafVDDeltaHatD*dDeltaHatDdT+DDeltafVDDeltaHatQ*dDeltaHatQdT+DDeltafVDDeltaHatP*dDeltaHatPdT;
+ *DfVDT = STensorOperation::doubledot(DfVDkcorpr,*DKcorprDT)+ DfVDDeltaHatD*dDeltaHatDdT+DfVDDeltaHatQ*dDeltaHatQdT+DfVDDeltaHatP*dDeltaHatPdT;
double &DplEnergyDT=q1->getRefToDPlasticEnergyDT();
DplEnergyDT=(1.-f0)*(H*dDeltaHatPdT+dRdT)*DeltaHatP+(1.-f0)*yield*dDeltaHatPdT;
}
// compute DfVDEpr
- DfVDEpr = DDeltafVDEpr;
- DfVDEpr.daxpy(dDeltaHatDdEpr,DDeltafVDDeltaHatD);
- DfVDEpr.daxpy(dDeltaHatQdEpr,DDeltafVDDeltaHatQ);
- DfVDEpr.daxpy(dDeltaHatPdEpr,DDeltafVDDeltaHatP);
+ STensorOperation::multSTensor3STensor43(DfVDkcorpr,DKcorprDEpr,DfVDEpr);
+ DfVDEpr.daxpy(dDeltaHatDdEpr,DfVDDeltaHatD);
+ DfVDEpr.daxpy(dDeltaHatQdEpr,DfVDDeltaHatQ);
+ DfVDEpr.daxpy(dDeltaHatPdEpr,DfVDDeltaHatP);
// compute DfVDtildefV
- DfVDtildefV = DDeltafVDDeltaHatD*dDeltaHatDdtildefV+DDeltafVDDeltaHatQ*dDeltaHatQdtildefV+DDeltafVDDeltaHatP*dDeltaHatPdtildefV;
+ DfVDtildefV = DfVDDeltaHatD*dDeltaHatDdtildefV+DfVDDeltaHatQ*dDeltaHatQdtildefV+DfVDDeltaHatP*dDeltaHatPdtildefV;
computePlasticTangentNonLocal(q1,q0,devKcorpr,kcorEqpr,pcorpr,DGNp,devDGNp,trDGNp,ENp,
DeltaHatD,DeltaHatQ,DeltaHatP,
@@ -1542,9 +1534,11 @@ void mlawNonLocalPorosity::predictorCorrectorLocal(
_j2IH->hardening(ipvprev->getMatrixPlasticStrain(),ipvprev->getConstRefToIPJ2IsotropicHardening(), ipvcur->getMatrixPlasticStrain(),ipvcur->getRefToIPJ2IsotropicHardening());
}
- // nucleation from previous
- // cannot copy because nucleation void part must be set to zero initially
- nucleation(ipvcur->getMatrixPlasticStrain(),ipvprev->getMatrixPlasticStrain(),ipvprev->getLocalPorosity(),ipvcur,ipvprev);
+ // no nucleation
+ for (int i=0; i< _gdnLawContainer.size(); i++)
+ {
+ ipvcur->getRefToIPNucleation(i).resetNucleation(ipvprev->getConstRefToIPNucleation(i));
+ }
// reset to previous state
ipvcur->getRefToIPCoalescence() = (ipvprev->getConstRefToIPCoalescence());
@@ -1764,7 +1758,36 @@ void mlawNonLocalPorosity::predictorCorrectorLocal(
+void mlawNonLocalPorosity::voidCharacteristicsEvoltionLocal(const STensor3& kcor,
+ const IPNonLocalPorosity *q0, IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
+ const double DtildefVstarDDeltafV,
+ double& DyieldfVDDeltafV,
+ double& DDeltafVDDeltaHatD, double& DDeltafVDDeltaHatQ, double& DDeltafVDDeltaHatP,
+ STensor3& DDeltafVDKcor) const
+{
+ double hatPn = q0->getMatrixPlasticStrain();
+ // nucleatino
+ nucleation(hatPn+DeltaHatP,hatPn,q0->getLocalPorosity(),q1,q0);
+ // local void growth
+ localPorosityGrowth(q1, DeltaHatD,DeltaHatQ,DeltaHatP,
+ kcor,q0,true,&DDeltafVDDeltaHatD,&DDeltafVDDeltaHatQ,&DDeltafVDDeltaHatP,&DDeltafVDKcor);
+
+ // yield porosity
+ double & tildefVstar = q1->getRefToCorrectedPorosity();
+ double& yieldfV = q1->getRefToYieldPorosity();
+
+ tildefVstar = q0->getCorrectedPorosity() + DtildefVstarDDeltafV*(q1->getLocalPorosity() - q0->getLocalPorosity());
+ yieldfV = _correctedRegularizedFunction->getVal(tildefVstar);
+ DyieldfVDDeltafV = _correctedRegularizedFunction->getDiff(tildefVstar)*DtildefVstarDDeltafV;
+
+ // update void
+ _voidEvolutionLaw->evolve(q0->getYieldPorosity(),q1->getYieldPorosity(),DeltaHatD,DeltaHatQ,DeltaHatP,
+ q0,q1, q0->getConstRefToIPVoidState(), q1->getRefToIPVoidState());
+ _coalescenceLaw->computeIPCoalescence(q0,q1,&q0->getConstRefToIPCoalescence(),&q1->getRefToIPCoalescence());
+
+};
@@ -1789,37 +1812,31 @@ bool mlawNonLocalPorosity::plasticCorrectorLocal(const STensor3& F1, const doubl
double DeltaHatD = 1.e-8;
double DeltaHatQ = 0.;
double DeltaHatP = 1.e-8;
- double DeltafV = 0.;
+
double DDeltafVDDeltaHatD = 0.;
double DDeltafVDDeltaHatQ = 0.;
- double DDeltafVDDeltaHatP = 0.;
-
-
+ double DDeltafVDDeltaHatP = 0.;
static STensor3 DDeltafVDKcorpr;
+ double DyieldfVDDeltafV = 1.;
+
+ voidCharacteristicsEvoltionLocal(Kcorpr,q0,q1,DeltaHatD,DeltaHatQ,DeltaHatP,
+ DtildefVstarDDeltafV,
+ DyieldfVDDeltafV, DDeltafVDDeltaHatD, DDeltafVDDeltaHatQ, DDeltafVDDeltaHatP,DDeltafVDKcorpr);
double kcorEqpr = sqrt(1.5*devKcorpr.dotprod());
- // get derivatives
- localPorosityGrowth(DeltafV, DeltaHatD,DeltaHatQ,DeltaHatP,
- Kcorpr,q0,q1,
- true,&DDeltafVDDeltaHatD,&DDeltafVDDeltaHatQ,&DDeltafVDDeltaHatP,&DDeltafVDKcorpr);
-
+
// First residual
static SVector3 res;
static SVector3 sol;
static STensor3 J, invJ;
-
- double & yield = q1->getRefToCurrentViscoplasticYieldStress();
+
+ double & yield = q1->getRefToCurrentViscoplasticYieldStress();
yield = q1->getConstRefToIPJ2IsotropicHardening().getR();
double H = q1->getConstRefToIPJ2IsotropicHardening().getDR();
-
- double & tildefVstar = q1->getRefToCorrectedPorosity();
- tildefVstar = tildefVstarPrev + DtildefVstarDDeltafV*DeltafV;
- double& yieldfV = q1->getRefToYieldPorosity();
- yieldfV = _correctedRegularizedFunction->getVal(tildefVstar);
- double DyieldfVDDeltafV = _correctedRegularizedFunction->getDiff(tildefVstar)*DtildefVstarDDeltafV;
-
- computeResidual(res,J,kcorEqpr,pcorpr,yield,H,yieldfV,DyieldfVDDeltafV,q0,q1,
- DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,
+
+ double yieldfV = q1->getYieldPorosity();
+ computeResidualLocal(res,J,kcorEqpr,pcorpr,yield,H,yieldfV,DyieldfVDDeltafV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,
DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
double muT = _mu;
@@ -1863,24 +1880,11 @@ bool mlawNonLocalPorosity::plasticCorrectorLocal(const STensor3& F1, const doubl
else
DeltaHatP -= sol(2);
+ voidCharacteristicsEvoltionLocal(Kcorpr,q0,q1,DeltaHatD,DeltaHatQ,DeltaHatP,
+ DtildefVstarDDeltafV,
+ DyieldfVDDeltafV, DDeltafVDDeltaHatD, DDeltafVDDeltaHatQ, DDeltafVDDeltaHatP,DDeltafVDKcorpr);
-
- nucleation(hatPn+DeltaHatP,hatPn,q0->getLocalPorosity(),q1,q0);
-
- localPorosityGrowth(DeltafV, DeltaHatD,DeltaHatQ,DeltaHatP,
- Kcorpr,q0,q1,
- true, &DDeltafVDDeltaHatD,&DDeltafVDDeltaHatQ,&DDeltafVDDeltaHatP,&DDeltafVDKcorpr);
-
- updateInternalVariables(q1,q0,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
-
- tildefVstar = tildefVstarPrev + DtildefVstarDDeltafV*DeltafV;
- yieldfV = _correctedRegularizedFunction->getVal(tildefVstar);
- DyieldfVDDeltafV = _correctedRegularizedFunction->getDiff(tildefVstar)*DtildefVstarDDeltafV;
-
- // update void
- _voidEvolutionLaw->evolve(q0->getYieldPorosity(),q1->getYieldPorosity(),DeltaHatD,DeltaHatQ,DeltaHatP,
- q0,q1, q0->getConstRefToIPVoidState(), q1->getRefToIPVoidState());
- _coalescenceLaw->computeIPCoalescence(q0,q1,&q0->getConstRefToIPCoalescence(),&q1->getRefToIPCoalescence());
+ updateInternalVariables(q1,q0,DeltaHatD,DeltaHatQ,DeltaHatP,T);
yield = q1->getConstRefToIPJ2IsotropicHardening().getR();
H = q1->getConstRefToIPJ2IsotropicHardening().getDR();
@@ -1893,9 +1897,10 @@ bool mlawNonLocalPorosity::plasticCorrectorLocal(const STensor3& F1, const doubl
yield += R;
H += dR/delta_t;
}
-
- computeResidual(res,J,kcorEqpr,pcorpr,yield,H,yieldfV,DyieldfVDDeltafV,q0,q1,
- DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,
+
+ yieldfV = q1->getYieldPorosity();
+ computeResidualLocal(res,J,kcorEqpr,pcorpr,yield,H,yieldfV,DyieldfVDDeltafV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,
DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
f = res.norm();
@@ -1906,12 +1911,12 @@ bool mlawNonLocalPorosity::plasticCorrectorLocal(const STensor3& F1, const doubl
if((ite > _maxite) or STensorOperation::isnan(f))
{
if (q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag()){
- Msg::Warning("No convergence in LOCAL plast. correct. ite= %d: active coales, norm= %e, DeltaHatD= %e, DeltaHatQ= %e, DeltaHatP= %e, DeltafV= %e, fv_y= %e,Cft= %e ",
- ite,f,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV, q1->getYieldPorosity(),q1->getConstRefToIPCoalescence().getConcentrationFactor());
+ Msg::Warning("No convergence in LOCAL plast. correct. ite= %d: active coales, norm= %e, DeltaHatD= %e, DeltaHatQ= %e, DeltaHatP= %e, fv_y= %e,Cft= %e ",
+ ite,f,DeltaHatD,DeltaHatQ,DeltaHatP, q1->getYieldPorosity(),q1->getConstRefToIPCoalescence().getConcentrationFactor());
}
else{
- Msg::Warning("No convergence in LOCAL plast. correct. ite= %d: inactive coa., norm= %e, DeltaHatD= %e, DeltaHatQ= %e, DeltaHatP= %e, DeltafV= %e, fv_y= %e,Cft= %e ",
- ite,f,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV, q1->getYieldPorosity(),q1->getConstRefToIPCoalescence().getConcentrationFactor());
+ Msg::Warning("No convergence in LOCAL plast. correct. ite= %d: inactive coa., norm= %e, DeltaHatD= %e, DeltaHatQ= %e, DeltaHatP= %e, fv_y= %e,Cft= %e ",
+ ite,f,DeltaHatD,DeltaHatQ,DeltaHatP, q1->getYieldPorosity(),q1->getConstRefToIPCoalescence().getConcentrationFactor());
double J3 = STensorOperation::determinantSTensor3(devKcorpr);
Msg::Warning("res=[%e; %e; %e;] fVloc= %e, tauEqpr= %e , ppr = %e, J3= %e, DyieldfVDDeltafV= %e, DDeltafVDDeltaHatV= [%e %e %e], yieldfG = %e",
res(0),res(1),res(2),q0->getLocalPorosity(), kcorEqpr/_j2IH->getYield0(), pcorpr/_j2IH->getYield0(), 27.*J3/(2.*kcorEqpr*kcorEqpr*kcorEqpr),
@@ -1925,7 +1930,7 @@ bool mlawNonLocalPorosity::plasticCorrectorLocal(const STensor3& F1, const doubl
static STensor3 DGNp, devDGNp;
double trDGNp;
static STensor43 ENp;
- updatePlasticState(q1,q0,F1,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,devKcorpr,kcorEqpr,pcorpr,
+ updatePlasticState(q1,q0,F1,DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,
yield,Fe,Ce,Ee,Le,dLe,DGNp,devDGNp,trDGNp,ENp,T);
double f0=q0->getInitialPorosity();
q1->getRefToDPlasticEnergyDMatrixPlasticStrain() = (1.-f0)*H*DeltaHatP + (1.-f0)*yield;
@@ -3543,18 +3548,36 @@ void mlawNonLocalPorosity::constitutive_NonLocalLogarithmPorosity(
}
};
-void mlawNonLocalPorosity::localPorosityGrowth(double& DeltafV,
+void mlawNonLocalPorosity::nucleation(const double p1, const double p0, const double f0, IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const{
+ // Update nucleation state
+ for (int i = 0; i < _gdnLawContainer.size(); i++){
+ // Get IpNucl
+ const IPNucleation& ipvnucprev = q0->getConstRefToIPNucleation(i);
+ IPNucleation& ipvnuc = q1->getRefToIPNucleation(i);
+ // Compute...
+ if (q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag())
+ {
+ ipvnuc.resetNucleation(ipvnucprev);
+ }
+ else
+ {
+ _gdnLawContainer[i]->nucleate(p1,p0,f0,ipvnucprev,ipvnuc);
+
+ }
+ }
+};
+
+void mlawNonLocalPorosity::localPorosityGrowth(IPNonLocalPorosity *q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const STensor3& Kcor,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity *q1,
- bool stiff,
- double* DDeltafVDDeltaHatD, double* DDeltafVDDeltaHatQ, double* DDeltafVDDeltaHatP, STensor3* DDeltafVDKcor) const
+ const IPNonLocalPorosity *q0,
+ bool stiff, double* DDeltafVDDeltaHatD, double* DDeltafVDDeltaHatQ, double* DDeltafVDDeltaHatP, STensor3* DDeltafVDKcor) const
{
// Compute local porosity growth if damage is not blocked (otherwise : no more growth)
if (q1->dissipationIsBlocked())
{
// No more growth as dissipation is blocked
- DeltafV =0.;
+ q1->getRefToLocalPorosity() = q0->getLocalPorosity();
if (stiff){
*DDeltafVDDeltaHatD = 0.;
*DDeltafVDDeltaHatQ = 0.;
@@ -3606,20 +3629,23 @@ void mlawNonLocalPorosity::localPorosityGrowth(double& DeltafV,
if (kw > 0. and DeltaHatD > 0. and kcorEq > 0)
{
-
J3 = STensorOperation::determinantSTensor3(devKcor);
-
Xi = 27.*J3/(2.*kcorEq*kcorEq*kcorEq);
u += kw*(1.- Xi*Xi)*fVn*DeltaHatD;
v -= kw*(1.- Xi*Xi)*DeltaHatD;
}
// get local porosity
- DeltafV = u/v;
+ double DeltafV = u/v;
// Regularise porosity growth and rates
- double ffPrev = _localRegularizedFunction->inverse(fVn);
- double ff = ffPrev+DeltafV;
- DeltafV = _localRegularizedFunction->getVal(ff) - fVn;
+ double regularizedRate = 1.;
+ if (fVn + DeltafV > _localfVFailure){
+ q1->getRefToLocalPorosity() = _localfVFailure;
+ regularizedRate = (_localfVFailure - q1->getLocalPorosity())/DeltafV;
+ }
+ else{
+ q1->getRefToLocalPorosity() = q0->getLocalPorosity()+ DeltafV;
+ }
if (stiff){
// Compute derivatives
@@ -3667,7 +3693,6 @@ void mlawNonLocalPorosity::localPorosityGrowth(double& DeltafV,
STensorOperation::multSTensor3STensor43(DDeltafVDdevKcor,_I4dev,*DDeltafVDKcor);
}
- double regularizedRate = _localRegularizedFunction->getDiff(ff);
*DDeltafVDDeltaHatD *= regularizedRate;
*DDeltafVDDeltaHatQ *= regularizedRate;
*DDeltafVDDeltaHatP *= regularizedRate;
@@ -3936,7 +3961,8 @@ void mlawNonLocalPorosity::constitutive(
}
};
-bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STensor3& F1, const double& yieldfV,const double& DyieldfVDDeltafV,
+bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STensor3& F1,
+ const double& yieldfV,const std::vector<double>& DyieldfVDnonlocalVar, const STensor3& DyieldfVDKcorpr,
const STensor3& Kcorpr, const STensor3& devKcorpr, const double& pcorpr,
const IPNonLocalPorosity *q0, IPNonLocalPorosity *q1,
STensor3& Fe, STensor3& Ce, STensor3& Ee, STensor43& Le, STensor63& dLe,
@@ -3949,28 +3975,12 @@ bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STen
const STensor3* DKcorprDT, const STensor3* DdevKcorprDT, const double* DpcorprDT,
STensor3* DKcorDT, STensor3* DFpDT, fullMatrix<double>* DlocalVarDT
) const{
-
+
+ // unknowns
double DeltaHatD = 0.;
double DeltaHatQ = 0.;
double DeltaHatP = 0.;
- // deformation
- double DeltafV = 0.;
-
- double DeltaHatDNonLocal = q1->getNonLocalDeviatoricPlasticStrain() - q0->getNonLocalDeviatoricPlasticStrain();
- double DeltaHatQNonLocal = q1->getNonLocalVolumetricPlasticStrain() - q0->getNonLocalVolumetricPlasticStrain();
- double DeltaHatPNonLocal = q1->getNonLocalMatrixPlasticStrain() - q0->getNonLocalMatrixPlasticStrain();
-
- double DDeltafVDDeltaHatDNonLocal = 0.;
- double DDeltafVDDeltaHatQNonLocal = 0.;
- double DDeltafVDDeltaHatPNonLocal = 0.;
-
- static STensor3 DDeltafVDKcorpr;
-
- double kcorEqpr = sqrt(1.5*devKcorpr.dotprod());
- localPorosityGrowth(DeltafV, DeltaHatDNonLocal,DeltaHatQNonLocal,DeltaHatPNonLocal,
- Kcorpr,q0,q1,
- true, &DDeltafVDDeltaHatDNonLocal,&DDeltafVDDeltaHatQNonLocal,&DDeltafVDDeltaHatPNonLocal,&DDeltafVDKcorpr);
// First residual
static SVector3 res;
@@ -3981,9 +3991,10 @@ bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STen
double & yield = q1->getRefToCurrentViscoplasticYieldStress();
yield = q1->getConstRefToIPJ2IsotropicHardening().getR();
double H = q1->getConstRefToIPJ2IsotropicHardening().getDR();
-
- computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,yield,H,yieldfV,DyieldfVDDeltafV,q0,q1,
- DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
+
+ double kcorEqpr = sqrt(1.5*STensorOperation::doubledot(devKcorpr,devKcorpr));
+ computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,yield,H,yieldfV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,T);
int ite = 0;
double f = res.norm();
@@ -4017,7 +4028,7 @@ bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STen
else
DeltaHatP -= sol(2);
- updateInternalVariables(q1,q0,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
+ updateInternalVariables(q1,q0,DeltaHatD,DeltaHatQ,DeltaHatP,T);
yield = q1->getConstRefToIPJ2IsotropicHardening().getR();
H = q1->getConstRefToIPJ2IsotropicHardening().getDR();
// Add strain rate effects
@@ -4030,8 +4041,8 @@ bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STen
H += dR/delta_t;
}
- computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,yield,H,yieldfV,DyieldfVDDeltafV,q0,q1,
- DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
+ computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,yield,H,yieldfV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,T);
f = res.norm();
//printf("ite = %d maxite = %d norm = %e ,DeltaHatD = %e,DeltaHatQ = %e,DeltaHatP = %e,DeltafV =%e !! \n",ite,_maxite,f,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV);
@@ -4040,7 +4051,7 @@ bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STen
ite++;
if((ite > _maxite) or STensorOperation::isnan(f))
{
- printf("No convergence for plastic correction with ite = %d maxite = %d norm = %e ,DeltaHatD = %e,DeltaHatQ = %e,DeltaHatP = %e,DeltafV =%e !!\n",ite,_maxite,f,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV);
+ printf("No convergence for plastic correction with ite = %d maxite = %d norm = %e ,DeltaHatD = %e,DeltaHatQ = %e,DeltaHatP = %e!!\n",ite,_maxite,f,DeltaHatD,DeltaHatQ,DeltaHatP);
return false;
}
}
@@ -4050,9 +4061,9 @@ bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STen
static STensor43 ENp;
// update internal value first
- updateInternalVariables(q1,q0,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
+ updateInternalVariables(q1,q0,DeltaHatD,DeltaHatQ,DeltaHatP,T);
- updatePlasticState(q1,q0,F1,DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,devKcorpr,kcorEqpr,pcorpr,
+ updatePlasticState(q1,q0,F1,DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,
yield,Fe,Ce,Ee,Le,dLe,DGNp,devDGNp,trDGNp,ENp,T);
double DplasticEnergyDHatP = (1-q0->getInitialPorosity())*(yield +H*DeltaHatP);
@@ -4061,39 +4072,24 @@ bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STen
/* Stiffness computation*/
if(stiff)
{
- static STensor3 DDeltafVDEpr; // partial derivatives
- STensorOperation::multSTensor3STensor43(DDeltafVDKcorpr,DKcorprDEpr,DDeltafVDEpr);
- double DDeltafVDT =0.; // partial derivatives
double dRdT= 0.;
if (isThermomechanicallyCoupled()){
dRdT = q1->getConstRefToIPJ2IsotropicHardening().getDRDT();
- DDeltafVDT = STensorOperation::doubledot(DDeltafVDKcorpr,*DKcorprDT);
}
-
/* Compute internal variables derivatives from residual derivatives */
static STensor3 dres0dEpr, dres1dEpr, dres2dEpr;
- static std::vector<double> Dres0DNonLocalVar, Dres1DNonLocalVar, Dres2DNonLocalVar, DyieldfVDNonLocalVar;
- if (Dres0DNonLocalVar.size() != _numNonLocalVar){
- Dres0DNonLocalVar.resize(_numNonLocalVar);
- Dres1DNonLocalVar.resize(_numNonLocalVar);
- Dres2DNonLocalVar.resize(_numNonLocalVar);
- DyieldfVDNonLocalVar.resize(_numNonLocalVar);
- }
-
- DyieldfVDNonLocalVar[0] = DyieldfVDDeltafV*DDeltafVDDeltaHatQNonLocal;
- DyieldfVDNonLocalVar[1] = DyieldfVDDeltafV*DDeltafVDDeltaHatPNonLocal;
- DyieldfVDNonLocalVar[2] = DyieldfVDDeltafV*DDeltafVDDeltaHatDNonLocal;
-
+ static std::vector<double> Dres0DNonLocalVar(_numNonLocalVar,0.), Dres1DNonLocalVar(_numNonLocalVar,0.), Dres2DNonLocalVar(_numNonLocalVar,0.);
+
double dres0dT, dres1dT,dres2dT;
computeDResidual_multipleNonLocalVariables(dres0dEpr,dres1dEpr,dres2dEpr,
Dres0DNonLocalVar,Dres1DNonLocalVar,Dres2DNonLocalVar,
devKcorpr,kcorEqpr,pcorpr,yield,
- yieldfV,DyieldfVDDeltafV,DyieldfVDNonLocalVar,
- DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
+ yieldfV,DyieldfVDKcorpr,DyieldfVDnonlocalVar,
+ DKcorprDEpr,DdevKcorprDEpr,DpcorprDEpr,
q0,q1,
DeltaHatD,DeltaHatQ,DeltaHatP,
- T,&dRdT,DdevKcorprDT,DpcorprDT,&DDeltafVDT,
+ T,&dRdT,DdevKcorprDT,DpcorprDT,DKcorprDT,
&dres0dT,&dres1dT,&dres2dT);
static STensor3 dDeltaHatDdEpr, dDeltaHatQdEpr, dDeltaHatPdEpr;
@@ -4118,12 +4114,7 @@ bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STen
DlocalVarDEpr[1] = dDeltaHatPdEpr;
DlocalVarDEpr[2] = dDeltaHatDdEpr;
- static std::vector<double> dDeltaHatDDNonLocalVar, dDeltaHatQDNonLocalVar, dDeltaHatPDNonLocalVar;
- if (dDeltaHatDDNonLocalVar.size() != _numNonLocalVar){
- dDeltaHatDDNonLocalVar.resize(_numNonLocalVar);
- dDeltaHatQDNonLocalVar.resize(_numNonLocalVar);
- dDeltaHatPDNonLocalVar.resize(_numNonLocalVar);
- }
+ static std::vector<double> dDeltaHatDDNonLocalVar(_numNonLocalVar,0.), dDeltaHatQDNonLocalVar(_numNonLocalVar,0.), dDeltaHatPDNonLocalVar(_numNonLocalVar,0.);
for (int i=0; i< _numNonLocalVar; i++){
// dtilde_fv
@@ -4186,6 +4177,70 @@ bool mlawNonLocalPorosity::plasticCorrector_multipleNonLocalVariables(const STen
return true;
};
+void mlawNonLocalPorosity::voidCharacteristicsEvoltion_multipleNonLocalVariables(const STensor3& kcor,
+ const IPNonLocalPorosity *q0, IPNonLocalPorosity* q1,
+ std::vector<double>& DyieldfVDNonlocalVars,
+ STensor3& DyieldfVDKcor, const bool blockVoidGrowth
+ ) const
+{
+ if (blockVoidGrowth)
+ {
+ // purely plastic problem without void growth
+ q1->getRefToLocalPorosity() = q0->getLocalPorosity();
+ q1->getRefToNonLocalPorosity() = q0->getNonLocalPorosity();
+ q1->getRefToCorrectedPorosity() = q0->getCorrectedPorosity();
+ q1->getRefToLocalLogarithmPorosity() = q0->getLocalLogarithmPorosity();
+ q1->getRefToYieldPorosity() = q0->getYieldPorosity();
+ //
+ for (int i=0; i< _numNonLocalVar; i++)
+ {
+ DyieldfVDNonlocalVars[i] = 0.;
+ }
+ STensorOperation::zero(DyieldfVDKcor);
+ }
+ else
+ {
+ // note the computation order
+ // nucleation--> porosities (local, nonlocal, yield ...) -> void state -> coalescence
+
+ // Update nucleation state first
+ nucleation(q1->getNonLocalMatrixPlasticStrain(),q0->getNonLocalMatrixPlasticStrain(),q0->getYieldPorosity(),q1,q0);
+
+ // update internal data with nonlocal variable
+ // nonlocal increment
+ double DeltaHatQNonLocal = q1->getNonLocalVolumetricPlasticStrain() - q0->getNonLocalVolumetricPlasticStrain();
+ double DeltaHatPNonLocal = q1->getNonLocalMatrixPlasticStrain() - q0->getNonLocalMatrixPlasticStrain();
+ double DeltaHatDNonLocal = q1->getNonLocalDeviatoricPlasticStrain() - q0->getNonLocalDeviatoricPlasticStrain();
+
+ // note that nonlocal var = [hatQ hatP hatD]
+ // void growth --> yield porosity = local porosity
+ localPorosityGrowth(q1,DeltaHatDNonLocal,DeltaHatQNonLocal,DeltaHatPNonLocal,
+ kcor,q0,true,
+ &DyieldfVDNonlocalVars[2],&DyieldfVDNonlocalVars[0],&DyieldfVDNonlocalVars[1],&DyieldfVDKcor);
+ // nonlocal
+ q1->getRefToNonLocalPorosity() = q1->getLocalPorosity();
+ // corrected porosity
+ double & tildefVstar = q1->getRefToCorrectedPorosity();
+ double DfVtildeStarDtildefV = q1->getConstRefToIPCoalescence().getAccelerateRate(); // current rate
+ tildefVstar = q0->getCorrectedPorosity() + DfVtildeStarDtildefV*(q1->getNonLocalPorosity() - q0->getNonLocalPorosity());
+ // yield
+ q1->getRefToYieldPorosity() = _correctedRegularizedFunction->getVal(tildefVstar);
+
+ // derivatives of yield porosity
+ double DyieldfVDDtildefV = _correctedRegularizedFunction->getDiff(tildefVstar)*DfVtildeStarDtildefV;
+ for (int i=0; i< 3; i++)
+ {
+ DyieldfVDNonlocalVars[i] *= DyieldfVDDtildefV;
+ }
+ DyieldfVDKcor *= (DyieldfVDDtildefV);
+
+ // update void state
+ _voidEvolutionLaw->evolve(q0->getYieldPorosity(),q1->getYieldPorosity(),DeltaHatDNonLocal,DeltaHatQNonLocal,DeltaHatPNonLocal,
+ q0,q1, q0->getConstRefToIPVoidState(),q1->getRefToIPVoidState());
+ _coalescenceLaw->computeIPCoalescence(q0,q1,&q0->getConstRefToIPCoalescence(),&q1->getRefToIPCoalescence());
+ }
+};
+
void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const STensor3& F0, // initial deformation gradient (input @ time n)
const STensor3& F1, // updated deformation gradient (input @ time n+1)
STensor3 &P, // updated 1st Piola-Kirchhoff stress tensor (output)
@@ -4248,7 +4303,10 @@ void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const ST
ipvcur->getRefToIPCoalescence() = (ipvprev->getConstRefToIPCoalescence());
//
// no nucleation
- this->nucleation(ipvcur->getNonLocalMatrixPlasticStrain(),ipvprev->getNonLocalMatrixPlasticStrain(), ipvprev->getYieldPorosity(),ipvcur,ipvprev);
+ for (int i=0; i< _gdnLawContainer.size(); i++)
+ {
+ ipvcur->getRefToIPNucleation(i).resetNucleation(ipvprev->getConstRefToIPNucleation(i));
+ }
//
const STensor3* Fp0 = &(ipvprev->getConstRefToFp());
@@ -4438,51 +4496,17 @@ void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const ST
}
if (!elasticFollowedBlocked){
- //Non local
- // update internal data with nonlocal variable
- // nonlocal in
- double DeltaHatQNonLocal = ipvcur->getNonLocalVolumetricPlasticStrain() - ipvprev->getNonLocalVolumetricPlasticStrain();
- double DeltaHatDNonLocal = ipvcur->getNonLocalDeviatoricPlasticStrain() - ipvprev->getNonLocalDeviatoricPlasticStrain();
- double DeltaHatPNonLocal = ipvcur->getNonLocalMatrixPlasticStrain() - ipvprev->getNonLocalMatrixPlasticStrain();
- double DeltafV = 0.;
-
- // Update nucleation state
- nucleation(ipvcur->getNonLocalMatrixPlasticStrain(),ipvprev->getNonLocalMatrixPlasticStrain(),ipvprev->getYieldPorosity(),ipvcur,ipvprev);
-
- // void growth
- localPorosityGrowth(DeltafV,DeltaHatDNonLocal,DeltaHatQNonLocal,DeltaHatPNonLocal,
- corKirpr,ipvprev,ipvcur);
-
- // internal data
- updateInternalVariables(ipvcur,ipvprev,0.,0.,0.,DeltafV,&T);
-
- // update yield porosity as it depends on nonlocal variable
- double & tildefVstar = ipvcur->getRefToCorrectedPorosity();
- double DfVtildeStarDtildefVPrev = ipvprev->getConstRefToIPCoalescence().getAccelerateRate();
- tildefVstar = ipvprev->getCorrectedPorosity() + DfVtildeStarDtildefVPrev*DeltafV;
-
- double& yieldfV = ipvcur->getRefToYieldPorosity(); // value
- double DyieldfVDDtildefV; // rate
- if (plasticFollowedBlocked){
- yieldfV = ipvprev->getYieldPorosity();
- DyieldfVDDtildefV = 0.;
- }
- else{
- yieldfV = _correctedRegularizedFunction->getVal(tildefVstar);
- DyieldfVDDtildefV = _correctedRegularizedFunction->getDiff(tildefVstar)*DfVtildeStarDtildefVPrev;
- }
-
- _voidEvolutionLaw->evolve(ipvprev->getYieldPorosity(),ipvcur->getYieldPorosity(),DeltaHatDNonLocal,DeltaHatQNonLocal,DeltaHatPNonLocal,
- ipvprev,ipvcur, ipvprev->getConstRefToIPVoidState(), ipvcur->getRefToIPVoidState());
- _coalescenceLaw->computeIPCoalescence(ipvprev,ipvcur,&ipvprev->getConstRefToIPCoalescence(),&ipvcur->getRefToIPCoalescence());
-
-
+ static std::vector<double> DyieldfVDNonlocalVar(3,0.);
+ static STensor3 DyieldfVDKcorpr;
+ voidCharacteristicsEvoltion_multipleNonLocalVariables(corKirpr,ipvprev,ipvcur,DyieldfVDNonlocalVar,DyieldfVDKcorpr,plasticFollowedBlocked);
// check with previous value
+ double yieldfV = ipvcur->getYieldPorosity(); // does not change during corrector
double yieldF = yieldFunction(kCorPrEq,ppr,R,yieldfV,ipvprev,ipvcur,&T);
plastic = (yieldF > _tol);
-
if (plastic){
- correctorOK = plasticCorrector_multipleNonLocalVariables(F1,yieldfV,DyieldfVDDtildefV,corKirpr,kcorprDev,ppr,ipvprev,ipvcur,
+ correctorOK = plasticCorrector_multipleNonLocalVariables(F1,
+ yieldfV,DyieldfVDNonlocalVar,DyieldfVDKcorpr,
+ corKirpr,kcorprDev,ppr,ipvprev,ipvcur,
Fe,Ce,Ee,Le,dLe,stiff,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,
DcorKirDEpr,DcorKirDNonLocalVar,
DFpDEpr,DFpDNonLocalVar,
@@ -4494,17 +4518,21 @@ void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const ST
// substepping devide by 0 untile converge
int numSubStep = 2.;
int numAttempt = 0;
- static STensor3 dF, Fcur;
+
+ static STensor3 dF;
dF = F1;
dF -= F0;
-
- static IPNonLocalPorosity ipvTemp(*ipvprev);
- double dT,Tcur;//---------added for substepping temperature
- dT = T - T0;
+ double DeltaHatQNonLocal = ipvcur->getNonLocalVolumetricPlasticStrain() - ipvprev->getNonLocalVolumetricPlasticStrain();
+ double DeltaHatPNonLocal = ipvcur->getNonLocalMatrixPlasticStrain() - ipvprev->getNonLocalMatrixPlasticStrain();
+ double DeltaHatDNonLocal = ipvcur->getNonLocalDeviatoricPlasticStrain() - ipvprev->getNonLocalDeviatoricPlasticStrain();
+ double dT = T - T0;
+
+ static IPNonLocalPorosity* ipvTemp = dynamic_cast<IPNonLocalPorosity*>(ipvprev->clone()); // do not known type of ip a priori
+
while (true){
bool success = true;
int iter = 0;
- ipvTemp.operator =(*dynamic_cast<const IPVariable*>(ipvprev)); // take privious value
+ ipvTemp->operator =(*dynamic_cast<const IPVariable*>(ipvprev)); // take privious value
while (iter < numSubStep){
@@ -4516,61 +4544,33 @@ void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const ST
estimateStiff = stiff;
}
- //
+ //update current kinematic variables
+ static STensor3 Fcur;
Fcur = F0;
Fcur.daxpy(dF,fact);
ipvcur->getRefToNonLocalVolumetricPlasticStrain() = ipvprev->getNonLocalVolumetricPlasticStrain() + fact*DeltaHatQNonLocal;
ipvcur->getRefToNonLocalDeviatoricPlasticStrain() = ipvprev->getNonLocalDeviatoricPlasticStrain() + fact*DeltaHatDNonLocal;
ipvcur->getRefToNonLocalMatrixPlasticStrain() = ipvprev->getNonLocalMatrixPlasticStrain() + fact*DeltaHatPNonLocal;
- Tcur = T0 + fact*dT;//------for temperature substepping
-
+ double Tcur = T0 + fact*dT;//------for temperature substepping
+ ipvcur->setTemperature(Tcur); // some dependence on temperature must be followed
+
//elastic predictor
- Fp0 = &ipvTemp.getConstRefToFp();
- // Fp
- ipvcur->getRefToFp() = ipvTemp.getConstRefToFp();
-
+ Fp0 = &ipvTemp->getConstRefToFp();
elasticPredictor(Fcur,*Fp0,corKirpr,kcorprDev,ppr,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,Fepr,Cepr,Epr,Lepr,dLepr,estimateStiff,
&Tcur,&DcorKirprDT,&DkcorprDevDT,&dpprDT);
-
//check yield condition
kCorPrEq = sqrt(1.5*kcorprDev.dotprod());
- R = ipvTemp.getConstRefToIPJ2IsotropicHardening().getR();
-
- // void growth
- double subDeltaHatDNonLocal = DeltaHatDNonLocal/(double)numSubStep;
- double subDeltaHatQNonLocal = DeltaHatQNonLocal/(double)numSubStep;
- double subDeltaHatPNonLocal = DeltaHatPNonLocal/(double)numSubStep;
-
-
- nucleation(ipvcur->getNonLocalMatrixPlasticStrain(),ipvTemp.getNonLocalMatrixPlasticStrain(),ipvprev->getYieldPorosity(),ipvcur,&ipvTemp);
-
- localPorosityGrowth(DeltafV,subDeltaHatDNonLocal,subDeltaHatQNonLocal,subDeltaHatPNonLocal,
- corKirpr,&ipvTemp,ipvcur);
- // internal data
- updateInternalVariables(ipvcur,&ipvTemp,0.,0.,0.,DeltafV,&Tcur);
-
- // update yield porosity as it depends on nonlocal variable
- tildefVstar = ipvTemp.getCorrectedPorosity() + DfVtildeStarDtildefVPrev*DeltafV;
- if (plasticFollowedBlocked){
- yieldfV = ipvprev->getYieldPorosity();
- DyieldfVDDtildefV = 0.;
- }
- else{
- yieldfV = _correctedRegularizedFunction->getVal(tildefVstar);
- DyieldfVDDtildefV = _correctedRegularizedFunction->getDiff(tildefVstar)*DfVtildeStarDtildefVPrev;
- }
+ R = ipvTemp->getConstRefToIPJ2IsotropicHardening().getR();
+ voidCharacteristicsEvoltion_multipleNonLocalVariables(corKirpr,ipvTemp,ipvcur,DyieldfVDNonlocalVar,DyieldfVDKcorpr,plasticFollowedBlocked);
- _voidEvolutionLaw->evolve(ipvTemp.getYieldPorosity(),ipvcur->getYieldPorosity(),subDeltaHatDNonLocal,subDeltaHatQNonLocal,subDeltaHatPNonLocal,
- &ipvTemp,ipvcur, ipvTemp.getConstRefToIPVoidState(), ipvcur->getRefToIPVoidState());
- _coalescenceLaw->computeIPCoalescence(&ipvTemp,ipvcur,&ipvTemp.getConstRefToIPCoalescence(),&ipvcur->getRefToIPCoalescence());
-
-
- yieldF = yieldFunction(kCorPrEq,ppr,R,yieldfV,&ipvTemp,ipvcur,&Tcur);
+ yieldfV = ipvcur->getYieldPorosity(); // does not change during corrector
+ yieldF = yieldFunction(kCorPrEq,ppr,R,yieldfV,ipvTemp,ipvcur,&Tcur);
plastic = (yieldF > _tol);
if (plastic){
- correctorOK = plasticCorrector_multipleNonLocalVariables(Fcur,yieldfV,DyieldfVDDtildefV,
- corKirpr,kcorprDev,ppr,&ipvTemp,ipvcur,
+ correctorOK = plasticCorrector_multipleNonLocalVariables(Fcur,
+ yieldfV,DyieldfVDNonlocalVar,DyieldfVDKcorpr,
+ corKirpr,kcorprDev,ppr,ipvTemp,ipvcur,
Fe,Ce,Ee,Le,dLe,estimateStiff,DcorKirprDEpr,DkcorprDevDEpr,DpprDEpr,
DcorKirDEpr,DcorKirDNonLocalVar,
DFpDEpr,DFpDNonLocalVar,
@@ -4582,10 +4582,9 @@ void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const ST
break;
}
}
-
// next step
if (iter < numSubStep){
- ipvTemp.operator =(*dynamic_cast<const IPVariable*>(ipvcur));
+ ipvTemp->operator =(*static_cast<const IPVariable*>(ipvcur));
}
}
if (!success){
@@ -4610,7 +4609,11 @@ void mlawNonLocalPorosity::predictorCorrector_multipleNonLocalVariables(const ST
}
else{
//reset nucleation
- nucleation(ipvprev->getNonLocalMatrixPlasticStrain(),ipvprev->getNonLocalMatrixPlasticStrain(),ipvprev->getYieldPorosity(),ipvcur,ipvprev);
+ // no nucleation
+ for (int i=0; i< _gdnLawContainer.size(); i++)
+ {
+ ipvcur->getRefToIPNucleation(i).resetNucleation(ipvprev->getConstRefToIPNucleation(i));
+ }
}
}
}
@@ -4743,24 +4746,7 @@ void mlawNonLocalPorosity::I1J2J3_updateInternalVariables(IPNonLocalPorosity *q1
STensorOperation::decomposeDevTr(DeltaEp, devDeltaEp, DeltaHatQ);
DeltaHatD = sqrt(2.*STensorOperation::doubledot(devDeltaEp,devDeltaEp)/3.);
-
- double& hatP = q1->getRefToMatrixPlasticStrain();
- hatP = hatPn + DeltaHatP;
-
- double& hatQ = q1->getRefToVolumetricPlasticStrain();
- hatQ = hatQn + DeltaHatQ;
-
- double& hatD = q1->getRefToDeviatoricPlasticStrain();
- hatD = hatDn + DeltaHatD;
-
- // update hardening
- if (isThermomechanicallyCoupled()){
- _j2IH->hardening(hatPn,q0->getConstRefToIPJ2IsotropicHardening(),hatP,*T, q1->getRefToIPJ2IsotropicHardening());
- }
- else{
- _j2IH->hardening(hatPn,q0->getConstRefToIPJ2IsotropicHardening(),hatP,q1->getRefToIPJ2IsotropicHardening());
- }
-
+ updateInternalVariables(q1,q0,DeltaHatD,DeltaHatQ,DeltaHatP);
};
void mlawNonLocalPorosity::I1J2J3_voidCharacteristicsEvoltion_NonLocal(const STensor3& kcor,
@@ -4777,7 +4763,6 @@ void mlawNonLocalPorosity::I1J2J3_voidCharacteristicsEvoltion_NonLocal(const STe
double DeltaHatPNonLocal = nonlocalVars[1] - nonlocalVarsPrev[1];
double DeltaHatDNonLocal = nonlocalVars[2] - nonlocalVarsPrev[2];
- double DeltafVtilde;
double DDeltafVtildeDDeltaHatDNonLocal, DDeltafVtildeDDeltaHatQNonLocal, DDeltafVtildeDDeltaHatPNonLocal;
static STensor3 DDeltafVtildeDKcor;
@@ -4785,25 +4770,24 @@ void mlawNonLocalPorosity::I1J2J3_voidCharacteristicsEvoltion_NonLocal(const STe
nucleation(q1->getNonLocalMatrixPlasticStrain(),q0->getNonLocalMatrixPlasticStrain(),q0->getYieldPorosity(),q1,q0);
// nonlocal void growth
- localPorosityGrowth(DeltafVtilde, DeltaHatDNonLocal,DeltaHatQNonLocal,DeltaHatPNonLocal,
- kcor,q0,q1, stiff,
+ localPorosityGrowth(q1, DeltaHatDNonLocal,DeltaHatQNonLocal,DeltaHatPNonLocal,
+ kcor,q0, stiff,
&DDeltafVtildeDDeltaHatDNonLocal, &DDeltafVtildeDDeltaHatQNonLocal,&DDeltafVtildeDDeltaHatPNonLocal,
&DDeltafVtildeDKcor);
//
- double& fV = q1->getRefToLocalPorosity();
+ double fV = q1->getLocalPorosity();
double& fVtilde = q1->getRefToNonLocalPorosity();
double& fVtildeStar = q1->getRefToCorrectedPorosity();
double& yieldfV = q1->getRefToYieldPorosity();
// update void
- fV = q0->getLocalPorosity() + DeltafVtilde;
- fVtilde = q0->getNonLocalPorosity() + DeltafVtilde;
+ fVtilde = fV;
// update yield porosity as it depends on nonlocal variable
- double DfVtildeStarDtildefVPrev = q0->getConstRefToIPCoalescence().getAccelerateRate();
- fVtildeStar = q0->getCorrectedPorosity() + DfVtildeStarDtildefVPrev*DeltafVtilde;
+ double DfVtildeStarDtildefV = q1->getConstRefToIPCoalescence().getAccelerateRate();
+ fVtildeStar = q0->getCorrectedPorosity() + DfVtildeStarDtildefV*(q1->getLocalPorosity() - q0->getLocalPorosity());
yieldfV = _correctedRegularizedFunction->getVal(fVtildeStar);
- double DyieldfVDDfVtilde = _correctedRegularizedFunction->getDiff(fVtildeStar)*DfVtildeStarDtildefVPrev;
+ double DyieldfVDDfVtilde = _correctedRegularizedFunction->getDiff(fVtildeStar)*DfVtildeStarDtildefV;
// update void state
_voidEvolutionLaw->evolve(q0->getYieldPorosity(),q1->getYieldPorosity(),DeltaHatDNonLocal,DeltaHatQNonLocal,DeltaHatPNonLocal,
@@ -5432,8 +5416,11 @@ void mlawNonLocalPorosity::I1J2J3_predictorCorrector_NonLocal(const STensor3& F0
// coalescence from previous state
ipvcur->getRefToIPCoalescence() = (ipvprev->getConstRefToIPCoalescence());
//
- // no nucleation
- this->nucleation(ipvcur->getNonLocalMatrixPlasticStrain(),ipvprev->getNonLocalMatrixPlasticStrain(), ipvprev->getYieldPorosity(),ipvcur,ipvprev);
+ // no nucleation
+ for (int i=0; i< _gdnLawContainer.size(); i++)
+ {
+ ipvcur->getRefToIPNucleation(i).resetNucleation(ipvprev->getConstRefToIPNucleation(i));
+ }
//
const STensor3* Fp0 = &(ipvprev->getConstRefToFp());
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorous.h b/NonLinearSolver/materialLaw/mlawNonLocalPorous.h
index 835c6280073f9d745b7bbef17a9875883198fa23..8b8e5c376f946b599e818652ff738f72d57e1155 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorous.h
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorous.h
@@ -89,7 +89,6 @@ protected:
// Failure management
double _fV_failure; // max porosity in yield surface (value for which fG=0)
double _localfVFailure; // max local porosity at failure (almost 1)
- scalarFunction* _localRegularizedFunction; // Regularisation function for local values
scalarFunction* _correctedRegularizedFunction; // Regularisation function for yield/effective values
double _failedTol; // tolerance applied on max values at which the IP is considered as failed.
@@ -135,7 +134,6 @@ protected:
virtual void setShearPorosityGrowthFactor(const double k);
virtual void setStressTriaxialityFunction_kw(const scalarFunction& fct);
virtual void setFailureTolerance(const double NewFailureTol);
- virtual void setLocalRegularizedFunction(const scalarFunction& fct);
virtual void setCorrectedRegularizedFunction(const scalarFunction& fct);
virtual void setPostBlockedDissipationBehavior(const int method);
virtual void clearCLengthLaw();
@@ -368,27 +366,36 @@ protected:
// Case-dependent functions
virtual double yieldFunction(const double kcorEq, const double pcor, const double R, const double yieldfV,
const IPNonLocalPorosity* q0, const IPNonLocalPorosity* q1, const double* T = NULL) const = 0;
-
- virtual void computeResidual(SVector3& res, STensor3& J,
+
+ // when considereing nonlocal porosity,
+ virtual void computeResidual_NonLocalPorosity(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
+ const double yieldfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
const double* T = NULL
) const = 0;
- virtual void computeDResidual(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
+ virtual void computeDResidual_NonLocalPorosity(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const double DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDKcorpr, const double DyieldfVDtildefV,
+ const STensor43& DKcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const STensor3* DKcorprDT = NULL,
double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
) const = 0;
+
+ virtual void computeResidualLocal(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV, const double DyieldfVDfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
+ const double* T = NULL
+ ) const = 0;
virtual void computeDResidualLocal(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
@@ -403,40 +410,41 @@ protected:
// Multiple nonlocal version
virtual void computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J,
- const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfVtilde,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double* T = NULL
- ) const = 0;
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double* T = NULL
+ ) const = 0;
virtual void computeDResidual_multipleNonLocalVariables(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- std::vector<double> &dres0dtildefV, std::vector<double> &dres1dtildefV, std::vector<double> &dres2dtildefV,
+ std::vector<double> &dres0dNonlocalVars, std::vector<double> &dres1dNonlocalVars, std::vector<double> &dres2dNonlocalVars,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const std::vector<double>& DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDkcorpr, const std::vector<double>& DyieldfVDNonlocalVars,
+ const STensor43& DkcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const STensor3* DKcorprVDT = NULL,
double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
) const = 0;
// Internal state and failure management
virtual double getOnsetCriterion(const IPNonLocalPorosity* q1) const = 0;
- virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const;
- virtual void nucleation(const double p1, const double p0, const double f0, IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const=0;
virtual void checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const = 0;
virtual void checkFailed(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const = 0;
+ virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const;
+ virtual void nucleation(const double p1, const double p0, const double f0, IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const;
+
// void expansion from growth, nucleation and shear
// in the integrated form, it is a function of
// plastic defos as DeltaHatD=deviatoric plastic defo, DeltaHatQ = volumetric plastic defor, DeltaHatP=matrixPlasticDefo
// and kcor
- virtual void localPorosityGrowth(double& DeltafV,
+ virtual void localPorosityGrowth(IPNonLocalPorosity *q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
- const STensor3& Kcor, const IPNonLocalPorosity *q0, const IPNonLocalPorosity *q1,
+ const STensor3& Kcor, const IPNonLocalPorosity *q0,
bool stiff = false,
double* DDeltafVDDeltaHatD=NULL, double* DDeltafVDDeltaHatQ=NULL, double* DDeltafVDDeltaHatP=NULL, STensor3* DDeltafVDKcor=NULL) const;
@@ -506,12 +514,12 @@ protected:
// Update internal variables during corrector iterations
virtual void updateInternalVariables(IPNonLocalPorosity *q1, const IPNonLocalPorosity* q0,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, const double DeltafV,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T = NULL) const;
// Update plastic state at the end of a plastic corrector step (and check for coalescence)
virtual void updatePlasticState(IPNonLocalPorosity *q1, const IPNonLocalPorosity* q0, const STensor3& F1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, const double DeltafV,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double yield,
STensor3& Fe1, STensor3& Ce, STensor3& Ee, STensor43& Le, STensor63& dLe,
STensor3& DGNp, STensor3& devDGNp, double& trDGNp, STensor43& Dexp,
@@ -555,8 +563,15 @@ protected:
double* dLocalPorosityDT = NULL,
const double T0 = 0.,const double T = 0.
) const;
-
- virtual bool plasticCorrector_NonLocalPorosity(const STensor3& F1, const double& tildefVstar,const double& DtildefVstarDtildefV,
+
+ virtual void voidCharacteristicsEvoltion_NonLocalPorosity(const STensor3& kcor,
+ const IPNonLocalPorosity *q0, IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
+ double& DfVDDeltaHatD, double& DfVDDeltaHatQ, double& DfVDDeltaHatP,
+ STensor3& DfVDkcorpr) const;
+
+ virtual bool plasticCorrector_NonLocalPorosity(const STensor3& F1,
+ const double& yieldfV,const double& DyieldfVDtildefV,
const STensor3& Kcorpr, const STensor3& devKcorpr, const double& pcorpr,
const IPNonLocalPorosity *q0, IPNonLocalPorosity *q1,
STensor3& Fe, STensor3& Ce, STensor3& Ee, STensor43& Le, STensor63& dLe,
@@ -571,7 +586,13 @@ protected:
) const;
-
+ virtual void voidCharacteristicsEvoltionLocal(const STensor3& kcor,
+ const IPNonLocalPorosity *q0, IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
+ const double DtildefVstarDDeltafV,
+ double& DyieldfVDDeltafV,
+ double& DDeltafVDDeltaHatD, double& DDeltafVDDeltaHatQ, double& DDeltafVDDeltaHatP,
+ STensor3& DDeltafVDKcor) const;
// Predictor-corrector and plastic corrector for local case
virtual void predictorCorrectorLocal(const STensor3& F0, // initial deformation gradient (input @ time n)
const STensor3& Fn, // updated deformation gradient (input @ time n+1)
@@ -743,8 +764,13 @@ protected:
STensor3* DKcorDT= NULL, STensor3* DFpDT= NULL
) const;
-
- virtual bool plasticCorrector_multipleNonLocalVariables(const STensor3& F1, const double& yieldfV,const double& DyieldfVDDeltafV,
+ virtual void voidCharacteristicsEvoltion_multipleNonLocalVariables(const STensor3& kcor,
+ const IPNonLocalPorosity *q0, IPNonLocalPorosity* q1,
+ std::vector<double>& DyieldfVDNonlocalVars,
+ STensor3& DyieldfVDKcor, const bool blockVoidGrowth
+ ) const;
+ virtual bool plasticCorrector_multipleNonLocalVariables(const STensor3& F1,
+ const double& yieldfV,const std::vector<double>& DyieldfVDNonlocalVar, const STensor3& DyieldfVDKcorpr,
const STensor3& Kcorpr, const STensor3& devKcorpr, const double& pcorpr,
const IPNonLocalPorosity *q0, IPNonLocalPorosity *q1,
STensor3& Fe, STensor3& Ce, STensor3& Ee, STensor43& Le, STensor63& dLe,
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.cpp b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.cpp
index fc9c1aaa599451392de1499b52df752f15c532cc..cd7b9b0069552d7561ebcd5945001cc8b6690260 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.cpp
@@ -127,9 +127,6 @@ double mlawNonLocalPorousThomasonLaw::getOnsetCriterion(const IPNonLocalPorosity
return q1->getConstRefToIPVoidState().getVoidLigamentRatio();
};
-void mlawNonLocalPorousThomasonLaw::nucleation(const double p1, const double p0, const double f0, IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const{
- // do nothing
-};
void mlawNonLocalPorousThomasonLaw::checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const{
// do nothing
@@ -424,12 +421,11 @@ void mlawNonLocalPorousThomasonLaw::plasticFlowDerivatives(fullVector<double>& g
};
-void mlawNonLocalPorousThomasonLaw::computeResidual(SVector3& res, STensor3& J,
+void mlawNonLocalPorousThomasonLaw::computeResidual_NonLocalPorosity(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
+ const double yieldfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
const double* T
) const
{
@@ -439,19 +435,19 @@ void mlawNonLocalPorousThomasonLaw::computeResidual(SVector3& res, STensor3& J,
// Get values
double Chi = voidState->getVoidLigamentRatio();
- double DChiDDeltafV = voidState->getDVoidLigamentRatioDYieldPorosity()*DyieldfVDfV;
+ double DChiDyieldfV = voidState->getDVoidLigamentRatioDYieldPorosity();
double DChiDDeltaHatD = voidState->getDVoidLigamentRatioDDeviatoricPlasticDeformation();
double DChiDDeltaHatQ = voidState->getDVoidLigamentRatioDVolumetricPlasticDeformation();
double DChiDDeltaHatP = voidState->getDVoidLigamentRatioDMatrixPlasticDeformation();
double W = voidState->getVoidAspectRatio();
- double DWDDeltafV = voidState->getDVoidAspectRatioDYieldPorosity()*DyieldfVDfV;
+ double DWDyieldfV = voidState->getDVoidAspectRatioDYieldPorosity();
double DWDDeltaHatD = voidState->getDVoidAspectRatioDDeviatoricPlasticDeformation();
double DWDDeltaHatQ = voidState->getDVoidAspectRatioDVolumetricPlasticDeformation();
double DWDDeltaHatP = voidState->getDVoidAspectRatioDMatrixPlasticDeformation();
double gamma = voidState->getVoidShapeFactor();
- double DgammaDDeltafV = voidState->getDVoidShapeFactorDYieldPorosity()*DyieldfVDfV;
+ double DgammaDyieldfV = voidState->getDVoidShapeFactorDYieldPorosity();
double DgammaDDeltaHatD = voidState->getDVoidShapeFactorDDeviatoricPlasticDeformation();
double DgammaDDeltaHatQ = voidState->getDVoidShapeFactorDVolumetricPlasticDeformation();
double DgammaDDeltaHatP = voidState->getDVoidShapeFactorDMatrixPlasticDeformation();
@@ -506,34 +502,25 @@ void mlawNonLocalPorousThomasonLaw::computeResidual(SVector3& res, STensor3& J,
computeYieldDerivatives(gradf,kcorEq,pcor,R,Chi,W,q0,q1,T);
plasticFlowDerivatives(gradNs,gradNv,kcorEq,pcor,R,Chi,W,q0,q1,T);
- double Dres0DfV = gradf(3)*DChiDDeltafV + gradf(4)*DWDDeltafV;
-
// Dres0DDeltaHatD
- J(0,0) = gradf(0)*DkcorEqDDeltaHatD+ gradf(3)*DChiDDeltaHatD + gradf(4)*DWDDeltaHatD + Dres0DfV*DDeltafVDDeltaHatD;
+ J(0,0) = gradf(0)*DkcorEqDDeltaHatD+ gradf(3)*DChiDDeltaHatD + gradf(4)*DWDDeltaHatD;
// Dres0DDeltaHatQ
- J(0,1) = gradf(1)*DpcorDDeltaHatQ + gradf(3)*DChiDDeltaHatQ + gradf(4)*DWDDeltaHatQ+ Dres0DfV*DDeltafVDDeltaHatQ;
+ J(0,1) = gradf(1)*DpcorDDeltaHatQ + gradf(3)*DChiDDeltaHatQ + gradf(4)*DWDDeltaHatQ;
// Dres0DDeltaHatP
- J(0,2) = gradf(2)*H+ + gradf(3)*DChiDDeltaHatP + gradf(4)*DWDDeltaHatP + Dres0DfV*DDeltafVDDeltaHatP;
-
- double Dres1DfV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDDeltafV +
- R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltafV;
-
+ J(0,2) = gradf(2)*H+ + gradf(3)*DChiDDeltaHatP + gradf(4)*DWDDeltaHatP;
// Dres1DDeltaHatD
J(1,0) = R0*(Nv + DeltaHatD*gradNv(0)*DkcorEqDDeltaHatD - DeltaHatQ*gradNs(0)*DkcorEqDDeltaHatD) +
R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDDeltaHatD +
- R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltaHatD +
- Dres1DfV*DDeltafVDDeltaHatD;
+ R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltaHatD;
// Dres1DDeltaHatQ
J(1,1) = R0*(DeltaHatD*gradNv(1)*DpcorDDeltaHatQ- Ns - DeltaHatQ*gradNs(1)*DpcorDDeltaHatQ) +
R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDDeltaHatQ +
- R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltaHatQ +
- Dres1DfV*DDeltafVDDeltaHatQ;
+ R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltaHatQ;
// Dres1DDeltaHatP
- J(1,2) = R0*(DeltaHatD*gradNv(2) - DeltaHatQ*gradNs(2))*H +
+ J(1,2) = R0*(DeltaHatD*gradNv(2) - DeltaHatQ*gradNs(2))*H+
R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDDeltaHatP+
- R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltaHatP+
- Dres1DfV*DDeltafVDDeltaHatP;
+ R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltaHatP;
if (_balanceEquationType == 1){
// Dres2DDeltaHatd
@@ -544,13 +531,12 @@ void mlawNonLocalPorousThomasonLaw::computeResidual(SVector3& res, STensor3& J,
J(2,2) = (1-f0)*(R + H*DeltaHatP)/R0;
}
else if (_balanceEquationType == 2){
- double Dres2DfV =(2.*DgammaDDeltafV*yieldfV*DeltaHatP - 2.*Chi*DChiDDeltafV*DeltaHatD + 2.*gamma*DeltaHatP*DyieldfVDfV)/f0;
// Dres2DDeltaHatD
- J(2,0) = (2.*DgammaDDeltaHatD*yieldfV*DeltaHatP - Chi*Chi- 2.*Chi*DChiDDeltaHatD*DeltaHatD)/f0 + Dres2DfV*DDeltafVDDeltaHatD;
+ J(2,0) = (2.*DgammaDDeltaHatD*yieldfV*DeltaHatP - Chi*Chi- 2.*Chi*DChiDDeltaHatD*DeltaHatD)/f0;
// Dres2DDeltaHatQ
- J(2,1) = (2.*DgammaDDeltaHatQ*yieldfV*DeltaHatP - 2.*Chi*DChiDDeltaHatQ*DeltaHatD)/f0 + Dres2DfV*DDeltafVDDeltaHatQ;
+ J(2,1) = (2.*DgammaDDeltaHatQ*yieldfV*DeltaHatP - 2.*Chi*DChiDDeltaHatQ*DeltaHatD)/f0;
// Dres2DDeltaHatP
- J(2,2) = (2.*DgammaDDeltaHatP*yieldfV*DeltaHatP+2.*gamma*yieldfV - 2.*Chi*DChiDDeltaHatP*DeltaHatD)/f0 + Dres2DfV*DDeltafVDDeltaHatP;
+ J(2,2) = (2.*DgammaDDeltaHatP*yieldfV*DeltaHatP+2.*gamma*yieldfV - 2.*Chi*DChiDDeltaHatP*DeltaHatD)/f0;
}
else{
Msg::Fatal("balance equation is not defined");
@@ -559,31 +545,28 @@ void mlawNonLocalPorousThomasonLaw::computeResidual(SVector3& res, STensor3& J,
};
-void mlawNonLocalPorousThomasonLaw::computeDResidual(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
+void mlawNonLocalPorousThomasonLaw::computeDResidual_NonLocalPorosity(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const double DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDKcorpr, const double DyieldfVDtildefV,
+ const STensor43& DKcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T, const double* DRDT,
- const STensor3* DdevKcorprDT , const double* DpcorprDT, const double* DDeltafVDT,
+ const STensor3* DdevKcorprDT , const double* DpcorprDT, const STensor3* DKcorprDT,
double *dres0dT, double *dres1dT, double *dres2dT
) const{
const IPVoidState* voidState = &q1->getConstRefToIPVoidState();
double Chi = voidState->getVoidLigamentRatio();
- double DChiDDeltafV = voidState->getDVoidLigamentRatioDYieldPorosity()*DyieldfVDfV;
- double DChiDtildefV = voidState->getDVoidLigamentRatioDYieldPorosity()*DyieldfVDtildefV;
+ double DChiDyieldfV = voidState->getDVoidLigamentRatioDYieldPorosity();
double W = voidState->getVoidAspectRatio();
- double DWDDeltafV = voidState->getDVoidAspectRatioDYieldPorosity()*DyieldfVDfV;
- double DWDtildefV = voidState->getDVoidAspectRatioDYieldPorosity()*DyieldfVDtildefV;
+ double DWDyieldfV = voidState->getDVoidAspectRatioDYieldPorosity();
double gamma = voidState->getVoidShapeFactor();
- double DgammaDDeltafV = voidState->getDVoidShapeFactorDYieldPorosity()*DyieldfVDfV;
- double DgammaDtildefV = voidState->getDVoidShapeFactorDYieldPorosity()*DyieldfVDtildefV;
+ double DgammaDyieldfV = voidState->getDVoidShapeFactorDYieldPorosity();
double f0 = q1->getInitialPorosity();
double R0 = _j2IH->getYield0(); // reference value
@@ -594,7 +577,10 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual(STensor3 &dres0dEpr, STenso
KT=bulkModulus(*T);
muT=shearModulus(*T);
}
-
+
+ static STensor3 DyieldfVDEpr;
+ STensorOperation::multSTensor3STensor43(DyieldfVDKcorpr,DKcorprDEpr,DyieldfVDEpr);
+
static STensor3 DkcorEqprDEpr;
STensorOperation::multSTensor3STensor43(devKcorpr,DdevKcorprDEpr,DkcorEqprDEpr);
DkcorEqprDEpr *= (1.5/kcorEqpr);
@@ -614,29 +600,29 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual(STensor3 &dres0dEpr, STenso
// compute Dres0DEpr
double Dres0DkcorEq = gradf(0);
double Dres0Dpcor = gradf(1);
- double Dres0DDeltafV = gradf(3)*DChiDDeltafV+gradf(4)*DWDDeltafV;
+ double Dres0DyieldfV = gradf(3)*DChiDyieldfV+gradf(4)*DWDyieldfV;
dres0dEpr = DkcorEqDEpr;
dres0dEpr*= Dres0DkcorEq;
dres0dEpr.daxpy(DpcorDEpr,Dres0Dpcor);
- dres0dEpr.daxpy(DDeltafVDEpr,Dres0DDeltafV);
+ dres0dEpr.daxpy(DyieldfVDEpr,Dres0DyieldfV);
// compute Dres1DEpr
double Dres1DkcorEq = R0*(DeltaHatD*gradNv(0)- DeltaHatQ*gradNs(0));
double Dres1Dpcor = R0*(DeltaHatD*gradNv(1)- DeltaHatQ*gradNs(1));
- double Dres1DDeltafV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDDeltafV +
- R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltafV;
+ double Dres1DyieldfV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDyieldfV +
+ R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDyieldfV;
dres1dEpr = DkcorEqDEpr;
dres1dEpr *= Dres1DkcorEq;
dres1dEpr.daxpy(DpcorDEpr,Dres1Dpcor);
- dres1dEpr.daxpy(DDeltafVDEpr,Dres1DDeltafV);
+ dres1dEpr.daxpy(DyieldfVDEpr,Dres1DyieldfV);
// compute Dres2DEpr
double Dres2DkcorEq = 0.;
double Dres2Dpcor = 0.;
double Dres2DChi = 0.;
- double Dres2DDeltafV = 0.;
+ double Dres2DyieldfV = 0.;
if (_balanceEquationType == 1){
Dres2DkcorEq = (-DeltaHatD)/R0;
@@ -647,23 +633,22 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual(STensor3 &dres0dEpr, STenso
dres2dEpr.daxpy(DpcorDEpr,Dres2Dpcor);
}
else if (_balanceEquationType == 2){
- Dres2DDeltafV = (2.*DgammaDDeltafV*yieldfV*DeltaHatP - 2.*Chi*DChiDDeltafV*DeltaHatD + 2.*gamma*DeltaHatP*DyieldfVDfV)/f0;
- dres2dEpr = DDeltafVDEpr;
- dres2dEpr *= Dres2DDeltafV;
+ Dres2DyieldfV = (2.*DgammaDyieldfV*yieldfV*DeltaHatP - 2.*Chi*DChiDyieldfV*DeltaHatD + 2.*gamma*DeltaHatP)/f0;
+ dres2dEpr = DyieldfVDEpr;
+ dres2dEpr *= Dres2DyieldfV;
}
else{
Msg::Fatal("balance equation is not defined");
}
// compute nonlocal derivative
- dres0dtildefV = gradf(3)*DChiDtildefV+gradf(4)*DWDtildefV;
- dres1dtildefV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDtildefV +
- R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDtildefV;
+ dres0dtildefV = (gradf(3)*DChiDyieldfV+gradf(4)*DWDyieldfV)*DyieldfVDtildefV;
+ dres1dtildefV = (R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDyieldfV + R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDyieldfV)*DyieldfVDtildefV;
if (_balanceEquationType == 1){
dres2dtildefV = 0.;
}
else if (_balanceEquationType == 2){
- dres2dtildefV = (2.*DgammaDtildefV*yieldfV*DeltaHatP - 2.*Chi*DChiDtildefV*DeltaHatD + 2.*gamma*DeltaHatP*DyieldfVDtildefV)/f0;
+ dres2dtildefV = (2.*DgammaDyieldfV*yieldfV*DeltaHatP - 2.*Chi*DChiDyieldfV*DeltaHatD + 2.*gamma*DeltaHatP)*DyieldfVDtildefV/f0;
}
else{
Msg::Fatal("balance equation is not defined");
@@ -673,11 +658,11 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual(STensor3 &dres0dEpr, STenso
double dKdT=bulkModulusDerivative(*T);
double dmudT=shearModulusDerivative(*T);
+
static double DkcorEqDT, DpcorDT;
DkcorEqDT = STensorOperation::doubledot(devKcorpr,*DdevKcorprDT)*1.5/kcorEqpr - 3.*dmudT*DeltaHatD;
DpcorDT = *DpcorprDT -dKdT* DeltaHatQ;
-
double partDres0DT = gradf(2)*(*DRDT) + gradf(5); // yield temperature-dependent
*dres0dT = Dres0DkcorEq*DkcorEqDT + Dres0Dpcor*DpcorDT + partDres0DT;
@@ -689,7 +674,8 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual(STensor3 &dres0dEpr, STenso
*dres2dT = Dres2DkcorEq*DkcorEqDT + Dres2Dpcor*DpcorDT+ partDres2DT;
}
else if (_balanceEquationType == 2){
- *dres2dT = Dres2DDeltafV*(*DDeltafVDT);
+ double DyieldfVDT = STensorOperation::doubledot(DyieldfVDKcorpr,*DKcorprDT);
+ *dres2dT = Dres2DyieldfV*DyieldfVDT;
}
else{
Msg::Fatal("balance equation is not defined");
@@ -697,6 +683,140 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual(STensor3 &dres0dEpr, STenso
}
};
+void mlawNonLocalPorousThomasonLaw::computeResidualLocal(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV, const double DyieldfVDfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
+ const double* T
+ ) const
+{
+ // Get ipv
+ const IPVoidState* voidStatePrev = &q0->getConstRefToIPVoidState();
+ const IPVoidState* voidState = &q1->getConstRefToIPVoidState();
+
+ // Get values
+ double Chi = voidState->getVoidLigamentRatio();
+ double DChiDDeltafV = voidState->getDVoidLigamentRatioDYieldPorosity()*DyieldfVDfV;
+ double DChiDDeltaHatD = voidState->getDVoidLigamentRatioDDeviatoricPlasticDeformation();
+ double DChiDDeltaHatQ = voidState->getDVoidLigamentRatioDVolumetricPlasticDeformation();
+ double DChiDDeltaHatP = voidState->getDVoidLigamentRatioDMatrixPlasticDeformation();
+
+ double W = voidState->getVoidAspectRatio();
+ double DWDDeltafV = voidState->getDVoidAspectRatioDYieldPorosity()*DyieldfVDfV;
+ double DWDDeltaHatD = voidState->getDVoidAspectRatioDDeviatoricPlasticDeformation();
+ double DWDDeltaHatQ = voidState->getDVoidAspectRatioDVolumetricPlasticDeformation();
+ double DWDDeltaHatP = voidState->getDVoidAspectRatioDMatrixPlasticDeformation();
+
+ double gamma = voidState->getVoidShapeFactor();
+ double DgammaDDeltafV = voidState->getDVoidShapeFactorDYieldPorosity()*DyieldfVDfV;
+ double DgammaDDeltaHatD = voidState->getDVoidShapeFactorDDeviatoricPlasticDeformation();
+ double DgammaDDeltaHatQ = voidState->getDVoidShapeFactorDVolumetricPlasticDeformation();
+ double DgammaDDeltaHatP = voidState->getDVoidShapeFactorDMatrixPlasticDeformation();
+
+
+ double yieldfVprev = q0->getYieldPorosity();
+ double Chiprev = voidStatePrev->getVoidLigamentRatio();
+ double gammprev = voidStatePrev->getVoidShapeFactor();
+
+ double f0 = q1->getInitialPorosity();
+ double R0 = _j2IH->getYield0();
+
+
+ // Get elastic parameters
+ double KT= _K;
+ double muT = _mu;
+ if (isThermomechanicallyCoupled()){
+ KT = bulkModulus(*T);
+ muT =shearModulus(*T);
+ };
+
+ // Corrected stress invariants
+ double kcorEq = kcorEqpr - 3.*muT*DeltaHatD;
+ double pcor = pcorpr - KT*DeltaHatQ;
+
+
+ // Residues
+ // 1. yield function
+ res(0) = yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
+
+
+ // 2. plastic flow normality
+ double Ns, Nv; // plastic flow components
+ plasticFlow(Ns,Nv,kcorEq,pcor,R,Chi,W,q0,q1,T);
+ res(1) = R0*(DeltaHatD*Nv- DeltaHatQ*Ns);
+
+ // plastic energy balance
+ if (_balanceEquationType == 1){
+ res(2) = ((1-f0)*R*DeltaHatP - kcorEq*DeltaHatD - pcor*DeltaHatQ)/R0;
+ }
+ else if (_balanceEquationType == 2){
+ res(2) = (2.*gamma*yieldfV*DeltaHatP - (Chi*Chi)*DeltaHatD)/f0;
+ }
+ else{
+ Msg::Fatal("balance equation is not defined");
+ }
+
+ double DkcorEqDDeltaHatD = -3.*muT;
+ double DpcorDDeltaHatQ = -KT;
+
+ static fullVector<double> gradf, gradNs, gradNv;
+ computeYieldDerivatives(gradf,kcorEq,pcor,R,Chi,W,q0,q1,T);
+ plasticFlowDerivatives(gradNs,gradNv,kcorEq,pcor,R,Chi,W,q0,q1,T);
+
+ double Dres0DfV = gradf(3)*DChiDDeltafV + gradf(4)*DWDDeltafV;
+
+ // Dres0DDeltaHatD
+ J(0,0) = gradf(0)*DkcorEqDDeltaHatD+ gradf(3)*DChiDDeltaHatD + gradf(4)*DWDDeltaHatD + Dres0DfV*DDeltafVDDeltaHatD;
+ // Dres0DDeltaHatQ
+ J(0,1) = gradf(1)*DpcorDDeltaHatQ + gradf(3)*DChiDDeltaHatQ + gradf(4)*DWDDeltaHatQ+ Dres0DfV*DDeltafVDDeltaHatQ;
+ // Dres0DDeltaHatP
+ J(0,2) = gradf(2)*H+ + gradf(3)*DChiDDeltaHatP + gradf(4)*DWDDeltaHatP + Dres0DfV*DDeltafVDDeltaHatP;
+
+ double Dres1DfV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDDeltafV +
+ R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltafV;
+
+
+ // Dres1DDeltaHatD
+ J(1,0) = R0*(Nv + DeltaHatD*gradNv(0)*DkcorEqDDeltaHatD - DeltaHatQ*gradNs(0)*DkcorEqDDeltaHatD) +
+ R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDDeltaHatD +
+ R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltaHatD +
+ Dres1DfV*DDeltafVDDeltaHatD;
+ // Dres1DDeltaHatQ
+ J(1,1) = R0*(DeltaHatD*gradNv(1)*DpcorDDeltaHatQ- Ns - DeltaHatQ*gradNs(1)*DpcorDDeltaHatQ) +
+ R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDDeltaHatQ +
+ R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltaHatQ +
+ Dres1DfV*DDeltafVDDeltaHatQ;
+ // Dres1DDeltaHatP
+ J(1,2) = R0*(DeltaHatD*gradNv(2) - DeltaHatQ*gradNs(2))*H +
+ R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDDeltaHatP+
+ R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltaHatP+
+ Dres1DfV*DDeltafVDDeltaHatP;
+
+ if (_balanceEquationType == 1){
+ // Dres2DDeltaHatd
+ J(2,0) = (-kcorEq - DkcorEqDDeltaHatD*DeltaHatD)/R0;
+ // Dres2DDeltaHatq
+ J(2,1) = (-pcor - DpcorDDeltaHatQ*DeltaHatQ)/R0;
+ // Dres2DDeltaHatp
+ J(2,2) = (1-f0)*(R + H*DeltaHatP)/R0;
+ }
+ else if (_balanceEquationType == 2){
+ double Dres2DfV =(2.*DgammaDDeltafV*yieldfV*DeltaHatP - 2.*Chi*DChiDDeltafV*DeltaHatD + 2.*gamma*DeltaHatP*DyieldfVDfV)/f0;
+ // Dres2DDeltaHatD
+ J(2,0) = (2.*DgammaDDeltaHatD*yieldfV*DeltaHatP - Chi*Chi- 2.*Chi*DChiDDeltaHatD*DeltaHatD)/f0 + Dres2DfV*DDeltafVDDeltaHatD;
+ // Dres2DDeltaHatQ
+ J(2,1) = (2.*DgammaDDeltaHatQ*yieldfV*DeltaHatP - 2.*Chi*DChiDDeltaHatQ*DeltaHatD)/f0 + Dres2DfV*DDeltafVDDeltaHatQ;
+ // Dres2DDeltaHatP
+ J(2,2) = (2.*DgammaDDeltaHatP*yieldfV*DeltaHatP+2.*gamma*yieldfV - 2.*Chi*DChiDDeltaHatP*DeltaHatD)/f0 + Dres2DfV*DDeltafVDDeltaHatP;
+ }
+ else{
+ Msg::Fatal("balance equation is not defined");
+ }
+
+
+};
void mlawNonLocalPorousThomasonLaw::computeDResidualLocal(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
@@ -821,10 +941,10 @@ void mlawNonLocalPorousThomasonLaw::computeDResidualLocal(STensor3 &dres0dEpr, S
void mlawNonLocalPorousThomasonLaw::computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
+ const double yieldfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double* T
+ const double* T
) const{
// note that Chi, W, gamma are nonlocal
@@ -919,27 +1039,27 @@ void mlawNonLocalPorousThomasonLaw::computeResidual_multipleNonLocalVariables(SV
void mlawNonLocalPorousThomasonLaw::computeDResidual_multipleNonLocalVariables(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- std::vector<double> &dres0dtildefV, std::vector<double> &dres1dtildefV, std::vector<double> &dres2dtildefV,
+ std::vector<double> &dres0dNonlocalVars, std::vector<double> &dres1dNonlocalVars, std::vector<double> &dres2dNonlocalVars,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const std::vector<double>& DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDkcorpr, const std::vector<double>& DyieldfVDNonlocalVars,
+ const STensor43& DkcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T, const double* DRDT,
- const STensor3* DdevKcorprDT , const double* DpcorprDT, const double* DDeltafVDT,
+ const STensor3* DdevKcorprDT , const double* DpcorprDT, const STensor3* DKcorprVDT,
double *dres0dT, double *dres1dT, double *dres2dT
) const{
const IPVoidState* voidState = &q1->getConstRefToIPVoidState();
double Chi = voidState->getVoidLigamentRatio();
- double DChiDDeltafV = voidState->getDVoidLigamentRatioDYieldPorosity()*DyieldfVDfV;
+ double DChiDyieldfV = voidState->getDVoidLigamentRatioDYieldPorosity();
double W = voidState->getVoidAspectRatio();
- double DWDDeltafV = voidState->getDVoidAspectRatioDYieldPorosity()*DyieldfVDfV;
+ double DWDyieldfV = voidState->getDVoidAspectRatioDYieldPorosity();
double gamma = voidState->getVoidShapeFactor();
- double DgammaDDeltafV = voidState->getDVoidShapeFactorDYieldPorosity()*DyieldfVDfV;
+ double DgammaDyieldfV = voidState->getDVoidShapeFactorDYieldPorosity();
double f0 = q1->getInitialPorosity();
double R0 = _j2IH->getYield0(); // reference value
@@ -950,6 +1070,9 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual_multipleNonLocalVariables(S
KT=bulkModulus(*T);
muT=shearModulus(*T);
}
+
+ static STensor3 DyieldfVDEpr;
+ STensorOperation::multSTensor3STensor43(DyieldfVDkcorpr,DkcorprDEpr,DyieldfVDEpr);
static STensor3 DkcorEqprDEpr;
STensorOperation::multSTensor3STensor43(devKcorpr,DdevKcorprDEpr,DkcorEqprDEpr);
@@ -970,31 +1093,32 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual_multipleNonLocalVariables(S
// compute Dres0DEpr
double Dres0DkcorEq = gradf(0);
double Dres0Dpcor = gradf(1);
- double Dres0DDeltafV = gradf(3)*DChiDDeltafV+gradf(4)*DWDDeltafV;
+ double Dres0DyieldfV = gradf(3)*DChiDyieldfV+gradf(4)*DWDyieldfV;
dres0dEpr = DkcorEqDEpr;
dres0dEpr*= Dres0DkcorEq;
dres0dEpr.daxpy(DpcorDEpr,Dres0Dpcor);
- dres0dEpr.daxpy(DDeltafVDEpr,Dres0DDeltafV);
+ dres0dEpr.daxpy(DyieldfVDEpr,Dres0DyieldfV);
// compute Dres1DEpr
double Dres1DkcorEq = R0*(DeltaHatD*gradNv(0)- DeltaHatQ*gradNs(0));
double Dres1Dpcor = R0*(DeltaHatD*gradNv(1)- DeltaHatQ*gradNs(1));
- double Dres1DDeltafV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDDeltafV +
- R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDDeltafV;
+ double Dres1DyieldfV = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDyieldfV +
+ R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDyieldfV;
dres1dEpr = DkcorEqDEpr;
dres1dEpr *= Dres1DkcorEq;
dres1dEpr.daxpy(DpcorDEpr,Dres1Dpcor);
- dres1dEpr.daxpy(DDeltafVDEpr,Dres1DDeltafV);
+ dres1dEpr.daxpy(DyieldfVDEpr,Dres1DyieldfV);
// compute Dres2DEpr
double Dres2DkcorEq = 0.;
double Dres2Dpcor = 0.;
double Dres2DChi = 0.;
- double Dres2DDeltafV = 0.;
+ double Dres2DyieldfV = 0.;
if (_balanceEquationType == 1){
+ // ((1-f0)*R*DeltaHatP - kcorEq*DeltaHatD - pcor*DeltaHatQ)/R0;
Dres2DkcorEq = (-DeltaHatD)/R0;
Dres2Dpcor = (-DeltaHatQ)/R0;
@@ -1003,21 +1127,17 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual_multipleNonLocalVariables(S
dres2dEpr.daxpy(DpcorDEpr,Dres2Dpcor);
}
else if (_balanceEquationType == 2){
- Dres2DDeltafV = (2.*DgammaDDeltafV*yieldfV*DeltaHatP - 2.*Chi*DChiDDeltafV*DeltaHatD+ 2.*gamma*DeltaHatP*DyieldfVDfV)/f0;
+ // (2.*gamma*yieldfV*DeltaHatP - (Chi*Chi)*DeltaHatD)/f0;
+ Dres2DyieldfV = (2.*DgammaDyieldfV*yieldfV*DeltaHatP - 2.*Chi*DChiDyieldfV*DeltaHatD+ 2.*gamma*DeltaHatP)/f0;
- dres2dEpr= DDeltafVDEpr;
- dres2dEpr *= Dres2DDeltafV;
+ dres2dEpr= DyieldfVDEpr;
+ dres2dEpr *= Dres2DyieldfV;
}
else{
Msg::Fatal("balance equation is not defined");
}
- static std::vector<double> DChiDtildefV, DWDtildefV;
- if (DChiDtildefV.size() != _numNonLocalVar){
- DChiDtildefV.resize(_numNonLocalVar);
- DWDtildefV.resize(_numNonLocalVar);
- }
-
+ static std::vector<double> DChiDNonlocalVars(3,0.), DWDNonlocalVars(3,0.);
double DChiDDeltaHatDNonLocal = voidState->getDVoidLigamentRatioDDeviatoricPlasticDeformation();
double DChiDDeltaHatQNonLocal = voidState->getDVoidLigamentRatioDVolumetricPlasticDeformation();
@@ -1027,29 +1147,29 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual_multipleNonLocalVariables(S
double DWDDeltaHatQNonLocal = voidState->getDVoidAspectRatioDVolumetricPlasticDeformation();
double DWDDeltaHatPNonLocal = voidState->getDVoidAspectRatioDMatrixPlasticDeformation();
- DChiDtildefV[0] = DChiDDeltaHatQNonLocal;
- DChiDtildefV[1] = DChiDDeltaHatPNonLocal;
- DChiDtildefV[2] = DChiDDeltaHatDNonLocal;
+ DChiDNonlocalVars[0] = DChiDDeltaHatQNonLocal;
+ DChiDNonlocalVars[1] = DChiDDeltaHatPNonLocal;
+ DChiDNonlocalVars[2] = DChiDDeltaHatDNonLocal;
- DWDtildefV[0] = DWDDeltaHatQNonLocal;
- DWDtildefV[1] = DWDDeltaHatPNonLocal;
- DWDtildefV[2] = DWDDeltaHatDNonLocal;
+ DWDNonlocalVars[0] = DWDDeltaHatQNonLocal;
+ DWDNonlocalVars[1] = DWDDeltaHatPNonLocal;
+ DWDNonlocalVars[2] = DWDDeltaHatDNonLocal;
for (int i=0; i< _numNonLocalVar; i++){
- DChiDtildefV[i] += voidState->getDVoidLigamentRatioDYieldPorosity()*DyieldfVDtildefV[i];
- DWDtildefV[i] += voidState->getDVoidAspectRatioDYieldPorosity()*DyieldfVDtildefV[i];
+ DChiDNonlocalVars[i] += voidState->getDVoidLigamentRatioDYieldPorosity()*DyieldfVDNonlocalVars[i];
+ DWDNonlocalVars[i] += voidState->getDVoidAspectRatioDYieldPorosity()*DyieldfVDNonlocalVars[i];
}
for (int i=0; i< _numNonLocalVar; i++){
// compute nonlocal derivative
- dres0dtildefV[i] = gradf(3)*DChiDtildefV[i]+gradf(4)*DWDtildefV[i];
- dres1dtildefV[i] = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDtildefV[i] +
- R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDtildefV[i];
+ dres0dNonlocalVars[i] = gradf(3)*DChiDNonlocalVars[i]+gradf(4)*DWDNonlocalVars[i];
+ dres1dNonlocalVars[i] = R0*(DeltaHatD*gradNv(3)- DeltaHatQ*gradNs(3))*DChiDNonlocalVars[i] +
+ R0*(DeltaHatD*gradNv(4)- DeltaHatQ*gradNs(4))*DWDNonlocalVars[i];
if (_balanceEquationType == 1){
- dres2dtildefV[i] = 0.;
+ dres2dNonlocalVars[i] = 0.;
}
else if (_balanceEquationType == 2){
- dres2dtildefV[i] = Dres2DChi*DChiDtildefV[i] + 2.*gamma*DeltaHatP*DyieldfVDtildefV[i]/f0;
+ dres2dNonlocalVars[i] = Dres2DChi*DChiDNonlocalVars[i] + 2.*gamma*DeltaHatP*DyieldfVDNonlocalVars[i]/f0;
}
else{
Msg::Fatal("balance equation is not defined");
@@ -1077,7 +1197,8 @@ void mlawNonLocalPorousThomasonLaw::computeDResidual_multipleNonLocalVariables(S
*dres2dT = Dres2DkcorEq*DkcorEqDT + Dres2Dpcor*DpcorDT+ partDres2DT;
}
else if (_balanceEquationType == 2){
- *dres2dT = Dres2DDeltafV*(*DDeltafVDT);
+ double DyieldfVDT= STensorOperation::doubledot(DyieldfVDkcorpr,*DKcorprVDT);
+ *dres2dT = Dres2DyieldfV*DyieldfVDT;
}
else{
Msg::Fatal("balance equation is not defined");
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.h b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.h
index ce84311bba571921ddbaae3c98f62a5547a5821b..71cf120320c511bd507d893972052a1e64aa1901 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.h
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.h
@@ -75,28 +75,37 @@ public:
virtual double yieldFunction(const double kcorEq, const double pcor, const double R, const double yieldfV,
const IPNonLocalPorosity* q0, const IPNonLocalPorosity* q1, const double* T = NULL) const;
- virtual void computeResidual(SVector3& res, STensor3& J,
+ virtual void computeResidual_NonLocalPorosity(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
+ const double yieldfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
const double* T = NULL
) const;
- virtual void computeDResidual(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
+ virtual void computeDResidual_NonLocalPorosity(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const double DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDKcorpr, const double DyieldfVDtildefV,
+ const STensor43& DKcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const STensor3* DKcorprDT = NULL,
double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
) const;
-
+
+
+ virtual void computeResidualLocal(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV, const double DyieldfVDfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
+ const double* T = NULL
+ ) const;
+
virtual void computeDResidualLocal(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
const double yieldfV, const double DyieldfVDfV,
@@ -110,27 +119,26 @@ public:
virtual void computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J,
- const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double* T = NULL
- ) const;
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double* T = NULL
+ ) const;
virtual void computeDResidual_multipleNonLocalVariables(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- std::vector<double> &dres0dtildefV, std::vector<double> &dres1dtildefV, std::vector<double> &dres2dtildefV,
+ std::vector<double> &dres0dNonlocalVars, std::vector<double> &dres1dNonlocalVars, std::vector<double> &dres2dNonlocalVars,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const std::vector<double>& DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDkcorpr, const std::vector<double>& DyieldfVDtildefV,
+ const STensor43& DkcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const STensor3* DKcorprVDT = NULL,
double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
) const;
virtual double getOnsetCriterion(const IPNonLocalPorosity* q1) const;
- virtual void nucleation(const double p1, const double p0, const double f0, IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const;
virtual void checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const;
virtual void checkFailed(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const;
@@ -217,72 +225,83 @@ class mlawNonLocalPorousThomasonLawWithLodeEnhanced : public mlawNonLocalPorousT
virtual double yieldFunction(const double kcorEq, const double pcor, const double R, const double yieldfV,
const IPNonLocalPorosity* q0, const IPNonLocalPorosity* q1, const double* T = NULL) const
{
- Msg::Fatal("mlawNonLocalPorousThomasonLawWithLodeEnhanced cannot be used");
+ Msg::Fatal("mlawNonLocalPorousThomasonLawWithLodeEnhanced::yieldFunction cannot be used");
};
- virtual void computeResidual(SVector3& res, STensor3& J,
- const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
- const double* T = NULL
- ) const
+ virtual void computeResidual_NonLocalPorosity(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double* T = NULL
+ ) const
{
- Msg::Fatal("mlawNonLocalPorousThomasonLawWithLodeEnhanced cannot be used");
+ Msg::Fatal("mlawNonLocalPorousThomasonLawWithLodeEnhanced::computeResidual_NonLocalPorosity cannot be used");
};
- virtual void computeDResidual(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
- const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const double DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
- const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
- double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
- ) const
+ virtual void computeDResidual_NonLocalPorosity(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
+ double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
+ const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
+ const double yieldfV, const STensor3& DyieldfVDKcorpr, const double DyieldfVDtildefV,
+ const STensor43& DKcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
+ const double* T = NULL, const double* DRDT= NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const STensor3* DKcorprDT = NULL,
+ double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
+ ) const
{
- Msg::Fatal("mlawNonLocalPorousThomasonLawWithLodeEnhanced cannot be used");
+ Msg::Fatal("mlawNonLocalPorousThomasonLawWithLodeEnhanced::computeDResidual_NonLocalPorosity cannot be used");
+ };
+
+ virtual void computeResidualLocal(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV, const double DyieldfVDfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
+ const double* T = NULL
+ ) const
+ {
+ Msg::Fatal("mlawNonLocalPorousThomasonLawWithLodeEnhanced::computeResidualLocal cannot be used");
};
virtual void computeDResidualLocal(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
- const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
- double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
- ) const
+ const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
+ const double yieldfV, const double DyieldfVDfV,
+ const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
+ const double* T = NULL, const double* DRDT= NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
+ double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
+ ) const
{
- Msg::Fatal("mlawNonLocalPorousThomasonLawWithLodeEnhanced cannot be used");
+ Msg::Fatal("mlawNonLocalPorousThomasonLawWithLodeEnhanced::computeDResidualLocal cannot be used");
};
// Multiple nonlocal version
- virtual void computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J,
- const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfVtilde,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double* T = NULL
- ) const
+ void computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double* T = NULL
+ ) const
{
Msg::Fatal("mlawNonLocalPorousThomasonLawWithLodeEnhanced cannot be used");
};
virtual void computeDResidual_multipleNonLocalVariables(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- std::vector<double> &dres0dtildefV, std::vector<double> &dres1dtildefV, std::vector<double> &dres2dtildefV,
+ std::vector<double> &dres0dNonlocalVars, std::vector<double> &dres1dNonlocalVars, std::vector<double> &dres2dNonlocalVars,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const std::vector<double>& DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDkcorpr, const std::vector<double>& DyieldfVDtildefV,
+ const STensor43& DkcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const STensor3* DKcorprVDT = NULL,
double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
) const
{
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp
index bfe96ff2dcc0a7355ce5de6424dd9d16084570c2..a3e06082fedb20544914b34d98cec38ebe907a3d 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp
@@ -188,11 +188,6 @@ void mlawNonLocalPorousCoupledLaw::setFailureTolerance(const double NewFailureTo
_mlawCoales->setFailureTolerance(NewFailureTol);
};
-void mlawNonLocalPorousCoupledLaw::setLocalRegularizedFunction(const scalarFunction& fct){
- mlawNonLocalPorosity::setLocalRegularizedFunction(fct);
- _mlawGrowth->setLocalRegularizedFunction(fct);
- _mlawCoales->setLocalRegularizedFunction(fct);
-};
void mlawNonLocalPorousCoupledLaw::setCorrectedRegularizedFunction(const scalarFunction& fct){
mlawNonLocalPorosity::setCorrectedRegularizedFunction(fct);
_mlawGrowth->setCorrectedRegularizedFunction(fct);
@@ -320,72 +315,93 @@ double mlawNonLocalPorousCoupledLaw::yieldFunction(const double kcorEq, const do
else return _mlawGrowth->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
};
-
-void mlawNonLocalPorousCoupledLaw::computeResidual(SVector3& res, STensor3& J,
+void mlawNonLocalPorousCoupledLaw::computeResidual_NonLocalPorosity(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
+ const double yieldfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
const double* T
) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
- _mlawCoales->computeResidual(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
- DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,
- DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
+ _mlawCoales->computeResidual_NonLocalPorosity(res,J,kcorEqpr,pcorpr,R,H,yieldfV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,T);
}
else {
- _mlawGrowth->computeResidual(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
- DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,
- DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
+ _mlawGrowth->computeResidual_NonLocalPorosity(res,J,kcorEqpr,pcorpr,R,H,yieldfV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,T);
}
}
-void mlawNonLocalPorousCoupledLaw::computeDResidual(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
+
+void mlawNonLocalPorousCoupledLaw::computeDResidual_NonLocalPorosity(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const double DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDKcorpr, const double DyieldfVDtildefV,
+ const STensor43& DKcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T, const double* DRDT,
- const STensor3* DdevKcorprDT , const double* DpcorprDT, const double* DDeltafVDT,
+ const STensor3* DdevKcorprDT , const double* DpcorprDT, const STensor3* DKcorprDT,
double *dres0dT, double *dres1dT, double *dres2dT
) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
- _mlawCoales->computeDResidual(dres0dEpr,dres1dEpr,dres2dEpr,
+ _mlawCoales->computeDResidual_NonLocalPorosity(dres0dEpr,dres1dEpr,dres2dEpr,
dres0dtildefV,dres1dtildefV,dres2dtildefV,
devKcorpr,kcorEqpr,pcorpr,R,
- yieldfV,DyieldfVDfV,DyieldfVDtildefV,
- DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
+ yieldfV,DyieldfVDKcorpr,DyieldfVDtildefV,
+ DKcorprDEpr,DdevKcorprDEpr,DpcorprDEpr,
q0, q1,
DeltaHatD,DeltaHatQ,DeltaHatP,
T,DRDT,
- DdevKcorprDT ,DpcorprDT,DDeltafVDT,
+ DdevKcorprDT ,DpcorprDT,DKcorprDT,
dres0dT,dres1dT,dres2dT);
}
else {
- _mlawGrowth->computeDResidual(dres0dEpr,dres1dEpr,dres2dEpr,
+ _mlawGrowth->computeDResidual_NonLocalPorosity(dres0dEpr,dres1dEpr,dres2dEpr,
dres0dtildefV,dres1dtildefV,dres2dtildefV,
devKcorpr,kcorEqpr,pcorpr,R,
- yieldfV,DyieldfVDfV,DyieldfVDtildefV,
- DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
+ yieldfV,DyieldfVDKcorpr,DyieldfVDtildefV,
+ DKcorprDEpr,DdevKcorprDEpr,DpcorprDEpr,
q0, q1,
DeltaHatD,DeltaHatQ,DeltaHatP,
T,DRDT,
- DdevKcorprDT ,DpcorprDT,DDeltafVDT,
+ DdevKcorprDT ,DpcorprDT,DKcorprDT,
dres0dT,dres1dT,dres2dT);
}
};
+void mlawNonLocalPorousCoupledLaw::computeResidualLocal(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV, const double DyieldfVDfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
+ const double* T
+ ) const{
+
+ bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
+
+ if (coalesFlag){
+ _mlawCoales->computeResidualLocal(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,
+ DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
+ }
+ else {
+ _mlawGrowth->computeResidualLocal(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,
+ DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
+ }
+}
+
+
void mlawNonLocalPorousCoupledLaw::computeDResidualLocal(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
const double yieldfV, const double DyieldfVDfV,
@@ -753,56 +769,56 @@ void mlawNonLocalPorousCoupledLaw::checkFailed(IPNonLocalPorosity* q1, const IPN
void mlawNonLocalPorousCoupledLaw::computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
+ const double yieldfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double* T
+ const double* T
) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
- _mlawCoales->computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
- DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
+ _mlawCoales->computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,R,H,yieldfV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,T);
}
else {
- _mlawGrowth->computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
- DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
+ _mlawGrowth->computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,R,H,yieldfV,q0,q1,
+ DeltaHatD,DeltaHatQ,DeltaHatP,T);
}
};
void mlawNonLocalPorousCoupledLaw::computeDResidual_multipleNonLocalVariables(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- std::vector<double> &dres0dtildefV, std::vector<double> &dres1dtildefV, std::vector<double> &dres2dtildefV,
- const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const std::vector<double>& DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
- const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
- const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
- const double* T, const double* DRDT,
- const STensor3* DdevKcorprDT, const double* DpcorprDT, const double* DDeltafVDT,
- double *dres0dT, double *dres1dT, double *dres2dT
+ std::vector<double> &dres0dNonlocalVars, std::vector<double> &dres1dNonlocalVars, std::vector<double> &dres2dNonlocalVars,
+ const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
+ const double yieldfV, const STensor3& DyieldfVDkcorpr, const std::vector<double>& DyieldfVDNonlocalVars,
+ const STensor43& DkcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
+ const double* T, const double* DRDT,
+ const STensor3* DdevKcorprDT, const double* DpcorprDT, const STensor3* DKcorprVDT,
+ double *dres0dT, double *dres1dT, double *dres2dT
) const{
bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
if (coalesFlag){
_mlawCoales->computeDResidual_multipleNonLocalVariables(dres0dEpr,dres1dEpr,dres2dEpr,
- dres0dtildefV,dres1dtildefV,dres2dtildefV,
+ dres0dNonlocalVars,dres1dNonlocalVars,dres2dNonlocalVars,
devKcorpr,kcorEqpr,pcorpr,R,
- yieldfV,DyieldfVDfV,DyieldfVDtildefV,
- DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
+ yieldfV,DyieldfVDkcorpr,DyieldfVDNonlocalVars,
+ DkcorprDEpr,DdevKcorprDEpr,DpcorprDEpr,
q0, q1,
DeltaHatD,DeltaHatQ,DeltaHatP,
T,DRDT,
- DdevKcorprDT ,DpcorprDT,DDeltafVDT,
+ DdevKcorprDT ,DpcorprDT,DKcorprVDT,
dres0dT,dres1dT,dres2dT);
}
else {
_mlawGrowth->computeDResidual_multipleNonLocalVariables(dres0dEpr,dres1dEpr,dres2dEpr,
- dres0dtildefV,dres1dtildefV,dres2dtildefV,
+ dres0dNonlocalVars,dres1dNonlocalVars,dres2dNonlocalVars,
devKcorpr,kcorEqpr,pcorpr,R,
- yieldfV,DyieldfVDfV,DyieldfVDtildefV,
- DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
+ yieldfV,DyieldfVDkcorpr,DyieldfVDNonlocalVars,
+ DkcorprDEpr,DdevKcorprDEpr,DpcorprDEpr,
q0, q1,
DeltaHatD,DeltaHatQ,DeltaHatP,
T,DRDT,
- DdevKcorprDT ,DpcorprDT,DDeltafVDT,
+ DdevKcorprDT ,DpcorprDT,DKcorprVDT,
dres0dT,dres1dT,dres2dT);
}
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.h b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.h
index 374dca8f9525213922fd58d628fdfdc47616024e..c55d6bdff6b15ba357a2c171250faa4a686be806 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.h
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.h
@@ -58,7 +58,6 @@ class mlawNonLocalPorousCoupledLaw : public mlawNonLocalPorosity
virtual void setLodeOffset(const bool fl);
virtual void setShearPorosityGrowthFactor(const double k);
- virtual void setLocalRegularizedFunction(const scalarFunction& fct);
virtual void setCorrectedRegularizedFunction(const scalarFunction& fct);
virtual void setFailureTolerance(const double NewFailureTol);
virtual void setPostBlockedDissipationBehavior(const int method);
@@ -112,27 +111,34 @@ class mlawNonLocalPorousCoupledLaw : public mlawNonLocalPorosity
virtual double yieldFunction(const double kcorEq, const double pcor, const double R, const double yieldfV,
const IPNonLocalPorosity* q0, const IPNonLocalPorosity* q1, const double* T = NULL) const;
- virtual void computeResidual(SVector3& res, STensor3& J,
+ virtual void computeResidual_NonLocalPorosity(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
+ const double yieldfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
const double* T = NULL
) const;
- virtual void computeDResidual(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
+ virtual void computeDResidual_NonLocalPorosity(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const double DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDKcorpr , const double DyieldfVDtildefV,
+ const STensor43& DKcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const STensor3* DKcorprDT = NULL,
double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
) const;
-
+
+ virtual void computeResidualLocal(SVector3& res, STensor3& J,
+ const double kcorEqpr, const double pcorpr, const double R, const double H,
+ const double yieldfV, const double DyieldfVDfV,
+ const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
+ const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
+ const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
+ const double* T = NULL
+ ) const;
virtual void computeDResidualLocal(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
const double yieldfV, const double DyieldfVDfV,
@@ -147,21 +153,21 @@ class mlawNonLocalPorousCoupledLaw : public mlawNonLocalPorosity
virtual void computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J,
const double kcorEqpr, const double pcorpr, const double R, const double H,
- const double yieldfV, const double DyieldfVDfV,
+ const double yieldfV,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
- const double DeltafV, const double* T = NULL
+ const double* T = NULL
) const;
virtual void computeDResidual_multipleNonLocalVariables(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
- std::vector<double> &dres0dtildefV, std::vector<double> &dres1dtildefV, std::vector<double> &dres2dtildefV,
+ std::vector<double> &dres0dNonlocalVars, std::vector<double> &dres1dNonlocalVars, std::vector<double> &dres2dNonlocalVars,
const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
- const double yieldfV, const double DyieldfVDfV, const std::vector<double>& DyieldfVDtildefV,
- const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
+ const double yieldfV, const STensor3& DyieldfVDkcorpr, const std::vector<double>& DyieldfVDNonlocalVars,
+ const STensor43& DkcorprDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
const double* T = NULL, const double* DRDT= NULL,
- const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const double* DDeltafVDT = NULL,
+ const STensor3* DdevKcorprDT = NULL, const double* DpcorprDT = NULL, const STensor3* DKcorprVDT = NULL,
double *dres0dT= NULL, double *dres1dT= NULL, double *dres2dT = NULL
) const;
diff --git a/NonLinearSolver/materialLaw/voidStateEvolutionLaw.cpp b/NonLinearSolver/materialLaw/voidStateEvolutionLaw.cpp
index e6a364508f89219555c47aa37603eb2704ffdf36..4f2ac57efa7fc9ae38ac89b5923851f25698fd50 100644
--- a/NonLinearSolver/materialLaw/voidStateEvolutionLaw.cpp
+++ b/NonLinearSolver/materialLaw/voidStateEvolutionLaw.cpp
@@ -58,10 +58,12 @@ void SphericalVoidStateEvolutionLawWithMatrixPlasticDeformation::evolve(const do
const IPVoidState& sPrev, IPVoidState& sCur) const{
double lambdaPrev = sPrev.getVoidSpacingRatio();
+ double roughChiPrev = sPrev.getRoughVoidLigamentRatio();
double chiPrev = sPrev.getVoidLigamentRatio();
double& lambda = sCur.getRefToVoidSpacingRatio();
// Update Chi and lambda from previous state and then, increment
+ double& roughChi = sCur.getRefToRoughVoidLigamentRatio();
double& Chi = sCur.getRefToVoidLigamentRatio();
// new lambda
double DlambdaDHatP = 0;
@@ -72,11 +74,8 @@ void SphericalVoidStateEvolutionLawWithMatrixPlasticDeformation::evolve(const do
}
// new regularised chi
double fact = yieldfV1*lambda/(yieldfV0*lambdaPrev);
- double ChiComputed = chiPrev*pow(fact,1./3.);
-
- double ffPrev = _ChiRegularizedFunction->inverse(chiPrev);
- double ff = ffPrev+ChiComputed-chiPrev;
- Chi = _ChiRegularizedFunction->getVal(ff);
+ roughChi = roughChiPrev*pow(fact,1./3.);
+ Chi = _ChiRegularizedFunction->getVal(roughChi);
double& dChidf = sCur.getRefToDVoidLigamentRatioDYieldPorosity();
@@ -84,12 +83,12 @@ void SphericalVoidStateEvolutionLawWithMatrixPlasticDeformation::evolve(const do
double& dChidhatD = sCur.getRefToDVoidLigamentRatioDDeviatoricPlasticDeformation();
double& dChidhatP = sCur.getRefToDVoidLigamentRatioDMatrixPlasticDeformation();
- double DiffChiComputed = _ChiRegularizedFunction->getDiff(ff);
+ double DiffRoughChi = _ChiRegularizedFunction->getDiff(roughChi);
// new Chi rate
- dChidf = DiffChiComputed*ChiComputed/(3.*yieldfV1);
+ dChidf = DiffRoughChi*roughChi/(3.*yieldfV1);
dChidhatQ = 0.;
dChidhatD = 0.;
- dChidhatP = DiffChiComputed*DlambdaDHatP*ChiComputed/(3.*lambda);
+ dChidhatP = DiffRoughChi*DlambdaDHatP*roughChi/(3.*lambda);
};
@@ -105,10 +104,12 @@ void SphericalVoidStateEvolutionLawWithDeviatoricPlasticDeformation::evolve(cons
const IPVoidState& sPrev, IPVoidState& sCur) const{
double lambdaPrev = sPrev.getVoidSpacingRatio();
+ double roughChiPrev = sPrev.getRoughVoidLigamentRatio();
double chiPrev = sPrev.getVoidLigamentRatio();
double& lambda = sCur.getRefToVoidSpacingRatio();
// Update Chi and lambda from previous state and then, increment
+ double& roughChi = sCur.getRefToRoughVoidLigamentRatio();
double& Chi = sCur.getRefToVoidLigamentRatio();
// new lambda
double DlambdaDHatD = 0;
@@ -119,21 +120,19 @@ void SphericalVoidStateEvolutionLawWithDeviatoricPlasticDeformation::evolve(cons
}
// new regularised chi
double fact = yieldfV1*lambda/(yieldfV0*lambdaPrev);
- double ChiComputed = chiPrev*pow(fact,1./3.);
- double ffPrev = _ChiRegularizedFunction->inverse(chiPrev);
- double ff = ffPrev+ChiComputed-chiPrev;
- Chi = _ChiRegularizedFunction->getVal(ff);
+ roughChi = roughChiPrev*pow(fact,1./3.);
+ Chi = _ChiRegularizedFunction->getVal(roughChi);
double& dChidf = sCur.getRefToDVoidLigamentRatioDYieldPorosity();
double& dChidhatQ = sCur.getRefToDVoidLigamentRatioDVolumetricPlasticDeformation();
double& dChidhatD = sCur.getRefToDVoidLigamentRatioDDeviatoricPlasticDeformation();
double& dChidhatP = sCur.getRefToDVoidLigamentRatioDMatrixPlasticDeformation();
- double DiffChiComputed = _ChiRegularizedFunction->getDiff(ff);
+ double DiffChiComputed = _ChiRegularizedFunction->getDiff(roughChi);
// new Chi rate
- dChidf = DiffChiComputed*ChiComputed/(3.*yieldfV1);
+ dChidf = DiffChiComputed*roughChi/(3.*yieldfV1);
dChidhatQ = 0.;
- dChidhatD = DiffChiComputed*DlambdaDHatD*ChiComputed/(3.*lambda);
+ dChidhatD = DiffChiComputed*DlambdaDHatD*roughChi/(3.*lambda);
dChidhatP = 0.;
};
@@ -181,7 +180,8 @@ void SpheroidalVoidStateEvolutionLawWithAspectRatioEvolution::evolve(const doubl
double& dWdHatQ = sCur.getRefToDVoidAspectRatioDVolumetricPlasticDeformation();
double& dWdHatD = sCur.getRefToDVoidAspectRatioDDeviatoricPlasticDeformation();
double& dWdHatP = sCur.getRefToDVoidAspectRatioDMatrixPlasticDeformation();
-
+
+ double& roughX = sCur.getRefToRoughVoidLigamentRatio();
double& X = sCur.getRefToVoidLigamentRatio();
double& dXdf = sCur.getRefToDVoidLigamentRatioDYieldPorosity();
double& dXdHatQ = sCur.getRefToDVoidLigamentRatioDVolumetricPlasticDeformation();
@@ -201,6 +201,7 @@ void SpheroidalVoidStateEvolutionLawWithAspectRatioEvolution::evolve(const doubl
dgammadHatP = 0.;
// previpis value
+ double roughXPrev = sPrev.getRoughVoidLigamentRatio();
double lambdaPrev = sPrev.getVoidSpacingRatio();
double XPrev = sPrev.getVoidLigamentRatio();
double WPrev = sPrev.getVoidAspectRatio();
@@ -224,46 +225,42 @@ void SpheroidalVoidStateEvolutionLawWithAspectRatioEvolution::evolve(const doubl
// new regularised chi
double fact = yieldfV1*lambda/(yieldfV0*lambdaPrev);
- double XComputed = XPrev*pow(fact,1./3.);
- double ffPrev = _ChiRegularizedFunction->inverse(XPrev);
- double ff = ffPrev+XComputed-XPrev;
- X = _ChiRegularizedFunction->getVal(ff);
- double DXComputed = _ChiRegularizedFunction->getDiff(ff);
+ roughX = roughXPrev*pow(fact,1./3.);
+ X = _ChiRegularizedFunction->getVal(roughX);
+ double DXComputed = _ChiRegularizedFunction->getDiff(roughX);
// new Chi rate
- dXdf = DXComputed*XComputed/(3.*yieldfV1);
+ dXdf = DXComputed*roughX/(3.*yieldfV1);
dXdHatQ = 0.;
- dXdHatD = DXComputed*DlambdaDhatD*XComputed/(3.*lambda);
+ dXdHatD = DXComputed*DlambdaDhatD*roughX/(3.*lambda);
dXdHatP = 0.;
}
else{
if (DeltahatD > 0){
- double Zprev = XPrev*WPrev;
+ double Zprev = roughXPrev*WPrev;
double DeltaZ = 1.5*(lambda-lambdaPrev)/_kappa;
double Z = Zprev+DeltaZ;
- double Xcomputed = sqrt(1.5-(1.5-XPrev*XPrev)*Zprev/Z);
- double dXcomputeddhatD = (3./4.)*(lambda*Xcomputed/Z)*(1.5/Xcomputed/Xcomputed-1.);
-
- double ffPrev = _ChiRegularizedFunction->inverse(XPrev);
+ roughX = sqrt(1.5-(1.5-roughXPrev*roughXPrev)*Zprev/Z);
+ double droughXdhatD = (3./4.)*(lambda*roughX/Z)*(1.5/roughX/roughX-1.);
- double ff = ffPrev+Xcomputed-XPrev;
- X = _ChiRegularizedFunction->getVal(ff);
- double DXDXcomputed = _ChiRegularizedFunction->getDiff(ff);
+ X = _ChiRegularizedFunction->getVal(roughX);
+ double DXDroughX = _ChiRegularizedFunction->getDiff(roughX);
dXdf =0.;
dXdHatQ = 0.;
- dXdHatD = DXDXcomputed*dXcomputeddhatD;
+ dXdHatD = DXDroughX*droughXdhatD;
dXdHatP = 0.;
- W = Z/Xcomputed;
+ W = Z/roughX;
dWdf = 0.;
dWdHatQ = 0.;
- dWdHatD = (9./4.)*(lambda/Xcomputed)*(1.-0.5/Xcomputed*Xcomputed);
+ dWdHatD = (9./4.)*(lambda/roughX)*(1.-0.5/roughX*roughX);
dWdHatP = 0.;
}
else{
+ roughX = roughXPrev;
X = XPrev;
dXdf =0.;
dXdHatQ = 0.;
diff --git a/NonLinearSolver/nlsolver/scalarFunction.h b/NonLinearSolver/nlsolver/scalarFunction.h
index 9ccc85f5a3bbf0197fcedc43d383230d48dfef50..faae3957749e2d9152fa91422aff7ddcd0ea4959 100644
--- a/NonLinearSolver/nlsolver/scalarFunction.h
+++ b/NonLinearSolver/nlsolver/scalarFunction.h
@@ -25,7 +25,6 @@ class scalarFunction{
virtual ~scalarFunction(){};
virtual void setElement(MElement* el) const{}
virtual double getVal(const double x) const = 0; // function value at x
- virtual double inverse(const double y) const = 0; // resolve function y = f(x) to get x
virtual double getDiff(const double x) const = 0; // function derivative at x
virtual double getVal(const std::vector<double>& input) const {
return getVal(input[0]);
@@ -59,11 +58,6 @@ class constantScalarFunction : public scalarFunction{
virtual ~constantScalarFunction(){};
virtual double getVal(const double x) const {return _f;};
- virtual double inverse(const double y) const {
- // resolve function y = f(x) to get x
- Msg::Error("there is infinite solution in constantScalarFunction::inverse");
- return 0.;
- };
virtual double getDiff(const double x) const {return 0.;};
virtual scalarFunction* clone() const {return new constantScalarFunction(*this);};
#endif // SWIG
@@ -96,13 +90,7 @@ class elementWiseScalarFunction : public constantScalarFunction{
std::map<int,double>::const_iterator it = _constantOnElements.find(_curELe->getNum());
if (it != _constantOnElements.end()) return it->second;
return _f;
- };
- virtual double inverse(const double y) const {
- // resolve function y = f(x) to get x
- Msg::Error("there is infinite solution in elementWiseScalarFunction::inverse");
- return 0.;
- };
-
+ };
virtual double getDiff(const double x) const {
Msg::Error("elementWiseScalarFunction::getDiff should not be used");
return 0.;
@@ -136,11 +124,6 @@ class twoValuesFunction : public scalarFunction{
if (x < _x01 and x > _x00) return _A0;
else return _A1;
};
- virtual double inverse(const double y) const {
- // resolve function y = f(x) to get x
- Msg::Error("there is infinite solution in twoValuesFunction::inverse");
- return 0.;
- };
virtual double getDiff(const double x) const {return 0.;};
virtual scalarFunction* clone() const {return new twoValuesFunction(*this);};
#endif // SWIG
@@ -167,13 +150,6 @@ class linearScalarFunction : public scalarFunction{
virtual ~linearScalarFunction(){};
virtual double getVal(const double x) const {return _a*(x-_x0)+_f0;};
- virtual double inverse(const double y) const {
- if (_a == 0.){
- Msg::Error("twoValuesFunction::inverse cannot be used");
- return 0.;
- }
- return _x0+(y-_f0)/_a;
- };
virtual double getDiff(const double x) const {return _a;};
virtual scalarFunction* clone() const {return new linearScalarFunction(*this);};
#endif // SWIG
@@ -217,9 +193,6 @@ protected:
}
};
- virtual double inverse(const double y) const {
- Msg::Fatal("piecewiseScalarFunction::inverse has not been implemented yet");
- };
virtual double getDiff(const double x) const {
int n = _var.size();
if (x < _var[0]) return 0.;
@@ -276,24 +249,6 @@ class twoVariableExponentialSaturationScalarFunction : public scalarFunction{
}
}
};
- virtual double inverse(const double y) const {
- if(_upperBound){
- if (y >= _f2) Msg::Fatal("no solution in twoVariableExponentialSaturationScalarFunction::inverse upperBound");
- if (y < _f1) return y;
- else{
- double fact = (_f2-y)/(_f2-_f1);
- return _f1 - (_f2-_f1)*log(fact);
- }
- }
- else{
- if (y <= _f2) Msg::Fatal("no solution in twoVariableExponentialSaturationScalarFunction::inverse lowerBound");
- if (y > _f1) return y;
- else{
- double fact = (_f2-y)/(_f2-_f1);
- return _f1 - (_f2-_f1)*log(fact);
- }
- }
- }
virtual double getDiff(const double x) const {
if (_upperBound){
if (x < _f1)
@@ -370,20 +325,6 @@ protected:
}
};
- virtual double inverse(const double y) const {
- if ( y >= _f2max or y <= _f2min) Msg::Fatal("fourVariableExponentialSaturationScalarFunction::inverse has no solution with val = %e",y);
- if (y < _f1min){
- double fact = (_f2min-y)/(_f2min-_f1min);
- return _f1min - (_f2min-_f1min)*log(fact);
- }
- else if (y > _f1max){
- double fact = (_f2max-y)/(_f2max-_f1max);
- return _f1max - (_f2max-_f1max)*log(fact);
- }
- else{
- return y;
- }
- };
virtual double getDiff(const double x) const {
if (x < _f1min)
@@ -443,9 +384,6 @@ class polynomialScalarFunction : public scalarFunction{
else if (_order == 3) return _coeffs[0]+_coeffs[1]*x+_coeffs[2]*x*x+_coeffs[3]*x*x*x;
else Msg::Fatal("polynomialScalarFunction with order %d has not been implemented");
};
- virtual double inverse(const double y) const {
- Msg::Fatal("polynomialScalarFunction::inverse is not implemented!!!");
- };
virtual double getDiff(const double x) const {
if (_order == 0) return 0.;
else if (_order == 1) return _coeffs[1];
@@ -477,10 +415,6 @@ class asinwt : public scalarFunction{
virtual ~asinwt(){};
virtual double getVal(const double x) const {return _a*sin(_omega*x);};
- virtual double inverse(const double y) const{
- if (fabs(y) > fabs(_a)) Msg::Fatal("asinwt::inverse has not solution");
- else return (asin(y/_a))/_omega;
- }
virtual double getDiff(const double x) const {return _a*_omega*cos(_omega*x);};
virtual scalarFunction* clone() const {return new asinwt(*this);};
#endif // SWIG
diff --git a/dG3D/benchmarks/CMakeLists.txt b/dG3D/benchmarks/CMakeLists.txt
index 40952cdaa489cc6d43eb031cfb828c0aa945e271..4fdfc5eaa531d289972d83154bad476857fd2bfd 100644
--- a/dG3D/benchmarks/CMakeLists.txt
+++ b/dG3D/benchmarks/CMakeLists.txt
@@ -176,3 +176,4 @@ add_subdirectory(Thomason_cube_I1J2J3)
add_subdirectory(GursonThomason_cube_I1J2J3)
add_subdirectory(GursonThomason_cube_withLodeEnhanced)
add_subdirectory(Thomason_planeStrain_withLodeEnhanced)
+add_subdirectory(LocalGurson_PathFollowing)
diff --git a/dG3D/benchmarks/LocalGurson_PathFollowing/NodalDisplacement18comp1.csv b/dG3D/benchmarks/LocalGurson_PathFollowing/NodalDisplacement18comp1.csv
deleted file mode 100644
index 36f0a6347036841273c338576e68e81037c07767..0000000000000000000000000000000000000000
--- a/dG3D/benchmarks/LocalGurson_PathFollowing/NodalDisplacement18comp1.csv
+++ /dev/null
@@ -1,213 +0,0 @@
-3.333333e-03;5.231572e-04
-6.666667e-03;-3.608092e-03
-1.000000e-02;-9.666361e-03
-1.333333e-02;-1.596241e-02
-1.666667e-02;-2.235176e-02
-2.000000e-02;-2.873383e-02
-2.333333e-02;-3.502402e-02
-2.666667e-02;-4.336159e-02
-3.000000e-02;-5.166774e-02
-3.333333e-02;-6.270994e-02
-3.666667e-02;-7.738285e-02
-4.000000e-02;-9.689127e-02
-4.333333e-02;-1.228686e-01
-4.666667e-02;-1.575335e-01
-5.000000e-02;-2.040408e-01
-5.333333e-02;-2.588823e-01
-5.666667e-02;-3.144424e-01
-6.000000e-02;-3.708322e-01
-6.333333e-02;-4.280457e-01
-6.666667e-02;-4.744352e-01
-7.000000e-02;-5.215554e-01
-7.333333e-02;-5.695508e-01
-7.666667e-02;-6.185366e-01
-8.000000e-02;-6.813001e-01
-8.333333e-02;-7.458442e-01
-8.666667e-02;-8.122454e-01
-9.000000e-02;-8.667188e-01
-9.333333e-02;-9.224292e-01
-9.666667e-02;-9.793877e-01
-1.000000e-01;-1.037639e+00
-1.033333e-01;-1.097254e+00
-1.066667e-01;-1.146043e+00
-1.100000e-01;-1.185856e+00
-1.133333e-01;-1.218368e+00
-1.166667e-01;-1.244935e+00
-1.200000e-01;-1.272028e+00
-1.233333e-01;-1.299722e+00
-1.266667e-01;-1.328087e+00
-1.300000e-01;-1.357252e+00
-1.333333e-01;-1.386753e+00
-1.366667e-01;-1.398605e+00
-1.400000e-01;-1.408076e+00
-1.433333e-01;-1.420623e+00
-1.466667e-01;-1.437200e+00
-1.500000e-01;-1.443835e+00
-1.533333e-01;-1.450458e+00
-1.566667e-01;-1.459252e+00
-1.600000e-01;-1.466274e+00
-1.633333e-01;-1.475567e+00
-1.666667e-01;-1.481728e+00
-1.700000e-01;-1.487864e+00
-1.733333e-01;-1.495993e+00
-1.766667e-01;-1.504129e+00
-1.800000e-01;-1.514934e+00
-1.833333e-01;-1.518235e+00
-1.866667e-01;-1.520856e+00
-1.900000e-01;-1.524308e+00
-1.933333e-01;-1.528824e+00
-1.966667e-01;-1.534671e+00
-2.000000e-01;-1.542136e+00
-2.033333e-01;-1.551469e+00
-2.066667e-01;-1.560065e+00
-2.100000e-01;-1.570364e+00
-2.133333e-01;-1.579410e+00
-2.166667e-01;-1.587281e+00
-2.200000e-01;-1.596165e+00
-2.233333e-01;-1.603320e+00
-2.266667e-01;-1.609001e+00
-2.300000e-01;-1.613379e+00
-2.333333e-01;-1.616687e+00
-2.366667e-01;-1.619111e+00
-2.400000e-01;-1.620807e+00
-2.433333e-01;-1.621847e+00
-2.466667e-01;-1.622240e+00
-2.500000e-01;-1.622109e+00
-2.533333e-01;-1.621336e+00
-2.566667e-01;-1.619946e+00
-2.600000e-01;-1.617432e+00
-2.633333e-01;-1.614212e+00
-2.666667e-01;-1.610391e+00
-2.700000e-01;-1.606272e+00
-2.716667e-01;-1.604967e+00
-2.733333e-01;-1.604307e+00
-2.750000e-01;-1.603443e+00
-2.766667e-01;-1.602343e+00
-2.783333e-01;-1.601089e+00
-2.800000e-01;-1.600345e+00
-2.816667e-01;-1.600868e+00
-2.818750e-01;-1.600909e+00
-2.820833e-01;-1.600944e+00
-2.822917e-01;-1.600993e+00
-2.825000e-01;-1.601059e+00
-2.827083e-01;-1.601145e+00
-2.829167e-01;-1.601253e+00
-2.831250e-01;-1.601380e+00
-2.833333e-01;-1.601511e+00
-2.835417e-01;-1.601580e+00
-2.837500e-01;-1.601452e+00
-2.839583e-01;-1.601377e+00
-2.841667e-01;-1.601324e+00
-2.843750e-01;-1.601266e+00
-2.845833e-01;-1.601206e+00
-2.847917e-01;-1.601143e+00
-2.850000e-01;-1.601074e+00
-2.852083e-01;-1.600996e+00
-2.854167e-01;-1.600907e+00
-2.856250e-01;-1.600805e+00
-2.858333e-01;-1.600748e+00
-2.860417e-01;-1.600681e+00
-2.862500e-01;-1.600604e+00
-2.864583e-01;-1.600560e+00
-2.866667e-01;-1.600506e+00
-2.868750e-01;-1.600468e+00
-2.870833e-01;-1.600421e+00
-2.872917e-01;-1.600391e+00
-2.875000e-01;-1.600369e+00
-2.877083e-01;-1.600346e+00
-2.879167e-01;-1.600326e+00
-2.881250e-01;-1.600303e+00
-2.883333e-01;-1.600294e+00
-2.885417e-01;-1.600284e+00
-2.887500e-01;-1.600274e+00
-2.889583e-01;-1.600268e+00
-2.891667e-01;-1.600262e+00
-2.893750e-01;-1.600259e+00
-2.895833e-01;-1.600258e+00
-2.897917e-01;-1.600259e+00
-2.900000e-01;-1.600259e+00
-2.902083e-01;-1.600261e+00
-2.904167e-01;-1.600264e+00
-2.906250e-01;-1.600269e+00
-2.908333e-01;-1.600275e+00
-2.910417e-01;-1.600285e+00
-2.912500e-01;-1.600298e+00
-2.914583e-01;-1.600310e+00
-2.916667e-01;-1.600328e+00
-2.918750e-01;-1.600333e+00
-2.920833e-01;-1.600340e+00
-2.922917e-01;-1.600350e+00
-2.925000e-01;-1.600361e+00
-2.927083e-01;-1.600375e+00
-2.929167e-01;-1.600394e+00
-2.931250e-01;-1.600402e+00
-2.933333e-01;-1.600410e+00
-2.935417e-01;-1.600418e+00
-2.937500e-01;-1.600427e+00
-2.939583e-01;-1.600440e+00
-2.941667e-01;-1.600460e+00
-2.943750e-01;-1.600485e+00
-2.945833e-01;-1.600512e+00
-2.947917e-01;-1.600521e+00
-2.950000e-01;-1.600532e+00
-2.952083e-01;-1.600543e+00
-2.954167e-01;-1.600553e+00
-2.956250e-01;-1.600563e+00
-2.958333e-01;-1.600570e+00
-2.960417e-01;-1.600575e+00
-2.962500e-01;-1.600577e+00
-2.964583e-01;-1.600579e+00
-2.965625e-01;-1.600580e+00
-2.966146e-01;-1.600580e+00
-2.966667e-01;-1.600580e+00
-2.967187e-01;-1.600580e+00
-2.967708e-01;-1.600580e+00
-2.968229e-01;-1.600580e+00
-2.968750e-01;-1.600580e+00
-2.969271e-01;-1.600580e+00
-2.969792e-01;-1.600580e+00
-2.970312e-01;-1.600580e+00
-2.970833e-01;-1.600580e+00
-2.971354e-01;-1.600581e+00
-2.971875e-01;-1.600581e+00
-2.972396e-01;-1.600581e+00
-2.972917e-01;-1.600581e+00
-2.973437e-01;-1.600581e+00
-2.973698e-01;-1.600581e+00
-2.973958e-01;-1.600581e+00
-2.974219e-01;-1.600581e+00
-2.974479e-01;-1.600581e+00
-2.974740e-01;-1.600581e+00
-2.975000e-01;-1.600581e+00
-2.975260e-01;-1.600581e+00
-2.975521e-01;-1.600581e+00
-2.975781e-01;-1.600581e+00
-2.976042e-01;-1.600581e+00
-2.976302e-01;-1.600581e+00
-2.976562e-01;-1.600581e+00
-2.976823e-01;-1.600581e+00
-2.977083e-01;-1.600581e+00
-2.977344e-01;-1.600581e+00
-2.977604e-01;-1.600581e+00
-2.977865e-01;-1.600581e+00
-2.978125e-01;-1.600581e+00
-2.978385e-01;-1.600581e+00
-2.978646e-01;-1.600581e+00
-2.978906e-01;-1.600581e+00
-2.979167e-01;-1.600581e+00
-2.979427e-01;-1.600581e+00
-2.979687e-01;-1.600581e+00
-2.979948e-01;-1.600581e+00
-2.980208e-01;-1.600581e+00
-2.980469e-01;-1.600581e+00
-2.980729e-01;-1.600581e+00
-2.980990e-01;-1.600581e+00
-2.981250e-01;-1.600581e+00
-2.981510e-01;-1.600581e+00
-2.981771e-01;-1.600581e+00
-2.982031e-01;-1.600581e+00
-2.982292e-01;-1.600581e+00
-2.982552e-01;-1.600581e+00
-2.982812e-01;-1.600581e+00
-2.983073e-01;-1.600581e+00
-2.983333e-01;-1.600581e+00
diff --git a/dG3D/benchmarks/LocalGurson_PathFollowing/NodalDisplacement19comp1.csv b/dG3D/benchmarks/LocalGurson_PathFollowing/NodalDisplacement19comp1.csv
deleted file mode 100644
index 93219e10b5a0c83afbe682839b12e1505998a9dd..0000000000000000000000000000000000000000
--- a/dG3D/benchmarks/LocalGurson_PathFollowing/NodalDisplacement19comp1.csv
+++ /dev/null
@@ -1,213 +0,0 @@
-3.333333e-03;2.322882e-03
-6.666667e-03;1.005008e-02
-1.000000e-02;1.655122e-02
-1.333333e-02;2.314741e-02
-1.666667e-02;2.976523e-02
-2.000000e-02;3.637185e-02
-2.333333e-02;4.293414e-02
-2.666667e-02;5.167398e-02
-3.000000e-02;6.041818e-02
-3.333333e-02;7.208282e-02
-3.666667e-02;8.763348e-02
-4.000000e-02;1.083510e-01
-4.333333e-02;1.359063e-01
-4.666667e-02;1.725148e-01
-5.000000e-02;2.211113e-01
-5.333333e-02;2.776734e-01
-5.666667e-02;3.341823e-01
-6.000000e-02;3.908783e-01
-6.333333e-02;4.480852e-01
-6.666667e-02;4.944498e-01
-7.000000e-02;5.415301e-01
-7.333333e-02;5.894724e-01
-7.666667e-02;6.383926e-01
-8.000000e-02;7.010575e-01
-8.333333e-02;7.654848e-01
-8.666667e-02;8.317506e-01
-9.000000e-02;8.861017e-01
-9.333333e-02;9.416763e-01
-9.666667e-02;9.984845e-01
-1.000000e-01;1.056569e+00
-1.033333e-01;1.115998e+00
-1.066667e-01;1.164623e+00
-1.100000e-01;1.204293e+00
-1.133333e-01;1.236672e+00
-1.166667e-01;1.263109e+00
-1.200000e-01;1.290063e+00
-1.233333e-01;1.317605e+00
-1.266667e-01;1.345794e+00
-1.300000e-01;1.374713e+00
-1.333333e-01;1.403895e+00
-1.366667e-01;1.415615e+00
-1.400000e-01;1.424972e+00
-1.433333e-01;1.437357e+00
-1.466667e-01;1.453698e+00
-1.500000e-01;1.460238e+00
-1.533333e-01;1.466763e+00
-1.566667e-01;1.475421e+00
-1.600000e-01;1.482331e+00
-1.633333e-01;1.491469e+00
-1.666667e-01;1.497527e+00
-1.700000e-01;1.503554e+00
-1.733333e-01;1.511538e+00
-1.766667e-01;1.519530e+00
-1.800000e-01;1.530164e+00
-1.833333e-01;1.533417e+00
-1.866667e-01;1.536004e+00
-1.900000e-01;1.539415e+00
-1.933333e-01;1.543883e+00
-1.966667e-01;1.549681e+00
-2.000000e-01;1.557099e+00
-2.033333e-01;1.566395e+00
-2.066667e-01;1.574968e+00
-2.100000e-01;1.585246e+00
-2.133333e-01;1.594274e+00
-2.166667e-01;1.602128e+00
-2.200000e-01;1.610979e+00
-2.233333e-01;1.618090e+00
-2.266667e-01;1.623722e+00
-2.300000e-01;1.628046e+00
-2.333333e-01;1.631296e+00
-2.366667e-01;1.633661e+00
-2.400000e-01;1.635296e+00
-2.433333e-01;1.636273e+00
-2.466667e-01;1.636601e+00
-2.500000e-01;1.636406e+00
-2.533333e-01;1.635549e+00
-2.566667e-01;1.634079e+00
-2.600000e-01;1.631460e+00
-2.633333e-01;1.628128e+00
-2.666667e-01;1.624277e+00
-2.700000e-01;1.619805e+00
-2.716667e-01;1.618343e+00
-2.733333e-01;1.617587e+00
-2.750000e-01;1.616550e+00
-2.766667e-01;1.615089e+00
-2.783333e-01;1.612840e+00
-2.800000e-01;1.609783e+00
-2.816667e-01;1.604521e+00
-2.818750e-01;1.604128e+00
-2.820833e-01;1.603809e+00
-2.822917e-01;1.603380e+00
-2.825000e-01;1.602804e+00
-2.827083e-01;1.602033e+00
-2.829167e-01;1.601000e+00
-2.831250e-01;1.599615e+00
-2.833333e-01;1.597745e+00
-2.835417e-01;1.595184e+00
-2.837500e-01;1.591546e+00
-2.839583e-01;1.589417e+00
-2.841667e-01;1.587697e+00
-2.843750e-01;1.585392e+00
-2.845833e-01;1.582312e+00
-2.847917e-01;1.578205e+00
-2.850000e-01;1.572756e+00
-2.852083e-01;1.565551e+00
-2.854167e-01;1.556064e+00
-2.856250e-01;1.543576e+00
-2.858333e-01;1.535256e+00
-2.860417e-01;1.524303e+00
-2.862500e-01;1.509881e+00
-2.864583e-01;1.500271e+00
-2.866667e-01;1.487624e+00
-2.868750e-01;1.477526e+00
-2.870833e-01;1.464249e+00
-2.872917e-01;1.453651e+00
-2.875000e-01;1.443116e+00
-2.877083e-01;1.429291e+00
-2.879167e-01;1.415617e+00
-2.881250e-01;1.397752e+00
-2.883333e-01;1.388791e+00
-2.885417e-01;1.377054e+00
-2.887500e-01;1.361806e+00
-2.889583e-01;1.349748e+00
-2.891667e-01;1.334165e+00
-2.893750e-01;1.319002e+00
-2.895833e-01;1.299571e+00
-2.897917e-01;1.291004e+00
-2.900000e-01;1.282588e+00
-2.902083e-01;1.271708e+00
-2.904167e-01;1.261107e+00
-2.906250e-01;1.247558e+00
-2.908333e-01;1.234513e+00
-2.910417e-01;1.218566e+00
-2.912500e-01;1.204233e+00
-2.914583e-01;1.193852e+00
-2.916667e-01;1.182211e+00
-2.918750e-01;1.178989e+00
-2.920833e-01;1.175052e+00
-2.922917e-01;1.170444e+00
-2.925000e-01;1.166435e+00
-2.927083e-01;1.162202e+00
-2.929167e-01;1.158372e+00
-2.931250e-01;1.156904e+00
-2.933333e-01;1.155734e+00
-2.935417e-01;1.154859e+00
-2.937500e-01;1.154306e+00
-2.939583e-01;1.153958e+00
-2.941667e-01;1.153807e+00
-2.943750e-01;1.153508e+00
-2.945833e-01;1.152336e+00
-2.947917e-01;1.151581e+00
-2.950000e-01;1.150327e+00
-2.952083e-01;1.148340e+00
-2.954167e-01;1.146025e+00
-2.956250e-01;1.142607e+00
-2.958333e-01;1.138813e+00
-2.960417e-01;1.135511e+00
-2.962500e-01;1.133174e+00
-2.964583e-01;1.130787e+00
-2.965625e-01;1.129956e+00
-2.966146e-01;1.129770e+00
-2.966667e-01;1.129583e+00
-2.967187e-01;1.129396e+00
-2.967708e-01;1.129209e+00
-2.968229e-01;1.129022e+00
-2.968750e-01;1.128835e+00
-2.969271e-01;1.128648e+00
-2.969792e-01;1.128461e+00
-2.970312e-01;1.128273e+00
-2.970833e-01;1.128086e+00
-2.971354e-01;1.127898e+00
-2.971875e-01;1.127710e+00
-2.972396e-01;1.127460e+00
-2.972917e-01;1.127126e+00
-2.973437e-01;1.126992e+00
-2.973698e-01;1.126965e+00
-2.973958e-01;1.126944e+00
-2.974219e-01;1.126923e+00
-2.974479e-01;1.126901e+00
-2.974740e-01;1.126880e+00
-2.975000e-01;1.126859e+00
-2.975260e-01;1.126838e+00
-2.975521e-01;1.126821e+00
-2.975781e-01;1.126804e+00
-2.976042e-01;1.126781e+00
-2.976302e-01;1.126763e+00
-2.976562e-01;1.126749e+00
-2.976823e-01;1.126729e+00
-2.977083e-01;1.126703e+00
-2.977344e-01;1.126686e+00
-2.977604e-01;1.126673e+00
-2.977865e-01;1.126659e+00
-2.978125e-01;1.126645e+00
-2.978385e-01;1.126632e+00
-2.978646e-01;1.126618e+00
-2.978906e-01;1.126605e+00
-2.979167e-01;1.126591e+00
-2.979427e-01;1.126580e+00
-2.979687e-01;1.126565e+00
-2.979948e-01;1.126551e+00
-2.980208e-01;1.126531e+00
-2.980469e-01;1.126505e+00
-2.980729e-01;1.126488e+00
-2.980990e-01;1.126474e+00
-2.981250e-01;1.126460e+00
-2.981510e-01;1.126447e+00
-2.981771e-01;1.126433e+00
-2.982031e-01;1.126419e+00
-2.982292e-01;1.126405e+00
-2.982552e-01;1.126391e+00
-2.982812e-01;1.126378e+00
-2.983073e-01;1.126364e+00
-2.983333e-01;1.126350e+00
diff --git a/dG3D/benchmarks/LocalGurson_PathFollowing/Plane_notch.py b/dG3D/benchmarks/LocalGurson_PathFollowing/Plane_notch.py
index 86c8fb5ac38d597e10d8a3b0d8e65b97ea78cb47..bec947abda75e09e277a52ecb58cbb3f7db190b4 100644
--- a/dG3D/benchmarks/LocalGurson_PathFollowing/Plane_notch.py
+++ b/dG3D/benchmarks/LocalGurson_PathFollowing/Plane_notch.py
@@ -39,11 +39,10 @@ Gdn1 = LinearNucleationLaw(NucleationLawNum1, 0.2, 0.005, 1.)
# material law creation
-BulkLaw1 = LocalDamageGursonDG3DMaterialLaw(BulkLawNum1, rho, young, nu, q1,q2,q3,fVinitial, Harden1,1e-6,False,1e-8)
+BulkLaw1 = NonLocalDamageGursonDG3DMaterialLaw(BulkLawNum1, rho, young, nu, q1,q2,q3,fVinitial, Harden1,1e-6,False,1e-8)
+BulkLaw1.setNonLocalMethod(0)
BulkLaw1.setOrderForLogExp(3)
BulkLaw1.setNucleationLaw(Gdn1)
-BulkLaw1.setPredictorCorrectorParameters(2,100,1.,2)
-BulkLaw1.setFailureTolerance(0.9,1e-1,0.9)
BulkLaw1.setSubStepping(True,3)
# Solver parameters - shift
@@ -54,7 +53,7 @@ r_impl = 1.0
# Solver parameters - implicit
# ===================================================================================
soltype = 1 # StaticLinear=0 (default, StaticNonLinear=1, Explicit=2, # Multi=3, Implicit=4, Eigen=5)
-nstep = 300 # Number of step
+nstep = 70 # Number of step
tol=1.e-6 # Relative tolerance for NR scheme
tolAbs = 1.e-8 # Absolute tolerance for NR scheme
nstepArchIP=1 # Number of step between 2 archiving
@@ -150,14 +149,14 @@ mysolver.internalPointBuildView("local_fV",IPField.LOCAL_POROSITY,1,1)
mysolver.internalPointBuildView("local_fV_max",IPField.LOCAL_POROSITY,1,IPField.MAX_VALUE)
mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, 1)
mysolver.internalPointBuildView("triaxiality",IPField.STRESS_TRIAXIALITY, 1, 1)
-
+mysolver.internalPointBuildView("yield_fV",IPField.YIELD_POROSITY,1,1)
mysolver.archivingForceOnPhysicalGroup("Edge", 1400, 0, nstepArchForce)
mysolver.archivingNodeDisplacement(401,0, nstepArchForce)
mysolver.archivingNodeDisplacementOnPhysical(0,100,1,nstepArchForce)
-mysolver.archivingIPOnPhysicalGroup("Face", 700, IPField.DAMAGE,IPField.MAX_VALUE,nstepArchForce)
+mysolver.archivingIPOnPhysicalGroup("Face", 700, IPField.YIELD_POROSITY,IPField.MAX_VALUE,nstepArchForce)
mysolver.archivingIPOnPhysicalGroup("Face", 700, IPField.LOCAL_POROSITY,IPField.MAX_VALUE,nstepArchForce)
mysolver.solve()
check = TestCheck()
-check.equal(-8.171800e+02,mysolver.getArchivedForceOnPhysicalGroup("Edge", 1400, 0),1.e-6)
+check.equal(-4.319212e+02,mysolver.getArchivedForceOnPhysicalGroup("Edge", 1400, 0),1.e-6)
diff --git a/dG3D/benchmarks/LocalGurson_PathFollowing/force1400comp0.csv b/dG3D/benchmarks/LocalGurson_PathFollowing/force1400comp0.csv
deleted file mode 100644
index 277cb7a360cd8967b806a6889ac60aec6baedf6a..0000000000000000000000000000000000000000
--- a/dG3D/benchmarks/LocalGurson_PathFollowing/force1400comp0.csv
+++ /dev/null
@@ -1,213 +0,0 @@
-3.333333e-03;-9.990288e+02
-6.666667e-03;-2.253347e+03
-1.000000e-02;-2.406483e+03
-1.333333e-02;-2.510439e+03
-1.666667e-02;-2.589631e+03
-2.000000e-02;-2.653479e+03
-2.333333e-02;-2.706600e+03
-2.666667e-02;-2.766120e+03
-3.000000e-02;-2.816298e+03
-3.333333e-02;-2.872671e+03
-3.666667e-02;-2.934286e+03
-4.000000e-02;-2.999600e+03
-4.333333e-02;-3.066359e+03
-4.666667e-02;-3.131415e+03
-5.000000e-02;-3.190458e+03
-5.333333e-02;-3.233093e+03
-5.666667e-02;-3.256273e+03
-6.000000e-02;-3.265241e+03
-6.333333e-02;-3.263152e+03
-6.666667e-02;-3.254917e+03
-7.000000e-02;-3.241615e+03
-7.333333e-02;-3.223713e+03
-7.666667e-02;-3.201507e+03
-8.000000e-02;-3.167938e+03
-8.333333e-02;-3.127975e+03
-8.666667e-02;-3.081460e+03
-9.000000e-02;-3.039278e+03
-9.333333e-02;-2.992348e+03
-9.666667e-02;-2.940222e+03
-1.000000e-01;-2.882295e+03
-1.033333e-01;-2.817692e+03
-1.066667e-01;-2.760451e+03
-1.100000e-01;-2.710519e+03
-1.133333e-01;-2.664079e+03
-1.166667e-01;-2.619062e+03
-1.200000e-01;-2.570328e+03
-1.233333e-01;-2.517161e+03
-1.266667e-01;-2.455886e+03
-1.300000e-01;-2.369640e+03
-1.333333e-01;-2.258439e+03
-1.366667e-01;-2.212056e+03
-1.400000e-01;-2.172446e+03
-1.433333e-01;-2.115643e+03
-1.466667e-01;-2.033212e+03
-1.500000e-01;-1.999961e+03
-1.533333e-01;-1.965554e+03
-1.566667e-01;-1.917985e+03
-1.600000e-01;-1.879003e+03
-1.633333e-01;-1.825132e+03
-1.666667e-01;-1.788842e+03
-1.700000e-01;-1.751013e+03
-1.733333e-01;-1.700120e+03
-1.766667e-01;-1.649680e+03
-1.800000e-01;-1.579215e+03
-1.833333e-01;-1.558744e+03
-1.866667e-01;-1.542636e+03
-1.900000e-01;-1.521617e+03
-1.933333e-01;-1.494462e+03
-1.966667e-01;-1.460229e+03
-2.000000e-01;-1.417281e+03
-2.033333e-01;-1.363543e+03
-2.066667e-01;-1.312287e+03
-2.100000e-01;-1.247775e+03
-2.133333e-01;-1.186553e+03
-2.166667e-01;-1.130106e+03
-2.200000e-01;-1.058946e+03
-2.233333e-01;-9.924668e+02
-2.266667e-01;-9.311516e+02
-2.300000e-01;-8.744794e+02
-2.333333e-01;-8.220825e+02
-2.366667e-01;-7.741716e+02
-2.400000e-01;-7.296713e+02
-2.433333e-01;-6.873558e+02
-2.466667e-01;-6.473365e+02
-2.500000e-01;-6.103359e+02
-2.533333e-01;-5.650414e+02
-2.566667e-01;-5.238882e+02
-2.600000e-01;-4.745466e+02
-2.633333e-01;-4.307740e+02
-2.666667e-01;-3.917970e+02
-2.700000e-01;-3.569718e+02
-2.716667e-01;-3.477019e+02
-2.733333e-01;-3.431687e+02
-2.750000e-01;-3.372253e+02
-2.766667e-01;-3.294609e+02
-2.783333e-01;-3.192801e+02
-2.800000e-01;-3.088291e+02
-2.816667e-01;-2.959601e+02
-2.818750e-01;-2.949540e+02
-2.820833e-01;-2.941399e+02
-2.822917e-01;-2.930447e+02
-2.825000e-01;-2.915751e+02
-2.827083e-01;-2.896095e+02
-2.829167e-01;-2.869859e+02
-2.831250e-01;-2.834855e+02
-2.833333e-01;-2.787956e+02
-2.835417e-01;-2.723949e+02
-2.837500e-01;-2.632183e+02
-2.839583e-01;-2.578561e+02
-2.841667e-01;-2.536204e+02
-2.843750e-01;-2.480821e+02
-2.845833e-01;-2.409074e+02
-2.847917e-01;-2.317143e+02
-2.850000e-01;-2.200928e+02
-2.852083e-01;-2.056425e+02
-2.854167e-01;-1.880488e+02
-2.856250e-01;-1.671412e+02
-2.858333e-01;-1.546578e+02
-2.860417e-01;-1.395861e+02
-2.862500e-01;-1.218771e+02
-2.864583e-01;-1.114226e+02
-2.866667e-01;-9.892048e+01
-2.868750e-01;-8.996644e+01
-2.870833e-01;-7.934254e+01
-2.872917e-01;-7.180597e+01
-2.875000e-01;-6.502859e+01
-2.877083e-01;-5.707496e+01
-2.879167e-01;-5.024944e+01
-2.881250e-01;-4.265987e+01
-2.883333e-01;-3.942225e+01
-2.885417e-01;-3.560078e+01
-2.887500e-01;-3.125190e+01
-2.889583e-01;-2.828454e+01
-2.891667e-01;-2.490900e+01
-2.893750e-01;-2.206602e+01
-2.895833e-01;-1.892144e+01
-2.897917e-01;-1.772313e+01
-2.900000e-01;-1.662312e+01
-2.902083e-01;-1.529443e+01
-2.904167e-01;-1.410870e+01
-2.906250e-01;-1.271117e+01
-2.908333e-01;-1.148678e+01
-2.910417e-01;-1.006616e+01
-2.912500e-01;-8.844003e+00
-2.914583e-01;-7.993207e+00
-2.916667e-01;-7.027279e+00
-2.918750e-01;-6.784877e+00
-2.920833e-01;-6.479999e+00
-2.922917e-01;-6.103859e+00
-2.925000e-01;-5.759598e+00
-2.927083e-01;-5.345105e+00
-2.929167e-01;-4.861076e+00
-2.931250e-01;-4.669663e+00
-2.933333e-01;-4.490884e+00
-2.935417e-01;-4.323930e+00
-2.937500e-01;-4.167580e+00
-2.939583e-01;-3.978710e+00
-2.941667e-01;-3.768032e+00
-2.943750e-01;-3.543317e+00
-2.945833e-01;-3.310794e+00
-2.947917e-01;-3.233855e+00
-2.950000e-01;-3.142202e+00
-2.952083e-01;-3.036225e+00
-2.954167e-01;-2.944997e+00
-2.956250e-01;-2.841528e+00
-2.958333e-01;-2.754239e+00
-2.960417e-01;-2.693644e+00
-2.962500e-01;-2.657116e+00
-2.964583e-01;-2.623217e+00
-2.965625e-01;-2.612659e+00
-2.966146e-01;-2.610602e+00
-2.966667e-01;-2.608468e+00
-2.967187e-01;-2.606328e+00
-2.967708e-01;-2.604182e+00
-2.968229e-01;-2.602042e+00
-2.968750e-01;-2.599891e+00
-2.969271e-01;-2.597741e+00
-2.969792e-01;-2.595593e+00
-2.970312e-01;-2.593448e+00
-2.970833e-01;-2.591306e+00
-2.971354e-01;-2.589168e+00
-2.971875e-01;-2.587035e+00
-2.972396e-01;-2.584353e+00
-2.972917e-01;-2.580275e+00
-2.973437e-01;-2.578953e+00
-2.973698e-01;-2.579799e+00
-2.973958e-01;-2.580210e+00
-2.974219e-01;-2.580063e+00
-2.974479e-01;-2.579902e+00
-2.974740e-01;-2.579767e+00
-2.975000e-01;-2.579565e+00
-2.975260e-01;-2.579422e+00
-2.975521e-01;-2.579629e+00
-2.975781e-01;-2.579456e+00
-2.976042e-01;-2.578519e+00
-2.976302e-01;-2.579502e+00
-2.976562e-01;-2.578716e+00
-2.976823e-01;-2.580032e+00
-2.977083e-01;-2.577659e+00
-2.977344e-01;-2.578499e+00
-2.977604e-01;-2.578933e+00
-2.977865e-01;-2.578612e+00
-2.978125e-01;-2.578726e+00
-2.978385e-01;-2.578530e+00
-2.978646e-01;-2.578528e+00
-2.978906e-01;-2.578291e+00
-2.979167e-01;-2.578286e+00
-2.979427e-01;-2.578688e+00
-2.979687e-01;-2.576984e+00
-2.979948e-01;-2.576947e+00
-2.980208e-01;-2.578507e+00
-2.980469e-01;-2.575899e+00
-2.980729e-01;-2.576730e+00
-2.980990e-01;-2.576844e+00
-2.981250e-01;-2.576685e+00
-2.981510e-01;-2.576600e+00
-2.981771e-01;-2.576421e+00
-2.982031e-01;-2.576333e+00
-2.982292e-01;-2.576188e+00
-2.982552e-01;-2.576099e+00
-2.982812e-01;-2.575917e+00
-2.983073e-01;-2.575826e+00
-2.983333e-01;-2.575689e+00
diff --git a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
index e89c52643eba9737d5c73d6204ed47aed646e926..de13c7714c8a431be79952a41f60e2f448bc585a 100644
--- a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
+++ b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
@@ -924,10 +924,6 @@ void NonLocalPorousDuctileDamageDG3DMaterialLaw::setCorrectedRegularizedFunction
_nldPorous->setCorrectedRegularizedFunction(fct);
};
-void NonLocalPorousDuctileDamageDG3DMaterialLaw::setLocalRegularizedFunction(const scalarFunction& fct){
- _nldPorous->setLocalRegularizedFunction(fct);
-};
-
void NonLocalPorousDuctileDamageDG3DMaterialLaw::setNonLocalMethod(const int i) {
_nldPorous->setNonLocalMethod(i);
};
diff --git a/dG3D/src/nonLocalDamageDG3DMaterialLaw.h b/dG3D/src/nonLocalDamageDG3DMaterialLaw.h
index d0f74bdc351c9e890c344f686051b8d18d0df5b5..cd0ebf3745e03cd8697e326d8333a6d10371e920 100644
--- a/dG3D/src/nonLocalDamageDG3DMaterialLaw.h
+++ b/dG3D/src/nonLocalDamageDG3DMaterialLaw.h
@@ -345,7 +345,6 @@ class NonLocalPorousDuctileDamageDG3DMaterialLaw : public dG3DMaterialLaw{
void setPostBlockedDissipationBehavior(const int method);
void setFailureTolerance(const double NewFailureTol);
void setCorrectedRegularizedFunction(const scalarFunction& fct);
- void setLocalRegularizedFunction(const scalarFunction& fct);
void setNonLocalMethod(const int i);
void setShearPorosityGrowthFactor(const double k);
void setStressTriaxialityFunction_kw(const scalarFunction& fct);