lsmastermatgrad.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 "lsmastermatgrad.h"
#include "Materials/isolinearelasticmaterial.h"
#include "gausspoint.h"
#include "floatmatrix.h"
#include "floatarray.h"
#include "intarray.h"
#include "stressvector.h"
#include "strainvector.h"
#include "mathfem.h"
#include "contextioerr.h"
#include "datastream.h"
#include "classfactory.h"
 
namespace oofem {
REGISTER_Material(LargeStrainMasterMaterialGrad);
 
// constructor
LargeStrainMasterMaterialGrad :: LargeStrainMasterMaterialGrad(int n, Domain *d) : LargeStrainMasterMaterial(n, d), GradientDamageMaterialExtensionInterface(d)
{
}
 
 
// specifies whether a given material mode is supported by this model
bool
LargeStrainMasterMaterialGrad :: hasMaterialModeCapability(MaterialMode mode) const
{
    return mode == _3dMat;
}
 
// creates a new material status  corresponding to this class
std::unique_ptr<MaterialStatus> 
LargeStrainMasterMaterialGrad :: CreateStatus(GaussPoint *gp) const
{
    return std::make_unique<LargeStrainMasterMaterialStatus>(gp, this->giveDomain(), slaveMat);
}
 
 
void
LargeStrainMasterMaterialGrad :: giveStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep) const
//
// Returns characteristic material stiffness matrix of the receiver
//
{
    OOFEM_ERROR("shouldn't be called");
}
 
 
void
LargeStrainMasterMaterialGrad :: giveGradientDamageStiffnessMatrix_uu(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    MaterialMode mMode = gp->giveMaterialMode();
    switch ( mMode ) {
    case _3dMat:
        answer = LargeStrainMasterMaterial :: give3dMaterialStiffnessMatrix(mode, gp, tStep);
        break;
    default:
        OOFEM_ERROR( "unknown mode (%s)", __MaterialModeToString(mMode) );
    }
}
 
void
LargeStrainMasterMaterialGrad :: giveGradientDamageStiffnessMatrix_du(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    MaterialMode mMode = gp->giveMaterialMode();
    switch ( mMode ) {
    case _3dMat:
        give3dKappaMatrix(answer, mode, gp, tStep);
        break;
    default:
        OOFEM_ERROR( "unknown mode (%s)", __MaterialModeToString(mMode) );
    }
}
 
void
LargeStrainMasterMaterialGrad :: giveGradientDamageStiffnessMatrix_ud(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    MaterialMode mMode = gp->giveMaterialMode();
    switch ( mMode ) {
    case _3dMat:
        give3dGprime(answer, mode, gp, tStep);
        break;
    default:
        OOFEM_ERROR( "unknown mode (%s)", __MaterialModeToString(mMode) );
    }
}
 
 
void
LargeStrainMasterMaterialGrad :: giveGradientDamageStiffnessMatrix_dd_BB(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{}
 
 
void
LargeStrainMasterMaterialGrad :: give3dGprime(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    LargeStrainMasterMaterialStatus *status = static_cast< LargeStrainMasterMaterialStatus * >( this->giveStatus(gp) );
    this->initTempStatus(gp);
    FloatMatrix gPrime;
    GradientDamageMaterialExtensionInterface *graddpmat = dynamic_cast< GradientDamageMaterialExtensionInterface * >( domain->giveMaterial(slaveMat)->giveInterface(GradientDamageMaterialExtensionInterfaceType) );
    if ( graddpmat == NULL ) {
        OOFEM_WARNING("material %d has no Structural support", slaveMat);
        return;
    }
 
    graddpmat->giveGradientDamageStiffnessMatrix_ud(gPrime, mode, gp, tStep);
    gPrime.at(4, 1) =  2. * gPrime.at(4, 1);
    gPrime.at(5, 1) =  2. * gPrime.at(5, 1);
    gPrime.at(6, 1) =  2. * gPrime.at(6, 1);
    answer.beProductOf(status->givePmatrix(), gPrime);
}
 
 
void
LargeStrainMasterMaterialGrad :: give3dKappaMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    LargeStrainMasterMaterialStatus *status = static_cast< LargeStrainMasterMaterialStatus * >( this->giveStatus(gp) );
    this->initTempStatus(gp);
    FloatMatrix kappaMatrix;
    GradientDamageMaterialExtensionInterface *graddpmat = dynamic_cast< GradientDamageMaterialExtensionInterface * >( domain->giveMaterial(slaveMat)->giveInterface(GradientDamageMaterialExtensionInterfaceType) );
    if ( graddpmat == NULL ) {
        OOFEM_WARNING("material %d has no Structural support", slaveMat);
        return;
    }
 
    graddpmat->giveGradientDamageStiffnessMatrix_du(kappaMatrix, mode, gp, tStep);
    kappaMatrix.at(1, 4) = 2. * kappaMatrix.at(1, 4);
    kappaMatrix.at(1, 5) = 2. * kappaMatrix.at(1, 5);
    kappaMatrix.at(1, 6) = 2. * kappaMatrix.at(1, 6);
    answer.beProductTOf( kappaMatrix, status->givePmatrix() );
}
 
 
 
 
void
LargeStrainMasterMaterialGrad :: giveNonlocalInternalForces_N_factor(double &answer, double nlDamageDrivingVariable, GaussPoint *gp, TimeStep *tStep)
{
    answer = nlDamageDrivingVariable;
}
 
