linearelasticmaterial.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 "linearelasticmaterial.h"
#include "gausspoint.h"
#include "sm/CrossSections/simplecrosssection.h"
#include "sm/Materials/structuralms.h"
#include "dynamicinputrecord.h"
#include "floatarrayf.h"
#include "floatmatrixf.h"
 
namespace oofem {
void
LinearElasticMaterial :: initializeFrom(InputRecord &ir)
{
    StructuralMaterial :: initializeFrom(ir);
 
    preCastStiffnessReduction = 0.99999999;
    IR_GIVE_OPTIONAL_FIELD(ir, preCastStiffnessReduction, _IFT_LinearElasticMaterial_preCastStiffRed);
}
 
 
void
LinearElasticMaterial :: giveInputRecord(DynamicInputRecord &input)
{
    StructuralMaterial :: giveInputRecord(input);
    input.setField(this->preCastStiffnessReduction, _IFT_LinearElasticMaterial_preCastStiffRed);
}
 
void
LinearElasticMaterial :: computesSubTangents()
{
    tangentPlaneStrain = {
        tangent(0, 0), tangent(1, 0), tangent(2, 0), tangent(3, 0),
        tangent(0, 1), tangent(1, 1), tangent(2, 1), tangent(3, 1),
        tangent(0, 2), tangent(1, 2), tangent(2, 2), tangent(3, 2),
        tangent(0, 3), tangent(1, 3), tangent(2, 3), tangent(3, 3),
    };
 
    auto c = inv(tangent);
    FloatMatrixF<3,3> reduced = {
        c(0,0), c(0,1), c(0,5),
        c(1,0), c(1,1), c(1,5),
        c(5,0), c(5,1), c(5,5),
    };
    tangentPlaneStress = inv(reduced);
}
 
 
FloatMatrixF<6,6>
LinearElasticMaterial :: give3dMaterialStiffnessMatrix(MatResponseMode mode,
                                                       GaussPoint *gp,
                                                       TimeStep *tStep) const
{
    if ( tStep->giveIntrinsicTime() < this->castingTime ) {
        return tangent *  (1. - this->preCastStiffnessReduction);
    } else {
        return tangent;
    }
}
 
 
FloatArrayF<6>
LinearElasticMaterial :: giveThermalDilatationVector(GaussPoint *gp, TimeStep *tStep) const
{
    return alpha;
}
 
 
FloatArrayF<6>
LinearElasticMaterial :: giveRealStressVector_3d(const FloatArrayF<6> &strain, GaussPoint *gp, TimeStep *tStep) const
{
    auto status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
 
    auto d = this->give3dMaterialStiffnessMatrix(TangentStiffness, gp, tStep);
 
    FloatArrayF<6> stress;
    if ( this->castingTime < 0. ) { // no changes in material stiffness ->> total formulation
        auto thermalStrain = this->computeStressIndependentStrainVector_3d(gp, tStep, VM_Total);
        auto strainVector = strain - thermalStrain;
        stress = dot(d, strainVector);
    } else { // changes in material stiffness ->> incremental formulation
        auto thermalStrain = this->computeStressIndependentStrainVector_3d(gp, tStep, VM_Incremental);
        auto strainIncrement = strain - thermalStrain - FloatArrayF<6>(status->giveStrainVector());
 
        stress = dot(d, strainIncrement) + status->giveStressVector();
    }
 
    // update gp
    status->letTempStrainVectorBe(strain);
    status->letTempStressVectorBe(stress);
    return stress;
}
 
 
FloatArrayF<6>
LinearElasticMaterial :: giveRealStressVector_3dDegeneratedShell(const FloatArrayF<6> &strain, GaussPoint *gp, TimeStep *tStep) const
{
    auto status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
 
    auto d = this->give3dMaterialStiffnessMatrix(TangentStiffness, gp, tStep);
 
    d.at(1, 1) -= d.at(1, 3) * d.at(3, 1) / d.at(3, 3);
    d.at(2, 1) -= d.at(2, 3) * d.at(3, 1) / d.at(3, 3);
    d.at(1, 2) -= d.at(1, 3) * d.at(3, 2) / d.at(3, 3);
    d.at(2, 2) -= d.at(2, 3) * d.at(3, 2) / d.at(3, 3);
 
    d.at(3, 1) = 0.0;
    d.at(3, 2) = 0.0;
    d.at(3, 3) = 0.0;
    d.at(2, 3) = 0.0;
    d.at(1, 3) = 0.0;
 
    FloatArrayF<6> stress, strainVector;
    if ( this->castingTime < 0. ) { // no changes in material stiffness ->> total formulation
        auto thermalStrain = this->computeStressIndependentStrainVector_3d(gp, tStep, VM_Total);
        strainVector = strain - thermalStrain;
        stress = dot(d, strainVector);
    } else { // changes in material stiffness ->> incremental formulation
        auto thermalStrain = this->computeStressIndependentStrainVector_3d(gp, tStep, VM_Incremental);
        auto strainIncrement = strain - thermalStrain - FloatArrayF<6>(status->giveStrainVector());
 
        stress = dot(d, strainIncrement) + status->giveStressVector();
    }
 
    stress = dot(d, strainVector);
 
    // update gp
    status->letTempStrainVectorBe(strain);
    status->letTempStressVectorBe(stress);
    return stress;
}
 
 
FloatArrayF<3>
LinearElasticMaterial :: giveRealStressVector_PlaneStress(const FloatArrayF<3> &reducedStrain, GaussPoint *gp, TimeStep *tStep) const
{
    auto status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
 
    auto d = this->givePlaneStressStiffMtrx(TangentStiffness, gp, tStep);
 
    FloatArray answer;
    FloatArray strainVector;
    if ( this->castingTime < 0. ) { // no changes in material stiffness ->> total formulation
        this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Total);
        answer.beProductOf(d, strainVector);
    } else { // changes in material stiffness ->> incremental formulation
        this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Incremental);
        auto strainIncrement = strainVector - status->giveStrainVector();
 
