Basic arithmetic¶

Addition¶

template<typename _TensorLHS, typename _TensorRHS, typename std::enable_if_t<tmech::is_tensor_type<typename std::decay<_TensorLHS>::type>::value>* = nullptr, typename std::enable_if_t<tmech::is_tensor_type<typename std::decay<_TensorRHS>::type>::value>* = nullptr> constexpr auto operator+(_TensorLHS &&__tensor_lhs, _TensorRHS &&__tensor_rhs)¶

Addition of two tensors of same rank and dimension

$C_{ijkl...} = A_{ijkl...} + B_{ijkl...}$

.

tmech::tensor<double, 3, 2> a,b,c;
a.randn(); b.randn();
c = a + b;

Template Parameters

_TensorLHS: Left hand side tensor object
_TensorRHS: Right hand side tensor object

Parameters

__tensor_lhs: Left hand side tensor expression.
__tensor_rhs: Right hand side tensor expression.

Subtraction¶

template<typename _TensorLHS, typename _TensorRHS, typename std::enable_if_t<tmech::is_tensor_type<typename std::decay<_TensorLHS>::type>::value>* = nullptr, typename std::enable_if_t<tmech::is_tensor_type<typename std::decay<_TensorRHS>::type>::value>* = nullptr> constexpr auto operator-(_TensorLHS &&__tensor_lhs, _TensorRHS &&__tensor_rhs)¶

Subtraction of two tensors of same rank and dimension.

$C_{ijkl...} = A_{ijkl...} - B_{ijkl...}$

.

tmech::tensor<double, 3, 2> a,b,c;
a.randn(); b.randn();
c = a - b;

Template Parameters

_TensorLHS: Left hand side tensor object
_TensorRHS: Right hand side tensor object

Parameters

__tensor_lhs: Left hand side tensor expression.
__tensor_rhs: Right hand side tensor expression.

Single contraction¶

template<typename _TensorLHS, typename _TensorRHS, typename std::enable_if_t<tmech::is_tensor_type<typename std::decay<_TensorLHS>::type>::value>* = nullptr, typename std::enable_if_t<tmech::is_tensor_type<typename std::decay<_TensorRHS>::type>::value>* = nullptr> constexpr auto operator*(_TensorLHS &&__tensor_lhs, _TensorRHS &&__tensor_rhs)¶

Single contraction of two tensors. Assume that the left hand side tensor is of rank m and the right hand side tensor of rank n. The resulting tensor is of rank m+n-2. The most right index of the left hand side tensor is contracted with the most left index of the right hand tensor

$C_{i_{1},...,i_{m+n-2}} = A_{i_{1},...,i_{m-1},k} * B_{k,i_{2},...,i_{n}}$

.

tmech::tensor<double, 3, 2> a,b,c;
a.randn(); b.randn();
c = a * b;

Template Parameters

_TensorLHS: Left hand side tensor object
_TensorRHS: Right hand side tensor object

Parameters

__tensor_lhs: Left hand side tensor expression.
__tensor_rhs: Right hand side tensor expression.

Scalar update¶

template<typename _Scalar, typename _Tensor, typename, typename> constexpr auto operator*(_Scalar &&__scalar, _Tensor &&__tensor)¶

Scalar update of a tensor

$C_{ijkl...} = a*C_{ijkl...}$

.

tmech::tensor<double, 3, 2> a,b;
a.randn();
b = 2.0*a;

Template Parameters

_Tensor: Tensor object
_Scalar: Scalar type. Must be a std::fundamental

Parameters

__scalar: Scalar
__tensor: Tensor expression.

template<typename _Scalar, typename _Tensor, typename std::enable_if<std::is_fundamental<typename std::decay<_Scalar>::type>::value>::type* = nullptr, typename std::enable_if_t<tmech::is_tensor_type<typename std::decay<_Tensor>::type>::value>* = nullptr> constexpr auto operator/(_Tensor &&__tensor, _Scalar &&__scalar)¶

Tensor division by a scalar

$C_{ijkl...} = C_{ijkl...}/a$

.

tmech::tensor<double, 3, 2> a,b;
a.randn();
b = a/2.0;

Template Parameters

_Tensor: Tensor object
_Scalar: Scalar type. Must be a std::fundamental

Parameters

__scalar: Scalar
__tensor: Tensor expression.

