Merge branch 'jl_gurs' into 'master'

Transfer jl_gurs updates to master

See merge request !36

NonLinearSolver/internalPoints/ipHardening.cpp

@@ -249,6 +249,38 @@ IPVariable * IPLinearFollowedByExponentialJ2IsotropicHardening::clone() const
 
 
 
+/* IPLinearFollowedByPowerLawJ2IsotropicHardening */
+IPLinearFollowedByPowerLawJ2IsotropicHardening::IPLinearFollowedByPowerLawJ2IsotropicHardening(): IPJ2IsotropicHardening()
+{
+}
+
+IPLinearFollowedByPowerLawJ2IsotropicHardening::IPLinearFollowedByPowerLawJ2IsotropicHardening(const IPLinearFollowedByPowerLawJ2IsotropicHardening &source) : IPJ2IsotropicHardening(source)
+{
+}
+
+IPLinearFollowedByPowerLawJ2IsotropicHardening &IPLinearFollowedByPowerLawJ2IsotropicHardening::operator=(const IPVariable &source)
+{
+  IPJ2IsotropicHardening::operator=(source);
+  const IPLinearFollowedByPowerLawJ2IsotropicHardening* src = dynamic_cast<const IPLinearFollowedByPowerLawJ2IsotropicHardening*>(&source);
+  if(src != NULL)
+  {
+  }
+  return *this;
+}
+
+void IPLinearFollowedByPowerLawJ2IsotropicHardening::restart()
+{
+  IPJ2IsotropicHardening::restart();
+}
+
+IPVariable * IPLinearFollowedByPowerLawJ2IsotropicHardening::clone() const
+{
+  return new IPLinearFollowedByPowerLawJ2IsotropicHardening(*this);
+}
+
+
+
+
 
 IPPolynomialJ2IsotropicHardening::IPPolynomialJ2IsotropicHardening(): IPJ2IsotropicHardening()
 {

NonLinearSolver/internalPoints/ipHardening.h

@@ -160,6 +160,7 @@ class IPLinearExponentialJ2IsotropicHardening : public IPJ2IsotropicHardening
 };
 
 
