rcsd.C

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
 */
 
#include "rcsd.h"
#include "gausspoint.h"
#include "floatmatrix.h"
#include "floatarray.h"
#include "mathfem.h"
#include "sm/Materials/isolinearelasticmaterial.h"
#include "datastream.h"
#include "contextioerr.h"
#include "classfactory.h"
 
#include <cstring>
 
namespace oofem {
REGISTER_Material(RCSDMaterial);
 
RCSDMaterial :: RCSDMaterial(int n, Domain *d) : RCM2Material(n, d)
{
    linearElasticMaterial = new IsotropicLinearElasticMaterial(n, d);
}
 
 
RCSDMaterial :: ~RCSDMaterial()
{
    delete linearElasticMaterial;
}
 
 
void
RCSDMaterial :: giveRealStressVector(FloatArray &answer, GaussPoint *gp,
                                     const FloatArray &totalStrain,
                                     TimeStep *tStep) const
//
// returns real stress vector in 3d stress space of receiver according to
// previous level of stress and current
// strain increment, the only way, how to correctly update gp records
//
{
    FloatMatrix Ds0;
    double equivStrain;
    FloatArray princStress, reducedAnswer;
    FloatArray reducedStrainVector, strainVector, principalStrain;
    FloatArray reducedSpaceStressVector;
    FloatMatrix tempCrackDirs;
    RCSDMaterialStatus *status = static_cast< RCSDMaterialStatus * >( this->giveStatus(gp) );
 
    this->initTempStatus(gp);
 
    // subtract stress independent part
    // note: eigenStrains (temperature) is not contained in mechanical strain stored in gp
    // therefore it is necessary to subtract always the total eigen strain value
    this->giveStressDependentPartOfStrainVector(reducedStrainVector, gp, totalStrain,
                                                tStep, VM_Total);
 
    StructuralMaterial :: giveFullSymVectorForm( strainVector, reducedStrainVector, gp->giveMaterialMode() );
 
    tempCrackDirs = status->giveTempCrackDirs();
    this->computePrincipalValDir(principalStrain, tempCrackDirs,
                                 strainVector,
                                 principal_strain);
 
 
    if ( status->giveTempMode() == RCSDMaterialStatus :: rcMode ) {
        // rotating crack mode
 
        this->giveRealPrincipalStressVector3d(princStress, gp, principalStrain, tempCrackDirs, tStep);
        princStress.resize(6);
        tempCrackDirs = status->giveTempCrackDirs();
        answer = this->transformStressVectorTo(tempCrackDirs, princStress, 1);
 
        StructuralMaterial :: giveReducedSymVectorForm( reducedSpaceStressVector, answer, gp->giveMaterialMode() );
        status->letTempStressVectorBe(reducedSpaceStressVector);
 
        this->updateCrackStatus(gp, status->giveCrackStrainVector());
 
        StructuralMaterial :: giveReducedSymVectorForm( reducedAnswer, answer, gp->giveMaterialMode() );
        answer = reducedAnswer;
 
        // test if transition to scalar damage mode take place
        double minSofteningPrincStress = this->Ft, dCoeff, CurrFt, E, ep, ef, damage;
        int ipos = 0;
        for ( int i = 1; i <= 3; i++ ) {
            if ( status->giveTempCrackStatus(i) == pscm_SOFTENING ) {
                if ( princStress.at(i) < minSofteningPrincStress ) {
                    minSofteningPrincStress = princStress.at(i);
                    ipos = i;
                }
            }
        }
 
        CurrFt = this->computeStrength( gp, status->giveCharLength(ipos) );
 
        if ( minSofteningPrincStress <= this->SDTransitionCoeff * CurrFt ) {
            // sd transition takes place
 
            E = linearElasticMaterial->give(Ex, gp);
            ep = CurrFt / E;
            ef = this->giveMinCrackStrainsForFullyOpenCrack(gp, ipos);
            dCoeff = ( E / ( princStress.at(ipos) / principalStrain.at(ipos) ) );
            this->giveMaterialStiffnessMatrix(Ds0, SecantStiffness, gp, tStep);
            // compute reached equivalent strain
            equivStrain = this->computeCurrEquivStrain(gp, reducedStrainVector, E, tStep);
            damage = this->computeDamageCoeff(equivStrain, dCoeff, ep, ef);
 
 
            status->setDamageEpsfCoeff(ef);
            status->setDamageEpspCoeff(ep);
            status->setDamageStiffCoeff(dCoeff);
            status->setDs0Matrix(Ds0);
            status->setTempMaxEquivStrain(equivStrain);
            status->setTempDamageCoeff(damage);
            status->setTempMode(RCSDMaterialStatus :: sdMode);
        }
    } else {
        // scalar damage mode
        double ep, ef, E, dCoeff;
        double damage;
        //int ipos;
 
        E = linearElasticMaterial->give(Ex, gp);
        equivStrain = this->computeCurrEquivStrain(gp, reducedStrainVector, E, tStep);
        equivStrain = max( equivStrain, status->giveTempMaxEquivStrain() );
        reducedSpaceStressVector.beProductOf(* status->giveDs0Matrix(), reducedStrainVector);
        dCoeff = status->giveDamageStiffCoeff();
        ef = status->giveDamageEpsfCoeff();
        ep = status->giveDamageEpspCoeff();
        damage = this->computeDamageCoeff(equivStrain, dCoeff, ep, ef);
        reducedSpaceStressVector.times(1.0 - damage);
 
        answer = reducedSpaceStressVector;
 
        status->letTempStressVectorBe(reducedSpaceStressVector);
 
        status->setTempMaxEquivStrain(equivStrain);
        status->setTempDamageCoeff(damage);
    }
 
    status->letTempStrainVectorBe(totalStrain);
}
 
 
void
RCSDMaterial :: giveEffectiveMaterialStiffnessMatrix(FloatMatrix &answer,
                                                     MatResponseMode rMode, GaussPoint *gp,
                                                     TimeStep *tStep) const
//
// returns effective material stiffness matrix in full form
// for gp stress strain mode
//
{
    RCSDMaterialStatus *status = static_cast< RCSDMaterialStatus * >( this->giveStatus(gp) );
 
    if ( status->giveTempMode() == RCSDMaterialStatus :: rcMode ) {
        // rotating crack mode
 
        RCM2Material :: giveEffectiveMaterialStiffnessMatrix(answer, rMode, gp, tStep);
        return;
    } else {
        // rcsd mode
 
        if ( ( rMode == TangentStiffness ) || ( rMode == SecantStiffness ) ) {
            FloatMatrix reducedAnswer;
            double dCoeff;
 
            reducedAnswer = * status->giveDs0Matrix();
            dCoeff = 1.0 - status->giveDamageCoeff();
            if ( dCoeff < RCSD_DAMAGE_EPS ) {
                dCoeff = RCSD_DAMAGE_EPS;
            }
 
            reducedAnswer.times(dCoeff);
 
            answer = reducedAnswer;
        } else if ( rMode == ElasticStiffness ) {
            this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, rMode, gp, tStep);
            return;
        } else {
            OOFEM_ERROR("unsupported mode");
        }
    }
}
 
 
double
RCSDMaterial :: computeDamageCoeff(double equivStrain, double dStiffCoeff, double ep, double ef) const
{
    double damage = 1. - dStiffCoeff * ( 1. - ef / equivStrain ) / ( 1. - ef / ep );
    if ( damage > 1.0 ) {
        return 1.0;
    }
 
    return damage;
}
 
 
double
RCSDMaterial :: computeCurrEquivStrain(GaussPoint *gp, const FloatArray &reducedTotalStrainVector, double e, TimeStep *tStep) const
{
    FloatArray effStress, princEffStress, fullEffStress;
    FloatMatrix De;
    double answer = 0.0;
 
    linearElasticMaterial->giveStiffnessMatrix(De, TangentStiffness, gp, tStep);
    effStress.beProductOf(De, reducedTotalStrainVector);
    StructuralMaterial :: giveFullSymVectorForm( fullEffStress, effStress, gp->giveMaterialMode() );
 
    this->computePrincipalValues(princEffStress, fullEffStress, principal_stress);
    for ( int i = 1; i <= 3; i++ ) {
        answer = max( answer, macbra( princEffStress.at(i) ) );
    }
 
    return answer / e;
}
 
 
void
RCSDMaterial :: initializeFrom(InputRecord &ir)
{
    RCM2Material :: initializeFrom(ir);
    IR_GIVE_FIELD(ir, SDTransitionCoeff, _IFT_RCSDMaterial_sdtransitioncoeff);
}
 
 
double
RCSDMaterial :: give(int aProperty, GaussPoint *gp) const
// Returns the value of the property aProperty (e.g. the Young's modulus
// 'E') of the receiver.
{
    if ( aProperty == pscm_SDTransitionCoeff ) {
        return this->SDTransitionCoeff;
    }
 
    return RCM2Material :: give(aProperty, gp);
}
 
 
int
RCSDMaterial :: checkSizeLimit(GaussPoint *gp, double charLength) const
//
// checks if element size (charLength) is too big
// so that tension strength must be reduced followed
// by sudden stress drop
//
{
    double Ee, Gf, Ft, LeCrit;
 
    Ee = this->give(pscm_Ee, gp);
    Gf = this->give(pscm_Gf, gp);
    Ft = this->give(pscm_Ft, gp);
 
    LeCrit = 2.0 * Gf * Ee / ( Ft * Ft );
    return ( charLength < LeCrit );
}
 
 
double
RCSDMaterial :: computeStrength(GaussPoint *gp, double charLength) const
//
// computes strength for given gp,
// which may be reduced according to length of "fracture process zone"
// to be energetically correct
//
{
    double Ee, Gf, Ft;
 
    Ee = this->give(pscm_Ee, gp);
    Gf = this->give(pscm_Gf, gp);
    Ft = this->give(pscm_Ft, gp);
 
    if ( !this->checkSizeLimit(gp, charLength) ) {
        // we reduce Ft and there is no softening but sudden drop
        Ft = sqrt(2. * Ee * Gf / charLength);
        //
        OOFEM_LOG_RELEVANT("Reducing Ft to %f in element %d, gp %d, Le %f\n",
                           Ft, gp->giveElement()->giveNumber(), gp->giveNumber(), charLength);
    }
 
    return Ft;
}
 
 
double
RCSDMaterial :: giveMinCrackStrainsForFullyOpenCrack(GaussPoint *gp, int i) const
//
// computes MinCrackStrainsForFullyOpenCrack for given gp and i-th crack
//
{
    RCM2MaterialStatus *status = static_cast< RCM2MaterialStatus * >( this->giveStatus(gp) );
    double Le, Gf, Ft;
 
    Le = status->giveCharLength(i);
    Gf = this->give(pscm_Gf, gp);
    Ft = this->computeStrength(gp, Le);
 
    return 2.0 * Gf / ( Le * Ft );
}
 
 
/*
 * void
 * RCSDMaterial :: updateStatusForNewCrack (GaussPoint* gp, int i, double Le)
 * //
 * // updates gp status when new crack-plane i is formed with charLength Le
 * // updates Le and computes and sets minEffStrainForFullyOpenCrack
 * //
 * {
 * RCM2MaterialStatus *status = (RCM2MaterialStatus*) this -> giveStatus (gp);
 *
 * if (Le <= 0) {
 * char errMsg [80];
 * sprintf (errMsg,"Element %d returned zero char length",
 *    gp->giveElement()->giveNumber());
 * OOFEM_ERROR(errMsg);
 * }
 *
 * status -> setCharLength(i, Le);
 * status ->
 * setMinCrackStrainsForFullyOpenCrack (i, this->giveMinCrackStrainsForFullyOpenCrack(gp,i));
 * }
 */
 
