MarmotMechanicsCore

template<size_t materialTangentSize, size_t nIntegrationDependentStateVars>
class AdaptiveSubstepper

#include <AdaptiveSubstepper.h>

Adaptive sub-stepper employing error estimation and Richardson extrapolation.

Splits a strain increment into sub-increments, performs a full-step and two half-steps for each sub-increment, and estimates the integration error via Richardson extrapolation. The sub-increment size is adapted automatically based on the error estimate.

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

Template Parameters:

• materialTangentSize – Size of the material tangent matrix.

• nIntegrationDependentStateVars – Number of integration-dependent state variables.

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, 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)

Construct an AdaptiveSubstepper.

Parameters:

• initialStepSize – Initial sub-step size as a fraction of the total increment (0, 1].

• minimumStepSize – Minimum allowable sub-step size; triggers a warning if reached.

• maxScaleUpFactor – Maximum factor by which the sub-step size may be increased.

• scaleDownFactor – Factor applied to the sub-step size when a sub-step is discarded.

• integrationErrorTolerance – Richardson-extrapolation error threshold for accepting a sub-step.

• nPassesToIncrease – Number of successful sub-steps before the step size is allowed to grow.

• Cel – Linear-elastic stiffness matrix (6×6).

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

#include <AdaptiveSubstepperExplicit.h>

Adaptive sub-stepper for semi-explicit (forward-Euler-based) return-mapping algorithms.

Analogous to AdaptiveSubstepper but adapted for explicit integration schemes. Error estimation is performed by Richardson extrapolation between a full-step and two half-steps.

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

Template Parameters:

• materialTangentSize – Size of the material tangent matrix.

• nIntegrationDependentStateVars – Number of integration-dependent state variables.

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

Auxiliary matrix relating material-state columns to strain columns (materialTangentSize × 6)

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

Vector to carry the internal state of a material.

Public Functions

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

Construct an AdaptiveSubstepperExplicit.

Parameters:

• initialStepSize – Initial sub-step size as a fraction of the total increment (0, 1].

• minimumStepSize – Minimum allowable sub-step size; triggers a warning if reached.

• maxScaleUpFactor – Maximum factor by which the sub-step size may be increased.

• scaleDownFactor – Factor applied to the sub-step size when a sub-step is discarded.

• integrationErrorTolerance – Richardson-extrapolation error threshold for accepting a sub-step.

• nPassesToIncrease – Number of successful sub-steps before the step size is allowed to grow.

• Cel – Linear-elastic stiffness matrix (6×6).

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

Square matrix type sized to the material tangent dimension.

Public Functions

DuvautLionsViscosity(double viscosity)

Constructs the viscosity model with the given viscosity parameter.

Parameters:

viscosity – Viscosity parameter \(\eta\) for the Duvaut-Lions regularisation.

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

enum factorial

Values:

enumerator value

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

Print this history entry to stdout.

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

#include <HughesWinget.h>

Hughes–Winget objective integration algorithm.

Given the deformation gradients at the beginning and end of an increment, this class computes:

• the symmetric strain increment in Voigt notation,

• the incremental rotation matrix,

• a rotated stress tensor, and

• consistent tangent contributions for mixed deformation-gradient/stress quantities.

Public Types

enum Formulation

Available integration formulations.

Values:

enumerator AbaqusLike

Formulation consistent with the Abaqus UMAT conventions.

Public Functions

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

Constructs a HughesWinget integrator from two successive deformation gradients.

Parameters:

• FOld – Deformation gradient at the start of the increment.

• FNew – Deformation gradient at the end of the increment.

• formulation – Integration formulation to use.

Marmot::Vector6d getStrainIncrement()

Returns the objective symmetric strain increment in Voigt notation.

Eigen::Matrix3d getRotationIncrement()

Returns the incremental rotation matrix \(\Delta\mathbf{R}\).

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

Rotates a Voigt stress (or strain) tensor with the incremental rotation.

Parameters:

tensor – Input tensor in Voigt notation.

Returns:

Rotated tensor in Voigt notation.

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

Compute the consistent tangent contribution \(\frac{d\boldsymbol{S}}{d\mathbf{F}}\).