+
 /* IPLinearFollowedByExponentialJ2IsotropicHardening */
 class IPLinearFollowedByExponentialJ2IsotropicHardening : public IPJ2IsotropicHardening
 {
@@ -177,6 +178,23 @@ class IPLinearFollowedByExponentialJ2IsotropicHardening : public IPJ2IsotropicHa
 
 
 
+/* IPLinearFollowedByPowerLawJ2IsotropicHardening */
+class IPLinearFollowedByPowerLawJ2IsotropicHardening : public IPJ2IsotropicHardening
+{
+
+ protected:
+ public:
+  IPLinearFollowedByPowerLawJ2IsotropicHardening();
+  IPLinearFollowedByPowerLawJ2IsotropicHardening(const IPLinearFollowedByPowerLawJ2IsotropicHardening &source);
+  virtual IPLinearFollowedByPowerLawJ2IsotropicHardening &operator=(const IPVariable &source);
+  virtual ~IPLinearFollowedByPowerLawJ2IsotropicHardening(){}
+  virtual void restart();
+  virtual IPVariable *clone() const;
+
+};
+
+
+
 class IPPolynomialJ2IsotropicHardening : public IPJ2IsotropicHardening{
   public:
     IPPolynomialJ2IsotropicHardening();

NonLinearSolver/internalPoints/ipNonLocalDamageGurson.cpp

@@ -24,7 +24,7 @@ IPNonLocalDamageGurson::IPNonLocalDamageGurson() : IPVariableMechanics(),  _elas
   _Fp(0,0)=1.;
   _Fp(1,1)=1.;
   _Fp(2,2)=1.;
-  _DFvStarDStrain=0;
+  _DFvStarDStrain=0.;
   _logsqrtCe = 0.; 
 };
 
@@ -32,8 +32,11 @@ IPNonLocalDamageGurson::IPNonLocalDamageGurson(double fVinitial, const J2Isotrop
                                                  const CLengthLaw *cll,
                                                  const std::vector<GursonDamageNucleation*> *gdnLawContainer) :
                                                  IPVariableMechanics(),  _elasticEnergy(0.),
-                                                 _eplmatrix(0.), _fV(fVinitial),  _fVstar(fVinitial),
-                                                 _nldfVstar(fVinitial), _DFvStarDNldfVstar(1.),_dissipationBlocked(false)
+                                                 _eplmatrix(0.), _fV(fVinitial),
+                                                 _fVstar(fVinitial),_nldfVstar(fVinitial),
+                                                 _DFvStarDNldfVstar(1.),
+                                                 _dissipationBlocked(false),
+                                                 _NonLocalToLocal(false)
 {
   ipvJ2IsotropicHardening=NULL;
   if(j2IH ==NULL) Msg::Error("IPNonLocalDamageGurson::IPNonLocalDamageGurson has no j2IH");
@@ -61,7 +64,7 @@ IPNonLocalDamageGurson::IPNonLocalDamageGurson(double fVinitial, const J2Isotrop
   
   _logsqrtCe = 0.;
   
-  _DFvStarDStrain=0;
+  _DFvStarDStrain = 0.;
   
 };
 
@@ -93,6 +96,7 @@ IPNonLocalDamageGurson::IPNonLocalDamageGurson(const IPNonLocalDamageGurson &sou
   _DFvStarDStrain=source._DFvStarDStrain;
   _DFvStarDNldfVstar=source._DFvStarDNldfVstar;
   _dissipationBlocked = source._dissipationBlocked;
+  _NonLocalToLocal = source._NonLocalToLocal;
 }
 
 IPNonLocalDamageGurson& IPNonLocalDamageGurson::operator=(const IPVariable &source)
@@ -154,26 +158,12 @@ IPNonLocalDamageGurson& IPNonLocalDamageGurson::operator=(const IPVariable &sour
     _DFvStarDStrain=src->_DFvStarDStrain;
     _DFvStarDNldfVstar=src->_DFvStarDNldfVstar;
     _dissipationBlocked = src->_dissipationBlocked;
+    _NonLocalToLocal = src->_NonLocalToLocal;
   }
   return *this;
 }
 
 
-double IPNonLocalDamageGurson::defoEnergy() const
-{
-  return _elasticEnergy;
-}
-
-
-double IPNonLocalDamageGurson::plasticEnergy() const
-{
-  if(ipvJ2IsotropicHardening != NULL)
-    return ipvJ2IsotropicHardening->getIntegR();
-  else
-    return 0.;
-}
-
-
 void IPNonLocalDamageGurson::restart()
 {
   IPVariableMechanics::restart();
@@ -190,6 +180,7 @@ void IPNonLocalDamageGurson::restart()
   restartManager::restart(_logsqrtCe);
   restartManager::restart(_DFvStarDStrain);
   restartManager::restart(_dissipationBlocked);
+  restartManager::restart(_NonLocalToLocal);
   return;
 }
 

NonLinearSolver/internalPoints/ipNonLocalDamageGurson.h

@@ -25,9 +25,6 @@ class IPNonLocalDamageGurson : public IPVariableMechanics
   IPJ2IsotropicHardening* ipvJ2IsotropicHardening; // Hardening law
   std::vector<IPGursonNucleation*> ipvgdnContainer; // Nucleation law vector
   
-  // Energies
-  double _elasticEnergy; // elastic energy stored
-
   // Internal variables
   double _eplmatrix;     // equivalent plastic strain in the matrix (hatp)
   double _fV;            // local porosity
@@ -42,7 +39,14 @@ class IPNonLocalDamageGurson : public IPVariableMechanics
   double _DFvStarDNldfVstar;  //derivative of the local corrected porosity with the non local corrected porosity
   
   // Damage and transition managing
-  bool _dissipationBlocked;
+  bool _dissipationBlocked; // True if dissipation is blocked at the IPv
+  bool _NonLocalToLocal;    // True if transition from local to non-local allowed
+  
+  // Energies
+  double _elasticEnergy; // elastic energy stored
+    // plastic energy : stored in hardening law
+
+
 
 
  public:
@@ -65,10 +69,7 @@ class IPNonLocalDamageGurson : public IPVariableMechanics
     }
   }
 
-
   // General functions
