mdm.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 mdm_h
#define mdm_h

#define MDM_USE_MMAClosestIPTransfer
//#define MDM_USE_MMAShapeFunctProjection
//#define MDM_USE_MMALeastSquareProjection
 
 
// Allows to select different mapping algorithms
#define MDM_MAPPING_DEBUG 1
 
#include "microplanematerial.h"
#include "structuralnonlocalmaterialext.h"
#include "matconst.h"
#include "sm/Materials/isolinearelasticmaterial.h"
#include "sm/Materials/structuralms.h"
#include "materialmapperinterface.h"
#include "mmaclosestiptransfer.h"
 
#ifdef MDM_MAPPING_DEBUG
 #include "mmashapefunctprojection.h"
 #include "mmaleastsquareprojection.h"
 
#else
 
 #ifdef MDM_USE_MMAShapeFunctProjection
  #include "mmashapefunctprojection.h"
 #endif
 #ifdef MDM_USE_MMALeastSquareProjection
  #include "mmaleastsquareprojection.h"
 #endif
#endif
 
#include "mmashapefunctprojection.h"
#include "mmaleastsquareprojection.h"


#define _IFT_MDM_Name "mdm"
#define _IFT_MDM_talpha "talpha"
#define _IFT_MDM_parmd "parmd"
#define _IFT_MDM_nonloc "nonloc"
#define _IFT_MDM_r "r"
#define _IFT_MDM_efp "efp"
#define _IFT_MDM_ep "ep"
#define _IFT_MDM_gf "gf"
#define _IFT_MDM_ft "ft"
#define _IFT_MDM_formulation "formulation"
#define _IFT_MDM_mode "mode"
#define _IFT_MDM_mapper "mapper"
#define _IFT_MDM_sourceRegionSet "sourceregset"
 
namespace oofem {

class MDMStatus : public StructuralMaterialStatus, public StructuralNonlocalMaterialStatusExtensionInterface
{
protected:
    FloatArray Psi, PsiTemp;
    FloatMatrix DamageTensor, DamageTensorTemp, localDamageTensor;
    FloatArray tempDamageTensorEigenValues, damageTensorEigenValues;
    FloatMatrix tempDamageTensorEigenVectors, damageTensorEigenVectors;
 
public:
    MDMStatus(GaussPoint * g, int nsd, int nmplanes);
 
    void setTempDamageTensorEigenVals(FloatArray src) { tempDamageTensorEigenValues = std :: move(src); }
    void setTempDamageTensorEigenVec(FloatMatrix src) { tempDamageTensorEigenVectors = std :: move(src); }
    const FloatArray &giveTempDamageTensorEigenVals() { return tempDamageTensorEigenValues; }
    const FloatArray &giveDamageTensorEigenVals() { return damageTensorEigenValues; }
    const FloatMatrix &giveTempDamageTensorEigenVec() { return tempDamageTensorEigenVectors; }
    const FloatMatrix &giveDamageTensorEigenVec() { return damageTensorEigenVectors; }
 
    double giveMicroplaneTempDamage(int m) { return PsiTemp.at(m); }
    double giveMicroplaneDamage(int m) { return Psi.at(m); }
    void setMicroplaneTempDamage(int m, double val) { PsiTemp.at(m) = val; }
    const FloatArray & giveMicroplaneDamageValues() { return Psi; }
    void setMicroplaneTempDamageValues(FloatArray src) { PsiTemp = std :: move(src); }
 
    const FloatMatrix &giveTempDamageTensor() { return DamageTensorTemp; }
    const FloatMatrix &giveDamageTensor() { return DamageTensor; }
    void setTempDamageTensor(FloatMatrix src) { DamageTensorTemp = std :: move(src); }
    void setLocalDamageTensorForAverage(FloatMatrix src) { localDamageTensor = std :: move(src); }
    const FloatMatrix &giveLocalDamageTensorForAverage() { return localDamageTensor; }
    const FloatMatrix *giveLocalDamageTensorForAveragePtr() { return & localDamageTensor; }
 
    const char *giveClassName() const override { return "MDMStatus"; }
 
    void printOutputAt(FILE *file, TimeStep *tStep) const override;
 
    void initTempStatus() override;
    void updateYourself(TimeStep *tStep) override;
 
    void saveContext(DataStream &stream, ContextMode mode) override;
    void restoreContext(DataStream &stream, ContextMode mode) override;

    Interface *giveInterface(InterfaceType it) override;
};
 

class MDM : public MicroplaneMaterial, public StructuralNonlocalMaterialExtensionInterface, public MaterialModelMapperInterface
{
public:
    enum MDMFormulatrionType { COMPLIANCE_DAMAGE, STIFFNESS_DAMAGE };
    enum MDMModeType { mdm_3d, mdm_2d };
 
protected:
    int ndc = 0; 
    int nsd = 0; 
    int type_dam = 0;
    int type_soft = 0;

    mutable double mdm_Ep = 0.;
    mutable double mdm_Efp = 0.; 
    double ParMd = 0.; 
    double tempDillatCoeff = 0.; 

    double Gf = 0.;
    double Ft = 0.;
 
    MDMFormulatrionType formulation = COMPLIANCE_DAMAGE;
    MDMModeType mdmMode = mdm_3d;

    IsotropicLinearElasticMaterial linearElasticMaterial;

    int nonlocal = 0;
    double R = 0.;

    std::unique_ptr<Set> sourceElemSet;
 
#ifdef MDM_MAPPING_DEBUG
    static MMAShapeFunctProjection mapperSFT;
    static MMALeastSquareProjection mapperLST;
 
