TPZPlasticStep.h

Go to the documentation of this file.

#ifndef TPZPLASTICSTEP_H
#define TPZPLASTICSTEP_H

#include "TPZTensor.h"
#include "pzreal.h"
#include "pzfmatrix.h"
#include "TPZPlasticState.h"
#include "TPZPlasticIntegrMem.h"
#include "TPZElasticResponse.h"

#include "pzlog.h"
#include "TPZPlasticBase.h"

#include <set>
#include <ostream>

enum EElastoPlastic
{
                EAuto = 0,
                EForceElastic = 1,
                EForcePlastic = 2
};

template <class YC_t, class TF_t, class ER_t>
class TPZPlasticStep: public TPZPlasticBase 
{
public:
                
virtual int ClassId() const override;
    TPZPlasticStep(REAL alpha=0.):fYC(), fTFA(), fER(), fResTol(1.e-12),
                fIntegrTol(1.e-4), fMaxNewton(30), fMinLambda(1.),
                fMinStepSize(1.e-3), fN(), fPlasticMem(),
                fMaterialTensionSign(1), fInterfaceTensionSign(1)
                { fN.m_hardening = alpha; }
                
    TPZPlasticStep(const TPZPlasticStep & source)
                {
                                fYC                                             = source.fYC;
                                fTFA                            = source.fTFA;
                                fER                                             = source.fER;
                                fResTol                         = source.fResTol;
                                fIntegrTol  = source.fIntegrTol;
                                fMaxNewton  = source.fMaxNewton;                
                                fMinLambda  = source.fMinLambda;
                                fMinStepSize= source.fMinStepSize;
                                fN                                              = source.fN;
                                fPlasticMem = source.fPlasticMem;
                                fMaterialTensionSign  = source.fMaterialTensionSign;
                                fInterfaceTensionSign = source.fInterfaceTensionSign;
                }
                
                TPZPlasticStep & operator=(const TPZPlasticStep & source)
                {
                                fYC                                             = source.fYC;
                                fTFA                            = source.fTFA;
                                fER                                             = source.fER;
                                fResTol                         = source.fResTol;
                                fIntegrTol  = source.fIntegrTol;
                                fMaxNewton  = source.fMaxNewton;                
                                fMinLambda  = source.fMinLambda;
                                fMinStepSize= source.fMinStepSize;
                                fN                                              = source.fN;
                                fPlasticMem = source.fPlasticMem;
                                fMaterialTensionSign  = source.fMaterialTensionSign;
                                fInterfaceTensionSign = source.fInterfaceTensionSign;
                                return *this;
                }
                
                virtual const char *Name() const override
                {
                                return "TPZPlasticStep generic class";          
                }
                
                virtual void Print(std::ostream & out) const override
                {
                                out << "\n" << this->Name();
                                out << "\n YC_t:";
                                fYC.Print(out);
                                out << "\n TF_t:";
                                fTFA.Print(out);                
                                out << "\n ER_t:";
                                fER.Print(out);
                                out << "\nTPZPlasticStep Internal members:";
                                out << "\n fResTol = " << fResTol;
                                out << "\n fMaxNewton = " << fMaxNewton;
                                out << "\n fMinLambda = " << fMinLambda;
                                out << "\n fPlasticMem.NElements() = " << fPlasticMem.NElements();
                                out << "\n fN = ";
                                fN.Print(out);
                }
                
                virtual void ApplyStrain(const TPZTensor<REAL> &epsTotal) override;
    
    
    void SetOutFile(string outfile);
                
                typedef YC_t fNYields;
                
                virtual void ApplyStrainComputeSigma(const TPZTensor<REAL> &epsTotal, TPZTensor<REAL> &sigma, TPZFMatrix<REAL> * tangent = NULL) override;
                
                virtual void ApplyStrainComputeDep(const TPZTensor<REAL> &epsTotal, TPZTensor<REAL> &sigma, TPZFMatrix<REAL> &Dep) override;
                
    virtual void ApplyLoad(const TPZTensor<REAL> & sigma, TPZTensor<REAL> &epsTotal) override;
                
    void ResetPlasticMem()
    {
        fPlasticMem.Resize(0);
    }
    
                void SetTensionSign(int s);
                
                int SignCorrection()const;
        
        TPZPlasticCriterion& GetYC() override{
            return fYC;
        }

                
protected:
                
    void ApplyStrain_Internal(const TPZTensor<REAL> &epsTotal);
                
    void ApplyStrainComputeDep_Internal(const TPZTensor<REAL> &epsTotal, TPZTensor<REAL> &sigma, TPZFMatrix<REAL> &Dep);
                
                void ApplyStrainComputeSigma_Internal(const TPZTensor<REAL> &epsTotal, TPZTensor<REAL> &sigma)
                {
                                ApplyStrain_Internal(epsTotal);
                                
                                REAL A;
                                
                                int n = fPlasticMem.NElements();
                                
                                ComputePlasticVars(fPlasticMem[n-1].m_elastoplastic_state, sigma, A);
                }
                
