diff --git a/NonLinearSolver/internalPoints/ipField.cpp b/NonLinearSolver/internalPoints/ipField.cpp
index 34b6c99cbf7a35bf9aea3dd550ce749fa0fea620..0bf49a6e404fe8291cc0e62051e2c145901b7e48 100644
--- a/NonLinearSolver/internalPoints/ipField.cpp
+++ b/NonLinearSolver/internalPoints/ipField.cpp
@@ -1061,6 +1061,7 @@ std::string IPField::ToString(const int i){
else if (i == PRESSURE) return "PRESSURE";
else if (i == viscousDissipatedEnergy) return "viscousDissipatedEnergy";
else if (i == mullinsDissipatedEnergy) return "mullinsDissipatedEnergy";
+ else if (i == yieldSurface) return "yieldSurface";
else if (i == USER1) return "USER1";
else if (i == USER2) return "USER2";
else if (i == USER3) return "USER3";
diff --git a/NonLinearSolver/internalPoints/ipField.h b/NonLinearSolver/internalPoints/ipField.h
index 9ef8b3b46042298a7c260003c73cd44a98419a4c..acc9a937a0119508fe5b1e8ef6e2da7240d1b27c 100644
--- a/NonLinearSolver/internalPoints/ipField.h
+++ b/NonLinearSolver/internalPoints/ipField.h
@@ -342,7 +342,7 @@ class IPField : public elementsField {
MAX_DEFO_ENERGY, MULLINS_DAMAGE,
PRESSURE,
// dissipated Energies
- viscousDissipatedEnergy,mullinsDissipatedEnergy,
+ viscousDissipatedEnergy,mullinsDissipatedEnergy, yieldSurface,
ONE, _EMEnergy,
UNDEFINED};
enum Operator { MEAN_VALUE=1, MIN_VALUE, MAX_VALUE, CRUDE_VALUE};
diff --git a/NonLinearSolver/internalPoints/ipHyperelastic.cpp b/NonLinearSolver/internalPoints/ipHyperelastic.cpp
index f654e0fa33f89cba40a1d032f72bab2f66160d59..e13a221894971b2e663ee921f9251eaa859f5784 100644
--- a/NonLinearSolver/internalPoints/ipHyperelastic.cpp
+++ b/NonLinearSolver/internalPoints/ipHyperelastic.cpp
@@ -15,7 +15,7 @@
IPHyperViscoElastic::IPHyperViscoElastic(const int N):IPVariableMechanics(),_N(N),_elasticEnergy(0.),_Ee(0.),_kirchhoff(0.),
_irreversibleEnergy(0.),_DirreversibleEnergyDF(0.), _viscousEnergyPart(0.), _dElasticEnergyPartdF(0.), _dViscousEnergyPartdF(0.), _elasticBulkPropertyScaleFactor(1.), _elasticShearPropertyScaleFactor(1.),
- _intA(0.), _intB(0.), _Jacobian(1.)
+ _intA(0.), _intB(0.), _Jacobian(1.), _yieldSurface(0.), _SVM(0.)
{
_A.clear();
_B.clear();
@@ -32,7 +32,7 @@ IPHyperViscoElastic::IPHyperViscoElastic(const IPHyperViscoElastic& src): IPVari
_irreversibleEnergy(src._irreversibleEnergy),_DirreversibleEnergyDF(src._DirreversibleEnergyDF),
_viscousEnergyPart(src._viscousEnergyPart), _dElasticEnergyPartdF(src._dElasticEnergyPartdF),
_dViscousEnergyPartdF(src._dViscousEnergyPartdF), _elasticBulkPropertyScaleFactor(src._elasticBulkPropertyScaleFactor), _elasticShearPropertyScaleFactor(src._elasticShearPropertyScaleFactor),
- _intA(src._intA),_intB(src._intB), _psi_branch(src._psi_branch), _Jacobian(src._Jacobian)
+ _intA(src._intA),_intB(src._intB), _psi_branch(src._psi_branch), _Jacobian(src._Jacobian), _yieldSurface(src._yieldSurface), _SVM(src._SVM)
{
};
@@ -58,6 +58,8 @@ IPHyperViscoElastic& IPHyperViscoElastic::operator =(const IPVariable& src){
_intB= psrc->_intB;
_psi_branch= psrc->_psi_branch;
_Jacobian=psrc->_Jacobian;
+ _yieldSurface=psrc->_yieldSurface;
+ _SVM=psrc->_SVM;
}
return *this;
};
@@ -81,6 +83,8 @@ void IPHyperViscoElastic::restart() {
restartManager::restart(_intB);
restartManager::restart(_psi_branch);
restartManager::restart(_Jacobian);
+ restartManager::restart(_yieldSurface);
+ restartManager::restart(_SVM);
};
void IPHyperViscoElastic::getViscoElasticStrain(int i, STensor3& Ev) const
diff --git a/NonLinearSolver/internalPoints/ipHyperelastic.h b/NonLinearSolver/internalPoints/ipHyperelastic.h
index d119c9071c88c419811346a64d4f8b597a8b825f..624467dcc0c88306e54f1348b8b74ceac8bc9f63 100644
--- a/NonLinearSolver/internalPoints/ipHyperelastic.h
+++ b/NonLinearSolver/internalPoints/ipHyperelastic.h
@@ -36,6 +36,8 @@ class IPHyperViscoElastic : public IPVariableMechanics{
std::vector<double> _psi_branch; // N elements
double _Jacobian; // NEW100
+ double _yieldSurface; // NEW100
+ double _SVM; // NEW100 von mises
protected:
//For energy sources
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
index a4f98f54443c93fbd61dc5356ae8e3baaf8d86ec..5996e752926c25808c5cd783a228f7b44635365c 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
@@ -762,6 +762,11 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
}
Ee *= 0.5;
+ // update yieldSurface in IP -> NEW100
+ q1->_yieldSurface = f;
+ q1->_SVM = PhiEq;
+ q1->_pressure = ptilde;
+
// Corrected corKir, mandel, backStress are already in the IP due to getIterated_DPhi
/*
// Correct A, B
@@ -2013,6 +2018,8 @@ void mlawNonLinearTVENonLinearTVP::getDho(const IPNonLinearTVP *q1, const double
for (int p=0; p<3; p++)
for (int q=0; q<3; q++)
dMeDphi(i,j,k,l) += dMedGN(i,j,p,q)*Gamma*(3.*_Idev(p,q,k,l) + 2.*_b/3 * 1./3 * _I(p,q)*_I(k,l));
+
+ dMeDphi *= 1/q1->_Jacobian; // NEW100
// Get dXdPhi
static STensor43 dXdPhi;
@@ -2027,7 +2034,6 @@ void mlawNonLinearTVENonLinearTVP::getDho(const IPNonLinearTVP *q1, const double
static STensor43 Dho2;
Dho2 = _I4;
Dho2 -= dMeDphi;
- Dho2 *= 1/q1->_Jacobian; // NEW100
Dho2 += dXdPhi;
STensorOperation::inverseSTensor43(Dho2,Dho2inv);
@@ -2051,9 +2057,10 @@ void mlawNonLinearTVENonLinearTVP::getDho(const IPNonLinearTVP *q1, const double
STensorOperation::zero(Ge_TrEeTensor_temp);
Ge_TrEeTensor_temp = (*Ge_TrEeTensor);
Ge_TrEeTensor_temp *= (1.5*Gamma);
+ Ge_TrEeTensor_temp *= 1/q1->_Jacobian; // NEW100
Dho2_u += Ge_TrEeTensor_temp;
}
- Dho2_u *= 1/q1->_Jacobian; // NEW100
+
STensorOperation::inverseSTensor43(Dho2_u,Dho2_u_inv);
// Get Dho2_v, Dho2_v_inv
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index c5c3875ee33b5b80b2d4da250943bc42333ba580..5f83b5992b5100ee99d0ba087735423d15949ec5 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -6240,6 +6240,7 @@ double NonLinearTVPDG3DIPVariable::get(const int i) const{
else if (i == IPField::Dev_corKir_Inf_YY){return _ipViscoElastoPlastic->_devCorKirinf_TVE(1,1);}
else if (i == IPField::Dev_corKir_Inf_YZ){return _ipViscoElastoPlastic->_devCorKirinf_TVE(1,2);}
else if (i == IPField::Dev_corKir_Inf_ZZ){return _ipViscoElastoPlastic->_devCorKirinf_TVE(2,2);}
+ else if (i == IPField::PRESSURE){return _ipViscoElastoPlastic->_pressure;}
else if (i == IPField::viscousDissipatedEnergy){return _ipViscoElastoPlastic->_viscousDissipatedEnergy;}
else if (i == IPField::mullinsDissipatedEnergy){return _ipViscoElastoPlastic->_mullinsDissipatedEnergy;}
else if (i == IPField::Eve1_XX){
@@ -6379,6 +6380,14 @@ double NonLinearTVPDG3DIPVariable::get(const int i) const{
{
return _ipViscoElastoPlastic->_mullinsDamage;
}
+ else if (i == IPField::yieldSurface)
+ {
+ return _ipViscoElastoPlastic->_yieldSurface;
+ }
+ else if (i == IPField::SVM)
+ {
+ return _ipViscoElastoPlastic->_SVM;
+ }
else
{
return dG3DIPVariable::get(i);
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 6ed9c64645d0c5bec195a60db135b56d76f3c012..0483ec109bc3066d9711853e0cf6ac739952f007 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -247,7 +247,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
_stressLaw->stress(ipvcur,ipvprev,false,checkfrac,dTangent);
-#if 0
+#if 1
_stressLaw->stress(ipvcur,ipvprev,true,checkfrac,dTangent);
ipvcur->getConstRefToDeformationGradient().print("F Analytique"); // FLE
ipvcur->getConstRefToFirstPiolaKirchhoffStress().print("P Analytique"); // FLE
@@ -622,7 +622,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
}
}
-#if 0
+#if 1
ipvcur->getConstRefToDeformationGradient().print("F Numerique"); // FLE
ipvcur->getConstRefToFirstPiolaKirchhoffStress().print("P Numerique"); // FLE
ipvcur->getConstRefToTangentModuli().print("dPdE Numerique");