        answer.beProductOf(d, strainIncrement);
        answer.add( status->giveStressVector() );
    }
 
    // update gp
    status->letTempStrainVectorBe(reducedStrain);
    status->letTempStressVectorBe(answer);
    return answer;
}
 
 
FloatArrayF<1>
LinearElasticMaterial :: giveRealStressVector_1d(const FloatArrayF<1> &reducedStrain, GaussPoint *gp, TimeStep *tStep) const
{
    auto status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
    auto d = this->give1dStressStiffMtrx(TangentStiffness, gp, tStep);
 
    FloatArray answer;
    FloatArray strainVector;
    if ( this->castingTime < 0. ) {  // no changes in material stiffness ->> total formulation
        this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Total);
        answer.beProductOf(d, strainVector);
    } else { // changes in material stiffness ->> incremental formulation
        this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Incremental);
        auto strainIncrement = strainVector - status->giveStrainVector();
 
        answer.beProductOf(d, strainIncrement);
        answer.add( status->giveStressVector() );
    }
 
    // update gp
    status->letTempStrainVectorBe(reducedStrain);
    status->letTempStressVectorBe(answer);
    return answer;
}
 
 
FloatArrayF<2>
LinearElasticMaterial :: giveRealStressVector_Warping(const FloatArrayF<2> &reducedStrain, GaussPoint *gp, TimeStep *tStep) const
{
    // reducedStrain contains the stress components tau_zy and tau_zx, computed as    //        tau_zy = G * theta * ( x + dPsi/dy )
    //        tau_zx = G * theta * (-y + dPsi/dx )
    // where x and y are the global coordinates of the Gauss point (the origin must be at the centroid of the cross section)
    //       G is the shear modulus of elasticity and theta is the relative twist (dPhi_z/dz)
    
    
    auto status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
    double G = this->giveShearModulus();
    FloatArray gcoords;
    Element *elem = gp->giveElement();
    elem->computeGlobalCoordinates( gcoords, gp->giveNaturalCoordinates() );
    FloatArrayF<2> answer;
 
    if ( this->castingTime < 0. ) { // no changes in material stiffness ->> total formulation
        answer.at(1) = reducedStrain.at(1);
        answer.at(2) = reducedStrain.at(2);
 
        answer *= G;
    } else { // changes in material stiffness ->> incremental formulation
        FloatArray strainIncrement;
        strainIncrement.beDifferenceOf( reducedStrain, status->giveStrainVector() );
 
        answer.at(1) = strainIncrement.at(1);
        answer.at(2) = strainIncrement.at(2);
        answer *= G;
 
        if ( ( tStep->giveIntrinsicTime() < this->castingTime ) ) {
            answer *= (1. - this->preCastStiffnessReduction);
        }
 
        answer += FloatArrayF<2>(status->giveStressVector());
    }
 
    // update gp
 
    status->letTempStrainVectorBe(reducedStrain);
    status->letTempStressVectorBe(answer);
    return answer;
}
 
 
FloatArrayF<2>
LinearElasticMaterial :: giveRealStressVector_2dBeamLayer(const FloatArrayF<2> &reducedStrain, GaussPoint *gp, TimeStep *tStep) const
{
    auto status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
 
    auto d = this->give2dBeamLayerStiffMtrx(TangentStiffness, gp, tStep);
 
    FloatArray answer;
    FloatArray strainVector, strainIncrement;
    if ( this->castingTime < 0. ) { // no changes in material stiffness ->> total formulation
        this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Total);
        answer.beProductOf(d, strainVector);
    } else { // changes in material stiffness ->> incremental formulation
        this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Incremental);
        strainIncrement.beDifferenceOf( strainVector, status->giveStrainVector() );
 
        answer.beProductOf(d, strainIncrement);
        answer.add( status->giveStressVector() );
    }
 
    // update gp
    status->letTempStrainVectorBe(reducedStrain);
    status->letTempStressVectorBe(answer);
    return answer;
}
 
 
FloatArrayF<5>
LinearElasticMaterial :: giveRealStressVector_PlateLayer(const FloatArrayF<5> &reducedStrain, GaussPoint *gp, TimeStep *tStep) const
{
    auto status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
 
    auto d = this->givePlateLayerStiffMtrx(TangentStiffness, gp, tStep);
 
    FloatArray answer;
    FloatArray strainVector, strainIncrement;
    if ( this->castingTime < 0. ) { // no changes in material stiffness ->> total formulation
        this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Total);
        answer.beProductOf(d, strainVector);
    } else { // changes in material stiffness ->> incremental formulation
        this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Incremental);
        strainIncrement.beDifferenceOf( strainVector, status->giveStrainVector() );
 
        answer.beProductOf(d, strainIncrement);
        answer.add( status->giveStressVector() );
    }
 
    // update gp
    status->letTempStrainVectorBe(reducedStrain);
    status->letTempStressVectorBe(answer);
    return answer;
}
 
 
FloatArrayF<3>
LinearElasticMaterial :: giveRealStressVector_Fiber(const FloatArrayF<3> &reducedStrain, GaussPoint *gp, TimeStep *tStep) const
{
    auto status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
 
    auto d = this->giveFiberStiffMtrx(TangentStiffness, gp, tStep);
 
    FloatArray answer;
    FloatArray strainVector;
    if ( this->castingTime < 0. ) { // no changes in material stiffness ->> total formulation
        this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Total);
        answer.beProductOf(d, strainVector);
    } else { // changes in material stiffness ->> incremental formulation
        this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Incremental);
        auto strainIncrement = strainVector - status->giveStrainVector();
 
        answer.beProductOf(d, strainIncrement);
        answer.add( status->giveStressVector() );
    }
 
    // update gp
    status->letTempStrainVectorBe(reducedStrain);
    status->letTempStressVectorBe(answer);
    return answer;
}
 
