Merge branch 'Coal_jl2' into 'master'

Coal jl2

See merge request !50

NonLinearSolver/internalPoints/CMakeLists.txt

@@ -19,7 +19,6 @@ set(SRC
   ipNucleation.cpp
   ipCoalescence.cpp
   ipNonLocalPorosity.cpp
-  ipNonLocalDamageGursonThermoMechanics.cpp
   ipTransverseIsotropic.cpp
   ipTransverseIsoCurvature.cpp
   ipTransverseIsoYarnB.cpp

NonLinearSolver/internalPoints/ipNonLocalDamageGursonThermoMechanics.cpp

-//
-// Description: storing class for Non Local Damage Gurson Thermo Mechanics
-//
-//
-// Author:  <L. Noels>, (C) 2014
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-#include "ipNonLocalDamageGursonThermoMechanics.h"
-
-
-IPNonLocalDamageGursonThermoMechanics::IPNonLocalDamageGursonThermoMechanics() : IPNonLocalPorosity(),_thermalEnergy(0.), _DirreversibleEnergyDT(0.), _ThermoElasticCoupling(0.)
-{};
-
-IPNonLocalDamageGursonThermoMechanics::IPNonLocalDamageGursonThermoMechanics(const IPNonLocalDamageGursonThermoMechanics &source) : IPNonLocalPorosity(source), _thermalEnergy(source._thermalEnergy), _DirreversibleEnergyDT(0.), _ThermoElasticCoupling(0.) {};
-
-IPNonLocalDamageGursonThermoMechanics::IPNonLocalDamageGursonThermoMechanics(double fVinitial, const J2IsotropicHardening *j2IH, const CLengthLaw *cll,
-                         const std::vector<NucleationLaw*> *ipvgdnContainer,const CoalescenceLaw *coaleslaw) :IPNonLocalPorosity(fVinitial,j2IH,cll,ipvgdnContainer,coaleslaw),  _thermalEnergy(0.), 
-                                                                                                              _DirreversibleEnergyDT(0.), _ThermoElasticCoupling(0.) {};
- 
-IPNonLocalDamageGursonThermoMechanics& IPNonLocalDamageGursonThermoMechanics::operator=(const IPVariable &source)
-{
-  IPNonLocalPorosity::operator=(source);
-  const IPNonLocalDamageGursonThermoMechanics* src = dynamic_cast<const IPNonLocalDamageGursonThermoMechanics*>(&source);
-  if(src)
-  {
-    _DirreversibleEnergyDT = src->_DirreversibleEnergyDT;
-    _ThermoElasticCoupling = src->_ThermoElasticCoupling;
-    _thermalEnergy=src->_thermalEnergy;
-  }
-  return *this;
-}
-
-double IPNonLocalDamageGursonThermoMechanics::defoEnergy() const
-{
-  return IPNonLocalPorosity::defoEnergy();
-
-}
-double IPNonLocalDamageGursonThermoMechanics::plasticEnergy() const
-{
-  return IPNonLocalPorosity::plasticEnergy();
-}
-
-
-void IPNonLocalDamageGursonThermoMechanics::restart(){
-  IPNonLocalPorosity::restart();
-  restartManager::restart(_thermalEnergy);
-  restartManager::restart(_DirreversibleEnergyDT);
-  restartManager::restart(_ThermoElasticCoupling);
-}

NonLinearSolver/internalPoints/ipNonLocalDamageGursonThermoMechanics.h

-//
-// Description: storing class for non local damage gurson thermomechanical law
-// Author:  <L. Noels>, (C) 2014
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-#ifndef IPNONLOCALDAMAGEGURSONTHERMOMECHANICS_H_
-#define IPNONLOCALDAMAGEGURSONTHERMOMECHANICS_H_
-#include "ipNonLocalPorosity.h"
-#include "STensor3.h"
-
-class IPNonLocalDamageGursonThermoMechanics : public IPNonLocalPorosity
-{
- public:
-  double _thermalEnergy; //
-  double _DirreversibleEnergyDT;
-  double _ThermoElasticCoupling;
-  
- 
-
- public:
-  IPNonLocalDamageGursonThermoMechanics();
-  IPNonLocalDamageGursonThermoMechanics(const IPNonLocalDamageGursonThermoMechanics &source);
-  IPNonLocalDamageGursonThermoMechanics& operator=(const IPVariable &source);
-  IPNonLocalDamageGursonThermoMechanics(double fVinitial, const J2IsotropicHardening *j2IH, const CLengthLaw *cll,
-                         const std::vector<NucleationLaw*> *ipvgdnContainer,const CoalescenceLaw *coaleslaw);
-  virtual const double& getConstRefToDIrreversibleEnergyDT() const{return _DirreversibleEnergyDT;};
-  virtual double& getRefToDIrreversibleEnergyDT() {return _DirreversibleEnergyDT;};
-
-  virtual const double& getConstRefToThermoelasticCoupling() const{return _ThermoElasticCoupling;};
-  virtual double& getRefToThermoelasticCoupling() {return _ThermoElasticCoupling;};
-  
-  virtual ~IPNonLocalDamageGursonThermoMechanics(){}
-  virtual double defoEnergy() const;
-  virtual double plasticEnergy() const;
-  virtual double getThermalEnergy() const {return _thermalEnergy;};
-  IPNonLocalDamageGursonThermoMechanics& operator=(const IPNonLocalPorosity &source);
-  virtual IPVariable* clone() const {return new IPNonLocalDamageGursonThermoMechanics(*this);};
-  virtual void restart();
-};
-
-#endif // IPNONLOCALDAMAGEGURSONTHERMOMECHANICS_H_

NonLinearSolver/internalPoints/ipNonLocalPorosity.cpp

@@ -15,7 +15,7 @@ IPNonLocalPorosity::IPNonLocalPorosity() : IPVariableMechanics(),  _elasticEnerg
                                                    _eplmatrix(0.), _fV(0),  _fVstar(0), _fVtilde(0.), _yieldfV(0.),
                                                    _dissipationBlocked(false),_NonLocalToLocal(false), _Fp(1.),
                                                    _Ee(0.),_DirreversibleEnergyDF(0.), _DirreversibleEnergyDtildefV(0.),
-						  _currentOutwardNormal(NULL),_referenceOutwardNormal(NULL),_corKir(0.), _yieldStress(0.),
+																									_currentOutwardNormal(NULL),_referenceOutwardNormal(NULL),_corKir(0.), _yieldStress(0.),
                                                   _failed(false), _hatq(0.),_nonlocalEplmatrix(0.),_nonlocalHatQ(0.), _fVinitial(0.)
 {
   ipvJ2IsotropicHardening=NULL;
@@ -34,8 +34,8 @@ IPNonLocalPorosity::IPNonLocalPorosity(double fVinitial, const J2IsotropicHarden
                                                  _fVstar(fVinitial),_fVtilde(fVinitial), _yieldfV(fVinitial),
                                                   _Fp(1.),_Ee(0.), _DirreversibleEnergyDF(0.), _DirreversibleEnergyDtildefV(0.),
                                                  _dissipationBlocked(false), _NonLocalToLocal(false),
-                                                 _currentOutwardNormal(NULL), _referenceOutwardNormal(NULL),_corKir(0.),
-                                                 _failed(false),_hatq(0.),_nonlocalEplmatrix(0.),_nonlocalHatQ(0.), _fVinitial(fVinitial)
+																								_currentOutwardNormal(NULL), _referenceOutwardNormal(NULL),_corKir(0.),
+                                                _failed(false),_hatq(0.),_nonlocalEplmatrix(0.),_nonlocalHatQ(0.), _fVinitial(fVinitial)
 {
   ipvJ2IsotropicHardening=NULL;
   if(j2IH ==NULL) Msg::Error("IPNonLocalPorosity::IPNonLocalPorosity has no j2IH");
@@ -59,7 +59,7 @@ IPNonLocalPorosity::IPNonLocalPorosity(double fVinitial, const J2IsotropicHarden
   ipvCoal=NULL;
   if(coaleslaw ==NULL) Msg::Error("IPNonLocalPorosity::IPNonLocalPorosity has no Coalescence law");
   coaleslaw->createIPVariable(ipvCoal);
-  
+  coaleslaw->computeInitialChi(ipvCoal,fVinitial);
 
   if(fVinitial <0. or fVinitial >1.) Msg::Error("IPNonLocalPorosity::IPNonLocalPorosity wrong initial porosity");
 };
@@ -226,14 +226,14 @@ void IPNonLocalPorosity::restart()
   restartManager::restart(_DirreversibleEnergyDF);
   restartManager::restart(_DirreversibleEnergyDtildefV);
   restartManager::restart(_corKir);
-  restartManager::restart(_yieldStress);
+	restartManager::restart(_yieldStress);
   restartManager::restart(_failed);
   restartManager::restart(_hatq);
   restartManager::restart(_nonlocalEplmatrix);
   restartManager::restart(_nonlocalHatQ);
   restartManager::restart(_fVinitial);
-  _currentOutwardNormal = NULL;
-  _referenceOutwardNormal = NULL;
+	_currentOutwardNormal = NULL;
+	_referenceOutwardNormal = NULL;
   return;
   
 }

NonLinearSolver/internalPoints/ipNonLocalPorosity.h

@@ -31,8 +31,8 @@ class IPNonLocalPorosity : public IPVariableMechanics
   std::vector<IPNucleation*> ipvgdnContainer; // Nucleation law vector
   IPCoalescence* ipvCoal; // Coalescence law
 
-  const SVector3* _currentOutwardNormal; // current normal at Gauss point  
-  const SVector3* _referenceOutwardNormal; // refrence nomral at Gauss point
+	const SVector3* _currentOutwardNormal; // current normal at Gauss point 
+	const SVector3* _referenceOutwardNormal; // refrence nomral at Gauss point
   
   // Internal variables
   double _fVinitial; // initial porosity
@@ -46,11 +46,11 @@ class IPNonLocalPorosity : public IPVariableMechanics
   double _hatq;  // volumetric plastic deformation = tr(Dp)
   double _nonlocalHatQ; // nonlocal volumetric plastic deformation 
 
-  double _yieldStress; // current viscoplastic yield stress = initial part + hardening part+ viscoplastic part
+	double _yieldStress; // current viscoplastic yield stress = initial part + hardening part+ viscoplastic part
 
   STensor3 _Fp;    // plastic part of the deformation gradient
   STensor3 _Ee; // elastic natural strain tensor
-  STensor3 _corKir; // corotational Kirchhoff stress
+	STensor3 _corKir; // corotational Kirchhoff stress
 
   // Damage and transition managing
   bool _dissipationBlocked; // True if dissipation is blocked at the IPv
@@ -63,6 +63,8 @@ class IPNonLocalPorosity : public IPVariableMechanics
   double _irreversibleEnergy;
   STensor3 _DirreversibleEnergyDF;
   double _DirreversibleEnergyDtildefV;
+
+
  public:
   // Constructor & destructor
   IPNonLocalPorosity();
@@ -137,7 +139,7 @@ class IPNonLocalPorosity : public IPVariableMechanics
   virtual double damageEnergy() const {return 0.;}; // dissipation by damage
   
   virtual double irreversibleEnergy() const {return _irreversibleEnergy;};
-  virtual double & getRefToIrreversibleEnergy() {return _irreversibleEnergy;};
+	virtual double & getRefToIrreversibleEnergy() {return _irreversibleEnergy;};
   
   virtual const STensor3& getConstRefToDIrreversibleEnergyDF() const{return _DirreversibleEnergyDF;};
   virtual STensor3& getRefToDIrreversibleEnergyDF() {return _DirreversibleEnergyDF;};
@@ -145,12 +147,6 @@ class IPNonLocalPorosity : public IPVariableMechanics
   virtual const double& getConstRefToDIrreversibleEnergyDtildefV() const{return _DirreversibleEnergyDtildefV;};
   virtual double& getRefToDIrreversibleEnergyDtildefV() {return _DirreversibleEnergyDtildefV;};
 
-  virtual const double& getConstRefToDIrreversibleEnergyDT() const{ Msg::Fatal("getRefToDIrreversibleEnergyDT only defined in thermoelasticity");};//-------added
-  virtual double& getRefToDIrreversibleEnergyDT() {Msg::Fatal("getRefToDIrreversibleEnergyDT only defined in thermoelasticity");};//-------added
-
-  virtual const double& getConstRefToThermoelasticCoupling() const{Msg::Fatal("getRefToDIrreversibleEnergyDT only defined in thermoelasticity");};//------added
-  virtual double& getRefToThermoelasticCoupling() {Msg::Fatal("getRefToDIrreversibleEnergyDT only defined in thermoelasticity");};//------added
-
   // Access functions - Internal variables
   virtual double getDamage() const{return getNonLocalPorosity();};
   
@@ -181,8 +177,8 @@ class IPNonLocalPorosity : public IPVariableMechanics
   virtual double getYieldPorosity() const {return _yieldfV;};
   virtual double& getRefToYieldPorosity() {return _yieldfV;};
   
-  virtual double getCurrentViscoplasticYieldStress() const {return _yieldStress;};
-  virtual double& getRefToCurrentViscoplasticYieldStress() {return _yieldStress;};
+	virtual double getCurrentViscoplasticYieldStress() const {return _yieldStress;};
+	virtual double& getRefToCurrentViscoplasticYieldStress() {return _yieldStress;};
       
   virtual const STensor3& getConstRefToFp() const {return _Fp;}
   virtual STensor3& getRefToFp(){return _Fp;}
@@ -200,7 +196,7 @@ class IPNonLocalPorosity : public IPVariableMechanics
   virtual bool getCoalescenceActiveFlag() const {return ipvCoal->getCoalescenceActiveFlag();};
   virtual void setCoalescenceActiveFlag(const bool fl){ipvCoal->getRefToCoalescenceActiveFlag() = fl;};
  
-  virtual const STensor3& getConstRefToCorotationalKirchhoffStress() const {return _corKir;};
+	virtual const STensor3& getConstRefToCorotationalKirchhoffStress() const {return _corKir;};
   virtual STensor3& getRefToCorotationalKirchhoffStress() {return _corKir;};
 
   // Access functions - Included IPVs

NonLinearSolver/internalPoints/ipPhenomenologicalSMP.cpp

@@ -12,12 +12,12 @@
 #include "restartManager.h"
 #include "STensorOperations.h"
 
-IPPhenomenologicalSMP::IPPhenomenologicalSMP():IPVariableMechanics(),_thermalEnergy(0.) , _eps(0.),_SMPEnergy(0.)
+IPPhenomenologicalSMP::IPPhenomenologicalSMP(double _zg):IPVariableMechanics(),_thermalEnergy(0.) , _eps(0.),_SMPEnergy(0.)
 {
   STensorOperation::unity(Fpg);
   STensorOperation::unity(Fp);
   STensorOperation::unity(Ff);
-  zg=0.;
+  zg=_zg;
 };
 
 

NonLinearSolver/internalPoints/ipPhenomenologicalSMP.h

@@ -12,25 +12,25 @@
 #include "ipLinearThermoMechanics.h"
 #include "STensor3.h"
 #include "ipvariable.h"
-#include "ipHardening.h"
-#include "j2IsotropicHardening.h"
+//#include "mlawSMP.h"
 
 class IPPhenomenologicalSMP : public IPVariableMechanics//: public IPLinearThermoMechanics  //, public mlawSMP
 {
 public:
-  //double _SMPEnergy;
+   double _SMPEnergy;
 protected:
-  IPJ2IsotropicHardening* ipvJ2IsotropicHardening;
   STensor3 Fpg,Fp,Ff;
   double zg;
-  double _SMPEnergy;
+  //double _SMPEnergy;
   double _thermalEnergy;
   double _eps; // equivalent plastic strain
  public:
-  IPPhenomenologicalSMP();  
-  IPPhenomenologicalSMP(const J2IsotropicHardening *j2IH);
+  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;};
 
@@ -47,18 +47,6 @@ protected:
 
   virtual double& getRefToEquivalentPlasticDefo() {return _eps;};
   virtual const double& getConstRefToEquivalentPlasticDefo() const {return _eps;};
-  virtual const IPJ2IsotropicHardening &getConstRefToIPJ2IsotropicHardening() const
-  {
-    if(ipvJ2IsotropicHardening==NULL)
-      Msg::Error("IPJ2linear: ipvJ2IsotropicHardening not initialized");
-    return *ipvJ2IsotropicHardening;
-  }
-  virtual IPJ2IsotropicHardening &getRefToIPJ2IsotropicHardening()
-  {
-    if(ipvJ2IsotropicHardening==NULL)
-      Msg::Error("IPJ2linear: ipvJ2IsotropicHardening not initialized");
-    return *ipvJ2IsotropicHardening;
-  }
   
   
 protected:

NonLinearSolver/materialLaw/CMakeLists.txt

@@ -68,7 +68,6 @@ set(SRC
   elasticPotential.cpp
   cohesiveElasticPotential.cpp
   mlawAnisotropicStoch.cpp
-  mlawNonLocalDamageGursonFullyCoupledThermoMechanics.cpp
   #  src/op_eshelby.cpp
    # Headers without cpp change this ??
   ID.h

NonLinearSolver/materialLaw/CoalescenceLaw.cpp

@@ -164,9 +164,17 @@ OriginalThomasonCoalescenceLaw::OriginalThomasonCoalescenceLaw(const OriginalTho
 void OriginalThomasonCoalescenceLaw::createIPVariable(IPCoalescence* &ipv) const
 {
   if(ipv != NULL) delete ipv;
+
   ipv = new IPOriginalThomasonCoalescence();
+  
+};
+
+void OriginalThomasonCoalescenceLaw::computeInitialChi(IPCoalescence* ipvCoal, double fVinitial) const
+{
+  ipvCoal->getRefToLigamentRatio() = pow(1.5* fVinitial *_initialVoidSpacingRatio,1./3.);
 };
 
+
 void OriginalThomasonCoalescenceLaw::checkCoalescence(const double fV, const IPNonLocalPorosity* ipv, const  IPCoalescence* ipvCoalPrev, IPCoalescence* ipvCoal) const{
 	const IPOriginalThomasonCoalescence* ipThomasonPrev = dynamic_cast<const IPOriginalThomasonCoalescence*>(ipvCoalPrev);
 	IPOriginalThomasonCoalescence* ipThomason = dynamic_cast<IPOriginalThomasonCoalescence*>(ipvCoal);

NonLinearSolver/materialLaw/CoalescenceLaw.h

@@ -46,6 +46,8 @@ class CoalescenceLaw
     virtual void createIPVariable(IPCoalescence* &ipv) const=0;
     virtual void checkCoalescence(const double fV, const IPNonLocalPorosity* ipv, const  IPCoalescence* ipvCoalPrev,  IPCoalescence* ipvCoal) const = 0;
 		virtual void forceCoalescence(const IPNonLocalPorosity* ipv, IPCoalescence* ipvCoal) = 0;
+    virtual void computeInitialChi(IPCoalescence* ipvCoal, double fVinitial) const = 0;
+
   
 #endif // SWIG
 };
@@ -72,6 +74,7 @@ class NoCoalescenceLaw : public CoalescenceLaw
 		virtual bool nonLocalCheck() const {return true;};
     virtual void checkCoalescence(const double fV, const IPNonLocalPorosity* ipv,const  IPCoalescence* ipvCoalPrev,  IPCoalescence* ipvCoal) const {};
 		virtual void forceCoalescence(const IPNonLocalPorosity* ipv, IPCoalescence* ipvCoal) {};
+    virtual void computeInitialChi(IPCoalescence* ipvCoal, double fVinitial) const {};
 #endif // SWIG
 };
 
@@ -102,6 +105,7 @@ class FstarCoalescenceLaw : public CoalescenceLaw
 		virtual bool nonLocalCheck() const {return _nonLocalCheck;};
     virtual void checkCoalescence(const double fV, const IPNonLocalPorosity* ipv, const  IPCoalescence* ipvCoalPrev,   IPCoalescence* ipvCoal) const;
 		virtual void forceCoalescence(const IPNonLocalPorosity* ipv, IPCoalescence* ipvCoal);
+    virtual void computeInitialChi(IPCoalescence* ipvCoal, double fVinitial) const {};
     #endif // SWIG
 };
 
@@ -133,6 +137,7 @@ class OriginalThomasonCoalescenceLaw : public CoalescenceLaw{
 		virtual bool nonLocalCheck() const {return _nonLocalCheck;};
     virtual void checkCoalescence(const double fV, const IPNonLocalPorosity* ipv,const  IPCoalescence* ipvCoalPrev, IPCoalescence* ipvCoal) const;
 		virtual void forceCoalescence(const IPNonLocalPorosity* ipv, IPCoalescence* ipvCoal);
+    virtual void computeInitialChi(IPCoalescence* ipvCoal, double fVinitial) const;
 		#endif //SWIG
 };
 

NonLinearSolver/materialLaw/j2IsotropicHardening.cpp

@@ -14,57 +14,11 @@
 J2IsotropicHardening::J2IsotropicHardening(const int num, double yield0, const bool init): _num(num), _initialized(init),
 											_yield0(yield0)
 {
-  _temFunc_Sy0= new constantScalarFunction(1.);
-  _temFunc_h= new constantScalarFunction(1.);
-  _temFunc_hexp= new constantScalarFunction(1.);
-  _temFunc_h1= new constantScalarFunction(1.);
-  _temFunc_h2= new constantScalarFunction(1.);
-  _temFunc_pexp= new constantScalarFunction(1.);
-  _temFunc_hexp2= new constantScalarFunction(1.);
-  _temFunc_K= new constantScalarFunction(1.);
-  _temFunc_p0= new constantScalarFunction(1.);
-  
   if(_yield0 < 0) Msg::Error("J2IsotropicHardening: negative yield stress");
 }
 
 J2IsotropicHardening::J2IsotropicHardening(const J2IsotropicHardening &source)
 {
-  _temFunc_Sy0 = NULL; 
-  if (source._temFunc_Sy0 != NULL){
-    _temFunc_Sy0 = source._temFunc_Sy0->clone();
-  }
-  _temFunc_h = NULL; 
-  if (source._temFunc_h != NULL){
-    _temFunc_h = source._temFunc_h->clone();
-  }
-  _temFunc_hexp = NULL; 
-  if (source._temFunc_hexp != NULL){
-    _temFunc_hexp = source._temFunc_hexp->clone();
-  }
-  _temFunc_h1 = NULL; 
-  if (source._temFunc_h1 != NULL){
-    _temFunc_h1 = source._temFunc_h1->clone();
-  }
-  _temFunc_h2 = NULL; 
-  if (source._temFunc_h2 != NULL){
-    _temFunc_h2 = source._temFunc_h2->clone();
-  }
-  _temFunc_pexp = NULL; 
-  if (source._temFunc_pexp != NULL){
-    _temFunc_pexp = source._temFunc_pexp->clone();
-  }
-  _temFunc_hexp2 = NULL; 
-  if (source._temFunc_hexp2 != NULL){
-    _temFunc_hexp2 = source._temFunc_hexp2->clone();
-  }
-  _temFunc_K = NULL; 
-  if (source._temFunc_K != NULL){
-    _temFunc_K = source._temFunc_K->clone();
-  }
-  _temFunc_p0 = NULL; 
-  if (source._temFunc_p0 != NULL){
-    _temFunc_p0 = source._temFunc_p0->clone();
-  }
   _num = source._num;
   _initialized = source._initialized;
   _yield0 = source._yield0;
@@ -72,44 +26,12 @@ J2IsotropicHardening::J2IsotropicHardening(const J2IsotropicHardening &source)
 
 J2IsotropicHardening& J2IsotropicHardening::operator=(const J2IsotropicHardening &source)
 {
-  _temFunc_Sy0 = NULL; 
-  if (source._temFunc_Sy0 != NULL){
-    _temFunc_Sy0 = source._temFunc_Sy0->clone();
-  }
   _num = source._num;
   _initialized = source._initialized;
   _yield0 = source._yield0;
   return *this;
 }
-//-----------------------added
-void J2IsotropicHardening::setTemperatureFunction_Sy0(const scalarFunction& Tfunc){
-      if (_temFunc_Sy0 != NULL) delete _temFunc_Sy0;
-      _temFunc_Sy0 = Tfunc.clone();}
-void J2IsotropicHardening::setTemperatureFunction_h(const scalarFunction& Tfunc){
-      if (_temFunc_h != NULL) delete _temFunc_h;
-      _temFunc_h = Tfunc.clone();}
-void J2IsotropicHardening::setTemperatureFunction_hexp(const scalarFunction& Tfunc){
-      if (_temFunc_hexp != NULL) delete _temFunc_hexp;
-      _temFunc_hexp = Tfunc.clone();}
-void J2IsotropicHardening::setTemperatureFunction_h1(const scalarFunction& Tfunc){
-      if (_temFunc_h1 != NULL) delete _temFunc_h1;
-      _temFunc_h1 = Tfunc.clone();}
-void J2IsotropicHardening::setTemperatureFunction_h2(const scalarFunction& Tfunc){
-      if (_temFunc_h2 != NULL) delete _temFunc_h2;
-      _temFunc_h2 = Tfunc.clone();}
-void J2IsotropicHardening::setTemperatureFunction_pexp(const scalarFunction& Tfunc){
-      if (_temFunc_pexp != NULL) delete _temFunc_pexp;
-      _temFunc_pexp = Tfunc.clone();}
-void J2IsotropicHardening::setTemperatureFunction_hexp2(const scalarFunction& Tfunc){
-      if (_temFunc_hexp2 != NULL) delete _temFunc_hexp2;
-      _temFunc_hexp2 = Tfunc.clone();}
-void J2IsotropicHardening::setTemperatureFunction_K(const scalarFunction& Tfunc){
-      if (_temFunc_K != NULL) delete _temFunc_K;
-      _temFunc_K = Tfunc.clone();}
-void J2IsotropicHardening::setTemperatureFunction_p0(const scalarFunction& Tfunc){
-      if (_temFunc_p0 != NULL) delete _temFunc_p0;
-      _temFunc_p0 = Tfunc.clone();}
-//------------------------------------------end
+
 PerfectlyPlasticJ2IsotropicHardening::PerfectlyPlasticJ2IsotropicHardening(const int num, double yield0, const bool init) :
 					J2IsotropicHardening(num,yield0,init)
 {
@@ -151,31 +73,10 @@ void PerfectlyPlasticJ2IsotropicHardening::hardening(double p, IPJ2IsotropicHard
   }
   ipv.set(R,dR,ddR,intR);
 }
-//-------------------------------------added
-void PerfectlyPlasticJ2IsotropicHardening::hardening(double p,double T, IPJ2IsotropicHardening &ipv) const
-{
-  double tol=1.e-15;
-  double dR=0, ddR=0, intR=0;
-  double R = getYield0(T);
-  if(p>0)
-  {
-    R=getYield0(T)+tol*p;
-    dR=tol;
-    ddR=0;
-    intR=0.5*tol*p*p+getYield0(T)*p;
-  }
-  ipv.set(R,dR,ddR,intR);
-}
-//-----------------------------------------end
 J2IsotropicHardening * PerfectlyPlasticJ2IsotropicHardening::clone() const
 {
   return new PerfectlyPlasticJ2IsotropicHardening(*this);
 }
-double PerfectlyPlasticJ2IsotropicHardening::DRDT(double p, double T) const
-{
-  double dRdT=getDYield0DT(T);
-  return dRdT;
-}
 
 PowerLawJ2IsotropicHardening::PowerLawJ2IsotropicHardening(const int num, double yield0, double h, double hexp, double pth, const bool init) :
 					J2IsotropicHardening(num,yield0,init), _h(h), _hexp(hexp), _pth(pth)
@@ -218,6 +119,17 @@ void PowerLawJ2IsotropicHardening::hardening(double p, IPJ2IsotropicHardening &i
   double R = getYield0();
   if(p< _pth)
   {
+  /*  R=getYield0();
+    if(_h<1.)
+    {
+      dR = 1e20;
+      ddR = -1e20;
+    }
+    else
+    {
+      dR = _h;
+      ddR = 0.;
+    }*/
    R=getYield0()+_h*pow(_pth,_hexp)/_pth*p;
    dR=_h*pow(_pth,_hexp)/_pth;
    ddR=0.;
@@ -234,67 +146,11 @@ void PowerLawJ2IsotropicHardening::hardening(double p, IPJ2IsotropicHardening &i
   }
   ipv.set(R,dR,ddR,intR);
 }
-//--------------------------------------------------added
-void PowerLawJ2IsotropicHardening::hardening(double p,double T, IPJ2IsotropicHardening &ipv) const
-{
-  double tol=1.e-16;
-  double dR=0, ddR=0, intR=0;
-  double R = getYield0(T);
-  double hT=_h*_temFunc_h->getVal(T);
-  double hexpT=_hexp*_temFunc_hexp->getVal(T);
-  if(p< tol)
-  {
-    R=getYield0(T);
-    if(hT<1.)
-    {
-      dR = 1e20;
-      ddR = -1e20;
-    }
-    else
-    {
-      dR = hT;
-      ddR = 0.;
-    }
-  }
-  else
-  {
-    R=getYield0(T)+hT*pow(p,hexpT);
-    dR = hT*hexpT*pow(p,hexpT)/p;
-    ddR = dR*(hexpT-1.)/p;
-    intR=getYield0(T)*p;
-    if(fabs(hexpT+1.)!=0)
-        intR+=hT*pow(p,hexpT+1.)/(hexpT+1.);
-  }
-  ipv.set(R,dR,ddR,intR);
-}
-//--------------------------------------------------------------------------end
 J2IsotropicHardening * PowerLawJ2IsotropicHardening::clone() const
 {
   return new PowerLawJ2IsotropicHardening(*this);
 }
-//--------------------------------added
-double PowerLawJ2IsotropicHardening::DRDT(double p, double T) const
-{
-  double tol=1.e-16;
-  double dhdT=_h*_temFunc_h->getDiff(T);
-  double dhexpdT=_hexp*_temFunc_hexp->getDiff(T);
-  double hT=_h*_temFunc_h->getVal(T);
-  double hexpT=_hexp*_temFunc_h->getVal(T);
-  double dRdT;
-  if(p<tol && dhexpdT != 0)
-  {
-    dRdT=-1e20;
-  }
-  else if (p<tol && dhexpdT == 0)
-  {
-    dRdT=getDYield0DT(T)+dhdT*pow(p,hexpT);
-  }
-  else
-  {
-    dRdT=getDYield0DT(T)+dhdT*pow(p,hexpT)+hT*dhexpdT*log(p)*pow(p,hexpT);
-  }
-}
-//-----------------------------------------end
+
 ExponentialJ2IsotropicHardening::ExponentialJ2IsotropicHardening(const int num, double yield0, double h, double hexp, const bool init) :
 					J2IsotropicHardening(num,yield0,init), _h(h), _hexp(hexp)
 {
@@ -347,42 +203,6 @@ void ExponentialJ2IsotropicHardening::hardening(double p, IPJ2IsotropicHardening
   }
   ipv.set(R,dR,ddR,intR);
 }
-//-----------------------------------------------added
-void ExponentialJ2IsotropicHardening::hardening(double p,double T, IPJ2IsotropicHardening &ipv) const
-{
-  double tol=1.e-16;
-  double dR=0, ddR=0, intR=0;
-  double R = getYield0(T);
-  double hT=_h*_temFunc_h->getVal(T);
-  double hexpT=_hexp*_temFunc_hexp->getVal(T);
-  if(p< tol)
-  {
-    R=getYield0(T);
-    dR = hT*hexpT;
-    ddR = -hexpT*(dR);
-  }
-  else
-  {
-    double tmp = exp(-hexpT*p);
-    R = getYield0(T)*+hT*(1.-tmp);
-    dR = hT*hexpT*tmp;
-    ddR = -hexpT*(dR);
-    intR=getYield0(T)*p+hT*p;
-    if(fabs(hexpT)!=0) intR+=hT*(exp(-hexpT*p)-exp(0))/(hexpT);
-  }
-  ipv.set(R,dR,ddR,intR);
-}
-
-double ExponentialJ2IsotropicHardening::DRDT(double p, double T) const
-{
-  double dhdT=_h*_temFunc_h->getDiff(T);
-  double dhexpdT=_hexp*_temFunc_hexp->getDiff(T);
-  double hexpT=_hexp*_temFunc_hexp->getVal(T);
-  double hT=_h*_temFunc_h->getVal(T);
-  double dRdT=getDYield0DT(T)+dhdT*(1-exp(-hexpT*p))+hT*(dhexpdT*p*exp(-hexpT*p));
-  return dRdT;
-}
-//----------------------------------------end
 J2IsotropicHardening * ExponentialJ2IsotropicHardening::clone() const
 {
   return new ExponentialJ2IsotropicHardening(*this);
@@ -439,47 +259,11 @@ void SwiftJ2IsotropicHardening::hardening(double p, IPJ2IsotropicHardening &ipv)
   }
   ipv.set(R,dR,ddR,intR);
 }
-//----------------------added
-void SwiftJ2IsotropicHardening::hardening(double p,double T, IPJ2IsotropicHardening &ipv) const
-{
-  double tol=1.e-16;
-  double dR=0, ddR=0, intR=0;
-  double R = getYield0(T);
-  double hT = _h*_temFunc_h->getVal(T);
-  double hexpT = _hexp*_temFunc_hexp->getVal(T);
-  if(p< tol)
-  {
-    R=getYield0(T);
-    dR = getYield0(T)*hT*hexpT;
-    ddR = dR*(hexpT-1.)*hT;
-  }
-  else
-  {
-    double tmp = 1.+hT*p;
-    R = getYield0(T)*pow(tmp,hexpT);
-    dR = getYield0(T)*hT*hexpT*pow(tmp,hexpT)/tmp;
-    ddR = (dR)*(hexpT-1.)*hT/tmp;
-    if(fabs(hexpT+1.)!=0 && fabs(hT)!=0.) intR=getYield0(T)*pow(1.+hT*p,hexpT+1.)/(hexpT+1.)/hT;
-  }
-  ipv.set(R,dR,ddR,intR);
-}
-
-double SwiftJ2IsotropicHardening::DRDT(double p,double T) const
-{
-  double dhdT=_h*_temFunc_h->getDiff(T);
-  double dhexpdT=_hexp*_temFunc_hexp->getDiff(T);
-  double hexpT=_hexp*_temFunc_hexp->getVal(T);
-  double hT=_h*_temFunc_h->getVal(T);
-  double dRdT=getDYield0DT(T)*pow(1+hT*p,hexpT)+getYield0(T)*(dhdT*hexpT/(1+hT*p)+log(1+hT*p)*dhexpdT)*pow(1+hT*p,hexpT);
-  return dRdT;
-}
-//------------------------------end
 J2IsotropicHardening * SwiftJ2IsotropicHardening::clone() const
 {
   return new SwiftJ2IsotropicHardening(*this);
 }
 
-//--------------------------------------------------------------------------------------------------------end
 LinearExponentialJ2IsotropicHardening::LinearExponentialJ2IsotropicHardening(const int num, double yield0, double h1,
 					double h2, double hexp, const bool init) :
 					J2IsotropicHardening(num,yield0,init), _h1(h1), _h2(h2), _hexp(hexp)
@@ -535,44 +319,7 @@ void LinearExponentialJ2IsotropicHardening::hardening(double p, IPJ2IsotropicHar
   }
   ipv.set(R,dR,ddR,intR);
 }
-//-----------------------------added
-void LinearExponentialJ2IsotropicHardening::hardening(double p,double T, IPJ2IsotropicHardening &ipv) const
-{
-  double tol=1.e-16;
-  double dR=0, ddR=0, intR=0;
-  double R = getYield0(T);
-  double h1T=_h1*_temFunc_h1->getVal(T);
-  double h2T=_h2*_temFunc_h2->getVal(T);
-  double hexpT=_hexp*_temFunc_hexp->getVal(T);
-  if(p< tol)
-  {
-    R=getYield0(T);
-    dR = h1T + h2T*hexpT;
-    ddR = -h2T*hexpT*hexpT;
-  }
-  else
-  {
-    double tmp = exp(-hexpT*p);
-    R = getYield0(T)+h1T*p + h2T*(1.-tmp);
-    dR = h1T + h2T*hexpT*tmp;
-    ddR = - h2T*hexpT*hexpT*tmp;
-    intR=getYield0(T)*p+ h1T*p*p/2.+h2T*p;
-    if(fabs(hexpT)!=0) intR+=h2T*(exp(-hexpT*p)-exp(0))/(hexpT);
-  }
-  ipv.set(R,dR,ddR,intR);
-}
-double LinearExponentialJ2IsotropicHardening::DRDT(double p,double T) const
-{
-  double dh1dT=_h1*_temFunc_h1->getDiff(T);
-  double dhexpdT=_hexp*_temFunc_hexp->getDiff(T);
-  double hexpT=_hexp*_temFunc_hexp->getVal(T);
-  double h1T=_h1*_temFunc_h1->getVal(T);
-  double dh2dT=_h2*_temFunc_h2->getDiff(T);
-  double h2T=_h2*_temFunc_h2->getVal(T);
-  double dRdT=getDYield0DT(T)+dh1dT*p+dh2dT*(1-exp(-hexpT*p))+h2T*(dhexpdT*p*exp(-hexpT*p));
-  return dRdT;
-}
-//----------------------------end
+
 J2IsotropicHardening * LinearExponentialJ2IsotropicHardening::clone() const
 {
   return new LinearExponentialJ2IsotropicHardening(*this);
@@ -647,62 +394,7 @@ void LinearFollowedByExponentialJ2IsotropicHardening::hardening(double p, IPJ2Is
   }
   ipv.set(R,dR,ddR,intR);
 }
-//-----------------------------added
-void LinearFollowedByExponentialJ2IsotropicHardening::hardening(double p,double T, IPJ2IsotropicHardening &ipv) const
-{
-  double R = getYield0(T);
-  double dR(0.), ddR(0.), intR(0.);
-  double h1T = _h1*_temFunc_h1->getVal(T);
-  double pexpT = _pexp*_temFunc_pexp->getVal(T);
-  double h2T = _h2*_temFunc_h2->getVal(T);
-  double hexp2T = _hexp2*_temFunc_hexp2->getVal(T);
-  if(p < 1.e-16) // Elastic case
-  {
-    dR = h1T;
-    ddR = 0.;
-  }
-  else if (p <= pexpT) // Plastic case: linear part
-  {
-    R += h1T*p;
-    dR = h1T;
-    ddR = 0.;
-    intR = getYield0(T)*p + 0.5*h1T*p*p;
-  }
-  else // Plastic case: exponentional (saturation) part 
-  {
-    double tmp = exp(-(p-pexpT)/hexp2T);
-    R += h1T*p + h2T*(1.-tmp);
-    dR = h1T + h2T * tmp /hexp2T;
-    ddR = -h2T * tmp /hexp2T /hexp2T;
-    
-    intR = getYield0(T)*p;
-    intR += 0.5*h1T*p*p;
-    intR += h2T*(p-pexpT) + h2T*(tmp-exp(0.)) *hexp2T;
-  }
-  ipv.set(R,dR,ddR,intR);
-}
-double LinearFollowedByExponentialJ2IsotropicHardening::DRDT(double p,double T) const
-{
-  double dh1dT=_h1*_temFunc_h1->getDiff(T);
-  double dhexp2dT=_hexp2*_temFunc_hexp2->getDiff(T);
-  double hexp2T=_hexp2*_temFunc_hexp2->getVal(T);
-  double h1T=_h1*_temFunc_h1->getVal(T);
-  double dh2dT=_h2*_temFunc_h2->getDiff(T);
-  double h2T=_h2*_temFunc_h2->getVal(T);
-  double pexpT=_pexp*_temFunc_pexp->getVal(T);
-  double dpexpdT=_pexp*_temFunc_pexp->getDiff(T);
-  double dRdT;
-  if (p<pexpT)
-  {
-     dRdT=getDYield0DT(T)+dh1dT*p;
-  }
-  else
-  {
-     dRdT=getDYield0DT(T)+dh1dT*p+dh2dT*(1-exp(-hexp2T*(p-pexpT)))+h2T*(dhexp2dT*(p-pexpT)-hexp2T*dpexpdT)*exp(-hexp2T*p);
-  }
-  return dRdT;
-}
-//----------------------------end
+
 J2IsotropicHardening* LinearFollowedByExponentialJ2IsotropicHardening::clone() const
 {
   return new LinearFollowedByExponentialJ2IsotropicHardening(*this);
@@ -787,70 +479,6 @@ void LinearFollowedByPowerLawJ2IsotropicHardening::hardening(double p, IPJ2Isotr
   }
   ipv.set(R,dR,ddR,intR);
 }
-//-----------------------------------------------------added
-void LinearFollowedByPowerLawJ2IsotropicHardening::hardening(double p,double T, IPJ2IsotropicHardening &ipv) const
-{
-  double tol=1.e-16;
-  double dR=0., ddR=0., intR=0.;
-  double R = getYield0(T);
-  double h1T=_h1*_temFunc_h1->getVal(T);
-  double pexpT=_pexp*_temFunc_pexp->getVal(T);
-  double hexpT=_hexp*_temFunc_hexp->getVal(T);
-  // elastic case
-  if(p < tol)
-  {
-    // why ????
-    if(h1T<1.)
-    {
-      dR = 1e20;
-      ddR = -1e20;
-    }
-    else
-    {
-      dR = h1T;
-      ddR = 0.;
-    }
-  }
-  // linear part 
-  else if (p <= pexpT)
-  {
-    R = getYield0(T) + h1T*p;
-    dR = h1T;
-    ddR = 0.;
-    intR = getYield0(T)*p + 0.5*h1T*p*p;
-  }
-  // power part
-  else
-  {
-    R = (getYield0(T)+h1T*pexpT) * pow(p/pexpT,hexpT);
-    dR = R*hexpT/p;
-    ddR = dR*(hexpT-1.)/p;
-    
-    intR = getYield0(T)*pexpT + 0.5*h1T*pexpT*pexpT;
-    intR += R*p /(hexpT+1.) -((getYield0(T)+h1T*pexpT)*pexpT/(hexpT+1.));
-  }
-  ipv.set(R,dR,ddR,intR);
-}
-double LinearFollowedByPowerLawJ2IsotropicHardening::DRDT(double p,double T) const
-{
-  double dh1dT=_h1*_temFunc_h1->getDiff(T);
-  double dhexpdT=_hexp*_temFunc_hexp->getDiff(T);
-  double hexpT=_hexp*_temFunc_hexp->getVal(T);
-  double h1T=_h1*_temFunc_h1->getVal(T);
-  double pexpT=_pexp*_temFunc_pexp->getVal(T);
-  double dpexpdT=_pexp*_temFunc_pexp->getDiff(T);
-  double dRdT;
-  if(p<pexpT)
-  {
-    dRdT=getDYield0DT(T)+dh1dT*p;
-  }
-  else
-  {
-    dRdT=getDYield0DT(T)*(1+h1T*pexpT)*pow(p/pexpT,hexpT)+getYield0(T)*(1+dh1dT*pexpT+h1T*dpexpdT)*pow(p/pexpT,hexpT)+getYield0(T)*(1+h1T*pexpT)*((-dpexpdT/pexpT)*hexpT+log(p/pexpT)*dhexpdT)*pow(p/    pexpT,hexpT);
-  }
-  return dRdT;
-}
-//--------------------------------end
 J2IsotropicHardening * LinearFollowedByPowerLawJ2IsotropicHardening::clone() const
 {
   return new LinearFollowedByPowerLawJ2IsotropicHardening(*this);
@@ -958,26 +586,7 @@ void TwoExpJ2IsotropicHaderning::hardening(double p, IPJ2IsotropicHardening &ipv
   double intR = _yield0*p +_K*(exp(p)+0.5*exp(-2.*p));
   ipv.set(R,dR,ddR,intR);
 }
-//-------------------------------added
-void TwoExpJ2IsotropicHaderning::hardening(double p,double T, IPJ2IsotropicHardening &ipv) const
-{
-  double yield0T=getYield0(T);
-  double KT=_K*_temFunc_K->getVal(T);
-  double R =yield0T+ KT*(exp(p)-exp(-2.*p));
 
-  double dR = KT*(exp(p)+2.*exp(-2.*p));
-  double ddR = KT*(exp(p)-4.*exp(-2.*p));
-  double intR = yield0T*p +KT*(exp(p)+0.5*exp(-2.*p));
-  ipv.set(R,dR,ddR,intR);
-}
-double TwoExpJ2IsotropicHaderning::DRDT(double p,double T) const
-{
-  double KT=_K*_temFunc_K->getVal(T);
-  double dKdT=_K*_temFunc_K->getDiff(T);
-  double dRdT=getDYield0DT(T)+dKdT*(exp(p)-exp(-2.*p));
-  return dRdT;
-}
-//-------------------------------------end
 J2IsotropicHardening * TwoExpJ2IsotropicHaderning::clone() const
 {
   return new TwoExpJ2IsotropicHaderning(*this);
@@ -1016,30 +625,7 @@ void TanhJ2IsotropicHardening::hardening(double p, IPJ2IsotropicHardening &ipv)
   double intR = _yield0*p +_K*_p0*log(cosh(p/_p0));
   ipv.set(R,dR,ddR,intR);
 }
-//------------------------------added
-void TanhJ2IsotropicHardening::hardening(double p,double T, IPJ2IsotropicHardening &ipv) const
-{ 
-  double p0T=_p0*_temFunc_p0->getVal(T);
-  double yield0T=getYield0(T);
-  double KT=_K*_temFunc_K->getVal(T);
-  double th = tanh(p/p0T);
-  double R =yield0T+ KT*th;
-  double dR = KT*(1.-th*th)/p0T;
-  double ddR = KT*2.*th*(th*th-1.)/(p0T*p0T);
-  double intR = yield0T*p +KT*p0T*log(cosh(p/p0T));
-  ipv.set(R,dR,ddR,intR);
-}
-double TanhJ2IsotropicHardening::DRDT(double p,double T) const
-{
-  double KT= _K*_temFunc_K->getVal(T);
-  double dKdT=_K*_temFunc_K->getDiff(T);
-  double p0T= _p0*_temFunc_p0->getVal(T);
-  double dp0dT=_p0*_temFunc_p0->getDiff(T);
-  double th=tanh(p/p0T);
-  double dRdT=getDYield0DT(T)+dKdT*tanh(p/p0T)+KT*(1-th*th)*(-p*dp0dT/p0T/p0T);
-  return dRdT;
-}
-//------------------------------------------end
+
 J2IsotropicHardening * TanhJ2IsotropicHardening::clone() const
 {
   return new TanhJ2IsotropicHardening(*this);

NonLinearSolver/materialLaw/j2IsotropicHardening.h

@@ -15,7 +15,6 @@
 #include "ipstate.h"
 #include "MElement.h"
 #include "ipHardening.h"
-#include "scalarFunction.h"
 #endif
 class J2IsotropicHardening{
  public :
@@ -28,17 +27,6 @@ class J2IsotropicHardening{
   int _num; // number of law (must be unique !)
   bool _initialized; // to initialize law
   double _yield0;
-//-----------------------------added
-  scalarFunction* _temFunc_Sy0;
-  scalarFunction* _temFunc_h;
-  scalarFunction* _temFunc_hexp;
-  scalarFunction* _temFunc_h1;
-  scalarFunction* _temFunc_h2;
-  scalarFunction* _temFunc_pexp;
-  scalarFunction* _temFunc_hexp2;
-  scalarFunction* _temFunc_K;
-  scalarFunction* _temFunc_p0;
-//-------------------------------end
  public:
   // constructor
 #ifndef SWIG
@@ -54,21 +42,9 @@ class J2IsotropicHardening{
   virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const=0;
   virtual double getYield0() const {return _yield0;}
   virtual J2IsotropicHardening * clone() const=0;
-  virtual double getYield0(double T) const {return _yield0*_temFunc_Sy0->getVal(T);}
-  virtual double getDYield0DT(double T) const {return _yield0*_temFunc_Sy0->getDiff(T);}
-  virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const=0;
-  virtual double DRDT(double p, double T) const=0;
-  void setTemperatureFunction_Sy0(const scalarFunction& Tfunc);
-  void setTemperatureFunction_h(const scalarFunction& Tfunc);
-  void setTemperatureFunction_hexp(const scalarFunction& Tfunc);
-  void setTemperatureFunction_h1(const scalarFunction& Tfunc);
-  void setTemperatureFunction_h2(const scalarFunction& Tfunc);
-  void setTemperatureFunction_pexp(const scalarFunction& Tfunc);
-  void setTemperatureFunction_hexp2(const scalarFunction& Tfunc);
-  void setTemperatureFunction_K(const scalarFunction& Tfunc);
-  void setTemperatureFunction_p0(const scalarFunction& Tfunc);
 #endif
 };
+
 class PerfectlyPlasticJ2IsotropicHardening : public J2IsotropicHardening
 {
  public :
@@ -82,10 +58,8 @@ class PerfectlyPlasticJ2IsotropicHardening : public J2IsotropicHardening
   virtual hardeningname getType() const{return J2IsotropicHardening::perfectlyPlasticJ2IsotropicHardening; };
   virtual void createIPVariable(IPJ2IsotropicHardening* &ipv) const;
   virtual void initLaws(const std::map<int,J2IsotropicHardening*> &maplaw) {}; //nothing now, we will see if we use the mapping
-  virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
+   virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
   virtual J2IsotropicHardening * clone() const;
-  virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const;
-  virtual double DRDT(double p, double T) const;
 #endif
 };
 
@@ -95,6 +69,7 @@ class PowerLawJ2IsotropicHardening : public J2IsotropicHardening
  // R = yield0 + (h pth^hexp)/pth p if p<=pth
  protected :
   double _h, _hexp, _pth;
+
  public:
   // constructor
   PowerLawJ2IsotropicHardening(const int num,  double yield0, double h, double hexp, double pth=1.e-16, bool init=true);
@@ -108,8 +83,6 @@ class PowerLawJ2IsotropicHardening : public J2IsotropicHardening
   virtual void initLaws(const std::map<int,J2IsotropicHardening*> &maplaw) {}; //nothing now, we will see if we use the mapping
   virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
   virtual J2IsotropicHardening * clone() const;
-  virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const;
-  virtual double DRDT(double p, double T) const;
 #endif
 };
 
@@ -118,6 +91,7 @@ class ExponentialJ2IsotropicHardening : public J2IsotropicHardening
  // R = yield0 +h * (1-exp (-hexp *p))
  protected :
   double _h, _hexp;
+
  public:
   // constructor
   ExponentialJ2IsotropicHardening(const int num,  double yield0, double h, double hexp, bool init=true);
@@ -131,8 +105,6 @@ class ExponentialJ2IsotropicHardening : public J2IsotropicHardening
   virtual void initLaws(const std::map<int,J2IsotropicHardening*> &maplaw) {}; //nothing now, we will see if we use the mapping
   virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
   virtual J2IsotropicHardening * clone() const;
-  virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const;
-  virtual double DRDT(double p, double T) const;
 #endif
 };
 
@@ -141,6 +113,7 @@ class SwiftJ2IsotropicHardening : public J2IsotropicHardening
  // R = yield0 * (1+ h*p)^hexp
  protected :
   double _h, _hexp;
+
  public:
   // constructor
   SwiftJ2IsotropicHardening(const int num,  double yield0, double h, double hexp, bool init=true);
@@ -154,16 +127,15 @@ class SwiftJ2IsotropicHardening : public J2IsotropicHardening
   virtual void initLaws(const std::map<int,J2IsotropicHardening*> &maplaw) {}; //nothing now, we will see if we use the mapping
   virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
   virtual J2IsotropicHardening * clone() const;
-  virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const;
-  virtual double DRDT(double p,double T) const;
 #endif
 };
-//---------------------------------------------------------------------------end
+
 class LinearExponentialJ2IsotropicHardening : public J2IsotropicHardening
 {
  // R = yield0 +h1 * p + h2* (1- exp(-hexp*p))
  protected :
   double _h1, _h2, _hexp;
+
  public:
   // constructor
   LinearExponentialJ2IsotropicHardening(const int num,  double yield0, double h1, double h2, double hexp, bool init=true);
@@ -177,8 +149,6 @@ class LinearExponentialJ2IsotropicHardening : public J2IsotropicHardening
   virtual void initLaws(const std::map<int,J2IsotropicHardening*> &maplaw) {}; //nothing now, we will see if we use the mapping
   virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
   virtual J2IsotropicHardening * clone() const;
-  virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const;
-  virtual double DRDT(double p,double T) const;
 #endif
 };
 
@@ -204,8 +174,6 @@ class LinearFollowedByExponentialJ2IsotropicHardening : public J2IsotropicHarden
   virtual void initLaws(const std::map<int,J2IsotropicHardening*> &maplaw) {}; //nothing now, we will see if we use the mapping
   virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
   virtual J2IsotropicHardening * clone() const;
-  virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const;
-  virtual double DRDT(double p,double T) const;
 #endif
 };
 
@@ -231,8 +199,6 @@ class LinearFollowedByPowerLawJ2IsotropicHardening : public J2IsotropicHardening
   virtual void initLaws(const std::map<int,J2IsotropicHardening*> &maplaw) {}; //nothing now, we will see if we use the mapping
   virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
   virtual J2IsotropicHardening * clone() const;
-  virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const;
-  virtual double DRDT(double p,double T) const;
 #endif
 };
 
@@ -257,8 +223,6 @@ class PolynomialJ2IsotropicHardening : public J2IsotropicHardening{
   virtual void initLaws(const std::map<int,J2IsotropicHardening*> &maplaw) {}; //nothing now, we will see if we use the mapping
   virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
   virtual J2IsotropicHardening * clone() const;
-  virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const {Msg::Error("temperature dependent hardening is not implemented for PolynomialJ2IsotropicHardening");} ;
-  virtual double DRDT(double p,double T) const{Msg::Error("DRDT is not implemented for PolynomialJ2IsotropicHardening");};
   #endif // SWIG
 };
 
@@ -280,8 +244,6 @@ class TwoExpJ2IsotropicHaderning : public J2IsotropicHardening{
     virtual void initLaws(const std::map<int,J2IsotropicHardening*> &maplaw) {}; //nothing now, we will see if we use the mapping
     virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
     virtual J2IsotropicHardening * clone() const;
-    virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const;
-    virtual double DRDT(double p,double T) const;
   #endif // SWIG
 };
 
@@ -303,8 +265,6 @@ class TanhJ2IsotropicHardening : public J2IsotropicHardening{
     virtual void initLaws(const std::map<int,J2IsotropicHardening*> &maplaw) {}; //nothing now, we will see if we use the mapping
     virtual void hardening(double p, IPJ2IsotropicHardening &ipv) const;
     virtual J2IsotropicHardening * clone() const;
-    virtual void hardening(double p, double T, IPJ2IsotropicHardening &ipv) const;
-    virtual double DRDT(double p,double T) const;
   #endif // SWIG
 };
 

NonLinearSolver/materialLaw/mlaw.h

@@ -30,7 +30,7 @@ class materialLaw{
                  LinearThermoMechanics,J2ThermoMechanics,SMP,PhenomenologicalSMP,LinearElecTherMech,AnIsotropicElecTherMech,
                  hyperelastic, powerYieldLaw, powerYieldLawWithFailure,nonLocalDamagePowerYieldHyper,
                  localDamagePowerYieldHyperWithFailure,nonLocalDamagePowerYieldHyperWithFailure,ElecSMP,
-                 ThermalConducter,AnIsotropicTherMech, localDamageJ2Hyper,linearElastic,nonLocalDamageGursonThermoMechanics};
+                 ThermalConducter,AnIsotropicTherMech, localDamageJ2Hyper,linearElastic};
 
 
  protected :

NonLinearSolver/materialLaw/mlawNonLocalDamageGursonFullyCoupledThermoMechanics.h

-//
-// C++ Interface: material law
-//
-// Description: j2 thermo elasto-plastic Gurson law
-
-// Author:  <P. Harik>, (C) 2017
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-#ifndef MLAWNONLOCALGURSONFULLYCOUPLEDTHERMOMECHANICS_H_
-#define MLAWNONLOCALGURSONFULLYCOUPLEDTHERMOMECHANICS_H_
-
-#include "mlawNonLocalDamageGurson.h"
-#include "scalarFunction.h"
-#include "ipNonLocalDamageGursonThermoMechanics.h"
-
-class mlawNonLocalDamageGursonFullyCoupledThermoMechanics : public mlawNonLocalDamageGurson
-{
-  #ifndef SWIG
-
-  protected:
- 
-
-    double _TaylorQuineyFactor; // for heat conversion
-
-    double _ThermalConductivity;
-    double _cp; // heat capacity per unit field
-    scalarFunction* _temFunc_ThermalConductivity;
-    scalarFunction* _temFunc_cp;
-
-    bool _thermalEstimationPreviousConfig;
-    //double _Tref;	
-    STensor3 linearK;
-    STensor3 Stiff_alphaDilatation;
-    #endif //SWIG
-
-  public:
-    mlawNonLocalDamageGursonFullyCoupledThermoMechanics(const int num,const double E,const double nu, const double rho, 
-                           const double q1, const double q2, const double q3, const double fVinitial,
-                           const J2IsotropicHardening &j2IH,const CLengthLaw &cLLaw,
-                           const double tol=1.e-8, const bool matrixbyPerturbation = false, const double pert = 1e-8);
-
-    mlawNonLocalDamageGursonFullyCoupledThermoMechanics(const mlawNonLocalDamageGursonFullyCoupledThermoMechanics& src);
-    virtual ~mlawNonLocalDamageGursonFullyCoupledThermoMechanics();
-		
-    virtual materialLaw* clone() const {return new mlawNonLocalDamageGursonFullyCoupledThermoMechanics(*this);};
-    virtual bool withEnergyDissipation() const {return true;};
-		
-    const STensor3& getConductivityTensor() const {return linearK;};
-    const STensor3& getStiff_alphaDilatation() const {return Stiff_alphaDilatation;};
-    
-    //void setReferenceTemperature(const double Tr){_Tref = Tr;};
- 
-    void setReferenceThermalExpansionCoefficient(const double al){ 
-			_alp = al;
-			Stiff_alphaDilatation*= 0.;
-			Stiff_alphaDilatation(0,0) = al;
-			Stiff_alphaDilatation(1,1) = al;
-			Stiff_alphaDilatation(2,2) = al;
-		};
-
-
-    void setReferenceThermalConductivity(const double KK){ 
-			_ThermalConductivity = KK;
-			linearK*= 0.;
-			linearK(0,0) = KK;
-			linearK(1,1) = KK;
-			linearK(2,2) = KK;
-		};
-    void setTemperatureFunction_ThermalConductivity(const scalarFunction& Tfunc){
-      if (_temFunc_ThermalConductivity != NULL) delete _temFunc_ThermalConductivity;
-      _temFunc_ThermalConductivity = Tfunc.clone();
-    }
-
-    void setReferenceCp(const double Cp) {_cp = Cp;};
-    void setTemperatureFunction_Cp(const scalarFunction& Tfunc){
-      if (_temFunc_cp != NULL) delete _temFunc_cp;
-      _temFunc_cp = Tfunc.clone();
-    }
-
-    void setThermalEstimationPreviousConfig(const bool flag){
-      _thermalEstimationPreviousConfig = flag;
-    }
-
-    void setTaylorQuineyFactor(const double f){_TaylorQuineyFactor = f;};
-
-    virtual matname getType() const{return materialLaw::nonLocalDamageGursonThermoMechanics;}
-    //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 initLaws(const std::map<int,materialLaw*> &maplaw){};
- 
-    virtual double getInitialExtraDofStoredEnergyPerUnitField() const {return _cp;};
-    virtual double getExtraDofStoredEnergyPerUnitField(double T) const {return _cp*_temFunc_cp->getVal(T);};
-
-    // predictor-corrector estimation
-    /*virtual void constitutive(const STensor3& F0,         // previous deformation gradient (input @ time n)
-                            const STensor3& Fn,         // current deformation gradient (input @ time n+1)
-                            STensor3 &P,                // current 1st Piola-Kirchhoff stress tensor (output)
-                            const IPNonLocalDamageGursonThermoMechanics*q0,       // array of previous internal variables
-                            IPNonLocalDamageGursonThermoMechanics *q1,             // current array of internal variable (in ipvcur on output),
-                            STensor43 &Tangent,         // tangents (output)
-                            const bool stiff            // if true compute the tangents
-                           ) const;*/
-    virtual void constitutive(const STensor3& F0,         // initial deformation gradient (input @ time n)
-                            const STensor3& F1,         // updated deformation gradient (input @ time n+1)
-                            STensor3 &P1,                // updated 1st Piola-Kirchhoff stress tensor (output)
-                            const IPNonLocalDamageGursonThermoMechanics *q0,       // array of initial internal variable
-                            IPNonLocalDamageGursonThermoMechanics *q1,             // updated array of internal variable (in ipvcur on output),
-                            STensor43 &Tangent,         // mechanical tangents (output)
-                            //STensor3  &dLocalCorrectedPorosityDStrain,
-                            //STensor3  &dStressDNonLocalCorrectedPorosity,
-                            //double    &dLocalCorrectedPorosityDNonLocalCorrectedPorosity,
-                            const double& T0, // previous temperature
-                            const double& T, // temperature
-                            const SVector3 &gradT0, // previoustemeprature gradient
-  	                    const SVector3 &gradT, // temeprature gradient
-                            SVector3 &fluxT, // temperature flux
-                            STensor3 &dPdT, // mechanical-thermal coupling
-                            STensor3 &dfluxTdgradT, // thermal tengent
-                            SVector3 &dfluxTdT,
-                            STensor33 &dfluxTdF, // thermal-mechanical coupling
-                            double &thermalSource,   // - cp*dTdt
-                            double &dthermalSourcedT, // thermal source
-                            double &dthermalSourcedtildefVstar,//-----added
-                            STensor3 &dthermalSourcedF,
-                            double& mechanicalSource, // mechanical source--> convert to heat
-                            double & dmechanicalSourcedT,
-                            double & dmechanicalSourcedtildefVstar,//----added
-                            STensor3 & dmechanicalSourceF,
-                            STensor3  &dLocalPorosityDStrain,
-                            STensor3  &dStressDNonLocalPorosity,
-                            double    &dLocalPorosityDNonLocalPorosity,
-                            double    &dLocalPorosityDT,
-                            const bool stiff
-                            ) const;
-
-  protected:
-    double deformationEnergy(const STensor3& Ce, const double& T) const;
-    //unifromiser avec Gurson-> ajouter ce qui est relatif a ftilde
-    void predictorCorector(const STensor3& F0,         // initial deformation gradient (input @ time n)
-                            const STensor3& F,         // updated deformation gradient (input @ time n+1)
-                            STensor3 &P,                // updated 1st Piola-Kirchhoff stress tensor (output)
-                            const IPNonLocalDamageGursonThermoMechanics *q0,       // array of initial internal variable
-                            IPNonLocalDamageGursonThermoMechanics *q,             // updated array of internal variable (in ipvcur on output),
-                            const double& T0, // previous temperature
-			                      const double& T, // temperature
-			                      const SVector3 &gradT0, // previous temeprature gradient
-  			                    const SVector3 &gradT, // temeprature gradient
-                            SVector3 &fluxT, // temperature flux)
-                            double &thermalSource,
-                            double& mechanicalSource
-                            ) const;
-
-    void predictorCorector(const STensor3& F0,         // initial deformation gradient (input @ time n)
-                            const STensor3& F1,         // updated deformation gradient (input @ time n+1)
-                            STensor3 &P1,                // updated 1st Piola-Kirchhoff stress tensor (output)
-                            const IPNonLocalDamageGursonThermoMechanics *q0,       // array of initial internal variable
-                            IPNonLocalDamageGursonThermoMechanics *q1,             // updated array of internal variable (in ipvcur on output),
-                            STensor43 &Tangent,         // mechanical tangents (output)
-                           // STensor3  &dLocalCorrectedPorosityDStrain,
-                            //STensor3  &dStressDNonLocalCorrectedPorosity,
-                            //double    &dLocalCorrectedPorosityDNonLocalCorrectedPorosity,
-                            STensor43 &dFpdF, // plastic tangent
-                            const double& T0, // previous temperature
-			                      const double& T, // temperature
-			                      const SVector3 &gradT0, // previoustemeprature gradient
-  			                    const SVector3 &gradT, // temeprature gradient
-                            SVector3 &fluxT, // temperature flux
-                            STensor3 &dPdT, // mechanical-thermal coupling
-                            STensor3 &dfluxTdgradT, // thermal tengent
-                            SVector3 &dfluxTdT,
-                            STensor33 &dfluxTdF, // thermal-mechanical coupling
-                            double &thermalSource,   // - cp*dTdt
-                            double &dthermalSourcedT, // thermal source
-                            double &dthermalSourcedtildefVstar,//-----added
-                            STensor3 &dthermalSourcedF,
-                            double& mechanicalSource, // mechanical source--> convert to heat
-                            double & dmechanicalSourcedT,
-                            double & dmechanicalSourcedtildefVstar,//----added
-                            STensor3 & dmechanicalSourceF,
-                            STensor3  &dLocalPorosityDStrain,
-                            STensor3  &dStressDNonLocalPorosity,
-                            double    &dLocalPorosityDNonLocalPorosity,
-                            double    &dLocalPorosityDT,
-                            const bool stiff
-                            ) const; // tangent is used by pointer, if NULL->no compute
-
-  virtual bool isThermomechanicallyCoupled() const {return true;}
- 
-};
-
-#endif // MLAWNONLOCALDAMGEGURSONFULLYCOUPLEDTHERMOMECHANICS_H_