    // determines the mapping algorithm to be used
    enum MDMMapperType { mdm_cpt=0, mdm_sft=1, mdm_lst=2 };
    MDMMapperType mapperType = mdm_cpt;
#else
 
 #ifdef MDM_USE_MMAClosestIPTransfer
    static MMAClosestIPTransfer mapper;
 #endif
 #ifdef MDM_USE_MMAShapeFunctProjection
    static MMAShapeFunctProjection mapper;
 #endif
 #ifdef MDM_USE_MMALeastSquareProjection
    static MMALeastSquareProjection mapper;
 #endif
#endif

    static MMAClosestIPTransfer mapper2;
 
public:
    MDM(int n, Domain * d) : MicroplaneMaterial(n, d), 
        StructuralNonlocalMaterialExtensionInterface(d),
        MaterialModelMapperInterface(),
        linearElasticMaterial(n, d)
    {
        mdm_Ep = mdm_Efp = -1.0;
    }
 
    bool hasMaterialModeCapability(MaterialMode mode) const override;
 
    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<MDM*>(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<MDM*>(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<MDM*>(this)->giveRealStressVector(answer, gp, strain, tStep);
        return answer;
    }
 
    FloatMatrixF<6,6> give3dMaterialStiffnessMatrix(MatResponseMode mode, GaussPoint *gp, TimeStep *tStep) const override;
 
    int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep) override;
 
    FloatArrayF<6> giveThermalDilatationVector(GaussPoint *gp, TimeStep *tStep) const override;
 
    void initializeFrom(InputRecord &ir) override;
    void giveInputRecord(DynamicInputRecord &input) override;
 
    const char *giveInputRecordName() const override { return _IFT_MDM_Name; }
    const char *giveClassName() const override { return "MDM"; }
 
    void initializeData(int numberOfMicroplanes) override;
 
    void updateBeforeNonlocAverage(const FloatArray &strainVector, GaussPoint *gp, TimeStep *tStep) const override;
 
  double computeWeightFunction(const double cl, const FloatArray &src, const FloatArray &coord) const override;
 
    int hasBoundedSupport() const override { return 1; }
    virtual void giveSupportRadius(double &radius) { radius = this->R; }
 
    Interface *giveInterface(InterfaceType it) override;
 
    int MMI_map(GaussPoint *gp, Domain *oldd, TimeStep *tStep) override;
    int MMI_update(GaussPoint *gp, TimeStep *tStep, FloatArray *estrain = nullptr) override;
    int MMI_finish(TimeStep *tStep) override;
 
    int packUnknowns(DataStream &buff, TimeStep *tStep, GaussPoint *ip) override;
    int unpackAndUpdateUnknowns(DataStream &buff, TimeStep *tStep, GaussPoint *ip) override;
    int estimatePackSize(DataStream &buff, GaussPoint *ip) override;
    double predictRelativeComputationalCost(GaussPoint *gp) override;
    double predictRelativeRedistributionCost(GaussPoint *gp) override { return 1.0; }
 
    std::unique_ptr<MaterialStatus> CreateStatus(GaussPoint *gp) const override;
 
protected:
    void computeDamageTensor(FloatMatrix &tempDamageTensor, const FloatArray &totalStrain,
                             GaussPoint *gp, TimeStep *tStep) const;
    void computeLocalDamageTensor(FloatMatrix &tempDamageTensor, const FloatArray &totalStrain,
                                  GaussPoint *gp, TimeStep *tStep) const;
    double computeDamageOnPlane(GaussPoint *gp, int mnumber, const FloatArray &strain) const;
    void computePDC(FloatMatrix &tempDamageTensor, FloatArray &tempDamageTensorEigenVals,
                    FloatMatrix &tempDamageTensorEigenVec) const;
    void transformStrainToPDC(FloatArray &answer, FloatArray &strain,
                              FloatMatrix &t, GaussPoint *gp) const;
    void applyDamageTranformation(FloatArray &strainPDC, const FloatArray &tempDamageTensorEigenVals) const;
    void transformStressFromPDC(FloatArray &answer, const FloatArray &stressPDC, const FloatMatrix &t, GaussPoint *gp) const;
    void computeEffectiveStress(FloatArray &stressPDC, const FloatArray &strainPDC,
                                GaussPoint *gp, TimeStep *tStep) const;
    void giveMaterialStiffnessMatrix(FloatMatrix &answer,
                                     MatResponseMode mode, GaussPoint *gp,
                                     TimeStep *tStep) const;
    void applyDamageToStiffness(FloatMatrix &d, GaussPoint *gp) const;
    void transformStiffnessfromPDC(FloatMatrix &de, const FloatMatrix &t) const;
 
    FloatMatrixF<3,3> givePlaneStressStiffMtrx(MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep) const override;
 
    FloatMatrixF<4,4> givePlaneStrainStiffMtrx(MatResponseMode mmode, GaussPoint *gp, TimeStep *tStep) const override;
 
    void rotateTensor4(FloatMatrix &Dlocal, const FloatMatrix &t) const;
    void formTransformationMatrix(FloatMatrix &answer, const FloatMatrix &t, int n) const;
 
    //double giveParameterEfp(const FloatArray &reducedStrain, GaussPoint *gp);
    void giveRawMDMParameters(double &Efp, double &Ep, const FloatArray &reducedStrain, GaussPoint *gp) const;
};
} // end namespace oofem
#endif // mdm_h