Merge branch 'vdg-cm3'

NonLinearSolver/materialLaw/failureCriterion.cpp

@@ -140,8 +140,134 @@ bool CritialPlasticDeformationCriteration::isFailed(const IPVariable* ipvprev, c
   else{
     double Phif = this->getFailureCriterion(ipvprev,ipv);
     double r = viscoipv->getFailureOnset();  
-    double Ff =  Phif  - r;
-    if (Ff >=0) return true;
+    if (Phif >= r) return true;
+    else return false;
+  }
+};
+
+void ThreeCriticalStressCriterion::getFactor(const double sigt, const double sigc, const double sigs, double& aa, double& bb, double& cc) const{
+  double pc = -sigc/3.;
+  double pt = sigt/3.;
+  
+  double mat[2][2];
+  mat[0][0] = pc;
+  mat[0][1] = sigc-sigs;
+  mat[1][0] = pt;
+  mat[1][1] = sigt-sigs;
+  
+  double det = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
+  double inv[2][2];
+  if(fabs(det) > 0){
+    double ud = 1. / det;
+    inv[0][0] =  mat[1][1] * ud;
+    inv[1][0] = -mat[1][0] * ud;
+    inv[0][1] = -mat[0][1] * ud;
+    inv[1][1] =  mat[0][0] * ud;
+  }
+  else{
+    Msg::Fatal("Singular matrix 3x3 ThreeCriticalStressCriterion::ThreeCriticalStressCriterion");
+  }
+  
+  double vec[2];
+  vec[0] = (sigs-sigc)*pc;
+  vec[1] = (sigs-sigt)*pt;
+  
+  double sol[2];
+  sol[0] = inv[0][0]*vec[0]+inv[0][1]*vec[1];
+  sol[1] = inv[1][0]*vec[0]+inv[1][1]*vec[1];
+  
+  // compute factor a, b, c
+  aa = sol[0]*sol[1];
+  bb = sol[1];
+  cc = sigs - aa/bb;
+  
+  //printf(" aa = %e bb = %e cc = %e\n",aa,bb,cc);
+};
+
+ThreeCriticalStressCriterion::ThreeCriticalStressCriterion(const int num, const RateFailureLaw& trac, const RateFailureLaw& comp, const RateFailureLaw& shear): FailureCriterionBase(num){
+  _compFailure = comp.clone();
+  _tracFailure = trac.clone();
+  _shearFailure = shear.clone();
+  _isConstant = (comp.isConstant() && trac.isConstant() && shear.isConstant());
+  // esitmated a,b, c
+  if (_isConstant){
+    double sigc,sigt,sigs,dR;
+    _compFailure->get(0,sigc,dR);
+    _tracFailure->get(0,sigt,dR);
+    _shearFailure->get(0,sigs,dR);
+  
+    getFactor(sigt,sigc,sigs,a,b,c);
+  }
+};
+
+  
+ThreeCriticalStressCriterion::ThreeCriticalStressCriterion(const ThreeCriticalStressCriterion& src): FailureCriterionBase(src), a(src.a),b(src.b),c(src.c), _isConstant(src._isConstant){
+  _compFailure = NULL;
+  if (src._compFailure){
+    _compFailure= src._compFailure->clone();
+  }
+  _tracFailure = NULL;
+  if (src._tracFailure){
+    _tracFailure = src._tracFailure->clone();
+  }
+  _shearFailure = NULL;
+  if (src._shearFailure){
+    _shearFailure = src._shearFailure->clone();
+  }
+};
+
+ThreeCriticalStressCriterion::~ThreeCriticalStressCriterion(){
+  if (_compFailure){
+    delete _compFailure;
+    _compFailure = NULL;
+  }
+  if (_tracFailure){
+    delete _tracFailure;
+    _tracFailure = NULL;
+  }
+  if (_shearFailure){
+    delete _shearFailure;
+    _shearFailure = NULL;
+  }
+}
+
+double ThreeCriticalStressCriterion::getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const{
+  const IPHyperViscoElastoPlastic* viscoipv = dynamic_cast<const IPHyperViscoElastoPlastic*>(ipv);
+  if (viscoipv == NULL){
+    Msg::Error("PowerFailureCriterion::isFailed is implemented for IPHyperViscoElastoPlastic only");
+    return 1e10;
+  }
+  else{
+    if (!_isConstant){
+      double sigc,sigt,sigs,dR;
+      const IPHyperViscoElastoPlastic* viscoipvprev = static_cast<const IPHyperViscoElastoPlastic*>(ipvprev);
+      _compFailure->get(viscoipvprev->_DgammaDt,sigc,dR);
+      _tracFailure->get(viscoipvprev->_DgammaDt,sigt,dR);
+      _shearFailure->get(viscoipvprev->_DgammaDt,sigs,dR);
+      getFactor(sigt,sigc,sigs,a,b,c);
+    }
+    
+    // compute failure r and DrDF
+    const STensor3& K = viscoipv->_kirchhoff;
+    double pres = K.trace()/3.;
+    STensor3 devK = K.dev();
+    double Keq = sqrt(1.5*devK.dotprod());
+    // get coeficinet
+    //K.print("Kir stess");
+    //Msg::Info(" p = %e  KEq = %e, cr = %e",pres,Keq,Keq - a/(pres+b) -c);
+    return (Keq - a/(pres+b) -c);
+  }
+};
+bool ThreeCriticalStressCriterion::isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const{
+  const IPHyperViscoElastoPlastic* viscoipv = dynamic_cast<const IPHyperViscoElastoPlastic*>(ipv);
+  if (viscoipv == NULL){
+    Msg::Error("PowerFailureCriterion::isFailed is implemented for IPHyperViscoElastoPlastic only");
+    return false;
+  }
+  else{
+    double Phif = this->getFailureCriterion(ipvprev,ipv);
+    double r = viscoipv->getFailureOnset();  
+    if (Phif >= r) return true;
     else return false;
   }
 };
