add new hardening law - linear part followed by power law

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/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)
@@ -390,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 !)
@@ -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: