MarmotMathCore

struct ExponentialMapFailed : public std::exception

Inheritance diagram for Marmot::ContinuumMechanics::TensorUtility::TensorExponential::ExponentialMapFailed:

Collaboration diagram for Marmot::ContinuumMechanics::TensorUtility::TensorExponential::ExponentialMapFailed:

class NewtonConvergenceChecker

#include <NewtonConvergenceChecker.h>

A class to check convergence of Newton-Raphson iterations. This class provides methods to evaluate the convergence of Newton-Raphson iterations based on residual norms and relative increments.

Public Functions

NewtonConvergenceChecker(const Eigen::VectorXd &residualScaleVector, int nMaxNewtonCycles, int nMaxNewtonCyclesAlt, double newtonTol, double newtonRTol, double newtonTolAlt, double newtonRTolAlt)

Constructor for NewtonConvergenceChecker.

Parameters:

• residualScaleVector – A vector to scale the residual for convergence checks.

• nMaxNewtonCycles – Maximum number of Newton cycles allowed.

• nMaxNewtonCyclesAlt – Maximum number of alternative Newton cycles allowed.

• newtonTol – Tolerance for the Newton convergence based on absolute residual norm.

• newtonRTol – Tolerance for the Newton convergence based on relative increment norm.

• newtonTolAlt – Alternative tolerance for the Newton convergence based on absolute residual norm.

• newtonRTolAlt – Alternative tolerance for the Newton convergence based on relative increment norm.

double relativeNorm(const Eigen::VectorXd &increment, const Eigen::VectorXd &reference)

Compute the relative norm of the increment with respect to a reference vector.

The relative norm is computed as

\[\begin{split}\displaystyle r_{\rm rel} = \begin{cases} || \Delta \boldsymbol{X} ||_2 & \text{if } || \Delta \boldsymbol{X} ||_2 < 10^{-14} \\ 0.0 & \text{if } || \boldsymbol{X}_{\rm ref} ||_2 < 10^{-12} \\ \frac{ || \Delta \boldsymbol{X} ||_2 }{ || \boldsymbol{X}_{\rm ref} ||_2 } & \text{else} \end{cases} \end{split}\]

where \( || \square ||_2 \) denotes the Euclidean norm.

Parameters:

• increment – The increment vector.

• reference – The reference vector.

Returns:

The relative norm of the increment.

double residualNorm(const Eigen::VectorXd &Residual)

Compute the norm of the residual vector.

The residual is first scaled by the residualScaleVector before computing the norm.

This function computes

\[r = \sqrt{ R^{\rm s}_i R^{\rm s}_i } = || \boldsymbol{R}^{\rm s} ||_2 \]

where \( R^{\rm s}_i = s_i R_i \) (no summation over i) are the scaled residual components and \( s_i \) are the components of the residualScaleVector.

Parameters:

Residual – The residual vector.

Returns:

The norm of the scaled residual.

bool iterationFinished(const Eigen::VectorXd &residual, const Eigen::VectorXd &X, const Eigen::VectorXd &dX, int numberOfIterations)

Check if the iteration has finished based on residuals and increments.

The iteration is considered finished if either the solution has converged or the maximum number of alternative Newton cycles has been excee:ed:

if ( isConverged( residual, X, dX, numberOfIterations ) || numberOfIterations > nMaxNewtonCyclesAlt )
  return true;
else
 return false;

Parameters:

• residual – The current residual vector.

• X – The current solution vector.

• dX – The current increment vector.

• numberOfIterations – The current iteration count.

Returns:

True if the iteration has finished, false otherwise.

bool isConverged(const Eigen::VectorXd &residual, const Eigen::VectorXd &X, const Eigen::VectorXd &dX, int numberOfIterations)

Check if the solution has converged based on residuals and increments.

The solution is considered converged if both the absolute residual norm and the relative increment norm is below the specified tolerances. Alternative tolerances are used if the number of iterations exceeds nMaxNewtonCycles.

if ( numberOfIterations <= nMaxNewtonCycles ) {
  if ( relativeNorm( dX, X ) <= newtonRTol && residualNorm( residual ) <= newtonTol )
    return true;
} else if ( numberOfIterations <= nMaxNewtonCyclesAlt + 1 ){
  if ( relativeNorm( dX, X ) <= newtonRTolAlt && residualNorm( residual ) <= newtonTolAlt )
    return true;
}
else
  return false;

Parameters:

• residual – The current residual vector.

• X – The current solution vector.

• dX – The current increment vector.

• numberOfIterations – The current iteration count.

Returns:

True if the solution has converged, false otherwise.

Private Members

const Eigen::VectorXd residualScaleVector

A vector to scale the residual for convergence checks.

const int nMaxNewtonCycles

Maximum number of Newton cycles allowed.

const int nMaxNewtonCyclesAlt

Maximum number of alternative Newton cycles allowed.

const double newtonTol

Tolerance for the Newton convergence based on absolute residual norm.

const double newtonRTol

Tolerance for the Newton convergence based on relative increment norm.

const double newtonTolAlt

Alternative tolerance for the Newton convergence based on absolute residual norm.

const double newtonRTolAlt

Alternative tolerance for the Newton convergence based on relative increment norm.

template<typename T, std::size_t tensorSize>
struct TensorExponentialResult

Public Members

Fastor::Tensor<T, tensorSize, tensorSize> theExponential

Fastor::Tensor<T, tensorSize, tensorSize, tensorSize, tensorSize> theExponentialDerivative

namespace autodiff

namespace Eigen

namespace Marmot

Typedefs

typedef Eigen::Matrix<double, 6, 6> Matrix6d

typedef Eigen::Matrix<double, 6, 9> Matrix69d

typedef Eigen::Matrix<double, 9, 9> Matrix99d

typedef Eigen::Matrix<double, 3, 4> Matrix34d

typedef Eigen::Map<Matrix6d> mMatrix6d

typedef Eigen::Matrix<double, 3, 3> Matrix3d

typedef Eigen::Matrix<double, 3, 1> Vector3d

typedef Eigen::Matrix<double, 4, 1> Vector4d

typedef Eigen::Matrix<double, 6, 1> Vector6d

typedef Eigen::Matrix<double, 7, 1> Vector7d

typedef Eigen::Matrix<double, 8, 1> Vector8d

typedef Eigen::Matrix<double, 9, 1> Vector9d

typedef Eigen::Matrix<int, 8, 1> Vector8i

typedef Eigen::Matrix<double, 1, 6> RowVector6d

typedef Eigen::Map<Vector6d> mVector6d

typedef Eigen::Map<Eigen::VectorXd> mVectorXd

typedef Eigen::Map<const Marmot::Vector6d> mConstVector6d

typedef Eigen::Matrix<double, 3, 6> Matrix36d

typedef Eigen::Matrix<double, 3, 6> Matrix36

typedef Eigen::Matrix<double, 6, 3> Matrix63d

typedef Eigen::Matrix<double, 9, 9> Matrix9d

typedef std::complex<double> complexDouble

typedef Eigen::Matrix<complexDouble, 6, 1> Vector6cd

template<typename T>
using Vector6t = Eigen::Matrix<T, 6, 1>

template<typename T>
using VectorXt = Eigen::Matrix<T, -1, 1>

Enums

enum DimensionType

Dimension types for reduceTo2D and expandTo3D.

Values:

enumerator U

displacement dimension

enumerator W

rotation dimension

Functions

template<typename T, size_t nRows, size_t nCols, typename = std::enable_if<!std::is_const<std::remove_reference<T>>::value>>
inline auto mapEigenToFastor(const Fastor::Tensor<T, nRows, nCols> &fastor)

Map a Fastor Tensor (const) to an Eigen Map (row-major, const)

Note

This function works for second rank tensors (matrices) of size nRows x nCols.

Template Parameters:

• T – scalar type

• nRows – number of rows

• nCols – number of columns

Parameters:

fastor – a Fastor Tensor

Returns:

an Eigen Map

template<typename T, size_t nRows, size_t nCols, typename = void>
inline auto mapEigenToFastor(Fastor::Tensor<T, nRows, nCols> &fastor)

Map a Fastor Tensor to an Eigen Map (row-major)

Note

This function works for second rank tensors (matrices) of size nRows x nCols.

Template Parameters:

• T – scalar type

• nRows – number of rows

• nCols – number of columns

Parameters:

fastor – a Fastor Tensor

Returns:

an Eigen Map

template<typename T, size_t nRows, typename = void>
inline auto mapEigenToFastor(const Fastor::Tensor<T, nRows> &fastor)

Map a Fastor Tensor to an Eigen Map (column vector)

Note

This function works for first rank tensors (vectors) of size nRows.

Template Parameters:

• T – scalar type

• nRows – number of rows

Parameters:

fastor – a Fastor Tensor

Returns:

an Eigen Map

template<typename T, size_t nRows, typename = void>
inline auto mapEigenToFastor(const Fastor::TensorMap<T, nRows> &fastor)

Map a Fastor TensorMap to an Eigen Map (column vector)

Note

This function works only for first rank tensors (vectors) of size nRows.

Template Parameters:

• T – scalar type

• nRows – number of rows

Parameters:

fastor – a Fastor TensorMap

Returns:

an Eigen Map

template<typename T, size_t nRows, size_t nCols, typename = void>
inline auto mapEigenToFastor(const Fastor::TensorMap<T, nRows, nCols> &fastor)

Map a Fastor TensorMap to an Eigen Map (row-major)

Note

This function works only for second rank tensors (matrices) of size nRows x nCols.

Template Parameters:

• T – scalar type

• nRows – number of rows

• nCols – number of columns

Parameters:

fastor – a Fastor TensorMap

Returns:

an Eigen Map

template<template<typename, size_t...> class TensorType, typename T, size_t... Rest>
inline void copyFastorToColumnMajor(T *target, const TensorType<T, Rest...> &source)

Copy a Fastor tensor to a column-major array.

Template Parameters:

• TensorType – a Fastor tensor type

• T – scalar type

• Rest – a pack of sizes specifying the dimensions of the tensor

Parameters:

• target – pointer to the target array

• source – a Fastor tensor

template<DimensionType... dims, typename T, size_t... dims3D>
inline auto reduceTo2D(const Fastor::Tensor<T, dims3D...> &theTensor3D)

Reduce a 3D Fastor tensor to a 2D Fastor tensor.

This function is used to reduce a displacement and/or rotation tensor in 3D to a displacement and/or rotation tensor in 2D. The conversion is specified by the pack of DimensionType. For displacements (DimensionType U) the dimension is reduced from 3 to 2, for rotations (DimensionType W) the dimension is reduced from 3 to 1.

Template Parameters:

• dims – a pack of DimensionType specifying the conversion from 3D to 2D

• dims3D – a pack of sizes specifying the dimensions of the 3D tensor

Parameters:

• T – scalar type

• theTensor3D – a 3D input tensor

Returns:

the reduced 2D tensor

template<typename Derived, size_t order>
inline auto reduceTo2D(const Fastor::AbstractTensor<Derived, order> &theTensor3D)

Overload of reduceTo2d for AbstractTensor inputs.

Template Parameters:

• Derived – derived scalar type

• order – order of the tensor

Parameters:

theTensor3D – a 3D input tensor

Returns:

the reduced 2D tensor

int const3(size_t x)

Helper function to return 3 for any input size_t.

Parameters:

x – input size_t

Returns:

always returns 3

template<typename T, size_t... dims2D>
inline auto expandTo3D(const Fastor::Tensor<T, dims2D...> &theTensor2D)

Expand a 2D Fastor tensor to a 3D Fastor tensor.

This function is used to expand a displacement and/or rotation tensor in 2D to a displacement and/or rotation tensor in 3D. The conversion is specified by the pack of DimensionType. For displacements (DimensionType U) the dimension is expanded from 2 to 3, for rotations (DimensionType W) the dimension is expanded from 1 to 3.

Template Parameters:

• dims – a pack of DimensionType specifying the conversion from 2D to 3D

• dims2D – a pack of sizes specifying the dimensions of the 2D Tensor

Parameters:

• T – scalar type

• theTensor2D – a 2D input tensor

Returns:

the expanded 3D tensor

template<typename Derived, size_t order>
inline auto expandTo3D(const Fastor::AbstractTensor<Derived, order> &theTensor2D)

Overload of expandTo3D for AbstractTensor inputs.

Template Parameters:

• Derived – derived scalar type

• order – order of the tensor

Parameters:

theTensor2D – a 2D input tensor

Returns:

the expanded 3D tensor

template<typename T, size_t... Rest>
Fastor::Tensor<T, Rest...> multiplyFastorTensorWithScalar(Fastor::Tensor<T, Rest...> tensor, T scalar)

Multiply a Fastor tensor with a scalar.

Note

This function is a workaround for a bug in Fastor (issue #149).

Template Parameters:

• T – scalar type

• Rest – a pack of sizes specifying the dimensions of the tensor

Parameters:

• tensor – a Fastor tensor

• scalar – a scalar

Returns:

the resulting Fastor tensor

template<typename T>
T einsum_ij_ij_hardcoded(const FastorStandardTensors::Tensor33t<T> &A, const FastorStandardTensors::Tensor33t<T> &B)

Hardcoded implementation of the Einstein summation for two second order tensors.

This function computes the Einstein summation for two second order tensors A and B: \( C = A_{ij} B_{ij} \)

Note

This function is a workaround for cases in which Fastor::to_scalar does not work. This is the case for example when T is autodiff::dual.

Template Parameters:

T – scalar type

Parameters:

• A – first second order tensor

• B – second second order tensor

Returns:

the resulting scalar

template<typename T, size_t... Rest>
Fastor::Tensor<T, Rest...> fastorTensorFromDoubleTensor(const Fastor::Tensor<double, Rest...> &in)

Construct a arbitrary type Fastor tensor from a Fastor double tensor.

Note

This function is a workaround for lack of casting of pointer types e.g. no cast available for autodiff::dual* and double*

Template Parameters:

T – target scalar type

Parameters:

• Rest – a pack of sizes specifying the dimensions of the tensor

• in – a Fastor tensor of type double

Returns:

a Fastor tensor of type T

template<typename T, size_t... Rest>
Fastor::Tensor<T, Rest...> fastorTensorFromDoubleTensorMap(const Fastor::TensorMap<double, Rest...> &in)

Construct a arbitrary type Fastor tensor from a Fastor double tensor map.

Note

This function is a workaround for lack of casting of pointer types e.g. no cast available for autodiff::dual* and double*

Template Parameters:

T – target scalar type

Parameters:

• Rest – a pack of sizes specifying the dimensions of the tensor

• in – a Fastor tensor map of type double

Returns:

a Fastor tensor of type T

template<typename T, size_t... Rest>
Fastor::Tensor<double, Rest...> makeReal(const Fastor::Tensor<T, Rest...> &in)

Extract the real part of a Fastor tensor of arbitrary type.

Template Parameters:

T – source scalar type

Parameters:

• Rest – a pack of sizes specifying the dimensions of the tensor

• in – a Fastor tensor of type T

Returns:

a Fastor tensor of type double

template<typename T, size_t... Rest>
Fastor::Tensor<autodiff::dual, Rest...> makeDual(const Fastor::Tensor<T, Rest...> &in)

Construct a Fastor tensor of autodiff::dual from a Fastor tensor of arbitrary type.

Template Parameters:

T – source scalar type

Parameters:

• Rest – a pack of sizes specifying the dimensions of the tensor

• in – a Fastor tensor of type T

Returns:

a Fastor tensor of type autodiff::dual

template<size_t order, size_t... Rest>
Fastor::Tensor<autodiff::HigherOrderDual<order, double>, Rest...> makeHigherOrderDual(const Fastor::Tensor<double, Rest...> &in)

Construct a Fastor tensor of autodiff::HigherOrderDual from a Fastor double tensor.

Template Parameters:

order – order of the HigherOrderDual

Parameters:

• Rest – a pack of sizes specifying the dimensions of the tensor

• in – a Fastor tensor of type double

Returns:

a Fastor tensor of type autodiff::HigherOrderDual

template<typename T>
FastorStandardTensors::Tensor33t<T> symmetric(const FastorStandardTensors::Tensor33t<T> &t)

Compute the symmetric part of a second order 3D Fastor tensor.

It actually computes \( \text{sym}\, t_ij = 0.5 ( t_{ij} + t_{ji} ) \)

Template Parameters:

T – scalar type

Parameters:

t – a second order Fastor tensor

Returns:

the symmetric part of the tensor

template<typename T>
FastorStandardTensors::Tensor33t<T> deviatoric(const FastorStandardTensors::Tensor33t<T> &t)

Compute the deviatoric part of a second order 3D Fastor tensor.

It actually computes \( \text{dev} \, t_{ij} = t_{ij} - \frac{1}{3}t_{kk}\delta_{ij}\)

Template Parameters:

T – scalar type

Parameters:

t – a second order Fastor tensor

Returns:

the deviatoric part of the tensor

namespace AutomaticDifferentiation

Typedefs

template<size_t... Rest>
using tensor_to_scalar_function_type = std::function<dual(const Fastor::Tensor<dual, Rest...> &T)>

Alias for a function mapping a tensor to a scalar with dual numbers.

Template Parameters:

Rest – dimensions of the tensor

template<size_t order, size_t... Rest>
using tensor_to_scalar_function_type_arbitrary_dual_order = std::function<HigherOrderDual<order, double>(const Fastor::Tensor<HigherOrderDual<order, double>, Rest...> &T)>

A type alias for a function mapping a tensor to a scalar with dual numbers of arbitrary order.

Template Parameters:

• order – current order of the dual numbers in the tensor

• Rest – dimensions of the tensor

using scalar_to_scalar_function_type = std::function<dual(const dual&)>

A type alias for a scalar-to-scalar function that takes and returns dual numbers.

using scalar_to_scalar_function_type_2nd = std::function<dual2nd(const dual2nd&)>

A type alias for a scalar-to-scalar function that takes and returns second order dual numbers.

using vector_to_vector_function_type_dual = std::function<VectorXdual(const VectorXdual &X)>

A type alias for a vector-to-vector function that takes and returns dual-valued vectors.

using vector_to_vector_function_type_dual2nd = std::function<VectorXdual2nd(const VectorXdual2nd &X)>

A type alias for a vector-to-vector function that takes and returns second order dual-valued vectors.

Functions

template<size_t order, size_t... Rest>
Fastor::Tensor<HigherOrderDual<order + 1, double>, Rest...> increaseDualOrderWithShift(const Fastor::Tensor<HigherOrderDual<order, double>, Rest...> &in)

Increase the order of an arbitray order dual-valued Fastor tensor by one and shift the derivatives.

Template Parameters:

• order – current order of the dual numers in the tensor

• Rest – dimensions of the tensor

Parameters:

in – input tensor

Returns:

output tensor with increased order dual numbers

template<size_t... Rest>
Fastor::Tensor<double, Rest...> df_dT(const tensor_to_scalar_function_type<Rest...> &f, const Fastor::Tensor<double, Rest...> &T)

Compute the gradient of a tensor-to-scalar function.

Template Parameters:

Rest – dimensions of the tensor

Parameters:

• f – function mapping a tensor to a scalar

• T – input tensor at which the gradient is evaluated

Returns:

gradient of f with respect to T, same shape as T

template<size_t order, size_t... Rest>
std::pair<HigherOrderDual<order, double>, Fastor::Tensor<HigherOrderDual<order, double>, Rest...>> df_dT(const tensor_to_scalar_function_type_arbitrary_dual_order<order + 1, Rest...> &f, const Fastor::Tensor<HigherOrderDual<order, double>, Rest...> &T)

Compute the gradient of a tensor-to-scalar function where the with dual numbers of arbitrary order.

Template Parameters:

• order – current order of the dual numers in the tensor

• Rest – dimensions of the tensor

Parameters:

• f – function mapping a tensor to a scalar

• T – input tensor at which the gradient is evaluated

Returns:

pair of function value and gradient of f with respect to T, same shape as T

template<size_t... RestF, size_t... RestT>
std::pair<Fastor::Tensor<double, RestF...>, Fastor::Tensor<double, RestF..., RestT...>> dF_dT(std::function<Fastor::Tensor<dual, RestF...>(const Fastor::Tensor<dual, RestT...>&)> &F, const Fastor::Tensor<double, RestT...> &T)

Compute the gradient of a tensor-to-tensor function.

Template Parameters:

• RestF – dimensions of the output tensor

• RestT – dimensions of the input tensor

Parameters:

• F – function mapping a tensor to a tensor

• T – input tensor at which the gradient is evaluated

Returns:

pair of function value and gradient of F with respect to T, gradient has shape (RestF…, RestT…)

dual2nd shiftTo2ndOrderDual(const dual &x)

Creates a 2nd order hyper-dual number from a dual number.

This function takes a dual number as input and returns a 2nd order hyper-dual number. This is done by keeping the real part and shifting the first derivative part to the second derivative part.

Parameters:

x – Input dual number

Returns:

2nd order hyper-dual number

VectorXdual2nd shiftTo2ndOrderDual(const VectorXdual &X)

Creates a vector of 2nd order hyper-dual numbers from a vector of dual numbers.

This function takes a vector of dual numbers as input and returns a vector of 2nd order hyper-dual numbers. This is done by keeping the real part and shifting the first derivative part to the second derivative part for each element in the vector.

Parameters:

X – Input vector of dual numbers

Returns:

Vector of 2nd order hyper-dual numbers

template<size_t order, typename T, typename G>
const auto &valnode(const Dual<T, G> &dual)

Access the value node of a dual number at a specific order.

This function recursively accesses the value node of a dual number at a specific order. It uses template metaprogramming to ensure that the order is valid for the given dual number type. If the order is 0, it returns the value of the dual number. If the order is greater than 0, it recursively calls itself on the value of the dual number until it reaches the desired order.

Template Parameters:

• order – The order of the dual number to access

• T – The type of the value stored in the dual number

• G – The type of the gradient stored in the dual number

Parameters:

dual – The dual number to access

Returns:

A reference to the value node at the specified order

template<size_t order, typename T, typename G>
auto &valnode(Dual<T, G> &dual)

template<size_t order>
autodiff::HigherOrderDual<order + 1, double> increaseDualOrderWithShift(const autodiff::HigherOrderDual<order, double> &in)

Increases the order of a hyper-dual number by 1.

This function takes a dual number of a given order and creates a new dual number with the order increased by 1. The derivative parts are elevated (shifted) one order.

Template Parameters:

order – The current order of the dual number

Parameters:

in – The input dual number

Returns:

A new dual number with order increased by 1

template<size_t order>
autodiff::HigherOrderDual<order - 1, double> decreaseDualOrder(autodiff::HigherOrderDual<order, double> &in)

Decreases the order of a hyper-dual number by 1.

This function takes a dual number of a given order and creates a new dual number with the order decreased by 1. The derivative parts are copied to the lower order dual number. The highest derivative part of the input dual number is discarded.

Template Parameters:

order – The current order of the dual number

Parameters:

in – The input dual number

Returns:

A new dual number with order decreased by 1

template<size_t order>
autodiff::HigherOrderDual<order - 1, double> decreaseDualOrderWithShift(autodiff::HigherOrderDual<order, double> &in)

Decreases the order of a hyper-dual number by 1 with a shift.

This function takes a dual number of a given order and creates a new dual number with the order decreased by 1. The derivative parts are shifted down by one order. The first derivative part of the input dual number is discarded.

Template Parameters:

order – The current order of the dual number

Parameters:

in – The input dual number

Returns:

A new dual number with order decreased by 1

template<size_t order>
Vector<HigherOrderDual<order + 1, double>, -1> increaseDualOrderWithShift(const Vector<HigherOrderDual<order, double>, -1> &in)

Increases the order of a vector of hyper-dual numbers by 1.

This function takes a vector of dual numbers of a given order and creates a new vector of dual numbers with the order increased by 1. The derivative parts are elevated (shifted) one order for each element in the vector.

Template Parameters:

order – The current order of the dual numbers in the vector

Parameters:

in – The input vector of dual numbers

Returns:

A new vector of dual numbers with order increased by 1

double df_dx(const scalar_to_scalar_function_type &f, const double &x)

Computes the derivative of a scalar-to-scalar function at a given point using automatic differentiation.

This function takes a scalar-to-scalar function and a point as input, and returns the derivative of the function at that point. It uses automatic differentiation to compute the derivative accurately.

Parameters:

• f – The scalar-to-scalar function to differentiate

• x – The point at which to evaluate the derivative

Returns:

The derivative of the function at the given point

dual df_dx(const scalar_to_scalar_function_type_2nd &f, const dual &x)

Computes the derivative of a dual-valued scalar-to-scalar function at a given point using automatic differentiation.

This function takes a 2nd order dual-valued scalar-to-scalar function and a point as input, and returns the derivative of the function at that point. Note that the function is 2nd order dual-valued, but the input is only 1st order dual-valued.

Parameters:

• f – The scalar-to-scalar function to differentiate

• x – The point at which to evaluate the derivative

Returns:

The derivative of the function at the given point

std::pair<VectorXd, MatrixXd> dF_dX(const vector_to_vector_function_type_dual &F, const VectorXd &X)

Computes the Jacobian of a vector-to-vector function at a given point using automatic differentiation.

This function takes a vector-to-vector function and a point as input, and returns a pair containing the function value and the Jacobian matrix at that point. It uses automatic differentiation to compute the Jacobian accurately.

Parameters:

• F – The vector-to-vector function to differentiate

• X – The point at which to evaluate the Jacobian

Returns:

A pair containing the function value and the Jacobian matrix at the given point

std::pair<VectorXdual, MatrixXdual> dF_dX_2nd(const vector_to_vector_function_type_dual2nd &F, const VectorXdual &X)

Computes the Jacobian of a dual-valued vector-to-vector function at a given point using automatic differentiation.

This function takes a 2nd order dual-valued vector-to-vector function and a point as input, and returns a dual-valued pair containing the function value and the Jacobian matrix at that point. Note that the function is 2nd order dual-valued, but the input is only 1st order dual-valued.

Parameters:

• F – The vector-to-vector function to differentiate

• X – The point at which to evaluate the Jacobian

Returns:

A pair containing the function value and the Jacobian matrix at the given point

namespace SecondOrder

Typedefs

template<size_t dim>
using tensor_to_scalar_function_type = std::function<dual2nd(const Fastor::Tensor<dual2nd, dim, dim> &T)>

Alias for a function mapping a tensor to a scalar with second order dual numbers.

Note

The implementation is currently limited to second rank (dim x dim) tensors

Template Parameters:

dim – dimension of the input tensor (dim x dim)

template<size_t dim>
using tensor_and_scalar_to_scalar_function_type = std::function<dual2nd(const Fastor::Tensor<dual2nd, dim, dim> &T, const dual2nd scalar)>

Alias for a function mapping a tensor and a scalar to a scalar with second order dual numbers.

Note

The implementation is currently limited to second rank (dim x dim) tensors

Template Parameters:

dim – dimension of the input tensor (dim x dim)

Functions

template<size_t dim>
std::tuple<double, Fastor::Tensor<double, dim, dim>, Fastor::Tensor<double, dim, dim, dim, dim>> d2f_dT2(const tensor_to_scalar_function_type<dim> &F, const Fastor::Tensor<double, dim, dim> &T)

Computes the second derivative of a function that maps a tensor to a scalar.

Note

The implementation is currently limited to second rank (dim x dim) tensors

Template Parameters:

dim – dimension of the input tensor (dim x dim)

Parameters:

• F – function mapping a (dim x dim) tensor to a scalar

• T – input tensor at which the derivative is evaluated

Returns:

tuple of function value, first derivative and second derivative of F with respect to T first derivative has shape (dim, dim) second derivative has shape (dim, dim, dim, dim)

template<size_t dim>
Fastor::Tensor<double, dim, dim> d2f_dTensor_dScalar(const tensor_and_scalar_to_scalar_function_type<dim> &F, const Fastor::Tensor<double, dim, dim> &T, const double scalar)

Computes the mixed second derivative of a function that maps a tensor and a scalar to a scalar.

Note

The implementation is currently limited to second rank (dim x dim) tensors

Template Parameters:

dim – dimension of the input tensor (dim x dim)

Parameters:

• F – function mapping a (dim x dim) tensor and a scalar to a scalar

• T – input tensor at which the derivative is evaluated

• scalar – input scalar at which the derivative is evaluated

Returns:

mixed second derivative of F with respect to T and the scalar, has shape (dim, dim)

namespace Constants

Functions

inline double cubicRootEps()

Inlined function to get the cubic root of machine epsilon for double precision.

This function computes and returns the cubic root of the machine epsilon for double precision floating-point numbers. Machine epsilon is the smallest positive number that, when added to 1.0, results in a value different from 1.0 due to rounding errors.

Returns:

The cubic root of machine epsilon for double precision.

inline double squareRootEps()

Inlined function to get the square root of machine epsilon for double precision.

This function computes and returns the square root of the machine epsilon for double precision floating-point numbers. Machine epsilon is the smallest positive number that, when added to 1.0, results in a value different from 1.0 due to rounding errors.

Returns:

The square root of machine epsilon for double precision.

Variables

double Pi = 3.141592653589793238463

Constant expression for $\pi$.

double numZeroPos = 1e-16

Constant expression for 1e-16.

double GoldenRatio = 1.618033988749895

The golden ratio, approximately 1.618033988749895.

double GoldenAngle = 2.399963229728653

The golden angle in radians, approximately 2.399963229728653.

const double SquareRootEps = squareRootEps()

const double CubicRootEps = cubicRootEps()

double sqrt3_8 = 0.61237243569579452454932101867647

Constant expression for \(\sqrt{3/8}\).

double sqrt2_3 = 0.8164965809277260327324280249019

Constant expression for \(\sqrt{2/3}\).

double sqrt3_2 = 1.2247448713915890490986420373529

Constant expression for \(\sqrt{3/2}\).

double sqrt2 = 1.4142135623730950488016887242097

Constant expression for \(\sqrt{2}\).

double sqrt3 = 1.7320508075688772935274463415059

Constant expression for \(\sqrt{3}\).

double sqrt6 = 2.4494897427831780981972840747059

Constant expression for \(\sqrt{6}\).

namespace ContinuumMechanics

namespace CommonTensors

Functions

EigenTensors::Tensor3333d Initialize_I2xI2()

Initializes the fourth-order tensor \(I_{ijkl} = \delta_{ij}\delta_{kl}\).

EigenTensors::Tensor3333d Initialize_Isym()

Initializes the symmetric fourth-order identity tensor \(I_{ijkl}^{sym}=\frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}) \).

