concretedpm2.h

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2013   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */

#ifndef ConcreteDPM2_h
#define ConcreteDPM2_h

#include "../sm/Materials/structuralmaterial.h"
#include "floatarray.h"
#include "floatmatrix.h"
#include "cltypes.h"
#include "../sm/Materials/structuralms.h"
#include "Materials/isolinearelasticmaterial.h"
#include "gausspoint.h"
#include "mathfem.h"

#define DYNCON_TOL 1.e-6
#define keep_track_of_dissipated_energy

#define _IFT_ConcreteDPM2_Name "con2dpm"
#define _IFT_ConcreteDPM2_fc "fc"
#define _IFT_ConcreteDPM2_fcZero "fczero"
#define _IFT_ConcreteDPM2_ft "ft"
#define _IFT_ConcreteDPM2_ecc "ecc"
#define _IFT_ConcreteDPM2_kinit "kinit"
#define _IFT_ConcreteDPM2_ahard "ahard"
#define _IFT_ConcreteDPM2_bhard "bhard"
#define _IFT_ConcreteDPM2_chard "chard"
#define _IFT_ConcreteDPM2_dhard "dhard"
#define _IFT_ConcreteDPM2_dilation "dilation"
#define _IFT_ConcreteDPM2_asoft "asoft"
#define _IFT_ConcreteDPM2_hp "hp"
#define _IFT_ConcreteDPM2_yieldtol "yieldtol"
#define _IFT_ConcreteDPM2_newtoniter "newtoniter"
#define _IFT_ConcreteDPM2_wf "wf"
#define _IFT_ConcreteDPM2_efc "efc"
#define _IFT_ConcreteDPM2_softeningType "stype"
#define _IFT_ConcreteDPM2_ftOne "ft1"
#define _IFT_ConcreteDPM2_wfOne "wf1"
#define _IFT_ConcreteDPM2_rateFlag "rateflag"
#define _IFT_ConcreteDPM2_timeFactor "timefactor"
#define _IFT_ConcreteDPM2_helem "helem"
#define _IFT_ConcreteDPM2_isoflag "isoflag"


namespace oofem {
class ConcreteDPM2Status : public StructuralMaterialStatus
{
public:
    enum state_flag_values {
        ConcreteDPM2_Elastic,
        ConcreteDPM2_Unloading,
        ConcreteDPM2_Plastic,
        ConcreteDPM2_Damage,
        ConcreteDPM2_PlasticDamage,
        ConcreteDPM2_VertexCompression,
        ConcreteDPM2_VertexTension,
        ConcreteDPM2_VertexCompressionDamage,
        ConcreteDPM2_VertexTensionDamage
    };
protected:

    FloatArray plasticStrain;
    FloatArray tempPlasticStrain;

    double dFDKappa;
    double deltaLambda;


    double kappaP;
    double tempKappaP;

    double kappaPPeak;

    double le;

    double alpha;
    double tempAlpha;

    double equivStrain;
    double tempEquivStrain;

    double equivStrainTension;
    double tempEquivStrainTension;

    double equivStrainCompression;
    double tempEquivStrainCompression;

    double kappaDTension;
    double tempKappaDTension;

    double kappaDCompression;
    double tempKappaDCompression;

    double kappaDTensionOne;
    double tempKappaDTensionOne;

    double kappaDCompressionOne;
    double tempKappaDCompressionOne;

    double kappaDTensionTwo;
    double tempKappaDTensionTwo;

    double kappaDCompressionTwo;
    double tempKappaDCompressionTwo;

    double damageTension;
    double tempDamageTension;

    double damageCompression;
    double tempDamageCompression;

    double deltaEquivStrain;

    double rateFactor;
    double tempRateFactor;

    double rateStrain;
    double tempRateStrain;

    int state_flag;
    int temp_state_flag;
#ifdef keep_track_of_dissipated_energy
    double stressWork;
    double tempStressWork;
    double dissWork;
    double tempDissWork;
#endif

public:
    ConcreteDPM2Status(int n, Domain *d, GaussPoint *gp);

    virtual ~ConcreteDPM2Status();