-  virtual double defoEnergy() const;
-  virtual double plasticEnergy() const;
   virtual void restart();
   virtual IPVariable* clone() const {return new IPNonLocalDamageGurson(*this);};
   
@@ -80,21 +81,29 @@ class IPNonLocalDamageGurson : public IPVariableMechanics
   #endif // HAVE_MPI
   
  
+ 
   // Damage managing
   virtual void blockDissipation(const bool fl){_dissipationBlocked = fl;};
   virtual bool dissipationIsBlocked() const {return _dissipationBlocked;};
   
 
-  // Access functions - members values
+  
+  // Access functions - Energy
   virtual double &getRefToElasticEnergy() {return _elasticEnergy;};
+  virtual double defoEnergy() const { return _elasticEnergy;};
+  virtual double plasticEnergy() const
+  {
+    if(ipvJ2IsotropicHardening != NULL) {return ipvJ2IsotropicHardening->getIntegR();}
+    else {return 0.;}
+  };
+  // Those need to be implemented ???
+  virtual double damageEnergy() const {return 0.;}; // dissipation by damage
+  virtual double irreversibleEnergy() const {return 0.;} // irreversible energy for path following
 
-  virtual double getDamage() const{return getLocalCorrectedPorosity();};
 
-  virtual const STensor3& getConstRefToFp() const {return _Fp;}
-  virtual STensor3& getRefToFp(){return _Fp;}
 
-  virtual const STensor3& getConstRefToLogsqrtCe() const{return _logsqrtCe;}
-  virtual STensor3& getRefToLogsqrtCe(){return _logsqrtCe;}
+  // Access functions - Internal variables
+  virtual double getDamage() const{return getLocalCorrectedPorosity();};
   
   virtual double getMatrixPlasticStrain() const { return _eplmatrix;}
   virtual double &getRefToMatrixPlasticStrain() { return _eplmatrix;}
@@ -108,6 +117,15 @@ class IPNonLocalDamageGurson : public IPVariableMechanics
   virtual double getNonLocalCorrectedPorosity() const { return _nldfVstar;}
   virtual double &getRefToNonLocalCorrectedPorosity() { return _nldfVstar;}
   
+  virtual const STensor3& getConstRefToFp() const {return _Fp;}
+  virtual STensor3& getRefToFp(){return _Fp;}
+  
+  virtual const STensor3& getConstRefToLogsqrtCe() const{return _logsqrtCe;}
+  virtual STensor3& getRefToLogsqrtCe(){return _logsqrtCe;}
+  
+  
+  
+  // Access functions - Derivatives
   virtual double &getRefToDLocalCorrectedPorosityDNonLocalCorrectedPorosity(){return _DFvStarDNldfVstar;}
   virtual double getDLocalCorrectedPorosityDNonLocalCorrectedPorosity() const{return _DFvStarDNldfVstar;}
   
@@ -115,6 +133,7 @@ class IPNonLocalDamageGurson : public IPVariableMechanics
   virtual STensor3 & getRefToDLocalCorrectedPorosityDStrain(){return _DFvStarDStrain;}
   
   
+  
   // Access functions - Included IPVs
   virtual const IPJ2IsotropicHardening &getConstRefToIPJ2IsotropicHardening() const
   {
@@ -176,7 +195,6 @@ class IPNonLocalDamageGurson : public IPVariableMechanics
   {
     return ipvgdnContainer;
   }
-  
 };
 
 #endif // IPNONLOCALDAMAGEGURSON_H_

NonLinearSolver/materialLaw/j2IsotropicHardening.cpp

@@ -315,7 +315,7 @@ J2IsotropicHardening * LinearExponentialJ2IsotropicHardening::clone() const
 
 
 
-/* LinearExponentialJ2IsotropicHardening */
+/* LinearFollowedByExponentialJ2IsotropicHardening */
 LinearFollowedByExponentialJ2IsotropicHardening::LinearFollowedByExponentialJ2IsotropicHardening(const int num,  double yield0,
     double h1, double pexp, double h2, double hexp2, bool init) :
     J2IsotropicHardening(num,yield0,init), _h1(h1), _pexp(pexp), _h2(h2), _hexp2(hexp2)