EigenTensors::Tensor3333d Initialize_Iskew()

Initializes the skew-symmetric fourth-order identity tensor \( I_{ijkl}^{skew}=\frac{1}{2}(\delta_{ik}\delta_{jl}-\delta_{il}\delta_{jk}) \).

EigenTensors::Tensor3333d Initialize_IFourthOrder()

Initializes the fourth-order identity tensor \( I_{ijkl} = \delta_{ik}\delta_{jl} \).

EigenTensors::Tensor3333d Initialize_IFourthOrderTranspose()

Initializes the transposed fourth-order identity tensor \( I_{ijkl}^{T} = \delta_{il}\delta_{jk} \).

EigenTensors::Tensor3333d Initialize_dDeviatoricStress_dStress()

Initializes the derivative tensor of deviatoric stress w.r.t. stress \( \frac{\partial s_{ij}}{\partial\sigma_{kl}} = \delta_{ik}\delta_{jl} - \frac{1}{3} \delta_{ij}\delta_{kl} \).

EigenTensors::Tensor333d Initialize_LeviCivita3D()

Initializes the 3D Levi-Civita permutation tensor \(E_{ijk}\).

A fully antisymmetric third-order tensor defined as

\[\begin{split} E_{ijk} = \begin{cases} +1 & \text{if } (i,j,k) \text{ is an even permutation of } (1,2,3), \\ -1 & \text{if } (i,j,k) \text{ is an odd permutation of } (1,2,3), \\ 0 & \text{if any two indices are equal.} \end{cases} \end{split}\]

EigenTensors::Tensor122d Initialize_LeviCivita2D()

2D Levi-Civita permutation tensor \(E_{ij}\).

A fully antisymmetric second-order tensor defined as

\[\begin{split} E_{ij} = \begin{cases} +1 & \text{if } (i,j) = (1,2), \\ -1 & \text{if } (i,j) = (2,1), \\ 0 & \text{if } i=j. \end{cases} \end{split}\]

Commonly used to represent 2D cross products and rotations in tensor notation.

EigenTensors::Tensor33d Initialize_I2()

Initializes the second-order identity tensor.

\( I_{ij} = \delta_{ij} \).

int getNumberOfDofForRotation(int nDim)

Returns the number of rotational DOFs for the given dimension.

Parameters:

nDim – Problem dimension (2 or 3).

Returns:

1 in 2D, 3 in 3D.