    virtual void initTempStatus();
    virtual void updateYourself(TimeStep *tStep);
    virtual void printOutputAt(FILE *file, TimeStep *tStep);
    virtual contextIOResultType saveContext(DataStream &, ContextMode mode, void *obj = NULL);
    virtual contextIOResultType restoreContext(DataStream &, ContextMode mode, void *obj = NULL);
    virtual const char *giveClassName() const { return "ConcreteDPM2Status"; }

    // Inline functions for access to state variables

    const FloatArray &givePlasticStrain() const { return plasticStrain; }


    double giveDeviatoricPlasticStrainNorm()
    {
        FloatArray dev;
        StructuralMaterial :: computeDeviatoricVolumetricSplit(dev, plasticStrain);
        return sqrt( .5 * ( 2. * dev [ 0 ] * dev [ 0 ] + 2. * dev [ 1 ] * dev [ 1 ] + 2. * dev [ 2 ] * dev [ 2 ] +
                    dev [ 3 ] * dev [ 3 ] + dev [ 4 ] * dev [ 4 ] + dev [ 5 ] * dev [ 5 ] ) );
    }

    double giveVolumetricPlasticStrain() const
    {
        return 1. / 3. * ( plasticStrain(0) + plasticStrain(1) + plasticStrain(2) );
    }

    double giveKappaP() const
    { return kappaP; }

    double giveKappaDTensionOne() const
    { return kappaDTensionOne; }

    double giveKappaDCompressionOne() const
    { return kappaDCompressionOne; }


    double giveKappaDTensionTwo() const
    { return kappaDTensionTwo; }


    double giveKappaDCompressionTwo() const
    { return kappaDCompressionTwo; }


    double giveEquivStrain() const
    { return equivStrain; }

    double giveEquivStrainTension() const
    { return equivStrainTension; }


    double giveEquivStrainCompression() const
    { return equivStrainCompression; }

    double giveDamageTension() const
    { return damageTension; }


    double giveDamageCompression() const
    { return damageCompression; }



    double giveRateFactor() const
    { return rateFactor; }

    double giveTempRateFactor() const
    { return tempRateFactor; }


    double giveRateStrain() const
    { return rateStrain; }

    void letTempRateStrainBe(double v)
    { tempRateStrain = v; }

    void letTempAlphaBe(double v)
    { tempAlpha = v; }

    int giveStateFlag() const
    { return state_flag; }


    // giveTemp:

    // Functions used to access the temp variables.
    const FloatArray &giveTempPlasticStrain() const { return tempPlasticStrain; }

    double giveTempVolumetricPlasticStrain() const
    { return 1. / 3. * ( tempPlasticStrain(0) + tempPlasticStrain(1) + tempPlasticStrain(2) ); }


    double giveTempKappaP() const
    { return tempKappaP; }


    double giveKappaDTension() const
    { return kappaDTension; }

    double giveAlpha() const
    { return alpha; }


    double giveKappaDCompression() const
    { return kappaDCompression; }


    double giveTempDamageTension() const
    { return tempDamageTension; }


    double giveTempDamageCompression() const
    { return tempDamageCompression; }

    double giveDeltaEquivStrain() const
    { return deltaEquivStrain; }


    int giveTempStateFlag() const
    { return temp_state_flag; }


    // letTemp...be :
    // Functions used by the material to assign a new value to a temp variable.
    void letTempPlasticStrainBe(const FloatArray &v)
    { tempPlasticStrain = v; }

    void letDeltaLambdaBe(double v)
    { deltaLambda = v; }

    void letTempKappaPBe(double v)
    { tempKappaP = v; }

    void letTempKappaDTensionBe(double v)
    { tempKappaDTension = v; }

    void letTempKappaDCompressionBe(double v)
    { tempKappaDCompression = v; }

    void letTempKappaDTensionOneBe(double v)
    { tempKappaDTensionOne = v; }

    void letTempKappaDCompressionOneBe(double v)
    { tempKappaDCompressionOne = v; }

    void letTempKappaDTensionTwoBe(double v)
    { tempKappaDTensionTwo = v; }

    void letTempKappaDCompressionTwoBe(double v)
    { tempKappaDCompressionTwo = v; }

    void letTempDamageTensionBe(double v)
    { tempDamageTension = v; }

    void letTempDamageCompressionBe(double v)
    { tempDamageCompression = v; }

    void letTempRateFactorBe(double v)
    { tempRateFactor = v; }

    void letTempEquivStrainBe(double v)
    { tempEquivStrain = v; }