Basis rearrangement¶

General basis rearrangement¶

template<typename _Sequence, typename _Tensor, typename, typename> constexpr auto tmech::basis_change(_Tensor &&__tensor)¶

General basis rearrangement. Controlled by template parameter _Sequence, which contains the new order of bases.

//Basis 1,2,3,4 is swaped to 3,4,1,2.
tmech::tensor<double, 3, 4> A, B;
A.randn();
B = tmech::basis_change<tmech::sequence<3,4,1,2>>(A);
B = tmech::basis_change<tmech::sequence<3,4,1,2>>(A+2*A);

Self assignment

//Basis 1,2,3,4 is swaped to 3,4,1,2.
tmech::tensor<double, 3, 4> A;
A.randn();
A = tmech::eval(tmech::basis_change<tmech::sequence<3,4,1,2>>(A));
A = tmech::eval(tmech::basis_change<tmech::sequence<3,4,1,2>>(A+2*A));

Template Parameters

_Sequence: A tmech::sequence<>, which contains the new order of bases
_Derived: Tensor object

Parameters

__tensor: Tensor expression from which the basis rearrangement is to be formed

template<typename _Sequence, typename _Tensor, typename, typename> constexpr auto tmech::basis_change(_Tensor &&__tensor)

General basis rearrangement. Controlled by template parameter _Sequence, which contains the new order of bases.

//Basis 1,2,3,4 is swaped to 3,4,1,2.
tmech::tensor<double, 3, 4> A, B;
A.randn();
B = tmech::basis_change<tmech::sequence<3,4,1,2>>(A);
B = tmech::basis_change<tmech::sequence<3,4,1,2>>(A+2*A);

Self assignment

//Basis 1,2,3,4 is swaped to 3,4,1,2.
tmech::tensor<double, 3, 4> A;
A.randn();
A = tmech::eval(tmech::basis_change<tmech::sequence<3,4,1,2>>(A));
A = tmech::eval(tmech::basis_change<tmech::sequence<3,4,1,2>>(A+2*A));

Template Parameters

_Sequence: A tmech::sequence<>, which contains the new order of bases
_Derived: Tensor object

Parameters

__tensor: Tensor expression from which the basis rearrangement is to be formed

Transposition¶

template<typename _Tensor, typename> constexpr auto tmech::trans(_Tensor &&__tensor)¶

Wrapper function for transposition for a second-order tensor and major transposition of a fourth-order tensor.

The transpose of a second-order tensor in index notation $(A_{ij})^T = A_{ji}$ .\ The transpose of a fourth-order tensor in index notation $(A_{ijkl})^T = A_{klij}$ .

tmech::tensor<double,3,2> A, B;
A.randn();
B = tmech::trans(A);

tmech::tensor<double,3,2> A;
A.randn();
A = tmech::eval(tmech::trans(A));

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression from which the transposition is to be formed.

Left basis transposition¶

template<typename _Tensor, typename> constexpr auto tmech::transl(_Tensor &&__tensor)¶

Transposition of the left base a fourth-order tensor $(A_{ijkl})^{l} = A_{jikl}$ .

tmech::tensor<double,3,4> A, B;
A.randn();
B = tmech::transl(A);

tmech::tensor<double,3,4> A;
A.randn();
A = tmech::eval(tmech::transl(A));

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression from which the left transpose part is to be formed.

Outer product¶

General outer product¶

template<typename _SequenceLHS, typename _SequenceRHS, typename _TensorLHS, typename _TensorRHS, typename, typename> constexpr auto tmech::outer_product(_TensorLHS &&__tensor_lhs, _TensorRHS &&__tensor_rhs)¶

General outer product. Controlled by template parameters _SequenceLHS and _SequenceRHS. Bases contained in sequence _SequenceLHS are used for ordered element acces in __lhs tensor expression. Bases contained in sequence _SequenceRHS are used for ordered element acces in __rhs tensor expression. Both togther peform outer product.

using SeqL = tmech::sequence<1,2>;
using SeqR = tmech::sequence<3,4>;
//Bases 1,2 of the new tensor C
//are given by a and bases 3,4 are
//given by b.
tmech::tensor<double, 3, 2> a, b;
tmech::tensor<double, 3, 4> C;
C = tmech::outer_product<SeqL, SeqR>(a,b);

Template Parameters

_SequenceLHS: Left hand side sequence, which contains numbers of bases to peform the outer product.
_SequenceRHS: Right hand side sequence, which contains numbers of bases to peform the outer product.
_TensorLHS: Left hand side tensor object
_TensorRHS: Right hand side tensor object

Parameters

__tensor_lhs: Left hand side tensor expression.
__tensor_rhs: Reft hand side tensor expression.

Dyadic product¶

template<typename _TensorLHS, typename _TensorRHS, typename, typename> constexpr auto tmech::otimes(_TensorLHS &&__tensor_lhs, _TensorRHS &&__tensor_rhs)¶

The dyadic product between a tensor object of order n and a tensor object of order m is a tensor of order m + n. The dyadic product between two first-order tensors is defined as

$A_{i_{1},...,i_{n},i_{n+1},...,i_{m}} = a_{i_{i},...,i_{n}} b_{i_{n+1},...,i_{m}}$

.

tmech::tensor<double,3,2> A, B;
tmech::tensor<double,3,4> C{tmech::otimes(A,B)};
tmech::tensor<double,3,6> D{tmech::otimes(C,B)};
tmech::tensor<double,3,8> E{tmech::otimes(C,C)};

Template Parameters

_TensorLHS: Left hand side tensor object
_TensorRHS: Right hand side tensor object

Parameters

__tensor_lhs: Left hand side tensor expression.
__tensor_rhs: Reft hand side tensor expression.

Lower dyadic product¶

template<typename _TensorLHS, typename _TensorRHS, typename, typename> constexpr auto tmech::otimesl(_TensorLHS &&__tensor_lhs, _TensorRHS &&__tensor_rhs)¶

The $\underline{\otimes}$ product between two second-order tensors is defined as

$\FourthT{C} = \SecondT{A} \underline{\otimes} \SecondT{B} = A_{ij} B_{kl} \SecondT{e}_i \otimes \SecondT{e}_l \otimes \SecondT{e}_j \otimes \SecondT{e}_k.$

.

tmech::tensor<double,3,2> A, B;
tmech::tensor<double,3,4> C{tmech::otimesl(A,B)};

Template Parameters

_TensorLHS: Left hand side tensor object
_TensorRHS: Right hand side tensor object

Parameters

__tensor_lhs: Left hand side tensor expression.
__tensor_rhs: Reft hand side tensor expression.

Upper dyadic product¶

template<typename _TensorLHS, typename _TensorRHS, typename, typename> constexpr auto tmech::otimesu(_TensorLHS &&__tensor_lhs, _TensorRHS &&__tensor_rhs)¶

The $\overline{\otimes}$ product between two second-order tensors is defined as

$\FourthT{C} = \SecondT{A} \overline{\otimes} \SecondT{B} = A_{ij} B_{kl} \SecondT{e}_i \otimes \SecondT{e}_k \otimes \SecondT{e}_j \otimes \SecondT{e}_l$

.

tmech::tensor<double,3,2> A, B;
tmech::tensor<double,3,4> C{tmech::otimesu(A,B)};

Template Parameters

_TensorLHS: Left hand side tensor object
_TensorRHS: Right hand side tensor object

Parameters

__tensor_lhs: Left hand side tensor expression.
__tensor_rhs: Reft hand side tensor expression.

Inner product¶

General inner product¶

template<typename _SequenceLHS, typename _SequenceRHS, typename _TensorLHS, typename _TensorRHS, typename, typename> constexpr auto tmech::inner_product(_TensorLHS &&__tensor_lhs, _TensorRHS &&__tensor_rhs)¶

General inner product. Controlled by template parameters _SequenceLHS and _SequenceRHS. Bases contained in sequence _SequenceLHS are contracted with bases contained in sequence _SequenceRHS.

using SeqL = tmech::sequence<3,4>;
using SeqR = tmech::sequence<1,2>;
//Double contraction of two 4th order tensors
tmech::tensor<double, 3, 4> A, B, C;
C = tmech::inner_product<SeqL, SeqR>(A,B);

//Double contraction of a 4th and a 2th order tensor
tmech::tensor<double, 3, 2> a, c;
c = tmech::inner_product<SeqL, SeqR>(C,a);

Template Parameters