    void ApplyLoad_Internal(const TPZTensor<REAL> & sigma, TPZTensor<REAL> &epsTotal);
                
                
public:
    virtual void Phi(const TPZTensor<REAL> &epsTotal, TPZVec<REAL> &phi) const override;
                
                
public:
                
    virtual void Phi_Internal(const TPZTensor<REAL> &epsTotal, TPZVec<REAL> &phi) const;
                
    bool IsStrainElastic(const TPZPlasticState<REAL> &state)const;
                
    virtual void SetElasticResponse(TPZElasticResponse &ER) override;

    virtual TPZElasticResponse GetElasticResponse() const override;
    virtual void SetState(const TPZPlasticState<REAL> &state) override;
    
    virtual void SetUp(const TPZTensor<REAL> & epsTotal);
                
    virtual TPZPlasticState<REAL> GetState()const override;
                
                virtual int IntegrationSteps()const override;
                
                virtual void SetIntegrTol(REAL integrTol)
                {
                                fIntegrTol = integrTol;
                }
                
                virtual void SetMinStepSize(REAL minStepSize)
                {
                                fMinStepSize = minStepSize;
                }
                
                virtual REAL IntegrationOverview(TPZVec<REAL> & plastifLen);
                
protected:      
                
    virtual void SetState_Internal(const TPZPlasticState<REAL> &state);
                
public:
    virtual TPZPlasticState<REAL> GetState_Internal()const;
                
protected:

    virtual void ProcessStrain(const TPZTensor<REAL> &epsTotal, const EElastoPlastic ep = EAuto);

    virtual void ProcessStrainNoSubIncrement(const TPZTensor<REAL> &epsTotal, const EElastoPlastic ep = EAuto);
                
    virtual void ProcessLoad(const TPZTensor<REAL> &sigma, const EElastoPlastic ep = EAuto);

    virtual void ComputeDep(TPZTensor<REAL> & sigma, TPZFMatrix<REAL> &Dep);
//    virtual void ComputeDep2(TPZTensor<REAL> & sigma, TPZFMatrix<REAL> &Dep);
                
                template<class T>
    void ComputePlasticVars(const TPZPlasticState<T> & state_T,
                                                                                                                TPZTensor<T> & sigma_T,
                                                                                                                T& A_T)const;
                
                REAL FindPointAtYield(const TPZTensor<REAL>& epsTotalNp1,
                                                                                                  TPZPlasticState<REAL> & stateAtYield)const ;
                
    template <class T1, class T2>
    void PlasticResidual(const TPZPlasticState<T1> &N_T1,
                                                                                                 TPZPlasticState<T2> &Np1_T2,
                                                                                                 const TPZVec<T2> &delGamma_T2,
                                                                                                       TPZVec<T2> &res_T2,
                               REAL &normEpsPErr,
                               int silent = 0)const; 
    

    template <class T1, class T2>
    void PlasticResidualRK(
                                                             const TPZPlasticState<T1> &N_T1,
                                                             TPZPlasticState<T2> &Np1_T2,
                                                             const TPZVec<T2> &delGamma_T2,
                                                             TPZVec<T2> &res_T2,
                                                             REAL &normEpsPErr,
                                                             int silent=0)const;
                                
    template<class T1, class T2, class TVECTOR>
    REAL UpdatePlasticVars(const TPZPlasticState<T1> &N_T1,
                                                                                                   TPZPlasticState<T2> &Np1_T2,
                                                                                                   TPZVec<T2> &delGamma_T2,
                                                                                                   TPZVec<T2> &res_T2,
                                                                                                   TVECTOR & Sol_TVECTOR,
                                                                                                   TPZVec<int> & validEqs,
                                                                                                   int updateVars = 1)const;   
                
                template<class T>
                void InitializePlasticFAD(const TPZPlasticState<REAL> & state, const TPZVec<REAL> &delGamma,
                                                                                                                  TPZPlasticState<T> &state_T,         TPZVec<T> &delGamma_T,
                                                                                                                  const int nVars = 7 + YC_t::NYield)const;
                
                template <class T>
                int InitializeValidEqs(TPZVec<T> &res_T, TPZVec<int> & validEqs);
                
    template <class T>
    int RemoveInvalidEqs(TPZVec<T> & delGamma_T, TPZVec<T> & res_T, TPZVec<int> &validEqs);
                
    void InitialGuess(
                      const TPZPlasticState<REAL> &N,
                      TPZPlasticState<REAL> &Np1,
                      TPZVec<REAL> &delGamma,
                      TPZVec<int> &validEqs
                      );
    int PlasticLoop(const TPZPlasticState<REAL> &N,
                                                                                TPZPlasticState<REAL> &Np1,
                                                                                TPZVec<REAL> &delGamma,
                                                                                REAL &normEpsPErr,
                                                                                REAL &lambda,
                                                                                TPZVec<int> & validEqs);
                
