concretedpm2.h

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2025   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 "sm/Materials/isolinearelasticmaterial.h"
#include "gausspoint.h"
#include "mathfem.h"
 
#define CDPM2_TOL 1.e-6
#define keep_track_of_dissipated_energy

#define _IFT_ConcreteDPM2_Name "con2dpm"
#define _IFT_ConcreteDPM2_fc "fc"
#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_strengthratetype "sratetype"
#define _IFT_ConcreteDPM2_energyratetype "eratetype"
#define _IFT_ConcreteDPM2_deltatime "deltat"
#define _IFT_ConcreteDPM2_helem "helem"
#define _IFT_ConcreteDPM2_damflag "damflag"
 
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:


    FloatArrayF< 6 >plasticStrain;
    FloatArrayF< 6 >tempPlasticStrain;
 
    FloatArrayF< 6 >reducedStrain;
    FloatArrayF< 6 >tempReducedStrain;
 
    FloatArrayF< 6 >effectiveStress;
    FloatArrayF< 6 >tempEffectiveStress;


    double kappaP = 0.;
    double tempKappaP = 0.;
 
    double kappaPPeak = 0.;
 
    double deltaLambda = 0.;
 
    double le = 0.;
 
    double alpha = 0.;
    double tempAlpha = 0.;
 
    double equivStrain = 0.;
    double tempEquivStrain = 0.;
 
    double equivStrainTension = 0.;
    double tempEquivStrainTension = 0.;
 
    double equivStrainCompression = 0.;
    double tempEquivStrainCompression = 0.;
 
    double kappaDTension = 0.;
    double tempKappaDTension = 0.;
 
    double kappaDCompression = 0.;
    double tempKappaDCompression = 0.;
 
    double kappaDTensionOne = 0.;
    double tempKappaDTensionOne = 0.;
 
    double kappaDCompressionOne = 0.;
    double tempKappaDCompressionOne = 0.;
 
    double kappaDTensionTwo = 0.;
    double tempKappaDTensionTwo = 0.;
 
    double kappaDCompressionTwo = 0.;
    double tempKappaDCompressionTwo = 0.;
 
    double damageTension = 0.;
    double tempDamageTension = 0.;
 
    double damageCompression = 0.;
    double tempDamageCompression = 0.;
 
    double deltaEquivStrain = 0.;
 
    double rateFactor = 1.;
    double tempRateFactor = 0.;

    double rateStrain = 0.;
    double tempRateStrain = 0.;

    int state_flag = ConcreteDPM2Status::ConcreteDPM2_Elastic;
    int temp_state_flag = ConcreteDPM2Status::ConcreteDPM2_Elastic;
 
 
#ifdef keep_track_of_dissipated_energy
    double stressWork = 0.;
    double tempStressWork = 0.;
    double dissWork = 0.;
    double tempDissWork = 0.;
#endif
 
