TPZYCModifiedMohrCoulomb.h

Go to the documentation of this file.

/* Generated by Together */// $Id: TPZYCModifiedMohrCoulomb.h,v 1.9 2010-12-13 19:34:58 diogo Exp $

#ifndef TPZYCModifiedMohrCoulomb_H
#define TPZYCModifiedMohrCoulomb_H

#include "TPZTensor.h"
#include "pzlog.h"
#include "TPZPlasticCriterion.h"
#ifndef WIN32
#include <fenv.h>//NAN DETECTOR
#endif

#ifdef LOG4CXX_PLASTICITY
//static LoggerPtr logMohr(Logger::getLogger("TPZYCMohrOriginal"));
#endif

class TPZYCModifiedMohrCoulomb : public TPZPlasticCriterion {
    
public:
    
    enum {NYield = 1};
    
    virtual int ClassId() const override;

    
    const char * Name() const
    {
        return "TPZYCModifiedMohrCoulomb";
    }
    
    void Print(std::ostream & out) const override
    {
        out << Name();
    }
    
    int GetForceYield()
    {
        return 0; // nothing to be done in this yield criterium
    }
    
    void SetForceYield(const int forceYield)
    {
        // nothing to be done in this yield criterium
    }
    
    void Write(TPZStream& buf, int withclassid) const override{
        buf.Write(&fPhi);
        buf.Write(&fCoesion);
        buf.Write(&fPi);
    }
    
    void Read(TPZStream& buf, void* context) override{
        buf.Read(&fPhi);
        buf.Read(&fCoesion);
        buf.Read(&fPi);
    }
    
    void SetYieldStatusMode(const TPZTensor<REAL> & sigma, const REAL & A)
    {
        // nothing to be done in this yield criterium
    }
    
    
    void SetUp(const REAL & phi)
    {
        fPhi = phi;
        fPi = M_PI;
        fCoesion = 9.2376;
    }
    
    
    template < class T>
    void Compute(const TPZTensor<T> & sigma, const T & A, TPZVec<T> &result, int checkForcedYield = 0) const;
    
    template <class T>
    void N(const TPZTensor<T> & sigma,const T & A,  TPZVec<TPZTensor<T> > & Ndir, int checkForcedYield = 0) const;
    
    template <class T>
    void H(const TPZTensor<T> & sigma,const T & A,  TPZVec<T> & h, int checkForcedYield = 0) const;

    template <class T>
    void AlphaMultiplier(const T &A, T &multiplier) const
    {
        multiplier = T(1.);
    }

    void YieldFunction(const TPZVec<STATE>& sigma, STATE kprev, TPZVec<STATE>& yield) const override{
        TPZTensor<STATE> sigmaTensor;
        sigmaTensor.XX() = sigma[0];
        sigmaTensor.YY() = sigma[1];
        sigmaTensor.ZZ() = sigma[2];
        Compute(sigmaTensor, kprev, yield, 0);
    }

    virtual int GetNYield() const override{
        return as_integer(NYield);
    }
    
public:
    
    REAL fPhi;
    
    
protected:
    
    REAL fCoesion,fPi;
public:

    
};

static bool fFlag,fFlag1,fFlag2;