void
LargeStrainMasterMaterialGrad :: giveNonlocalInternalForces_B_factor(FloatArray &answer, const FloatArray &nlDamageDrivingVariable_grad, GaussPoint *gp, TimeStep *tStep)
{
    answer = nlDamageDrivingVariable_grad;
    answer.times(internalLength * internalLength);
}
 
 
void
LargeStrainMasterMaterialGrad :: computeLocalDamageDrivingVariable(double &answer, GaussPoint *gp, TimeStep *tStep)
{
    // LargeStrainMasterMaterialStatus *status = static_cast< LargeStrainMasterMaterialStatus * >( this->giveStatus(gp) );
    // @todo: solve this
    //  answer = status->giveTempKappa();w
    answer = 1;
}
 
void
LargeStrainMasterMaterialGrad :: giveFirstPKStressVectorGrad(FloatArray &answer1, double &answer2, GaussPoint *gp, const FloatArray &vF, double nonlocalCumulatedStrain, TimeStep *tStep)
{
    LargeStrainMasterMaterialStatus *status = static_cast< LargeStrainMasterMaterialStatus * >( this->giveStatus(gp) );
    this->initTempStatus(gp);
    MaterialMode mode = gp->giveMaterialMode();
    if  ( mode == _3dMat ) {
        Material *mat;
        StructuralMaterial *sMat;
        mat = domain->giveMaterial(slaveMat);
        sMat = dynamic_cast< StructuralMaterial * >(mat);
        if ( sMat == NULL ) {
            OOFEM_WARNING("material %d has no Structural support", slaveMat);
            return;
        }
 
        GradientDamageMaterialExtensionInterface *dpmat = static_cast< GradientDamageMaterialExtensionInterface * >( sMat->giveInterface(GradientDamageMaterialExtensionInterfaceType) );
        if ( !dpmat ) {
            OOFEM_ERROR("Material doesn't implement the required DpGrad interface!");
        }
 
 
        double lambda1, lambda2, lambda3, E1, E2, E3;
        FloatArray eVals, SethHillStrainVector, stressVector, stressM;
        FloatMatrix F, C, eVecs, SethHillStrain;
        FloatMatrix L1, L2, T;
        //store of deformation gradient into 3x3 matrix
        F.beMatrixForm(vF);
        //compute right Cauchy-Green tensor(C), its eigenvalues and eigenvectors
        C.beTProductOf(F, F);
        // compute eigen values and eigen vectors of C
        C.jaco_(eVals, eVecs, 40);
        // compute Seth - Hill's strain measure, it depends on mParameter
        lambda1 = eVals.at(1);
        lambda2 = eVals.at(2);
        lambda3 = eVals.at(3);
        if ( m == 0 ) {
            E1 = 1. / 2. * log(lambda1);
            E2 = 1. / 2. * log(lambda2);
            E3 = 1. / 2. * log(lambda3);
        } else {
            E1 = 1. / ( 2. * m ) * ( pow(lambda1, m) - 1. );
            E2 = 1. / ( 2. * m ) * ( pow(lambda2, m) - 1. );
            E3 = 1. / ( 2. * m ) * ( pow(lambda3, m) - 1. );
        }
 
        SethHillStrain.resize(3, 3);
        for ( int i = 1; i < 4; i++ ) {
            for ( int j = 1; j < 4; j++ ) {
                SethHillStrain.at(i, j) = E1 * eVecs.at(i, 1) * eVecs.at(j, 1) + E2 *eVecs.at(i, 2) * eVecs.at(j, 2) + E3 *eVecs.at(i, 3) * eVecs.at(j, 3);
            }
        }
 
 
 
        SethHillStrainVector.beSymVectorFormOfStrain(SethHillStrain);
        dpmat->giveRealStressVectorGradientDamage(stressVector, answer2, gp, SethHillStrainVector, nonlocalCumulatedStrain, tStep);
        //        this->constructTransformationMatrix(T, eVecs);
 
        stressVector.at(4) = 2 * stressVector.at(4);
        stressVector.at(5) = 2 * stressVector.at(5);
        stressVector.at(6) = 2 * stressVector.at(6);
 
 
        stressM.beProductOf(T, stressVector);
        stressM.at(4) = 1. / 2. *  stressM.at(4);
        stressM.at(5) = 1. / 2. *  stressM.at(5);
        stressM.at(6) = 1. / 2. *  stressM.at(6);
 
        //        this->constructL1L2TransformationMatrices(L1, L2, eVals, stressM, E1, E2, E3);
 
        FloatMatrix junk, P, TL;
        FloatArray secondPK;
        junk.beProductOf(L1, T);
        P.beTProductOf(T, junk);
        //transformation of the stress to the 2PK stress and then to 1PK
        stressVector.at(4) = 0.5 * stressVector.at(4);
        stressVector.at(5) = 0.5 * stressVector.at(5);
        stressVector.at(6) = 0.5 * stressVector.at(6);
        secondPK.beProductOf(P, stressVector);
        answer1.beProductOf(F, secondPK); // P = F*S
        junk.zero();
        junk.beProductOf(L2, T);
        TL.beTProductOf(T, junk);
 
        status->setPmatrix(P);
        status->setTLmatrix(TL);
        status->letTempStressVectorBe(answer1);
    } else {
        OOFEM_ERROR("Unknown material mode.");
    }
}
 
 
void
LargeStrainMasterMaterialGrad :: initializeFrom(InputRecord &ir)
{
    LargeStrainMasterMaterial :: initializeFrom(ir);
}
} // end namespace oofem