MarmotMechanicsCore

template<size_t materialTangentSize, size_t nIntegrationDependentStateVars>
class AdaptiveSubstepper

Public Types

typedef Eigen::Matrix<double, materialTangentSize, materialTangentSize> TangentSizedMatrix

Matrix for describing the nonlinear equation system of the return mapping algorithm.

Adaptive Substepper, employing an error estimation and Richardson Extrapolation for an Implicit Return Mapping algorithm

Matthias Neuner (2016), developed within the DK-CIM collaboration

typedef Eigen::Matrix<double, nIntegrationDependentStateVars, 1> IntegrationStateVector

Vector to carry the internal state of a material.

Public Functions

AdaptiveSubstepper(double initialStepSize, double minimumStepSize, double maxScaleUpFactor, double scaleDownFactor, double integrationErrorTolerance, int nPassesToIncrease, const Matrix6d &Cel)

void setConvergedProgress(const Marmot::Vector6d &stressOld, const IntegrationStateVector &stateVarsOld)

Set the current converged state.

bool isFinished()

Check if the subincrementation is finished.

double getNextSubstep()

get the next subincrement size

int getNumberOfSubsteps()

get the finished number of subincrements

bool finishSubstep(const Marmot::Vector6d &resultStress, const TangentSizedMatrix &dXdY, const IntegrationStateVector &stateVars)

finish a subincrement with plasticity

void finishElasticSubstep(const Marmot::Vector6d &resultStress)

finish a subincrement with elasticity only

bool discardSubstep()

discard the current subincrement

bool repeatSubstep(double decrementationFactor)

repeat the current subincrement with a smaller step

void getConvergedProgress(Marmot::Vector6d &stress, IntegrationStateVector &stateVars)

get the last converged state

Matrix6d getCurrentTangentOperator()

get the algorithmic tangent for the current state

void getResults(Marmot::Vector6d &stress, Matrix6d &consistentTangent, IntegrationStateVector &stateVars)

Write the current results.

Private Types

enum SubsteppingState

Values:

enumerator FullStep

enumerator FirstHalfStep

enumerator SecondHalfStep

Private Functions

bool acceptSubstepWithFullStepOnly()

bool splitCurrentSubstep()

Private Members

const double initialStepSize

const double minimumStepSize

const double maxScaleUpFactor

const double scaleDownFactor

const double integrationErrorTolerance

const int nPassesToIncrease

const bool ignoreErrorToleranceOnMinimumStepSize

double currentProgress

double currentSubstepSize

int passedSubsteps

int substepIndex

Marmot::Vector6d stressProgress

internal storages for the progress of the total increment

IntegrationStateVector stateProgress

TangentSizedMatrix consistentTangentProgress

Marmot::Vector6d stressProgressHalfTemp

Marmot::Vector6d stressProgressFullTemp

IntegrationStateVector stateProgressHalfTemp

IntegrationStateVector stateProgressFullTemp

TangentSizedMatrix consistentTangentProgressHalfTemp

TangentSizedMatrix consistentTangentProgressFullTemp

TangentSizedMatrix elasticTangent

SubsteppingState currentState

template<size_t materialTangentSize, size_t nIntegrationDependentStateVars>
class AdaptiveSubstepperExplicit

Public Types

typedef Eigen::Matrix<double, materialTangentSize, materialTangentSize> TangentSizedMatrix

Matrix for describing the nonlinear equation system of the return mapping algorithm.

typedef Eigen::Matrix<double, materialTangentSize, 6> MatrixStateStrain

typedef Eigen::Matrix<double, nIntegrationDependentStateVars, 1> IntegrationStateVector

Public Functions

AdaptiveSubstepperExplicit(double initialStepSize, double minimumStepSize, double maxScaleUpFactor, double scaleDownFactor, double integrationErrorTolerance, int nPassesToIncrease, const Matrix6d &Cel)

void setConvergedProgress(const Marmot::Vector6d &stressOld, const IntegrationStateVector &stateVarsOld)

Set the current converged state.

bool isFinished()

Check if the subincrementation is finished.

double getNextSubstep()

get the next subincrement size

int getNumberOfSubsteps()

get the finished number of subincrements

int getNumberDiscardedSubstepsDueToError()

get the number of discarded subincrements

bool finishSubstep(const Marmot::Vector6d &resultStress, const TangentSizedMatrix &dXdY, const TangentSizedMatrix &dYdXOld, const IntegrationStateVector &stateVars)

finish a subincrement with plasticity

void finishElasticSubstep(const Marmot::Vector6d &resultStress)

finish a subincrement with elasticity only

bool discardSubstep()

discard the current subincrement

bool repeatSubstep(double decrementationFactor)

repeat the current subincrement with a smaller step

void getConvergedProgress(Marmot::Vector6d &stress, IntegrationStateVector &stateVars)

get the last converged state

Matrix6d getCurrentTangentOperator()

get the algorithmic tangent for the current state

void getResults(Marmot::Vector6d &stress, Matrix6d &consistentTangent, IntegrationStateVector &stateVars)

Write the current results.

Private Types

enum SubsteppingState

Values:

enumerator FullStep

enumerator FirstHalfStep

enumerator SecondHalfStep

Private Functions

const TangentSizedMatrix &IXX()

const MatrixStateStrain &IX6()

bool acceptSubstepWithFullStepOnly()

bool splitCurrentSubstep()

Private Members

const double initialStepSize

const double minimumStepSize

const double maxScaleUpFactor

const double scaleDownFactor

const double integrationErrorTolerance

const int nPassesToIncrease

const bool ignoreErrorToleranceOnMinimumStepSize

const Matrix6d &Cel

double currentProgress

double currentSubstepSize

int passedSubsteps

int substepIndex

int discardedDueToError

Marmot::Vector6d stressProgress

internal storages for the progress of the total increment

IntegrationStateVector stateProgress

internal storages for the progress of the total increment

MatrixStateStrain consistentTangentProgress

internal storages for the progress of the total increment

Marmot::Vector6d stressProgressHalfTemp

temporal storages, which are used until a cycle full/half/half has finished successfully

Marmot::Vector6d stressProgressFullTemp

IntegrationStateVector stateProgressHalfTemp

temporal storages, which are used until a cycle full/half/half has finished successfully

IntegrationStateVector stateProgressFullTemp

MatrixStateStrain consistentTangentProgressHalfTemp

temporal storages, which are used until a cycle full/half/half has finished successfully

MatrixStateStrain consistentTangentProgressFullTemp

SubsteppingState currentState

template<int nMatTangentSize>
class DuvautLionsViscosity

#include <DuvautLionsViscosity.h>

Implementation of Duvaut-Lions viscosity for a material with nMatTangentSize internal degrees of freedom.

Todo:

: Update member names to more descriptive ones

Public Types

typedef Eigen::Matrix<double, nMatTangentSize, nMatTangentSize> TangentSizedMatrix

Public Functions

DuvautLionsViscosity(double viscosity)

double applyViscosityOnStateVar(double stateVarTrial, double StateVarInf, double dT)

Apply viscosity on a scalar internal variable depending on the extremal solutions for t=0 (trialState) and t= \(\infty\), and timestep dt

Marmot::Vector6d applyViscosityOnStress(const Marmot::Vector6d &trialStress, const Marmot::Vector6d &stressInf, double dT)

Apply viscosity on voigt sized rank two tensor depending on the extremal solutions for t=0 (trialState) and t= \(\infty\), and timestep dt

TangentSizedMatrix applyViscosityOnMatTangent(const TangentSizedMatrix &matTangentInv, double dT)

Apply viscosity on the inverse (algorithmic) material tangent depending on the extremal solutions for t=0 (trialState) and t= \(\infty\), and timestep dt

Todo:

: Check if application to inverse can be replaced by application to non-inverse in general

Private Members

const double viscosity

template<int N>
struct Factorial

#include <MarmotKelvinChain.h>

Compile-time factorial.

Recursive template structure computing the factorial of N at compile time.

Example:

int f = Factorial<5>::value; // f = 120

Template Parameters:

N – Non-negative integer whose factorial is to be computed.

Public Types

Values:

enumerator value

template<>
struct Factorial<0>

#include <MarmotKelvinChain.h>

Template specialization of Factorial for 0.

Defines the base case of the factorial recursion, setting 0! = 1.

Public Types

Values:

enumerator value

Values:

template<typename T = double>
struct HaighWestergaardCoordinates

#include <HaighWestergaard.h>

Aggregate of the Haigh-Westergaard coordinates (invariants) \(\xi\), \(\rho\), \(\theta\).

Public Members

T xi

Hydrostatic component \(\xi\).

T rho

Deviatoric radius \(\rho\).

T theta

Lode angle \(\theta\) specified in radian.

struct HistoryEntry

#include <MarmotMaterialPointSolverHypoElastic.h>

Struct to record the history of the simulation.

Each entry contains time, stress, strain, and state variables.

Public Functions

inline void print() const

Public Members

double time

Time at the history entry.

Marmot::Vector6d stress

Stress at the history entry.

Marmot::Vector6d strain

Strain at the history entry.

Marmot::Matrix6d dStressdStrain

Material tangent at the history entry.

Eigen::VectorXd stateVars

State variables at the history entry.

class HughesWinget

Public Types

enum Formulation

Values:

enumerator AbaqusLike

Public Functions

HughesWinget(const Eigen::Matrix3d &FOld, const Eigen::Matrix3d &FNew, Formulation formulation)

Marmot::Vector6d getStrainIncrement()

Eigen::Matrix3d getRotationIncrement()

Marmot::Vector6d rotateTensor(const Marmot::Vector6d &tensor)

Marmot::EigenTensors::Tensor633d compute_dS_dF(const Marmot::Vector6d &stress, const Eigen::Matrix3d &FInv, const Marmot::Matrix6d &dChauchyDEps)

Eigen::Matrix3d compute_dScalar_dF(const Eigen::Matrix3d &FInv, const Marmot::Vector6d &dScalarDEps)

Private Members

Formulation theFormulation

Eigen::Matrix3d l

Eigen::Matrix3d dOmega

Eigen::Matrix3d dR

Marmot::Vector6d dEps

struct Increment

#include <MarmotMaterialPointSolverHypoElastic.h>

Struct to define a loading increment.

Each increment contains strain and stress increments, control flags, time information, and iteration limits.

Public Members

Marmot::Vector6d strainIncrement

Strain increment for the increment.

Marmot::Vector6d stressIncrement

Stress increment for the increment.

Eigen::Vector<bool, 6> isStrainComponentControlled

Flags to indicate which strain components are controlled.

Eigen::Vector<bool, 6> isStressComponentControlled

Flags to indicate which stress components are controlled.

double timeOld

Old time at the beginning of the increment.

double dT

Time step size for the increment.

class MarmotMaterialGradientEnhancedHypoElastic : public MarmotMaterialGradientEnhancedMechanical

#include <MarmotMaterialGradientEnhancedHypoElastic.h>

Inheritance diagram for MarmotMaterialGradientEnhancedHypoElastic:

Collaboration diagram for MarmotMaterialGradientEnhancedHypoElastic:

Derived abstract base class for gradient-enhanced hypo-elastic materials expressed purely in rate form.

In general, the nominal stress rate tensor \( \sigRate \) can be written as a function of the nominal stress tensor \( \sig \), the stretching rate tensor \( \epsRate \) the local variable \(\kappa\), the nonlocal variable \(\bar{\kappa}\), internal variables \(\boldsymbol{x}\) and the time \( t \)

\[ \displaystyle \sigRate = f( \sig, \kappa, \epsRate, \dot{\bar{\kappa}}, \boldsymbol{x}, t, ...) \]

In course of numerical time integration, this relation will be formulated incrementally as

\[ \displaystyle \Delta \sig = f ( \sig_n, \kappa_n, \Delta\eps, \Delta\bar{\kappa}, \boldsymbol{x}_n, \Delta t, t_n, ...) \]

with

\[ \displaystyle \Delta\eps = \epsRate\, \Delta t \]

\[ \displaystyle \Delta\bar{\kappa} = \dot{\bar{\kappa}}\, \Delta t \]

and the algorithmic tangents

\[ \displaystyle \frac{d \sig }{d \eps } = \frac{d \Delta \sig }{d \Delta \eps } \]

\[ \displaystyle \frac{d \sig }{d \bar{\kappa}} = \frac{d \Delta \sig }{d \Delta \bar{\kappa}} \]

\[ \displaystyle \frac{d \kappa }{d \eps }= \frac{d \Delta \kappa }{d \Delta \eps } \]

\[ \displaystyle \frac{d \kappa }{d \bar{\kappa}}= \frac{d \Delta \kappa }{d \Delta \bar{\kappa}} = l_{\text{nonlocal}} \]

Public Functions

virtual void computeStress(double *stress, double &KLocal, double &nonLocalRadius, double *dStress_dDeformationGradient, double *dKLocal_dDeformationGradient, double *dStress_dK, const double *FOld, const double *FNew, const double KOld, const double dK, const double *timeOld, const double dT, double &pNewDT) override

For a given deformation gradient and a nonlocal variable at the old and the current time, compute the Cauchy stress and the local variable and the algorithmic tangents.

Parameters:

• stress – [inout] Cauchy stress

• K_local – [inout] local variable

• nonLocalRadius – [inout] nonlocal radius representing the tangent \(\frac{d \Delta \kappa }{d \Delta \bar{\kappa}}\)

• dStressDDDeformationGradient – [inout] Derivative of the Cauchy stress tensor with respect to the deformation gradient \(\boldsymbol{F}\)

• dK_localDDeformationGradient – [inout] Derivative of the local variable with respect to the deformation gradient \(\boldsymbol{F}\)

• dStressDK – [inout] Derivative of the Cauchy stress tensor with respect to the nonlocal variable

• FOld – [in] Deformation gradient at the old (pseudo-)time

• FNew – [in] Deformation gradient at the current (pseudo-)time

• KOld – [in] Local variable at the old (pseudo-)time

• dK – [in] Increment of the local variable

• timeOld – [in] Old (pseudo-)time