template<int sizeI, int sizeJ>
Eigen::Matrix<double, sizeI * sizeJ, sizeI * sizeJ> makeIndexSwapTensor()

Constructs an aux. matrix, which helps to swap indices in Eigen::Matrices abused as higher order Tensors by multiplication.

The transformation is defined as:

\[ T_{(ij)(kl)} \cdot P_{(kl)(lk)} = T_{(ij)(lk)} \]

Eigen::Matrix3d t;
t << 1,2,3,
     4,5,6,
     7,8,9;
const auto P = makeIndexSwapTensor<3,3>();
std::cout << t.reshaped().transpose() << std::endl;
std::cout << t.reshaped().transpose() * P << std::endl;

// Output:
// 1,4,7,2,5,8,3,6,9
// 1,2,3,4,5,6,7,8,9

Template Parameters:

• sizeI – Number of rows in the tensor index block.

• sizeJ – Number of columns in the tensor index block.

Returns:

Eigen::Matrix<double, sizeI * sizeJ, sizeI * sizeJ> The index swap tensor.

template<int nDim>
Eigen::TensorFixedSize<double, Eigen::Sizes<getNumberOfDofForRotation(nDim), nDim, nDim>> getReferenceToCorrectLeviCivita()

Provides the reference Levi-Civita permutation tensor for the given spatial dimension.

Template Parameters:

nDim – Spatial dimension (2 or 3).

Returns:

The Levi-Civita tensor:

• In 2D: \(E_{ij}\), a second-order antisymmetric tensor with components A fully antisymmetric second-order tensor defined as

  \[\begin{split} E_{ij} = \begin{cases} +1 & \text{if } (i,j) = (1,2), \\ -1 & \text{if } (i,j) = (2,1), \\ 0 & \text{if } i=j. \end{cases} \end{split}\]
  Commonly used to represent 2D cross products and rotations in tensor notation.

• In 3D: \(E_{ijk}\), a third-order antisymmetric tensor with components A fully antisymmetric third-order tensor defined as

  \[\begin{split} E_{ijk} = \begin{cases} +1 & \text{if } (i,j,k) \text{ is an even permutation of } (1,2,3), \\ -1 & \text{if } (i,j,k) \text{ is an odd permutation of } (1,2,3), \\ 0 & \text{if any two indices are equal.} \end{cases} \end{split}\]

Variables

const EigenTensors::Tensor3333d I2xI2 = Initialize_I2xI2()

Fourth-order tensor \(I_{ijkl} = \delta_{ij}\delta_{kl}\).

const EigenTensors::Tensor3333d Isym = Initialize_Isym()

Symmetric fourth-order identity tensor \( I_{ijkl}^{sym}=\frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}) \).

const EigenTensors::Tensor3333d Iskew = Initialize_Iskew()

Skew-symmetric part of the fourth-order identity tensor \(I_{ijkl}^{skew}=\frac{1}{2}(\delta_{ik}\delta_{jl}-\delta_{il}\delta_{jk}) \).

const EigenTensors::Tensor3333d IFourthOrder = Initialize_IFourthOrder()

Fourth-order identity tensor \( I_{ijkl} = \delta_{ik}\delta_{jl} \).

const EigenTensors::Tensor3333d IFourthOrderTranspose = Initialize_IFourthOrderTranspose()

Transposed fourth-order identity tensor \( I_{ijkl}^{T} = \delta_{il}\delta_{jk} \).

const EigenTensors::Tensor3333d dDeviatoricStress_dStress = Initialize_dDeviatoricStress_dStress()

Derivative of the deviatoric stress with respect to stress \( \frac{\partial s_{ij}}{\partial\sigma_{kl}} = \delta_{ik}\delta_{jl} - \frac{1}{3} \delta_{ij}\delta_{kl} \).

const EigenTensors::Tensor333d LeviCivita3D = Initialize_LeviCivita3D()

3D Levi-Civita permutation tensor \(E_{ijk}\).

A fully antisymmetric third-order tensor defined as

\[\begin{split} E_{ijk} = \begin{cases} +1 & \text{if } (i,j,k) \text{ is an even permutation of } (1,2,3), \\ -1 & \text{if } (i,j,k) \text{ is an odd permutation of } (1,2,3), \\ 0 & \text{if any two indices are equal.} \end{cases} \end{split}\]

const EigenTensors::Tensor122d LeviCivita2D = Initialize_LeviCivita2D()

2D Levi-Civita permutation tensor \(\varepsilon_{ij}\).

A fully antisymmetric second-order tensor defined as

\[\begin{split} E_{ij} = \begin{cases} +1 & \text{if } (i,j) = (1,2), \\ -1 & \text{if } (i,j) = (2,1), \\ 0 & \text{if } i=j. \end{cases} \end{split}\]

Commonly used to represent 2D cross products and rotations in tensor notation.

const EigenTensors::Tensor33d I2 = Initialize_I2()

Second-order identity tensor.

\( I_{ij} = \delta_{ij} \).

namespace TensorUtility

Functions

int d(int a, int b)

Kronecker delta function \( \delta_{ab} \).

Returns:

1 if a == b, otherwise 0.

template<int x, int y, typename T, typename = std::enable_if<!std::is_const<std::remove_reference<T>>::value>>
auto as(T &t)

Map an object’s raw data as a fixed-size Eigen matrix.

Creates an Eigen::Map view of t reinterpreted as an x-by-y matrix. No data is copied; the map directly references the storage returned by t.data().

Note

Caller must ensure t.data() has at least x*y elements and matches Eigen’s default (column-major) layout.

Template Parameters:

• x – Number of rows (compile-time).

• y – Number of columns (compile-time).

• T – Object type, must define Scalar and provide data().

Parameters:

t – Input object.

Returns:

Eigen::Map<Eigen::Matrix<typename T::Scalar, x, y>>.

template<int x, int y, typename T, typename = void>
auto as(const T &t)

Map an object’s raw data as a fixed-size Eigen matrix (const version).

Creates an Eigen::Map view of t reinterpreted as an x-by-y matrix. No data is copied; the map directly references the storage returned by t.data().

Note

Caller must ensure t.data() has at least x*y elements and matches Eigen’s default (column-major) layout.

Template Parameters:

• x – Number of rows (compile-time).

• y – Number of columns (compile-time).

• T – Object type, must define Scalar and provide data().

Parameters:

t – Input object.

Returns:

Eigen::Map<Eigen::Matrix<typename T::Scalar, x, y>>.

template<typename Derived, typename = std::enable_if<!std::is_const<std::remove_reference<Derived>>::value>>
auto flatten(Derived &t)

Flattens an Eigen object into a 1D column vector map.

Creates an Eigen::Map that reinterprets the data of t as a column vector of size RowsAtCompileTime * ColsAtCompileTime. No copy is made; modifications through the map affect t.

Template Parameters:

Derived – Eigen matrix/array type (non-const).

Parameters:

t – Eigen object to be flattened.

Returns:

Eigen::Map<Eigen::Matrix<Scalar, Rows*Cols, 1>> referencing the data of t.

template<typename Derived, typename = void>
auto flatten(const Derived &t)

Flattens an Eigen object into a 1D column vector map (const version).

template<int... Pairs>
auto contractionDims()

Build contraction dimension pairs for Eigen tensor operations.

Generates an Eigen::array<Eigen::IndexPair<int>, N> from the compile-time parameter pack Pairs, where each consecutive pair defines a contraction index mapping. Evaluated entirely at compile time.

Template Parameters:

Pairs – Sequence of integers; must contain an even number of values forming index pairs.

Returns:

Array of index pairs usable in Eigen tensor contraction.

Eigen::Matrix3d dyadicProduct(const Eigen::Vector3d &vector1, const Eigen::Vector3d &vector2)

Compute the dyadic (outer) product of two 3D vectors.

Forms a 3x3 matrix where each entry is given by \( a_{ij} = b_i c_j \).

Parameters:

• vector1 – First 3D vector.

• vector2 – Second 3D vector.

Returns:

3x3 matrix representing the dyadic product.

namespace IndexNotation

Functions

template<int nDim>
std::pair<int, int> fromVoigt(int ij)

Convert a Voigt index to tensor indices.

Maps a Voigt notation index ij to the corresponding (i, j) tensor indices for a given dimension nDim.

Template Parameters:

nDim – Problem dimension (1, 2, or 3).

Parameters:

ij – Voigt index.

Throws:

std::invalid_argument – if nDim or ij is invalid.

Returns:

Pair of tensor indices (i, j).

template<int nDim>
int toVoigt(int i, int j)

Maps tensor indices (i, j) to the corresponding Voigt notation index for a given dimension nDim.

Template Parameters:

nDim – Problem dimension (1, 2, or 3).

Parameters:

• i – Row index of the tensor.

• j – Column index of the tensor.

Throws:

std::invalid_argument – if nDim is invalid.

Returns:

Voigt index corresponding to (i, j).

template<int nDim>
Eigen::TensorFixedSize<double, Eigen::Sizes<VOIGTFROMDIM(nDim), nDim, nDim>> voigtMap()

Construct the Voigt mapping tensor.

Creates a 3rd-order tensor that maps tensor indices (i, j) to their corresponding Voigt index. Each entry is 1 at (toVoigt<nDim>(i, j), i, j) and 0 elsewhere.

Template Parameters:

nDim – Problem dimension (1, 2, or 3).

Returns:

A tensor of shape (VoigtSize, nDim, nDim) encoding the Voigt mapping.

namespace TensorExponential

Functions

template<typename T, size_t tensorSize>
Fastor::Tensor<T, tensorSize, tensorSize> computeTensorExponential(const Fastor::Tensor<T, tensorSize, tensorSize> &theTensor, int maxIterations, double tolerance, double alternativeTolerance)

Computes the exponential of a square second rank tensor.

Evaluates exp(theTensor) using a series expansion up to maxIterations or until the terms converge below tolerance.

Template Parameters:

• T – Element type of the tensor.

• tensorSize – Dimension of the square tensor.

Parameters:

• theTensor – Input tensor.

• maxIterations – Maximum number of series terms to compute.

• tolerance – Convergence threshold for the series terms.

• alternativeTolerance – Secondary threshold to detect failure.

Throws:

ExponentialMapFailed – If convergence is not achieved.

Returns:

Tensor representing the exponential of theTensor.

namespace FirstOrderDerived

Functions

template<typename T, size_t tensorSize>
TensorExponentialResult<T, tensorSize> computeTensorExponential(const Fastor::Tensor<T, tensorSize, tensorSize> &theTensor, int maxIterations, double tolerance)

Computes the exponential of a square second rank tensor and its first-order derivative.

See the documentation of computeTensorExponential without derivative for additional details.

namespace EigenTensors

Typedefs

typedef Eigen::TensorFixedSize<double, Eigen::Sizes<6, 3, 3>> Tensor633d

typedef Eigen::TensorFixedSize<double, Eigen::Sizes<3, 2, 2>> Tensor322d

typedef Eigen::TensorFixedSize<double, Eigen::Sizes<3, 3, 3, 3>> Tensor3333d

typedef Eigen::TensorFixedSize<double, Eigen::Sizes<3, 3, 3>> Tensor333d