template < class T>
void TPZYCModifiedMohrCoulomb::Compute(const TPZTensor<T> & sigma, const T & A,TPZVec<T> &res, int checkForcedYield) const
{
    
    TPZTensor<T> eigenval,dSigma1,dSigma2,dSigma3;
    sigma.Eigenvalue(eigenval, dSigma1, dSigma2, dSigma3);
    T sigma1,sigma2,sigma3;
    sigma1 = eigenval.XX();
    sigma2 = eigenval.YY();
    sigma3 = eigenval.ZZ();
    T res0,res1,res2;
    //
   // T debug = T(2.)*A*cos(fPhi);
    res0=(sigma1 - sigma3) + (sigma1 + sigma3)*sin(fPhi) - T(2.)*A*cos(fPhi);
    res1=(sigma2 - sigma3) + (sigma2 + sigma3)*sin(fPhi) - T(2.)*A*cos(fPhi);
    res2=(sigma1 - sigma2) + (sigma1 + sigma2)*sin(fPhi) - T(2.)*A*cos(fPhi);
   // cout << " res0 " <<res0 << endl;

    fFlag=false;
    fFlag1=false;
    fFlag2=false;
    
    if(TPZExtractVal::val(res0) >= 0. && TPZExtractVal::val(res1) < 0. && TPZExtractVal::val(res2) < 0.)
    {
        res[0]=(sigma1 - sigma3) + (sigma1 + sigma3)*sin(fPhi) - T(2.)*A*cos(fPhi);
        fFlag =true;
        return;
    }
    if(TPZExtractVal::val(res0) >=0. && TPZExtractVal::val(res1) < 0. && TPZExtractVal::val(res2) > -1.)
    {
        T Eta = T(6. * sin(fPhi)/(sqrt(3.)*(3.+sin(fPhi))));//INNER
        T Ksi = T(6. * cos(fPhi)/(sqrt(3.) * (3.+sin(fPhi))));//INNER
        T I1,J2,p;
        I1 = sigma.I1();
        J2 = sigma.J2();
        p = I1 / T(3.);
        res[0] =  sqrt( J2 ) + p * T(Eta) - A * T(Ksi);
        fFlag1 =true;
        return;
    }
    
    if(TPZExtractVal::val(res0) >=0. && TPZExtractVal::val(res1) > -1. && TPZExtractVal::val(res2)< 0.)
    {
        T Eta = T(6. * sin(fPhi)/(sqrt(3.)*(3.-sin(fPhi))));//OUTER
        T Ksi = T(6. * cos(fPhi)/(sqrt(3.) * (3.-sin(fPhi))));//OUTER
        T I1,J2,p;
        I1 = sigma.I1();
        J2 = sigma.J2();
        p = I1 / T(3.);
        res[0] =  sqrt( J2 ) + p * T(Eta) - A * T(Ksi);
        fFlag2 =true;
        return;
    }
    
    if(TPZExtractVal::val(res0) > TPZExtractVal::val(res1) && TPZExtractVal::val(res0) > TPZExtractVal::val(res2))
    {
        res[0]=res0;
        fFlag =true;
        return;
    }
    
    if(TPZExtractVal::val(res1) > TPZExtractVal::val(res0) && TPZExtractVal::val(res1) > TPZExtractVal::val(res2))
    {
        T Eta = T(6. * sin(fPhi)/(sqrt(3.)*(3.+sin(fPhi))));//INNER
        T Ksi = T(6. * cos(fPhi)/(sqrt(3.) * (3.+sin(fPhi))));//INNER
        T I1,J2,p;
        I1 = sigma.I1();
        J2 = sigma.J2();
        p = I1 / T(3.);
        res[0] =  sqrt( J2 ) + p * T(Eta) - A * T(Ksi);
        fFlag1 =true;
        return;
    }
    if(TPZExtractVal::val(res2) > TPZExtractVal::val(res0) &&  TPZExtractVal::val(res2) >TPZExtractVal::val(res1))
    {
        T Eta = T(6. * sin(fPhi)/(sqrt(3.)*(3.-sin(fPhi))));//OUTER
        T Ksi = T(6. * cos(fPhi)/(sqrt(3.) * (3.-sin(fPhi))));//OUTER
        T I1,J2,p;
        I1 = sigma.I1();
        J2 = sigma.J2();
        p = I1 / T(3.);
        res[0] =  sqrt( J2 ) + p * T(Eta) - A * T(Ksi);
        fFlag2 =true;
        return;

    }
    
    //res[0]=res0;
    //fFlag =true;
#ifdef MACOS
    feclearexcept(FE_ALL_EXCEPT);
                int Res = fetestexcept(FE_ALL_EXCEPT);
                if(Res)
                {
                                std::cout << " \n " << __PRETTY_FUNCTION__ <<"\n NAN DETECTED \n";
                                DebugStop();
                }
                
                Res = fetestexcept(FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW);
                if(Res)
                {
                                std::cout << " \n " << __PRETTY_FUNCTION__ <<"\n NAN DETECTED \n";
                                DebugStop();
                }
#endif 
}