void
LinearElasticMaterial :: giveEshelbyStressVector_PlaneStrain(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedF, TimeStep *tStep)
{
    FloatArray fullFv;
    StructuralMaterial :: giveFullVectorFormF(fullFv, reducedF, _PlaneStrain);
 
    FloatMatrix H;
    H.beMatrixForm(fullFv);
 
    H(0, 0) -= 1.0;
    H(1, 1) -= 1.0;
    H(2, 2) -= 1.0;
 
 
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
    const FloatArray &stressV = status->giveTempStressVector();
 
    FloatArray fullStressV;
    StructuralMaterial :: giveFullSymVectorForm(fullStressV, stressV, _PlaneStrain);
 
    FloatMatrix stress;
    stress.beMatrixFormOfStress(fullStressV);
 
 
    double energyDens = giveEnergyDensity(gp, tStep);
 
 
    FloatMatrix eshelbyStress(3, 3);
    eshelbyStress.beTProductOf(H, stress);
    eshelbyStress.negated();
 
    eshelbyStress(0, 0) += energyDens;
    eshelbyStress(1, 1) += energyDens;
    eshelbyStress(2, 2) += energyDens;
 
 
    FloatArray eshelbyStressV;
    eshelbyStressV.beVectorForm(eshelbyStress);
    StructuralMaterial :: giveReducedVectorForm(answer, eshelbyStressV, _PlaneStrain);
}
 
double
LinearElasticMaterial :: giveEnergyDensity(GaussPoint *gp, TimeStep *tStep)
{
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
    const FloatArray &strain = status->giveTempStrainVector();
    const FloatArray &stress = status->giveTempStressVector();
 
    return 0.5 * stress.dotProduct(strain);
}
 
int
LinearElasticMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
    if ( type == IST_ElasticStrainTensor ) {
        this->giveStressDependentPartOfStrainVector(answer, gp, status->giveStrainVector(), tStep, VM_Total);
        return 1;
    } else if ( type == IST_ThermalStrainTensor ) {
        answer = this->computeStressIndependentStrainVector_3d(gp, tStep, VM_Total);
        return 1;
    } else if ( type == IST_CreepStrainTensor ) {
        answer.resize(6);
        answer.zero();
        return 1;
    } else if ( type == IST_FreeEnergyDensity ) {
        answer.resize(1);
        answer.at(1) = this->giveEnergyDensity(gp,tStep);
        return 1;
    } else {
        return StructuralMaterial :: giveIPValue(answer, gp, type, tStep);
    }
}
 
 
std::unique_ptr<MaterialStatus> 
LinearElasticMaterial :: CreateStatus(GaussPoint *gp) const
{
    return std::make_unique<StructuralMaterialStatus>(gp);
}
} // end namespace oofem