rheoChM.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 rheochm_h
#define rheochm_h

#define keep_track_of_strains
 
#include "sm/Materials/structuralmaterial.h"
//#include "sm/Materials/linearelasticmaterial.h"
#include "floatarray.h"
#include "floatmatrix.h"
 
#include "matconst.h"
#include "sm/Elements/structuralelement.h"
#include "sm/Materials/structuralms.h"


#define _IFT_RheoChainMaterial_n "n"
#define _IFT_RheoChainMaterial_alphaOne "a1"
#define _IFT_RheoChainMaterial_alphaTwo "a2"
#define _IFT_RheoChainMaterial_lattice "lattice"
#define _IFT_RheoChainMaterial_relmatage "relmatage"
#define _IFT_RheoChainMaterial_begoftimeofinterest "begoftimeofinterest"
#define _IFT_RheoChainMaterial_endoftimeofinterest "endoftimeofinterest"
#define _IFT_RheoChainMaterial_timefactor "timefactor"
#define _IFT_RheoChainMaterial_talpha "talpha"
 
namespace oofem {
#define MNC_NPOINTS 30
#define TIME_DIFF   1.e-10

class RheoChainMaterialStatus : public StructuralMaterialStatus
{
protected:
    int nUnits = 0;
    std :: vector< FloatArray >hiddenVars;
    std :: vector< FloatArray >tempHiddenVars;

    FloatArray shrinkageStrain;

    double currentTime = 0.;
 
#ifdef keep_track_of_strains
    double thermalStrain = 0.;
    double tempThermalStrain = 0.;
#endif
 
public:
    RheoChainMaterialStatus(GaussPoint *g, int nunits);
 
    void printOutputAt(FILE *file, TimeStep *tStep) const override;
 
    virtual const FloatArray &giveViscoelasticStressVector() const { return stressVector; }
 
    FloatArray &giveHiddenVarsVector(int i) { return hiddenVars [ i - 1 ]; }
    FloatArray &giveTempHiddenVarsVector(int i) { return tempHiddenVars [ i - 1 ]; }
    FloatArray *letHiddenVarsVectorBe(int i, FloatArray *);
    void letTempHiddenVarsVectorBe(int i, FloatArray &valueArray);
 
    FloatArray *giveShrinkageStrainVector() { return & shrinkageStrain; }
    void setShrinkageStrainVector(FloatArray src) { shrinkageStrain = std :: move(src); }

    double giveCurrentTime() { return currentTime; }
    void setCurrentTime(double src) { currentTime = src; }
 
    void initTempStatus() override;
    void updateYourself(TimeStep *tStep) override;
 
    void saveContext(DataStream &stream, ContextMode mode) override;
    void restoreContext(DataStream &stream, ContextMode mode) override;
 
 
#ifdef keep_track_of_strains
    void setTempThermalStrain(double src) { tempThermalStrain = src; }
    double giveTempThermalStrain(void) { return tempThermalStrain; }
    double giveThermalStrain(void) { return thermalStrain; }
#endif
 
    // definition
    const char *giveClassName() const override { return "RheoChainMaterialStatus"; }
};
 

class RheoChainMaterial : public StructuralMaterial
{
protected:
    double talpha = 0.;
    int nUnits = 0;
    double relMatAge = 0.;
 
    bool lattice = false;
    double nu = 0.;
    double alphaOne = 0., alphaTwo = 0.;
    mutable double EparValTime = -1.;

    double begOfTimeOfInterest = 0.; // local one or taken from e-model
    double endOfTimeOfInterest = 0.; // local one or taken from e-model
    //    LinearElasticMaterial *linearElasticMaterial = nullptr;
    StructuralMaterial *linearElasticMaterial = nullptr;

    mutable FloatArray EparVal;
    //FloatArray relaxationTimes;
    mutable FloatArray charTimes;
    FloatArray discreteTimeScale;

    double timeFactor = 0.;
 
public:
    RheoChainMaterial(int n, Domain *d);
    virtual ~RheoChainMaterial();
 
    void giveRealStressVector(FloatArray &answer, GaussPoint *gp,
                              const FloatArray &reducedStrain, TimeStep *tStep) const override;
 
    FloatArrayF< 6 >giveRealStressVector_3d(const FloatArrayF< 6 > &strain, GaussPoint *gp, TimeStep *tStep) const override
    {
        FloatArray answer;
        const_cast< RheoChainMaterial * >( this )->giveRealStressVector(answer, gp, strain, tStep);
        return answer;
    }
    FloatArrayF< 4 >giveRealStressVector_PlaneStrain(const FloatArrayF< 4 > &strain, GaussPoint *gp, TimeStep *tStep) const override
    {
        FloatArray answer;
        const_cast< RheoChainMaterial * >( this )->giveRealStressVector(answer, gp, strain, tStep);
        return answer;
    }
    FloatArrayF< 3 >giveRealStressVector_PlaneStress(const FloatArrayF< 3 > &strain, GaussPoint *gp, TimeStep *tStep) const override
    {
        FloatArray answer;
        const_cast< RheoChainMaterial * >( this )->giveRealStressVector(answer, gp, strain, tStep);
        return answer;
    }
    FloatArrayF< 1 >giveRealStressVector_1d(const FloatArrayF< 1 > &strain, GaussPoint *gp, TimeStep *tStep) const override
    {
        FloatArray answer;
        const_cast< RheoChainMaterial * >( this )->giveRealStressVector(answer, gp, strain, tStep);
        return answer;
    }
    FloatArrayF< 2 >giveRealStressVector_2dBeamLayer(const FloatArrayF< 2 > &strain, GaussPoint *gp, TimeStep *tStep) const override
    {
        FloatArray answer;
        const_cast< RheoChainMaterial * >( this )->giveRealStressVector(answer, gp, strain, tStep);
        return answer;
    }
    FloatArrayF< 5 >giveRealStressVector_PlateLayer(const FloatArrayF< 5 > &strain, GaussPoint *gp, TimeStep *tStep) const override
    {
        FloatArray answer;
        const_cast< RheoChainMaterial * >( this )->giveRealStressVector(answer, gp, strain, tStep);
        return answer;
    }
 