double
RCSDMaterial :: giveCrackingModulus(MatResponseMode rMode, GaussPoint *gp,
                                    double crackStrain, int i) const
//
// returns current cracking modulus according to crackStrain for i-th
// crackplane
// now linear softening is implemented
// see also CreateStatus () function.
// softening modulus represents a relation between the normal crack strain
// rate and the normal stress rate.
//
{
    //double Ee, Gf;
    double Cf, Ft, Le, minEffStrainForFullyOpenCrack;
    RCM2MaterialStatus *status = static_cast< RCM2MaterialStatus * >( this->giveStatus(gp) );
 
    //
    // now we have to set proper reduced strength and softening modulus Et
    // in order to obtain energetically correct solution using the concept
    // of fracture energy Gf, which is a material constant.
    //
    //Ee = this->give(pscm_Ee);
    //Gf = this->give(pscm_Gf);
    Ft = this->give(pscm_Ft, gp);
    Le = status->giveCharLength(i);
 
    Ft = this->computeStrength(gp, Le); 
    minEffStrainForFullyOpenCrack = this->giveMinCrackStrainsForFullyOpenCrack(gp, i);
 
    if ( rMode == TangentStiffness ) {
        if ( this->checkSizeLimit(gp, Le) ) {
            if ( ( crackStrain >= minEffStrainForFullyOpenCrack ) ||
                ( status->giveTempMaxCrackStrain(i) >= minEffStrainForFullyOpenCrack ) ) {
                // fully open crack - no stiffness
                Cf = 0.;
            } else if ( crackStrain >= status->giveTempMaxCrackStrain(i) ) {
                // further softening
                Cf = -Ft / minEffStrainForFullyOpenCrack;
            } else {
                // unloading or reloading regime
                Cf = Ft * ( minEffStrainForFullyOpenCrack - status->giveTempMaxCrackStrain(i) ) /
                     ( status->giveTempMaxCrackStrain(i) * minEffStrainForFullyOpenCrack );
            }
        } else {
            Cf = 0.;
        }
    } else {
        if ( this->checkSizeLimit(gp, Le) ) {
            if ( ( crackStrain >= minEffStrainForFullyOpenCrack ) ||
                ( status->giveTempMaxCrackStrain(i) >= minEffStrainForFullyOpenCrack ) ) {
                // fully open crack - no stiffness
                Cf = 0.;
            } else {
                // unloading or reloading regime
                Cf = Ft * ( minEffStrainForFullyOpenCrack - status->giveTempMaxCrackStrain(i) ) /
                     ( status->giveTempMaxCrackStrain(i) * minEffStrainForFullyOpenCrack );
            }
        } else {
            Cf = 0.;
        }
    }
 
    return Cf;
}
 
 
 