typedef Eigen::TensorFixedSize<double, Eigen::Sizes<1, 2, 2>> Tensor122d

typedef Eigen::TensorFixedSize<double, Eigen::Sizes<2, 2, 2, 2>> Tensor2222d

typedef Eigen::TensorFixedSize<double, Eigen::Sizes<2, 2, 1, 2>> Tensor2212d

typedef Eigen::TensorFixedSize<double, Eigen::Sizes<2, 1, 2, 2>> Tensor2122d

typedef Eigen::TensorFixedSize<double, Eigen::Sizes<2, 1, 1, 2>> Tensor2112d

typedef Eigen::TensorFixedSize<double, Eigen::Sizes<3, 3>> Tensor33d

namespace FastorIndices

Typedefs

using A = Fastor::Index<A_>

using Ai = Fastor::Index<A_, i_>

using AB = Fastor::Index<A_, B_>

using B = Fastor::Index<B_>

using IJ = Fastor::Index<I_, J_>

using IJKL = Fastor::Index<I_, J_, K_, L_>

using IJML = Fastor::Index<I_, J_, M_, L_>

using IK = Fastor::Index<I_, K_>

using Ii = Fastor::Index<I_, i_>

using IikK = Fastor::Index<I_, i_, k_, K_>

using Ik = Fastor::Index<I_, k_>

using IL = Fastor::Index<I_, L_>

using Im = Fastor::Index<I_, m_>

using JI = Fastor::Index<J_, I_>

using JK = Fastor::Index<J_, K_>

using JL = Fastor::Index<J_, L_>

using Jk = Fastor::Index<J_, k_>

using KI = Fastor::Index<K_, I_>

using KJ = Fastor::Index<K_, J_>

using KJN = Fastor::Index<K_, J_, N_>

using KL = Fastor::Index<K_, L_>

using KLMN = Fastor::Index<K_, L_, M_, N_>

using KLMP = Fastor::Index<K_, L_, M_, P_>

using KLNM = Fastor::Index<K_, L_, N_, M_>

using KLPM = Fastor::Index<K_, L_, P_, M_>

using KLm = Fastor::Index<K_, L_, m_>

using KMJ = Fastor::Index<K_, M_, J_>

using KMN = Fastor::Index<K_, M_, N_>

using Ki = Fastor::Index<K_, i_>

using Kk = Fastor::Index<K_, k_>

using L = Fastor::Index<L_>

using LI = Fastor::Index<L_, I_>

using LK = Fastor::Index<L_, K_>

using LN = Fastor::Index<L_, N_>

using Lm = Fastor::Index<L_, m_>

using LmN = Fastor::Index<L_, m_, N_>

using MJKL = Fastor::Index<M_, J_, K_, L_>

using MK = Fastor::Index<M_, K_>

using ML = Fastor::Index<M_, L_>

using MNL = Fastor::Index<M_, N_, L_>

using MPm = Fastor::Index<M_, P_, m_>

using Mi = Fastor::Index<M_, i_>

using NLJl = Fastor::Index<N_, L_, J_, l_>

using Nm = Fastor::Index<N_, m_>

using Pm = Fastor::Index<P_, m_>

using i = Fastor::Index<i_>

using iA = Fastor::Index<i_, A_>

using iAkB = Fastor::Index<i_, A_, k_, B_>

using iB = Fastor::Index<i_, B_>

using iI = Fastor::Index<i_, I_>

using iIKL = Fastor::Index<i_, I_, K_, L_>

using iIjJ = Fastor::Index<i_, I_, j_, J_>

using iIkK = Fastor::Index<i_, I_, k_, K_>

using iIkL = Fastor::Index<i_, I_, k_, L_>

using iImn = Fastor::Index<i_, I_, m_, n_>

using iJ = Fastor::Index<i_, J_>

using iJKL = Fastor::Index<i_, J_, K_, L_>

using iJLl = Fastor::Index<i_, J_, L_, l_>

using iJkL = Fastor::Index<i_, J_, k_, L_>

using iJl = Fastor::Index<i_, J_, l_>

using iK = Fastor::Index<i_, K_>

using iKjL = Fastor::Index<i_, K_, j_, L_>

using iL = Fastor::Index<i_, L_>

using iM = Fastor::Index<i_, M_>

using iN = Fastor::Index<i_, N_>

using iNL = Fastor::Index<i_, N_, L_>

using ij = Fastor::Index<i_, j_>

using ijB = Fastor::Index<i_, j_, B_>

using ijKJ = Fastor::Index<i_, j_, K_, J_>

using ijKL = Fastor::Index<i_, j_, K_, L_>

using ijL = Fastor::Index<i_, j_, L_>

using ijLm = Fastor::Index<i_, j_, L_, m_>

using ijk = Fastor::Index<i_, j_, k_>

using ijkB = Fastor::Index<i_, j_, k_, B_>

using ijkK = Fastor::Index<i_, j_, k_, K_>

using ijkl = Fastor::Index<i_, j_, k_, l_>

using ijl = Fastor::Index<i_, j_, l_>

using ijm = Fastor::Index<i_, j_, m_>

using ijmn = Fastor::Index<i_, j_, m_, n_>

using ijnk = Fastor::Index<i_, j_, n_, k_>

using ijnm = Fastor::Index<i_, j_, n_, m_>

using ik = Fastor::Index<i_, k_>

using im = Fastor::Index<i_, m_>

using imk = Fastor::Index<i_, m_, k_>

using imkl = Fastor::Index<i_, m_, k_, l_>

using imL = Fastor::Index<i_, m_, L_>

using imLk = Fastor::Index<i_, m_, L_, k_>

using in = Fastor::Index<i_, n_>

using inB = Fastor::Index<i_, n_, B_>

using inkB = Fastor::Index<i_, n_, k_, B_>

using j = Fastor::Index<j_>

using jA = Fastor::Index<j_, A_>

using jB = Fastor::Index<j_, B_>

using jJ = Fastor::Index<j_, J_>

using jL = Fastor::Index<j_, L_>

using jLm = Fastor::Index<j_, L_, m_>

using ji = Fastor::Index<j_, i_>

using jin = Fastor::Index<j_, i_, n_>

using jk = Fastor::Index<j_, k_>

using jK = Fastor::Index<j_, K_>

using jkB = Fastor::Index<j_, k_, B_>

using jkl = Fastor::Index<j_, k_, l_>

using jl = Fastor::Index<j_, l_>

using k = Fastor::Index<k_>

using kA = Fastor::Index<k_, A_>

using kB = Fastor::Index<k_, B_>

using kI = Fastor::Index<k_, I_>

using kK = Fastor::Index<k_, K_>

using kL = Fastor::Index<k_, L_>

using kM = Fastor::Index<k_, M_>

using kNL = Fastor::Index<k_, N_, L_>

using kj = Fastor::Index<k_, j_>

using kJ = Fastor::Index<k_, J_>

using kl = Fastor::Index<k_, l_>

using km = Fastor::Index<k_, m_>

using l = Fastor::Index<l_>

using lB = Fastor::Index<l_, B_>

using lm = Fastor::Index<l_, m_>

using m = Fastor::Index<m_>

using mK = Fastor::Index<m_, K_>

using mLl = Fastor::Index<m_, L_, l_>

using mj = Fastor::Index<m_, j_>

using mjL = Fastor::Index<m_, j_, L_>

using mn = Fastor::Index<m_, n_>

using mnKL = Fastor::Index<m_, n_, K_, L_>

using mnij = Fastor::Index<m_, n_, i_, j_>

using mnkB = Fastor::Index<m_, n_, k_, B_>

using mnkL = Fastor::Index<m_, n_, k_, L_>

using nB = Fastor::Index<n_, B_>

using to_IJKL = Fastor::OIndex<I_, J_, K_, L_>

using to_IJkK = Fastor::OIndex<I_, J_, k_, K_>

using to_IJkL = Fastor::OIndex<I_, J_, k_, L_>

using to_IikK = Fastor::OIndex<I_, i_, k_, K_>

using to_IjkK = Fastor::OIndex<I_, j_, k_, K_>

using to_Ii = Fastor::OIndex<I_, i_>

using to_NLJl = Fastor::OIndex<N_, L_, J_, l_>

using to_iIKL = Fastor::OIndex<i_, I_, K_, L_>

using to_iIjJ = Fastor::OIndex<i_, I_, j_, J_>

using to_iImn = Fastor::OIndex<i_, I_, m_, n_>

using to_ij = Fastor::OIndex<i_, j_>

using to_ijIJ = Fastor::OIndex<i_, j_, I_, J_>

using to_ijKL = Fastor::OIndex<i_, j_, K_, L_>

using to_ijL = Fastor::OIndex<i_, j_, L_>

using to_ijLk = Fastor::OIndex<i_, j_, L_, k_>

using to_ijLm = Fastor::OIndex<i_, j_, L_, m_>

using to_ijk = Fastor::OIndex<i_, j_, k_>

using to_ijkK = Fastor::OIndex<i_, j_, k_, K_>

using to_ijkL = Fastor::OIndex<i_, j_, k_, L_>

using to_ijKl = Fastor::OIndex<i_, j_, K_, l_>

using to_ijkl = Fastor::OIndex<i_, j_, k_, l_>

using to_ijm = Fastor::OIndex<i_, j_, m_>

using to_ijmM = Fastor::OIndex<i_, j_, m_, M_>

using to_jAB = Fastor::OIndex<j_, A_, B_>

using to_jAkB = Fastor::OIndex<j_, A_, k_, B_>

using to_ji = Fastor::OIndex<j_, i_>

using to_jikL = Fastor::OIndex<j_, i_, k_, L_>

using to_jikl = Fastor::OIndex<j_, i_, k_, l_>

using to_jkiB = Fastor::OIndex<j_, k_, i_, B_>

using to_kK = Fastor::OIndex<k_, K_>

using to_kL = Fastor::OIndex<k_, L_>

Enums

Values:

enumerator i_

enumerator j_

enumerator k_

enumerator l_

enumerator m_

enumerator n_

enumerator A_

enumerator B_

enumerator I_

enumerator J_

enumerator K_

enumerator L_

enumerator M_

enumerator N_

enumerator P_

namespace FastorStandardTensors

Typedefs

using Tensor3d = Fastor::Tensor<double, 3>

using Tensor33d = Fastor::Tensor<double, 3, 3>

using Tensor333d = Fastor::Tensor<double, 3, 3, 3>

using Tensor3333d = Fastor::Tensor<double, 3, 3, 3, 3>

template<typename T>
using Tensor3t = Fastor::Tensor<T, 3>

template<typename T>
using Tensor33t = Fastor::Tensor<T, 3, 3>

template<typename T>
using Tensor333t = Fastor::Tensor<T, 3, 3, 3>

template<typename T>
using Tensor3333t = Fastor::Tensor<T, 3, 3, 3, 3>

using TensorMap3d = Fastor::TensorMap<double, 3>

using TensorMap33d = Fastor::TensorMap<double, 3, 3>

using TensorMap333d = Fastor::TensorMap<double, 3, 3, 3>