template <class T>
void TPZYCModifiedMohrCoulomb::N(const TPZTensor<T> & sigma,const T & A,  TPZVec<TPZTensor<T> > & Ndir, int checkForcedYield) const
{
    //Flow vectors feitos no mathematica baseados nas yield functions do chen
    TPZTensor<T> eigenval,dSigma1,dSigma2,dSigma3;
    sigma.Eigenvalue(eigenval, dSigma1, dSigma2, dSigma3);
    //    Ndir.Resize(3);
    
    //    TPZTensor<T> dSigma1Copy(dSigma1);
    //    TPZTensor<T> dSigma2Copy(dSigma2);
    //    TPZTensor<T> dSigma3Copy(dSigma3);
    //    dSigma1.Add(dSigma3,-1);
    //    dSigma1Copy.Add(dSigma3Copy,1);
    //    dSigma1Copy.Multiply(sin(fPhi),1);
    //    dSigma1.Add(dSigma1Copy,1);
    //
    //    Ndir[0]=dSigma1;
    //
    //    sigma.Eigenvalue(eigenval, dSigma1, dSigma2, dSigma3);
    //    dSigma1Copy = dSigma1;
    //    dSigma2Copy = dSigma2;
    //    dSigma3Copy = dSigma3;
    //    dSigma2.Add(dSigma3,-1);
    //    dSigma2Copy.Add(dSigma3Copy,1);
    //    dSigma2Copy.Multiply(sin(fPhi),1);
    //    dSigma2.Add(dSigma2Copy,1);
    //
    //    Ndir[1]=dSigma2;
    //
    //    sigma.Eigenvalue(eigenval, dSigma1, dSigma2, dSigma3);
    //    dSigma1Copy = dSigma1;
    //    dSigma2Copy = dSigma2;
    //    dSigma3Copy = dSigma3;
    //    dSigma1.Add(dSigma2,-1);
    //    dSigma1Copy.Add(dSigma2Copy,1);
    //    dSigma1Copy.Multiply(sin(fPhi),1);
    //    dSigma1.Add(dSigma1Copy,1);
    //
    //    Ndir[2]=dSigma1;
    
    if(fFlag == true)
    {
        //        TPZTensor<T> dSigma1Copy(dSigma1);
        //        TPZTensor<T> dSigma2Copy(dSigma2);
        //        TPZTensor<T> dSigma3Copy(dSigma3);
        //        dSigma1.Add(dSigma3,-1);
        //        dSigma1Copy.Add(dSigma3Copy,1);
        //        dSigma1Copy.Multiply(sin(fPhi),1);
        //        dSigma1.Add(dSigma1Copy,1);
        //        Ndir[0]=dSigma1;
        TPZTensor<T> e1e1,e2e2,e3e3;
        e1e1.XX()=1.;
        e3e3.ZZ()=1.;
        
        e1e1*=T(1.)+T(sin(fPhi));
        e3e3*=T(1.)-T(sin(fPhi));
        e1e1.Add(e3e3,-1.);
        
        Ndir[0]=e1e1;
        
    }
    
    
    if(fFlag1 == true)
    {
        T Eta = T(6. * sin(fPhi)/(sqrt(3.)*(3.+sin(fPhi))));//INNER
        T J2 = sigma.J2();
        TPZTensor<T> s;
        sigma.S(s);
        s *= T(0.5) / sqrt(J2);
        //Hydrostatic part
        T EtaOver3 = Eta/T(3.);
        s.XX() += EtaOver3;
        s.YY() += EtaOver3;
        s.ZZ() += EtaOver3;
        Ndir[0] = s;
    }
    
    if(fFlag2 == true)
    {
        T Eta = T(6. * sin(fPhi)/(sqrt(3.)*(3.-sin(fPhi))));//OUTER
        T J2 = sigma.J2();
        TPZTensor<T> s;
        sigma.S(s);
        s *= T(0.5) / sqrt(J2);
        //Hydrostatic part
        T EtaOver3 = Eta/T(3.);
        s.XX() += EtaOver3;
        s.YY() += EtaOver3;
        s.ZZ() += EtaOver3;
        Ndir[0] = s;
    }
}

template <class T>
void TPZYCModifiedMohrCoulomb::H(const TPZTensor<T> & sigma,const T & A,  TPZVec<T> & h, int checkForcedYield) const
{
    //h.Resize(3);
    if(fFlag == true)
    {
        h[0] = 2.* cos(fPhi);
    }
    if(fFlag1 == true)
    {
        T Ksi = T(6. * cos(fPhi)/(sqrt(3.) * (3.+sin(fPhi))));//INNER
        h[0] = Ksi;
    }
    if(fFlag2 == true)
    {
        T Ksi = T(6. * cos(fPhi)/(sqrt(3.) * (3.-sin(fPhi))));//OUTER
        h[0] = Ksi;
    }
    //h[1] = 2.* cos(fPhi);
    //h[2] = 2.* cos(fPhi);
}

#endif