double
RCSDMaterial :: giveNormalCrackingStress(GaussPoint *gp, double crackStrain, int i) const
//
// returns receivers Normal Stress in crack i  for given cracking strain
//
{
    double Cf, Ft, Le, answer, minEffStrainForFullyOpenCrack;
    RCM2MaterialStatus *status = static_cast< RCM2MaterialStatus * >( this->giveStatus(gp) );
    minEffStrainForFullyOpenCrack = this->giveMinCrackStrainsForFullyOpenCrack(gp, i);
 
    Cf = this->giveCrackingModulus(TangentStiffness, gp, crackStrain, i); // < 0
    Le = status->giveCharLength(i);
    Ft = this->computeStrength(gp, Le);
 
    if ( this->checkSizeLimit(gp, Le) ) {
        if ( ( crackStrain >= minEffStrainForFullyOpenCrack ) ||
            ( status->giveTempMaxCrackStrain(i) >= minEffStrainForFullyOpenCrack ) ) {
            // fully open crack - no stiffness
            answer = 0.;
        } else if ( crackStrain >= status->giveTempMaxCrackStrain(i) ) {
            // further softening
            answer = Ft + Cf * crackStrain;
        } else if ( crackStrain <= 0. ) {
            // crack closing
            // no stress due to cracking
            answer = 0.;
        } else {
            // unloading or reloading regime
            answer = crackStrain * Ft *
                     ( minEffStrainForFullyOpenCrack - status->giveTempMaxCrackStrain(i) ) /
                     ( status->giveTempMaxCrackStrain(i) * minEffStrainForFullyOpenCrack );
        }
    } else {
        answer = 0.;
    }
 
    return answer;
}
 
 
 
 
RCSDMaterialStatus :: RCSDMaterialStatus(GaussPoint *g) :
    RCM2MaterialStatus(g)
{}
 
 
void
RCSDMaterialStatus :: printOutputAt(FILE *file, TimeStep *tStep) const
{
    char s [ 11 ];
 
    StructuralMaterialStatus :: printOutputAt(file, tStep);
    fprintf(file, "status { ");
    if ( this->giveTempMode() == rcMode ) {
        fprintf(file, "mode :rc ");
 
        if ( this->giveTempAlreadyCrack() ) {
            for ( int i = 1; i <= 3; i++ ) {
                switch ( crackStatuses.at(i) ) {
                case pscm_NONE:
                    strcpy(s, "NONE");
                    break;
                case pscm_OPEN:
                    strcpy(s, "OPEN");
                    break;
                case pscm_SOFTENING:
                    strcpy(s, "SOFTENING");
                    break;
                case pscm_RELOADING:
                    strcpy(s, "RELOADING");
                    break;
                case pscm_UNLOADING:
                    strcpy(s, "UNLOADING");
                    break;
                default:
                    strcpy(s, "UNKNOWN");
                    break;
                }
 
                fprintf( file, "crack %d {status %s, normal to crackplane { %f %f %f }} ",
                        i, s, crackDirs.at(1, i), crackDirs.at(2, i), crackDirs.at(3, i) );
            }
        }
    } else {
        // sd mode
        fprintf( file, "mode :sd, damageCoeff = %f ", this->giveTempDamageCoeff() );
    }
 
    fprintf(file, "}\n");
}
 
 
void
RCSDMaterialStatus :: initTempStatus()
//
// initializes temp variables according to variables form previous equlibrium state.
// builds new crackMap
//
{
    RCM2MaterialStatus :: initTempStatus();
 
    tempMaxEquivStrain = maxEquivStrain;
    tempDamageCoeff = damageCoeff;
    tempMode = mode;
}
 
 
void
RCSDMaterialStatus :: updateYourself(TimeStep *tStep)
//
// updates variables (nonTemp variables describing situation at previous equilibrium state)
// after a new equilibrium state has been reached
// temporary variables are having values corresponding to newly reched equilibrium.
//
{
    RCM2MaterialStatus :: updateYourself(tStep);
 
    maxEquivStrain = tempMaxEquivStrain;
    damageCoeff = tempDamageCoeff;
    mode = tempMode;
}
 
 
void
RCSDMaterialStatus :: saveContext(DataStream &stream, ContextMode mode)
{
    RCM2MaterialStatus :: saveContext(stream, mode);
 
    if ( !stream.write(maxEquivStrain) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( !stream.write(damageCoeff) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( !stream.write(mode) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    contextIOResultType iores;
    if ( ( iores = Ds0.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
}
 
 
void
RCSDMaterialStatus :: restoreContext(DataStream &stream, ContextMode mode)
{
    RCM2MaterialStatus :: restoreContext(stream, mode);
 
    contextIOResultType iores;
    if ( !stream.read(maxEquivStrain) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( !stream.read(damageCoeff) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( !stream.read(mode) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( ( iores = Ds0.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
}
} // end namespace oofem