using TensorMap3333d = Fastor::TensorMap<double, 3, 3, 3, 3>

namespace Spatial3D

Variables

const Tensor33d I = Tensor33d((Eigen::Matrix3d() << Eigen::Matrix3d::Identity()).finished().data(), Fastor::ColumnMajor)

const Tensor333d LeviCivita = Tensor333d(Marmot::ContinuumMechanics::CommonTensors::LeviCivita3D.data(), Fastor::ColumnMajor)

const Tensor3333d IHyd = Tensor3333d(Marmot::ContinuumMechanics::CommonTensors::I2xI2.data(), Fastor::ColumnMajor)

const Tensor3333d ISymm = Tensor3333d(Marmot::ContinuumMechanics::CommonTensors::Isym.data(), Fastor::ColumnMajor)

const Tensor3333d ISkew = Tensor3333d(Marmot::ContinuumMechanics::CommonTensors::Iskew.data(), Fastor::ColumnMajor)

const Tensor3333d I4 = Tensor3333d(Marmot::ContinuumMechanics::CommonTensors::IFourthOrder.data(), Fastor::ColumnMajor)

const Tensor3333d ITranspose = Tensor3333d(Marmot::ContinuumMechanics::CommonTensors::IFourthOrderTranspose.data(), Fastor::ColumnMajor)

const Tensor3333d Deviatoric = I4 - 1. / 3 * IHyd

const Tensor3333d DeviatoricTranspose = Fastor::transpose(DeviatoricTranspose)

const Tensor3333d DeviatoricSymmetric = ISymm - 1. / 3 * IHyd

namespace Math

Functions

template<typename T>
bool isNaN(T x)

Checks if a scalar is NaN.

Checks if value x is NaN (not a number). This is done by checking if x is unequal to itself, which is only true for NaN values.

Parameters:

x – scalar to be checked

Returns:

true if x is NaN, false otherwise

double linearInterpolation(double x, double x0, double x1, double y0, double y1)

Performs linear interpolation.

Performs linear interpolation to find the value at point x given two known points (x0, y0) and (x1, y1). If x is outside the range [x0, x1], the function will perform linear extrapolation.

Parameters:

• x – Point at which to interpolate

• x0 – First known point

• x1 – Second known point

• y0 – Value at first known point

• y1 – Value at second known point

Returns:

Interpolated value at point x

double exp(double x)

Computes the exponential of value x with numerical limits check.

If x is larger than the maximum limit of double precision floating point numbers, the maximum limit is returned. If x is smaller than the minimum limit, the minimum limit is returned.

Parameters:

x – Exponent to which e is raised

Returns:

exponential of x

unsigned long factorial(unsigned int n)

Runtime factorial.

Parameters:

n – An unsigned integer

Returns:

The factorial of n

int getExponentPowerTen(const double x)

Extracts the exponent to the power of ten from a floating point number.

This function extracts the exponent to the power of ten from a given floating point number. For example, for an input of 3e5, the function will return 5. If the input number is very close to zero (between -1e-16 and 1e-16), the function returns 0.

Parameters:

x – Input floating point number

Returns:

Exponent to the power of ten

double radToDeg(const double alpha)

Convert angle from radiant to degree.

Converts angle alpha given in radiant to degree.

Parameters:

alpha – Angle in radiant

Returns:

Angle in degree

double degToRad(const double alpha)

Convert angle from degree to radiant.

Converts angle alpha given in degree to radiant.

Parameters:

alpha – Angle in degree

Returns:

Angle in radiant

double macauly(double scalar)

Macaulay function applied to a scalar.

The Macaulay function, also as positive part operator, is defined as:

\[\begin{split}\langle x \rangle = \begin{cases} x, & x \geq 0 \\ 0, & x < 0 \end{cases} \end{split}\]

Parameters:

scalar – Input value

Returns:

Positive part of scalar

int heaviside(double scalar)

Heaviside step function applied to a scalar.

The Heaviside step function is defined as:

\[\begin{split}H(x) = \begin{cases} 1, & x \geq 0 \\ 0, & x < 0 \end{cases} \end{split}\]

Parameters:

scalar – Input value

Returns:

1 if scalar >= 0, 0 otherwise

int heavisideExclude0(double scalar)

Heaviside step function excluding zero applied to a scalar.

The Heaviside step function excluding zero is defined as:

\[\begin{split}H(x) = \begin{cases} 1, & x > 0 \\ 0, & x \leq 0 \end{cases} \end{split}\]

Parameters:

scalar – Input value

Returns:

1 if scalar > 0, 0 otherwise

template<typename T>
int sgn(const T &val) noexcept

Sign function applied to a scalar.

The sign function is defined as:

\[\begin{split} \text{sign}(x) = \begin{cases} 1, & x > 0 \\ 0, & x = 0 \\ -1, & x < 0 \end{cases} \end{split}\]

Parameters:

val – Input value

Returns:

1 if val > 0, -1 if val < 0, 0 if val == 0

double makeReal(const double &value)

Converts various scalar types to double precision floating point numbers.

This function converts input values of different numeric types, including:

• double

• std::complex<double>

• autodiff::real

• autodiff::dual

• Eigen::Matrix with elements of any type convertible to double

The function is overloaded and templated to handle these different types appropriately.

Parameters:

value – Input value of various numeric types

Returns:

Converted value as double

double makeReal(const std::complex<double> &value)

double makeReal(const autodiff::real &value)

template<typename T, typename G>
double makeReal(const autodiff::detail::Dual<T, G> &number)

Converts autodiff::dual numbers to double precision floating point numbers.

This function converts an autodiff::dual number to a double precision floating point number. It is templated to handle dual numbers with different underlying types.

Template Parameters:

T – Underlying type of the autodiff::dual number

Parameters:

• G – Gradient type of the autodiff::dual number

• number – Input autodiff::dual number

Returns:

Converted value as double

template<typename T, int... Rest>
Eigen::Matrix<double, Rest...> makeReal(const Eigen::Matrix<T, Rest...> mat)

Extracts the real part of a arbitrary scalartype-valued Matrix.

Template Parameters:

• T – scalar type

• Rest... – parameter pack for additional matrix information, e.g, dimensions

Parameters:

mat – T-valued matrix

Returns:

double-valued matrix

template<typename T>
Eigen::VectorXd makeReal(Eigen::Vector<T, Eigen::Dynamic> in)

Extracts the real part of a arbitrary scalartype-valued Vector.

Template Parameters:

T – scalar type

Parameters:

in – T-valued vector

Returns:

double-valued vector

template<int nRows, int nCols>
Eigen::Matrix<double, nRows, nCols> macaulyMatrix(const Eigen::Matrix<double, nRows, nCols> &mat)

Apply Macaulay function to a matrix.

Applies the Macaulay function element-wise to the input matrix.

Todo:

: Can be replaced easily with Eigen’s array() functionality ???

Template Parameters:

• nRows – Number of rows in the matrix

• nCols – Number of columns in the matrix

Parameters:

mat – Input matrix

Returns:

Matrix with Macaulay function applied element-wise

template<typename functionType, typename yType, typename ...Args>
yType explicitEuler(yType yN, const double dt, functionType fRate, Args&&... fRateArgs)

Explicit Euler time integrator.

This function computes one single time step using the explicit Euler method:

\[ \boldsymbol{y}_{n+1} = \boldsymbol{y}_n + f(\boldsymbol{y}_n) * \Delta t \]

where \( \boldsymbol{y}_n \) is the current value, \( \Delta t \) is the time step size, and \( f(\boldsymbol{y}_n) \) is the rate of change.

Parameters:

• yN – Current value

• dt – Time step size

• fRate – Function that computes the rate of change

• fRateArgs – Additional arguments for the rate function

Returns:

Updated value after time step

template<int ySize, typename functionType, typename ...Args>
Eigen::Matrix<double, ySize, 1> semiImplicitEuler(Eigen::Matrix<double, ySize, 1> yN, const double dt, functionType fRate, Args&&... fRateArgs)

Implicit Euler time integrator for a linear vector-valued rate equation.

This function computes one single time step using the implicit Euler method for a linear vector-valued rate equation:

\[ \boldsymbol{y}_{n+1} = \boldsymbol{y}_n + f(\boldsymbol{y}_{n+1}) \Delta t \]

where \(\boldsymbol{y}_n \) is the current value, \( \Delta t \) is the time step size, and \( f(\boldsymbol{y}) \) is the linear rate of change evaluated at the next time step.

Todo:

: Is this really semi-implicit? It looks like a fully implicit Euler step for linear systems.

: Use external central difference function?

Parameters:

• yN – Current value

• dt – Time step size

• fRate – Function that computes the rate of change

• fRateArgs – Additional arguments for the rate function

Returns:

Updated value after time step

template<typename functionType, typename yType, typename ...Args>
yType explicitEulerRichardson(yType yN, const double dt, functionType fRate, Args&&... fRateArgs)

Explicit Euler integration with error estimation based on Richardson extrapolation.

This function computes one single time step using the explicit Euler method with Richardson extrapolation for error estimation:

\[ \boldsymbol{y}_{n+1} = 2 \left( \boldsymbol{y}_n + f\left(\boldsymbol{y}_n + \frac{f(\boldsymbol{y}_n) \Delta t}{2}\right) \frac{\Delta t}{2} \right) - \left( \boldsymbol{y}_n + f(\boldsymbol{y}_n) \Delta t \right) \]

where \( \boldsymbol{y}_n \) is the current value, \( \Delta t \) is the time step size, and \( f(\boldsymbol{y}) \) is the rate of change.

Template Parameters:

• functionType – Type of the rate function

• yType – Type of the current value

• Args – Additional argument types for the rate function

Parameters:

• yN – Current value

• dt – Time step size

• fRate – Function that computes the rate of change

• fRateArgs – Additional arguments for the rate function

Returns:

Updated value after time step

template<int ySize, typename functionType, typename ...Args>
std::tuple<Eigen::Matrix<double, ySize, 1>, double> explicitEulerRichardsonWithErrorEstimator(Eigen::Matrix<double, ySize, 1> yN, const double dt, const double TOL, functionType fRate, Args&&... fRateArgs)

Explicit Euler integration based on Richardson extrapolation with error estimation and time step estimation.

This function computes one single time step using the explicit Euler method with Richardson extrapolation for error estimation and adaptive time stepping:

\[ \boldsymbol{y}_{n+1} = 2 \left( \boldsymbol{y}_n + f\left(\boldsymbol{y}_n + \frac{f(\boldsymbol{y}_n) \Delta t}{2}\right) \frac{\Delta t}{2} \right) - \left( \boldsymbol{y}_n + f(\boldsymbol{y}_n) \Delta t \right) \]

where \( \boldsymbol{y}_n \) is the current value, \( \Delta t \) is the current time step size, and \( f(\boldsymbol{y}) \) is the rate of change.

The function also estimates the error of the time step and adjusts the time step size for the next iteration based on the desired tolerance TOL. The new time step size is computed as:

\[ \Delta t_{\text{new}} = \Delta t \cdot \min\left(2, \max\left(0.2, 0.9 \sqrt{\frac{TOL}{EST}}\right)\right) \]

