diff --git a/NonLinearSolver/materialLaw/failureCriterion.cpp b/NonLinearSolver/materialLaw/failureCriterion.cpp
index d77d7e27c9a40db6ddafd84a1d5c3735528b3f54..1513e7b2cce6a626a1a812851a7666e15431d87a 100644
--- a/NonLinearSolver/materialLaw/failureCriterion.cpp
+++ b/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
diff --git a/NonLinearSolver/materialLaw/failureCriterion.h b/NonLinearSolver/materialLaw/failureCriterion.h
index 4687df58678f9293df0e05f8f6a74282bea6256f..fcf979c6eabf53d5dd6183e97fc15263a88c39bd 100644
--- a/NonLinearSolver/materialLaw/failureCriterion.h
+++ b/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_
diff --git a/NonLinearSolver/materialLaw/materialStrengthLaw.h b/NonLinearSolver/materialLaw/materialStrengthLaw.h
index 49df951467885ea8c2614d76d21db913d3e4d944..6b68992cc9374432105728dfbdbd2dfa9b738474 100644
--- a/NonLinearSolver/materialLaw/materialStrengthLaw.h
+++ b/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
};