_SequenceLHS: Left hand side sequence, which contains the numbers of bases used for contraction.
_SequenceRHS: Right hand side sequence, which contains the numbers of bases used for contraction.
_TensorLHS: Left hand side tensor object
_TensorRHS: Right hand side tensor object

Parameters

__tensor_lhs: Left hand side tensor expression.
__tensor_rhs: Reft hand side tensor expression.

Double contraction¶

template<typename _TensorLHS, typename _TensorRHS, typename, typename> constexpr auto tmech::dcontract(_TensorLHS &&__tensor_lhs, _TensorRHS &&__tensor_rhs)¶

A double contraction between two tensors objects contracts the two most right and left indices. The result of a double contraction between a tensor of order n and a tensor of order m is a tensor of order m + n - 4.

tmech::tensor<double,3,4> A;
tmech::tensor<double,3,2> B;
A.randn(); B.randn();
auto c{tmech::dcontract(A,B)};

Template Parameters

_TensorLHS: Left hand side tensor object
_TensorRHS: Right hand side tensor object

Parameters

__tensor_lhs: Left hand side tensor expression.
__tensor_rhs: Reft hand side tensor expression.

Fourth contraction¶

template<typename _TensorLHS, typename _TensorRHS, typename, typename> constexpr auto tmech::ddcontract(_TensorLHS &&__tensor_lhs, _TensorRHS &&__tensor_rhs)¶

A fourth contraction between two tensors objects contracts the fourth most right and left indices. The result of a double contraction between a tensor of order n and a tensor of order m is a tensor of order m + n - 8.

tmech::tensor<double,3,4> A;
tmech::tensor<double,3,2> B;
A.randn(); B.randn();
auto c{tmech::ddcontract(A,B)};

Template Parameters

_TensorLHS: Left hand side tensor object
_TensorRHS: Right hand side tensor object

Parameters

__tensor_lhs: Left hand side tensor expression.
__tensor_rhs: Reft hand side tensor expression.

Invers of a tensor¶

Second order tensors and fourth order tensors with minior-symmetry¶

template<typename ..._Sequences, typename _Tensor, typename> constexpr auto tmech::inv(_Tensor &&__tensor)¶

Inverse of a second- and a fourth-order tensor. The inverse of a second order tensor is defined by $\SecondT{A}^{-1}\SecondT{A} = \SecondT{I}$ , where $\SecondT{I}$ is the second order identity tensor.

tmech::tensor<double, 3, 2> A, B;
A.randn();
B = tmech::inv(A);

The inverse of a fourth order tensor depends on the minior-symmetry. Using the following differentiation rule $\FourthT{C}_{ijkl} = \frac{\partial \SecondT{a}_{ij} }{\partial \SecondT{b}_{kl}}$ yields a minior-symmetry in the first and in the second pair of indicies $i,j$ and $k,l$ , respectively. In this case the inverse is defined by $(\FourthT{A}^{-1})_{ijmn}\FourthT{A}_{mnkl} =\frac{1}{2}\left(\SecondT{I}_{ik}\SecondT{I}_{jl} + \SecondT{I}_{il}\SecondT{I}_{jk} \right)$ .

tmech::tensor<double, 3, 4> A, B;
A.randn();
//the sequences are indicating the pairs of minior-symmetry
B = tmech::inv<tmech::sequence<1,2>, tmech::sequence<3,4>>(A);
// the rule above is the most used one, therefore the following is sufficient.
// the correct indicies are set inside the function.
B = tmech::inv(A);

Using the following differentiation rule $\FourthT{C}_{iklj} = \frac{\partial \SecondT{a}_{ij} }{\partial \SecondT{b}_{kl}}$ , yields a minior symmetry in the outer and in the inner pair of indicies $i,j$ and $k,l$ , respectively.

tmech::tensor<double, 3, 4> A, B;
A.randn();
//the sequences are indicating the pairs of minior-symmetry
B = tmech::inv<tmech::sequence<1,4>, tmech::sequence<2,3>>(A);

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression from which the inverse is to be formed.

Fully anisotropic tensors¶

template<typename ..._Sequences, typename _Tensor, typename> constexpr auto tmech::invf(_Tensor &&__tensor)¶

Inverse of a second- and a full anisotropic fourth-order tensor. The inverse of a second order tensor is defined by $\SecondT{A}^{-1}\SecondT{A} = \SecondT{I}$ , where $\SecondT{I}$ is the second order identity tensor. In the case of a second oder tensor the ouput is the same as from the function inv.