    int PlasticIntegrate(
                                                                                                 const TPZPlasticState<REAL> &N,
                                                                                                 TPZPlasticState<REAL> &Np1,
                                                                                                 const REAL &TolEpsPErr);
                
                
    template<class T1, class T_VECTOR, class T_MATRIX>//T1:input residual fad type (FAD FAD or FAD), T_MATRIX: output matrix &vector of type (FAD or REAL, respectvly)
    void ExtractTangent(
                                                                                                const TPZVec<T1> & epsRes_FAD,
                                                                                                T_VECTOR & ResVal, //TPZFMatrix<REAL> for the T1=fad<real> type
                                                                                                REAL & resnorm, // REAL 
                                                                                                T_MATRIX & tangent, // TPZFMatrix<REAL> for T1=fad<real> type
                        TPZVec<int> & validEqs,
                        TPZVec<REAL> &pivots,
                                                                                                const int precond = 1,
                                                                                                const int resetInvalidEqs = 1);
                
                template <int N>
                void PushPlasticMem(const TPZPlasticState<REAL> & state,
                                                                                                const REAL & k,
                                                                                                const REAL & lambda,
                                                                                                const TPZManVector<REAL,N> & delGamma,
                                                                                                const TPZVec<int> & validEqs,
                                                                                                const int forceYield);
                
public:
                void SetMaterialElasticOrPlastic(int choice=1/*plastic*/);
    
    void SetResidualTolerance(STATE tol)
    {
        fResTol = tol;
    }
                
    void Write(TPZStream& buf, int withclassid) const override;
    void Read(TPZStream& buf, void* context) override;
public:
                
    YC_t fYC;
    TF_t fTFA;
    ER_t fER;
                
                
protected:
                
                REAL fResTol;
                
public:
                REAL fIntegrTol;
                
                int fMaxNewton; 
                
                REAL fMinLambda;
                
                REAL fMinStepSize;
                
protected:
                
    TPZPlasticState<REAL> fN;
                
                TPZStack< TPZPlasticIntegrMem<REAL, YC_t::NYield> > fPlasticMem;
                
                int fMaterialTensionSign;
                
                int fInterfaceTensionSign;
                
    
    
    //ofstream fOutfile(string &str);
                
public:
                
                
    int NumCases() 
    {
                                return 1;
    }
                
    static TPZTensor<REAL> gRefDeform;
                
    void LoadState(TPZFMatrix<REAL> &state)
    {
                                int i;
                                for(i=0; i<6; i++) gRefDeform[i] = state(i,0);
    }
                
    void ComputeTangent(TPZFMatrix<REAL> &tangent, TPZVec<REAL> &coefs, int icase)
    {
#ifdef LOG4CXX_PLASTICITY
                                LoggerPtr logger(Logger::getLogger("plasticity.plasticstep"));
#endif
                                TPZTensor<REAL> sigma;
                                switch(icase)
                                {
                                                case 0:
                                                                TPZTensor<REAL> sig;
                                                                tangent.Redim(6,6);
                                                {
                                                                ProcessStrain(gRefDeform);
                                                                ComputeDep(sig, tangent);
                                                }
                                                                //Sigma(gRefDeform,sig,tangent);
#ifdef LOG4CXX_PLASTICITY
                                                                std::stringstream sout;
                                                                sout << "matriz tangent for checkconv " << tangent;
                                                                LOGPZ_DEBUG(logger,sout.str().c_str());
#endif
                                                                break;
                                }
    }
                
    void Residual(TPZFMatrix<REAL> &res,int icase)
    {
#ifdef LOG4CXX_PLASTICITY
                                LoggerPtr logger(Logger::getLogger("plasticity.plasticstep"));
#endif
                                TPZTensor<REAL> sigma;
                                switch(icase)
                                {
                                                case 0:
                                                                TPZTensor<REAL> sig;
                                                                TPZFMatrix<REAL> tangent(6,6);
                                                                res.Redim(6,1);
                                                {
                                                                ProcessStrain(gRefDeform);
                                                                ComputeDep(sig, tangent);
                                                }
                                                                //Sigma(gRefDeform,sig,tangent);
#ifdef LOG4CXX_PLASTICITY
                                                                std::stringstream sout;
                                                                sout << "sigma for residual " << sig;
                                                                LOGPZ_DEBUG(logger,sout.str().c_str());
#endif
                                                                int i;
                                                                for(i=0; i<6; i++)
                                                                {
                                                                                res(i,0) = sig[i];
                                                                }
                                                                break;
                                }
    }
                

};

template <class YC_t, class TF_t, class ER_t>
int TPZPlasticStep<YC_t, TF_t, ER_t>::ClassId() const{
    return Hash("TPZPlasticStep") ^ TPZPlasticBase::ClassId() << 1 ^ YC_t().ClassId() << 2 ^  TF_t().ClassId() << 3 ^ ER_t().ClassId() << 4;
}

#endif //TPZPLASTICSTEP_H