    FloatArrayF< 6 >giveThermalDilatationVector(GaussPoint *gp, TimeStep *tStep) const override;
 
    /*    virtual void giveThermalDilatationVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
     *    { answer.clear(); }*/

    virtual double giveEModulus(GaussPoint *gp, TimeStep *tStep) const = 0;
 
    /*    virtual double giveIncrementalModulus(GaussPoint *gp, TimeStep *tStep) {
     * if ( (tStep->giveIntrinsicTime() < this->castingTime) && ( this->zeroStiffness > 0. ) ) {
     *  return this->zeroStiffness;
     * } else {
     *  return this->giveEModulus(gp, tStep);
     * }
     * }*/

    virtual FloatArray computeCharCoefficients(double tPrime, GaussPoint *gp, TimeStep *tStep) const = 0;
 
    // identification and auxiliary functions
    bool hasMaterialModeCapability(MaterialMode mode) const override;
 
    bool hasCastingTimeSupport() const override { return true; }
 
    const char *giveClassName() const override { return "RheoChainMaterial"; }
    void initializeFrom(InputRecord &ir) override;
 
    int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep) override;
 
    // store & restore context functions
    void saveIPContext(DataStream &stream, ContextMode mode, GaussPoint *gp) override;
    void restoreIPContext(DataStream &stream, ContextMode mode, GaussPoint *gp) override;

    /*    void giveStiffnessMatrix(FloatMatrix &answer,
     *                               MatResponseMode mode,
     *                               GaussPoint *gp,
     *                               TimeStep *tStep) override;
     */
 
    FloatMatrixF< 6, 6 >give3dMaterialStiffnessMatrix(MatResponseMode mode, GaussPoint *gp, TimeStep *tStep) const override;
    FloatMatrixF< 3, 3 >givePlaneStressStiffMtrx(MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep) const override;
    FloatMatrixF< 4, 4 >givePlaneStrainStiffMtrx(MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep) const override;
    FloatMatrixF< 1, 1 >give1dStressStiffMtrx(MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep) const override;
 
    // maybe not needed afterall?
    //  FloatMatrixF< 3, 3 >give2dLatticeStiffnessMatrix(MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep) const;
 
    //  FloatMatrixF< 6, 6 >give3dLatticeStiffnessMatrix(MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep) const;
 
    FloatArray computeStressIndependentStrainVector(GaussPoint *gp, TimeStep *tStep, ValueModeType mode) const override;

    virtual void giveShrinkageStrainVector(FloatArray &answer,
                                           GaussPoint *gp,
                                           TimeStep *tStep,
                                           ValueModeType mode) const
    { answer.clear(); }
 
    // Note: must take LoadResponseMode into account
    virtual void giveEigenStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep, ValueModeType mode) const { }
 
    std::unique_ptr<MaterialStatus> CreateStatus(GaussPoint *gp) const override;
 
    double giveAlphaOne() const { return this->alphaOne; }
    double giveAlphaTwo() const { return this->alphaTwo; }

    double givePoissonsRatio() const { return nu; }

    virtual double computeCreepFunction(double t, double t_prime, GaussPoint *gp, TimeStep *tStep) const = 0;

    bool isActivated(TimeStep *tStep) const override {
        if ( this->preCastingTimeMat > 0 ) {
            return true;
        } else {
            return Material :: isActivated(tStep);
        }
    }

    virtual double giveEquivalentTime(GaussPoint *gp, TimeStep *tStep) const
    { return ( tStep->giveTargetTime() - tStep->giveTimeIncrement() / 2 ); }
 
 
protected:
    virtual bool hasIncrementalShrinkageFormulation() const { return false; }

    static FloatArray generateLogTimeScale(double from, double to, int nsteps);
    const FloatArray &giveDiscreteTimes() const;

    void computeDiscreteRelaxationFunction(FloatArray &answer, const FloatArray &tSteps, double t0, double tr, GaussPoint *gp, TimeStep *tStep) const;

    void giveUnitComplianceMatrix(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep) const;
    void giveUnitStiffnessMatrix(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep) const;

    virtual void updateEparModuli(double tPrime, GaussPoint *gp, TimeStep *tStep) const;

    double giveEparModulus(int iChain) const;

    virtual void computeCharTimes();

    double giveCharTime(int) const;

    virtual double giveCharTimeExponent(int i) const { return 1.0; }

    //    LinearElasticMaterial *giveLinearElasticMaterial();
    StructuralMaterial *giveLinearElasticMaterial();

    double giveEndOfTimeOfInterest();

    void computeTrueStressIndependentStrainVector(FloatArray &answer, GaussPoint *gp,
                                                  TimeStep *tStep, ValueModeType mode) const;
};
} // end namespace oofem
#endif // rheochm_h