Given the current Cauchy stress, the inverse deformation gradient, and the Cauchy-stress/strain tangent, returns the third-order tensor mapping variations of the deformation gradient to variations of the Cauchy stress.

Parameters:

• stress – Current Cauchy stress in Voigt notation.

• FInv – Inverse of the deformation gradient at the end of the increment.

• dChauchyDEps – Algorithmic Cauchy-stress/strain tangent (6×6).

Returns:

The 6×3×3 tangent tensor \(\frac{d\boldsymbol{S}}{d\mathbf{F}}\).

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

Compute the consistent tangent \(\frac{d\phi}{d\mathbf{F}}\) for a scalar quantity.

Given the inverse deformation gradient and the derivative of a scalar (e.g., an internal variable) with respect to the strain in Voigt notation, returns the corresponding 3×3 matrix mapping variations of the deformation gradient.

Parameters:

• FInv – Inverse of the deformation gradient at the end of the increment.

• dScalarDEps – Derivative of the scalar quantity with respect to the strain, in Voigt notation.

Returns:

The 3×3 sensitivity matrix \(\frac{d\phi}{d\mathbf{F}}\).

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 MarmotMaterialHypoElastic

#include <MarmotMaterialHypoElastic.h>

Inheritance 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

inline MarmotMaterialHypoElastic(const double *matProperties_, int nMaterialProperties_, int materialNumber_)

Constructs the material with a given set of material properties and an identifier.

Parameters:

• matProperties_ – [in] Pointer to the array of material properties.

• nMaterialProperties_ – [in] Number of entries in matProperties_.

• materialNumber_ – [in] Integer identifying this material instance.

virtual ~MarmotMaterialHypoElastic() = default

Default destructor.

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(state3D &state, Marmot::Matrix6d &dStress_dStrain, const Marmot::Vector6d &dStrain, const timeInfo &timeInfo) const = 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:

• state – [inout] A state3D instance carrying stress, strain energy, and state variables

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

• dStrain – [in] linearized strain increment

• timeInfo – Structure carrying the current (pseudo-)time and the (pseudo-)time increment

inline virtual void computeStressExplicit(state3D &state, const Marmot::Vector6d &dStrain, const timeInfo &timeInfo) const

Explicit version of computeStress for use in explicit time integration schemes. The algorithmic tangent is not needed in explicit schemes and will therefore not be computed.

Note

The default implementation calls computeStress and ignores the algorithmic tangent.

Note

Derived classes may override this method for efficiency reasons.

Parameters:

• state – [inout] A state3D instance carrying stress, strain energy, and state variables

• dStrain – [in] linearized strain increment

• timeInfo – Structure carrying time information

virtual void computePlaneStress(state2D &stress2D, Marmot::Matrix3d &dStress_dStrain2D, const Marmot::Vector3d &dStrain2D, const timeInfo &timeInfo) const

Plane stress implementation of computeStress.

virtual void computeUniaxialStress(state1D &stress1D, double &dStress_dStrain1D, const double dStrain, const timeInfo &timeInfo) const

Uniaxial stress implementation of computeStress.

inline StateView getStateView(const std::string &stateName, double *stateVars) const

Get a view to the state variables.

Parameters:

• stateName – Name of the state variable

• stateVars – Pointer to the state variable array

Returns:

StatView to access the state variable

inline int getNumberOfRequiredStateVars() const

Get the total number of required state variables.

Returns:

Total number of required state variables

inline virtual void initializeYourself(double *stateVars, int nStateVars)

Initialize the state variables at a material point.

Note

The default implementation initializes all state variables to zero.

Parameters:

• stateVars – Pointer to the state variable array

• nStateVars – Number of state variables

virtual double getMaximumWaveSpeed(const state3D &state) const

Get the maximum wave speed of the material.

The default implementation computes the 3D algorithmic tangent and returns sqrt(max(C_ii) / rho) with C_ii from the Voigt tangent diagonal entries.

Parameters:

state – [in] Current material state

Returns:

Maximum wave speed

virtual double getDensity(const double *stateVars) const = 0

Get the mass density of the material.

Parameters:

stateVars – Pointer to the state variable array

Returns:

Mass density of the material

Public Members

const int materialNumber

Integer identifier for this material instance.

MarmotStateLayoutDynamic stateLayout

Layout of the state variables.

double characteristicElementLength

Characteristic element length.

Protected Attributes

const double *materialProperties

Pointer to the array of material properties.

const int nMaterialProperties

Number of material properties.

class MarmotMaterialHypoElasticAD : public MarmotMaterialHypoElastic

#include <MarmotMaterialHypoElasticAD.h>

Inheritance diagram for MarmotMaterialHypoElasticAD:

Collaboration diagram for MarmotMaterialHypoElasticAD:

Derived class of MarmotMaterialHypoElastic providing automatic differentiation support for computing the algorithmic tangent via dual numbers.

Public Functions

virtual void computeStressAD(state3DAD &state, const Marmot::Vector6dual &dStrain, const timeInfo &timeInfo) const = 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:

• state – [inout] State carrying the dual Cauchy stress, strain energy, and state variables

• dStrain – [in] linearized strain increment

• timeInfo – Old (pseudo-)time

virtual void computeStress(state3D &state, Marmot::Matrix6d &dStressDDStrain, const Marmot::Vector6d &dStrain, const timeInfo &timeInfo) const 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:

• state – [inout] A state3D instance carrying stress, strain energy, and state variables

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

• dStrain – [in] linearized strain increment

• timeInfo – Structure carrying the current (pseudo-)time and the (pseudo-)time increment

inline MarmotMaterialHypoElastic(const double *matProperties_, int nMaterialProperties_, int materialNumber_)

Inherits the base-class constructor; see MarmotMaterialHypoElastic::MarmotMaterialHypoElastic().

class MarmotMaterialHypoElasticFactory

#include <MarmotMaterialHypoElasticFactory.h>

Factory class for creating material instances.

This class provides a mechanism to register materials by their code and name, and to create material instances based on their properties. It allows for dynamic material creation without hardcoding specific material types.

Public Types

using materialFactoryFunction = std::function<MarmotMaterialHypoElastic*(const double *materialProperties, int nMaterialProperties, int materialNumber)>

Function signature for material factory functions registered in the map.

Public Functions

MarmotMaterialHypoElasticFactory() = delete

Public Static Functions

static MarmotMaterialHypoElastic *createMaterial(const std::string &materialName, const double *materialProperties, int nMaterialProperties, int materialNumber)

Create a material instance based on its name and properties.

Parameters:

• materialName – [in] Registered name of the material.

• materialProperties – [in] Array of material properties.

• nMaterialProperties – [in] Number of properties in the array.

• materialNumber – [in] Unique identifier for the material instance.

Returns:

Pointer to the created MarmotMaterial instance, or nullptr if creation failed.

template<class T>
static inline bool registerMaterial(const std::string &materialName)

Register a material with its name and an auto-generated factory function.

Parameters:

materialName – [in] Unique name under which the material is registered.

Returns:

True if registration was successful, false if the name already exists.

Private Types

using MaterialFactoryMap = std::unordered_map<std::string, materialFactoryFunction>

Private Static Functions

static MaterialFactoryMap &materialFactoryFunctionByName()

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

• options – Solver options controlling iteration limits and tolerances

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() 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

std::unique_ptr<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_)

Constructs the manager from an existing statevar array and a pre-built layout.

Parameters:

• theStateVars – Pointer to the first element of the statevar array.

• theLayout_ – Layout describing the offsets and lengths of each named entry.

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

Reduction of the generalized failure criterion to a specific type.

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

Menetrey Willam parameters.

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)

Compute the polar radius and its first and second derivative with respect to the Lode angle.

Template Parameters:

T – scalar type

Parameters:

• theta – Lode angle

• e – eccentricity parameter

Returns:

tuple of polar radius, its first and second derivative with respect to the Lode angle

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 <PerezFougetSubstepperMarkII.h>

Sub-stepper for linear-elastic/elastoplastic materials using an implicit return-mapping.

Implements the Pérez–Fouget (CPPM) sub-stepping strategy. The consistent algorithmic tangent is assembled incrementally via sub-step contributions.

Template Parameters:

nSizeMatTangent – Size of the square material tangent matrix.

Public Types

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

Square matrix of size nSizeMatTangent used for the incremental tangent assembly.

Public Functions

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

Construct a PerezFougetSubstepper.

Parameters:

• initialStepSize – Initial sub-step size as a fraction of the total increment (0, 1].

• minimumStepSize – Minimum allowable sub-step size; triggers a warning if reached.

• scaleUpFactor – Factor by which the sub-step size is increased after nPassesToIncrease successes.

• scaleDownFactor – Factor applied to the sub-step size when convergence fails.

• nPassesToIncrease – Number of successful sub-steps before the step size may be increased.

• Cel – Linear-elastic stiffness matrix (6×6).

bool isFinished()

check if the subincrement process has finished

double getNextSubstep()

get the next subincrement size

bool decreaseSubstepSize()

decrease 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

TangentSizedMatrix consistentTangent

template<int sizeMaterialState>
class PerezFougetSubstepperTime

#include <PerezFougetSubstepperTime.h>

Pérez–Fouget sub-stepper with time-dependent elastic stiffness.

Compared to PerezFougetSubstepper, this variant does not require the elastic stiffness Cel at construction time; instead it is passed explicitly to extendConsistentTangent at each sub-step.

Template Parameters:

sizeMaterialState – Total size of the material state vector.

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)

Construct a PerezFougetSubstepperTime.

Parameters:

• initialStepSize – Initial sub-step size as a fraction of the total increment (0, 1].

• minimumStepSize – Minimum allowable sub-step size; triggers a warning if reached.

• scaleUpFactor – Factor by which the sub-step size is increased after nPassesToIncrease successes.

• scaleDownFactor – Factor applied to the sub-step size when convergence fails.

• nPassesToIncrease – Number of successful sub-steps before the step size may be increased.

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)

Accumulate the consistent tangent for an elastic sub-step.

Parameters:

CelT – Elastic stiffness matrix (6×6) at the current (end-of-sub-step) time.

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

Accumulate the consistent tangent for an elastoplastic sub-step.

Parameters:

• CelT – Elastic stiffness matrix (6×6) at the current time.

• matTangent – Inverse algorithmic tangent (TangentSizedMatrix); left-applied to the current accumulation.

Matrix6d consistentStiffness()

Return the assembled consistent algorithmic stiffness after all sub-steps are finished.

Returns:

The 6×6 consistent tangent operator.

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 state1D

#include <MarmotMaterialHypoElastic.h>

Structure holding the material state at a material point for 1D uniaxial stress.

Public Functions

inline state1D()

Public Members

double stress

1D Cauchy stress

double elasticEnergyDensity

Elastic strain energy density.

double dissipation

Dissipation.

double *stateVars

Pointer to array of state variables.

struct state2D

#include <MarmotMaterialHypoElastic.h>

Structure holding the material state at a material point for 2D plane stress.

Public Functions

inline state2D()

Public Members

Marmot::Vector3d stress

2D Cauchy stress tensor in Voigt notation

double elasticEnergyDensity

Elastic strain energy density.

double dissipation

Dissipation.

double *stateVars

Pointer to array of state variables.

struct state3D

#include <MarmotMaterialHypoElastic.h>

Structure to hold the material state at a material point in 3D.

Public Functions

inline state3D()

Default constructor for state3D Initializes stress to zero, energy and dissipation to zero, and stateVars to nullptr.

inline state3D(Marmot::Vector6d stress_, double elasticEnergyDensity_, double dissipation_, double *stateVars_)

Constructor for initializing state3D.

Parameters:

• stress_ – Cauchy stress tensor in Voigt notation

• elasticEnergyDensity_ – Elastic strain energy density

• dissipation_ – Dissipation

• stateVars_ – Pointer to array of state variables

Public Members

Marmot::Vector6d stress

Cauchy stress tensor in Voigt notation.

double elasticEnergyDensity

Elastic strain energy density.

double dissipation

Dissipation.

double *stateVars

Pointer to array of state variables.

struct state3DAD

#include <MarmotMaterialHypoElasticAD.h>

Structure holding the material state for 3D using dual numbers (for automatic differentiation).

Public Functions

inline state3DAD()

