LocalIsotropicDamage.h

Go to the documentation of this file.

#pragma once

#include "MechanicsInterface.h"
#include "LinearElasticDamage.h"
#include "ModifiedMisesStrainNorm.h"
#include "damageLaws/DamageLawExponential.h"

namespace NuTo
{
namespace Laws
{
template <int TDim>
class EvolutionImplicit;

template <int TDim, typename TDamageLaw = Constitutive::DamageLawExponential,
          typename TEvolution = EvolutionImplicit<TDim>>
class LocalIsotropicDamage : public MechanicsInterface<TDim>
{
public:
    LocalIsotropicDamage(LinearElasticDamage<TDim> elasticDamage, TDamageLaw damageLaw, TEvolution evolution)
        : mElasticDamage(elasticDamage)
        , mDamageLaw(damageLaw)
        , mEvolution(evolution)
    {
    }

    LocalIsotropicDamage(Material::Softening m, eDamageApplication damageApplication = eDamageApplication::FULL)
        : mElasticDamage(m, damageApplication)
        , mDamageLaw(m)
        , mEvolution(m)
    {
    }

    EngineeringStress<TDim> Stress(EngineeringStrain<TDim> strain, double deltaT, CellIds ids) const override
    {
        double kappa = mEvolution.Kappa(strain, deltaT, ids);
        double omega = mDamageLaw.Damage(kappa);
        return mElasticDamage.Stress(strain, omega);
    }

    EngineeringTangent<TDim> Tangent(EngineeringStrain<TDim> strain, double deltaT, CellIds ids) const override
    {
        double kappa = mEvolution.Kappa(strain, deltaT, ids);
        double omega = mDamageLaw.Damage(kappa);
        double dOmegadKappa = mDamageLaw.Derivative(kappa);

        return mElasticDamage.DstressDstrain(strain, omega) +
               mElasticDamage.DstressDomega(strain, omega) * dOmegadKappa *
                       mEvolution.DkappaDstrain(strain, deltaT, ids);
    }

    void Update(EngineeringStrain<TDim> strain, double deltaT, CellIds ids)
    {
        mEvolution.Update(strain, deltaT, ids);
    }

public:
    // Intentionally made public. I assume that we know what we do. And this avoids lots of code, either getters
    // (possibly nonconst) or forwarding functions like `double Kappa(...), double Damage(...)`
    LinearElasticDamage<TDim> mElasticDamage;
    TDamageLaw mDamageLaw;
    TEvolution mEvolution;
};

template <int TDim>
class EvolutionImplicit
{

public:
    EvolutionImplicit(Constitutive::ModifiedMisesStrainNorm<TDim> strainNorm, size_t numCells = 1,
                      size_t numIpsPerCell = 1)
        : mStrainNorm(strainNorm)
    {
        ResizeHistoryData(numCells, numIpsPerCell);
    }

    EvolutionImplicit(Material::Softening m, size_t numCells = 1, size_t numIpsPerCell = 1)
        : mStrainNorm(m)
    {
        ResizeHistoryData(numCells, numIpsPerCell);
    }


    double Kappa(EngineeringStrain<TDim> strain, double, CellIds ids) const
    {
        return std::max(mStrainNorm.Value(strain), mKappas(ids.cellId, ids.ipId));
    }

    Eigen::Matrix<double, 1, Voigt::Dim(TDim)> DkappaDstrain(EngineeringStrain<TDim> strain, double, CellIds ids) const
    {
        if (mStrainNorm.Value(strain) >= mKappas(ids.cellId, ids.ipId))
            return mStrainNorm.Derivative(strain).transpose();
        return Eigen::Matrix<double, 1, Voigt::Dim(TDim)>::Zero();
    }

    void Update(EngineeringStrain<TDim> strain, double deltaT, CellIds ids)
    {
        mKappas(ids.cellId, ids.ipId) = Kappa(strain, deltaT, ids);
    }

    void ResizeHistoryData(size_t numCells, size_t numIpsPerCell)
    {
        mKappas.setZero(numCells, numIpsPerCell);
    }

public:
    Constitutive::ModifiedMisesStrainNorm<TDim> mStrainNorm;
    Eigen::MatrixXd mKappas;
};

} /* Laws */
} /* NuTo */