tmech::tensor<double, 3, 2> A, B;
A.randn();
B = tmech::invf(A);
//or
B = tmech::inv(A);

This function assumes, that a fourth order tensor has no minior-symmetry and is therefore fully anisotropic. In this case the inverse is defined by $(\FourthT{A}^{-1})_{ijmn}\FourthT{A}_{mnkl} =\SecondT{I}_{ij}\SecondT{I}_{kl}$ .

tmech::tensor<double, 3, 4> A, B;
A.randn();
B = tmech::invf<tmech::sequence<1,2,3,4>>(A);
// the rule above is the most used one, therefore the following is sufficient.
// the correct indicies are set inside the function.
B = tmech::invf(A);

Using other kind of differentiation rules as $\FourthT{C}_{ijkl} = \frac{\partial \SecondT{a}_{ij} }{\partial \SecondT{b}_{kl}}$ must be indicated, due to the internal storage scheme. For example using $\FourthT{C}_{iklj} = \frac{\partial \SecondT{a}_{ij} }{\partial \SecondT{b}_{kl}}$ .

tmech::tensor<double, 3, 4> A, B;
A.randn();
B = tmech::invf<tmech::sequence<1,3,4,2>>(A);

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression from which the inverse is to be formed.

Volumetric and deviatoric parts¶

Volumetric part¶

template<typename _Tensor, typename> constexpr auto tmech::vol(_Tensor &&__tensor)¶

Volumetric part of a second-order tensor $\text{vol}(\SecondT{A}) = \frac{1}{d}\text{trace}(\SecondT{A})\SecondT{I}$ , where $d$ is the dimension.

tmech::tensor<double,3,2> A, B;
A.randn();
B = tmech::vol(A);

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression from which the volumetric part is to be formed.

Deviatoric part¶

template<typename _Tensor, typename> constexpr auto tmech::dev(_Tensor &&__tensor)¶

Deviatoric part of a second-order tensor $\SecondT{A} = \text{vol}(\SecondT{A}) + \text{dev}(\SecondT{A})$ .

tmech::tensor<double,3,2> A, B;
A.randn();
B = tmech::dev(A);

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression from which the deviatoric part is to be formed.

Symmetric and skew-symmetric parts¶

Symmetric part¶

template<typename _Tensor, typename> constexpr auto tmech::sym(_Tensor &&__tensor)¶

Symmetric part of a second-order tensor $\text{sym}\left(\SecondT{A}\right) = \frac{1}{2}\left(\SecondT{A} + \SecondT{A}^T \right)$ .

tmech::tensor<double,3,2> A, B;
A.randn();
B = tmech::sym(A);

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression from which the symmetric part is to be formed.

Skew-symmetric part¶

template<typename _Tensor, typename> constexpr auto tmech::skew(_Tensor &&__tensor)¶

Skew-symmetric part of a second-order tensor $\text{skew}\left(\SecondT{A}\right) = \frac{1}{2}\left(\SecondT{A} - \SecondT{A}^T \right)$ .

tmech::tensor<double, 3, 2> A, B;
A.randn();
B = tmech::skew(A);

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression from which the skew-symmetric part is to be formed.

Isotropic tensor functions¶

General function¶

template<typename _Function, typename _Tensor, typename> constexpr auto tmech::isotropic_tensor_function(_Tensor &&__tensor)¶

Isotropic tensor function of a symmetric second-order tensor.

$\SecondT{Y} = \sum_{i=1}^m y_i \SecondT{E}_i$

struct sqrt{
template<typename  T>
static constexpr inline auto apply(T const& value){return std::sqrt(value);}

template<typename  T>
static constexpr inline auto derivative(T const& value){return 1./(2.*std::sqrt(value));}
};

tmech::tensor<double, 3, 2> A, B;
A.randn();
B = tmech::isotropic_tensor_function<sqrt>(A);

Template Parameters

_Function:
_Tensor: Tensor object

Parameters

__tensor: Tensor expression

Square root¶

template<typename _Tensor, typename> constexpr auto tmech::sqrt(_Tensor &&__tensor)¶

Square root of a positive semi-definite symmetric second-order tensor $\text{sqrt}(\SecondT{A}) = \sqrt{\SecondT{A}}$ . The square is given by spectral decomposition