    void letTempEquivStrainTensionBe(double v)
    { tempEquivStrainTension = v; }

    void letTempEquivStrainCompressionBe(double v)
    { tempEquivStrainCompression = v; }


    double giveLe() { return le; }

    void setLe(double ls)
    { le = ls; }
    void letTempStateFlagBe(const int v)
    { temp_state_flag = v; }


    void letKappaPPeakBe(double kappa)
    { kappaPPeak = kappa; }
#ifdef keep_track_of_dissipated_energy
    double giveStressWork() { return stressWork; }
    double giveTempStressWork() { return tempStressWork; }
    void setTempStressWork(double w) { tempStressWork = w; }
    double giveDissWork() { return dissWork; }
    double giveTempDissWork() { return tempDissWork; }
    void setTempDissWork(double w) { tempDissWork = w; }
    void computeWork(GaussPoint *gp, double ft);
#endif
};


//   ********************************
//   *** CLASS DYNAMIC CONCRETE   ***
//   ********************************

class ConcreteDPM2 : public StructuralMaterial
{
public:

protected:
    enum ConcreteDPM2_ReturnType { RT_Regular, RT_Tension, RT_Compression, RT_Auxiliary };
    ConcreteDPM2_ReturnType returnType;

    enum ConcreteDPM2_ReturnResult { RR_NotConverged, RR_Converged };
    ConcreteDPM2_ReturnResult returnResult;

    double fc, ft, ecc;

    int isotropicFlag;

    double e0;

    double AHard;
    double BHard;
    double CHard;
    double DHard;

    double hardeningModulus;

    double ASoft;

    double yieldHardPrimePeak;

    double yieldHardInitial;

    double dilationConst;

#if 1
    double sig;

    double rho;

    double thetaTrial;
#endif

    double m;

    double mQ;

    double helem;

    LinearElasticMaterial *linearElasticMaterial;

    double eM;
    double gM;
    double kM;
    double nu;

    double efCompression;

    double wf;

    double wfOne;

    double ftOne;

    double yieldTol;

    double yieldTolDamage;

    int newtonIter;

    int softeningType;

    double timeFactor;

    double fcZero;

    int strainRateFlag;


public:
    ConcreteDPM2(int n, Domain *d);
    virtual ~ConcreteDPM2();
    virtual IRResultType initializeFrom(InputRecord *ir);

    virtual const char *giveClassName() const { return "ConcreteDPM2"; }
    virtual const char *giveInputRecordName() const { return _IFT_ConcreteDPM2_Name; }

    LinearElasticMaterial *giveLinearElasticMaterial()
    { return linearElasticMaterial; }

    virtual void giveRealStressVector_1d(FloatArray &answer, GaussPoint *gp, const FloatArray &totalStrain, TimeStep *tStep);

    virtual void giveRealStressVector_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &strainVector, TimeStep *tStep);

    void performPlasticityReturn(GaussPoint *gp,
                                 const FloatMatrix &D,
                                 const FloatArray &strain);

    bool checkForVertexCase(double &answer,
                            double sig,
                            double tempKappa,
                            bool mode1d);

    double performRegularReturn(FloatArray &stress,
                                double kappaP,
                                GaussPoint *gp);

    void compute1dJacobian(FloatMatrix &answer,
                           double totalsigma,
                           double tempKappa,
                           double deltaLambda,
                           GaussPoint *gp);
    void computeJacobian(FloatMatrix &answer,
                         double sig,
                         double rho,
                         double tempKappa,
                         double deltaLambda,
                         GaussPoint *gp);

    double performVertexReturn(FloatArray &stress,
                               double apexStress,
                               double tempKappaP,
                               GaussPoint *gp);

    double computeYieldValue(double sig,
                             double rho,
                             double theta,
                             double tempKappa) const;

    double computeHardeningOne(double tempKappa) const;

    double computeHardeningOnePrime(double tempKappa) const;

    double computeHardeningTwo(double tempKappa) const;

    double computeHardeningTwoPrime(double tempKappa) const;

    double computeDFDKappa(double sig,
                           double rho,
                           double tempKappa,
                           bool mode1d);


    double computeDKappaDDeltaLambda(double sig, double rho, double tempKappa);
    double computeDKappaDDeltaLambda1d(double sig, double tempKappa);


