diff --git a/NonLinearSolver/internalPoints/CMakeLists.txt b/NonLinearSolver/internalPoints/CMakeLists.txt
index 3deba2244d968abcc638694c913fd0b276023aca..242c455c6aeca1baf5c039b5cc7c8999bc8936d4 100644
--- a/NonLinearSolver/internalPoints/CMakeLists.txt
+++ b/NonLinearSolver/internalPoints/CMakeLists.txt
@@ -31,8 +31,9 @@ set(SRC
ipJ2SmallStrains.cpp
ipFiniteStrain.cpp
ipSMP.cpp
+ ipPhenomenologicalSMP.cpp
ipViscoelastic.cpp
- ipvariable.cpp
+ ipvariable.cpp
# inline
ipvariable.h
ipEOS.cpp
@@ -40,7 +41,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..f99a75b246a35b75652ac605da7b477fb06ec5c7
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.cpp
@@ -0,0 +1,59 @@
+//
+// 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():IPVariableMechanics(),_thermalEnergy(0.) , _eps(0.),_SMPEnergy(0.), xi(0.)
+{
+ STensorOperation::unity(Fp);
+
+};
+
+
+IPPhenomenologicalSMP::IPPhenomenologicalSMP(const IPPhenomenologicalSMP &source) :IPVariableMechanics()//: IPLinearThermoMechanics(source) //n
+{
+ _SMPEnergy=source._SMPEnergy;//n
+ Fp = source.Fp;
+ xi = source.xi;
+ _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)
+ {
+ Fp = src->Fp;
+ xi = src->xi;
+ _thermalEnergy = src->_thermalEnergy;
+ _eps=src->_eps;
+ _SMPEnergy=src->_SMPEnergy;
+ }
+ return *this;
+}
+
+
+void IPPhenomenologicalSMP::restart()
+{
+ IPVariableMechanics::restart();
+ restartManager::restart(Fp);
+ restartManager::restart(xi);
+ 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..a79bf434672ee1b9f7acd98291c4b69d231bbed5
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipPhenomenologicalSMP.h
@@ -0,0 +1,52 @@
+//
+// 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 Fp;
+ double xi;
+ //double _SMPEnergy;
+ double _thermalEnergy;
+ double _eps; // equivalent plastic strain
+ public:
+ IPPhenomenologicalSMP();
+
+ //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 &getRefToFp() {return Fp;}
+ virtual const STensor3 &getConstRefToFp() const {return Fp;}
+ virtual double &getRefToXi() {return xi;}
+ virtual double getXi() const{return xi;}
+ 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 1de59a4b96761b52aba2a5c361d3d2e3d79ab2e4..0ee59de6b6ceb20e5db3e7057f11b71933576cc6 100644
--- a/NonLinearSolver/materialLaw/CMakeLists.txt
+++ b/NonLinearSolver/materialLaw/CMakeLists.txt
@@ -43,6 +43,7 @@ set(SRC
mlawLinearElecTherMech.cpp
mlawAnIsotropicElecTherMech.cpp
mlawSMP.cpp
+ mlawPhenomenologicalSMP.cpp
mlawEOS.cpp
mlawViscoelastic.cpp
mlawHyperelastic.cpp
@@ -56,9 +57,9 @@ set(SRC
mlawLocalDamageJ2Hyper.cpp
mlawAnIsotropicTherMech.cpp
mlawJ2VMSmallStrain.cpp
- cohesiveDamageLaw.cpp
+ cohesiveDamageLaw.cpp
mlawHyperelasticDamage.cpp
- elasticPotential.cpp
+ elasticPotential.cpp
cohesiveElasticPotential.cpp
mlawAnisotropicStoch.cpp
# src/op_eshelby.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..ef59746dc45523ef8f97e37dfbb5e65918c45a75
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.cpp
@@ -0,0 +1,334 @@
+//
+// 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 "mlawThermalConducter.h"
+#include "MInterfaceElement.h"
+#include "ipPhenomenologicalSMP.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 Kx,const double Ky,const double Kz)
+ :mlawThermalConducter(num,alpha, beta, gamma, Kx, Ky,Kz) ,_t0(t0),_order(-1)
+
+
+{
+ _k0=_k;
+ 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;
+ double _alphagl1=1.; //to be updated
+ alphaDilatation(0,0)= _alphagl1;
+ alphaDilatation(1,1)= _alphagl1;
+ alphaDilatation(2,2)= _alphagl1;
+ 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.3; // TO BE UPDATED
+ double mu = 1e6; // TO BE UPDATED
+ double E = 2.*mu*(1.+nu);
+ 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);
+ }
+ }
+ }
+ }
+}
+
+mlawPhenomenologicalSMP::mlawPhenomenologicalSMP(const mlawPhenomenologicalSMP &source) : mlawThermalConducter(source),_t0(source._t0), _k0(source._k0),
+ _Stiff_alphaDilatation(source._Stiff_alphaDilatation),_alphaDilatation(source._alphaDilatation), _order(source._order)
+
+{
+
+}
+
+
+//mlawSMP& mlawSMP::operator=(const mlawLinearThermoMechanics &source)
+mlawPhenomenologicalSMP& mlawPhenomenologicalSMP::operator=(const materialLaw &source)
+{
+ mlawThermalConducter::operator=(source);
+ const mlawPhenomenologicalSMP* src =static_cast<const mlawPhenomenologicalSMP*>(&source);
+ _k0=src->_k0;
+ _t0=src->_t0;
+ _Stiff_alphaDilatation=src->_Stiff_alphaDilatation;
+ _alphaDilatation=src->_alphaDilatation;
+ _order =src->_order;
+}
+
+double mlawPhenomenologicalSMP::getExtraDofStoredEnergyPerUnitField(double T) const
+{
+ //here is an approximation for constant
+// double Tg=_Tr;
+ double specificheat=0.; //_c0;
+ //if(T<Tg)
+ //{
+ // specificheat=_c0-_c1*(T-Tg);
+ //}
+ 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();
+ IPVariable* ipv1 = new IPPhenomenologicalSMP();
+ IPVariable* ipv2 = new IPPhenomenologicalSMP();
+ 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 &w,
+ double &dwdT,
+ STensor3 &dwdF
+ ) const
+{
+ w=0.;dwdT=0.; STensorOperation::zero(dwdF); STensorOperation::zero(dqdF); STensorOperation::zero(dqdT); STensorOperation::zero(dqdgradT); STensorOperation::zero(dPdT); STensorOperation::zero(fluxT);
+ // _k0=_k;
+ double Tg=0.;;
+ // const double pi=3.14159;
+ STensorOperation::zero(P);
+ STensorOperation::zero(Tangent);
+
+ static STensor3 Cedis3;
+ STensorOperation::zero(Cedis3);
+
+ double specificheat=0.;
+
+ double dh;
+ double dcpdt=0.;
+ dh = getTimeStep();
+#if 0
+ //cp is not constant
+ if(T<Tg)
+ {
+ specificheat=_c0-_c1*(T-Tg);
+ }
+ else
+ specificheat=_c0;
+
+ if(dh>0.)
+ {
+ w=-specificheat*(T-T0)/dh;// the - sign become from the term of thermal source
+ w+=(depsitau1)*_wp/dh;
+ w+=(depsitau2)*_wp/dh;
+ }
+ else
+ w=0.;
+ // q1->_thermalEnergy = specificheat*T; //(T-_t0);
+
+ if (T<Tg)
+ {
+ dcpdt=-_c1;
+ dcpdF=dTgdF;
+ dcpdF*=_c1;
+ }
+ else
+ {
+ dcpdt=0.;
+ STensorOperation::zero(dcpdF);
+ }
+
+ if(dh>0.)
+ {
+ dwdT=-dcpdt*(T-T0)/dh-specificheat/dh;
+ dwdT+=(dtauepsilon1dt)*_wp/dh;
+ dwdT+=(dtauepsilon2dt)*_wp/dh;
+
+
+ for(int K = 0; K< 3; K++){
+ for(int i = 0; i< 3; i++){
+ dwdF(K,i)=-dcpdF(K,i)*(T-T0)/dh;
+ dwdF(K,i)+=(dtauepsilon1dF(K,i))*_wp/dh;
+ dwdF(K,i)+=(dtauepsilon2dF(K,i))*_wp/dh;
+ }
+ }
+ }
+ else
+ {
+ dwdT=0.;
+ STensorOperation::zero(dwdF);
+ }
+#endif
+// P =P3;
+ //P.print("p");
+ //put flux in reference configuration
+ static STensor3 Fninv;
+ STensorOperation::inverseSTensor3(Fn, Fninv);
+ double J= Fn.determinant();
+ static STensor3 _Kref;
+ STensorOperation::zero(_Kref);
+ 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);//*J; //put in reference configuration
+ }
+ }
+ }
+ }
+ STensorOperation::multSTensor3SVector3(_Kref, gradT, fluxT);
+ fluxT*=J;
+
+ static STensor3 dkdT;
+ STensorOperation::zero(dkdT);
+ if(stiff)
+ {
+
+ //dPdT=dPdT3;
+ //Tangent=Tangent3;
+ 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=");
+}
+
+
diff --git a/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.h b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.h
new file mode 100644
index 0000000000000000000000000000000000000000..b056e1ed0e5985b88f23116bb7f2bdd4380189a1
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mlawPhenomenologicalSMP.h
@@ -0,0 +1,89 @@
+//
+// 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 "mlawThermalConducter.h"
+#include "../internalPoints/ipPhenomenologicalSMP.h"
+#include "mlaw.h"
+#include "STensor3.h"
+#include "STensor43.h"
+#include "STensor63.h"
+#include "ipSMP.h"
+
+
+class mlawPhenomenologicalSMP : public mlawThermalConducter
+{
+ protected:
+ double _t0;
+
+ STensor3 _k0;
+ STensor3 _alphaDilatation;
+ STensor3 _Stiff_alphaDilatation;
+
+ int _order;
+
+ public:
+ mlawPhenomenologicalSMP(const int num, const double rho, const double alpha, const double beta, const double gamma ,const double t0, const double Kx,const double Ky, const double Kz);
+#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 T) const;
+ virtual double getInitialExtraDofStoredEnergyPerUnitField() const {getExtraDofStoredEnergyPerUnitField(_t0);}
+ 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 &w,
+ double &dwdt,
+ STensor3 &dwdf
+ ) const;
+
+ const STensor3& getInitialConductivityTensor() const { return _k0;};
+ virtual const STensor3& getStiff_alphaDilatation()const { return _Stiff_alphaDilatation;};
+
+
+ protected:
+
+ #endif // SWIG
+};
+
+#endif // MLAWPHENOMENOLOGICALSMP_H_
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 8dfb384f896385ce00c28bfe05e43df565bac798..c2e533c45b62bf6b8ea2adc78e37211f55c568b7 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -2079,6 +2079,75 @@ void SMPDG3DIPVariable::restart()
restartManager::restart(_TMIPSMP);
return;
}
+//
+PhenomenologicalSMPDG3DIPVariable::PhenomenologicalSMPDG3DIPVariable( double tinitial,const mlawPhenomenologicalSMP &_lawSMP, const bool oninter) :
+ ThermoMechanicsDG3DIPVariableBase(oninter)
+{
+ _TMIPSMP = new IPPhenomenologicalSMP();
+ 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 55a76c13b47eaff9a584823e70c374b85742caf3..477fae6a607f9756b079c109c47e4c3fcf077e30 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -42,6 +42,7 @@
#include "mlawLinearThermoMechanics.h"
#include "mlawNonLocalDamageGurson.h"
#include "mlawSMP.h"
+#include "mlawPhenomenologicalSMP.h"
#include "mlawViscoelastic.h"
#include "mlawLinearElecTherMech.h"
#include "mlawAnIsotropicElecTherMech.h"
@@ -1511,6 +1512,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 fce62a2ae640b03b2a8fe90b3797fb18ae2435b8..b12c4faf5966c9d1caa5b19612cb7d8700218367 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -2905,7 +2905,118 @@ void SMPDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const b
ipvcur->setRefToStiff_alphaDilatation(this->Stiff_alphaDilatation);
}
+//
+PhenomenologicalSMPDG3DMaterialLaw::PhenomenologicalSMPDG3DMaterialLaw(const int num,const double rho,const double alpha, const double beta, const double gamma ,
+ const double t0, const double Kx,const double Ky, const double Kz):
+ dG3DMaterialLaw (num,rho, true)
+{
+ _lawTMSMP = new mlawPhenomenologicalSMP( num,rho,alpha,beta,gamma,t0,Kx,Ky,Kz);
+
+ 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 T) const
+{
+ return _lawTMSMP->getExtraDofStoredEnergyPerUnitField(T);
+}
+
+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->setRefToElasticTangentModuli(this->elasticStiffness);
+ ipvcur->setRefTolinearK(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 663f96b8c39a2cd7c5580fbd04fc5f19f78155b4..65798665de01f3395ea71553513cf9676346de42 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -24,6 +24,7 @@
#include "mlawLinearThermoMechanics.h"
#include "mlawJ2VMSmallStrain.h"
#include "mlawSMP.h"
+#include "mlawPhenomenologicalSMP.h"
#include "mlawEOS.h"
#include "mlawViscoelastic.h"
#include "mlawLinearElecTherMech.h"
@@ -1254,6 +1255,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 Kx,const double Ky, const double Kz);
+
+ 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
{