$\sqrt{\SecondT{A}} = \sum_{i}^m \sqrt{\lambda_i} \SecondT{E}_i,$

where $m$ is the number of non repeated eigenvalues $\lambda_i$ and $\SecondT{E}_i$ is the corresponding eigenbase.

tmech::tensor<double, 3, 2> A, B;
A = tmech::sym(tmech::randn<double,3,2>());
B = tmech::sqrt(A);

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression

Exponential map¶

template<typename _Tensor, typename> constexpr auto tmech::exp_sym(_Tensor &&__tensor)¶

Exponential map $\text{exp}(\SecondT{A})$ of a positive semi-definite symmetric second-order tensor. The exponential map is given by spectral decomposition

$\sqrt{\SecondT{A}} = \sum_{i}^m \exp{\lambda_i} \SecondT{E}_i,$

where $m$ is the number of non repeated eigenvalues $\lambda_i$ and $\SecondT{E}_i$ is the corresponding eigenbase.

tmech::tensor<double, 3, 2> A, B;
A = tmech::sym(tmech::randn<double,3,2>());
B = tmech::exp_sym(A);

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression

Logarithmic map¶

template<typename _Tensor, typename> constexpr auto tmech::log(_Tensor &&__tensor)¶

Logarithmic map $\text{log}(\SecondT{A})$ of a positive semi-definite symmetric second-order tensor. The Logarithmic map is given by spectral decomposition

$\sqrt{\SecondT{A}} = \sum_{i}^m \log{\lambda_i} \SecondT{E}_i,$

where $m$ is the number of non repeated eigenvalues $\lambda_i$ and $\SecondT{E}_i$ is the corresponding eigenbase.

tmech::tensor<double, 3, 2> A, B;
A = tmech::sym(tmech::randn<double,3,2>());
B = tmech::log(A);

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression

Positive-negative decomposition¶

General method¶

template<typename _Tensor, typename> constexpr auto tmech::positive_negative_decomposition(_Tensor &&__tensor)¶

Positive and negative decomposition $\SecondT{A} = \SecondT{A}^{+} + \SecondT{A}^{-}$ of a positive semi-definite symmetric second-order tensor is given by spectral decomposition

$\SecondT{A}^{+} = \sum_{i}^m \text{pos}(\lambda_i) \SecondT{E}_i,\quad \text{with} \quad \text{pos}(\lambda_i)= \begin{cases} \lambda_i ~\text{if}~ \lambda_i > 0\\ 0 ~\text{if}~ \lambda_i < 0 \end{cases}\\ {\SecondT{A}}^{-} = \sum_{i}^m \text{neg}(\lambda_i) \SecondT{E}_i,\quad \text{with} \quad \text{neg}(\lambda_i)= \begin{cases} \lambda_i ~\text{if}~ \lambda_i < 0\\ 0 ~\text{if}~ \lambda_i > 0 \end{cases}$

where $m$ is the number of non repeated eigenvalues $\lambda_i$ and $\SecondT{E}_i$ is the corresponding eigenbase.

tmech::tensor<double, 3, 2> A, B;
A = tmech::sym(tmech::randn<double,3,2>());
auto A_pos_neg = tmech::positive_negative_decomposition(A);
auto A_pos = A_pos_neg.positive();
auto A_neg = A_pos_neg.negative();

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression

Only positive part¶

template<typename _Tensor, typename> constexpr auto tmech::positive(_Tensor &&__tensor)¶

The positive part $\SecondT{A}^{+} = \SecondT{A} - \SecondT{A}^{-}$ of a positive semi-definite symmetric second-order tensor is given by spectral decomposition

$\SecondT{A}^{+} = \sum_{i}^m \text{pos}(\lambda_i) \SecondT{E}_i,\quad \text{with} \quad \text{pos}(\lambda_i)= \begin{cases} \lambda_i ~\text{if}~ \lambda_i > 0\\ 0 ~\text{if}~ \lambda_i < 0 \end{cases}$

where $m$ is the number of non repeated eigenvalues $\lambda_i$ and $\SecondT{E}_i$ is the corresponding eigenbase.

tmech::tensor<double, 3, 2> A, B;
A = tmech::sym(tmech::randn<double,3,2>());
B = tmech::positive(A);

