diff --git a/NonLinearSolver/internalPoints/ipField.cpp b/NonLinearSolver/internalPoints/ipField.cpp
index 1e232cc6b7caae0655245948a59f4e9a174b5e2e..f653c6f9a7aa1a520305f97e89dca3a0ed71f24e 100644
--- a/NonLinearSolver/internalPoints/ipField.cpp
+++ b/NonLinearSolver/internalPoints/ipField.cpp
@@ -112,6 +112,16 @@ std::string IPField::ToString(const int i){
else if (i == P_ZX) return "P_ZX";
else if (i == P_ZY) return "P_ZY";
else if (i == P_ZZ) return "P_ZZ";
+
+ else if (i == P_XX_EFF) return "P_XX_EFF";
+ else if (i == P_XY_EFF) return "P_XY_EFF";
+ else if (i == P_XZ_EFF) return "P_XZ_EFF";
+ else if (i == P_YX_EFF) return "P_YX_EFF";
+ else if (i == P_YY_EFF) return "P_YY_EFF";
+ else if (i == P_YZ_EFF) return "P_YZ_EFF";
+ else if (i == P_ZX_EFF) return "P_ZX_EFF";
+ else if (i == P_ZY_EFF) return "P_ZY_EFF";
+ else if (i == P_ZZ_EFF) return "P_ZZ_EFF";
// for VT case (VM or VLM)
diff --git a/NonLinearSolver/internalPoints/ipField.h b/NonLinearSolver/internalPoints/ipField.h
index a9e276e0391293854ae7ff4d7ce0fd7d056a993e..7fa320ea097e5cb372309c9542c2d041ff250087 100644
--- a/NonLinearSolver/internalPoints/ipField.h
+++ b/NonLinearSolver/internalPoints/ipField.h
@@ -122,7 +122,8 @@ class IPField : public elementsField {
U_XX, U_XY, U_XZ, U_YY, U_YZ, U_ZZ,
GL_XX, GL_XY, GL_XZ, GL_YY, GL_YZ, GL_ZZ,
F_XX, F_XY, F_XZ, F_YX, F_YY,F_YZ, F_ZX, F_ZY, F_ZZ,
- P_XX, P_XY, P_XZ,P_YX, P_YY, P_YZ, P_ZX, P_ZY, P_ZZ,
+ P_XX, P_XY, P_XZ,P_YX, P_YY, P_YZ, P_ZX, P_ZY, P_ZZ,
+ P_XX_EFF, P_XY_EFF, P_XZ_EFF,P_YX_EFF, P_YY_EFF, P_YZ_EFF, P_ZX_EFF, P_ZY_EFF, P_ZZ_EFF,
// for VT case (VM or VLM)
// for each VT phase
diff --git a/NonLinearSolver/materialLaw/STensorOperations.h b/NonLinearSolver/materialLaw/STensorOperations.h
index 2caa31c3011f2eb54069648a5b967c09c9746dd1..fc7251f13a2fc4162903fa9ff1f6c9cd9795777a 100644
--- a/NonLinearSolver/materialLaw/STensorOperations.h
+++ b/NonLinearSolver/materialLaw/STensorOperations.h
@@ -342,7 +342,7 @@ namespace STensorOperation{
//here some with 1/2
};
- inline void inverseSTensor3(const STensor3 &a, STensor3 &ainv)
+ inline void inverseSTensor3(const STensor3 &a, STensor3 &ainv, STensor43* DainvDa=NULL, bool sym=true)
{
double udet = 1./determinantSTensor3(a);
ainv(0,0) = (a(1,1) * a(2,2) - a(1,2) * a(2,1))* udet;
@@ -354,6 +354,31 @@ namespace STensorOperation{
ainv(0,2) = (a(0,1) * a(1,2) - a(0,2) * a(1,1))* udet;
ainv(1,2) = -(a(0,0) * a(1,2) - a(0,2) * a(1,0))* udet;
ainv(2,2) = (a(0,0) * a(1,1) - a(0,1) * a(1,0))* udet;
+
+ if (DainvDa != NULL)
+ {
+ for (int i=0; i< 3; i++)
+ {
+ for (int j=0; j< 3; j++)
+ {
+ for (int k=0; k< 3; k++)
+ {
+ for (int l=0; l< 3; l++)
+ {
+ if (sym)
+ {
+ (*DainvDa)(i,j,k,l) = -0.5*(ainv(i,k)*ainv(j,l)+ ainv(i,l)*ainv(j,k));
+ }
+ else
+ {
+ (*DainvDa)(i,j,k,l) = -ainv(i,k)*ainv(l,j);
+ }
+ }
+ }
+ }
+ }
+ }
+
};
inline void inverseAndTransposeSTensor3(const STensor3 &a, STensor3 &ainvT)
diff --git a/NonLinearSolver/materialLaw/elasticPotential.cpp b/NonLinearSolver/materialLaw/elasticPotential.cpp
index a82d698775ffd506940c002d599048fc7fdd4676..b34a1f3f78df8996e06df64b68f566244ccd9f7a 100644
--- a/NonLinearSolver/materialLaw/elasticPotential.cpp
+++ b/NonLinearSolver/materialLaw/elasticPotential.cpp
@@ -1617,6 +1617,135 @@ void NeoHookeanElasticPotential::getLocalPositiveValForNonLocalEquation(const ST
}
};
+
+CompressiveNeoHookeanPotential::CompressiveNeoHookeanPotential(const int num, const double E, const double nu) :
+ elasticPotential(num), _E(E), _nu(nu){
+ double K = _E/(3.*(1.-2.*_nu));
+ _mu = _E/(2.*(1.+_nu));
+ _lambda = K- 2.*_mu/3.;
+ static const STensor3 I(1.);
+ STensorOperation::zero(_Cel);
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ for (int k=0; k<3; k++){
+ for (int l=0; l<3; l++){
+ _Cel(i,j,k,l) = _lambda*I(i,j)*I(k,l)+_mu*(I(i,k)*I(j,l)+I(i,l)*I(j,k));
+ }
+ }
+ }
+ }
+};
+
+CompressiveNeoHookeanPotential::CompressiveNeoHookeanPotential(const CompressiveNeoHookeanPotential& src) : elasticPotential(src),
+ _E(src._E),_nu(src._nu),_lambda(src._lambda),_mu(src._mu), _Cel(src._Cel){};
+
+double CompressiveNeoHookeanPotential::get(const STensor3& F) const{
+ double J = F.determinant();
+ static STensor3 C, invC;
+ STensorOperation::multSTensor3FirstTranspose(F,F,C);
+ STensorOperation::inverseSTensor3(C,invC);
+ double traceC = C.trace();
+ return 0.5*_lambda*log(J)*log(J) - _mu*log(J) + 0.5*_mu*(traceC-3);
+};
+void CompressiveNeoHookeanPotential::constitutive(const STensor3& F, STensor3& P, const bool stiff, STensor43* DPDF, const bool dTangent, STensor63* dCalgdeps) const{
+ double J = F.determinant();
+ double logJ = log(J);
+ static STensor3 C, invC, S;
+ STensorOperation::multSTensor3FirstTranspose(F,F,C);
+ static STensor43 DinvCDC;
+ if (stiff)
+ {
+ STensorOperation::inverseSTensor3(C,invC,&DinvCDC,true);
+ }
+ else
+ {
+ STensorOperation::inverseSTensor3(C,invC);
+ }
+
+ //S=mu*(I-invC)+lambda*(logJ)*invC
+ STensorOperation::diag(S,_mu);
+ STensorOperation::axpy(S,invC,-_mu+_lambda*logJ);
+ STensorOperation::multSTensor3(F, S, P);
+ if (stiff)
+ {
+ static STensor43 DSDC;
+ DSDC = DinvCDC;
+ DSDC *= (_lambda*logJ-_mu);
+ STensorOperation::prodAdd(invC,invC,0.5*_lambda,DSDC);
+ STensorOperation::zero(*DPDF);
+ for (int i = 0; i < 3; i++) {
+ for (int J = 0; J < 3; J++) {
+ for (int k = 0; k < 3; k++) {
+ for (int L = 0; L < 3; L++) {
+ for (int K = 0; K < 3; K++) {
+ for (int I = 0; I < 3; I++) {
+ (*DPDF)(i, J, k, L) += 2.*DSDC(I, J, K, L) * F(i,I) * F(k,K);
+ }
+ }
+ if (i == k) {
+ (*DPDF)(i, J, k, L) += S(J,L);
+ }
+ }
+ }
+ }
+ }
+ }
+ if (dTangent)
+ {
+ Msg::Error("compute dCalgdeps is not implemented in CompressiveNeoHookeanPotential");
+ }
+};
+void CompressiveNeoHookeanPotential::getLocalValForNonLocalEquation(const STensor3& F, double& val, const bool stiff, STensor3* DvalDF) const{
+ double Y = get(F);
+ double E = getYoungModulus();
+ val = sqrt(2*Y/E);
+ if (stiff){
+ STensorOperation::zero(*DvalDF);
+ this->constitutive(F,*DvalDF,false,NULL);
+ if (val > 1e-10){
+ (*DvalDF) *= (1./(E*val));
+ }
+ else{
+ (*DvalDF) *= (1./(E*1e-10));
+ }
+ }
+};
+
+
+// for anisotropic
+double CompressiveNeoHookeanPotential::getPositiveValue(const STensor3& F) const {
+ return get(F);
+};
+double CompressiveNeoHookeanPotential::getNegativeValue(const STensor3& F) const {
+ return 0;
+};
+void CompressiveNeoHookeanPotential::constitutivePositive(const STensor3& F, STensor3& P, const bool stiff, STensor43* L) const {
+ constitutive(F,P,stiff,L);
+};
+void CompressiveNeoHookeanPotential::constitutiveNegative(const STensor3& F, STensor3& P, const bool stiff, STensor43* L) const {
+ STensorOperation::zero(P);
+ if (stiff){
+ STensorOperation::zero(*L);
+ }
+};
+void CompressiveNeoHookeanPotential::getLocalPositiveValForNonLocalEquation(const STensor3& F, double& val, const bool stiff, STensor3* DvalDF) const {
+ double Y = getPositiveValue(F);
+ double E = getYoungModulus();
+ val = sqrt(2*Y/E);
+ if (stiff){
+ STensorOperation::zero(*DvalDF);
+ this->constitutivePositive(F,*DvalDF,false,NULL);
+ if (val > 1e-10){
+ (*DvalDF) *= (1./(E*val));
+ }
+ else{
+ (*DvalDF) *= (1./(E*1e-10));
+ }
+
+ }
+};
+
+
HyperElasticPolynomialPotential::HyperElasticPolynomialPotential(const int num, const double E, const double nu, const bool pert, const double tol):
elasticPotential(num), _E(E), _nu(nu), _ByPerturbation(pert), _tol(tol)
{
diff --git a/NonLinearSolver/materialLaw/elasticPotential.h b/NonLinearSolver/materialLaw/elasticPotential.h
index 8ea2239aa8e4eaa613a825689d914c83cd561bc5..98b3202f5f8173ee8f6bea0fcd568cb6e853e5d4 100644
--- a/NonLinearSolver/materialLaw/elasticPotential.h
+++ b/NonLinearSolver/materialLaw/elasticPotential.h
@@ -13,6 +13,7 @@
#include "STensor3.h"
#include "STensor43.h"
#include "STensorOperations.h"
+#include <limits>
class elasticPotential{
#ifndef SWIG
@@ -236,6 +237,35 @@ class NeoHookeanElasticPotential : public elasticPotential{
#endif //SWIG
};
+class CompressiveNeoHookeanPotential : public elasticPotential{
+ protected:
+ double _E, _nu, _lambda, _mu;
+ STensor43 _Cel;
+
+ public:
+ CompressiveNeoHookeanPotential(const int num, const double E, const double nu);
+ #ifndef SWIG
+ CompressiveNeoHookeanPotential(const CompressiveNeoHookeanPotential& src);
+ virtual ~CompressiveNeoHookeanPotential(){}
+ virtual const STensor43& getElasticTensor() const {return _Cel;};
+ virtual double getYoungModulus() const {return _E;};
+ virtual double getPoissonRatio() const {return _nu;};
+ virtual double get(const STensor3& F) const; // elastic potential
+ virtual void constitutive(const STensor3& F, STensor3& P, const bool sitff, STensor43* L,
+ const bool dTangent=false, STensor63* dCalgdeps = NULL) const; // constitutive law
+ virtual void getLocalValForNonLocalEquation(const STensor3& F, double& val, const bool stiff, STensor3* DvalDF) const;
+ virtual elasticPotential* clone() const {
+ return new CompressiveNeoHookeanPotential(*this);
+ }
+ // for anisotropic
+ virtual double getPositiveValue(const STensor3& F) const;
+ virtual double getNegativeValue(const STensor3& F) const;
+ virtual void constitutivePositive(const STensor3& F, STensor3& P, const bool sitff, STensor43* L) const;
+ virtual void constitutiveNegative(const STensor3& F, STensor3& P, const bool sitff, STensor43* L) const;
+ virtual void getLocalPositiveValForNonLocalEquation(const STensor3& F, double& val, const bool stiff, STensor3* DvalDF) const;
+ #endif //SWIG
+};
+
// See https://en.wikipedia.org/wiki/Polynomial_(hyperelastic_model)
class HyperElasticPolynomialPotential: public elasticPotential
{
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
index a4b3b492597d67431faac86e9960415a3c25a9dd..7a5c788b19900f5c6ae52013e86a2c2057ca23e9 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
@@ -1207,8 +1207,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
return;
}
};
-
- //std::cout << "############################################### "<<std::endl;
+
q1->_DgammaDt = Dgamma/this->getTimeStep();
// update effective stress tensor
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index a4cb803e49fcbb1acc43db32ca8a18aa98591ad7..9bca682e9997a4771ab4037db084140edf5626f2 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -8332,15 +8332,41 @@ void nonLinearMechSolver::archivingNodeIPOnPhysicalGroup(const int dim, const in
}
};
+/**
+ * @brief Archiving a value IP on an interface element. The result is stored under a csv file.
+ * @param elenumMinus, elenumPlus: interface between two bulk elements
+ * @param ipval: a quantity specified in IPField::Ouput
+ * @param elemval: an operation to perform: 1 - mean value, 2 - min value, 3 - max value, and aribtrary another value - value at a specific Gausspoint specifyied by numip
+ * @param elemval: an operation to perform: 1 - mean value, 2 - min value, 3 - max value, and aribtrary another value - value at a specific Gausspoint specifyied by numip
+ * @param numip (0 by default): specifying a IP to archive, ranging from 0 to npts-1. This optiton is considered if elemval is not 1, 2, or 3
+ * @param nstep (1 by default): archiving frequence after each nstep time steps
+ */
void nonLinearMechSolver::archivingInterfaceElementIP(const int elenumMinus,const int elenumPlus, const int ipval,const int elemval, const int numip, const int nstep){
vaip.push_back(IPField::ip2archive(-1, elenumMinus, elenumPlus,ipval,elemval,nstep,numip));
};
+
+/**
+ * @brief Archiving a value IP on a bulk element. . The result is stored under a csv file.
+ * @param elenum: element number
+ * @param ipval: a quantity specified in IPField::Ouput
+ * @param elemval: an operation to perform between 1 - mean value, 2 - min value, 3 - max value, and aribtrary another value - the value at a specific Gausspoint specifyied by numip
+ * @param numip (0 by default): specifying a IP to archive, ranging from 0 to npts-1. This optiton is considered if elemval is not 1, 2, or 3
+ * @param nstep (1 by default): archiving frequence after each nstep time steps
+ */
void nonLinearMechSolver::archivingElementIP(const int elenum,const int ipval,const int elemval, const int numip, const int nstep){
// save data on numip on elemement elenum
vaip.push_back(IPField::ip2archive(-1, elenum, elenum,ipval,elemval,nstep,numip));
};
+/**
+ * @brief Archiving a value on a physical (only mean, min, max over this physical are available).
+ * @param onwhat: entity identification between Node, Edge, Face, or Volume
+ * @param numphys: physical number of the entity
+ * @param ipval: a quantity specified in IPField::Ouput
+ * @param elemval: an operation to perform between 1 - mean value, 2 - min value, 3 - max value
+ * @param nstep (1 by default): archiving frequence after each nstep time steps
+ */
void nonLinearMechSolver::archivingIPOnPhysicalGroup(std::string onwhat, const int numphys, const int ipval,const int elemval, const int nstep){
const std::string node("Node");
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 2cad40102a6a1c0efddf19adc50e88132a963684..8ce50726a688ff37116e673383591de6f96ec008 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -849,6 +849,44 @@ double dG3DIPVariable::get(const int comp) const
else if (comp == IPField::P_ZZ){
return P(2,2);
}
+ // effective P
+ else if (comp == IPField::P_XX_EFF){
+ double Dam = get(IPField::DAMAGE);
+ return P(0,0)/(1.0-Dam);
+ }
+ else if (comp == IPField::P_XY_EFF){
+ double Dam = get(IPField::DAMAGE);
+ return P(0,1)/(1.0-Dam);
+ }
+ else if (comp == IPField::P_XZ_EFF){
+ double Dam = get(IPField::DAMAGE);
+ return P(0,2)/(1.0-Dam);
+ }
+ else if (comp == IPField::P_YX_EFF){
+ double Dam = get(IPField::DAMAGE);
+ return P(1,0)/(1.0-Dam);
+ }
+ else if (comp == IPField::P_YY_EFF){
+ double Dam = get(IPField::DAMAGE);
+ return P(1,1)/(1.0-Dam);
+ }
+ else if (comp == IPField::P_YZ_EFF){
+ double Dam = get(IPField::DAMAGE);
+ return P(1,2)/(1.0-Dam);
+ }
+ else if (comp == IPField::P_ZX_EFF){
+ double Dam = get(IPField::DAMAGE);
+ return P(2,0)/(1.0-Dam);
+ }
+ else if (comp == IPField::P_ZY_EFF){
+ double Dam = get(IPField::DAMAGE);
+ return P(2,1)/(1.0-Dam);
+ }
+ else if (comp == IPField::P_ZZ_EFF){
+ double Dam = get(IPField::DAMAGE);
+ return P(2,2)/(1.0-Dam);
+ }
+
else if (comp == IPField::S_XX)
{
return secondPK(0,0);
@@ -873,30 +911,30 @@ double dG3DIPVariable::get(const int comp) const
{
return secondPK(2,2);
}
- else if (comp == IPField::JACOBIAN){
+ else if (comp == IPField::JACOBIAN){
return getJ();
- }
- else if (comp == IPField::PRESSION){
+ }
+ else if (comp == IPField::PRESSION){
return (cauchy(0,0)+cauchy(1,1)+cauchy(2,2))/3.;
- }
- else if (comp == IPField::FIRST_PRINCIPAL_STRESS){
+ }
+ else if (comp == IPField::FIRST_PRINCIPAL_STRESS){
static fullMatrix<double> eigVec(3,3);
static fullVector<double> eigVal(3);
cauchy.eig(eigVec,eigVal,true);
return eigVal(0);
- }
- else if (comp == IPField::SECOND_PRINCIPAL_STRESS){
+ }
+ else if (comp == IPField::SECOND_PRINCIPAL_STRESS){
static fullMatrix<double> eigVec(3,3);
static fullVector<double> eigVal(3);
cauchy.eig(eigVec,eigVal,true);
return eigVal(1);
- }
- else if (comp == IPField::THIRD_PRINCIPAL_STRESS){
+ }
+ else if (comp == IPField::THIRD_PRINCIPAL_STRESS){
static fullMatrix<double> eigVec(3,3);
static fullVector<double> eigVal(3);
cauchy.eig(eigVec,eigVal,true);
return eigVal(2);
- }
+ }
else if (comp == IPField::LODE_PARAMETER){
static fullMatrix<double> eigVec(3,3);
static fullVector<double> eigVal(3);
@@ -914,21 +952,21 @@ double dG3DIPVariable::get(const int comp) const
const double pi = 3.14159265358979323846;
return (acos(xi))/3.*(180./pi);
}
- else if (comp == IPField::DISP_JUMP_X){
+ else if (comp == IPField::DISP_JUMP_X){
return getConstRefToJump()(0);
- }
- else if (comp == IPField::DISP_JUMP_Y){
+ }
+ else if (comp == IPField::DISP_JUMP_Y){
return getConstRefToJump()(1);
- }
- else if (comp == IPField::DISP_JUMP_Z){
+ }
+ else if (comp == IPField::DISP_JUMP_Z){
return getConstRefToJump()(2);
- }
- else if (comp == IPField::STRESS_TRIAXIALITY){
+ }
+ else if (comp == IPField::STRESS_TRIAXIALITY){
double p = (cauchy(0,0)+cauchy(1,1)+cauchy(2,2))/3.;
double svm = vonMises();
if (p ==0.) return 0.;
else return p/svm;
- }
+ }
else if (comp == IPField::INCOMPATIBLE_STRAIN_X){
const SVector3& istrain = this->getConstRefToIncompatibleJump();
return istrain(0);
@@ -941,24 +979,24 @@ double dG3DIPVariable::get(const int comp) const
const SVector3& istrain = this->getConstRefToIncompatibleJump();
return istrain(2);
}
- else if (comp == IPField::DAMAGE_IS_BLOCKED){
- return dissipationIsBlocked();
- }
+ else if (comp == IPField::DAMAGE_IS_BLOCKED){
+ return dissipationIsBlocked();
+ }
else if (comp == IPField::ACTIVE_DISSIPATION){
return this->dissipationIsActive();
}
- else if (comp == IPField::DEFO_ENERGY){
- return this->defoEnergy();
- }
- else if (comp == IPField::PLASTIC_ENERGY){
- return this->plasticEnergy();
- }
+ else if (comp == IPField::DEFO_ENERGY){
+ return this->defoEnergy();
+ }
+ else if (comp == IPField::PLASTIC_ENERGY){
+ return this->plasticEnergy();
+ }
else if (comp == IPField::DAMAGE_ENERGY){
return this->damageEnergy();
}
- else if (comp == IPField::IRREVERSIBLE_ENERGY){
- return this->irreversibleEnergy();
- }
+ else if (comp == IPField::IRREVERSIBLE_ENERGY){
+ return this->irreversibleEnergy();
+ }
else if (comp == IPField::DELETED){
if (this->isDeleted()) return 1.;
else return 0.;
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index 016efa98814b6fc3917c75e598fb7a560cd2ddc0..23ddec85816a0bfc478f85b9c0d1d8a371a6c970 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -2289,6 +2289,33 @@ class LinearElasticDG3DIPVariable :public dG3DIPVariable{
virtual IPVariable* clone() const {return new LinearElasticDG3DIPVariable(*this);};
};
+class HyperElasticDG3DIPVariable :public dG3DIPVariable{
+ protected:
+ double _defoEnergy;
+ public:
+ HyperElasticDG3DIPVariable(const bool createBodyForceHO, const bool oninter):dG3DIPVariable(createBodyForceHO,oninter), _defoEnergy(0){}; // A empty constructor is mandatory to build FractureCohesive3DIPVariable OK like this as the fracture is on interface??
+ HyperElasticDG3DIPVariable(const HyperElasticDG3DIPVariable &source):dG3DIPVariable(source), _defoEnergy(source._defoEnergy){};
+ virtual HyperElasticDG3DIPVariable &operator = (const IPVariable &_source){
+ dG3DIPVariable::operator =(_source);
+ const HyperElasticDG3DIPVariable* psrc = dynamic_cast<const HyperElasticDG3DIPVariable*>(&_source);
+ if (psrc != NULL)
+ {
+ _defoEnergy = psrc->_defoEnergy;
+ }
+ return *this;
+ };
+ virtual ~HyperElasticDG3DIPVariable()
+ {
+
+ }
+ virtual IPVariable* getInternalData() {return NULL;};
+ virtual const IPVariable* getInternalData() const {return NULL;};
+ virtual double defoEnergy() const {return _defoEnergy;};
+ virtual double& getRefToDefoEnergy() {return _defoEnergy;};
+ virtual IPVariable* clone() const {return new HyperElasticDG3DIPVariable(*this);};
+};
+
+
class dG3DIPVariableVoid :public dG3DIPVariable{
protected:
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 455a815151cb3087d1c95f274e5022c52d777f7e..3fa7e33a237cbaaacca069955cbbb9b75440661b 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -1046,9 +1046,9 @@ void dG3DHyperelasticMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyF
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
- IPVariable* ipvi = new LinearElasticDG3DIPVariable(hasBodyForce,inter);
- IPVariable* ipv1 = new LinearElasticDG3DIPVariable(hasBodyForce,inter);
- IPVariable* ipv2 = new LinearElasticDG3DIPVariable(hasBodyForce,inter);
+ IPVariable* ipvi = new HyperElasticDG3DIPVariable(hasBodyForce,inter);
+ IPVariable* ipv1 = new HyperElasticDG3DIPVariable(hasBodyForce,inter);
+ IPVariable* ipv2 = new HyperElasticDG3DIPVariable(hasBodyForce,inter);
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
}
@@ -1061,20 +1061,21 @@ void dG3DHyperelasticMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBod
if(iele == NULL){
inter=false;
}
- ipv = new LinearElasticDG3DIPVariable(hasBodyForce,inter);
+ ipv = new HyperElasticDG3DIPVariable(hasBodyForce,inter);
}
void dG3DHyperelasticMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac, const bool dTangent)
{
/* get ipvariable */
- dG3DIPVariableBase* ipvcur = static_cast<dG3DIPVariableBase*>(ipv);;
- const dG3DIPVariableBase* ipvprev = static_cast<const dG3DIPVariableBase*>(ipvp);;
+ HyperElasticDG3DIPVariable* ipvcur = static_cast<HyperElasticDG3DIPVariable*>(ipv);;
+ const HyperElasticDG3DIPVariable* ipvprev = static_cast<const HyperElasticDG3DIPVariable*>(ipvp);;
/* strain xyz */
const STensor3& F = ipvcur->getConstRefToDeformationGradient();
STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress();
STensor43& H = ipvcur->getRefToTangentModuli();
_elasticLaw->constitutive(F,P,stiff,&H);
+ ipvcur->getRefToDefoEnergy() = _elasticLaw->get(F);
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}