• dt – [in] (Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

• pNewDT – [inout] Suggestion for a new time increment

virtual void computeStress(double *stress, double &K_local, double &nonLocalRadius, double *dStressDDStrain, double *dK_localDDStrain, double *dStressDK, const double *dStrain, double KOld, double dK, const double *timeOld, const double dT, double &pNewDT) = 0

For a given deformation gradient and a nonlocal variable at the old and the current time, compute the Cauchy stress and the local variable and the algorithmic tangents.

Parameters:

• stress – [inout] Cauchy stress

• K_local – [inout] local variable

• nonLocalRadius – [inout] nonlocal radius representing the tangent \(\frac{d \Delta \kappa }{d \Delta \bar{\kappa}}\)

• dStressDDstrain – [inout] Algorithmic tangent representing the derivative of the Cauchy stress tensor with respect to the linearized strain

• dK_localDDStrain – [inout] Derivative of the local variable with respect to the linearized strain

• dStressDK – [inout] Derivative of the Cauchy stress tensor with respect to the nonlocal variable

• KOld – [in] Local variable at the old (pseudo-)time

• dK – [in] Increment of the local variable

• timeOld – [in] Old (pseudo-)time

• dt – [in] (Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

• pNewDT – [inout] Suggestion for a new time increment

virtual void computePlaneStress(double *stress2D, double &KLocal2D, double &nonLocalRadius, double *dStress_DStrain2D, double *dKLocal_dStrain2D, double *dStress_dK2D, const double *dStrain2D, double KOld, double dK, const double *timeOld, const double dT, double &pNewDT)

Plane stress implementation of computeStress.

inline virtual void computePlaneStress(double *stress, double &K_local, double &nonLocalRadius, double *dStressDDFNew, double *dK_localDDFNew, double *dStressDK, const double *FOld, const double *FNew, double KOld, double dK, const double *timeOld, const double dT, double &pNewDT)

Plane stress implementation of computeStress.

Todo:

why is dStrain an in-out parameter?

class MarmotMaterialGradientEnhancedMechanical : public MarmotMaterial

#include <MarmotMaterialGradientEnhancedMechanical.h>

Inheritance diagram for MarmotMaterialGradientEnhancedMechanical:

Collaboration diagram for MarmotMaterialGradientEnhancedMechanical:

Abstract basic class for Mechanical materials with scalar nonlocal interaction.

Formulated incrementally as σ_np, K_local_np = f (σ_n, dxdX_n, dxdX_np, Δt, t_n, , Kn, ΔK, K_local_n.. )

Algorithmic tangents: dσdF = d σ_np d (dxdX_np) dK_LocaldF = d K_local_np d (dxdX_np) dσdK = d σ_np d ΔK

Subclassed by MarmotMaterialGradientEnhancedHypoElastic

Public Functions

virtual void computeStress(double *stress, double &K_local, double &nonLocalRadius, double *dStressDDFNew, double *dK_localDDFNew, double *dStressDK, const double *FOld, const double *FNew, const double KOld, const double dK, const double *timeOld, const double dT, double &pNewDT) = 0

inline virtual void computePlaneStress(double *stress, double &K_local, double &nonLocalRadius, double *dStressDDFNew, double *dK_localDDFNew, double *dStressDK, const double *FOld, const double *FNew, double KOld, double dK, const double *timeOld, const double dT, double &pNewDT)

inline virtual void computeUniaxialStress(double *stress, double &K_local, double &nonLocalRadius, double *dStressDDFNew, double *dK_localDDFNew, double *dStressDK, const double *FOld, const double *FNew, double KOld, double dK, const double *timeOld, const double dT, double &pNewDT)

class MarmotMaterialHyperElastic : public MarmotMaterialMechanical

#include <MarmotMaterialHyperElastic.h>

Inheritance diagram for MarmotMaterialHyperElastic:

Collaboration diagram for MarmotMaterialHyperElastic:

Derived abstract base class for simple, purely hyperelastic materials to be used for finite elements based on the total lagrangian kinematic description (TL elements). The second Piola - Kirchhoff stress tensor \( S \) will be derived by

\[ \displaystyle S = \frac{\partial f(\boldsymbol{E},t )}{\partial \boldsymbol{E}} \]

with the Green - Lagrange strain tensor \( \boldsymbol{E} \)

\[ \displaystyle E = \frac{1}{2}\,\left(\boldsymbol{F}^T\cdot \boldsymbol{F} - \boldsymbol{I} \right) \]

as work conjugated measure and the variable \( \boldsymbol{F} \) denoting the deformation gradient. The algorithmic tangent will be calculated by

\[ \displaystyle \frac{d \boldsymbol{S}}{d \boldsymbol{E}} \]

Public Functions

virtual void computeStress(double *S, double *dSdE, const double *FOld, const double *FNew, const double *timeOld, const double dT, double &pNewDT) override

For a given deformation gradient at the old and the current time, compute the 2nd Piola-Kirchhoff stress and the algorithmic tangent \(\frac{\partial\boldsymbol{S}^{(n+1)}}{\partial\boldsymbol{E}^{(n+1)}}\).

Todo:

A default implementation is provided.

Parameters:

• S – [inout] 2nd Piola-Kirchhoff stress

• dSdE – [inout] Algorithmic tangent representing the derivative of the 2nd Piola-Kirchhoff stress tensor with respect to the Green-Lagrange strain tensor \(\boldsymbol{E}\).

• FOld – [in] Deformation gradient at the old (pseudo-)time

• FNew – [in] Deformation gradient at the current (pseudo-)time

• timeOld – [in] Old (pseudo-)time

• dt – [in] (Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

• pNewDT – [inout] Suggestion for a new time increment

virtual void computeStressPK2(double *S, double *dSdE, const double *E, const double *timeOld, const double dT, double &pNewDT) = 0

For a given Green-Lagrange strain, compute the 2nd Piola-Kirchhoff stress and the algorithmic tangent \(\frac{\partial\boldsymbol{S}^{(n+1)}}{\partial\boldsymbol{E}^{(n+1)}}\).

Todo:

Should we use function overloading in this case and simple also use computeStress for the function name?

Parameters:

• S – [inout] 2nd Piola-Kirchhoff stress

• dSdE – [inout] Algorithmic tangent representing the derivative of the 2nd Piola-Kirchhoff stress tensor with respect to the Green-Lagrange strain tensor \(\boldsymbol{E}\).

• deltaE – [in] Green-Lagrange strain increment

• timeOld – [in] Old (pseudo-)time

• dt – [in] (Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

• pNewDT – [inout] Suggestion for a new time increment

virtual void computePlaneStressPK2(double *S2D, double *dSdE2D, const double *E2D, const double *timeOld, const double dT, double &pNewDT)

Plane stress implementation of computeStressPK2.

virtual void computeUniaxialStressPK2(double *S1D, double *dSdE1D, const double *E1D, const double *timeOld, const double dT, double &pNewDT)

Uniaxial stress implementation of computeStressPK2.

class MarmotMaterialHypoElastic : public MarmotMaterialMechanical

#include <MarmotMaterialHypoElastic.h>

Inheritance diagram for MarmotMaterialHypoElastic:

Collaboration diagram for MarmotMaterialHypoElastic:

Derived abstract base class for elastic materials expressed purely in rate form.

In general, the nominal stress rate tensor \( \sigRate \) can be written as a function of the nominal stress tensor \( \sig \), the stretching rate tensor \( \epsRate \) and the time \( t \).

\[ \displaystyle \sigRate = f( \sig, \epsRate, t, ...) \]

In course of numerical time integration, this relation will be formulated incrementally as

\[ \displaystyle \Delta \sig = f ( \sig_n, \Delta\eps, \Delta t, t_n, ...) \]

with

\[ \displaystyle \Delta\eps = \epsRate\, \Delta t \]

and the algorithmic tangent

\[ \displaystyle \frac{d \sig }{d \eps } = \frac{d \Delta \sig }{d \Delta \eps } \]

This formulation is compatible with an Abaqus interface.

Subclassed by MarmotMaterialHypoElasticAD

Public Functions

void setCharacteristicElementLength(double length)

Set the characteristic element length at the considered quadrature point. It is needed for the regularization of materials with softening behavior based on the mesh-adjusted softening modulus.

Parameters:

length – [in] characteristic length; will be assigned to characteristicElementLength

virtual void computeStress(double *stress, double *dStressDDStrain, const double *FOld, const double *FNew, const double *timeOld, const double dT, double &pNewDT) override

For a given deformation gradient at the old and the current time, compute the Cauchy stress and the algorithmic tangent \(\frac{\partial\boldsymbol{\sigma}^{(n+1)}}{\partial\boldsymbol{F}^{(n+1)}}\).

Todo:

A default implementation is provided.

Parameters:

• stress – [inout] Cauchy stress

• dSdE – [inout] Algorithmic tangent representing the derivative of the Cauchy stress tensor with respect to the deformation gradient \(\boldsymbol{F}\)

• FOld – [in] Deformation gradient at the old (pseudo-)time

• FNew – [in] Deformation gradient at the current (pseudo-)time

• timeOld – [in] Old (pseudo-)time

• dt – [in] (Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

• pNewDT – [inout] Suggestion for a new time increment

virtual void computeStress(double *stress, double *dStressDDStrain, const double *dStrain, const double *timeOld, const double dT, double &pNewDT) = 0

For a given linearized strain increment \(\Delta\boldsymbol{\varepsilon}\) at the old and the current time, compute the Cauchy stress and the algorithmic tangent \(\frac{\partial\boldsymbol{\sigma}^{(n+1)}}{\partial\boldsymbol{\varepsilon}^{(n+1)}}\).

Parameters:

• stress – [inout] Cauchy stress

• dStressDDstrain – [inout] Algorithmic tangent representing the derivative of the Cauchy stress tensor with respect to the linearized strain

• dStrain – [in] linearized strain increment

• timeOld – [in] Old (pseudo-)time

• dt – [in] (Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

• pNewDT – [inout] Suggestion for a new time increment

virtual void computePlaneStress(double *stress2D, double *dStress_dStrain2D, const double *dStrain2D, const double *timeOld, const double dT, double &pNewDT)

virtual void computeUniaxialStress(double *stress1D, double *dStress_dStrain1D, const double *dStrain, const double *timeOld, const double dT, double &pNewDT)

virtual void computePlaneStress(double *stress2D, double *dStress_dF2DNew, const double *FOld2D, const double *FNew2D, const double *timeOld, const double dT, double &pNewDT)

Plane stress implementation of computeStress.

virtual void computeUniaxialStress(double *stress1D, double *dStress1D_dF1DNew, const double *F1DOld, const double *F1DNew, const double *timeOld, const double dT, double &pNewDT)

Uniaxial stress implementation of computeStress.

Public Members

double characteristicElementLength

Characteristic element length.

class MarmotMaterialHypoElasticAD : public MarmotMaterialHypoElastic

Inheritance diagram for MarmotMaterialHypoElasticAD:

Collaboration diagram for MarmotMaterialHypoElasticAD:

Public Functions

virtual void computeStressAD(autodiff::dual *stress, const autodiff::dual *dStrain, const double *timeOld, const double dT, double &pNewDT) = 0

Compute the Cauchy stress tensor \(\boldsymbol{\sigma}\) given an increment of the linearized strain tensor \(\Delta\boldsymbol{\varepsilon}\).

Dual numbers are used for both the Cauchy stress tensor and the linearized strain increment, such that the algorithmic tangent operator \(\frac{\partial\boldsymbol{\sigma}^{(n+1)}}{\partial\boldsymbol{\varepsilon}^{(n+1)}}\) can be obtained by means of automatic differentiation.

Parameters:

• stress – [inout] Cauchy stress tensor

• dStrain – [in] linearized strain increment

• timeOld – [in] Old (pseudo-)time

• dT – [in] (Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

• pNewDT – [inout] Suggestion for a new time increment

virtual void computeStress(double *stress, double *dStressDDStrain, const double *dStrain, const double *timeOld, const double dT, double &pNewDT) override

For a given linearized strain increment \(\Delta\boldsymbol{\varepsilon}\) at the old and the current time, compute the Cauchy stress and the algorithmic tangent \(\frac{\partial\boldsymbol{\sigma}^{(n+1)}}{\partial\boldsymbol{\varepsilon}^{(n+1)}}\).

Parameters:

• stress – [inout] Cauchy stress

• dStressDDstrain – [inout] Algorithmic tangent representing the derivative of the Cauchy stress tensor with respect to the linearized strain

• dStrain – [in] linearized strain increment

• timeOld – [in] Old (pseudo-)time

• dt – [in] (Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

• pNewDT – [inout] Suggestion for a new time increment

class MarmotMaterialMechanical : public MarmotMaterial

#include <MarmotMaterialMechanical.h>

Inheritance diagram for MarmotMaterialMechanical:

Collaboration diagram for MarmotMaterialMechanical:

Abstract basic class for Mechanical materials. ‘Mechanical’ is meant in the ‘most general sense’, i.e., any material which describes a mechanical (cauchy) stress - deformation relationship (e.g, hyperelastic, hypoelastic, elasto-plastic, visco-elastic materials)

σ = f (σ, dxdX, t, .. ),

formulated incrementally as σ_np = f (σ_n, dxdX_n, dxdX_np, Δt, t_n, .. )

Algorithmic tangent: dσdF = d σ_np d (dxdX_np)

Format:

/ d σ_11 d F_00, d σ_11 d F_10, d σ_11 d F_20, d σ_11 d F_01, \ | d σ_22 d F_00, | | d σ_33 d F_00, | | … | | … | \ … /

such that it can be interpreted as a column major 6x3x3 tensor (4th order, left voigt tensor)

Subclassed by MarmotMaterialHyperElastic, MarmotMaterialHypoElastic

Public Functions

virtual void computeStress(double *stress, double *dStress_dFNew, const double *FOld, const double *FNew, const double *timeOld, const double dT, double &pNewDT) = 0

virtual void computePlaneStress(double *stress2D, double *dStress_dF2DNew, const double *FOld2D, const double *FNew2D, const double *timeOld, const double dT, double &pNewDT)

virtual void computeUniaxialStress(double *stress1D, double *dStress1D_dF1DNew, const double *F1DOld, const double *F1DNew, const double *timeOld, const double dT, double &pNewDT)

class MarmotMaterialPointSolverHypoElastic

#include <MarmotMaterialPointSolverHypoElastic.h>

Solver for material point problems with hypo-elastic materials.

This class implements a solver for material point problems using hypo-elastic material models. It supports loading steps with controlled strain and stress components, adaptive time stepping, and history recording.

Public Functions

MarmotMaterialPointSolverHypoElastic(std::string &materialName, double *materialProperties, int nMaterialProperties, const SolverOptions &options)

Constructor for the MarmotMaterialPointSolverHypoElastic class.

Parameters:

• materialName – Name of the hypo-elastic material model

• materialProperties – Array of material properties

• nMaterialProperties – Number of material properties

void addStep(const Step &step)

Add a loading step to the solver.

Parameters:

step – The Step to be added

inline std::vector<Step> getSteps() const

Get the list of added loading steps.

Returns:

A vector of Step containing the added steps

inline void clearSteps()

Clear all added loading steps.

void setInitialState(const Marmot::Vector6d &initialStress, const Eigen::VectorXd &initialStateVars)

Set the initial state of the material model.

Parameters:

• initialStress – The initial stress in Voigt notation

• initialStateVars – The initial state variables

inline int getNumberOfStateVariables(int &nStateVarsOut) const

Get the number of state variables in the material model.

Returns:

The number of state variables

void resetToInitialState()

Reset the solver to the initial state.

This function resets the initial stress and state variables of the material model.

void solve()

Solve the material point problem for all added steps.

inline std::vector<HistoryEntry> getHistory() const

Get the recorded history of the simulation.

Returns:

A vector of HistoryEntry containing the recorded history

inline void clearHistory()

Clear the recorded history.

void printHistory()

Print the recorded history to the console.

void exportHistoryToCSV(const std::string &filename)

Export the recorded history to a CSV file.

Parameters:

filename – The name of the CSV file to export to

Private Functions

void solveStep(const Step &step)

Solve a single loading step.

This function iterates over the increments defined in the step, calling solveIncrement for each increment. It manages time stepping and ensures that the entire step is covered.

Parameters:

step – The Step to be solved

void solveIncrement(const Increment &increment)

Solve a single increment within a loading step.

This function implements a Newton-Raphson iterative solver to compute the stress and strain state for the given increment. It updates the material state variables and records the history after convergence.

Parameters:

increment – The Increment to be solved

Throws:

std::runtime_error – if the solver does not converge

Marmot::Vector6d computeResidual(const Marmot::Vector6d &stressIncrement, const Marmot::Vector6d &target, const Increment &increment)

Compute the residual for the current increment.

This function calculates the residual vector based on the difference between the computed stress increment and the target increment, taking into account which components are controlled by strain or stress.

Parameters:

• stressIncrement – The computed stress increment

• target – The target (mixed) stress/strain increment

• increment – The Increment containing control information

Returns:

The computed residual vector

void modifyTangent(Eigen::Matrix<double, 6, 6> &tangent, const Increment &increment)

Modify the material tangent matrix based on control type.

This function adjusts the tangent matrix to account for the components that are controlled by strain or stress, ensuring that the solver correctly handles mixed control scenarios. This is done by zeroing out rows corresponding to strain-controlled components and setting their diagonal entries to one.

Parameters:

• tangent – The material tangent matrix to be modified

• increment – The Increment containing control information

Private Members

MarmotMaterialHypoElastic *material

The hypo-elastic material model.

int nStateVars

Number of state variables in the material model.

Eigen::VectorXd stateVars

Current state variables.

Eigen::VectorXd _initialStateVars

Initial state variables.

Eigen::VectorXd stateVarsTemp

Temporary state variables for computations.

Marmot::Vector6d stress = Marmot::Vector6d::Zero()

The stress in Voigt notation.

Marmot::Vector6d _initialStress = Marmot::Vector6d::Zero()

The initial stress in Voigt notation.

Marmot::Vector6d strain = Marmot::Vector6d::Zero()

The strain in Voigt notation.

Marmot::Matrix6d dStressDStrain = Marmot::Matrix6d::Zero()

The material tangent matrix in Voigt notation.

std::vector<Step> steps

List of loading steps.

std::vector<HistoryEntry> history

History of the simulation.

const SolverOptions options

Solver options.

class MarmotStateVarVectorManager

#include <MarmotStateVarVectorManager.h>

Collaboration diagram for MarmotStateVarVectorManager:

A convenience auxiliary class for managing multiple statevars with arbitrary length in a single consecutive double array.

Public Functions

inline StateView getStateView(const std::string &name) const

get a StateView for a statevar entry

inline double &find(const std::string &name) const

get the reference to the first array element of an entry in the statevar vector

inline bool contains(const std::string &name) const

check if the entry with name is managed

Protected Functions

inline MarmotStateVarVectorManager(double *theStateVars, const StateVarVectorLayout &theLayout_)

Protected Attributes

double *theStateVars

pointer to the first element in the statevar vector

const StateVarVectorLayout &theLayout

a const reference to the respective layout

Protected Static Functions

static inline StateVarVectorLayout makeLayout(const std::vector<StateVarEntryDefinition> &theEntries)

generate the statevar vector layout from a list of entries, defined by name and length

class MenetreyWillam

#include <MenetreyWillam.h>

Collaboration diagram for Marmot::ContinuumMechanics::CommonConstitutiveModels::MenetreyWillam:

Generalized failure criterion proposed by Menétrey and Willam.

Implementation of the generalized failure criterion proposed by Menétrey and Willam

\[ f(\xi,\rho,\theta) = \left(A_f\,\rho\right)^2+m\left(B_f\,\rho\,r(\theta,e)+C_f\,\xi\right) - 1 \]

with the Haigh-Westergaard coordinates \(\xi,\,\rho,\,\theta\) and the parameters \(A_f,\,B_f,\,C_f,\,m,\,e\) which are automatically evaluated dependent on the desired formulation (supported types see MenetreyWillamType).

Example

For a given pair of principal stress components in the Haigh-Westergaard stress space, evaluate the Mohr-Coulomb failure criterion:

#include <Marmot/MenetreyWillam.h>
#include <Marmot/HaighWestergaard.h>

using namespace Marmot;

void main()
{
  const double ft = 1.0;
  const double fc = 10.0;
  auto mw = MenetreyWillam( ft , MenetreyWillamType::MohrCoulomb, fc );
  const HaighWestergaard::HaighWestergaardCoordinates hw = { 1.0, -10, 0 };
  if ( mw.yieldCriterion( hw ) >= 0 )
    std::cout << "Material is yielding!" << std::endl;
  else
    std::cout << "Material is elastic!" << std::endl;
}

Public Types

enum class MenetreyWillamType

Values:

enumerator Mises

von-Mises failure criterion; only the tensile strength ft has to be specified.

enumerator Rankine

Rankine failure criterion; only the tensile strength ft has to be specified.

Note

In deviatoric sections, compared to the original Rankine criterion the vertices are slightly rounded (specified by the eccentricity parameter \(e\)).

enumerator DruckerPrager

Drucker-Prager failure criterion; both the tensile strength ft as well as the compressive strength fc have to be specified.

enumerator MohrCoulomb

Drucker-Prager failure criterion; both the tensile strength ft as well as the compressive strength fc have to be specified.

Note

In deviatoric sections, compared to the original Mohr-Coulomb criterion the vertices are slightly rounded (specified by the eccentricity parameter \(e\)).

Public Functions

MenetreyWillam(const double ft, const MenetreyWillamType &type = MenetreyWillamType::Rankine, const double fc = 0)

Constructor that takes the uniaxial tensile strength \(f_t\) and two optional arguments consisting of the specific type of failure criterion MenetreyWillamType and the uniaxial compressive strength \(f_c\). The call of the constructor automatically fills the corresponding Menetrey-Willam param.

void setParameters(const double ft, const double fc, const MenetreyWillamType &type)

This function can be used to reset the type of the specified failure criterion by entering the uniaxial compressive strength fc, the uniaxial tensile strength ft and the MenetreyWillamType.

template<typename T>
inline T polarRadius(const double &theta) const

Compute the polar radius \(r\) from the Lode angle \(\theta\). The eccentricity parameter will be used from the chosen Menetrey-Willam parameters param.

template<typename T>
inline std::pair<T, T> dPolarRadius_dTheta(const T &theta) const

Compute the polar radius \(r\) and its derivative \(\frac{dr}{d\theta}\) from the Lode angle \(\theta\). The eccentricity parameter will be used from the param struct.

template<typename T>
inline T yieldFunction(const ContinuumMechanics::HaighWestergaard::HaighWestergaardCoordinates<T> &hw, const double varEps = 0.0) const

Evaluate the yield function \(f\) depending on the Haigh-Westergaard stress coordinates hw. \(f<0\) means no yielding while \(f\geq0\) means yielding. When the optional fillet parameter varEps is given, a potential vertex along the hydrostatic axis is rounded and thus a smooth failure criterion is obtained.

Note

The yield function can be also used as plastic potential function if needed.

template<typename T>
inline std::tuple<T, T, T> dYieldFunction_dHaighWestergaard(const ContinuumMechanics::HaighWestergaard::HaighWestergaardCoordinates<T> &hw, const double varEps = 0.0) const

Evaluate the derivatives of the yield function with respect to the Haigh-Westergaard stress coordinates, i.e. \(df/d\xi,\,df/d\rho,\,df/d\theta\). When the optional fillet parameter varEps is given, a potential vertex along the hydrostatic axis is rounded. Thus, a smooth failure criterion is obtained, where derivatives can be calculated uniquely.

Public Members

struct Marmot::ContinuumMechanics::CommonConstitutiveModels::MenetreyWillam::MenetreyWillamParameters param

Public Static Functions

template<typename T>
static inline T polarRadius(const T &theta, const double &e)

Static version for computing the polar radius \(r\) from the Lode angle \(\theta\) and a specified value for the eccentricity parameter \(e\).

template<typename T>
static inline std::pair<T, T> dPolarRadius_dTheta(const T &theta, const double &e)

Static version for computing the polar radius \(r$\) and its derivative \(\frac{dr}{d\theta}\) from the Lode angle \(\theta\) and a specified value for the eccentricity parameter \(e\).

template<typename T>
static inline std::tuple<T, T, T> d2PolarRadius_dTheta2(const T &theta, const double &e)

static inline double abaqusMohrCoulombPotentialVarEpsToMenetreyWillam(const double varEps, const double psi)

Compute a fillet parameter for the vertex of the yield surface along the hydrostatic axis in the same way as Abaqus does. This parameter is only relevant in the case of the Drucker-Prager or the Mohr-Coulomb criterion. The calculated smoothing value can then be used as the optional input argument varEps in yieldFunction and dYieldFunction_dHaighWestergaard.

static inline double e(const double fc, const double ft)

Compute the eccentricity parameter based on a given compressive strength fc and tensile strength ft.

static inline double c(const double fc, const double ft)

Compute the cohesion based on a given compressive strength fc and tensile strength ft.

static inline double phi(const double fc, const double ft)

Compute the friction angle in radian based on a given compressive strength fc and tensile strength ft.

static inline double ft(const double c, const double phi)

Compute the tensile strength based on a given cohesion c and a friction angle phi.

static inline double fc(const double c, const double phi)

Compute the compressive strength based on a given cohesion c and a friction angle phi.

struct MenetreyWillamParameters

#include <MenetreyWillam.h>

Aggregate of the parameters for the generalized strength criterion proposed by Menétrey and Willam.

Public Members

double Af

additional parameter to capture the influence of the deviatoric stress invariant \(\rho\)

double Bf

addditional parameter to capture the influence of the deviatoric stress invariant \(\rho\) and the lode angle

double Cf

addditional parameter to capture the influence of the hydrostatic stress invariant \(\xi\)

double m

Friction parameter \(m\)

double e

Eccentricity parameter \(e\); to obtain a smooth and convex surface $e$ has to be in the range of \(0.5\leq e \leq 1\)

template<int nSizeMatTangent>
class PerezFougetSubstepper

#include <PerezFougetSubstepperExplicitMarkII.h>

Substepper for (linear elastic) elastoplastic materials, semi-explicit return mapping version

Substepper for (linear elastic) elastoplastic materials, implicit return mapping version

Public Types

typedef Eigen::Matrix<double, nSizeMatTangent, nSizeMatTangent> TangentSizedMatrix

typedef Eigen::Matrix<double, nSizeMatTangent, 6> MatrixStateStrain

typedef Eigen::Matrix<double, nSizeMatTangent, nSizeMatTangent> TangentSizedMatrix

Public Functions

PerezFougetSubstepper(double initialStepSize, double minimumStepSize, double scaleUpFactor, double scaleDownFactor, int nPassesToIncrease, const Matrix6d &Cel)

bool isFinished()

check if subincrementation process has finished

double getNextSubstep()

get the next subincrement size

bool decreaseSubstepSize()

decrease the next subincrement

void finishElasticSubstep()

finish an elastic only subincrement

void finishSubstep(const TangentSizedMatrix &dXdY, const TangentSizedMatrix &dYdXOld)

finish an elastoplastic subincremnt

Parameters:

• dXdY – inverse of the material tangent

• dXdYOld – inverse of the previous material tangent

Matrix6d consistentStiffness()

get the consistent algorithmic tangent

PerezFougetSubstepper(double initialStepSize, double minimumStepSize, double scaleUpFactor, double scaleDownFactor, int nPassesToIncrease, const Matrix6d &Cel)

bool isFinished()

check if the subincrement process has finished

double getNextSubstep()

get the next subincrement size

bool decreaseSubstepSize()

decreas the subincrement size

void finishElasticSubstep()

finish an elastic only subincrement

void finishSubstep(const TangentSizedMatrix &dXdY)

finish a subincrement, pass the inverse material tangent

Matrix6d consistentStiffness()

get the overall consistent algorithmic tangent

Private Members

const double initialStepSize

const double minimumStepSize

const double scaleUpFactor

const double scaleDownFactor

const int nPassesToIncrease

double currentProgress

double currentSubstepSize

int passedSubsteps

const Matrix6d &Cel

MatrixStateStrain I76

TangentSizedMatrix I77

MatrixStateStrain consistentTangent

TangentSizedMatrix consistentTangent

template<int sizeMaterialState>
class PerezFougetSubstepperTime

#include <PerezFougetSubstepperTime.h>

Modified Version of the Perez-Fouget Substepper, to account for time-variant elastic Stiffness Tensor Cel(t_n+1) NO changes in algorithmic formulation!

modifications:

• elastic substep consistent tangent matrix update needs current Cel(t) for update

• No need for Cel for object construction anymore

Public Types

typedef Eigen::Matrix<double, sizeMaterialState, sizeMaterialState> TangentSizedMatrix

Matrix to carry the Jacobian of a material state.

Public Functions

PerezFougetSubstepperTime(double initialStepSize, double minimumStepSize, double scaleUpFactor, double scaleDownFactor, int nPassesToIncrease)

bool isFinished()

Check if the subincrementation has finished.

double getNextSubstep()

get the next subincrement size

double getFinishedProgress()

get the total finished progress of the subincrementation process

bool decreaseSubstepSize()

decrease the next subincrement

void extendConsistentTangent(const Matrix6d &CelT)

void extendConsistentTangent(const Matrix6d &CelT, const TangentSizedMatrix &matTangent)

Matrix6d consistentStiffness()

Private Members

const double initialStepSize

const double minimumStepSize

const double scaleUpFactor

const double scaleDownFactor

const int nPassesToIncrease

double currentProgress

double currentSubstepSize

int passedSubsteps

TangentSizedMatrix elasticTangent

TangentSizedMatrix consistentTangent

struct Properties

#include <MarmotPronySeries.h>

Material properties for a generalized Maxwell (Prony series) model.

This structure stores the parameters required for evaluating the viscoelastic stress response of a material represented by a Prony series expansion.

Public Members

size_t nPronyTerms

Number of Prony series terms.

Corresponds to the number of Maxwell branches in the generalized Maxwell model.

Matrix6d ultimateStiffnessMatrix

Ultimate (equilibrium) stiffness matrix.

Represents the stiffness contribution that persists at long times ( \(t \to \infty\)), outside the viscoelastic relaxation terms.

Matrix<double, 6, -1> pronyStiffnesses

Stiffness matrices of each Prony series term.

Stored as a block matrix of size \(6 \times (6 \cdot nPronyTerms)\), with each \(6\times 6\) block corresponding to the stiffness matrix of one Maxwell branch.

Matrix<double, 6, -1> pronyRelaxationTimes

Relaxation times of each Prony series term.

Stored as a block matrix of size \(6 \times (6 \cdot nPronyTerms)\), with each \(6\times 6\) block containing the relaxation times associated with one Maxwell branch.

struct SolverOptions

#include <MarmotMaterialPointSolverHypoElastic.h>

Struct to define solver options.

Contains parameters for controlling the solver’s behavior.

Public Members

int maxIterations = 25

Maximum number of iterations per increment.

double residualTolerance = 1e-10

Convergence tolerance.

double correctionTolerance = 1e-10

Correction tolerance.

struct StateVarEntryDefinition

#include <MarmotStateVarVectorManager.h>

An entry in the statevar vector consists of the name and a certain length.

Public Members

std::string name

int length

struct StateVarEntryLocation

#include <MarmotStateVarVectorManager.h>

The location in the statevar vector consists of the index and its certain length.

Public Members

int index

int length

struct StateVarVectorLayout

#include <MarmotStateVarVectorManager.h>

The layout is defined by a map of names to Locations, and the resulting required total length of the statevar vector

Public Members

std::unordered_map<std::string, StateVarEntryLocation> entries

int nRequiredStateVars

struct Step

#include <MarmotMaterialPointSolverHypoElastic.h>

Struct to define a loading step.

Each step contains target strain and stress states, time information and time step control parameters.

Public Functions

inline void checkControl() const

Check that for each component, either strain or stress is controlled.

Throws:

std::runtime_error – if the condition is not met

Public Members

Marmot::Vector6d strainIncrementTarget

Target strain increment for the step.

Marmot::Vector6d stressIncrementTarget

Target stress increment for the step.

Eigen::Vector<bool, 6> isStrainComponentControlled

Flags to indicate which strain components are controlled.

Eigen::Vector<bool, 6> isStressComponentControlled

Flags to indicate which stress components are controlled.

double timeStart = 0.0

Start time of the step.

double timeEnd = 1.0

End time of the step.

double dTStart = 0.1

Initial time step size.

double dTMin = 1e-6

Minimum time step size.

double dTMax = 0.5

Maximum time step size.

int maxIncrements = 100

Maximum number of increments in the step.

template<int nYieldSurfaces>
class YieldSurfaceCombinationManager

#include <YieldSurfaceCombinationManager.h>

Manager for yield surface combinations for multisurface plasticity: Try different yield surface combinations and track already used combinations

Public Types

typedef Eigen::Array<bool, 1, nYieldSurfaces> YieldSurfFlagArr

An array to carry active/nonactive states of yield surfaces.

typedef Eigen::Array<double, 1, nYieldSurfaces> YieldSurfResArr

An array to carry values of yield functions.

Public Functions

YieldSurfaceCombinationManager()

bool getAnotherYieldFlagCombination(YieldSurfFlagArr &activeSurfaces)

get another unused combination of yield surfaces

void markYieldFlagCombinationAsUsed(const YieldSurfFlagArr &activeSurfaces)

set the current combination as used

void resetUsedYieldFlagCombinations()

reset all yieldsurfaces as unused

Public Members

Eigen::Array<bool, (1 << nYieldSurfaces) - 1, (nYieldSurfaces + 1)> yieldSurfaceCombinations

An array which contains every possible (reasonable) combination of yield surfaces.

Private Members

const int idxUsedFlag

Column in the YieldSurfFlagArr for marking a combination as used.

namespace ContinuumMechanics

namespace VoigtNotation

This namespace includes functions needed for calculations with stress and strain tensors written in voigt notation.

namespace Eigen

namespace Marmot

Functions

void discardTheIncrement(double &pNewDT, double value, const std::string &message)

namespace ContinuumMechanics

namespace CommonConstitutiveModels

namespace Elasticity

namespace Isotropic

Functions for the description of isotropic elastic behavior.

Functions

double constexpr E(const double K, const double G)

Computes the isotropic Young’s modulus \( E \) from the compression modulus \( K \) and shear modulus \( G \).

\[ \displaystyle E = \frac{9\,K\,G}{3\,K + G} \]

Parameters:

• K – Compression modulus \( K \).

• G – Shear modulus \( G \).

Returns:

Young’s modulus \( E \).

double constexpr nu(const double K, const double G)

Computes the isotropic Poisson’s ratio \( \nu \) from the compression modulus \( K \) and shear modulus \( G \).

\[ \displaystyle \nu = \frac{3\,K - 2\,G}{6\,K + 2\,G} \]

Parameters:

• K – Compression modulus \( K \).

• G – Shear modulus \( G \).

Returns:

Poisson’s ratio \( \nu \).

double constexpr shearModulus(const double E, const double nu)

Computes the isotropic shear modulus \( G \) from the Young’s modulus \( E \) and Poisson’s ratio \( \nu \).

\[ \displaystyle G = \frac{E}{2\,(1 + \nu)} \]

Parameters:

• E – Young’s modulus \( E \).

• nu – Poisson’s ratio \( \nu \).

Returns:

Shear modulus \( G \).

double constexpr lameParameter(const double E, const double nu)

Computes the isotropic Lamé parameter \( \lambda \) from the Young’s modulus \( E \) and Poisson’s ratio \( \nu \).

\[ \displaystyle \lambda = \frac{E\,\nu}{(1 + \nu)(1 - 2\,\nu)} \]

Parameters:

• E – Young’s modulus \( E \).

• nu – Poisson’s ratio \( \nu \).

Returns:

Lamé parameter \( \lambda \).

Matrix6d complianceTensor(const double E, const double nu)

Computes the isotropic compliance tensor \(\mathbb{C}^{-1}\).

\[\begin{split} \displaystyle \mathbb{ C }^{-1} = \begin{bmatrix} \frac{1}{E} & \frac{-\nu}{E} & \frac{-\nu}{E} & 0 & 0 & 0 \\ \frac{-\nu}{E} & \frac{1}{E} & \frac{-\nu}{E} & 0 & 0 & 0 \\ \frac{-\nu}{E} & \frac{-\nu}{E} & \frac{1}{E} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G} \end{bmatrix} \end{split}\]

with

\[ G = \frac{E}{2\,(1 + \nu)}. \]

Parameters:

• E – Young’s modulus \(E\).

• nu – Poisson’s ratio \(\nu\).

Returns:

Compliance tensor \(\mathbb{C}^{-1}\) as a 6x6 matrix.

Matrix6d stiffnessTensor(const double E, const double nu)

Computes the isotropic stiffness tensor \(\mathbb{C}\) from Young’s modulus \(E\) and Poisson’s ratio \(\nu\).

\[\begin{split}\displaystyle \mathbb{ C } = \frac{E\,(1-\nu)}{(1+\nu)(1-2\,\nu)} \begin{bmatrix} 1 & \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 0 & 0 & 0 \\ \frac{\nu}{1-\nu} & 1 & \frac{\nu}{1-\nu} & 0 & 0 & 0 \\ \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\,\nu}{2\,(1-\nu)} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\,\nu}{2\,(1-\nu)} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\,\nu}{2\,(1-\nu)} \end{bmatrix} \end{split}\]

Parameters:

• E – Young’s modulus \(E\).

• nu – Poisson’s ratio \(\nu\).

Returns:

Stiffness tensor \(\mathbb{C}\) as a 6x6 matrix.

Matrix6d stiffnessTensorKG(const double K, const double G)

Computes the isotropic stiffness tensor \(\mathbb{C}\) from the bulk modulus \(K\) and shear modulus \(G\).

Parameters:

• K – Bulk modulus \(K\).

• G – Shear modulus \(G\).

Returns:

Stiffness tensor \(\mathbb{C}\) as a 6x6 matrix.

namespace Orthotropic

Functions for the description of orthotropic elastic behavior.

Functions

Matrix6d complianceTensor(const double E1, const double E2, const double E3, const double nu12, const double nu23, const double nu13, const double G12, const double G23, const double G31)

Computes the orthotropic compliance tensor \(\mathbb{C}^{-1}\), defined in the principal material directions \(x_1\), \(x_2\), and \(x_3\).

\[\begin{split} \displaystyle \mathbb{ C }^{-1} = \begin{bmatrix} \frac{1}{E_1} & \frac{-\nu_{12}}{E_2} & \frac{-\nu_{13}}{E_3} & 0 & 0 & 0 \\ \frac{-\nu_{12}}{E_2} & \frac{1}{E_2} & \frac{-\nu_{23}}{E_3} & 0 & 0 & 0 \\ \frac{-\nu_{13}}{E_3} & \frac{-\nu_{23}}{E_3} & \frac{1}{E_3} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G_{12}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G_{13}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G_{23}} \end{bmatrix} \end{split}\]

Parameters:

• E1 – Young’s modulus \(E_1\) in \( x_1 \) - direction.

• E2 – Young’s modulus \(E_2\) in \( x_2 \) - direction.

• E3 – Young’s modulus \(E_3\) in \( x_3 \) - direction.

• nu12 – Poisson’s ratio \(\nu_{12}\) effective between \( x_1 \) and \( x_2 \) - direction.

• nu23 – Poisson’s ratio \(\nu_{23}\) effective between \( x_2 \) and \( x_3 \) - direction.

• nu13 – Poisson’s ratio \(\nu_{13}\) effective between \( x_1 \) and \( x_3 \) - direction.

• G12 – Shear modulus \(G_{12}\) effective between \( x_1 \) and \( x_2 \) - direction.

• G23 – Shear modulus \(G_{23}\) effective between \( x_2 \) and \( x_3 \) - direction.

• G31 – Shear modulus \(G_{31}\) effective between \( x_1 \) and \( x_3 \) - direction.

Returns:

Compliance tensor \(\mathbb{C}^{-1}\) as a 6x6 matrix.

Matrix6d stiffnessTensor(const double E1, const double E2, const double E3, const double nu12, const double nu23, const double nu13, const double G12, const double G23, const double G31)

Computes the orthotropic stiffness tensor \( \mathbb{C} \) as inverse of the orthotropic compliance tensor \( \mathbb{C}^{-1} \).

Parameters:

• E1 – Young’s modulus \(E_1\) in \( x_1 \) - direction.

• E2 – Young’s modulus \(E_2\) in \( x_2 \) - direction.

• E3 – Young’s modulus \(E_3\) in \( x_3 \) - direction.

• nu12 – Poisson’s ratio \(\nu_{12}\) effective between \( x_1 \) and \( x_2 \) - direction.

• nu23 – Poisson’s ratio \(\nu_{23}\) effective between \( x_2 \) and \( x_3 \) - direction.

• nu13 – Poisson’s ratio \(\nu_{13}\) effective between \( x_1 \) and \( x_3 \) - direction.

• G12 – Shear modulus \(G_{12}\) effective between \( x_1 \) and \( x_2 \) - direction.

• G23 – Shear modulus \(G_{23}\) effective between \( x_2 \) and \( x_3 \) - direction.

• G31 – Shear modulus \(G_{31}\) effective between \( x_1 \) and \( x_3 \) - direction.

Returns:

Stiffness tensor \(\mathbb{C}\) as a 6x6 matrix.

namespace TransverseIsotropic

Functions for the description of transversely isotropic elastic behavior.

Functions

Matrix6d complianceTensor(const double E1, const double E2, const double nu12, const double nu23, const double G12)

Computes the transversely isotropic compliance tensor \(\mathbb{C}^{-1}\).

\[\begin{split} \displaystyle \mathbb{ C }^{-1} = \begin{bmatrix} \frac{1}{E_1} & \frac{-\nu_{12}}{E_1} & \frac{-\nu_{12}}{E_1} & 0 & 0 & 0 \\ \frac{-\nu_{12}}{E_1} & \frac{1}{E_2} & \frac{-\nu_{23}}{E_2} & 0 & 0 & 0 \\ \frac{-\nu_{12}}{E_1} & \frac{-\nu_{23}}{E_2} & \frac{1}{E_2} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G_{12}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G_{12}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G_{23}} \end{bmatrix} \end{split}\]

The isotropic plane is defined with respect to the \(x_2\)– \(x_3\) axes of a local coordinate system. The tensor is expressed in terms of the out-of-plane Young’s modulus \(E_1\), shear modulus \(G_{12}\), and Poisson’s ratio \(\nu_{12}\), together with the in-plane Young’s modulus \(E_2\) and Poisson’s ratio \(\nu_{23}\).

The in-plane shear modulus \( G_{23} \) is given by

\[ \displaystyle G_{23} = \frac{E_2}{2\,(1 + \nu_{23})} \]

Parameters:

• E1 – Out-of-plane Young’s modulus \(E_1\).

• E2 – In-plane Young’s modulus \(E_2\).

• nu12 – Poisson’s ratio \(\nu_{12}\).

• nu23 – Poisson’s ratio \(\nu_{23}\).

• G12 – shear modulus \(G_{12}\).

Returns:

Compliance tensor \(\mathbb{C}^{-1}\) as a 6x6 matrix.

Matrix6d stiffnessTensor(const double E1, const double E2, const double nu12, const double nu23, const double G12)

Computes the transversely isotropic stiffness tensor \( \mathbb{C} \) as inverse of the transversely isotropic compliance tensor \( \mathbb{C}^{-1} \).

Parameters:

• E1 – Out-of-plane Young’s modulus \(E_1\).

• E2 – In-plane Young’s modulus \(E_2\).

• nu12 – Poisson’s ratio \(\nu_{12}\).

• nu23 – Poisson’s ratio \(\nu_{23}\).

• G12 – shear modulus \(G_{12}\).

Returns:

Stiffness tensor \(\mathbb{C}\) as a 6x6 matrix.

namespace HaighWestergaard

Functions

template<typename T = double>
HaighWestergaardCoordinates<T> haighWestergaard(const Eigen::Matrix<T, 6, 1> &stress)

Computes the stress coordinates in the Haigh-Westergaard space.

\[\xi=I_1/\sqrt{3}\]

\[\rho=\sqrt{2\,J_2}\]

\[\theta=\frac{1}{3}\,\arccos{\left(\frac{3\,\sqrt{3}}{2}\,\frac{J_3}{\sqrt{J_2^3}}\right)}\]

Note

The stress coordinates are computed from the invariants \(I_1,\,J_2,\,J_3\) of the stress tensor as follows:

Parameters:

stress – Stress tensor \(\sig\) given in Voigt notation.

HaighWestergaardCoordinates<double> haighWestergaardFromStrain(const Marmot::Vector6d &strain)

Computes the strain coordinates in the Haigh-Westergaard space.

Note

The computation is equal to haighWestergaard by replacing the stress invariants with the strain invariants.

Parameters:

strain – Strain tensor \(\eps\) given in Voigt notation.

namespace Kinematics

namespace DeformationGradient

Functions

template<int nDim>
Eigen::Matrix3d make3D(const Eigen::Ref<const Eigen::Matrix<double, nDim, nDim>> &tensor)

Embeds a lower-dimensional square tensor into a 3x3 tensor.

This function takes an \(n \times n\) tensor (with \(n=1,2,3\)) and returns a 3x3 tensor:

• For \(n=1\): the scalar is placed in the (0,0) entry of a 3x3 identity matrix.

• For \(n=2\): the 2x2 tensor is placed in the top-left block of a 3x3 identity matrix.

• For \(n=3\): the tensor is returned unchanged.

Template Parameters:

nDim – Dimension of the input tensor (1, 2, or 3).

Parameters:

tensor – An \(n \times n\) tensor represented as Eigen::Matrix<double,nDim,nDim>.

Returns:

A 3x3 tensor embedding the input according to the rules above.

namespace Strain

Functions

Marmot::Vector6d GreenLagrange(const Eigen::Matrix3d &F)

Computes the Green-Lagrange strain tensor from the deformation gradient.

Parameters:

F – The deformation gradient tensor as a 3x3 matrix.

Returns:

The Green-Lagrange strain tensor represented as a 6D vector in Voigt notation.

namespace VelocityGradient

Variables

const Eigen::TensorFixedSize<double, Eigen::Sizes<3, 3, 3, 3>> dOmega_dVelocityGradient

A 4th-order tensor representing the derivative of the spin tensor \( \Omega \) with respect to the velocity gradient tensor.

const Eigen::TensorFixedSize<double, Eigen::Sizes<6, 3, 3>> dStretchingRate_dVelocityGradient

Tensor representing the derivative of the stretching rate with respect to the velocity gradient.

namespace LocalizationAnalysis

Functions

Marmot::Matrix3d computeAcousticTensor(const Marmot::Matrix6d &materialTangent, const Marmot::Vector3d &normalVector)

bool localizationChecker(const Marmot::Matrix3d &acousticTensor)

Marmot::Vector3d computeNormalVector(double alpha, double beta)

double minimumDeterminantAcousticTensor(const Marmot::Matrix6d &materialTangent)

namespace PlaneStrain

Functions

Eigen::Matrix3d getPlaneStrainTangent(const Matrix6d &C)

Extracts the plane strain stiffness matrix from a 3D stiffness matrix.

Reduces the full 3D stiffness matrix \(\mathbb{C}\) to the plane strain stiffness matrix corresponding to the in-plane components.

Parameters:

C – 3D stiffness matrix ( \(6\times6\)) in Voigt notation.

Returns:

Plane strain stiffness matrix ( \(3\times3\)) in Voigt notation.

EigenTensors::Tensor322d reduce3D_dStress_dDeformationGradient(const EigenTensors::Tensor633d &dStressdDeformationGradient3D)

Extracts the plane strain derivative of the stress with respect to the deformation gradient.

Reduces the derivative of the 3D stress tensor in Voigt notation with respect to the deformation gradient \(F_{ij}\) to the plane strain setting.

Parameters:

dStressdDeformationGradient3D – Derivative of 3D stress with respect to the deformation gradient ( \(6\times3\times3\) tensor).

Returns:

Reduced derivative of stress with respect to deformation gradient ( \(3\times2\times2\) tensor) in the plane strain setting.

Eigen::Matrix<double, 6, 3> dStrainDStrainPlaneStrain()

Computes the derivative of the 3D strain tensor with respect to the plane strain tensor.

Returns:

Matrix ( \(6\times3\)) representing \(\frac{\partial \varepsilon_{3D}}{\partial \varepsilon_{plane}}\).

namespace PlaneStress

Functions

Eigen::Matrix3d getPlaneStressTangent(const Matrix6d &C)

Extracts the plane stress stiffness matrix from a 3D stiffness matrix.

Parameters:

C – 3D stiffness matrix ( \(6\times6\)) in Voigt notation.

Returns:

Plane stress stiffness matrix ( \(3\times3\)) in Voigt notation.

EigenTensors::Tensor322d compute_dStress_dDeformationGradient(const EigenTensors::Tensor633d &dStressdDeformationGradient3D)

Computes the plane stress derivative of the stress with respect to the deformation gradient.

Reduces the 3D derivative of the stress tensor with respect to the deformation gradient \( F_{ij} \) to a plane stress setting.

Parameters:

dStressdDeformationGradient3D – Derivative of 3D stress w.r.t. deformation gradient ( \(6 \times 3 \times 3\) tensor).

Returns:

Reduced derivative of stress w.r.t. deformation gradient ( \(3 \times 2 \times 2\) tensor) for plane stress.

Marmot::Vector6d planeStressCompensationStrain(const Marmot::Vector6d &elasticStrain, double nu)

Computes the out-of-plane strain correction component for plane stress conditions.

Determines the correction \(\varepsilon_{33}^{\text{corr}}\) such that the resulting stress satisfies plane stress conditions. \( \boldsymbol{\varepsilon}^{\text{el}} \), so that the resulting stress satisfies plane stress:

\[ \sigma^\text{plane} = \mathbb{C} : (\boldsymbol{\varepsilon}^{\text{el}} + \boldsymbol{\varepsilon}^{\text{corr}}) \]

Parameters:

• elasticStrain – The in-plane elastic strain vector \( \boldsymbol{\varepsilon}^{\text{el}} \) (6 components, Voigt notation).

• nu – Poisson’s ratio of the material.

Returns:

A 6-component strain vector \( \boldsymbol{\varepsilon}^{\text{corr}} \) containing the computed out-of-plane correction component \(\varepsilon_{33}^{\text{corr}}\) (Voigt notation).

Matrix6d planeStressTangentTransformationMatrix(const Matrix6d &tangent)

Computes the transformation matrix for the plane stress tangent.

Returns the 6x6 transformation matrix \(\mathbf{T}\) relating arbitrary strain increments to plane stress increments using the material stiffness.

Parameters:

tangent – 6x6 material tangent stiffness matrix ( \(\mathbf{C}\)).

Returns:

6x6 transformation matrix ( \(\mathbf{T}\)) for plane stress correction.

Eigen::Matrix<double, 6, 3> dStrainDStrainPlaneStress(const Matrix6d &tangent)

Computes the derivative of the 3D strain tensor with respect to the plane stress strain tensor.

Parameters:

tangent – 6x6 material tangent stiffness matrix ( \(\mathbf{C}\)).

Returns:

6x3 matrix representing \(\frac{\partial \varepsilon_{3D}}{\partial \varepsilon_{\sigma_{plane}}}\).

Eigen::Matrix<double, 3, 6> dStressPlaneStressDStress()

Computes the derivative of the plane stress tensor with respect to the 3D stress tensor.

Returns:

3x6 matrix representing \(\frac{\partial \sigma_{plane}}{\partial \sigma_{3D}}\).

namespace UniaxialStress

Functions

double getUniaxialStressTangent(const Eigen::Ref<const Matrix6d> &C)

Extracts the uniaxial stress-strain tangent from a 3D stiffness matrix.

Computes the tangent \(\frac{\partial\sigma}{\partial\varepsilon}\) in the loading direction (uniaxial) based on the full 3D stiffness matrix \(\mathbb{C}\) provided in Voigt notation.

Parameters:

C – 3D stiffness matrix ( \(6\times6\)) in Voigt notation.

Returns:

Uniaxial stress-strain tangent ( \(\frac{\partial\sigma}{\partial\varepsilon}\)).

namespace Viscoelasticity

namespace ComplianceFunctions

Functions

template<typename T_>
T_ logPowerLaw(T_ tau, double m, double n)

Logarithmic power-law compliance function.

Computes the compliance function

\[ \Phi(\tau) = m \, \ln\!\left( 1 + \tau^n \right), \]

creep for long times.

Template Parameters:

T_ – Scalar or autodiff type of the argument.

Parameters:

• tau – [in] Retardation time or evaluation point.

• m – [in] Scaling factor.

• n – [in] Exponent controlling the growth rate.

Returns:

The compliance value \(\Phi(\tau)\).

template<typename T_>
T_ powerLaw(T_ tau, double m, double n)

Power-law compliance function.

Computes the compliance function

\[ \Phi(\tau) = m \, \tau^n, \]

creep for short times.

Template Parameters:

T_ – Scalar or autodiff type of the argument.

Parameters:

• tau – [in] Retardation time or evaluation point.

• m – [in] Scaling factor.

• n – [in] Exponent controlling the growth rate.

Returns:

The compliance value \(\Phi(\tau)\).

namespace PronySeries

Typedefs

typedef Eigen::Matrix<double, 6, Eigen::Dynamic> StateVarMatrix

Matrix of viscoelastic state variables.

Dynamic-size matrix of shape \([6 \times n]\), where each \(6\times 1\) block corresponds to the internal state variables of a Maxwell branch in the Prony series.

typedef Eigen::Map<StateVarMatrix> mapStateVarMatrix

Mapped view of a state variable matrix.

Provides an Eigen::Map wrapper around a contiguous block of memory that stores state variables in the same layout as StateVarMatrix, avoiding copies.

Functions

void evaluatePronySeries(const Properties &props, Vector6d &stress, Matrix6d &stiffness, Eigen::Ref<mapStateVarMatrix> stateVars, const Vector6d &dStrain, const double dT, const bool updateStateVars = false)

Evaluate the Prony series viscoelastic response.

Computes the stress and tangent stiffness matrix contributions of a generalized Maxwell (Prony series) viscoelastic material, and optionally updates the internal state variables.

In what follows we use indicial notation without the Einstein summation over repeated indices except explicitly stated. The stress is updated as

\[ \sigma^{k+1}_{ij} =\; \sum_{k,l}C^\infty_{ijkl} \Delta \varepsilon_{kl} \;+\; \sum_{\mu=1}^{M} \left[ \sum_{k,l}\eta^\mu_{ijkl}\left(1 - e^{-\Delta t / \tau^{\mu}_{ij}}\right)\frac{\Delta\varepsilon_{kl}}{\Delta t} \;-\; \sigma^{ve,\mu}_{ij}\left(1 - e^{-\Delta t/\tau^\mu_{ij}}\right) \right], \]

where \(C^\infty_{ijkl}\) is the ultimate stiffness, \(\tau^\mu_{ij}\) the relaxation times for each strain component, \(\eta^\mu_{ijkl} = \tau^\mu_{ij} C^\mu_{ijkl}\) the effective viscosities, and \(\sigma^{ve,k}_{kl}\) the internal state variables.

Parameters:

• props – [in] Prony series material properties (relaxation times, stiffnesses, number of terms, and ultimate stiffness).

• stress – [inout] 6-component stress vector to be updated.

• stiffness – [inout] 6×6 tangent stiffness matrix to be updated.

• stateVars – [inout] State variable matrix of size [6×(6·nPronyTerms)], storing internal viscoelastic strain-like variables.

• dStrain – [in] Incremental strain vector (6-components).

• dT – [in] Time increment.

• updateStateVars – [in] If true, update the state variables for the next step.

void updateStateVars(const Properties &props, Eigen::Ref<mapStateVarMatrix> stateVars, const Vector6d &dStrain, const double dT)

Update the Prony series state variables.

Propagates the internal viscoelastic state variables

In what follows we don’t use the Einstein summation over repeated indices, the products are to be taken as hadamard products except otherwise stated. \(\sigma^{ve,\mu,k}_{kl}\) from time step \(t^k\) to \(t^{k+1}\) according to

\[ \sigma^{ve,\mu,k+1}_{kl} = e^{-\Delta t / \tau^\mu_{ij}}\, \sigma^{ve,\mu,k}_{kl} + \sum_{k,l}\left( \frac{\eta^\mu_{ijkl}}{\Delta t}(1 - e^{-\Delta t / \tau^\mu_{ij}}) \right)\Delta\varepsilon_{kl}, \]

where \(\tau^\mu_{ij}\) are the relaxation times and \(\eta^\mu_{ij} = \tau^\mu_{ij} C^\mu_{ijkl}\).

Parameters:

• props – [in] Prony series material properties (relaxation times, stiffnesses, and number of terms).

• stateVars – [inout] State variable matrix of size [6×(6·nPronyTerms)], storing internal viscoelastic strain-like variables.

• dStrain – [in] Incremental strain vector (6-components).

• dT – [in] Time increment.

namespace VoigtNotation

Enums

enum VoigtSize

Values:

enumerator OneD

enumerator TwoD

enumerator ThreeD

enumerator Axial

Functions

VoigtSize voigtSizeFromDimension(int x)

Computes the size of the Voigt notation array based on the spatial dimension.

Note

The calculation is based on the formula: \( \frac{x * x + x}{2} \).

Parameters:

x – The spatial dimension (e.g., 2 for 2D, 3 for 3D).

Returns:

The size of the Voigt notation array as a VoigtSize.

Eigen::Vector3d voigtToPlaneVoigt(const Marmot::Vector6d &voigt)

Converts a 3D stress vector in Voigt notation to a plane stress vector.

Note

This function extracts only the in-plane components of a 3D stress vector.

Parameters:

voigt – 6-component 3D stress vector in Voigt notation:

\[\begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{33} & \sigma_{12} & \sigma_{13} & \sigma_{23} \end{bmatrix}^T \]

where \(\sigma_{33} = \sigma_{13} = \sigma_{23} = 0\) in plane stress conditions.

Returns:

3-component plane stress vector in Voigt notation:

\[\begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{12} \end{bmatrix}^T \]

Vector4d voigtToAxisymmetricVoigt(const Vector6d &voigt)

Marmot::Vector6d planeVoigtToVoigt(const Eigen::Vector3d &voigtPlane)

Converts a plane stress vector in Voigt notation back to a 3D Voigt vector.

Note

The out-of-plane components are set to zero, suitable for plane stress conditions.

Parameters:

voigtPlane – 3-component plane stress vector in Voigt notation:

\[\begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{12} \end{bmatrix}^T \]

Returns:

6-component 3D stress vector in Voigt notation:

\[\begin{bmatrix} \sigma_{11} & \sigma_{22} & 0 & \sigma_{12} & 0 & 0 \end{bmatrix}^T \]

Vector6d axisymmetricVoigtToVoigt(const Vector4d &voigtAxisymmetric)

template<enum VoigtSize voigtSize>
Eigen::Matrix<double, voigtSize, 1> reduce3DVoigt(const Marmot::Vector6d &Voigt3D)

Reduces a 3D Voigt vector to a lower-dimensional vector defined by voigtSize.

Note

Considered cases:

• voigtSize = 1: Returns the \(\sigma_{11}\) component.

• voigtSize = 3: Returns the plane stress components using voigtToPlaneVoigt().

• voigtSize = 6: Returns the input 3D vector unchanged.

Template Parameters:

voigtSize – Target Voigt dimension (1, 3, or 6).

Parameters:

Voigt3D – Input 6-component 3D Voigt vector.

Throws:

std::invalid_argument – if voigtSize is not 1, 3, or 6.

Returns:

Reduced Voigt vector of size voigtSize.

template<enum VoigtSize voigtSize>
Marmot::Vector6d make3DVoigt(const Eigen::Matrix<double, voigtSize, 1> &Voigt)

Constructs a 3D Voigt vector from a lower-dimensional Voigt vector.

Note

Considered cases:

• voigtSize = 1: Returns a 3D vector with only \(\sigma_{11}\) non-zero:

  \[\begin{bmatrix} \sigma_{11} & 0 & 0 & 0 & 0 & 0 \end{bmatrix}^T \]

• voigtSize = 3: Returns a 3D Voigt vector by calling planeVoigtToVoigt().

• voigtSize = 6: Returns the input vector unchanged.

Template Parameters:

voigtSize – Dimension of the input Voigt vector (1, 3, or 6).

Parameters:

Voigt – Input Voigt vector of size voigtSize.

Throws:

std::invalid_argument – if voigtSize is not 1, 3, or 6.

Returns:

6-component 3D Voigt vector.

template<typename T>
Eigen::Matrix<T, 3, 3> voigtToStrain(const Eigen::Matrix<T, 6, 1> &voigt)

Converts a 6-component Voigt strain vector to a 3x3 strain tensor.

Note

Shear components in Voigt notation are engineering strains and are halved in the resulting tensor.

Template Parameters:

T – Scalar type (e.g., double, float).

Parameters:

voigt – 6-component strain vector in Voigt notation:

\[\begin{bmatrix} \varepsilon_{11} & \varepsilon_{22} & \varepsilon_{33} & \gamma_{12} & \gamma_{13} & \gamma_{23} \end{bmatrix}^T \]

Returns:

3x3 strain tensor:

\[\begin{split}\begin{bmatrix} \varepsilon_{11} & \gamma_{12}/2 & \gamma_{13}/2 \\ \gamma_{12}/2 & \varepsilon_{22} & \gamma_{23}/2 \\ \gamma_{13}/2 & \gamma_{23}/2 & \varepsilon_{33} \\ \end{bmatrix} \end{split}\]

template<typename T>
Eigen::Matrix<T, 3, 3> voigtToStress(const Eigen::Matrix<T, 6, 1> &voigt)

Converts a 6-component Voigt stress vector to a 3x3 stress tensor.

Template Parameters:

T – Scalar type (e.g., double, float).

Parameters:

voigt – 6-component stress vector in Voigt notation:

\[\begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{33} & \sigma_{12} & \sigma_{13} & \sigma_{23} \end{bmatrix}^T \]

Returns:

3x3 stress tensor:

\[\begin{split}\begin{bmatrix} \sigma_{11} & \sigma_{12} & \sigma_{13} \\ \sigma_{12} & \sigma_{22} & \sigma_{23} \\ \sigma_{13} & \sigma_{23} & \sigma_{33} \\ \end{bmatrix} \end{split}\]

Marmot::Vector6d strainToVoigt(const Eigen::Matrix3d &strainTensor)

Converts a 3x3 strain tensor to a 6-component Voigt strain vector.

Note

Shear components in the Voigt vector ( \(\gamma_{ij}\) ) are engineering strains, which are twice the tensor shear components.

Parameters:

strainTensor – 3x3 strain tensor:

\[\begin{split}\boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_{11} & \varepsilon_{12} & \varepsilon_{13} \\ \varepsilon_{21} & \varepsilon_{22} & \varepsilon_{23} \\ \varepsilon_{31} & \varepsilon_{32} & \varepsilon_{33} \\ \end{bmatrix} \end{split}\]

Returns:

6-component strain vector in Voigt notation:

\[\begin{split}\begin{aligned} \boldsymbol{\varepsilon}^{\text{Voigt}} &= \begin{bmatrix} \varepsilon_{11} & \varepsilon_{22} & \varepsilon_{33} & \gamma_{12} & \gamma_{13} & \gamma_{23} \end{bmatrix}^T \\ &= \begin{bmatrix} \varepsilon_{11} & \varepsilon_{22} & \varepsilon_{33} & 2 \varepsilon_{12} & 2 \varepsilon_{13} & 2 \varepsilon_{23} \end{bmatrix}^T \end{aligned} \end{split}\]

template<typename T = double>
Eigen::Matrix<T, 6, 1> stressToVoigt(const Eigen::Matrix<T, 3, 3> &stressTensor)

Converts a 3x3 stress tensor to a 6-component Voigt stress vector.

Note

Shear components are copied directly; no factor of 2 is applied.

Template Parameters:

T – Scalar type (e.g., double, float).

Parameters:

stressTensor – 3x3 stress tensor:

\[\begin{split}\boldsymbol{\sigma} = \begin{bmatrix} \sigma_{11} & \sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_{22} & \sigma_{23} \\ \sigma_{31} & \sigma_{32} & \sigma_{33} \\ \end{bmatrix} \end{split}\]

Returns:

6-component stress vector in Voigt notation:

\[\boldsymbol{\sigma}^{\text{Voigt}} = \begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{33} & \sigma_{12} & \sigma_{13} & \sigma_{23} \end{bmatrix}^T \]

Eigen::Matrix<double, 6, 6> stiffnessToVoigt(const Eigen::Tensor<double, 4> &C)

Converts a fourth-order stiffness tensor to its Voigt notation representation ( \( 6 \times 6 \) matrix).

Note

The input tensor C must follow the symmetry properties of a stiffness tensor for the conversion to be valid (minor symmetry).

Parameters:

C – The fourth-order stiffness tensor represented as an Eigen::Tensor<double, 4>.

Returns:

An Eigen::Matrix<double, 6, 6> representing the stiffness tensor in Voigt notation.

Eigen::Tensor<double, 4> voigtToStiffness(const Eigen::Matrix<double, 6, 6> &voigtStiffness)

Converts a stiffness matrix in Voigt notation ( \( 6 \times 6 \) matrix) to a 4th-order stiffness tensor ( \( 3 \times 3 \times 3 \times 3 \) tensor).

Parameters:

voigtStiffness – The \( 6 \times 6 \) matrix representing the stiffness in Voigt notation.

Returns:

An Eigen::Tensor of rank 4 (4th-order tensor) representing the stiffness tensor. The dimensions of the tensor are \( 3 \times 3 \times 3 \times 3 \).

template<int nDim>
Eigen::Matrix<double, nDim, nDim> stressMatrixFromVoigt(const Eigen::Matrix<double, VOIGTFROMDIM(nDim), 1> &Voigt)

Converts a stress vector in Voigt notation to its corresponding tensor form.

This function handles 1D, 2D, and 3D stress vectors in Voigt notation and returns the full symmetric stress tensor.

Note

For 2D and 3D, the off-diagonal shear components in the tensor are taken directly from the Voigt vector (no factor of 2 applied).

Template Parameters:

nDim – Dimension of the stress tensor (1, 2, or 3).

Parameters:

Voigt – Input stress vector in Voigt notation, size VOIGTFROMDIM(nDim) x 1:

• 1D: \( \boldsymbol{\sigma}^{\text{Voigt}} = \begin{bmatrix} \sigma_{11} \end{bmatrix}^T \)

• 2D: \( \boldsymbol{\sigma}^{\text{Voigt}} = \begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{12} \end{bmatrix}^T \)

• 3D: \( \boldsymbol{\sigma}^{\text{Voigt}} = \begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{33} & \sigma_{12} & \sigma_{13} & \sigma_{23} \end{bmatrix}^T \)

Throws:

std::invalid_argument – If nDim is not 1, 2, or 3.

Returns:

Symmetric stress tensor of size \( nDim \times nDim \).

template<int nDim>
Eigen::Matrix<double, VOIGTFROMDIM(nDim), 1> voigtFromStrainMatrix(const Eigen::Matrix<double, nDim, nDim> &strain)

Converts a strain tensor to its corresponding Voigt notation vector.

This function supports 1D, 2D, and 3D strain tensors and returns the corresponding Voigt vector (column vector). Shear components are scaled appropriately for engineering strain.

Note

Axisymmetric voigt vectors are not supported.

Template Parameters:

nDim – Dimension of the strain tensor (1, 2, or 3).

Parameters:

strain – Input strain tensor of size nDim x nDim.

Throws:

std::invalid_argument – If nDim is not 1, 2, or 3.

Returns:

Column vector of size VOIGTFROMDIM(nDim) x 1 representing the strain in Voigt notation.

• 1D: \( \boldsymbol{\varepsilon}^{\text{Voigt}} = \begin{bmatrix} \varepsilon_{11} \end{bmatrix}^T \)

• 2D: \( \boldsymbol{\varepsilon}^{\text{Voigt}} = \begin{bmatrix} \varepsilon_{11} & \varepsilon_{22} & 2 \varepsilon_{12} \end{bmatrix}^T \)

• 3D: \( \boldsymbol{\varepsilon}^{\text{Voigt}} = \begin{bmatrix} \varepsilon_{11} & \varepsilon_{22} & \varepsilon_{33} & 2 \varepsilon_{12} & 2 \varepsilon_{13} & 2 \varepsilon_{23} \end{bmatrix}^T \)

Variables

const Marmot::Vector6d P

Predefined 6D vector with scaling factor (2x) for shear components in Voigt notation.

The vector contains scaling factors: {1, 1, 1, 2, 2, 2}.

const Marmot::Vector6d PInv

Predefined 6D vector with inverse scaling factor (0.5x) for shear components in Voigt notation.

The vector contains inverse scaling factors: {1, 1, 1, 0.5, 0.5, 0.5}.

const Marmot::Vector6d I

Predefined 6D Vector representing the identity tensor in Voigt notation.

The vector contains ones in all components: {1, 1, 1, 0, 0, 0}.

const Marmot::Vector6d IHyd

Predefined 6D Vector representing the hydrostatic projection tensor in Voigt notation.

The vector contains: {1/3, 1/3, 1/3, 0, 0, 0}.

const Matrix6d IDev

Deviatoric projection tensor in Voigt notation.

\( 6 \times 6 \) matrix.

namespace Derivatives

Functions

Vector6d dStressMean_dStress()

Computes the derivative of the mean stress with respect to the stress vector.

\[ \frac{\partial \sigma_m}{\partial \boldsymbol{\sigma}} \]

Returns:

6-component vector representing \( \tfrac{\partial \sigma_m}{\partial \boldsymbol{\sigma}} \) in Voigt notation.

template<typename T>
Eigen::Matrix<T, 6, 1> dRho_dStress(T rho, const Eigen::Matrix<T, 6, 1> &stress)

Computes the derivative of the Haigh–Westergaard deviatoric radius \( \rho \) with respect to the stress vector.

\[ \frac{\partial \rho}{\partial \boldsymbol{\sigma}} \]

Note

Returns zero if \( \rho \leq 10^{-16} \) to avoid division by zero.

Template Parameters:

T – Scalar type (e.g., double, float).

Parameters:

• rho – Precomputed Haigh–Westergaard coordinate \( \rho \).

• stress – 6-component stress vector in Voigt notation.

Returns:

6-component vector of partial derivatives \( \tfrac{\partial \rho}{\partial \boldsymbol{\sigma}} \) in Voigt notation.

template<typename T>
Eigen::Matrix<T, 6, 1> dRhoStrain_dStrain(T rhoStrain, const Eigen::Matrix<T, 6, 1> &strain)

Computes the derivative \( \frac{d\, \varepsilon_\rho}{d\, \boldsymbol{\varepsilon}}\) of the haigh westergaard coordinate \(\strain_\rho \) with respect to the voigt notated strain vector \( \boldsymbol{\varepsilon} \)

Marmot::Vector6d dTheta_dStress(double theta, const Marmot::Vector6d &stress)

Computes the derivative of the lode angle \( \theta \) with respect to the stress vector.

Note

Returns zero if \( \theta \leq 10^{-15} \) or if \( \theta \geq \frac{\pi}{3} - 10^{-15} \).

Parameters:

• theta – Lode angle \( \theta \).

• stress – 6-component stress vector in Voigt notation.

Returns:

6-component vector of partial derivatives \( \tfrac{\partial \theta}{\partial \boldsymbol{\sigma}} \) in Voigt notation.

double dTheta_dJ2(const Marmot::Vector6d &stress)

Computes the derivative \( \frac{\partial \theta}{\partial J_2}\) of the lode angle \( \theta \) with respect to the second deviatoric invariant \( J_2 \).

Parameters:

stress – 6-component stress vector in Voigt notation.

Returns:

Scalar derivative \( \tfrac{\partial \theta}{\partial J_2} \).

double dTheta_dJ3(const Marmot::Vector6d &stress)

Computes the derivative \( \frac{\partial \theta}{\partial J_3}\) of the Lode angle \( \theta \) with respect to the third deviatoric invariant \( J_3 \).

Parameters:

stress – Stress vector in Voigt notation ( \( \boldsymbol{\sigma} \)).

Returns:

Value of the derivative \( \frac{\partial \theta}{\partial J_3} \).

double dThetaStrain_dJ2Strain(const Marmot::Vector6d &strain)

Computes the derivative \( \frac{\partial \theta^{(\varepsilon)}}{\partial J^{(\varepsilon)}_2}\) of the Lode angle in strain space \( \theta^{(\varepsilon)} \) with respect to the second deviatoric invariant \( J^{(\varepsilon)}_2 \).

Parameters:

strain – Strain vector in Voigt notation ( \( \boldsymbol{\varepsilon} \)).

Returns:

Value of the derivative \( \frac{\partial \theta^{(\varepsilon)}}{\partial J^{(\varepsilon)}_2} \).

double dThetaStrain_dJ3Strain(const Marmot::Vector6d &strain)

Computes the derivative \( \frac{\partial \theta^{(\varepsilon)}}{\partial J^{(\varepsilon)}_3}\) of the Lode angle in strain space \( \theta^{(\varepsilon)} \) with respect to the third deviatoric invariant \( J^{(\varepsilon)}_3 \).

Parameters:

strain – Strain vector in Voigt notation ( \( \boldsymbol{\varepsilon} \)).

Returns:

Value of the derivative \( \frac{\partial \theta^{(\varepsilon)}}{\partial J^{(\varepsilon)}_3} \).

Marmot::Vector6d dJ2_dStress(const Marmot::Vector6d &stress)

Computes the derivative \( \frac{\partial J_2}{\partial \boldsymbol{\sigma}}\) of the second deviatoric invariant \( J_2 \) with respect to the stress vector in Voigt notation.

Parameters:

stress – Stress vector in Voigt notation ( \( \boldsymbol{\sigma} \)).

Returns:

Vector of partial derivatives \( \frac{\partial J_2}{\partial \boldsymbol{\sigma}} \) in Voigt notation.

Marmot::Vector6d dJ3_dStress(const Marmot::Vector6d &stress)

Computes the derivative \( \frac{\partial J_3}{\partial \boldsymbol{\sigma}}\) of the third deviatoric invariant \( J_3 \) with respect to the stress vector in Voigt notation.

Parameters:

stress – Stress vector in Voigt notation ( \( \boldsymbol{\sigma} \)).

Returns:

Vector of partial derivatives \( \frac{\partial J_3}{\partial \boldsymbol{\sigma}} \) in Voigt notation.

Marmot::Vector6d dJ2Strain_dStrain(const Marmot::Vector6d &strain)

Computes the derivative \( \frac{\partial J^{(\varepsilon)}_2}{\partial \boldsymbol{\varepsilon}} \) of the second deviatoric invariant \( J^{(\varepsilon)}_2 \) with respect to the strain vector in Voigt notation.

Parameters:

strain – Strain vector in Voigt notation ( \( \boldsymbol{\varepsilon} \)).

Returns:

Vector of partial derivatives \( \frac{\partial J^{(\varepsilon)}_2}{\partial \boldsymbol{\varepsilon}} \) in Voigt notation.

Marmot::Vector6d dJ3Strain_dStrain(const Marmot::Vector6d &strain)

Computes the derivative \( \frac{\partial J^{(\varepsilon)}_3}{\partial \boldsymbol{\varepsilon}} \) of the third deviatoric invariant \( J^{(\varepsilon)}_3 \) with respect to the strain vector in Voigt notation.

Parameters:

strain – Strain vector in Voigt notation ( \( \boldsymbol{\varepsilon} \)).

Returns:

Vector of partial derivatives \( \frac{\partial J^{(\varepsilon)}_3}{\partial \boldsymbol{\varepsilon}} \) in Voigt notation.

Marmot::Vector6d dThetaStrain_dStrain(const Marmot::Vector6d &strain)

Computes the derivative \( \frac{\partial \theta^{(\varepsilon)}}{\partial \boldsymbol{\varepsilon}} \) of the Lode angle in strain space \( \theta^{(\varepsilon)} \) with respect to the strain vector in Voigt notation.

Parameters:

strain – Strain vector in Voigt notation ( \( \boldsymbol{\varepsilon} \)).

Returns:

Vector of partial derivatives \( \frac{\partial \theta^{(\varepsilon)}}{\partial \boldsymbol{\varepsilon}} \) in Voigt notation.

Marmot::Matrix36 dStressPrincipals_dStress(const Marmot::Vector6d &stress)

Computes the derivative \( \frac{\partial \sigma_I}{\partial \boldsymbol{\sigma}} \) of the principal stresses \( \sigma_I \) with respect to the stress vector in Voigt notation.

Parameters:

stress – Stress vector in Voigt notation ( \( \boldsymbol{\sigma} \)).

Returns:

Matrix of partial derivatives \( \frac{\partial \sigma_I}{\partial \boldsymbol{\sigma}} \), where each row corresponds to the gradient of one principal stress with respect to \( \boldsymbol{\sigma} \).

Eigen::Vector3d dStrainVolumetricNegative_dStrainPrincipal(const Marmot::Vector6d &strain)

Computes the derivative \( \frac{\partial \varepsilon^{vol}_{\ominus}}{\partial \varepsilon_I} \) of the volumetric compressive strain \( \varepsilon^{vol}_{\ominus} \) with respect to the principal strains \( \varepsilon_I \).

Parameters:

strain – Principal/total strain provided in Voigt notation ( \( \boldsymbol{\varepsilon} \)), 6 components.

Returns:

A 3-component vector ( \( \partial \varepsilon^{vol}_{\ominus} / \partial \varepsilon_I \)), where each entry is the partial derivative of the corresponding volumetric-compressive component with respect to the principal strains \( \varepsilon_1, \varepsilon_2, \varepsilon_3 \).

Matrix6d dEp_dE(const Matrix6d &CelInv, const Matrix6d &Cep)

Computes the derivative \( \frac{\partial \boldsymbol{\varepsilon}^{p}}{\partial \boldsymbol{\varepsilon}} \) of the Voigt-notated plastic strain vector \( \boldsymbol{\varepsilon}^{p} \) with respect to the Voigt-notated total strain vector \( \boldsymbol{\varepsilon} \).

The relation used is

\[ \frac{\partial \boldsymbol{\varepsilon}^{p}}{\partial \boldsymbol{\varepsilon}} \;=\; \boldsymbol{I} \;-\; \mathbb{C}^{-1}\,\mathbb{C}^{(ep)} \]

where \( \mathbb{C}^{-1} \) is the elastic compliance (6×6) and \( \mathbb{C}^{(ep)} \) the elastoplastic stiffness (6×6).

Parameters:

• CelInv – Elastic compliance tensor \( \mathbb{C}^{-1} \) (6×6 matrix).

• Cep – Elastoplastic stiffness tensor \( \mathbb{C}^{(ep)} \) (6×6 matrix).

Returns:

6×6 matrix \( \partial \boldsymbol{\varepsilon}^{p} / \partial \boldsymbol{\varepsilon} \) in Voigt form.

RowVector6d dDeltaEpv_dE(const Matrix6d &CelInv, const Matrix6d &Cep)

Computes the derivative \( \frac{\partial \Delta \varepsilon^{p,vol}}{\partial \boldsymbol{\varepsilon}} \) of the volumetric plastic strain increment \( \Delta \varepsilon^{p,vol} \) with respect to the Voigt-notated total strain vector \( \boldsymbol{\varepsilon} \).

Parameters:

• CelInv – Elastic compliance tensor \( \mathbb{C}^{-1} \) (6×6 matrix).

• Cep – Elastoplastic stiffness tensor \( \mathbb{C}^{(ep)} \) (6×6 matrix).

Returns:

Row vector (1×6) containing \( \partial \Delta \varepsilon^{p,vol} / \partial \boldsymbol{\varepsilon} \) in Voigt notation.

Marmot::Matrix36 dSortedStrainPrincipal_dStrain(const Marmot::Vector6d &dEp)

Computes the derivative \( \frac{\partial \varepsilon_I}{\partial \boldsymbol{\varepsilon}} \) of the principal strains \( \varepsilon_I \) with respect to the Voigt-notated strain vector \(\boldsymbol{\varepsilon} \).

Parameters:

dEp – Input strain vector in Voigt notation (6 components).

Returns:

A 3×6 matrix where each row contains the partial derivatives of one principal strain \( \varepsilon_1, \varepsilon_2, \varepsilon_3 \) with respect to the 6 Voigt strain components.

RowVector6d dDeltaEpvneg_dE(const Marmot::Vector6d &dEp, const Matrix6d &CelInv, const Matrix6d &Cep)

Computes the derivative \( \frac{\partial \Delta \varepsilon^{p,vol}_{\ominus}} {\partial \boldsymbol{\varepsilon}} \) of the volumetric plastic strain increment in compression \( \Delta \varepsilon^{p,vol}_{\ominus} \) with respect to the strain vector in Voigt notation \( \boldsymbol{\varepsilon} \).

Parameters:

• dEp – Incremental plastic strain vector in Voigt notation ( \(\Delta\boldsymbol{\varepsilon^{p}} \)), 6 components.

• CelInv – Elastic compliance tensor \( \mathbb{C}^{-1} \) (6×6 matrix).

• Cep – Elastoplastic stiffness tensor \( \mathbb{C}^{(ep)} \) (6×6 matrix).

Returns:

Row vector (1×6) containing \( \partial \Delta \varepsilon^{p,vol}_{\ominus} / \partial \boldsymbol{\varepsilon} \) in Voigt notation.

namespace Invariants

Functions

Eigen::Vector3d principalStrains(const Marmot::Vector6d &strain)

Computes the principal strains by solving the eigenvalue problem.

This function computes the eigenvalues of the strain tensor reconstructed from the given Voigt strain vector. The eigenvalues correspond to the principal strains \( \varepsilon_1, \varepsilon_2, \varepsilon_3 \):

\[ \displaystyle |\varepsilon_{ij} - \lambda\, \delta_{ij}| = 0 \hspace{.5cm} \Rightarrow \hspace{.5cm} \lambda^{(1)},\lambda^{(2)},\lambda^{(3)}\hspace{0.3cm} \widehat{=}\hspace{0.3cm} \varepsilon_1,\, \varepsilon_2,\, \varepsilon_3 \]

Parameters:

strain – 6-component strain vector in Voigt notation.

Returns:

A 3D vector containing the principal strains (unsorted).

Eigen::Vector3d principalStresses(const Marmot::Vector6d &stress)

Computes the principal stresses by solving the eigenvalue problem.

This function computes the eigenvalues of the stress tensor reconstructed from the given Voigt stress vector. The eigenvalues correspond to the principal stresses \( \sigma_1, \sigma_2, \sigma_3 \):

\[ \displaystyle |\sigma_{ij} - \lambda\, \delta_{ij}| = 0 \hspace{.5cm} \Rightarrow \hspace{.5cm} \lambda^{(1)},\lambda^{(2)},\lambda^{(3)}\hspace{0.3cm} \widehat{=}\hspace{0.3cm}\sigma_1,\, \sigma_2,\, \sigma_3 \]

Parameters:

stress – 6-component stress vector in Voigt notation.

Returns:

A 3D vector containing the principal stresses (unsorted).

Eigen::Vector3d sortedPrincipalStrains(const Marmot::Vector6d &strain)

Computes the principal strains from Haigh–Westergaard coordinates.

This function transforms the Haigh–Westergaard coordinates \( (\xi, \rho, \theta) \) into the corresponding principal strains \( \varepsilon_1, \varepsilon_2, \varepsilon_3 \) using:

\[\begin{split} \begin{bmatrix} \varepsilon_{1} \\ \varepsilon_{2} \\ \varepsilon_{3} \end{bmatrix} = \frac{1}{\sqrt{3}} \begin{bmatrix} \xi \\ \xi \\ \xi \end{bmatrix} + \sqrt{\tfrac{2}{3}}\, \rho \, \begin{bmatrix} \cos(\theta) \\ -\sin\!\left(\tfrac{\pi}{6} - \theta\right) \\ -\sin\!\left(\tfrac{\pi}{6} + \theta\right) \end{bmatrix} \end{split}\]

Note

The computation of \( \xi \), \( \rho \), and \( \theta \) is performed in haighWestergaardFromStrain().

Parameters:

strain – 6-component strain vector in Voigt notation.

Returns:

A 3D vector of the principal strains, sorted by magnitude.

Eigen::Matrix3d principalStressesDirections(const Marmot::Vector6d &stress)

Computes the principal stress directions corresponding to the principal stresses.

This function solves the eigenvalue problem

\[ \left( \boldsymbol{\sigma} - \sigma_k \cdot \boldsymbol{I} \right) \cdot \boldsymbol{x}^{(k)} = 0 \]

to obtain the eigenvectors \( \boldsymbol{x}^{(k)} \) associated with the principal stresses \( \sigma_k \).

Parameters:

stress – 6-component stress vector in Voigt notation.

Returns:

A 3x3 matrix whose columns are the principal stress directions (unsorted).

double vonMisesEquivalentStress(const Marmot::Vector6d &stress)

Computes the equivalent von Mises stress.

Defined as:

\[ \sigma_\text{eq} = \sqrt{3 \cdot J_2} \]

where \( J_2 \) is the second invariant of the deviatoric stress tensor (see J2()).

Parameters:

stress – 6-component stress vector in Voigt notation.

Returns:

The von Mises equivalent stress.

double vonMisesEquivalentStrain(const Marmot::Vector6d &strain)

Computes the equivalent von Mises strain.

Based on the deviatoric strain tensor \( e_{ij} \):

\[ \varepsilon^\text{(eq)} = \sqrt{\tfrac{2}{3} \, e_{ij} e_{ij}} \]

Parameters:

strain – 6-component strain vector in Voigt notation.

Returns:

The von Mises equivalent strain.

template<typename T>
T normStrain(const Eigen::Matrix<T, 6, 1> &strain)

Computes the Euclidean norm of the strain tensor.

Defined as:

\[ ||\boldsymbol{\varepsilon}|| = \sqrt{ \varepsilon_{ij}\varepsilon_{ij}} \]

Template Parameters:

T – Scalar type (e.g., double, float).

Parameters:

strain – 6-component strain vector in Voigt notation.

Returns:

The Euclidean norm of the strain tensor.

double normStress(const Marmot::Vector6d &stress)

Computes the Euclidean norm of the stress tensor.

Defined as:

\[ ||\boldsymbol{\sigma}|| = \sqrt{\sigma_{ij}\sigma_{ij}} \]

Parameters:

stress – 6-component stress vector in Voigt notation.

Returns:

The Euclidean norm of the stress tensor.

double StrainVolumetricNegative(const Marmot::Vector6d &strain)

Computes the volumetric plastic strain in compression.

Defined as:

\[ \varepsilon^{\text{vol}}_{\ominus} = \sum_{i=1}^3 \left\langle -\varepsilon_i \right\rangle \]

where \( \varepsilon_i \) are the principal strains and \( \langle \bullet \rangle \) denotes the Macaulay brackets.

Parameters:

strain – 6-component strain vector in Voigt notation.

Returns:

The volumetric compressive strain.

template<typename T>
T I1(const Eigen::Matrix<T, 6, 1> &stress)

Computes the first invariant \( I_1 \) of the stress tensor.

The first invariant of the stress tensor \( \boldsymbol{\sigma} \) is the trace:

\[ I_1 = \operatorname{tr}(\boldsymbol{\sigma}) = \sigma_{11} + \sigma_{22} + \sigma_{33} \]

Template Parameters:

T – Scalar type (e.g., double, float).

Parameters:

stress – 6-component stress vector in Voigt notation: \( \boldsymbol{\sigma}^{\text{Voigt}} = \begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{33} & \sigma_{12} & \sigma_{13} & \sigma_{23} \end{bmatrix}^T \)

Returns:

Value of the first invariant \( I_1 \).

double I1Strain(const Marmot::Vector6d &strain)

Computes the first invariant \( I^{(\varepsilon)}_1 \) of the strain tensor.

The first invariant of the strain tensor \( \boldsymbol{\varepsilon} \) is the trace:

\[ I^{(\varepsilon)}_1 = \operatorname{tr}(\boldsymbol{\varepsilon}) = \varepsilon_{11} + \varepsilon_{22} + \varepsilon_{33} \]

Parameters:

strain – 6-component strain vector in Voigt notation: \( \boldsymbol{\varepsilon}^{\text{Voigt}} = \begin{bmatrix} \varepsilon_{11} & \varepsilon_{22} & \varepsilon_{33} & \gamma_{12} & \gamma_{13} & \gamma_{23} \end{bmatrix}^T \)

Returns:

Value of the first strain invariant \( I^{(\varepsilon)}_1 \).

template<typename T>
T I2(const Eigen::Matrix<T, 6, 1> &stress)

Computes the second invariant \( I_2 \) of the stress tensor.

The second invariant of the stress tensor \( \boldsymbol{\sigma} \) is defined as:

\[ I_2 = \tfrac{1}{2}\left( \operatorname{tr}(\boldsymbol{\sigma})^2 - \operatorname{tr}(\boldsymbol{\sigma}^2) \right) \]

In Voigt notation, this expands to:

\[ I_2 = \sigma_{11}\sigma_{22} + \sigma_{22}\sigma_{33} + \sigma_{33}\sigma_{11} - \sigma_{12}^2 - \sigma_{13}^2 - \sigma_{23}^2 \]

Template Parameters:

T – Scalar type (e.g., double, float).

Parameters:

stress – 6-component stress vector in Voigt notation: \( \boldsymbol{\sigma}^{\text{Voigt}} = \begin{bmatrix} \sigma_{11} & \sigma_{22} & \sigma_{33} & \sigma_{12} & \sigma_{13} & \sigma_{23} \end{bmatrix}^T \)

Returns:

Value of the second invariant \( I_2 \).

double I2Strain(const Marmot::Vector6d &strain)

Computes the second invariant \( I^{(\varepsilon)}_2 \) of the strain tensor.

For a strain tensor \( \boldsymbol{\varepsilon} \) in Voigt notation, the second invariant is:

\[ I^{(\varepsilon)}_2 = \varepsilon_{11}\varepsilon_{22} + \varepsilon_{22}\varepsilon_{33} + \varepsilon_{33}\varepsilon_{11} - \tfrac{1}{4}\left(\gamma_{12}^2 + \gamma_{13}^2 + \gamma_{23}^2\right) \]

Parameters:

strain – 6-component strain vector in Voigt notation: \( \boldsymbol{\varepsilon}^{\text{Voigt}} = [ \varepsilon_{11}, \varepsilon_{22}, \varepsilon_{33}, \gamma_{12}, \gamma_{13}, \gamma_{23} ]^T \)

Returns:

Value of the second invariant \( I^{(\varepsilon)}_2 \).

template<typename T>
T I3(const Eigen::Matrix<T, 6, 1> &stress)

Computes the third invariant \( I_3 \) of the stress tensor.

Defined as the determinant of the Cauchy stress tensor \( \boldsymbol{\sigma} \):

\[ I_3 = \det(\boldsymbol{\sigma}) \]

Expanded in Voigt notation:

\[ I_3 = \sigma_{11}\sigma_{22}\sigma_{33} + 2\sigma_{12}\sigma_{13}\sigma_{23} - \sigma_{11}\sigma_{23}^2 - \sigma_{22}\sigma_{13}^2 - \sigma_{33}\sigma_{12}^2 \]

Template Parameters:

T – Scalar type (e.g., double, float).

Parameters:

stress – 6-component stress vector in Voigt notation.

Returns:

Value of the third invariant \( I_3 \).

double I3Strain(const Marmot::Vector6d &strain)

Computes the third invariant \( I^{(\varepsilon)}_3 \) of the strain tensor.

The third invariant is given by the determinant of the strain tensor:

\[ I^{(\varepsilon)}_3 = \det(\boldsymbol{\varepsilon}) \]

Parameters:

strain – 6-component strain vector in Voigt notation.

Returns:

Value of the third strain invariant \( I^{(\varepsilon)}_3 \).

template<typename T>
T J2(const Eigen::Matrix<T, 6, 1> &stress)

Computes the second invariant \( J_2 \) of the deviatoric stress tensor.

For the deviatoric stress tensor \( \boldsymbol{s} \):

\[ J_2 = \tfrac{1}{3} I_1^2 - I_2 \]

Template Parameters:

T – Scalar type (e.g., double, float).

Parameters:

stress – 6-component stress vector in Voigt notation.

Returns:

Value of the second deviatoric invariant \( J_2 \) (clamped to non-negative).

double J2Strain(const Marmot::Vector6d &strain)

Computes the second invariant \( J^{(\varepsilon)}_2 \) of the deviatoric strain tensor.

For the deviatoric strain tensor \( \boldsymbol{e} \):

\[ J^{(\varepsilon)}_2 = \tfrac{1}{3} \left(I^{(\varepsilon)}_1\right)^2 - I^{(\varepsilon)}_2 \]

Parameters:

strain – 6-component strain vector in Voigt notation.

Returns:

Value of the second deviatoric strain invariant \( J^{(\varepsilon)}_2 \).

template<typename T>
T J3(const Eigen::Matrix<T, 6, 1> &stress)

Computes the third invariant \( J_3 \) of the deviatoric stress tensor.

For the deviatoric stress tensor \( \boldsymbol{s} \):

\[ J_3 = \tfrac{2}{27} I_1^3 - \tfrac{1}{3} I_1 I_2 + I_3 \]

Template Parameters:

T – Scalar type (e.g., double, float).

Parameters:

stress – 6-component stress vector in Voigt notation.

Returns:

Value of the third deviatoric stress invariant \( J_3 \).

double J3Strain(const Marmot::Vector6d &strain)

Computes the third invariant \( J^{(\varepsilon)}_3 \) of the deviatoric strain tensor \(\boldsymbol{e} \).

The third invariant is obtained by reconstructing the strain tensor from Voigt notation and computing its determinant:

\[ J^{(\varepsilon)}_3 = \det(\boldsymbol{e}) \]

Parameters:

strain – 6-component strain vector in Voigt notation.

Returns:

Value of the third deviatoric strain invariant \( J^{(\varepsilon)}_3 \).

std::pair<Eigen::Vector3d, Eigen::Matrix<double, 3, 6>> principalValuesAndDerivatives(const Eigen::Matrix<double, 6, 1> &S)

Computes principal values and their derivatives for a symmetric 3×3 matrix in Voigt notation.

This routine implements a fast algorithm to compute the principal values (eigenvalues) of a symmetric second-order tensor in Voigt notation (off-diagonal entries are expected without the factor of 2).

Parameters:

S – 6-component symmetric matrix in Voigt notation.

Returns:

A pair consisting of:

• A vector containing the (unsorted) principal values.

• An \( 3 \times 6 \) matrix containing the derivatives of the principal values with respect to the Voigt components of \( S \).

namespace Transformations

Functions

Matrix6d transformationMatrixStrainVoigt(const Matrix3d &transformedCoordinateSystem)

Computes the transformation matrix \( T_{\varepsilon} \) to rotate a strain vector in Voigt notation into another Cartesian coordinate system.

This function modifies the stress transformation matrix \(T_{\varepsilon}\) to account for the engineering shear convention used for strains in Voigt notation. Specifically:

• The upper-right 3×3 block is scaled by 0.5

• The lower-left 3×3 block is scaled by 2

This function builds the 6×6 transformation matrix from the direction cosines of the new basis:

\[ \boldsymbol{\varepsilon}'_{\text{voigt}} \;=\; T_{\varepsilon}\,\boldsymbol{\varepsilon}_{\text{voigt}}. \]

which is equivalent to the tensor transformation

\[ \boldsymbol{\varepsilon}' \;=\; \boldsymbol{Q}\,\boldsymbol{\varepsilon}\,\boldsymbol{Q}^{T} \]

Parameters:

transformedCoordinateSystem – 3×3 rotation matrix \(\boldsymbol{Q}\) that maps from the original to the target Cartesian basis.

Returns:

6×6 transformation matrix \(T_{\varepsilon}\) that transforms strain vectors in Voigt notation.

Matrix6d transformationMatrixStressVoigt(const Matrix3d &transformedCoordinateSystem)

Computes the transformation matrix \( T_{\sigma} \) to rotate a stress vector in Voigt notation into another Cartesian coordinate system.

This function builds the 6×6 transformation matrix from the direction cosines of the new basis:

\[ \boldsymbol{\sigma}'_{\text{voigt}} \;=\; T_{\sigma}\,\boldsymbol{\sigma}_{\text{voigt}}, \]

which is equivalent to the tensor transformation

\[ \boldsymbol{\sigma}' \;=\; \boldsymbol{Q}\,\boldsymbol{\sigma}\,\boldsymbol{Q}^{T}. \]

Parameters:

transformedCoordinateSystem – 3×3 rotation matrix \(\boldsymbol{Q}\) that maps from the original to the target Cartesian basis.

Returns:

6×6 transformation matrix \(T_{\sigma}\) that transforms stress vectors in Voigt notation.

Matrix36d projectVoigtStressToPlane(const Vector3d &normalVector)

Calculates the projection matrix for mapping Voigt stress to a traction vector.

This function returns the \(3 \times 6\) projection matrix \(\mathbf{P}^{(\mathbf{n})}\) that maps a Voigt-notated stress vector \(\boldsymbol{\sigma}_{\text{voigt}} \in \mathbb{R}^6\) to the traction vector \(\mathbf{t}^{(\mathbf{n})} \in \mathbb{R}^3\) on a plane with an outward unit normal vector \(\mathbf{n}\). The relationship is given by:

\[\mathbf{t}^{(\mathbf{n})} = \mathbf{P}^{(\mathbf{n})} \, \boldsymbol{\sigma}_{\text{voigt}}. \]

This expression is mathematically equivalent to Cauchy’s formula \( \boldsymbol{t}^{(n)}_i = \sigma_{ji} \, n_j \).

Parameters:

normalVector – The \(3 \times 1\) unit vector \(\mathbf{n}\) representing the normal of the plane.

Returns:

The \(3 \times 6\) projection matrix \(\mathbf{P}^{(\mathbf{n})}\).

Matrix36d projectVoigtStrainToPlane(const Vector3d &normalVector)

Calculates the projection matrix for mapping a Voigt strain vector to the strain on a plane.

This function is similar to projectVoigtStressToPlane but for strains in voigt notation. For more details, see the documentation of projectVoigtStressToPlane.

Parameters:

normalVector – The \(3 \times 1\) unit vector \(\mathbf{n}\) representing the normal of the plane.

Returns:

The \(3 \times 6\) projection matrix \(\mathbf{P}^{(\mathbf{n})}\).

Marmot::Vector6d rotateVoigtStress(const Eigen::Matrix3d &Q, const Marmot::Vector6d &stress)

Rotates a stress tensor in Voigt notation.

This function rotates a stress tensor \( \boldsymbol{\sigma} \) using a rotation matrix \( \boldsymbol{Q} \), while keeping the tensor in Voigt notation. The rotated stress \( \boldsymbol{\sigma}^{\prime} \) is computed as:

\[\boldsymbol{\sigma}^{\prime} = \boldsymbol{Q} \, \boldsymbol{\sigma} \, \boldsymbol{Q}^{T}. \]

Parameters:

• Q – A \(3 \times 3\) rotation matrix \( \boldsymbol{Q} \).

• stress – A Voigt-notated stress vector \( \boldsymbol{\sigma} \in \mathbb{R}^6 \).

Returns:

The rotated Voigt-notated stress vector \( \boldsymbol{\sigma}^{\prime} \in \mathbb{R}^6 \).

Marmot::Vector6d transformStressToLocalSystem(const Marmot::Vector6d &stress, const Matrix3d &transformedCoordinateSystem)

Transforms a stress tensor in Voigt notation from the global to a local coordinate system.

The transformation convention used is: \( \boldsymbol{\sigma}_{\text{local}} = Q^T \, \boldsymbol{\sigma}_{\text{global}} \, Q \), where \( Q \) is the transformation matrix.

Parameters:

• stress – The stress tensor in Voigt notation ( \( \mathbb{R}^6 \)) in the global system.

• transformedCoordinateSystem – A \(3 \times 3\) orthonormal matrix whose columns are the axes of the local coordinate system expressed in the global system. Typically, this is a rotation matrix.

Returns:

The stress tensor in Voigt notation ( \( \mathbb{R}^6 \)) in the local system.

Marmot::Vector6d transformStrainToLocalSystem(const Marmot::Vector6d &stress, const Matrix3d &transformedCoordinateSystem)

Transforms a strain tensor in Voigt notation from the global to a local coordinate system.

The transformation convention used is: \( \boldsymbol{\varepsilon}_{\text{local}} = Q^T \, \boldsymbol{\varepsilon}_{\text{global}} \, Q \), where \( Q \) is the transformation matrix.

Parameters:

• stress – The strain tensor in Voigt notation ( \( \mathbb{R}^6 \)) in the global system.

• transformedCoordinateSystem – A \(3 \times 3\) orthonormal matrix whose columns are the axes of the local coordinate system expressed in the global system. Typically, this is a rotation matrix.

Returns:

The strain tensor in Voigt notation ( \( \mathbb{R}^6 \)) in the local system.

Marmot::Vector6d transformStressToGlobalSystem(const Marmot::Vector6d &stress, const Matrix3d &transformedCoordinateSystem)

Transforms a stress tensor in Voigt notation from a local to the global coordinate system.

The transformation convention used is: \( \boldsymbol{\sigma}_{\text{global}} = Q \, \boldsymbol{\sigma}_{\text{local}} \, Q^T \), where \( Q \) is the transformation matrix.

Parameters:

• stress – The stress tensor in Voigt notation ( \( \mathbb{R}^6 \)) in the local system.

• transformedCoordinateSystem – A \(3 \times 3\) orthonormal matrix whose columns are the axes of the local coordinate system expressed in the global system. Typically, this is a rotation matrix.

Returns:

The stress tensor in Voigt notation ( \( \mathbb{R}^6 \)) in the global system.

Marmot::Vector6d transformStrainToGlobalSystem(const Marmot::Vector6d &stress, const Matrix3d &transformedCoordinateSystem)

Transforms a strain tensor in Voigt notation from a local to the global coordinate system.

The transformation convention used is: \( \boldsymbol{\varepsilon}_{\text{global}} = Q \, \boldsymbol{\varepsilon}_{\text{local}} \, Q^T \), where \( Q \) is the transformation matrix.

Parameters:

• stress – The strain tensor in Voigt notation ( \( \mathbb{R}^6 \)) in the local system.

• transformedCoordinateSystem – A \(3 \times 3\) orthonormal matrix whose columns are the axes of the local coordinate system expressed in the global system. Typically, this is a rotation matrix.

Returns:

The strain tensor in Voigt notation ( \( \mathbb{R}^6 \)) in the global system.

Marmot::Matrix6d transformStiffnessToGlobalSystem(const Marmot::Matrix6d &stiffness, const Matrix3d &transformedCoordinateSystem)

Transforms a stiffness tensor in Voigt notation from a local to the global coordinate system.

The transformation convention used is: \( \mathbf{C}_{\text{global}} = T \, \mathbf{C}_{\text{local}} \, T^T \), where \( T \) is the transformation matrix derived from \( Q \).

Parameters:

• stiffness – The stiffness tensor in Voigt notation ( \( \mathbb{R}^{6 \times 6} \)) in the local system.

• transformedCoordinateSystem – A \(3 \times 3\) orthonormal matrix whose columns are the axes of the local coordinate system expressed in the global system. Typically, this is a rotation matrix.

Returns:

The stiffness tensor in Voigt notation ( \( \mathbb{R}^{6 \times 6} \)) in the global system.

namespace GeostaticStress

Functions

std::tuple<double, double, double> getGeostaticStressFromLinearDistribution(const double *geostaticStressDefintion, double coordinate_y)

namespace Materials

namespace KelvinChain

Typedefs

typedef Eigen::VectorXd Properties

Vector of material properties.

This namespace contains the functions that derive the equivalent Kelvin chain for a given compliance function

Convenience typedef for an Eigen dynamic-size vector (Eigen::VectorXd) that holds the material or model parameters of the Kelvin chain.

typedef Eigen::Map<Properties> mapProperties

Mapped view of material properties.

Alias for Eigen::Map<Properties>, which allows mapping an existing contiguous memory block as a Properties vector without copying.

typedef Eigen::Matrix<double, 6, Eigen::Dynamic> StateVarMatrix

Matrix of state variables.

Alias for an Eigen matrix of shape 6 × N (Eigen::Matrix<double, 6, Eigen::Dynamic>), storing viscoelastic strain variables for each Kelvin unit.

typedef Eigen::Map<StateVarMatrix> mapStateVarMatrix

Mapped view of the state variable matrix.

Alias for Eigen::Map<StateVarMatrix>, allowing access to an existing state variable array as an Eigen matrix without copying.

Functions

template<int k>
double evaluatePostWidderFormula(std::function<autodiff::Real<k, double>(autodiff::Real<k, double>)> phi, double tau)

Evaluates the Post–Widder inversion formula to approximate the discrete retardance function \(L_k(\tau)\).

Evaluates the Post-Widder formula for deriving the equivalent retardance \(L_\kappa(\tau)\) of the discrete Kelvin Chain model, given the continuous description of the compliance

\[\Phi(t) = \int_{\tau=0}^\infty L(\tau)\left(1-e^{-\frac{t}{\tau}}\right)d(\ln\tau)\]

according to the formula

\[L_k(\tau)=(-1)^{k-1}\frac{(-tk)^k\Phi^{(k)(t)}}{(k-1)!}\]

in the sense that \(L_k(\tau)=\lim_{\kappa\to\infty}L_\kappa(\tau)\) in the frequency domain.

Template Parameters:

k – \(k\) integer defining the order of differentiation that corresponds to the \(L_k\) term in the sequence of the retardance approximation

Parameters:

• phi – [in] compliance function whose discrete approximation is sought

• tau – [in] the continuous variable in the frequency domain

Returns:

val the value of the \(L_k(\tau)\) approximation.

template<int k>
double approximateZerothCompliance(std::function<autodiff::Real<k, double>(autodiff::Real<k, double>)> phi, double tauMin, double spacing = 10.)

Computes the zeroth-order (instantaneous) elastic compliance by numerically integrating the retardance function.

Approximates the instantaneous elastic compliance \(\Phi_0=\frac{1}{E_0}\),

\[\frac{1}{E_0}=\int_{\tau=0}^{\tau}L(\tau)d\ln(\tau)\]

It is used in order to avoid integrating over the origin of the frequency domain where \(\int_0^{\tau_0}d\ln(\tau)\) is singular.

Note

The integral is evaluated numerically using Simpson’s rule via NumericalAlgorithms::Integration::integrateScalarFunction. The integration range \([\tau_{\min}/\sqrt{\text{spacing}},\,1\times10^{-14}]\) is chosen to avoid the logarithmic singularity at \(\tau=0\). As Simpson’s rule is sensitive to changes in the sign of the summands, integration should be handled with care. A Gauss quadrature approach involving Chebyshev polynomials could be implemented in the future.

Template Parameters:

k – \(k\) integer defining the order of differentiation for the \(k-th\) approximation of the retardance.

Parameters:

• phi – [in] compliance function whose discrete approximation is sought.

• tauMin – [in] the minimum retardation time that is required to be of the same order or smaller than the total time of the analysis.

• spacing – [in] the increment between two sunsequent retardation times

Returns:

val the value of the zeroth compliance term.

template<int k>
Properties computeElasticModuli(std::function<autodiff::Real<k, double>(autodiff::Real<k, double>)> phi, Properties retardationTimes, bool gaussQuadrature = false)

Computes the discrete elastic moduli of the equivalent Kelvin chain from the \(L_k\) retardance approximation.

Evaluates the elastic moduli based on the \(L_k\) approximation of the retardance,

\[ \frac{1}{E_\mu}=(\ln 10)L_k(\tau_\mu),\;\mu=1,2,...,M\]

Template Parameters:

k – \(k\) integer defining the order of differentiation for the \(k-th\) approximation of the retardance.

Parameters:

• phi – [in] compliance function whose discrete approximation is sought

• retardationTimes – [in] the minimum retardation time that is required to be of the same order or smaller than the total time of the analysis

• gaussQuadrature – [in] flag if on Gauss Quadrature is performed else the Post-Widder formula is applied.

Returns:

elasticModuli the discrete elastic moduli of the equivalt Kelvin chain.

Properties generateRetardationTimes(int n, double min, double spacing)

Generates a sequence of logarithmically spaced retardation times for the Kelvin chain model.

This function evaluates the retardation times.

Parameters:

• n – [in] the number of kelvin units in the equivalent Kelvin chain.

• min – [in] the first retardation time.

• spacing – [in] the spacing between the retardation times.

Returns:

retardationTimes a vector containing the retardation time for each Kelvin unit in the Kelvin chain.

void updateStateVarMatrix(const double dT, Properties elasticModuli, Properties retardationTimes, Eigen::Ref<StateVarMatrix> stateVars, const Marmot::Vector6d &dStress, const Marmot::Matrix6d &unitComplianceMatrix)

Updates the viscoelastic strain state variables for each Kelvin unit over a given time increment.

For the given time increment \(\Delta t_k\) this function updates the visco-elastic state variables according to the update rule:

\[\varepsilon^{ev,\mu,k+1}_{ij} = \frac{\lambda^{\mu k}}{E_\mu}C^\nu_{ijkl}\Delta\sigma_{kl}+\beta^{\mu ,k}\varepsilon^{ev,\mu, k}_{ij}\]

Parameters:

• dT – [in] the time increment.

• elasticModuli – [in] vector containing the elastic moduli of ech Kelvin unit in the Kelvin chain.

• retardationTimes – [in] vector containing the retardation time for each Kelvin unit in the Kelvin chain.

• stateVars – [inout] the \([6\times \mu]\) matrix that contains the viscoelastic strain update for each unit of the Kelvin chain.

• dStress – [in] the \([6\times 1]\) vector of the total stress increment.

• unitComplianceMatrix – [in] the [6\times 6] compliance matrix of the material with unit compliance and given Poisson coeffiscient.

void evaluateKelvinChain(const double dT, Properties elasticModuli, Properties retardationTimes, StateVarMatrix stateVars, double &uniaxialCompliance, Marmot::Vector6d &dStrain, const double factor)

Evaluates the viscoelastic response of the Kelvin chain over a time increment.

For the given time increment \(\Delta t_k\) this function updates the visco-elastic state variables according to the update rule:

\[\overline{J}^k=\sum_{\mu=1}^M\frac{1-\lambda^{\mu,k}}{E^\mu}\]

\[\Delta\varepsilon^{'',k}_{ij} = \sum_{\mu=1}^M\left(1-\beta^{\mu,k}\right)\varepsilon^{ev,\kappa}_{ij}\]

Parameters:

• dT – [in] the time increment.

• elasticModuli – [in] vector containing the elastic moduli of ech Kelvin unit in the Kelvin chain.

• retardationTimes – [in] vector containing the retardation time for each Kelvin unit in the Kelvin chain.

• stateVars – [in] the \([6\times \mu]\) matrix that contains the viscoelastic strain update for each unit of the Kelvin chain.

• uniaxialCompliance – [inout] number containing the \(\sum^M_{\mu=1}\frac{\lambda^{\mu, k}}{E_\mu}\).

• dStrain – [inout] the viscoelastic update of the strain based on the state variables of the Kelvin chain.

• factor – [in] the solidification factor (in non aging viscoelasticity set to 1).

void computeLambdaAndBeta(double dT, double tau, double &lambda, double &beta)

Computes the time-dependent relaxation factors \(\lambda\) and \(\beta\) for a given Kelvin unit.

This function evaluates the time parameters for each Kelvin chain.

Parameters:

• dT – [in] the time increment.

• tau – [in] the retardation time for a given Kelvin unit.

• lambda – [out] time dependent factor for each Kelvin unit.

• beta – [out] time dependemnt factor for each Kelvin unit.

namespace NumericalAlgorithms

namespace TimeIntegration

Functions

template<int nDim>
void newmarkBetaIntegration(const double *du, double *v, double *a, double dT, double newmarkBeta, double newmarkGamma, double *da_ddu)

Newmark-Beta time integration.

This function implements the Newmark-Beta time integration method for linear dynamic problems. It updates the velocity and acceleration vectors based on the displacement increment and the time step.

Note

The function assumes that the input vectors are of size nDim.

Note

The function also assumes that the time step is positive and that the Newmark-Beta parameter is not zero.

Note

The function uses a small value (1e-16) to avoid division by zero

Parameters:

• nDim – Number of dimensions (2D or 3D)

• du – Displacement increment vector

• v – Velocity vector

• a – Acceleration vector

• dT – Time step

• newmarkBeta – Newmark-Beta parameter (0.25 for linear problems)

• newmarkGamma – Newmark-Gamma parameter (0.5 for linear problems)

• da_ddu – Derivative of acceleration with respect to displacement

Returns:

void

file AdaptiveSubstepper.h

#include “Marmot/MarmotJournal.h”

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for AdaptiveSubstepper.h:

file AdaptiveSubstepperExplicit.h

#include “Marmot/MarmotJournal.h”

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for AdaptiveSubstepperExplicit.h:

file DuvautLionsViscosity.h

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for DuvautLionsViscosity.h:

file HaighWestergaard.h

#include “Marmot/MarmotMath.h”

#include “Marmot/MarmotTypedefs.h”

#include “Marmot/MarmotVoigt.h”

Include dependency graph for HaighWestergaard.h:

This graph shows which files directly or indirectly include HaighWestergaard.h:

file HughesWinget.h

#include “Marmot/MarmotVoigt.h”

Include dependency graph for HughesWinget.h:

file MarmotElasticity.h

#include “Marmot/MarmotVoigt.h”

Include dependency graph for MarmotElasticity.h:

file MarmotGeostaticStress.h

#include “Marmot/MarmotJournal.h”

#include <Eigen/Core>

#include <functional>

Include dependency graph for MarmotGeostaticStress.h:

file MarmotKelvinChain.h

#include “Marmot/MarmotNumericalIntegration.h”

#include “Marmot/MarmotTypedefs.h”

#include “autodiff/forward/real.hpp”

#include <functional>

Include dependency graph for MarmotKelvinChain.h:

file MarmotKinematics.h

#include “Marmot/MarmotJournal.h”

#include “Marmot/MarmotTypedefs.h”

#include “Marmot/MarmotVoigt.h”

Include dependency graph for MarmotKinematics.h:

file MarmotLocalization.h

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for MarmotLocalization.h:

file MarmotLowerDimensionalStress.h

#include “Marmot/MarmotVoigt.h”

Include dependency graph for MarmotLowerDimensionalStress.h:

file MarmotMaterialGradientEnhancedHypoElastic.h

#include “Marmot/MarmotMaterialGradientEnhancedMechanical.h”

Include dependency graph for MarmotMaterialGradientEnhancedHypoElastic.h:

file MarmotMaterialGradientEnhancedMechanical.h

#include “Marmot/MarmotMaterial.h”

Include dependency graph for MarmotMaterialGradientEnhancedMechanical.h:

This graph shows which files directly or indirectly include MarmotMaterialGradientEnhancedMechanical.h:

file MarmotMaterialHyperElastic.h

#include “Marmot/MarmotMaterialMechanical.h”

Include dependency graph for MarmotMaterialHyperElastic.h:

file MarmotMaterialHypoElastic.h

#include “Marmot/MarmotMaterialMechanical.h”

Include dependency graph for MarmotMaterialHypoElastic.h:

This graph shows which files directly or indirectly include MarmotMaterialHypoElastic.h:

file MarmotMaterialHypoElasticAD.h

#include “Marmot/MarmotMaterialHypoElastic.h”

#include “autodiff/forward/dual.hpp”

Include dependency graph for MarmotMaterialHypoElasticAD.h:

file MarmotMaterialMechanical.h

#include “Marmot/MarmotMaterial.h”

Include dependency graph for MarmotMaterialMechanical.h:

This graph shows which files directly or indirectly include MarmotMaterialMechanical.h:

file MarmotMaterialPointSolverHypoElastic.h

#include “Marmot/MarmotMaterialHypoElastic.h”

#include “Marmot/MarmotTypedefs.h”

#include <iostream>

Include dependency graph for MarmotMaterialPointSolverHypoElastic.h:

file MarmotPronySeries.h

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for MarmotPronySeries.h:

file MarmotStateVarVectorManager.h

#include “Marmot/MarmotUtils.h”

#include <memory>

#include <string>

#include <unordered_map>

#include <vector>

Include dependency graph for MarmotStateVarVectorManager.h:

file MarmotUtility.h

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for MarmotUtility.h:

file MarmotViscoelasticity.h

#include “Marmot/MarmotVoigt.h”

Include dependency graph for MarmotViscoelasticity.h:

file MarmotVoigt.h

#include “Marmot/MarmotJournal.h”

#include “Marmot/MarmotMath.h”

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for MarmotVoigt.h:

This graph shows which files directly or indirectly include MarmotVoigt.h:

Defines

VOIGTFROMDIM(x)

file MenetreyWillam.h

#include “Marmot/HaighWestergaard.h”

#include <utility>

Include dependency graph for MenetreyWillam.h:

file NewmarkBetaIntegrator.h

#include <algorithm>

Include dependency graph for NewmarkBetaIntegrator.h:

file PerezFougetSubstepperExplicitMarkII.h

#include “Marmot/MarmotJournal.h”

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for PerezFougetSubstepperExplicitMarkII.h:

file PerezFougetSubstepperMarkII.h

#include “Marmot/MarmotJournal.h”

#include “Marmot/MarmotMath.h”

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for PerezFougetSubstepperMarkII.h:

file PerezFougetSubstepperTime.h

#include “Marmot/MarmotTypedefs.h”

#include “Marmot/MarmotJournal.h”

Include dependency graph for PerezFougetSubstepperTime.h:

file YieldSurfaceCombinationManager.h

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for YieldSurfaceCombinationManager.h:

page Todo List

Class Marmot::ContinuumMechanics::CommonConstitutiveModels::DuvautLionsViscosity< nMatTangentSize >

: Update member names to more descriptive ones

Member Marmot::ContinuumMechanics::CommonConstitutiveModels::DuvautLionsViscosity< nMatTangentSize >::applyViscosityOnMatTangent  (const TangentSizedMatrix &matTangentInv, double dT)

: Check if application to inverse can be replaced by application to non-inverse in general

Member MarmotMaterialHyperElastic::computeStress  (double *S, double *dSdE, const double *FOld, const double *FNew, const double *timeOld, const double dT, double &pNewDT) override

A default implementation is provided.

Member MarmotMaterialHyperElastic::computeStressPK2  (double *S, double *dSdE, const double *E, const double *timeOld, const double dT, double &pNewDT)=0

Should we use function overloading in this case and simple also use computeStress for the function name?

Member MarmotMaterialHypoElastic::computeStress  (double *stress, double *dStressDDStrain, const double *FOld, const double *FNew, const double *timeOld, const double dT, double &pNewDT) override

A default implementation is provided.

dir /home/runner/work/Marmot/Marmot/modules/core

dir /home/runner/work/Marmot/Marmot/modules/core/MarmotMechanicsCore/include

dir /home/runner/work/Marmot/Marmot/modules/core/MarmotMechanicsCore/include/Marmot

dir /home/runner/work/Marmot/Marmot/modules/core/MarmotMechanicsCore

dir /home/runner/work/Marmot/Marmot/modules