If both parts $\SecondT{A}^{+}$ and $\SecondT{A}^{-}$ are needed use positive_negative_decomposition().

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression

Only negative part¶

template<typename _Tensor, typename> constexpr auto tmech::negative(_Tensor &&__tensor)¶

The negative part $\SecondT{A}^{-} = \SecondT{A} - \SecondT{A}^{+}$ of a positive semi-definite symmetric second-order tensor is given by spectral decomposition

${\SecondT{A}}^{-} = \sum_{i}^m \text{neg}(\lambda_i) \SecondT{E}_i,\quad \text{with} \quad \text{neg}(\lambda_i)= \begin{cases} \lambda_i ~\text{if}~ \lambda_i < 0\\ 0 ~\text{if}~ \lambda_i > 0 \end{cases}$

where $m$ is the number of non repeated eigenvalues $\lambda_i$ and $\SecondT{E}_i$ is the corresponding eigenbase.

tmech::tensor<double, 3, 2> A, B;
A = tmech::sym(tmech::randn<double,3,2>());
B = tmech::negative(A);

If both parts $\SecondT{A}^{+}$ and $\SecondT{A}^{-}$ are needed use positive_negative_decomposition().

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression

Eigen-decomposition¶

template<typename _Tensor, typename> constexpr auto tmech::eigen_decomposition(_Tensor &&__tensor)¶

Eigendecomposition of a positive semi-definite symmetric second-order tensor

$\SecondT{Y} = \sum_{i=1}^m \lambda_i \FirstT{e}_i \otimes \FirstT{e}_i = \sum_{i=1}^m \lambda_i \SecondT{E}_i,$

where $m$ is the number of non repeated eigenvalues, $\lambda_i$ are the corresponding eigenvalues, $\FirstT{e}_i$ eigenvectors and $\SecondT{E}_i$ eigenbasen.

Using eigenbases to get the inverse

tmech::tensor<double, 3, 2> A, A_inv;
A = tmech::sym(tmech::randn<double,3,2>());
auto A_eig = tmech::eigen_decomposition(A);
const auto[eigenvalues, eigenbasis]{A_eig.decompose_eigenbasis()};
const auto idx{A_eig.non_repeated_eigenvalues_index()};
for(std::size_t i{0}; i<A_eig.number_non_repeated_eigenvalues(); ++i){
A_inv += (1.0/eigenvalues[idx[i]])*eigenbasis[idx[i]];
}
std::cout<<std::boolalpha<<tmech::almost_equal(tmech::inv(A), A_inv, 5e-7)<<std::endl;

Using eigenvectors to get the inverse

tmech::tensor<double, 3, 2> A, A_inv;
A = tmech::sym(tmech::randn<double,3,2>());
auto A_eig = tmech::eigen_decomposition(A);
const auto[eigenvalues, eigenvectors]{A_eig.decompose()};
const auto idx{A_eig.non_repeated_eigenvalues_index()};
for(std::size_t i{0}; i<A_eig.number_non_repeated_eigenvalues(); ++i){
A_inv += (1.0/eigenvalues[idx[i]])*tmech::otimes(eigenvectors[idx[i]],eigenvectors[idx[i]]);
}
std::cout<<std::boolalpha<<tmech::almost_equal(tmech::inv(A), A_inv, 5e-7)<<std::endl;

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression, from which the eigen decomposition is to be computed.

Sign decomposition¶

template<typename _Tensor, typename> constexpr auto tmech::sign(_Tensor &&__tensor, typename std::decay<_Tensor>::type::value_type __eps = 5e-7, std::size_t __max_iter = 10)¶

Sign tensor decomposition of a second-order tensor

$\SecondT{A} = \SecondT{S}\SecondT{N}, \quad \SecondT{S} = \text{sign}(\SecondT{A}), \quad \SecondT{N} = \sqrt{\SecondT{A}\SecondT{A}}$

The decomposition is based on a Newton iteration as describet here.

$\SecondT{X}_{k+1} =\frac{1}{2}\left(\SecondT{X}_{k} + \SecondT{X}_{k}^{-1}\right), \quad \text{with} \quad \SecondT{X}_{0} = \SecondT{A}$

which converges quadratically to $\text{sign}(\SecondT{A})$ for any $\SecondT{A}$ having no pure imaginary eigenvalues.