public:
    ConcreteDPM2Status(GaussPoint *gp);
 
    void initTempStatus() override;
    void updateYourself(TimeStep *tStep) override;
    void printOutputAt(FILE *file, TimeStep *tStep) const override;
    void saveContext(DataStream &stream, ContextMode mode) override;
    void restoreContext(DataStream &stream, ContextMode mode) override;
    const char *giveClassName() const override { return "ConcreteDPM2Status"; }
 
 
    // Inline functions for access to state variables

    const FloatArrayF< 6 > &giveReducedStrain() const { return reducedStrain; }

    const FloatArrayF< 6 > &giveTempReducedStrain() const { return tempReducedStrain; }
 

    const FloatArrayF< 6 > &giveTempEffectiveStress() const { return tempEffectiveStress; }

    const FloatArrayF< 6 > &givePlasticStrain() const { return plasticStrain; }

    double giveDeviatoricPlasticStrainNorm() const
    {
        auto dev = StructuralMaterial::computeDeviator(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 FloatArrayF< 6 > &giveTempPlasticStrain() const { return tempPlasticStrain; }

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

    double giveTempKappaP() const
    { return tempKappaP; }
 
    double giveDeltaLambda() const
    { return deltaLambda; }

    double giveKappaDTension() const
    { return kappaDTension; }

    double giveAlpha() const
    { return alpha; }

    double giveTempAlpha() const
    { return tempAlpha; }
 
 

    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 FloatArrayF< 6 > &v)
    { tempPlasticStrain = v; }
 
 
    void letTempReducedStrainBe(const FloatArrayF< 6 > &v)
    { tempReducedStrain = v; }
 
    void letTempEffectiveStressBe(const FloatArrayF< 6 > &v)
    { tempEffectiveStress = v; }
 

    void letTempKappaPBe(double v)
    { tempKappaP = v; }
 
    void letDeltaLambdaBe(double v)
    { deltaLambda = 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() const { 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 CONCRETE DAMAGE PLASTICITY MODEL 2   ***
//   **************************************************

class ConcreteDPM2 : public StructuralMaterial
{
public:
 
    enum ConcreteDPM2_ReturnResult {
        RR_Unknown,
        RR_NotConverged,
        RR_Converged
    };
 
    enum ConcreteDPM2_ReturnType {
        RT_Unknown,
        RT_Regular,
        RT_Tension,
        RT_Compression,
        RT_Auxiliary
    };
 
 
protected:
    double fc = 0., ft = 0., ecc = 0.;

    int damageFlag = 0;
 
    double e0 = 0.;

    double AHard = 0.;
    double BHard = 0.;
    double CHard = 0.;
    double DHard = 0.;

    double hardeningModulus = 0.;

    double ASoft = 0.;

    double yieldHardPrimePeak = 0.;

    double yieldHardInitial = 0.;

    double dilationConst = 0.;

    double m = 0.;

    double mQ = 0.;

    double helem = 0.;

    IsotropicLinearElasticMaterial linearElasticMaterial;

    double eM = 0.;
    double gM = 0.;
    double kM = 0.;
    double nu = 0.;

    double efCompression = 0.;

    double wf = 0.;

    double wfOne = 0.;

    double ftOne = 0.;

    double yieldTol = 0.;

    double yieldTolDamage = 0.;

    int newtonIter = 0;

    int softeningType = 0;

    double deltaTime = 0.;

    int strengthRateType = 0;

    int energyRateType = 0;
 
public:
    ConcreteDPM2(int n, Domain *d);
 
    void initializeFrom(InputRecord &ir) override;
 
    const char *giveClassName() const override { return "ConcreteDPM2"; }
    const char *giveInputRecordName() const override { return _IFT_ConcreteDPM2_Name; }
 
    FloatArrayF< 6 >giveRealStressVector_3d(const FloatArrayF< 6 > &strain, GaussPoint *gp, TimeStep *tStep) const override;
 
    bool hasMaterialModeCapability(MaterialMode mode) const override;

    FloatArrayF< 6 >performPlasticityReturn(GaussPoint *gp,
                                            const FloatMatrixF< 6, 6 > &D,
                                            const FloatArrayF< 6 > &strain) const;

    void checkForVertexCase(double &answer,
                            ConcreteDPM2_ReturnType &returnType,
                            double sig,
                            double tempKappa,
                            GaussPoint *gp) const;

    double performRegularReturn(FloatArrayF< 6 > &stress,
                                ConcreteDPM2_ReturnResult &returnResult,
                                ConcreteDPM2_ReturnType &returnType,
                                double kappaP,
                                GaussPoint *gp,
                                double theta) const;

    FloatMatrixF< 4, 4 >computeJacobian(double sig,
                                        double rho,
                                        double theta,
                                        double tempKappa,
                                        double deltaLambda,
                                        GaussPoint *gp) const;

    double performVertexReturn(FloatArrayF< 6 > &stress,
                               ConcreteDPM2_ReturnResult &returnResult,
                               ConcreteDPM2_ReturnType &returnType,
                               double apexStress,
                               double tempKappaP,
                               GaussPoint *gp) const;

    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 theta,
                           double tempKappa) const;

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

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

    FloatArrayF< 2 >computeDDuctilityMeasureDInv(double sig,
                                                 double rho,
                                                 double theta,
                                                 double tempKappa) const;

    FloatArrayF< 2 >computeDGDInv(double sig,
                                  double rho,
                                  double tempKappa) const;

    double computeRatioPotential(double sig,
                                 double rho,
                                 double tempKappa) const;

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

    FloatMatrixF< 2, 2 >computeDDGDDInv(double sig,
                                        double rho,
                                        double tempKappa) const;
 

    FloatMatrixF< 6, 6 >computeDDGDDStress(const FloatArrayF< 6 > &stress,
                                           const double tempKappa) const;
 

    FloatArrayF< 6 >computeDCosThetaDStress(const FloatArrayF< 6 > &stress) const;

    FloatMatrixF< 6, 6 >computeDDCosThetaDDStress(const FloatArrayF< 6 > &stress) const;

    FloatArrayF< 2 >computeDDGDInvDKappa(double sig,
                                         double rho,
                                         double tempKappa) const;
 

    FloatArrayF< 2 >computeDDKappaDDeltaLambdaDInv(double sig,
                                                   double rho,
                                                   double theta,
                                                   double tempKappa) const;

    FloatArrayF< 6 >computeDDKappaDDeltaLambdaDStress(const FloatArrayF< 6 > &stress, double tempKappa) const;
 

    FloatArrayF< 6 >computeDDGDStressDKappa(const FloatArrayF< 6 > &stress, double tempKappa) const;

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

    FloatArrayF< 2 >computeDFDInv(double sig,
                                  double rho,
                                  double theta,
                                  double tempKappa) const;

    FloatArrayF< 6 >computeDFDStress(const FloatArrayF< 6 > &stress,
                                     const double tempKappa) const;

    FloatArrayF< 6 >computeDGDStress(const FloatArrayF< 6 > &stress,
                                     const double tempKappa) const;

    FloatMatrixF< 8, 8 >computeFullJacobian(const FloatArrayF< 6 > &stress,
                                            const double deltaLambda,
                                            GaussPoint *gp,
                                            TimeStep *atTime,
                                            const double tempKappa) const;

    double computeTempKappa(double kappaInitial,
                            double sigTrial,
                            double rhoTrial,
                            double sig) const;

    FloatArrayF< 2 >computeDamage(const FloatArrayF< 6 > &strain, const FloatMatrixF< 6, 6 > &D, double timeFactor, GaussPoint *gp, TimeStep *tStep, double alpha, const FloatArrayF< 6 > &effectiveStress) const;
 

    int checkForUnAndReloading(double &tempEquivStrain,
                               double &minEquivStrain,
                               const FloatMatrixF< 6, 6 > &D,
                               GaussPoint *gp) const;

    double computeAlpha(FloatArrayF< 6 > &effectiveStressTension, FloatArrayF< 6 > &effectiveStressCompression, const FloatArrayF< 6 > &effectiveStress) const;

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

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

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

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

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

    double computeDuctilityMeasureDamage(GaussPoint *gp, const double sig, const double rho) const;

    void initDamaged(double kappa,
                     const FloatArrayF< 6 > &strain,
                     GaussPoint *gp) const;

    void computeCoordinates(const FloatArrayF< 6 > &stress, double &sig, double &rho, double &theta) const;
 

    void assignStateFlag(GaussPoint *gp) const;

    FloatArrayF< 6 >computeDRhoDStress(const FloatArrayF< 6 > &stress) const;

    double computeDRDCosTheta(const double theta, const double ecc) const;

    double computeDDRDDCosTheta(const double theta, const double ecc) const;
 

    FloatArrayF< 6 >computeDSigDStress() const;

    FloatMatrixF< 6, 6 >computeDDRhoDDStress(const FloatArrayF< 6 > &stress) const;
 
    FloatMatrixF< 6, 6 >give3dMaterialStiffnessMatrix(MatResponseMode mode, GaussPoint *gp, TimeStep *tStep) const override;

    FloatMatrixF< 6, 6 >compute3dSecantStiffness(GaussPoint *gp, TimeStep *tStep) const;

    FloatMatrixF< 6, 6 >compute3dTangentStiffness(GaussPoint *gp, TimeStep *tStep) const;
 
    bool isCharacteristicMtrxSymmetric(MatResponseMode rMode) const override { return false; }
 
    int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep) override;
 
    void saveContext(DataStream &stream, ContextMode mode) override;
    void restoreContext(DataStream &stream, ContextMode mode) override;
 
 
protected:
    std::unique_ptr<MaterialStatus> CreateStatus(GaussPoint *gp) const override;
};
} //end namespace oofem
#endif