where \( EST \) is the estimated error.

Todo:

: reuse explicitEulerRichardson function?

Template Parameters:

• ySize – Size of the state vector

• yType – Type of the current value

• Args – Additional argument types for the rate function

Parameters:

• functionType – Type of the rate function

• Args – Additional argument types for the rate function

• yN – Current value

• dt – Current time step size

• TOL – Desired tolerance for the error estimation

• fRate – Function that computes the rate of change

• fRateArgs – Additional arguments for the rate function

Returns:

Tuple containing the updated value after time step and the new time step size

Matrix3d directionCosines(const Matrix3d &transformedCoordinateSystem)

Computes the direction cosines between a given coordinate system and the global coordinate system.

The direction cosines are computed as the dot products between the basis vectors of the transformed coordinate system and the global coordinate system.

Parameters:

transformedCoordinateSystem – 3x3 matrix representing the transformed coordinate system

Returns:

3x3 matrix of direction cosines

Matrix3d orthonormalCoordinateSystem(Vector3d &normalVector)

Constructs a orthonormal coordinate system with a given normal vector as \(x_1\)-axis.

The orthonormal coordinate system is constructed such that:

• The first column corresponds to the normalized input normal vector ( \(x_1\)-axis).

• The second column is a unit vector orthogonal to the first ( \(x_2\)-axis).

• The third column is the cross product of the first two columns ( \(x_3\)-axis).

Todo:

: Maybe remove this function completely to avoid mistakes?

Note

Do not use this function if you want to control the direction of the axes in the plane orthogonal to the normal vector.

Parameters:

normalVector – Input normal vector, which will be normalized and used as \(x_1\)-axis.

Returns:

Orthonormal coordinate system as a 3x3 matrix.

Matrix3d orthonormalCoordinateSystem(const Vector3d &n1, const Vector3d &n2)

Constructs an orthonormal coordinate system with two given normal vectors.

The orthonormal coordinate system is constructed such that:

• The first column corresponds to the normalized first input normal vector ( \(x_1\)-axis).

• The second colum corresponds to the normalized second input normal vector ( \(x_2\)-axis).

• The third column is the normalized cross product of the first two columns ( \(x_3\)-axis).

Parameters:

• n1 – First input vector.

• n2 – Second input vector.

Throws:

std::invalid_argument – if the input normal vectors are not orthogonal.

Returns:

Orthonormal coordinate system as a 3x3 matrix.

Matrix3d transformToLocalSystem(const Matrix3d &T, const Matrix3d &transformedCoordinateSystem)

Transforms a second-order tensor from the global coordinate system to a local coordinate system.

The transformation is performed as: \( T_{local} = Q^T \, T_{global} \, Q \), where \( Q \) is the transformation matrix.

Parameters:

• T – The tensor to be transformed, represented as a 3x3 matrix in the global coordinate system.

• transformedCoordinateSystem – The 3x3 transformation matrix representing the local coordinate system axes in global coordinates.

Returns:

The tensor represented in the local coordinate system as a 3x3 matrix.

Matrix3d transformToGlobalSystem(const Matrix3d &T, const Matrix3d &transformedCoordinateSystem)

Transforms a second-order tensor from a local coordinate system to the global coordinate system.

The transformation is performed as: \( T_{global} = Q \, T_{local} \, Q^T \), where \( Q \) is the transformation matrix.

Parameters:

• T – The tensor to be transformed, represented as a 3x3 matrix in the local coordinate system.

• transformedCoordinateSystem – The 3x3 transformation matrix representing the local coordinate system axes in global coordinates.

Returns:

The tensor represented in the global coordinate system as a 3x3 matrix.

namespace NumericalAlgorithms

namespace Differentiation

Typedefs

using scalar_to_scalar_function_type = std::function<double(const double x)>

Type definition for a function that maps a double to another double.

using vector_to_vector_function_type = std::function<Eigen::VectorXd(const Eigen::VectorXd &X)>

Type definition for a function that maps an Eigen vector to another Eigen vector.

Functions

double forwardDifference(const scalar_to_scalar_function_type &f, const double x)

Approximates the first derivative of a scalar function f at point x using the forward difference method.

It actually computes