tmech::tensor<double, 3, 2> A, B;
A.randn();
B = tmech::sign(A);

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression from which the inverse is to be formed.
__eps: Tolerance for Newton iteration
__max_iter: Maximum of Newton iterations

Polar decomposition¶

template<typename _Tensor, typename> constexpr auto tmech::polar_decomposition(_Tensor &&__tensor, bool const __newton_method, typename std::decay<_Tensor>::type::value_type const __tol, std::size_t const __max_steps)¶

Polar decomposition of a positive semi-definite symmetric second-order tensor

$\SecondT{F}=\SecondT{R}\SecondT{U} = \SecondT{V}\SecondT{R}$

where $\SecondT{R}$ is an orthogonal tensor also knwon as the rotation tensor, $\SecondT{U}$ and $\SecondT{V}$ are symmetric tensors called the right and the left stretch tensor, respectively. This function provides two different method to determine $\SecondT{U}$ , $\SecondT{V}$ and $\SecondT{R}$ . The frist method uses spectral decomposition

$\SecondT{U} = \sqrt{\SecondT{F}^T\SecondT{F}}, \quad \SecondT{R} = \SecondT{F}\SecondT{U}^{-1}, \quad \SecondT{V} = \SecondT{R}\SecondT{U}\SecondT{R}^T,$

.

//use the spectral decomposition
tmech::tensor<double, 3, 2> F;
F.randn();
auto F_polar = tmech::polar_decomposition(F);
F = F_polar.R()*F_polar.U();
F = F_polar.V()*F_polar.R();

The second one is based on a Newton iteration

$\SecondT{R}_{k+1} =\frac{1}{2} \left( \SecondT{R}_k + \SecondT{R}_k^{-T}\right), \quad \text{with}\quad \SecondT{R}_0 = \SecondT{F}$

//use the newton iteration
tmech::tensor<double, 3, 2> F;
F.randn();
auto F_polar = tmech::polar_decomposition(F, // tensor to decompose
                                          true, //use newton iteration
                                          5e-7, // tolerance
                                          5 // number of max steps
                                          );
F = F_polar.R()*F_polar.U();
F = F_polar.V()*F_polar.R();

The following derivatives are also important

$\frac{\partial \SecondT{U}}{\partial \SecondT{F}}, \quad \frac{\partial \SecondT{V}}{\partial \SecondT{F}}, \quad \frac{\partial \SecondT{R}}{\partial \SecondT{F}}$

explicit results are given here

//use the spectral decomposition
tmech::tensor<double, 3, 2> F;
F.randn();
auto F_polar = tmech::polar_decomposition(F);
auto dR = F_polar.R().derivative();
auto dU = F_polar.U().derivative();
auto dV = F_polar.V().derivative();

Template Parameters

_Tensor: Tensor object

Parameters

__tensor: Tensor expression, from which the polar decomposition is to be computed.
__newton_method: True if Newton iteration should be used
__tol: Tolerance for Newton iteration
__max_steps: Maximum number of Newton iterations

Exponential map¶

template<typename _Tensor, typename> constexpr auto tmech::exp(_Tensor &&__tensor)¶

Exponential map of a nonsymmetric second order tensor

$\text{exp}\SecondT{A} = \sum_{n=0}^\infty \frac{1}{n!} \SecondT{A}^n$

.

tmech::tensor<double, 3, 2> A{9.064107e-01, -4.649874e-01,  4.431378e-01,
                              1.557860e+00,  2.493285e-01, -3.458549e-01,
                              1.747649e+00, -3.824761e-01, -8.110930e-01};

tmech::tensor<double, 3, 2>Aexp{2.333496383103521e+00,  -9.424245337660712e-01,   6.174107765089896e-01,
                                2.458154759604529e+00,   7.615487934000856e-01,   1.403592316302926e-01,
                                1.701192794766057e+00,  -7.849258234991792e-01,   8.233614408596337e-01};

 std::cout<<std::boolalpha
          <<tmech::almost_equal(Aexp, tmech::exp(A), 5e-7)<<std::endl;

 //check the derivative
 const auto func_exp = [](auto const& tensor){return tmech::exp(tensor);};
 std::cout<<std::boolalpha
          <<tmech::almost_equal(tmech::num_diff_central<tmech::sequence<1,2,3,4>>(func_exp,A), tmech::exp(A).derivative(), 5e-7)<<std::endl;