\ No newline at end of file

NonLinearSolver/materialLaw/failureCriterion.h

@@ -97,4 +97,29 @@ class CritialPlasticDeformationCriteration : public FailureCriterionBase{
     #endif //SWIG
   
 };
+
+class ThreeCriticalStressCriterion : public FailureCriterionBase{
+   protected:
+    #ifndef SWIG
+    RateFailureLaw* _compFailure;
+    RateFailureLaw* _tracFailure;
+    RateFailureLaw* _shearFailure;
+    bool _isConstant;
+    mutable double a, b, c; // cache for surface sigVM = c + a/(p+b) , a,b,c estimated from sigc, sigt, sigs
+    
+private:
+    void getFactor(const double sigt, const double sigc, const double sigs, double& aa, double& bb, double& cc) const;
+    #endif //SWIG
+ 
+  public:
+    ThreeCriticalStressCriterion(const int num, const RateFailureLaw& trac, const RateFailureLaw& comp, const RateFailureLaw& shear);
+    #ifndef SWIG
+    ThreeCriticalStressCriterion(const ThreeCriticalStressCriterion& src);
+    virtual ~ThreeCriticalStressCriterion();
+    virtual double getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const;
+    virtual bool isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const;
+    virtual FailureCriterionBase* clone() const {return new ThreeCriticalStressCriterion(*this);};
+    #endif //SWIG
+};
+
 #endif // FAILURECRITERION_H_

NonLinearSolver/materialLaw/materialStrengthLaw.h

@@ -26,6 +26,7 @@ class RateFailureLaw {
     virtual double getRefStrength() const = 0;
     virtual void get(const double rate, double& R, double& dR) const = 0;
     virtual RateFailureLaw* clone() const  = 0;
+    virtual bool isConstant() const = 0;
     #endif // SWIG
 };
 
@@ -51,7 +52,7 @@ class ConstantFailureLaw : public RateFailureLaw{
       dR = 0.;
     };
     virtual RateFailureLaw* clone() const  {return new ConstantFailureLaw(*this);};
-
+    virtual bool isConstant() const {return true;};
     #endif // SWIG
 };
 
@@ -69,6 +70,7 @@ class PowerFailureLaw : public RateFailureLaw{
     virtual double getRefStrength() const{return _X0;}
     virtual void get(const double rate, double& R, double& dR) const ;
     virtual RateFailureLaw* clone() const  {return new PowerFailureLaw(*this);};
+    virtual bool isConstant() const {return false;};
     #endif // SWIG
 };
 
@@ -88,6 +90,7 @@ class LograrithmicFailureLaw : public RateFailureLaw{
     virtual double getRefStrength() const{return _Xr;}
     virtual void get(const double rate, double& R, double& dR) const ;
     virtual RateFailureLaw* clone() const  {return new LograrithmicFailureLaw(*this);};
+    virtual bool isConstant() const {return false;};
     #endif // SWIG
 };