Default constructor for initializing state3DAD Initializes stress to zero, energy and dissipation to zero, and stateVars to nullptr.

inline state3DAD(const state3D &state)

Constructor for initializing state3DAD.

Parameters:

state – A state3D object containing the initial values for stress, energy density, dissipation, and state variables.

Public Members

Marmot::Vector6dual stress

Cauchy stress tensor in Voigt notation.

double elasticEnergyDensity

Elastic strain energy density.

double dissipation

Dissipation.

double *stateVars

Pointer to array of state variables.

struct StateVarEntryDefinition

#include <MarmotStateVarVectorManager.h>

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

Public Members

std::string name

Unique name identifying the state variable entry.

int length

Number of consecutive doubles occupied by this entry.

struct StateVarEntryLocation

#include <MarmotStateVarVectorManager.h>

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

Public Members

int index

Zero-based offset into the statevar array where the entry begins.

int length

Number of consecutive doubles occupied by this entry.

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

Map from entry name to its location.

int nRequiredStateVars

Total number of doubles required.

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.

struct timeInfo

#include <MarmotMaterialHypoElastic.h>

Structure carrying (pseudo-)time information passed to the material routines.

Public Members

double time

Current (pseudo-)time.

double dT

(Pseudo-)time increment from the old (pseudo-)time to the current (pseudo-)time

template<int nYieldSurfaces>
class YieldSurfaceCombinationManager

#include <YieldSurfaceCombinationManager.h>

Tracks and iterates over yield-surface activity combinations for multi-surface plasticity.

Maintains a table of all \(2^n - 1\) non-empty subsets of nYieldSurfaces yield surfaces and provides an interface to iterate over them during the return-mapping algorithm, marking already-tried combinations as used.

Template Parameters:

nYieldSurfaces – Number of yield surfaces.

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)

Signals that the current increment should be discarded.

Sets pNewDT to the requested value and prints a warning message. Intended to be called inside material routines when convergence cannot be achieved, instructing the host code to restart with a smaller step.

Parameters:

• pNewDT – Output parameter: suggested new time-step size.

• value – Suggested value for the new time-step size.

• message – Human-readable reason for discarding the increment.

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.

Marmot::EigenTensors::Tensor633d dGreenLagrangedDeformationGradient(const Eigen::Matrix3d &F)

Computes the derivative of the Green-Lagrange strain tensor with respect to the deformation gradient.

Parameters:

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

Returns:

A \(6 \times 3 \times 3\) tensor representing the derivative of the Green-Lagrange strain tensor.

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)

Computes the acoustic tensor for a given material tangent and normal vector.

Parameters:

• materialTangent – Fourth-order material tangent in Voigt notation.

• normalVector – Candidate band normal \(\mathbf{n}\).

Returns:

Acoustic tensor \(\mathbf{Q} = \mathbf{n} \cdot \mathbb{C} \cdot \mathbf{n}\).

bool localizationChecker(const Marmot::Matrix3d &acousticTensor)

Checks whether the given acoustic tensor indicates localization.

Parameters:

acousticTensor – Acoustic tensor to test.

Returns:

true if \(\det(\mathbf{Q}) \le 0\) (localization detected).

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

Constructs a unit normal vector from spherical angles.

Parameters:

• alpha – Polar angle (rad).

• beta – Azimuthal angle (rad).

Returns:

Unit normal vector \(\mathbf{n}\).

double minimumDeterminantAcousticTensor(const Marmot::Matrix6d &materialTangent)

Computes the minimum determinant of the acoustic tensor over all band normals.

Parameters:

materialTangent – Fourth-order material tangent in Voigt notation.

Returns:

Minimum determinant of the acoustic tensor.

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 EigenTensors::Tensor3333d &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.

EigenTensors::Tensor3333d 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 \).

Marmot::FastorStandardTensors::Tensor3333d voigtToStiffnessFastor(const Marmot::Matrix6d &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:

a Fastor::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)

Computes the three principal geostatic stress components at a given depth.

The stress distribution is assumed to vary linearly with the vertical coordinate coordinate_y as specified by geostaticStressDefintion.

Parameters:

• geostaticStressDefintion – Pointer to the array defining the linear distribution.