@@ -361,7 +361,7 @@ void LinearFollowedByExponentialJ2IsotropicHardening::hardening(double p, IPJ2Is
     dR = _h1;
     ddR = 0.;
   }
-  else if (p < _pexp) // Plastic case: linear part
+  else if (p <= _pexp) // Plastic case: linear part
   {
     R += _h1*p;
     dR = _h1;
@@ -371,13 +371,12 @@ void LinearFollowedByExponentialJ2IsotropicHardening::hardening(double p, IPJ2Is
   else // Plastic case: exponentional (saturation) part 
   {
     double tmp = exp(-(p-_pexp)/_hexp2);
-    R += _h1*_pexp + _h2*(1.-tmp);
-    dR = _h2 * tmp /_hexp2;
-    ddR = - dR /_hexp2;
-    // ddR = _h2*exp(-(p-_pexp)/_hexp) * (p-_pexp)/_hexp * -(p-_pexp)/_hexp + _h2*exp(-(p-_pexp)/_hexp) * p/_hexp;
+    R += _h1*p + _h2*(1.-tmp);
+    dR = _h1 + _h2 * tmp /_hexp2;
+    ddR = -_h2 * tmp /_hexp2 /_hexp2;
     
     intR = getYield0()*p;
-    intR += 0.5*_h1*_pexp*_pexp + _h1*_pexp*(p-_pexp);
+    intR += 0.5*_h1*p*p;
     intR += _h2*(p-_pexp) + _h2*(tmp-exp(0.)) *_hexp2;
   }
   ipv.set(R,dR,ddR,intR);
@@ -391,6 +390,89 @@ J2IsotropicHardening* LinearFollowedByExponentialJ2IsotropicHardening::clone() c
 
 
 
+
+/* LinearFollowedByPowerLawJ2IsotropicHardening */
+LinearFollowedByPowerLawJ2IsotropicHardening::LinearFollowedByPowerLawJ2IsotropicHardening(const int num,  double yield0, double h, double pexp, double hexp, bool init):
+        J2IsotropicHardening(num,yield0,init), _h1(h), _pexp(pexp), _hexp(hexp)
+{
+  if(_h1 < 0. or _hexp < 0.) Msg::Error("LinearFollowedByPowerLawJ2IsotropicHardening:: negative hardening parameters");
+  if(_pexp <= 1.e-9 ) Msg::Error("LinearFollowedByPowerLawJ2IsotropicHardening:: negative or too small hardening parameters p_exp, use PowerLawJ2IsotropicHardening instead");
+}
+
+LinearFollowedByPowerLawJ2IsotropicHardening::LinearFollowedByPowerLawJ2IsotropicHardening(const LinearFollowedByPowerLawJ2IsotropicHardening &source) :
+					J2IsotropicHardening(source)
+{
+  _h1 = source._h1;
+  _pexp = source._pexp;
+  _hexp = source._hexp;
+}
+
+LinearFollowedByPowerLawJ2IsotropicHardening& LinearFollowedByPowerLawJ2IsotropicHardening::operator=(const J2IsotropicHardening &source)
+{
+  J2IsotropicHardening::operator=(source);
+  const LinearFollowedByPowerLawJ2IsotropicHardening* src = dynamic_cast<const LinearFollowedByPowerLawJ2IsotropicHardening*>(&source);
+  if(src != NULL)
+  {
+    _h1 = src->_h1;
+    _pexp = src->_pexp;
+    _hexp = src->_hexp;
+  }
+  return *this;
+}
+void LinearFollowedByPowerLawJ2IsotropicHardening::createIPVariable(IPJ2IsotropicHardening* &ipv) const
+{
+  if(ipv != NULL) delete ipv;
+  ipv = new IPLinearFollowedByPowerLawJ2IsotropicHardening();
+}
+
+void LinearFollowedByPowerLawJ2IsotropicHardening::hardening(double p, IPJ2IsotropicHardening &ipv) const
+{
+  double tol=1.e-16;
+  double dR=0., ddR=0., intR=0.;
+  double R = getYield0();
+  // elastic case
+  if(p < tol)
+  {
+    // why ????
+    if(_h1<1.)
+    {
+      dR = 1e20;
+      ddR = -1e20;
+    }
+    else
+    {
+      dR = _h1;
+      ddR = 0.;
+    }
+  }
+  // linear part 
+  else if (p <= _pexp)
+  {
+    R = getYield0() + _h1*p;
+    dR = _h1;
+    ddR = 0.;
+    intR = getYield0()*p + 0.5*_h1*p*p;
+  }
+  // power part
+  else
+  {
+    R = (getYield0()+_h1*_pexp) * pow(p/_pexp,_hexp);
+    dR = R*_hexp/p;
+    ddR = dR*(_hexp-1.)/p;
+    
+    intR = getYield0()*_pexp + 0.5*_h1*_pexp*_pexp;
+    intR += R*p /(_hexp+1.) -((getYield0()+_h1*_pexp)*_pexp/(_hexp+1.));
+  }
+  ipv.set(R,dR,ddR,intR);
+}
+J2IsotropicHardening * LinearFollowedByPowerLawJ2IsotropicHardening::clone() const
+{
+  return new LinearFollowedByPowerLawJ2IsotropicHardening(*this);
+}
+
+
+
+/* PolynomialJ2IsotropicHardening*/
 PolynomialJ2IsotropicHardening::PolynomialJ2IsotropicHardening(const int num, double yield0, int order,  bool init):
 J2IsotropicHardening(num,yield0,init),_order(order){
   _coefficients.resize(_order+1);

NonLinearSolver/materialLaw/j2IsotropicHardening.h

@@ -20,7 +20,8 @@ class J2IsotropicHardening{
  public :
   enum hardeningname{perfectlyPlasticJ2IsotropicHardening,powerLawJ2IsotropicHardening, exponentialJ2IsotropicHardening,
                      swiftJ2IsotropicHardening, linearExponentialJ2IsotropicHardening,
-                     linearFollowedByExponentialJ2IsotropicHardening, polynomialJ2IsotropicHardening,
+                     linearFollowedByExponentialJ2IsotropicHardening, linearFollowedByPowerLawJ2IsotropicHardening,
+                     polynomialJ2IsotropicHardening,
                      twoExpJ2IsotropicHaderning, tanhJ2IsotropicHaderning};
  protected :
   int _num; // number of law (must be unique !)
@@ -155,7 +156,7 @@ class LinearFollowedByExponentialJ2IsotropicHardening : public J2IsotropicHarden
 {
  // R = yield0 + h1 * p
  // Then when p > p_exp 
- // R = yield0 + h1 * pexp + h2* (1- exp(-hexp2*(p-pexp)))
+ // R = yield0 + h1 * p + h2* (1- exp(-hexp2*(p-pexp)))
  protected :
   double _h1, _pexp, _h2, _hexp2;
 
@@ -176,6 +177,31 @@ class LinearFollowedByExponentialJ2IsotropicHardening : public J2IsotropicHarden
 };
 
 
+class LinearFollowedByPowerLawJ2IsotropicHardening : public J2IsotropicHardening
+{
+ // R = yield0 + h1 * p
+ // Then when p > pexp 
+ // R = yield0 * (1 + h1 * pexp) * (p/pexp) ^ h_exp
+ protected :
+  double _h1, _pexp, _hexp;
+
+ public:
+  // constructor
+  LinearFollowedByPowerLawJ2IsotropicHardening(const int num,  double yield0, double h, double pexp, double hexp, bool init=true);
+#ifndef SWIG
+  virtual ~LinearFollowedByPowerLawJ2IsotropicHardening(){}
+  LinearFollowedByPowerLawJ2IsotropicHardening(const LinearFollowedByPowerLawJ2IsotropicHardening &source);
+  LinearFollowedByPowerLawJ2IsotropicHardening& operator=(const J2IsotropicHardening &source);
+  virtual int getNum() const{return _num;}
+  virtual hardeningname getType() const{return J2IsotropicHardening::linearFollowedByPowerLawJ2IsotropicHardening; };
+  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 J2IsotropicHardening * clone() const;
+#endif
+};
+
+
 class PolynomialJ2IsotropicHardening : public J2IsotropicHardening{
   #ifndef SWIG
   protected: