TPZSandlerExtended.h

Go to the documentation of this file.

//
//  pzsandlerextPV.h
//  PZ
//
//  Created by Diogo Cecilio on 9/3/13.
//
//

#ifndef __PZ__pzsandlerextPV__
#define __PZ__pzsandlerextPV__

#include <iostream>
#include "TPZTensor.h"
#include "TPZElasticResponse.h"
#include "TPZPlasticState.h"
#include "TPZPlasticCriterion.h"

class TPZSandlerExtended : public TPZPlasticCriterion {
public:

    enum {
        NYield = 2
    };

    TPZSandlerExtended(STATE A, STATE B, STATE C, STATE D, STATE K, STATE G, STATE W, STATE R, STATE Phi, STATE N, STATE Psi, STATE kappa_0);
    TPZSandlerExtended(const TPZSandlerExtended & copy);
    TPZSandlerExtended();

    TPZSandlerExtended & operator=(const TPZSandlerExtended & source) {
        ftol = source.ftol;
        fA = source.fA;
        fB = source.fB;
        fC = source.fC;
        fD = source.fD;
        fK = source.fK;
        fG = source.fG;
        fW = source.fW;
        fR = source.fR;
        fPhi = source.fPhi;
        fN = source.fN;
        fPsi = source.fPsi;
        fE = source.fE;
        fnu = source.fnu;
        fkappa_0 = source.fkappa_0;
        fElasticResponse = source.fElasticResponse;

        return *this;
    }

    virtual ~TPZSandlerExtended();

    STATE GetX(STATE k);

    void Firstk(STATE &epsp, STATE &k) const;

    TPZElasticResponse GetElasticResponse();

    void SetElasticResponse(const TPZElasticResponse &ER);

    virtual TPZElasticResponse GetElasticResponse() const;

    STATE GetR();

    virtual void YieldFunction(const TPZVec<STATE> &sigma, STATE kprev, TPZVec<STATE> &yield) const override;

    virtual int GetNYield() const override {
        return as_integer(NYield);
    }

    virtual void Print(std::ostream &out) const override;
    
    template<class T>
    T F(const T x) const;
    
    template<class T>
    T DF(const T x) const;

    STATE GetF(STATE x) const;
    
    virtual int ClassId() const override;

    void Read(TPZStream& buf, void* context) override;

    void Write(TPZStream& buf, int withclassid) const override;
    
private:


    template<class T>
    T X(const T k) const;
    
    template<class T>
    T EpsEqX(T X) const;
    template<class T>
    T EpsEqk(const T k) const;
    template<class T>
    T ResLF2(const TPZVec<T> &pt, T theta, T beta, T k, STATE kprev) const;
    
    template<class T>
    T ResLF2IJ(const TPZVec<T> &sigtrIJ, T theta, T k, STATE kprev) const;
    
    STATE NormalToF1(STATE I1, STATE I1_ref) const;
    
    STATE ResLF1(const TPZVec<STATE> &sigtrial, const TPZVec<STATE> &sigproj, const STATE k, const STATE kprev) const;
    STATE DistF1(const TPZVec<STATE> &pt, const STATE xi, const STATE beta) const;
    STATE DistF2(const TPZVec<STATE> &pt, const STATE theta, const STATE beta, const STATE k) const;

    STATE DistF2IJ(const TPZVec<STATE> &sigtrialIJ, STATE theta, STATE k) const;


    void DDistFunc1(const TPZVec<STATE> &pt, STATE xi, STATE beta, TPZFMatrix<STATE> &ddistf1) const;
    
    template<class T>
    void Res1(const TPZVec<T> &trial_stress, T i1, T beta, T k, T kprev, TPZVec<T> & residue_1) const;
    
    template<class T>
    void Res2(const TPZVec<T> &trial_stress, T theta, T beta, T k, T kprev, TPZVec<T> & residue_2) const;
    
    template<class T>
    void Res2CoVertex(const TPZVec<T> &trial_stress, T beta, T k, T kprev, TPZVec<T> & residue_covertex) const;
    
    template<class T>
    void Res2Vertex(const TPZVec<T> &trial_stress, T k, T kprev, T & residue_vertex) const;
    
    template<class T>
    void DDistF2IJ(TPZVec<T> &sigtrialIJ, T theta, T L, STATE Lprev, TPZVec<T> &ddistf2) const;


    void D2DistFunc1(const TPZVec<STATE> &pt, STATE xi, STATE beta, TPZFMatrix<STATE> &d2distf1) const;

    void Jacobianf1(const TPZVec<STATE> &trial_stress, STATE i1, STATE beta, STATE k, TPZFMatrix<STATE> &jacobianf1)const;
    
    void Jacobianf2(const TPZVec<STATE> &trial_stress, STATE theta, STATE beta, STATE k, TPZFMatrix<STATE> &jacobianf2)const;
    
    void Jacobianf2CoVertex(const TPZVec<STATE> &trial_stress, STATE beta, STATE k, TPZFMatrix<STATE> &jacobianf2_covertex)const;
    
    void Jacobianf2Vertex(const TPZVec<STATE> &trial_stress, STATE k, STATE &jacobianf2_vertex)const;

    void JacobianVertex(const TPZVec<STATE> &trial_stress, STATE k, STATE &jacobian_vertex)const;
    
    void JacobianCoVertex(const TPZVec<STATE> &trial_stress, STATE beta, STATE k, TPZFMatrix<STATE> &jacobian_covertex)const;

    // the derivative are given in terms of theta, beta and k
    void DResLF2(const TPZVec<STATE> &pt, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec<STATE> &dresl) const;
    // the derivative are given in terms of k
    STATE DResLF1(const TPZVec<STATE> &sigtrial, const TPZVec<STATE> &sigproj, const STATE k, const STATE kprev) const;
    
    template<class T>
    void FromThetaKToSigIJ(const T &theta, const T &K, TPZVec<T> &sigIJ) const;

    void DF1Cart(STATE xi, STATE beta, TPZFMatrix<STATE> &DF1) const;

    void DF2Cart(STATE theta, STATE beta, STATE k, TPZFMatrix<STATE> &DF1) const;


    void GradF1SigmaTrial(const TPZVec<STATE> &sigtrial, STATE xi, STATE beta, TPZFMatrix<STATE> &deriv) const;

    void GradF2SigmaTrial(const TPZVec<STATE> &sigtrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZFMatrix<STATE> &deriv) const;

    void F1Cyl(STATE xi, STATE beta, TPZVec<STATE> &f1cyl) const;

    void F2Cyl(const STATE theta, const STATE beta, const STATE k, TPZVec<STATE> &f2cyl) const;

public:
    
    REAL InitialDamage(const TPZVec<REAL> &stress_p) const;

    void Phi(TPZVec<REAL> sigma, STATE alpha, TPZVec<STATE> &phi)const;

    void SurfaceParamF1(TPZVec<STATE> &sigproj, STATE &xi, STATE &beta) const;

    void SurfaceParamF2(const TPZVec<STATE> &sigproj, const STATE k, STATE &theta, STATE &beta) const;

    STATE NormalFunctionToF1(STATE & I1, STATE & k) const;

    void ProjectApex(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const;
    
    void ProjectF1(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const;

    void ProjectF2(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const;
    
    void ProjectCapVertex(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const;
    
    void ProjectCapCoVertex(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const;

    void ProjectCoVertex(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const;
    
    void ProjectVertex(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const;
    
    void ProjectRing(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const;

    void ProjectBetaConstF2(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const;

    //          virtual void ApplyStrainComputeSigma(TPZPlasticState<STATE> &plasticstate, TPZVec<STATE> &sigma);

    //    void ApplyStrainComputeSigma(const TPZVec<STATE> &eps, STATE kprev, TPZVec<STATE> &sigma,STATE &kproj) const;
    void ApplyStrainComputeSigma(TPZVec<STATE> &epst, TPZVec<STATE> &epsp, STATE & kprev, TPZVec<STATE> &epspnext, TPZVec<STATE> &stressnext, STATE & knext, TPZFMatrix<REAL> * tangent = NULL) const;



    void ComputeI1(TPZVec<STATE> stress, STATE &I1)const;
    void ComputeJ2(TPZVec<STATE> stress, STATE &J2)const;
    void ApplyStrainComputeElasticStress(TPZVec<STATE> &strain, TPZVec<STATE> &stress)const;
    void ApplyStressComputeElasticStrain(TPZVec<STATE> &stress, TPZVec<STATE> &strain)const;

    void ProjectSigma(const TPZVec<STATE> &sigmatrial, STATE kprev, TPZVec<STATE> &sigmaproj, STATE &kproj, int &m_type, TPZFMatrix<REAL> * gradient = NULL) const;

    void ProjectSigmaDep(const TPZVec<STATE> &sigmatrial, STATE kprev, TPZVec<STATE> &sigmaproj, STATE &kproj, TPZFMatrix<STATE> &GradSigma) const;

    void TaylorCheckDistF1(const TPZVec<STATE> &sigmatrial, STATE xi, STATE beta, TPZVec<STATE> &xnorm,
            TPZVec<STATE> &errnorm) const;
    void TaylorCheckDDistF1(const TPZVec<STATE> &sigmatrial, STATE xi, STATE beta, TPZVec<STATE> &xnorm,
            TPZVec<STATE> &errnorm) const;
    void TaylorCheckDDistF1DSigtrial(const TPZVec<STATE> &sigmatrial, STATE xi, STATE beta, TPZVec<STATE> &xnorm,
            TPZVec<STATE> &errnorm) const;
    void TaylorCheckDistF2(const TPZVec<STATE> &sigmatrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec<STATE> &xnorm,
            TPZVec<STATE> &errnorm) const;
    void TaylorCheckDDistF2(const TPZVec<STATE> &sigmatrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec<STATE> &xnorm,
            TPZVec<STATE> &errnorm) const;
    void TaylorCheckDDistF2DSigtrial(const TPZVec<STATE> &sigmatrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec<STATE> &xnorm,
            TPZVec<STATE> &errnorm) const;
    void TaylorCheckDF1Cart(STATE xi, STATE beta, TPZVec<STATE> &xnorm,
            TPZVec<STATE> &errnorm) const;

    void TaylorCheckDF2Cart(STATE theta, STATE beta, STATE k, TPZVec<STATE> &xnorm,
            TPZVec<STATE> &errnorm) const;

    void TaylorCheckParamF1Sigtrial(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const;

    void TaylorCheckProjectF1(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const;

    void TaylorCheckDtbkDsigtrial(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const;

    void TaylorCheckProjectF2(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const;

    static void ConvergenceRate(TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm, TPZVec<STATE> &convergence);

    static void CheckCoordinateTransformation(TPZVec<STATE> &cart);
    
    void ComputeCapTangent(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj, TPZFMatrix<REAL> * gradient) const;
    
    void ComputeCapVertexTangent(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj, TPZFMatrix<REAL> * gradient) const;

    void ComputeCapCoVertexTangent(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj, TPZFMatrix<REAL> * gradient) const;
    
    void ComputeFailureTangent(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj, TPZFMatrix<REAL> * gradient) const;
    
    void TaylorCheckProjectSigma(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const;
    

public:
    
    // Set up charecteristic parameters
    void SetUp(STATE A, STATE B, STATE C, STATE D, STATE K, STATE G, STATE W, STATE R, STATE Phi, STATE N, STATE Psi);
    
    // Set up the initial damage
    void SetInitialDamage(STATE kappa_0);
    
    static void MCormicRanchSand(TPZSandlerExtended &mat);
    static void ReservoirSandstone(TPZSandlerExtended &mat);
    static void SalemLimestone(TPZSandlerExtended &mat);
    static void PreSMat(TPZSandlerExtended &mat); // em MPa

    STATE A() {
        return fA;
    }
    
    void SetA(STATE A)
    {
        fA = A;
    }

    STATE B() {
        return fB;
    }
    
    void SetB(STATE B)
    {
        fB = B;
    }

    STATE C() {
        return fC;
    }
    
    void SetC(STATE C)
    {
        fC = C;
    }

    STATE D() {
        return fD;
    }

    STATE W() {
        return fW;
    }

    STATE K() {
        return fK;
    }

    STATE R() {
        return fR;
    }

    STATE G() {
        return fG;
    }

    STATE E() {
        return fE;
    }

    STATE N() {
        return fN;
    }

    STATE Poisson() {
        return fnu;
    }
    
    STATE InitialDamage() {
        return fkappa_0;
    }
    
    STATE Apex() const {
        STATE apex = log(fA/fC)/fB;
        return apex;
    }
    
    STATE X_0() const {
        STATE X_0 = this->X(fkappa_0);
        return X_0;
    }
    
    STATE CPerturbation() const {
        STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
        STATE C_per = (fD*fC)/(3.0*CK);
        return C_per;
    }
    
    STATE ftol;

private:

    STATE fA, fB, fC, fD, fW, fK, fR, fG, fPhi, fN, fPsi, fE, fnu, fkappa_0;

    //    bool fIsonCap;
    TPZElasticResponse fElasticResponse;


};


#endif /* defined(__PZ__pzsandlerextPV__) */