• coordinate_y – Vertical coordinate of the evaluation point.

Returns:

Tuple of three stress components (vertical, horizontal-1, horizontal-2).

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 coefficient.

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 Solvers

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

Template Parameters:

nDim – Number of dimensions (2D or 3D)

Parameters:

• 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

namespace MarmotLibrary

file AdaptiveSubstepper.h

#include “Marmot/MarmotJournal.h”

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for AdaptiveSubstepper.h:

Adaptive sub-stepping with Richardson extrapolation for implicit return-mapping.

file AdaptiveSubstepperExplicit.h

#include “Marmot/MarmotJournal.h”

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for AdaptiveSubstepperExplicit.h:

Adaptive sub-stepping with Richardson extrapolation for semi-explicit return-mapping.

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:

Haigh–Westergaard stress and strain invariant coordinates.

Provides the HaighWestergaardCoordinates aggregate and functions for computing the hydrostatic component \(\xi\), deviatoric radius \(\rho\), and Lode angle \(\theta\) from stress or strain tensors in Voigt notation.

file HughesWinget.h

#include “Marmot/MarmotVoigt.h”

Include dependency graph for HughesWinget.h:

Hughes–Winget objective stress-rate integrator.

Provides the HughesWinget class that, given old and new deformation gradients, computes the objective strain increment, the incremental rotation tensor, and the tangent operator contributions required by hypoelastic constitutive updates.

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:

Utilities for evaluating geostatic (gravity-induced) initial stress states.

Provides functions for computing the in-situ stress components from a linearly-varying stress distribution defined along a vertical coordinate axis.

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:

Localization (strain-softening band) analysis via the acoustic tensor.

Functions for checking material instability by evaluating whether the acoustic tensor \(\mathbf{Q}(\mathbf{n}) = \mathbf{n} \cdot \mathbb{C} \cdot \mathbf{n}\) loses positive definiteness for some normal vector \(\mathbf{n}\).

file MarmotLowerDimensionalStress.h

#include “Marmot/MarmotVoigt.h”

Include dependency graph for MarmotLowerDimensionalStress.h:

file MarmotMaterialHypoElastic.h

#include “Marmot/MarmotStateHelpers.h”

#include “Marmot/MarmotTypedefs.h”

#include <algorithm>

#include <cmath>

#include <vector>

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 dependency graph for MarmotMaterialHypoElasticAD.h:

file MarmotMaterialHypoElasticFactory.h

#include “Marmot/MarmotMaterialHypoElastic.h”

#include <cassert>

#include <functional>

#include <string>

#include <unordered_map>

Include dependency graph for MarmotMaterialHypoElasticFactory.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:

Miscellaneous utility helpers for controlling increment sizes.

file MarmotViscoelasticity.h

#include “Marmot/MarmotVoigt.h”

Include dependency graph for MarmotViscoelasticity.h:

file MarmotVoigt.h

#include “Marmot/MarmotFastorTensorBasics.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:

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:

Sub-stepper for semi-explicit return-mapping algorithms (Pérez–Fouget algorithm, Mark II).

file PerezFougetSubstepperMarkII.h

#include “Marmot/MarmotJournal.h”

#include “Marmot/MarmotMath.h”

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for PerezFougetSubstepperMarkII.h:

Sub-stepper for implicit return-mapping algorithms (Pérez–Fouget algorithm, Mark II).

file PerezFougetSubstepperTime.h

#include “Marmot/MarmotTypedefs.h”

#include “Marmot/MarmotJournal.h”

Include dependency graph for PerezFougetSubstepperTime.h:

Time-variant variant of the Pérez–Fouget sub-stepper.

A modified version of the Pérez–Fouget sub-stepper that accounts for a time-variant elastic stiffness tensor \( \mathbb{C}^\mathrm{el}(t_{n+1}) \). The algorithmic formulation is unchanged; the only difference is that the current elastic tangent must be supplied when extending the consistent tangent.

file YieldSurfaceCombinationManager.h

#include “Marmot/MarmotTypedefs.h”

Include dependency graph for YieldSurfaceCombinationManager.h:

Manager for multi-surface yield-surface combinations in multi-surface plasticity.

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

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