    virtual double computeDuctilityMeasure(double sig,
                                           double rho,
                                           double theta);



    void computeDDuctilityMeasureDInv(FloatArray &answer,
                                      double sig,
                                      double rho,
                                      double tempKappa);
    double computeDDuctilityMeasureDInv1d(double sigma, double tempKappa); //Dimitris change 1d implementation

    void computeDGDInv(FloatArray &answer,
                       double sig,
                       double rho,
                       double tempKappa);
    double computeDGDInv1d(double sig, double tempKappa);

    double computeRatioPotential(double sig,
                                 double tempKappa);

    double computeRateFactor(double alpha,
                             double timeFactor,
                             GaussPoint *gp,
                             TimeStep *deltaTime);

    void computeDDGDDInv(FloatMatrix &answer,
                         double sig,
                         double rho,
                         double tempKappa);
    double computeDDGDDInv1d(double sigma, double tempKappa);
    void computeDDGDInvDKappa(FloatArray &answer,
                              double sig,
                              double rho,
                              double tempKappa);

    double computeDDGDInvDKappa1d(double sigma, double tempKappa);
    void computeDDKappaDDeltaLambdaDInv(FloatArray &answer,
                                        double sig,
                                        double rho,
                                        double tempKappa);

    double computeDDKappaDDeltaLambdaDInv1d(double sigma, double tempKappa);
    double computeDDKappaDDeltaLambdaDKappa(double sig, double rho, double tempKappa);
    double computeDDKappaDDeltaLambdaDKappa1d(double sig, double tempKappa);


    void computeDFDInv(FloatArray &answer,
                       double sig,
                       double rho,
                       double tempKappa) const;
    double computeDFDInv1d(double sigma, double tempKappa) const;
    double computeTempKappa(double kappaInitial,
                            double sigTrial,
                            double rhoTrial,
                            double sig);


    void  computeDamage(FloatArray &answer, const FloatArray &strain, const FloatMatrix &D, double timeFactor, GaussPoint *gp, TimeStep *tStep, double alpha);


    int checkForUnAndReloading(double &tempEquivStrain,
                               double &minEquivStrain,
                               const FloatMatrix &D,
                               GaussPoint *gp);

    double computeAlpha(FloatArray &effectiveStressTension, FloatArray &effectiveStressCompression, FloatArray &effectiveStress);

    virtual double computeDamageParamTension(double equivStrain, double kappaOne, double kappaTwo, double le, double omegaOld);

    virtual double computeDamageParamCompression(double equivStrain, double kappaOne, double kappaTwo, double omegaOld);

    double computeDeltaPlasticStrainNormTension(double tempKappaD, double kappaD, GaussPoint *gp);

    double computeDeltaPlasticStrainNormCompression(double tempAlpha, double tempKappaD, double kappaD, GaussPoint *gp);

    virtual double computeEquivalentStrain(double sig, double rho, double theta);


    double computeDuctilityMeasureDamage(const FloatArray &strain, GaussPoint *gp);

    void initDamaged(double kappa,
                     const FloatArray &strain,
                     GaussPoint *gp);


    void computeTrialCoordinates(const FloatArray &stress, double &sig, double &rho, double &theta);


    void assignStateFlag(GaussPoint *gp);

    void computeDRhoDStress(FloatArray &answer,
                            const FloatArray &stress) const;

    void computeDSigDStress(FloatArray &answer) const;

    void computeDDRhoDDStress(FloatMatrix &answer,
                              const FloatArray &stress) const;

    void computeDCosThetaDStress(FloatArray &answer,
                                 const FloatArray &stress) const;

    double computeDRDCosTheta(double theta, double ecc) const;

    virtual void give1dStressStiffMtrx(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep);

    virtual void give3dMaterialStiffnessMatrix(FloatMatrix &answer,
                                               MatResponseMode mode, GaussPoint *gp, TimeStep *tStep);

    void compute3dSecantStiffness(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep);



    virtual bool isCharacteristicMtrxSymmetric(MatResponseMode rMode) { return false; }

    virtual int giveIPValue(FloatArray &answer,
                            GaussPoint *gp,
                            InternalStateType type,
                            TimeStep *tStep);

protected:
    MaterialStatus *CreateStatus(GaussPoint *gp) const;
};
} //end namespace oofem
#endif