diff --git a/NonLinearSolver/internalPoints/CMakeLists.txt b/NonLinearSolver/internalPoints/CMakeLists.txt
index 16a62825ae717e028d14719bc1f9f9a9f3940ee5..357e3d0ad1991d84ef8c1b39ba5b48172be2c185 100644
--- a/NonLinearSolver/internalPoints/CMakeLists.txt
+++ b/NonLinearSolver/internalPoints/CMakeLists.txt
@@ -32,8 +32,9 @@ set(SRC
ipJ2SmallStrains.cpp
ipFiniteStrain.cpp
ipSMP.cpp
+ ipPhenomenologicalSMP.cpp
ipViscoelastic.cpp
- ipvariable.cpp
+ ipvariable.cpp
# inline
ipvariable.h
ipEOS.cpp
@@ -41,7 +42,7 @@ set(SRC
ipHyperelastic.cpp
ipNonLocalDamageHyperelastic.cpp
ipElecSMP.cpp
- ipKinematicHardening.cpp
+ ipKinematicHardening.cpp
ipAnIsotropicTherMech.cpp
)
diff --git a/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.cpp b/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..117d87bebaedca8f00a580cb31f74fec3e34d675
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.cpp
@@ -0,0 +1,67 @@
+//
+// Description: storing class for Phenomenological SMP
+//
+//
+// Author: <M. Pareja>, (C) 2017
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#include "ipPhenomenologicalSMP.h"
+#include "mlawSMP.h"
+#include "restartManager.h"
+#include "STensorOperations.h"
+
+IPPhenomenologicalSMP::IPPhenomenologicalSMP(double _zg):IPVariableMechanics(),_thermalEnergy(0.) , _eps(0.),_SMPEnergy(0.)
+{
+ STensorOperation::unity(Fpg);
+ STensorOperation::unity(Fp);
+ STensorOperation::unity(Ff);
+ zg=_zg;
+};
+
+
+IPPhenomenologicalSMP::IPPhenomenologicalSMP(const IPPhenomenologicalSMP &source) :IPVariableMechanics()//: IPLinearThermoMechanics(source) //n
+{
+ _SMPEnergy=source._SMPEnergy;//n
+ Fpg = source.Fpg;
+ Fp = source.Fp;
+ Ff = source.Ff;
+ zg = source.zg;
+ _thermalEnergy= source._thermalEnergy;
+ _eps= source._eps;
+ _SMPEnergy=source._SMPEnergy;
+}
+
+
+IPPhenomenologicalSMP& IPPhenomenologicalSMP::operator=(const IPVariable &source)
+{
+ IPVariableMechanics::operator=(source);
+
+ const IPPhenomenologicalSMP* src = dynamic_cast<const IPPhenomenologicalSMP*>(&source);
+ if(src)
+ {
+ Fpg = src->Fpg;
+ Fp = src->Fp;
+ Ff = src->Ff;
+ zg = src->zg;
+ _thermalEnergy = src->_thermalEnergy;
+ _eps=src->_eps;
+ _SMPEnergy=src->_SMPEnergy;
+ }
+ return *this;
+}
+
+
+void IPPhenomenologicalSMP::restart()
+{
+ IPVariableMechanics::restart();
+ restartManager::restart(Fpg);
+ restartManager::restart(Fp);
+ restartManager::restart(Ff);
+ restartManager::restart(zg);
+ restartManager::restart(_SMPEnergy);
+ restartManager::restart(_thermalEnergy);
+ restartManager::restart(_eps);
+ return;
+}
diff --git a/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.h b/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.h
new file mode 100644
index 0000000000000000000000000000000000000000..e096c069e505c25066f79c9c5f9173d49432db79
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.h
@@ -0,0 +1,56 @@
+//
+// Description: storing class for linear thermomechanical law it is a "virtual implementation" you have to define an ipvariable
+// in your project which derive from this one.
+// Author: <M. Pareja>, (C) 2017
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef IPPHENOMENOLOGICALSMP_H_
+#define IPPHENOMENOLOGICALSMP_H_
+#include "ipLinearThermoMechanics.h"
+#include "STensor3.h"
+#include "ipvariable.h"
+//#include "mlawSMP.h"
+
+class IPPhenomenologicalSMP : public IPVariableMechanics//: public IPLinearThermoMechanics //, public mlawSMP
+{
+public:
+ double _SMPEnergy;
+protected:
+ STensor3 Fpg,Fp,Ff;
+ double zg;
+ //double _SMPEnergy;
+ double _thermalEnergy;
+ double _eps; // equivalent plastic strain
+ public:
+ IPPhenomenologicalSMP(double _zg);
+
+ //IPSMP(const materialLaw* l = NULL);
+ IPPhenomenologicalSMP(const IPPhenomenologicalSMP &source);
+ IPPhenomenologicalSMP& operator=(const IPVariable &source);
+ //virtual double defoEnergy() const ;
+ virtual double defoEnergy() const{return _SMPEnergy;}
+ virtual double getThermalEnergy() const {return _thermalEnergy;};
+
+ virtual STensor3 &getRefToFpg() {return Fpg;}
+ virtual const STensor3 &getConstRefToFpg() const {return Fpg;}
+ virtual STensor3 &getRefToFp() {return Fp;}
+ virtual const STensor3 &getConstRefToFp() const {return Fp;}
+ virtual STensor3 &getRefToFf() {return Ff;}
+ virtual const STensor3 &getConstRefToFf() const {return Ff;}
+ virtual double &getRefTozg() {return zg;}
+ virtual double getzg() const {return zg;}
+ virtual void restart();
+ virtual IPVariable* clone() const {return new IPPhenomenologicalSMP(*this);};
+
+ virtual double& getRefToEquivalentPlasticDefo() {return _eps;};
+ virtual const double& getConstRefToEquivalentPlasticDefo() const {return _eps;};
+
+
+protected:
+
+};
+
+#endif // IPPHENOMENOLOGICALSMP_H_
diff --git a/NonLinearSolver/materialLaw/CMakeLists.txt b/NonLinearSolver/materialLaw/CMakeLists.txt
index 0c77be16c081e3d19fc8e450ad3a015a30d48ab4..5651734146888ea29df40647568e0b8559df43e4 100644
--- a/NonLinearSolver/materialLaw/CMakeLists.txt
+++ b/NonLinearSolver/materialLaw/CMakeLists.txt
@@ -48,6 +48,7 @@ set(SRC
mlawLinearElecTherMech.cpp
mlawAnIsotropicElecTherMech.cpp
mlawSMP.cpp
+ mlawPhenomenologicalSMP.cpp
mlawEOS.cpp
mlawViscoelastic.cpp
mlawHyperelastic.cpp
diff --git a/NonLinearSolver/materialLaw/mlaw.h b/NonLinearSolver/materialLaw/mlaw.h
index 1712d2641d87424164bd7e0c333a1541debf2439..710035148a0cc004f266b13526fd7a7bf8be4158 100644
--- a/NonLinearSolver/materialLaw/mlaw.h
+++ b/NonLinearSolver/materialLaw/mlaw.h
@@ -27,7 +27,7 @@ class materialLaw{
cohesiveLaw, fracture, nonLinearShell, j2linear, viscoelastic, EOS, transverseIsotropic, transverseIsoCurvature,
transverseIsoYarnB, Anisotropic, AnisotropicStoch, nonLocalDamage, vumat,
numeric, secondOrderElastic, j2smallstrain,nonLocalDamageGurson,nonLocalDamageJ2Hyper, nonLocalDamageIsotropicElasticity,
- LinearThermoMechanics,J2ThermoMechanics,SMP,LinearElecTherMech,AnIsotropicElecTherMech,
+ LinearThermoMechanics,J2ThermoMechanics,SMP,PhenomenologicalSMP,LinearElecTherMech,AnIsotropicElecTherMech,
hyperelastic, powerYieldLaw, powerYieldLawWithFailure,nonLocalDamagePowerYieldHyper,
localDamagePowerYieldHyperWithFailure,nonLocalDamagePowerYieldHyperWithFailure,ElecSMP,
ThermalConducter,AnIsotropicTherMech, localDamageJ2Hyper,linearElastic};
diff --git a/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.cpp b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..13d191ac86018dac4cc7326c927dc7d8b321b8f0
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.cpp
@@ -0,0 +1,1156 @@
+//
+// C++ Interface: material law
+//
+// Description: SMP law
+//
+//
+// Author: <M. Pareja>, (C) 2017
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#include "mlawPhenomenologicalSMP.h"
+#include <math.h>
+#include <iostream>
+#include "MInterfaceElement.h"
+#include "ipPhenomenologicalSMP.h"
+#include "STensorOperations.h"
+//#include "STensor63.h"
+
+#include <vector>
+
+
+mlawPhenomenologicalSMP::mlawPhenomenologicalSMP(const int num,const double rho,const double alpha, const double beta, const double gamma,const double t0,
+ const double Kxg,const double Kyg, const double Kzg, const double Kxr, const double Kyr, const double Kzr, const double w, const double TaylorQuiney, const double c, const double cp, const double cg, const double cr, const double Tg, const double DeltaTg,const double Eg, const double Er, const double nug, const double nur, const double alphag, const double alphar,const double zg0,const double tau_y0):materialLaw(num,true),_alpha(alpha),_beta(beta),_gamma(gamma),_Kxg(Kxg),_Kyg(Kyg),_Kzg(Kzg), _Kxr(Kxr),_Kyr(Kyr),_Kzr(Kzr),_t0(t0),_order(-1),_w(w),_TaylorQuiney(TaylorQuiney), _c(c), _cp(cp),_cg(cg),_cr(cr),_Tg(Tg),_DeltaTg(DeltaTg),_Eg(Eg),_Er(Er),_Gg(Eg/2./(1.+nug)),_Gr(Er/2./(1.+nur)),_nug(nug),_nur(nur),_lambdag((Eg*nug)/(1.+nug)/(1.-2.*nug)),_lambdar((Er*nur)/(1.+nur)/(1.-2.*nur)),_Kg(Eg/3./(1.-2.*nug)),_Kr(Er/3./(1.-2.*nur)),_alphag(alphag),_alphar(alphar),_zg0(zg0),_tau_y0(tau_y0),_I2(1.),_I4(1.,1.){
+
+ double _Gg2=_Gg+_Gg;
+ _C0g*=0.;
+ _C0g(0,0,0,0) = _lambdag + _Gg2;
+ _C0g(1,1,0,0) = _lambdag;
+ _C0g(2,2,0,0) = _lambdag;
+ _C0g(0,0,1,1) = _lambdag;
+ _C0g(1,1,1,1) = _lambdag + _Gg2;
+ _C0g(2,2,1,1) = _lambdag;
+ _C0g(0,0,2,2) = _lambdag;
+ _C0g(1,1,2,2) = _lambdag;
+ _C0g(2,2,2,2) = _lambdag + _Gg2;
+
+ _C0g(1,0,1,0) = _Gg;
+ _C0g(2,0,2,0) = _Gg;
+ _C0g(0,1,0,1) = _Gg;
+ _C0g(2,1,2,1) = _Gg;
+ _C0g(0,2,0,2) = _Gg;
+ _C0g(1,2,1,2) = _Gg;
+
+ _C0g(0,1,1,0) = _Gg;
+ _C0g(0,2,2,0) = _Gg;
+ _C0g(1,0,0,1) = _Gg;
+ _C0g(1,2,2,1) = _Gg;
+ _C0g(2,0,0,2) = _Gg;
+ _C0g(2,1,1,2) = _Gg;
+
+ double _Gr2=_Gr+_Gr;
+ _C0r*=0.;
+ _C0r(0,0,0,0) = _lambdar + _Gr2;
+ _C0r(1,1,0,0) = _lambdar;
+ _C0r(2,2,0,0) = _lambdar;
+ _C0r(0,0,1,1) = _lambdar;
+ _C0r(1,1,1,1) = _lambdar + _Gr2;
+ _C0r(2,2,1,1) = _lambdar;
+ _C0r(0,0,2,2) = _lambdar;
+ _C0r(1,1,2,2) = _lambdar;
+ _C0r(2,2,2,2) = _lambdar + _Gr2;
+
+ _C0r(1,0,1,0) = _Gr;
+ _C0r(2,0,2,0) = _Gr;
+ _C0r(0,1,0,1) = _Gr;
+ _C0r(2,1,2,1) = _Gr;
+ _C0r(0,2,0,2) = _Gr;
+ _C0r(1,2,1,2) = _Gr;
+
+ _C0r(0,1,1,0) = _Gr;
+ _C0r(0,2,2,0) = _Gr;
+ _C0r(1,0,0,1) = _Gr;
+ _C0r(1,2,2,1) = _Gr;
+ _C0r(2,0,0,2) = _Gr;
+ _C0r(2,1,1,2) = _Gr;
+
+
+ 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++){
+ _I4dev(i,j,k,l) = 0.5*(_I2(i,k)*_I2(j,l)+_I2(i,l)*_I2(j,k)) - _I2(i,j)*_I2(k,l)/3.;
+ }
+ }
+ }
+ }
+ STensor3 R; //3x3 rotation matrix
+
+ double co1,co2,co3,si1,si2,si3;
+ double s1c2, c1c2;
+ double _pi(3.14159265359);
+ double fpi = _pi/180.;
+
+ co1 = cos(_alpha*fpi);
+ si1 = sin(_alpha*fpi);
+
+ co2 = cos(_beta*fpi);
+ si2 = sin(_beta*fpi);
+
+ co3 = cos(_gamma*fpi);
+ si3 = sin(_gamma*fpi);
+
+ s1c2 = si1*co2;
+ c1c2 = co1*co2;
+
+ R(0,0) = co3*co1 - s1c2*si3;
+ R(0,1) = co3*si1 + c1c2*si3;
+ R(0,2) = si2*si3;
+
+ R(1,0) = -si3*co1 - s1c2*co3;
+ R(1,1) = -si3*si1 + c1c2*co3;
+ R(1,2) = si2*co3;
+
+ R(2,0) = si1*si2;
+ R(2,1) = -co1*si2;
+ R(2,2) = co2;
+
+ STensor3 alphaDilatation;
+ alphaDilatation(0,0)= _alphar;
+ alphaDilatation(1,1)= _alphar;
+ alphaDilatation(2,2)= _alphar;
+ for(int i=0;i<3;i++)
+ {
+ for(int j=0;j<3;j++)
+ {
+ _alphaDilatation(i,j)=0.;
+ for(int m=0;m<3;m++)
+ {
+ for(int n=0;n<3;n++)
+ {
+ _alphaDilatation(i,j)+=R(m,i)*R(n,j)*alphaDilatation(m,n);
+ }
+ }
+ }
+ }
+ double nu = 0.29;
+ double E = 771.e6;
+ double mu = _Gg;
+ double lambda = (E*nu)/(1.+nu)/(1.-2.*nu);
+ double twicemu = mu+mu;
+
+ STensor43 K_;
+ K_*=0.;
+ K_(0,0,0,0) = lambda + twicemu;
+ K_(1,1,0,0) = lambda;
+ K_(2,2,0,0) = lambda;
+ K_(0,0,1,1) = lambda;
+ K_(1,1,1,1) = lambda + twicemu;
+ K_(2,2,1,1) = lambda;
+ K_(0,0,2,2) = lambda;
+ K_(1,1,2,2) = lambda;
+ K_(2,2,2,2) = lambda+ twicemu;
+
+ if(lambda>=1000.*mu)
+ {
+ mu = lambda + mu;
+ }
+
+ K_(1,0,1,0) = mu;
+ K_(2,0,2,0) = mu;
+ K_(0,1,0,1) = mu;
+ K_(2,1,2,1) = mu;
+ K_(0,2,0,2) = mu;
+ K_(1,2,1,2) = mu;
+
+ K_(0,1,1,0) = mu;
+ K_(0,2,2,0) = mu;
+ K_(1,0,0,1) = mu;
+ K_(1,2,2,1) = mu;
+ K_(2,0,0,2) = mu;
+ K_(2,1,1,2) = mu;
+
+
+ for(int i=0;i<3;i++)
+ {
+ for(int j=0;j<3;j++)
+ {
+ _Stiff_alphaDilatation(i,j)=0.;
+ for(int k=0;k<3;k++)
+ {
+ for(int l=0;l<3;l++)
+ {
+ _Stiff_alphaDilatation(i,j)+=K_(i,j,k,l)*_alphaDilatation(k,l);
+ }
+ }
+ }
+ }
+ STensor3 kg, kr;
+ // to be unifom in DG3D q= k' gradT with k'=-k instead of q=-kgradT
+ kg(0,0) = -_Kxg;
+ kg(1,1) = -_Kyg;
+ kg(2,2) = -_Kzg;
+ kr(0,0) = -_Kxr;
+ kr(1,1) = -_Kyr;
+ kr(2,2) = -_Kzr;
+ for(int i=0;i<3;i++)
+ {
+ for(int j=0;j<3;j++)
+ {
+ _kg(i,j)=0.;
+ for(int m=0;m<3;m++)
+ {
+ for(int n=0;n<3;n++)
+ {
+ _kg(i,j)+=R(m,i)*R(n,j)*kg(m,n);
+ }
+ }
+ }
+ }
+ for(int i=0;i<3;i++)
+ {
+ for(int j=0;j<3;j++)
+ {
+ _kr(i,j)=0.;
+ for(int m=0;m<3;m++)
+ {
+ for(int n=0;n<3;n++)
+ {
+ _kr(i,j)+=R(m,i)*R(n,j)*kr(m,n);
+ }
+ }
+ }
+ }
+}
+
+mlawPhenomenologicalSMP::mlawPhenomenologicalSMP(const mlawPhenomenologicalSMP &source) : materialLaw(source),_alpha(source._alpha),_beta(source._beta),_gamma(source._gamma)
+ ,_Kxg(source. _Kxg),_Kyg(source._Kyg),_Kzg(source._Kzg),_Kxr(source._Kxr),_Kyr(source._Kyr),_Kzr(source._Kzr), _kg(source._kg)
+ , _kr(source._kr), _t0(source._t0), _k0(source._k0),
+ _Stiff_alphaDilatation(source._Stiff_alphaDilatation),_alphaDilatation(source._alphaDilatation), _order(source._order), _w(source._w),_TaylorQuiney(source._TaylorQuiney), _c(source._c), _cp(source._cp),_cg(source._cg),_cr(source._cr),_Tg(source._Tg),_DeltaTg(source._DeltaTg),_Eg(source._Eg),_Er(source._Er),_Gg(source._Gg),_Gr(source._Gr),_nug(source._nug),_nur(source._nur),_lambdag(source._lambdag),_lambdar(source._lambdar),_Kg(source._Kg),_Kr(source._Kr),_alphag(source._alphag),_alphar(source._alphar),_zg0(source._zg0),_tau_y0(source._tau_y0),_C0r(source._C0r),_C0g(source._C0g),_I2(1.),_I4(1.,1.)
+{
+
+}
+
+double mlawPhenomenologicalSMP::soundSpeed() const
+{
+ //use glassy
+ double K = _Kg;
+ double mu = _Gg;
+ double nu = (3.*K-2.*mu)/2./(3.*K+mu);
+ double E = 2.*mu*(1.+nu);
+
+ double factornu = (1.-nu)/((1.+nu)*(1.-2.*nu));
+ return sqrt(E*factornu/_rho);
+}
+mlawPhenomenologicalSMP& mlawPhenomenologicalSMP::operator=(const materialLaw &source)
+{
+ materialLaw::operator=(source);
+ const mlawPhenomenologicalSMP* src =static_cast<const mlawPhenomenologicalSMP*>(&source);
+ _t0=src->_t0;
+ _Stiff_alphaDilatation=src->_Stiff_alphaDilatation;
+ _alphaDilatation=src->_alphaDilatation;
+ _order=src->_order;
+ _w=src->_w;
+ _TaylorQuiney=src->_TaylorQuiney;
+ _c=src->_c;
+ _cp=src->_cp;
+ _cg=src->_cg;
+ _cr=src->_cr;
+ _Tg=src->_Tg;
+ _DeltaTg=src->_DeltaTg;
+ _Kg=src->_Kg;
+ _Kr=src->_Kr;
+ _Gg=src->_Gg;
+ _Gr=src->_Gr;
+ _Eg=src->_Eg;
+ _Er=src->_Er;
+ _nug=src->_nug;
+ _nur=src->_nur;
+ _lambdag=src->_lambdag;
+ _lambdar=src->_lambdar;
+ _alphag=src->_alphag;
+ _alphar=src->_alphar;
+ _zg0=src->_zg0;
+ _tau_y0=src->_tau_y0;
+ _alpha = src->_alpha;
+ _beta = src->_beta;
+ _gamma = src->_gamma;
+ _Kxg = src->_Kxg;
+ _Kyg = src->_Kyg;
+ _Kzg = src->_Kzg;
+ _Kxr = src->_Kxr;
+ _Kyr = src->_Kyr;
+ _Kzr = src->_Kzr;
+ _kg = src->_kg;
+ _kr = src->_kr;
+ _k0 = src->_k0;
+ _C0r = src->_C0r;
+ _C0g = src->_C0g;
+}
+
+double mlawPhenomenologicalSMP::getExtraDofStoredEnergyPerUnitField(double zg) const
+{
+ double specificheat=zg*_cg+(1.-zg)*_cr;
+ return specificheat;
+}
+
+void mlawPhenomenologicalSMP::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+ bool inter=true;
+ const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+ if(iele==NULL) inter=false;
+ IPVariable* ipvi = new IPPhenomenologicalSMP(getzg0());
+ IPVariable* ipv1 = new IPPhenomenologicalSMP(getzg0());
+ IPVariable* ipv2 = new IPPhenomenologicalSMP(getzg0());
+ if(ips != NULL) delete ips;
+ ips = new IP3State(state_,ipvi,ipv1,ipv2);
+}
+
+
+void mlawPhenomenologicalSMP::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P,
+ const IPPhenomenologicalSMP *q0, IPPhenomenologicalSMP *q1,STensor43 &Tangent,const bool stiff) const
+{
+
+
+}
+
+
+void mlawPhenomenologicalSMP::constitutive( const STensor3& F0, // initial deformation gradient (input @ time n)
+ const STensor3& Fn, // updated deformation gradient (input @ time n+1)
+ // updated 1st Piola-Kirchhoff stress tensor (output)
+ STensor3 &P, // contains the initial values on input
+ const IPPhenomenologicalSMP *q0, // array of initial internal variable
+ IPPhenomenologicalSMP *q1, // updated array of internal variable (in ipvcur on output),
+ STensor43 &Tangent, // constitutive tangents (output)
+ double T0,
+ double T,
+ const SVector3 &gradT,
+ SVector3 &fluxT,
+ STensor3 &dPdT,
+ STensor3 &dqdgradT,
+ SVector3 &dqdT,
+ STensor33 &dqdF,
+ const bool stiff,// if true compute the tangents
+ double &wth,
+ double &dwthdT,
+ STensor3 &dwthdF
+ ) const
+{
+ wth=0.;dwthdT=0.; STensorOperation::zero(dwthdF); STensorOperation::zero(dqdF); STensorOperation::zero(dqdT); STensorOperation::zero(dqdgradT); STensorOperation::zero(dPdT); STensorOperation::zero(fluxT);
+ static STensor3 Pg,Pr;
+ STensorOperation::zero(Pg);
+ STensorOperation::zero(Pr);
+ static STensor43 Tangentr,Tangentg;
+ static STensor3 dPrdT,dPgdT;
+ STensorOperation::zero(Tangentg); STensorOperation::zero(dPgdT);
+ STensorOperation::zero(Tangentr); STensorOperation::zero(dPrdT);
+
+ double &zg=q1->getRefTozg();
+ double dgammataueq=0.;
+ double dh;
+ double dcpdt=0.;
+ dh = getTimeStep();
+
+ if(T>=_Tg+_DeltaTg)
+ {
+ constitutive1(F0,Fn,q0,q1,Tangentr,dPrdT,T,stiff,Pr);
+ }
+ if((T<=_Tg+_DeltaTg)and(T>=_Tg-_DeltaTg))
+ {
+ constitutive2(F0,Fn,q0,q1,Tangentr,Tangentg,dPrdT,dPgdT,T0,T,stiff,Pr,Pg,dh,dgammataueq);
+ }
+ if(T<=_Tg-_DeltaTg)
+ {
+ constitutive3(F0,Fn,q0,q1,Tangentg,dPgdT,T,stiff,Pg,dh,dgammataueq);
+ }
+
+ P =Pr;
+ P+=Pg;
+
+ double specificheat=zg*_cg+(1.-zg)*_cr;
+
+ if(dh>0.)
+ {
+ wth=-specificheat*(T-T0)/dh;// the - sign become from the term of thermal source
+ wth+=(dgammataueq)*_TaylorQuiney/dh;
+ }
+ else
+ wth=0.;
+
+ /* if(dh>0.)
+ {
+ dwthdT=-dcpdt*(T-T0)/dh-specificheat/dh;
+ dwthdT+=(dtauepsilon1dt)*_TaylorQuiney/dh;
+ dwthdT+=(dtauepsilon2dt)*_TaylorQuiney/dh;
+
+
+ for(int K = 0; K< 3; K++){
+ for(int i = 0; i< 3; i++){
+ dwthdF(K,i)=-dcpdF(K,i)*(T-T0)/dh;
+ dwthdF(K,i)+=(dtauepsilon1dF(K,i))*_TaylorQuiney/dh;
+ }
+ }
+ }
+ else*/
+ //{
+ dwthdT=0.;
+ STensorOperation::zero(dwthdF);
+ //}
+
+ static STensor3 Fninv;
+ STensorOperation::inverseSTensor3(Fn, Fninv);
+ double J= Fn.determinant();
+ static STensor3 _Kref,_k;
+ STensorOperation::zero(_Kref);STensorOperation::zero(_k);
+ _k=zg*_kg;
+ _k+=(1-zg)*_kr;
+ for(int K = 0; K< 3; K++)
+ {
+ for(int L = 0; L< 3; L++)
+ {
+ for(int i = 0; i< 3; i++)
+ {
+ for(int j = 0; j< 3; j++)
+ {
+ _Kref(K,L) += Fninv(K,i)*_k(i,j)*Fninv(L,j);
+ }
+ }
+ }
+ }
+ STensorOperation::multSTensor3SVector3(_Kref, gradT, fluxT);
+ fluxT*=J;
+
+ static STensor3 dkdT;
+ STensorOperation::zero(dkdT);
+ if(stiff)
+ {
+ dPdT=dPgdT;
+ dPdT+=dPrdT;
+ Tangent=Tangentg;
+ Tangent+=Tangentr;
+
+ STensorOperation::multSTensor3SVector3(dkdT, gradT, dqdT);
+ dqdT*=J;
+
+ dqdgradT = _Kref;
+ dqdgradT*=J;
+
+
+ for(int K = 0; K< 3; K++)
+ {
+ for(int m = 0; m< 3; m++)
+ {
+ for(int N = 0; N< 3; N++)
+ {
+ for(int L = 0; L< 3; L++)
+ {
+ dqdF(K,m,N) -= Fninv(K,m)*_Kref(N,L)*gradT(L)*J;
+ dqdF(K,m,N) -= _Kref(K,N)*Fninv(L,m)*gradT(L)*J;
+ dqdF(K,m,N) += _Kref(K,L)*gradT(L)*Fninv(N,m)*J;
+ }
+ }
+ }
+ }
+ }
+ // Tangent.print("tangent=");
+}
+
+void mlawPhenomenologicalSMP::constitutive1(const STensor3& F0, const STensor3& Fn, const IPPhenomenologicalSMP *q0, IPPhenomenologicalSMP *q1,STensor43& Tangentr, STensor3& dPrdT, double T, bool stiff, STensor3& Pr) const
+{
+
+ const STensor3 &Fp0=q0->getConstRefToFp();
+ STensor3 &Fp=q1->getRefToFp();
+ STensor3 &Ff=q1->getRefToFf();
+ STensor3 &Fpg=q1->getRefToFpg();
+ double &zg=q1->getRefTozg();
+
+ //Rubbery phase//
+
+ //Internal variables evoluton
+
+ zg=0.;
+ Fp=Fp0;
+ STensor3 subs1=Fn;
+ subs1-=_I2;
+ Ff=subs1;
+ Ff*=_c;
+ Ff+=_I2;
+ STensorOperation::unity(Fpg);
+
+ //Resulting variables
+
+ static STensor3 Fer, Eer,tauer;
+ static STensor43 dlnCer;
+ static STensor63 ddlnCer;
+ STensorOperation::zero(dlnCer);STensorOperation::zero(ddlnCer);STensorOperation::zero(Eer);
+
+ computePr(Fn,Fp,Fer,Eer,tauer,dlnCer,ddlnCer,zg,T,Pr);
+
+ //Tangent operators//
+ if(stiff){
+ tangentrubbery(Fer,Fp,tauer,dlnCer,ddlnCer,Eer,T,zg,dPrdT,Tangentr);
+ }
+
+}
+
+void mlawPhenomenologicalSMP::constitutive2(const STensor3& F0, const STensor3& Fn, const IPPhenomenologicalSMP *q0, IPPhenomenologicalSMP *q1, STensor43& Tangentr, STensor43& Tangentg, STensor3& dPrdT, STensor3& dPgdT, double T0, double T, bool stiff, STensor3& Pr, STensor3& Pg, const double dh, double& dgammataueq) const
+{
+
+ const STensor3 &Fp0=q0->getConstRefToFp();
+ STensor3 &Fp=q1->getRefToFp();
+ const STensor3 &Ff0=q0->getConstRefToFf();
+ STensor3 &Ff=q1->getRefToFf();
+ const STensor3 &Fpg0=q0->getConstRefToFpg();
+ STensor3 &Fpg=q1->getRefToFpg();
+ double &zg=q1->getRefTozg();
+
+ //Mixed phase//
+
+ zg=1./(1.+exp(2.*_w*(T-_Tg)));
+ Fp=Fp0;
+ Ff=Ff0;
+
+//printf("T0= %f, T= %f\n",T0,T);
+
+ if(T>T0){
+ Fpg=Fpg0;
+ static STensor3 Feg, Eeg,taueg;
+ static STensor43 dlnCeg;
+ static STensor63 ddlnCeg;
+
+ computePg(Fn,Ff,Fpg,Feg,Eeg,taueg,dlnCeg,ddlnCeg,zg,T,Pg);
+
+ double Eegpreq=0;
+ static STensor3 N;
+ N = _I2;
+ double Delta_gamma=0.;
+ double H_gamma=0.; //If tau_y evolves with Delta_gamma then it is not 0
+ double H_T=0.; //If tau_y evolves with T then it is not 0
+ static STensor43 dexpDgammaN;
+ dexpDgammaN=_I4;
+
+ tangentglassy(Fn,q0,q1,dexpDgammaN, Feg, dlnCeg,Eeg, N, Eegpreq, Delta_gamma, zg, H_gamma, H_T, T,Feg,taueg,dlnCeg,ddlnCeg,dPgdT,Tangentg);
+ }
+ if(T<=T0){
+ double tau_y=zg*_tau_y0;
+ predictorCorrector(F0,Fn,q0,q1,tau_y,Pg,dh,T,dgammataueq,stiff,dPgdT, Tangentg);
+ }
+
+ //Resulting variables
+
+ static STensor3 Fer, Eer,tauer;
+ static STensor43 dlnCer;
+ static STensor63 ddlnCer;
+ STensorOperation::zero(dlnCer);STensorOperation::zero(ddlnCer);STensorOperation::zero(Eer);
+
+ computePr(Fn,Fp,Fer,Eer,tauer,dlnCer,ddlnCer,zg,T,Pr);
+
+ //Tangent operators//
+ if(stiff){
+ tangentrubbery(Fer,Fp,tauer,dlnCer,ddlnCer,Eer,T,zg,dPrdT,Tangentr);
+ }
+
+}
+
+void mlawPhenomenologicalSMP::constitutive3(const STensor3& F0, const STensor3& Fn, const IPPhenomenologicalSMP *q0, IPPhenomenologicalSMP *q1,STensor43& Tangentg, STensor3& dPgdT, double T, bool stiff, STensor3& Pg, const double dh, double& dgammataueq) const
+{
+
+ const STensor3 &Fp0=q0->getConstRefToFp();
+ STensor3 &Fp=q1->getRefToFp();
+ const STensor3 &Ff0=q0->getConstRefToFf();
+ STensor3 &Ff=q1->getRefToFf();
+ const STensor3 &Fpg0=q0->getConstRefToFpg();
+ STensor3 &Fpg=q1->getRefToFpg();
+ double &zg=q1->getRefTozg();
+
+
+ //Glassy phase//
+
+ //Internal variables evoluton
+
+ zg=1.;
+ Ff=Ff0;
+ double tau_y=_tau_y0;
+ static STensor43 dtauegdF, dFpgdF;
+ static STensor3 dtauegdT, dFpgdT;
+ static STensor63 ddlnCeg;
+ STensorOperation::zero(dtauegdF);STensorOperation::zero(dtauegdT);STensorOperation::zero(dFpgdF);STensorOperation::zero(dFpgdT);STensorOperation::zero(ddlnCeg);
+ predictorCorrector(F0,Fn,q0,q1,tau_y,Pg,dh,T,dgammataueq,stiff,dPgdT,Tangentg);
+ static STensor3 subs3=Ff;
+ subs3+=Fpg;
+ subs3-=2.*_I2;
+ Fp=_cp*subs3;
+ Fp+=_I2;
+ if(!checkDeltaEp(Fp,Fp0)){
+ Fp=Fp0;
+ }
+}
+
+bool mlawPhenomenologicalSMP::checkDeltaEp(const STensor3& Fp0, const STensor3& Fp) const
+{
+ STensor3 Cp, lnCp, Ep;
+ STensor43 dlnCp;
+ STensorOperation::zero(Cp);STensorOperation::zero(lnCp);STensorOperation::zero(Ep);STensorOperation::zero(dlnCp);
+ double Epnorm0;
+ double Epnorm1;
+ for(int i=0;i<2;i++){
+ if(i==0){
+ STensorOperation::multSTensor3FirstTranspose(Fp0,Fp0,Cp);
+ STensorOperation::logSTensor3(Cp,_order,lnCp,&dlnCp);
+ Ep=0.5*lnCp;
+ Epnorm0=Ep.norm2();
+ }
+ if(i==1){
+ STensorOperation::multSTensor3FirstTranspose(Fp,Fp,Cp);
+ STensorOperation::logSTensor3(Cp,_order,lnCp,&dlnCp);
+ Ep=0.5*lnCp;
+ Epnorm1=Ep.norm2();
+ }
+ }
+ if(Epnorm1>Epnorm0) return true;
+ else return false;
+}
+
+
+void mlawPhenomenologicalSMP::predictorCorrector(const STensor3& F0, const STensor3& Fn, const IPPhenomenologicalSMP *q0, IPPhenomenologicalSMP *q1, double tau_y,STensor3& Pg, const double dh,double T,double& dgammataueq, bool stiff,STensor3& dPgdT, STensor43& Tangentg) const
+{
+ const double &gamma0=q0->getConstRefToEquivalentPlasticDefo();
+ double &gamma1=q1->getRefToEquivalentPlasticDefo();
+
+ gamma1=gamma0;
+
+ STensor3 &Ff=q1->getRefToFf();
+ const STensor3 &Fpg0=q0->getConstRefToFpg();
+ STensor3 &Fpg=q1->getRefToFpg();
+ double &zg=q1->getRefTozg();
+
+ static STensor3 Fegpr, Eegpr,tauegpr;
+ static STensor43 dlnCegpr;
+ static STensor63 ddlnCegpr;
+
+ computePg(Fn,Ff,Fpg0,Fegpr,Eegpr,tauegpr,dlnCegpr,ddlnCegpr,zg,T,Pg);
+
+ double taueqpr=sqrt(1.5*tauegpr.dev().dotprod());
+
+ double Eegpreq;
+ Eegpreq=sqrt((2./3.)*Eegpr.dev().dotprod());
+
+ double tolNR=1.e-6;
+
+ static STensor3 Feg, Eeg, devtaueg, taueg;
+ static STensor43 dlnCeg;
+ static STensor63 ddlnCeg;
+ Feg=Fegpr;
+ Eeg=Eegpr;
+ taueg=tauegpr;
+ dlnCeg=dlnCegpr;
+ ddlnCeg=ddlnCegpr;
+
+ static STensor3 N;
+ N = Eegpr.dev();
+ N*= 1/Eegpreq;
+
+ double f=taueqpr-tau_y;
+
+ double Delta_gamma=0.;
+ double H_gamma=0.; //If tau_y evolves with Delta_gamma then it is not 0
+ double H_T=0.; //If tau_y evolves with T then it is not 0
+
+ static STensor43 dexpDgammaN;
+
+ if(f>0.){
+
+ //Newton-Raphson//
+
+ //Initialization of variables
+
+ int ite = 0, maxite = 1000;
+ static bool flag=false;
+
+ double eps_f=fabs(f/taueqpr);
+ double j1=1./((3.*zg*_Gg)+H_gamma);
+ double taueq;
+
+ //Iterative process
+
+ while(eps_f>tolNR or ite<1)
+ {
+ //Compute Jacobian: in this case it's constant (outside of the iterative process)
+ double deps=j1*f;
+ gamma1+=deps;
+ Eeg-=N*deps;
+ devtaueg = Eeg.dev();
+ devtaueg *= (2.*_Gg*zg);
+
+ taueq=sqrt(1.5*devtaueg.dotprod());
+
+ f=taueq-tau_y;
+ eps_f=fabs(f/taueqpr);
+ if(ite > maxite)
+ {
+ Msg::Error("SMP itereration %d, residual %f,time step %f ", ite, eps_f ,dh);
+ break;
+ }
+ ite++;
+ }
+ Delta_gamma=gamma1-gamma0;
+
+ static STensor3 DepsN;
+ DepsN = N;
+ DepsN *= Delta_gamma;
+ static STensor3 expDgammaN;
+ STensorOperation::expSTensor3(DepsN,_order,expDgammaN,&dexpDgammaN);
+
+ //Update of variables
+ STensorOperation::multSTensor3(expDgammaN,Fpg0,Fpg);
+ computePg(Fn,Ff,Fpg,Feg,Eeg,tauegpr,dlnCeg,ddlnCeg,zg,T,Pg);
+
+ //Plastic dissipation term
+ dgammataueq=taueq*Delta_gamma;
+ }
+ else{gamma1=gamma0;}
+ if(stiff){
+ //Tangent operators glassy phase
+ tangentglassy(Fn,q0,q1,dexpDgammaN, Fegpr, dlnCegpr,Eegpr, N, Eegpreq, Delta_gamma, zg, H_gamma, H_T, T,Feg,taueg,dlnCeg,ddlnCeg,dPgdT,Tangentg);
+ }
+
+}
+
+
+void mlawPhenomenologicalSMP::tangentrubbery(const STensor3& Fer,const STensor3& Fp,const STensor3& tauer,const STensor43& dlnCer,const STensor63& ddlnCer,const STensor3& Eer,const double T,const double zg, STensor3& dPrdT, STensor43& Tangentr) const {
+
+ static STensor3 Fpinv;
+ STensorOperation::inverseSTensor3(Fp,Fpinv);
+
+ static STensor43 dFerdF,dtauerdF;
+ static STensor3 dtauerdT;
+ STensorOperation::zero(dFerdF);STensorOperation::zero(dtauerdF);STensorOperation::zero(dtauerdT);
+
+ static double dzgdT=0.;
+ if((T<=_Tg+_DeltaTg)and(T>=_Tg-_DeltaTg)){
+ dzgdT=-2.*_w*exp(2.*_w*(T-_Tg));
+ dzgdT*=1./pow((1.+2.*_w*exp(2.*_w*(T-_Tg))),2);
+ }
+
+dFerdF *= 0.;
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ for (int k=0; k<3; k++){
+ dFerdF(i,j,i,k) += Fpinv(k,j);
+ }
+ }
+ }
+
+static STensor63 DLDF;
+ if (_order != 1){
+ 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 p=0; p<3; p++){
+ for (int q=0; q<3; q++){
+ DLDF(i,j,k,l,p,q) = 0.;
+ for (int r=0; r<3; r++){
+ for (int s=0; s<3; s++){
+ for (int a=0; a<3; a++){
+ DLDF(i,j,k,l,p,q) += ddlnCer(i,j,k,l,r,s)*2.*Fer(a,r)*dFerdF(a,s,p,q);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ else{
+ DLDF *= 0.;
+ }
+
+static STensor43 EerToF, dtauerdEer;
+ dtauerdEer=(_C0r);
+ dtauerdEer*=(1.-zg);
+ 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++){
+ EerToF(i,j,k,l) = 0.;
+ for (int p=0; p<3; p++){
+ for (int q=0; q<3; q++){
+ EerToF(i,j,k,l) += dlnCer(i,j,p,q)*Fer(k,p)*Fpinv(l,q);
+ }
+ }
+ }
+ }
+ }
+ }
+
+STensorOperation::multSTensor43(dtauerdEer,EerToF,dtauerdF);
+
+
+dtauerdT *= 0.;
+ for( int E=0; E<3; E++){
+ // dtauerdT(E,E) += -(1.-zg)*3.*_Kr*_alphar;
+ for( int Z=0; Z<3; Z++){
+ // dtauerdT(E,Z) += dzgdT*3.*_Kr*_alphar*(T-_t0)*_I2(E,Z);
+ for( int N=0; N<3; N++){
+ for( int O=0; O<3; O++){
+ dtauerdT(E,Z) -= dzgdT*_C0r(E,Z,N,O)*Eer(N,O);
+ }
+ }
+ }
+ }
+
+static STensor43 DSDF; // S = corKir:L
+static STensor3 DSDT;
+STensorOperation::zero(DSDT);
+
+ 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++){
+ DSDF(i,j,k,l) = 0.;
+ DSDT(i,j)+= dtauerdT(k,l)*dlnCer(k,l,i,j);
+ for (int r=0; r<3; r++){
+ for (int s=0; s<3; s++){
+ DSDF(i,j,k,l) += dtauerdF(r,s,k,l)*dlnCer(r,s,i,j) ;
+ if (_order != 1){
+ DSDF(i,j,k,l) += tauer(r,s)*DLDF(r,s,i,j,k,l);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+
+static STensor3 Ser;
+STensorOperation::multSTensor3STensor43(tauer,dlnCer,Ser);
+
+// compute mechanical & thermal tangent
+ Tangentr *= 0.;
+ dPrdT *= 0.;
+ 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++){
+ dPrdT(i,j) += (Fer(i,k)*DSDT(k,l)*Fpinv(j,l));
+ for (int p=0; p<3; p++){
+ for (int q=0; q<3; q++){
+ Tangentr(i,j,p,q) += (dFerdF(i,k,p,q)*Ser(k,l)*Fpinv(j,l));
+ Tangentr(i,j,p,q) += (Fer(i,k)*DSDF(k,l,p,q)*Fpinv(j,l));
+ }
+ }
+ }
+ }
+ }
+ }
+
+}
+
+
+void mlawPhenomenologicalSMP::tangentglassy(const STensor3& Fn,const IPPhenomenologicalSMP *q0, IPPhenomenologicalSMP *q1, const STensor43& dexpDelta_gammaNpr,const STensor3& Fegpr, STensor43& dlnCegpr,const STensor3& Eegpr,const STensor3& Npr,const double Eegpreq,const double Deltagamma,const double zg,const double H_gamma,const double H_T,const double T,const STensor3& Feg,const STensor3& taueg,const STensor43& dlnCeg,const STensor63& ddlnCeg, STensor3& dPgdT, STensor43& Tangentg) const {
+
+ const STensor3& Fpg0 = q0->getConstRefToFpg();
+ static STensor3 Fpg0inv;
+ STensorOperation::inverseSTensor3(Fpg0,Fpg0inv);
+
+ STensor3 &Ff=q1->getRefToFf();
+ STensor3 &Fpg=q1->getRefToFpg();
+ STensor3 Ffinv, Fpginv;
+ STensorOperation::inverseSTensor3(Ff,Ffinv);STensorOperation::inverseSTensor3(Fpg,Fpginv);
+
+ static STensor43 dtauegdF, dFpgdF;
+ static STensor3 dtauegdT, dFpgdT;
+
+ STensorOperation::zero(dtauegdF);STensorOperation::zero(dtauegdT);STensorOperation::zero(dFpgdF);STensorOperation::zero(dFpgdT);
+
+ static STensor43 dFpgdEegpr, dtauegdEegpr;
+ STensorOperation::zero(dFpgdEegpr);STensorOperation::zero(dtauegdEegpr);
+
+ dtauegdEegpr = _C0g;
+ dtauegdEegpr *= zg;
+
+ static double dzgdT=0.;
+ if((T<=_Tg+_DeltaTg)and(T>=_Tg-_DeltaTg)){
+ dzgdT=-2.*_w*exp(2.*_w*(T-_Tg));
+ dzgdT*=1./pow((1.+2.*_w*exp(2.*_w*(T-_Tg))),2);
+ }
+
+ dtauegdT *= 0.;
+ for( int E=0; E<3; E++){
+ // dtauegdT(E,E) += -3.*_Kg*_alphag*zg;
+ for( int Z=0; Z<3; Z++){
+ // dtauegdT(E,Z) -= dzgdT*3.*_Kg*_alphag*(T-_t0)*_I2(E,Z);
+ for( int N=0; N<3; N++){
+ for( int O=0; O<3; O++){
+ dtauegdT(E,Z) += dzgdT*_C0g(E,Z,N,O)*Eegpr(N,O);
+ }
+ }
+ }
+ }
+
+if (Deltagamma > 0.){
+ static STensor3 DgammaDEepr;
+ DgammaDEepr = (Npr);
+ DgammaDEepr *= ((2.*_Gg*zg)/(3.*_Gg*zg+H_gamma));
+ static double dDeltagammadT=-(H_T+3.*zg*_Gg*dzgdT);
+ dDeltagammadT*=1./(3.*zg*_Gg+H_gamma);
+
+ static STensor43 BE;
+ tensprod(Npr,Npr,BE);
+ BE *= (2.*_Gg*zg/(3.*_Gg*zg+H_gamma)- 2.*Deltagamma/3./Eegpreq);
+ double val=Deltagamma/Eegpreq;
+ BE.axpy(val,_I4dev);
+
+ static STensor3 BT;
+ BT = (Npr);
+ BT*= dDeltagammadT;
+
+ // update with plasticity corKir
+ dtauegdEegpr.axpy(-2.*_Gg*zg,BE);
+ dtauegdT.daxpy(BT,-2.*_Gg*zg);
+ dtauegdT.daxpy(Npr,-2.*_Gg*dzgdT*Deltagamma);
+
+ static STensor43 DexpABE;
+ STensorOperation::multSTensor43(dexpDelta_gammaNpr,BE,DexpABE);
+ // update with plasticity Fp
+ 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 p=0; p<3; p++){
+ dFpgdEegpr(i,j,k,l) += DexpABE(i,p,k,l)*Fpg0(p,j);
+ }
+ }
+ }
+ }
+ }
+ static STensor3 DexpABT;
+ STensorOperation::multSTensor43STensor3(dexpDelta_gammaNpr,BT,DexpABT);
+ STensorOperation::multSTensor3Add(DexpABT,Fpg0,dFpgdT);
+
+ }
+
+ static STensor43 EegprToF;
+ 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++){
+ EegprToF(i,j,k,l) = 0.;
+ for (int p=0; p<3; p++){
+ for (int q=0; q<3; q++){
+ for (int r=0; r<3; r++){
+ EegprToF(i,j,k,l) += dlnCegpr(i,j,p,q)*Fegpr(k,p)*Ffinv(l,r)*Fpg0inv(r,q);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ STensorOperation::multSTensor43(dtauegdEegpr,EegprToF,dtauegdF);
+ if (Deltagamma>0.){
+ STensorOperation::multSTensor43(dFpgdEegpr,EegprToF,dFpgdF);
+ }
+ else{
+ STensorOperation::zero(dFpgdF);
+ }
+
+ // done DcorKirDF, DcorKirDT, DFpDF, DFpDT
+
+ static STensor43 DinvFpgdF;
+ static STensor3 DinvFpgdT;
+ STensorOperation::zero(DinvFpgdF);
+ STensorOperation::zero(DinvFpgdT);
+ if (Deltagamma >0.){
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ for (int p=0; p<3; p++){
+ for (int q=0; q<3; q++){
+ DinvFpgdT(i,j) -= Fpginv(i,p)*dFpgdT(p,q)*Fpginv(q,j);
+ for (int k=0; k<3; k++){
+ for (int l=0; l<3; l++){
+ DinvFpgdF(i,j,k,l) -= Fpginv(i,p)*dFpgdF(p,q,k,l)*Fpginv(q,j);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ static STensor43 dFegdF;
+ static STensor3 dFegdT;
+ STensorOperation::zero(dFegdF);
+ STensorOperation::zero(dFegdT);
+ 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++){
+ dFegdF(i,j,i,k) += Ffinv(k,l)*Fpginv(l,j);
+ dFegdT(i,j) += Fn(i,l)*Ffinv(l,k)*DinvFpgdT(k,j);
+ if (Deltagamma> 0.){
+ for (int p=0; p<3; p++){
+ for (int r=0; r<3; r++){
+ dFegdF(i,j,k,l) += Fn(i,r)*Ffinv(r,p)*DinvFpgdF(p,j,k,l);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ static STensor63 DLDF;
+ static STensor43 DLDT;
+ STensorOperation::zero(DLDF);
+ STensorOperation::zero(DLDT);
+ 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 r=0; r<3; r++){
+ for (int s=0; s<3; s++){
+ for (int a=0; a<3; a++){
+ DLDT(i,j,k,l) += ddlnCeg(i,j,k,l,r,s)*2.*Feg(a,r)*dFegdT(a,s);
+ for (int p=0; p<3; p++){
+ for (int q=0; q<3; q++){
+ DLDF(i,j,k,l,p,q) += ddlnCeg(i,j,k,l,r,s)*2.*Feg(a,r)*dFegdF(a,s,p,q);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+
+ static STensor43 DSDF; // S = corKir:L
+ static STensor3 DSDT;
+ STensorOperation::zero(DSDF);
+ STensorOperation::zero(DSDT);
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ for (int r=0; r<3; r++){
+ for (int s=0; s<3; s++){
+ DSDT(i,j) += dtauegdT(r,s)*dlnCeg(r,s,i,j) + taueg(r,s)*DLDT(r,s,i,j);
+ for (int k=0; k<3; k++){
+ for (int l=0; l<3; l++){
+ DSDF(i,j,k,l) += dtauegdF(r,s,k,l)*dlnCeg(r,s,i,j) + taueg(r,s)*DLDF(r,s,i,j,k,l);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ static STensor3 Seg;
+ STensorOperation::multSTensor3STensor43(taueg,dlnCeg,Seg);
+
+ // compute mechanical tangent
+ STensorOperation::zero(Tangentg);
+ STensorOperation::zero(dPgdT);
+ 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 p=0; p<3; p++){
+ dPgdT(i,j) += (dFegdT(i,k)*Seg(k,l)*Fpginv(p,l)*Ffinv(j,p) + Feg(i,k)*DSDT(k,l)*Fpginv(p,l)*Ffinv(j,p)+Feg(i,k)*Seg(k,l)*DinvFpgdT(p,l)*Ffinv(j,p));
+ for (int q=0; q<3; q++){
+ for (int r=0; r<3; r++){
+ Tangentg(i,j,p,q) += (dFegdF(i,k,p,q)*Seg(k,l)*Fpginv(r,l)*Ffinv(j,r) + Feg(i,k)*DSDF(k,l,p,q)*Fpginv(r,l)*Ffinv(j,r)+Feg(i,k)*Seg(k,l)*DinvFpgdF(r,l,p,q)*Ffinv(j,r));
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+}
+
+
+void mlawPhenomenologicalSMP::computePr(const STensor3& Fn,const STensor3& Fp,STensor3& Fer,STensor3& Eer,STensor3& tauer,STensor43& dlnCer,STensor63& ddlnCer,const double zg,const double T,STensor3& Pr)const{
+
+ static STensor3 Fpinv, Cer;
+ STensorOperation::inverseSTensor3(Fp,Fpinv);
+ STensorOperation::multSTensor3(Fn,Fpinv,Fer); //Fer
+ STensorOperation::multSTensor3FirstTranspose(Fer,Fer,Cer);//Cer
+
+ STensorOperation::logSTensor3(Cer,_order,Eer,&dlnCer,&ddlnCer);
+ Eer*=0.5; //(Elastic) Logarithmic strain of the rubbery phase
+
+ static STensor3 per, devtauer,Ser;
+ devtauer=Eer.dev();
+ devtauer*=(2.*_Gr);
+ double p_bracket=_Kr*(Eer.trace());//-3.*_alphar*(T-_t0)
+ per=_I2;
+ per*=p_bracket;
+ tauer=per;
+ tauer+=devtauer;
+ tauer*=(1.-zg);//Corotational stress of the rubbery phase
+
+ /*printf("T: %f\n",T);
+ Fn.print("Fn: \n");
+ Fp.print("Fp: \n");
+ Fer.print("Fer: \n");
+ Eer.print("Eer: \n");
+ tauer.print("tauer: \n");*/
+
+ STensorOperation::multSTensor3STensor43(tauer,dlnCer,Ser);//Elastic 2nd Piola-Kirchhoff stress of the rubbery phase
+
+ static STensor3 Pr_1, FpinvT;
+ STensorOperation::transposeSTensor3(Fpinv,FpinvT);
+ STensorOperation::multSTensor3(Fer,Ser,Pr_1);
+ STensorOperation::multSTensor3(Pr_1,FpinvT,Pr); //First Piola-Kirchhoff stres of the rubbery phase
+
+}
+
+
+void mlawPhenomenologicalSMP::computePg(const STensor3& Fn,const STensor3& Ff,const STensor3& Fpg,STensor3& Feg,STensor3& Eeg,STensor3& taueg,STensor43& dlnCeg,STensor63& ddlnCeg,const double zg,const double T,STensor3& Pg)const{
+
+ static STensor3 Ffinv, Fpginv;
+ STensorOperation::inverseSTensor3(Ff,Ffinv);
+ STensorOperation::inverseSTensor3(Fpg,Fpginv);
+ static STensor3 Feg_int;
+ STensorOperation::multSTensor3(Fn, Ffinv, Feg_int);
+ STensorOperation::multSTensor3(Feg_int, Fpginv, Feg);
+
+ static STensor3 Ceg;
+ STensorOperation::multSTensor3FirstTranspose(Feg, Feg, Ceg);
+ STensorOperation::logSTensor3(Ceg,_order,Eeg,&dlnCeg,&ddlnCeg);
+ Eeg*=0.5;
+
+ static STensor3 devtaueg, peg,Seg;
+ devtaueg=Eeg.dev();
+ devtaueg*=(2.*_Gg*zg);
+ double peg_bracket=zg*_Kg*(Eeg.trace());//-3.*_alphag*(T-_t0)
+ peg=_I2;
+ peg*=peg_bracket;
+ taueg=peg;
+ taueg+=devtaueg;
+
+ /*printf("T: %f\n",T);
+ Fn.print("Fn: \n");
+ Ff.print("Ff: \n");
+ Feg.print("Feg: \n");
+ Eeg.print("Eeg: \n");
+ taueg.print("taueg: \n");*/
+
+ STensorOperation::multSTensor3STensor43(taueg,dlnCeg,Seg);
+
+ static STensor3 Pg_1, Pg_2, FpginvT, FfinvT;
+ STensorOperation::transposeSTensor3(Fpginv,FpginvT);
+ STensorOperation::transposeSTensor3(Ffinv,FfinvT);
+ STensorOperation::multSTensor3(Feg,Seg,Pg_1);
+ STensorOperation::multSTensor3(Pg_1,FpginvT,Pg_2);
+ STensorOperation::multSTensor3(Pg_2,FfinvT,Pg); //First Piola-Kirchhoff stres of the glassy phase
+
+}
diff --git a/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.h b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.h
new file mode 100644
index 0000000000000000000000000000000000000000..865b139776c4410705b6de3def2be904109d6e04
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.h
@@ -0,0 +1,154 @@
+//
+// C++ Interface: material law
+// Linear Thermo Mechanical law
+// Author: <M. Pareja>, (C) 2017
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#ifndef MLAWPHENOMENOLOGICALSMP_H_
+#define MLAWPHENOMENOLOGICALSMP_H_
+#include "../internalPoints/ipPhenomenologicalSMP.h"
+#include "mlaw.h"
+#include "STensor3.h"
+#include "STensor43.h"
+#include "STensor63.h"
+#include "ipSMP.h"
+
+
+class mlawPhenomenologicalSMP : public materialLaw
+{
+ protected:
+ double _t0;
+ STensor3 _I2;
+ STensor43 _I4,_I4dev;
+ STensor3 _alphaDilatation;
+ STensor3 _Stiff_alphaDilatation;
+ int _order;
+ double _alpha; // parameter of anisotropy
+ double _beta; // parameter of anisotropy
+ double _gamma; // parameter of anisotropy
+ double _rho;
+ double _Kxg;
+ double _Kyg;
+ double _Kzg;
+ double _Kxr;
+ double _Kyr;
+ double _Kzr;
+ STensor3 _k0;
+ STensor3 _kg;
+ STensor3 _kr;
+ double _w, _c, _cp, _cg, _cr, _TaylorQuiney, _Tg, _DeltaTg, _Kg, _Kr,_Eg,_Er,_nug,_nur,_lambdag,_lambdar,_Gg, _Gr,_alphag,_alphar,_zg0,_tau_y0;
+
+ STensor43 _C0r,_C0g;
+
+ public:
+ mlawPhenomenologicalSMP(const int num, const double rho, const double alpha, const double beta, const double gamma ,const double t0, const double Kxg,const double Kyg, const double Kzg, const double Kxr, const double Kyr, const double Kzr,const double w, const double TaylorQuiney, const double c, const double cp, const double cg, const double cr,const double Tg,const double DeltaTg, const double Eg, const double Er, const double nug, const double nur,const double alphag, const double alphar, const double zg0, const double tau_y0);
+#ifndef SWIG
+ mlawPhenomenologicalSMP(const mlawPhenomenologicalSMP &source);
+ mlawPhenomenologicalSMP& operator=(const materialLaw &source);
+ virtual ~mlawPhenomenologicalSMP(){};
+ virtual materialLaw* clone() const {return new mlawPhenomenologicalSMP(*this);};
+ virtual bool withEnergyDissipation() const {return false;};
+ // function of materialLaw
+ virtual matname getType() const{return materialLaw::PhenomenologicalSMP;}
+ virtual void createIPState(IPStateBase* &ips,const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt =0) const;
+
+
+public:
+ void setStrainOrder(const int order){_order = order;};
+ double getInitialTemperature() const {return _t0;};
+ double getExtraDofStoredEnergyPerUnitField(double zg) const;
+ virtual double getzg0() const { return _zg0;};
+ //virtual double getInitialExtraDofStoredEnergyPerUnitField() const {getExtraDofStoredEnergyPerUnitField(zg);}
+ virtual void constitutive(
+ const STensor3& F0, // initial deformation gradient (input @ time n)
+ const STensor3& Fn, // updated deformation gradient (input @ time n+1)
+ STensor3 &P, // updated 1st Piola-Kirchhoff stress tensor (output)
+ const IPPhenomenologicalSMP *q0, // array of initial internal variable
+ IPPhenomenologicalSMP *q1, // updated array of internal variable (in ipvcur on output),
+ STensor43 &Tangent, // constitutive tangents (output)
+ const bool stiff
+ ) const;
+ virtual void constitutive(
+ const STensor3& F0, // initial deformation gradient (input @ time n)
+ const STensor3& Fn, // updated deformation gradient (input @ time n+1)
+ STensor3 &P, // updated 1st Piola-Kirchhoff stress tensor (output)
+ // contains the initial values on input
+ const IPPhenomenologicalSMP *q0, // array of initial internal variable
+ IPPhenomenologicalSMP *q1, // updated array of internal variable (in ipvcur on output),
+ STensor43 &Tangent, // constitutive tangents (output)
+ double T0,
+ double T,
+ const SVector3 &gradT,
+ SVector3 &fluxT,
+ STensor3 &dPdT,
+ STensor3 &dqdgradT,
+ SVector3 &dqdT,
+ STensor33 &dqdF,
+ const bool stiff,
+ double &wth,
+ double &dwthdt,
+ STensor3 &dwthdf
+ ) const;
+ void constitutive1(
+ const STensor3& F0,
+ const STensor3& Fn,
+ const IPPhenomenologicalSMP *q0,
+ IPPhenomenologicalSMP *q1,
+ STensor43& Tangentr,
+ STensor3& dPrdT,
+ double T,
+ bool stiff,
+ STensor3& Pr
+ ) const;
+ void constitutive2(
+ const STensor3& F0,
+ const STensor3& Fn,
+ const IPPhenomenologicalSMP *q0,
+ IPPhenomenologicalSMP *q1,
+ STensor43& Tangentr,
+ STensor43& Tangentg,
+ STensor3& dPrdT,
+ STensor3& dPgdT,
+ double T0,
+ double T,
+ bool stiff,
+ STensor3& Pr,
+ STensor3& Pg,
+ const double dh,
+ double& dgammataueq
+ ) const;
+ void constitutive3(
+ const STensor3& F0,
+ const STensor3& Fn,
+ const IPPhenomenologicalSMP *q0,
+ IPPhenomenologicalSMP *q1,
+ STensor43& Tangentg,
+ STensor3& dPgdT,
+ double T,
+ bool stiff,
+ STensor3& Pg,
+ const double dh,
+ double& dgammataueq
+ ) const;
+ const STensor3& getInitialConductivityTensor() const { return _k0;};
+ virtual const STensor3& getStiff_alphaDilatation()const { return _Stiff_alphaDilatation;};
+
+ virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}; // this law is initialized so nothing to do
+ virtual double soundSpeed() const; // default but you can redefine it for your case
+
+ protected:
+
+ void predictorCorrector(const STensor3& F0, const STensor3& Fn, const IPPhenomenologicalSMP *q0, IPPhenomenologicalSMP *q1, double tau_y,STensor3& Pg, const double dh,double T,double& dgammataueq,bool stiff,STensor3& dPgdT,STensor43& Tangentg) const;
+ bool checkDeltaEp(const STensor3& Fp0, const STensor3& Fp) const;
+ void tangentrubbery(const STensor3& Fer,const STensor3& Fp,const STensor3& tauer,const STensor43& dlnCer,const STensor63& ddlnCer,const STensor3& Eer,const double T,const double zg, STensor3& dPrdT, STensor43& Tangentr) const;
+ void tangentglassy(const STensor3& Fn,const IPPhenomenologicalSMP *q0,IPPhenomenologicalSMP *q1, const STensor43& dexpDelta_gammaNpr,const STensor3& Fegpr, STensor43& dlnCegpr,const STensor3& Eegpr,const STensor3& Npr,const double Eegpreq,const double Deltagamma,const double zg,const double H_gamma,const double H_T,const double T,const STensor3& Feg,const STensor3& taueg,const STensor43& dlnCeg,const STensor63& ddlnCeg, STensor3& dPgdT, STensor43& Tangentg) const;
+ //void plasticitytangent(const IPPhenomenologicalSMP *q0,STensor43& dtauegdF, STensor3& dtauegdT, STensor43& dFpgdF, STensor3& dFpgdT,const STensor43& dexpDelta_gammaNpr,const STensor3& Fegpr, STensor43& dlnCegpr,const STensor3& Eegpr,const STensor3& Npr,const double Eegpreq,const double Deltagamma,const double zg,const double H_gamma,const double H_T,const double T) const;
+void computePr(const STensor3& Fn,const STensor3& Fp,STensor3& Fer,STensor3& Eer,STensor3& tauer,STensor43& dlnCer,STensor63& ddlnCer,const double zg,const double T,STensor3& Pr)const;
+void computePg(const STensor3& Fn,const STensor3& Ff,const STensor3& Fpg,STensor3& Feg,STensor3& Eeg,STensor3& taueg,STensor43& dlnCeg,STensor63& ddlnCeg,const double zg,const double T,STensor3& Pg)const;
+
+ #endif // SWIG
+};
+
+#endif // MLAWPHENOMENOLOGICALSMP_H_
diff --git a/NonLinearSolver/materialLaw/mlawSMP.cpp b/NonLinearSolver/materialLaw/mlawSMP.cpp
index b3d5dbfa5acf1a945e77b97e7f27b162bdb3533d..4482238a74be81f19eb984bec9ad21d90591d817 100644
--- a/NonLinearSolver/materialLaw/mlawSMP.cpp
+++ b/NonLinearSolver/materialLaw/mlawSMP.cpp
@@ -143,6 +143,15 @@ mlawSMP::mlawSMP(const int num,const double rho,const double alpha, const double
}
}
}
+double mlawSMP::soundSpeed() const
+{
+ double nu = _Mgl1;
+ double mu = _Mugl1;
+ double E = 2.*mu*(1.+nu);
+
+ double factornu = (1.-nu)/((1.+nu)*(1.-2.*nu));
+ return sqrt(E*factornu/_rho);
+}
mlawSMP::mlawSMP(const mlawSMP &source) : mlawThermalConducter(source),_mu_groundState3(source._mu_groundState3),_Im3(source._Im3),
_pi(source._pi), _k0(source._k0)
diff --git a/NonLinearSolver/materialLaw/mlawSMP.h b/NonLinearSolver/materialLaw/mlawSMP.h
index 615ca53e7473a6483480a73255e91bcc333b0b27..6b464b37ae1925e6920573a65e86ec2021caf459 100644
--- a/NonLinearSolver/materialLaw/mlawSMP.h
+++ b/NonLinearSolver/materialLaw/mlawSMP.h
@@ -20,8 +20,8 @@
class mlawSMP : public mlawThermalConducter
{
protected:
- const STensor43 I4;
- const STensor3 I2;
+ const STensor43 I4;
+ const STensor3 I2;
double _mu_groundState3, _Im3;
double _pi;
double _Tr,_mu_groundState2,_Nmu_groundState2, _Im2,_Sgl2,_Sr2,_Delta,_m2, _epsilonr,_n, _epsilonp02 ;
@@ -61,7 +61,7 @@ class mlawSMP : public mlawThermalConducter
public:
void setStrainOrder(const int order){_order = order;};
virtual void setMechanism2(bool mechanism2) {_mechanism2=mechanism2;};
-
+ double soundSpeed() const;
virtual double getSa0() const { return _Sa0;}
virtual double getPhia0() const { return _phia01;}
virtual double getSastar0() const { return -_b1*_phia01;}
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 061c655bd30760a1ceb94aef47d0ffc72278da55..b4da849462e421ed01d49a838f2c381d99f236ca 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -1758,6 +1758,75 @@ void SMPDG3DIPVariable::restart()
restartManager::restart(_TMIPSMP);
return;
}
+//
+PhenomenologicalSMPDG3DIPVariable::PhenomenologicalSMPDG3DIPVariable( double tinitial,const mlawPhenomenologicalSMP &_lawSMP, const bool oninter) :
+ ThermoMechanicsDG3DIPVariableBase(oninter)
+{
+ _TMIPSMP = new IPPhenomenologicalSMP(_lawSMP.getzg0());
+ this->getRefToTemperature()=tinitial;
+}
+
+PhenomenologicalSMPDG3DIPVariable::PhenomenologicalSMPDG3DIPVariable(const PhenomenologicalSMPDG3DIPVariable &source) :
+ ThermoMechanicsDG3DIPVariableBase(source)
+{
+ _TMIPSMP = NULL;
+ if (source._TMIPSMP != NULL){
+ _TMIPSMP = dynamic_cast<IPPhenomenologicalSMP*>(source._TMIPSMP->clone());
+ }
+
+}
+
+PhenomenologicalSMPDG3DIPVariable& PhenomenologicalSMPDG3DIPVariable::operator=(const IPVariable &source)
+{
+ ThermoMechanicsDG3DIPVariableBase::operator=(source);
+ const PhenomenologicalSMPDG3DIPVariable* src = dynamic_cast<const PhenomenologicalSMPDG3DIPVariable*>(&source);
+ if(src != NULL)
+ {
+ if (src->_TMIPSMP != NULL) {
+ if (_TMIPSMP!=NULL){
+ _TMIPSMP->operator=(*dynamic_cast<const IPVariable*>(src->_TMIPSMP));
+ }
+ else
+ _TMIPSMP = dynamic_cast<IPPhenomenologicalSMP*>(src->_TMIPSMP->clone());
+ }
+ else{
+ if (_TMIPSMP!=NULL){
+ delete _TMIPSMP;
+ _TMIPSMP = NULL;
+ }
+ }
+ }
+ return *this;
+}
+double PhenomenologicalSMPDG3DIPVariable::get(const int i) const
+{
+ if(i == IPField::PLASTICSTRAIN)
+ {
+ return _TMIPSMP->getRefToEquivalentPlasticDefo();
+ }
+ else
+ {
+ return ThermoMechanicsDG3DIPVariableBase::get(i);
+ }
+}
+double PhenomenologicalSMPDG3DIPVariable::defoEnergy() const
+{
+ return getIPPhenomenologicalSMP()->defoEnergy();
+}
+double PhenomenologicalSMPDG3DIPVariable::plasticEnergy() const
+{
+ return getIPPhenomenologicalSMP()->plasticEnergy();
+}
+
+void PhenomenologicalSMPDG3DIPVariable::restart()
+{
+ ThermoMechanicsDG3DIPVariableBase::restart();
+ if (_TMIPSMP!= NULL)
+ restartManager::restart(_TMIPSMP);
+ return;
+}
+
+//
ElecTherMechDG3DIPVariableBase::ElecTherMechDG3DIPVariableBase( const bool oninter) : dG3DIPVariable(oninter),
Temperature(0.), gradT(0.), fluxT(0.),dPdT(0.), dFluxTdGradT(0.), dFluxTdT(0.), dFluxTdGradV(0.), dFluxTdV(0.), TemperatureJump(0.),interfaceFluxT(0.),
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index 4c7dee491b24d3532a257a6fc087dd266d9eb32d..a09d8c832d5e6455fafacb80c46085967a316a49 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -36,6 +36,7 @@
#include "mlawLocalDamageJ2Hyper.h"
#include "mlawLinearThermoMechanics.h"
#include "mlawSMP.h"
+#include "mlawPhenomenologicalSMP.h"
#include "mlawViscoelastic.h"
#include "mlawLinearElecTherMech.h"
#include "mlawAnIsotropicElecTherMech.h"
@@ -1711,6 +1712,44 @@ class SMPDG3DIPVariable : public ThermoMechanicsDG3DIPVariableBase
virtual void restart();
};
+class PhenomenologicalSMPDG3DIPVariable : public ThermoMechanicsDG3DIPVariableBase
+{
+ protected:
+ IPPhenomenologicalSMP *_TMIPSMP;
+
+ public:
+ PhenomenologicalSMPDG3DIPVariable(const mlawPhenomenologicalSMP &_lawSMP, const bool oninter=false){Msg::Error("MPDG3DIPVariable: Cannot use the default constructor");}
+ PhenomenologicalSMPDG3DIPVariable(double tinitial,const mlawPhenomenologicalSMP &_lawSMP, const bool oninter=false);
+ PhenomenologicalSMPDG3DIPVariable(const PhenomenologicalSMPDG3DIPVariable &source);
+ PhenomenologicalSMPDG3DIPVariable& operator=(const IPVariable &source);
+
+
+ virtual void setLocation(const IPVariable::LOCATION loc) {
+ ThermoMechanicsDG3DIPVariableBase::setLocation(loc);
+ if (_TMIPSMP){
+ _TMIPSMP->setLocation(loc);
+ }
+ };
+
+ /* specific function */
+ IPPhenomenologicalSMP* getIPPhenomenologicalSMP(){return _TMIPSMP;}
+ const IPPhenomenologicalSMP* getIPPhenomenologicalSMP() const{return _TMIPSMP;}
+ virtual ~PhenomenologicalSMPDG3DIPVariable()
+ {
+ if(_TMIPSMP != NULL)
+ { delete _TMIPSMP;
+ _TMIPSMP = NULL;
+ }
+ }
+ virtual double get(const int i) const;
+ virtual double defoEnergy() const;
+ virtual double plasticEnergy() const;
+ virtual double damageEnergy() const {return _TMIPSMP->damageEnergy();};
+ virtual double getInternalEnergyExtraDofDiffusion() const {return _TMIPSMP->getThermalEnergy();};
+ virtual IPVariable* clone() const {return new PhenomenologicalSMPDG3DIPVariable(*this);};
+ virtual void restart();
+};
+
class ElecTherMechDG3DIPVariableBase : public dG3DIPVariable
{
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 96153263453451853cd947abaa2e3f5fe1407065..6d8bc5f22ace288535faf2e62984a71f2ec22296 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -16,6 +16,7 @@
#include "ipField.h"
#include "mlawLinearThermoMechanics.h"
#include "mlawSMP.h"
+#include "mlawPhenomenologicalSMP.h"
#include "mlaw.h"
#include "mlawLinearElecTherMech.h"
#include "mlawAnIsotropicElecTherMech.h"
@@ -2282,6 +2283,116 @@ void SMPDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const b
}
+PhenomenologicalSMPDG3DMaterialLaw::PhenomenologicalSMPDG3DMaterialLaw(const int num, const double rho, const double alpha, const double beta, const double gamma ,const double t0, const double Kxg,const double Kyg, const double Kzg, const double Kxr, const double Kyr, const double Kzr,const double w, const double TaylorQuiney, const double c, const double cp, const double cg, const double cr,const double Tg,const double DeltaTg, const double Eg, const double Er, const double nug, const double nur,const double alphag, const double alphar, const double zg0, const double tau_y0):
+ dG3DMaterialLaw (num,rho, true)
+{
+ _lawTMSMP = new mlawPhenomenologicalSMP(num,rho,alpha,beta,gamma,t0,Kxg,Kyg,Kzg,Kxr,Kyr,Kzr,w,TaylorQuiney,c,cp,cg,cr,Tg,DeltaTg,Eg,Er,nug,nur,alphag,alphar,zg0,tau_y0);
+
+ double nu = 0.3; //to be updated
+ double mu = 1.e6; //to be updated
+ double E = 2.*mu*(1.+nu); //Where does this come from? To Have the E corresponding to the Vmax, MUmax
+ //In the linear & homogeneous & isotropic case ?
+
+ linearK = _lawTMSMP->getInitialConductivityTensor();
+ Stiff_alphaDilatation=_lawTMSMP->getStiff_alphaDilatation();
+
+ fillElasticStiffness(E, nu, elasticStiffness);
+}
+
+PhenomenologicalSMPDG3DMaterialLaw::PhenomenologicalSMPDG3DMaterialLaw(const PhenomenologicalSMPDG3DMaterialLaw &source) :
+ dG3DMaterialLaw (source)// , dG3DMaterialLaw (source)
+{
+ _lawTMSMP = NULL;
+ if (source._lawTMSMP != NULL){
+ _lawTMSMP = new mlawPhenomenologicalSMP(*(source._lawTMSMP));
+ }
+ linearK = source.linearK;
+ Stiff_alphaDilatation = source.Stiff_alphaDilatation;
+}
+
+void PhenomenologicalSMPDG3DMaterialLaw::setStrainOrder(const int order ){
+ _lawTMSMP->setStrainOrder(order);
+};
+
+double PhenomenologicalSMPDG3DMaterialLaw::getExtraDofStoredEnergyPerUnitField(double zg) const
+{
+ return _lawTMSMP->getExtraDofStoredEnergyPerUnitField(zg);
+}
+
+void PhenomenologicalSMPDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+ // check interface element
+ bool inter=true;
+ const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+ if(iele==NULL) inter=false;
+ IPVariable* ipvi = new PhenomenologicalSMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
+ IPVariable* ipv1 = new PhenomenologicalSMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
+ IPVariable* ipv2 = new PhenomenologicalSMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
+ if(ips != NULL) delete ips;
+ ips = new IP3State(state_,ipvi,ipv1,ipv2);
+}
+
+void PhenomenologicalSMPDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+ if(ipv !=NULL) delete ipv;
+ bool inter=true;
+ const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+ if(iele == NULL){
+ inter=false;
+ }
+ ipv = new PhenomenologicalSMPDG3DIPVariable(_lawTMSMP->getInitialTemperature(),*_lawTMSMP,inter);
+}
+
+
+void PhenomenologicalSMPDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
+{
+ /* get ipvariable */
+ PhenomenologicalSMPDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
+ const PhenomenologicalSMPDG3DIPVariable* ipvprev; //= static_cast <const nonLocalDamageDG3DIPVariable *> (ipvp);
+
+ FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
+ if(ipvtmp !=NULL)
+ {
+
+ ipvcur = static_cast<PhenomenologicalSMPDG3DIPVariable*>(ipvtmp->getIPvBulk());
+ const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
+ ipvprev = static_cast<const PhenomenologicalSMPDG3DIPVariable*>(ipvtmp2->getIPvBulk());
+ }
+ else
+ {
+ ipvcur = static_cast<PhenomenologicalSMPDG3DIPVariable*>(ipv);
+ ipvprev = static_cast<const PhenomenologicalSMPDG3DIPVariable*>(ipvp);
+ }
+
+ /* strain xyz */
+ const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
+ const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
+ const double temperature = ipvcur->getConstRefToTemperature();
+ const double T0 = ipvprev->getConstRefToTemperature();
+ const SVector3& gradT= ipvcur->getConstRefToGradT();
+
+ STensor3& dPdT=ipvcur->getRefTodPdT();
+ SVector3& fluxT=ipvcur->getRefToThermalFlux();
+ STensor3& dqdgradT=ipvcur->getRefTodThermalFluxdGradT();
+ SVector3& dqdT=ipvcur->getRefTodThermalFluxdT();
+ STensor33& dqdF=ipvcur->getRefTodThermalFluxdF();
+
+ STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
+ STensor43& Tangent=ipvcur->getRefToTangentModuli();
+ double &w = ipvcur->getRefToThermalSource();
+ double &dwdt = ipvcur->getRefTodThermalSourcedField();
+ STensor3& dwdf = ipvcur->getRefTodThermalSourcedF();
+ /* data for J2 law */
+ IPPhenomenologicalSMP* q1 = ipvcur->getIPPhenomenologicalSMP();
+ const IPPhenomenologicalSMP* q0 = ipvprev->getIPPhenomenologicalSMP();
+ _lawTMSMP->constitutive(F0,F1, P,q0, q1,Tangent,T0,temperature, gradT,fluxT,dPdT,dqdgradT,dqdT,dqdF,stiff,w,dwdt,dwdf);
+
+
+ ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
+ ipvcur->setRefToLinearConductivity(this->linearK);
+ ipvcur->setRefToStiff_alphaDilatation(this->Stiff_alphaDilatation);
+}
+//
LinearElecTherMechDG3DMaterialLaw::LinearElecTherMechDG3DMaterialLaw(const int num,const double rho,const double Ex, const double Ey, const double Ez, const double Vxy,
const double Vxz, const double Vyz,const double MUxy, const double MUxz, const double MUyz, const double alpha, const double beta, const double gamma,
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 0764e98adce6d7d8589cc24be4f71f5424efbc52..fd8fe2d9bcd4f0b682ab6d943a566a74c27c83f6 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -23,6 +23,7 @@
#include "mlawLinearThermoMechanics.h"
#include "mlawJ2VMSmallStrain.h"
#include "mlawSMP.h"
+#include "mlawPhenomenologicalSMP.h"
#include "mlawEOS.h"
#include "mlawViscoelastic.h"
#include "mlawLinearElecTherMech.h"
@@ -1022,6 +1023,49 @@ class SMPDG3DMaterialLaw :public dG3DMaterialLaw // ,public mlawSMP
};
+class PhenomenologicalSMPDG3DMaterialLaw :public dG3DMaterialLaw // ,public mlawSMP
+{
+
+ protected:
+ mlawPhenomenologicalSMP *_lawTMSMP; // pointer to allow to choose between LLN style or Gmsh Style. The choice is performed by the constructor (2 constructors with != arguments)
+
+ STensor3 linearK;
+ STensor3 Stiff_alphaDilatation;
+ public:
+
+ PhenomenologicalSMPDG3DMaterialLaw(const int num, const double rho, const double alpha, const double beta, const double gamma ,const double t0, const double Kxg,const double Kyg, const double Kzg, const double Kxr, const double Kyr, const double Kzr,const double w, const double TaylorQuiney, const double c, const double cp, const double cg, const double cr,const double Tg,const double DeltaTg, const double Eg, const double Er, const double nug, const double nur,const double alphag, const double alphar, const double zg0, const double tau_y0);
+
+ void setStrainOrder(const int order);
+
+#ifndef SWIG
+ virtual ~PhenomenologicalSMPDG3DMaterialLaw() {
+ if (_lawTMSMP !=NULL) {delete _lawTMSMP; _lawTMSMP = NULL;};
+ };
+ PhenomenologicalSMPDG3DMaterialLaw(const PhenomenologicalSMPDG3DMaterialLaw &source);
+ // set the time of _nldlaw
+ virtual void setTime(const double t,const double dtime){
+ dG3DMaterialLaw::setTime(t,dtime);
+ _lawTMSMP->setTime(t,dtime);
+ }
+ virtual materialLaw::matname getType() const {return _lawTMSMP->getType();}
+ virtual void createIPState(IPStateBase* &ips,const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const;
+ virtual void createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const;
+ virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}
+ virtual double soundSpeed() const{return _lawTMSMP->soundSpeed();} // or change value ??
+ virtual void stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff=true, const bool checkfrac=true);
+ virtual int getNumberOfExtraDofsDiffusion() const {return 1;}
+ double getExtraDofStoredEnergyPerUnitField(double T) const;
+ virtual materialLaw* clone() const{ return new PhenomenologicalSMPDG3DMaterialLaw(*this);};
+ virtual bool withEnergyDissipation() const {return _lawTMSMP->withEnergyDissipation();};
+ virtual void setMacroSolver(const nonLinearMechSolver* sv){
+ dG3DMaterialLaw::setMacroSolver(sv);
+ if (_lawTMSMP!=NULL)
+ _lawTMSMP->setMacroSolver(sv);
+ }
+#endif
+};
+
+
class LinearElecTherMechDG3DMaterialLaw : public dG3DMaterialLaw // public materialLaw
{
diff --git a/dG3D/src/dG3DTerms.cpp b/dG3D/src/dG3DTerms.cpp
index ce7d9d855d9071dda73bd8d22ebc80701faa673c..4b3d48ccaa417b3b9a17dcf7602c07b31c418e42 100644
--- a/dG3D/src/dG3DTerms.cpp
+++ b/dG3D/src/dG3DTerms.cpp
@@ -1087,6 +1087,7 @@ void dG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<double>
}
}
}
+ //mStiff.print("analytical bulk");
}
}