\[ f'(x) \approx \frac{f(x + h) - f(x)}{h} \]

with \( h = \sqrt{\epsilon}\, \max(1,|x|) \) and \( \epsilon \) being the machine precision.

Parameters:

• f – The scalar function for which the derivative is to be approximated.

• x – The point at which the derivative is to be approximated.

Returns:

The approximated first derivative of f at point x.

double centralDifference(const scalar_to_scalar_function_type &f, const double x)

Approximates the first derivative of a scalar function f at point x using the central difference method.

It actually approximates the first derivative by

\[ f'(x) \approx \frac{f(x + h) - f(x - h)}{2h} \]

with \( h = \sqrt[3]{\epsilon}\,\max( 1, |x|)\) and \( \epsilon \) being the machine precision.

Parameters:

• f – The scalar function for which the derivative is to be approximated.

• x – The point at which the derivative is to be approximated.

Returns:

The approximated first derivative of f at point x.

Eigen::MatrixXd forwardDifference(const vector_to_vector_function_type &F, const Eigen::VectorXd &X)

Approximates the Jacobian matrix of a vector function F at point X using the forward difference method.

It actually computes

\[ J_{ij} \approx \frac{F_i(X + h e_j) - F_i(X)}{h} \]

with \( h = \sqrt{\epsilon} (1 + |X_j|) \), \( e_j \) being the unit vector in the j-th direction, and \(\epsilon \) being the machine precision.

Parameters:

• F – The vector function for which the Jacobian matrix is to be approximated.

• X – The point at which the Jacobian matrix is to be approximated.

Returns:

The approximated Jacobian matrix of F at point X.

Eigen::MatrixXd centralDifference(const vector_to_vector_function_type &F, const Eigen::VectorXd &X)

Approximates the Jacobian matrix of a vector function F at point X using the central difference method.

It actually approximates the Jacobian matrix by

\[ J_{ij} \approx \frac{F_i(\boldsymbol{X} + h e_j) - F_i(\boldsymbol{X} - h e_j)}{2h} \]

with \( h = \sqrt[3]{\epsilon} \max(1, ||\boldsymbol{X}||) \), \( e_j \) being the unit vector in the j-th direction, and \( \epsilon \) being the machine precision.

Parameters:

• F – The vector function for which the Jacobian matrix is to be approximated.

• X – The point at which the Jacobian matrix is to be approximated.

Returns:

The approximated Jacobian matrix of F at point X.

namespace Complex

Typedefs

using complexDouble = std::complex<double>

template<size_t dim>
using tensor_to_scalar_function_type = std::function<std::complex<double>(const Fastor::Tensor<std::complex<double>, dim, dim> &T)>

Type alias for a function that maps a tensor to a scalar using complex numbers.

Template Parameters:

dim – The dimensions of the input Tensor

using scalar_to_scalar_function_type = std::function<complexDouble(const complexDouble x)>

Type definition for a function that maps a complex double to another complex double.

using vector_to_vector_function_type = std::function<Eigen::VectorXcd(const Eigen::VectorXcd &X)>

Type definition for a function that maps an Eigen complex vector to another Eigen complex vector.

Functions

template<size_t dim>
Fastor::Tensor<double, dim, dim> forwardDifference(const tensor_to_scalar_function_type<dim> &F, const Fastor::Tensor<double, dim, dim> &T)

Approximates the derivative of a function mapping a tensor to a scalar using the complex step method.

Template Parameters:

dim – The dimension of the input Tensor (assumed to be square)

Parameters:

• F – The function mapping a tensor to a scalar with complex numbers

• T – The point at which the derivative is evaluated

Returns:

The derivative of the function F at the point T

double forwardDifference(const scalar_to_scalar_function_type &f, const double x)

Approximates the first derivative of a scalar function f at point x using the complex step method.

It actually computes

\[ f'(x) \approx \frac{\text{Im}(f(x + ih))}{h} \]

with \( h = 10^{-20} \).

Note

This method is highly accurate and does not suffer from subtractive cancellation errors, making it suitable for functions where high precision is required.

Parameters:

• f – The scalar function for which the derivative is to be approximated.

• x – The point at which the derivative is to be approximated.

Returns:

The approximated first derivative of f at point x.

std::tuple<Eigen::VectorXd, Eigen::MatrixXd> forwardDifference(const vector_to_vector_function_type &F, const Eigen::VectorXd &X)

Approximates the Jacobian matrix of a vector function F at point X using the complex step method.

It actually computes

\[ J_{ij} \approx \frac{\text{Im}(F_i(X + ih e_j))}{h} \]

with \( h = 10^{-20} \), and \( e_j \) being the unit vector in the j-th direction.

Note

This method is highly accurate and does not suffer from subtractive cancellation errors, making it suitable for functions where high precision is required.

Parameters:

• F – The vector function for which the Jacobian matrix is to be approximated.

• X – The point at which the Jacobian matrix is to be approximated.

Returns:

A tuple containing the function value at X and the approximated Jacobian matrix of F at point X.

Eigen::MatrixXd centralDifference(const vector_to_vector_function_type &F, const Eigen::VectorXd &X)

Approximates the Jacobian matrix of a vector function F at point X using the complex step method with central differences.

It actually approximates the Jacobian matrix by

\[ J_{ij} \approx \frac{\text{Im}(F_i(\boldsymbol{X} + Ih e_j) - F_i(\boldsymbol{X} - Ih e_j))}{2h} \]

with \( I = \sqrt{2}/2( 1+1i )\), \( h = \sqrt[3]{\epsilon} \max(1, ||\boldsymbol{X}||) \), and \( e_j \) being the unit vector in the j-th direction.

Note

The formula can be found in Lai et al. (2005) Equ. 19.

Parameters:

• F – The vector function for which the Jacobian matrix is to be approximated.

• X – The point at which the Jacobian matrix is to be approximated.

Returns:

The approximated Jacobian matrix of F at point X.

Eigen::MatrixXd fourthOrderAccurateDerivative(const vector_to_vector_function_type &F, const Eigen::VectorXd &X)

Approximates the Jacobian matrix of a vector function F at point X using a fourth-order accurate complex step method.

It actually approximates the Jacobian matrix by

\[J_{ij} \approx \frac{\text{Im}(8 [ F_i(\boldsymbol{X} + I/2h e_j) - F_i(\boldsymbol{X} - I/2h e_j)] - [F_i(\boldsymbol{X} + Ih e_j) + F_i(\boldsymbol{X} - Ih e_j)])}{3\sqrt{2}h} \]

with \( I = \sqrt{2}/2( 1+1i )\), \( h = \sqrt{\epsilon} \max(1, ||\boldsymbol{X}||) \), and \( e_j \) being the unit vector in the j-th direction.

Note

The formula can be found in Lai et al. (2005) Equ. 24.

Parameters:

• F – The vector function for which the Jacobian matrix is to be approximated.

• X – The point at which the Jacobian matrix is to be approximated.

Returns:

The approximated Jacobian matrix of F at point X.

Variables

static const double imaginaryPerturbationSize = 1e-20

static const complexDouble imaginaryPerturbation = imaginaryPerturbationSize * imaginaryUnit

static const std::complex<double> imaginaryUnit = {0, 1}

A variable representing \( 0 + 1i \).

static const std::complex<double> complexUnit = {1, 1}

A variable representing \( 1 + 1i \).

static const std::complex<double> i_ = Marmot::Constants::sqrt2 / 2. * complexUnit

A variable representing \( \frac{\sqrt{2}}{2} (1+1i) \).

namespace ScalarToTensor

Functions

template<size_t... Rest>
std::tuple<Fastor::Tensor<double, Rest...>, Fastor::Tensor<double, Rest...>> forwardDifference(const std::function<Fastor::Tensor<complexDouble, Rest...>(const complexDouble)> &F, const double x)

Approximates the derivative of a function mapping a scalar to a tensor using the complex step method.

Template Parameters:

Rest – The dimensions of the output Tensor

Parameters:

• F – The function mapping a scalar to a tensor with complex numbers

• x – The point at which the derivative is evaluated

Returns:

A tuple containing the real part of F(x) and the derivative of F at the point x

namespace TensorToScalar

Typedefs

template<size_t... Rest>
using tensor_to_scalar_function_type = std::function<complexDouble(const Fastor::Tensor<complexDouble, Rest...> &T)>

Type alias for a function that maps a tensor to a scalar using complex numbers.

Template Parameters:

Rest – The dimensions of the input Tensor

Functions

template<size_t dim>
Fastor::Tensor<double, dim, dim> forwardDifference(const tensor_to_scalar_function_type<dim, dim> &f, const Fastor::Tensor<double, dim, dim> &T)

namespace TensorToTensor

Functions

template<size_t... RestF, size_t... RestT>
Fastor::Tensor<double, RestF..., RestT...> forwardDifference(std::function<Fastor::Tensor<complexDouble, RestF...>(const Fastor::Tensor<complexDouble, RestT...>&)> &F, const Fastor::Tensor<double, RestT...> &T)

Approximates the derivative of a function mapping a tensor to a tensor using the complex step method.

Template Parameters:

• RestF – The dimensions of the output Tensor

• RestT – The dimensions of the input Tensor

Parameters:

• F – The function mapping a tensor to a tensor

• T – The point at which the derivative is evaluated

Returns:

The derivative of the function F at the point T

namespace ScalarToTensor

Functions

template<size_t... Rest>
Fastor::Tensor<double, Rest...> forwardDifference(const std::function<Fastor::Tensor<double, Rest...>(const double)> &F, const double x)

Approximates the derivative of a function mapping a scalar to a tensor using the forward difference method.

Template Parameters:

Rest – The dimensions of the output Tensor

Parameters:

• F – The function mapping a scalar to a tensor

• x – The point at which the derivative is evaluated

Returns:

The derivative of the function F at the point x

template<size_t... Rest>
Fastor::Tensor<double, Rest...> centralDifference(const std::function<Fastor::Tensor<double, Rest...>(const double)> &F, const double x)

Approximates the derivative of a function mapping a scalar to a tensor using the central difference method.

Template Parameters:

Rest – The dimensions of the output Tensor

Parameters:

• F – The function mapping a scalar to a tensor

• x – The point at which the derivative is evaluated

Returns:

The derivative of the function F at the point x

namespace TensorToScalar

Typedefs

template<size_t... Rest>
using tensor_to_scalar_function_type = std::function<double(const Fastor::Tensor<double, Rest...> &T)>

Type alias for a function that maps a tensor to a scalar.

Template Parameters:

Rest – The dimensions of the input Tensor

Functions

template<size_t dim>
Fastor::Tensor<double, dim, dim> forwardDifference(const tensor_to_scalar_function_type<dim, dim> &f, const Fastor::Tensor<double, dim, dim> &T)

Approximates the derivative of a function mapping a tensor to a scalar using the forward difference method.

Template Parameters:

dim – The dimension of the input Tensor (assumed to be square)

Parameters:

• f – The function mapping a tensor to a scalar

• T – The point at which the derivative is evaluated

Returns:

The derivative of the function f at the point T

template<size_t dim>
Fastor::Tensor<double, dim, dim> centralDifference(const tensor_to_scalar_function_type<dim, dim> &F, const Fastor::Tensor<double, dim, dim> &T)

Approximates the derivative of a function mapping a tensor to a scalar using the central difference method.

Template Parameters:

dim – The dimension of the input Tensor (assumed to be square)

Parameters:

• F – The function mapping a tensor to a scalar

• T – The point at which the derivative is evaluated

Returns:

The derivative of the function F at the point T

namespace TensorToTensor

Functions

template<size_t... RestF, size_t... RestT>
Fastor::Tensor<double, RestF..., RestT...> forwardDifference(const std::function<Fastor::Tensor<double, RestF...>(const Fastor::Tensor<double, RestT...>&)> &F, const Fastor::Tensor<double, RestT...> &T)

Approximates the derivative of a function mapping a tensor to a tensor using the forward difference method.

Template Parameters:

• RestF – The dimensions of the output Tensor

• RestT – The dimensions of the input Tensor

Parameters:

• F – The function mapping a tensor to a tensor

• T – The point at which the derivative is evaluated

Returns:

The derivative of the function F at the point T

template<size_t... Rest1, size_t... Rest2>
Fastor::Tensor<double, Rest1..., Rest2...> centralDifference(const std::function<Fastor::Tensor<double, Rest1...>(const Fastor::Tensor<double, Rest2...>&)> &F, const Fastor::Tensor<double, Rest2...> &T)

Approximates the derivative of a function mapping a tensor to a tensor using the central difference method.

Template Parameters:

• Rest1 – The dimensions of the output Tensor

• Rest2 – The dimensions of the input Tensor

Parameters:

• F – The function mapping a tensor to a tensor

• T – The point at which the derivative is evaluated

Returns:

The derivative of the function F at the point T

namespace Integration

Typedefs

using scalar_to_scalar_function_type = std::function<double(const double x)>

Type alias for a function that takes a double and returns a double.

Enums

enum integrationRule

Enumeration of available numerical integration rules.

This enum defines the different numerical integration methods that can be used in the integrateScalarFunction function.

Values:

enumerator midpoint

Midpoint rule for numerical integration.

enumerator trapezodial

Trapezoidal rule for numerical integration.

enumerator simpson

Simpson’s rule for numerical integration.

Functions

double integrateScalarFunction(scalar_to_scalar_function_type f, const std::tuple<double, double> integrationLimits, const int n, const integrationRule intRule)

Numerically integrates a scalar function over a specified interval using a chosen integration rule.

This function approximates the integral of the given scalar function f over the interval defined by integrationLimits using the specified numerical integration rule. The number of subintervals n determines the accuracy of the approximation; a larger n generally leads to a more accurate result.

Example usage:

auto f = [](double x) { return x * x; }; // Function to integrate
double result = integrateScalarFunction(f, std::make_tuple(0.0, 1.0), 100, integrationRule::simpson);

Parameters:

• f – The scalar function to be integrated.

• integrationLimits – A tuple containing the lower and upper limits of integration.

• n – The number of subintervals to use for the integration (must be positive).

• intRule – The numerical integration rule to use (midpoint, trapezoidal, or Simpson’s rule).

Returns:

The approximate value of the integral of f over the specified interval.

namespace std

file MarmotAutomaticDifferentiation.h

#include “autodiff/forward/dual/eigen.hpp”

#include <autodiff/forward/dual/dual.hpp>

#include <functional>

Include dependency graph for MarmotAutomaticDifferentiation.h:

This graph shows which files directly or indirectly include MarmotAutomaticDifferentiation.h:

file MarmotAutomaticDifferentiationForFastor.h

#include “Fastor/Fastor.h”

#include “Marmot/MarmotAutomaticDifferentiation.h”

#include “Marmot/MarmotFastorTensorBasics.h”

#include “autodiff/forward/dual.hpp”

#include “autodiff/forward/dual/eigen.hpp”

#include <autodiff/forward/dual/dual.hpp>

#include <functional>

Include dependency graph for MarmotAutomaticDifferentiationForFastor.h:

file MarmotConstants.h

#include <cmath>

Include dependency graph for MarmotConstants.h:

This graph shows which files directly or indirectly include MarmotConstants.h:

file MarmotFastorTensorBasics.h

#include “Eigen/Core”

#include “Fastor/Fastor.h”

#include “Marmot/MarmotTensor.h”

#include <autodiff/forward/dual/dual.hpp>

Include dependency graph for MarmotFastorTensorBasics.h:

This graph shows which files directly or indirectly include MarmotFastorTensorBasics.h:

file MarmotMath.h

#include “Marmot/MarmotConstants.h”

#include “Marmot/MarmotTypedefs.h”

#include “autodiff/forward/dual.hpp”

#include “autodiff/forward/real.hpp”

#include <algorithm>

#include <autodiff/forward/dual/dual.hpp>

#include <complex>

Include dependency graph for MarmotMath.h:

This graph shows which files directly or indirectly include MarmotMath.h:

file MarmotNumericalDifferentiation.h

#include “Marmot/MarmotMath.h”

#include “Marmot/MarmotTypedefs.h”

#include <functional>

Include dependency graph for MarmotNumericalDifferentiation.h:

This graph shows which files directly or indirectly include MarmotNumericalDifferentiation.h:

file MarmotNumericalDifferentiationForFastor.h

#include “Marmot/MarmotConstants.h”

#include “Marmot/MarmotFastorTensorBasics.h”

#include “Marmot/MarmotNumericalDifferentiation.h”

#include <Fastor/config/macros.h>

#include <Fastor/tensor/Tensor.h>

#include <functional>

Include dependency graph for MarmotNumericalDifferentiationForFastor.h:

file MarmotNumericalIntegration.h

#include <functional>

Include dependency graph for MarmotNumericalIntegration.h:

file MarmotTensor.h

#include “Marmot/MarmotJournal.h”

#include “Marmot/MarmotTypedefs.h”

#include “Marmot/MarmotVoigt.h”

#include <utility>

Include dependency graph for MarmotTensor.h:

This graph shows which files directly or indirectly include MarmotTensor.h:

file MarmotTensorExponential.h

#include “Fastor/Fastor.h”

#include “Marmot/MarmotMath.h”

#include <iostream>

Include dependency graph for MarmotTensorExponential.h:

file MarmotTypedefs.h

#include “Eigen/Core”

#include “Eigen/Dense”

#include “unsupported/Eigen/CXX11/Tensor”

Include dependency graph for MarmotTypedefs.h:

This graph shows which files directly or indirectly include MarmotTypedefs.h:

file NewtonConvergenceChecker.h

#include <Eigen/Core>

Include dependency graph for NewtonConvergenceChecker.h:

dir /home/runner/work/Marmot/Marmot/modules/core

dir /home/runner/work/Marmot/Marmot/modules/core/MarmotMathCore/include

dir /home/runner/work/Marmot/Marmot/modules/core/MarmotMathCore/include/Marmot

dir /home/runner/work/Marmot/Marmot/modules/core/MarmotMathCore

dir /home/runner/work/Marmot/Marmot/modules