contimech.hpp

Overview 

General continuum mechanics functions for tensor operations, transformations, and common mechanical quantities (deviatoric parts, traces, invariants, etc.).

API Reference 

arma::mat dev(const arma::mat &m)

Returns the deviatoric part of a 3x3 matrix.

The deviatoric part of a tensor \( \mathbf{m} \) is defined as:

\[ \mathbf{m}' = \mathbf{m} - \frac{1}{3} \textrm{tr}(\mathbf{m}) \mathbf{I} \]

where \( \textrm{tr}(\mathbf{m}) = m_{11} + m_{22} + m_{33} \) is the trace and \( \mathbf{I} \) is the identity tensor. Example:

mat m = randu(3,3);
mat m_dev = dev(m);

Parameters:: m – The input 3x3 matrix (arma::mat)
Returns:: The 3x3 deviatoric part of the matrix input (arma::mat)

arma::mat sph(const arma::mat &m)

Returns the spherical part of a 3x3 matrix.

The spherical (or hydrostatic) part of a tensor \( \mathbf{m} \) is defined as:

\[ \mathbf{m}^{sph} = \frac{1}{3} \textrm{tr}(\mathbf{m}) \mathbf{I} \]

where \( \textrm{tr}(\mathbf{m}) = m_{11} + m_{22} + m_{33} \) is the trace and \( \mathbf{I} \) is the identity tensor. Note that \( \mathbf{m} = \mathbf{m}' + \mathbf{m}^{sph} \). Example:

mat m = randu(3,3);
mat m_sph = sph(m);

Parameters:: m – The input 3x3 matrix (arma::mat)
Returns:: The 3x3 spherical part of the matrix input (arma::mat)

double tr(const arma::vec &v)

Returns the trace of a tensor v expressed in Voigt notation.

The trace of a second-order tensor \( \mathbf{m} \) written in Voigt notation as \( \mathbf{v} = (v_1, v_2, v_3, v_4, v_5, v_6)^T \) is:

\[ \textrm{tr}(\mathbf{m}) = v_1 + v_2 + v_3 = m_{11} + m_{22} + m_{33} \]

Example:

vec v = randu(6);
double trace = tr(v);

Parameters:: v – The input tensor (arma::vec)
Returns:: The trace of the input tensor (double)

arma::vec dev(const arma::vec &v)

Returns the deviatoric part of a tensor v expressed in Voigt notation.

The deviatoric part of a tensor \( \mathbf{m} \) expressed in Voigt notation is:

\[ \mathbf{v}' = \mathbf{v} - \frac{1}{3} \textrm{tr}(\mathbf{v}) \mathbf{I}_{th} \]

where \( \mathbf{I}_{th} = (1, 1, 1, 0, 0, 0)^T \) is the expansion vector. This gives \( v'_i = v_i - \frac{1}{3}(v_1 + v_2 + v_3) \) for \( i = 1,2,3 \) and \( v'_i = v_i \) for \( i = 4,5,6 \). Example:

vec v = randu(6);
vec v_dev = dev(v);

Parameters:: v – The input tensor (arma::vec)
Returns:: The 6 vector, deviatoric part (arma::vec)

double Mises_stress(const arma::vec &v)

Provides the Von Mises stress \( \sigma^{Mises} \) of a second order stress tensor written as a vector v in Voigt notation.

The Von Mises equivalent stress is defined as:

\[ \sigma^{Mises} = \sqrt{ \frac{3}{2} \mathbf{\sigma}' : \mathbf{\sigma}' } = \sqrt{3 \, J_2} \]

where \( \mathbf{\sigma}' = \textrm{dev}(\mathbf{\sigma}) \) is the deviatoric stress tensor and \( J_2 \) is the second deviatoric invariant. Example:

vec v = randu(6);
double Mises_sig = Mises_stress(v);

Parameters:: v – The input stress tensor (arma::vec)
Returns:: The Mises equivalent (double)

arma::vec eta_stress(const arma::vec &v)

Provides the stress flow \( \boldsymbol{\eta}_{stress} \) of a second order stress tensor written as a vector v in Voigt notation.

The stress flow direction is defined as:

\[ \boldsymbol{\eta}_{stress} = \frac{\frac{3}{2} \mathbf{\sigma}'}{\sigma^{Mises}} = \frac{\frac{3}{2} \textrm{dev}(\mathbf{\sigma})}{\sigma^{Mises}} \]

This corresponds to the normal to the Von Mises yield surface and is used in associated flow rules. Example:

vec v = randu(6);
vec v_flow = eta_stress(v);

Parameters:: v – The input stress tensor (arma::vec)
Returns:: The 6 vector, flow (arma::vec)

arma::mat deta_stress(const arma::vec &v)

Provides the derivative of the J2 stress flow w.r.t. stress, \( \partial \boldsymbol{\eta}_{stress} / \partial \boldsymbol{\sigma} \) (6x6 Voigt).

With \( \boldsymbol{\eta}_{stress} = \mathbf{P}_{Mises}\,\boldsymbol{\sigma} / \sigma^{Mises} \), the derivative is the normalized-normal Hessian

\[ \frac{\partial \boldsymbol{\eta}_{stress}}{\partial \boldsymbol{\sigma}} = \frac{1}{\sigma^{Mises}} \left( \mathbf{P}_{Mises} - \boldsymbol{\eta}_{stress}\otimes\boldsymbol{\eta}_{stress} \right), \]

where \( \mathbf{P}_{Mises} \) is the Voigt metric with shear diagonal 3 (the engineering shear factor is already embedded, consistently with eta_stress() and Mises_stress()). The result is symmetric and satisfies \( (\partial \boldsymbol{\eta}_{stress}/\partial \boldsymbol{\sigma})\,\boldsymbol{\sigma} = \mathbf{0} \) (degree-0 homogeneity). Returns \( \mathbf{0}_{6\times6} \) when \( \sigma^{Mises} < \iota \). It supplies \( \partial \boldsymbol{\Lambda}_\varepsilon / \partial \boldsymbol{\sigma} \) for the Simo-Hughes algorithmic tangent of J2-associated UMATs.

Note

Voigt type: compliance-like (maps stress -> strain-flow, \( \partial \boldsymbol{\varepsilon} / \partial \boldsymbol{\sigma} \) structure). Rows are strain-like (engineering, doubled shear — as eta_stress), columns are stress-like (single shear); it transforms as \( \mathbf{A}' = \mathbf{Q}_\varepsilon\,\mathbf{A}\,\mathbf{Q}_\sigma^{-1} \) (like a compliance \( \mathbf{M} \)), NOT as a stiffness. This matches its use as \( \mathbf{L}\,(\Delta s\,\partial\boldsymbol{\Lambda}/\partial\boldsymbol{\sigma}) \) inside assemble_algorithmic_tangent(), where the strain-like output meets the stiffness \( \mathbf{L} \).

Parameters:: v – The input stress tensor (arma::vec, 6)
Returns:: The 6x6 Hessian of the Mises equivalent stress (arma::mat)

arma::vec eta_norm_stress(const arma::vec &v)

Provides the strain flow (direction) from a stress tensor (Euclidian norm), according to the Voigt convention for strains.

The normalized flow direction based on the Euclidian norm is:

\[ \boldsymbol{\eta}_{norm} = \frac{\mathbf{\sigma}}{\|\mathbf{\sigma}\|} = \frac{\mathbf{\sigma}}{\sqrt{\mathbf{\sigma}:\mathbf{\sigma}}} \]

Note that the Euclidian norm is expressed as \( \|\mathbf{m}\| = \sqrt{\mathbf{m}:\mathbf{m}} \), considering that v is the vector representation of a stress matrix m.

vec v = randu(6);
vec v_flow = eta_norm_stress(v);

Parameters:: v – The input stress tensor (arma::vec)
Returns:: The 6 vector, flow (arma::vec)

arma::vec eta_norm_strain(const arma::vec &v)

Provides the strain flow (direction) from a strain tensor (Euclidian norm), according to the Voigt convention for strains.

The normalized flow direction based on the Euclidian norm is:

\[ \boldsymbol{\eta}_{norm} = \frac{\boldsymbol{\varepsilon}}{\|\boldsymbol{\varepsilon}\|} = \frac{\boldsymbol{\varepsilon}}{\sqrt{\boldsymbol{\varepsilon}:\boldsymbol{\varepsilon}}} \]

Note that the Euclidian norm is expressed as \( \|\mathbf{m}\| = \sqrt{\mathbf{m}:\mathbf{m}} \), considering that v is the vector representation of a strain matrix m.

vec v = randu(6);
vec v_flow = eta_norm_strain(v);

Parameters:: v – The input strain tensor (arma::vec)
Returns:: The 6 vector, flow (arma::vec)

double norm_stress(const arma::vec &v)

Provides the Euclidian norm of a stress tensor.

The Euclidian (Frobenius) norm of a stress tensor is:

\[ \|\mathbf{\sigma}\| = \sqrt{\mathbf{\sigma}:\mathbf{\sigma}} = \sqrt{\sigma_{ij}\,\sigma_{ij}} \]

In Voigt notation: \( \|\mathbf{\sigma}\| = \sqrt{v_1^2 + v_2^2 + v_3^2 + 2(v_4^2 + v_5^2 + v_6^2)} \)

vec v = randu(6);
double norm=norm_stress(v);

Parameters:: v – The input stress tensor (arma::vec)
Returns:: The norm (double)

double norm_strain(const arma::vec &v)

Provides the Euclidian norm of a strain tensor.

The Euclidian (Frobenius) norm of a strain tensor is:

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

In Voigt notation: \( \|\boldsymbol{\varepsilon}\| = \sqrt{v_1^2 + v_2^2 + v_3^2 + \frac{1}{2}(v_4^2 + v_5^2 + v_6^2)} \)

vec v = randu(6);
double norm=norm_strain(v);

Parameters:: v – The input strain tensor (arma::vec)
Returns:: The norm (double)

double Mises_strain(const arma::vec &v)

Provides the Von Mises strain \( \varepsilon^{Mises} \) of a second order strain tensor.

The Von Mises equivalent strain is defined as:

\[ \varepsilon^{Mises} = \sqrt{ \frac{2}{3} \boldsymbol{\varepsilon}' : \boldsymbol{\varepsilon}' } \]

where \( \boldsymbol{\varepsilon}' = \textrm{dev}(\boldsymbol{\varepsilon}) \) is the deviatoric strain tensor. This definition is conjugate to the Von Mises stress such that \( \sigma^{Mises} \cdot \varepsilon^{Mises} = \mathbf{\sigma}' : \boldsymbol{\varepsilon}' \).

vec v = randu(6);
double Mises_eps = Mises_strain(v);

Parameters:: v – The input strain tensor (arma::vec)
Returns:: The Mises norm (double)

arma::vec eta_strain(const arma::vec &v)

Provides the strain flow \( \boldsymbol{\eta}_{strain} \) of a second order strain tensor written as a vector v in Voigt notation.

The strain flow direction is defined as:

\[ \boldsymbol{\eta}_{strain} = \frac{\frac{2}{3} \boldsymbol{\varepsilon}'}{\varepsilon^{Mises}} = \frac{\frac{2}{3} \textrm{dev}(\boldsymbol{\varepsilon})}{\varepsilon^{Mises}} \]

This is conjugate to the stress flow \( \boldsymbol{\eta}_{stress} \) used in plasticity. Example:

vec v = randu(6);
vec eps_f = eta_strain(v);

Parameters:: v – The input strain tensor (arma::vec)
Returns:: The 6 vector, flow (arma::vec)

double J2_stress(const arma::vec &v)

Provides the second invariant of the deviatoric part of a second order stress tensor, written as a vector v in Voigt notation.

The second deviatoric stress invariant \( J_2 \) is defined as:

\[ J_2 = \frac{1}{2} \mathbf{\sigma}' : \mathbf{\sigma}' = \frac{1}{2} \sigma'_{ij} \, \sigma'_{ij} \]

It is related to the Von Mises stress by \( \sigma^{Mises} = \sqrt{3 \, J_2} \).

vec v = randu(6);
double J2 = J2_stress(v);

Parameters:: v – The input stress tensor (arma::vec)
Returns:: The second invariant (double)

double J2_strain(const arma::vec &v)

Provides the second invariant of the deviatoric part of a second order strain tensor, written as a vector v in Voigt notation.

The second deviatoric strain invariant \( J_2 \) is defined as:

\[ J_2 = \frac{1}{2} \boldsymbol{\varepsilon}' : \boldsymbol{\varepsilon}' = \frac{1}{2} \varepsilon'_{ij} \, \varepsilon'_{ij} \]

where \( \boldsymbol{\varepsilon}' = \textrm{dev}(\boldsymbol{\varepsilon}) \) is the deviatoric strain tensor.

vec v = randu(6);
double J2 = J2_strain(v);

Parameters:: v – The input strain tensor (arma::vec)
Returns:: The second invariant (double)

double J3_stress(const arma::vec &v)

Provides the third invariant of the deviatoric part of a second order stress tensor, written as a vector v in Voigt notation.

The third deviatoric stress invariant \( J_3 \) is defined as:

\[ J_3 = \det(\mathbf{\sigma}') = \frac{1}{3} \sigma'_{ij} \, \sigma'_{jk} \, \sigma'_{ki} \]

This invariant characterizes the Lode angle and is used in pressure-dependent yield criteria.

vec v = randu(6);
double J3 = J3_stress(v);

Parameters:: v – The input stress tensor (arma::vec)
Returns:: The third invariant (double)

double J3_strain(const arma::vec &v)

Provides the third invariant of the deviatoric part of a second order strain tensor, written as a vector v in Voigt notation.

The third deviatoric strain invariant \( J_3 \) is defined as:

\[ J_3 = \det(\boldsymbol{\varepsilon}') = \frac{1}{3} \varepsilon'_{ij} \, \varepsilon'_{jk} \, \varepsilon'_{ki} \]

where \( \boldsymbol{\varepsilon}' = \textrm{dev}(\boldsymbol{\varepsilon}) \) is the deviatoric strain tensor.

vec v = randu(6);
double J3 = J3_strain(v);

Parameters:: v – The input strain tensor (arma::vec)
Returns:: The third invariant (double)

double Macaulay_p(const double &d)

Provides the results of the MacCaulay brackets operator \(\langle d \rangle^+\).

The positive Macaulay bracket (ramp function) is defined as:

\[\begin{split} \langle d \rangle^+ = \begin{cases} d & \text{if } d > 0 \\ 0 & \text{if } d \leq 0 \end{cases} = \frac{d + |d|}{2} \end{split}\]

constexpr int MIN = -10;
constexpr int MAX = 100;
double d = MIN + (double)(rand()) / ((double)(RAND_MAX/(MAX - MIN)));
double d_MC = Macaulay_p(d);

Parameters:: d – The input value (double)
Returns:: The value of \(\langle d \rangle^+\) (double)

double Macaulay_n(const double &d)

Provides the results of the MacCaulay brackets operator \(\langle d \rangle^-\).

The negative Macaulay bracket is defined as:

\[\begin{split} \langle d \rangle^- = \begin{cases} 0 & \text{if } d \geq 0 \\ d & \text{if } d < 0 \end{cases} = \frac{d - |d|}{2} \end{split}\]

Note that \( d = \langle d \rangle^+ + \langle d \rangle^- \).

constexpr int MIN = -10;
constexpr int MAX = 100;
double d = MIN + (double)(rand()) / ((double)(RAND_MAX/(MAX - MIN)));
double d_MC = Macaulay_n(d);

Parameters:: d – The input value (double)
Returns:: The value of \(\langle d \rangle^-\) (double)

double sign(const double &d)

Provides the results of the sign operator.

The sign (signum) function is defined as:

\[\begin{split} \textrm{sign}(d) = \begin{cases} +1 & \text{if } d > \iota \\ 0 & \text{if } |d| \leq \iota \\ -1 & \text{if } d < -\iota \end{cases} \end{split}\]

where \( \iota \) is a small tolerance (set in the simcoon parameters, by default \( 10^{-12} \)).

constexpr int MIN = -10;
constexpr int MAX = 100;
double d = MIN + (double)(rand()) / ((double)(RAND_MAX/(MAX - MIN)));
double d_sign = sign(d);

Parameters:: d – The input value (double)
Returns:: The sign of the input value (double)

arma::vec normal_ellipsoid(const double &u, const double &v, const double &a1, const double &a2, const double &a3)

Provides the normalized vector normal to an ellipsoid with semi-principal axes of length a1, a2, a3. The direction of the normalized vector is set by angles u and v.

Provides the normalized vector to an ellipsoid with semi-principal axes of length a1, a2, a3. The direction of the normalized vector is set by angles u and v. These 2 angles correspond to the rotations in the plan defined by the center of the ellipsoid, a1 and a2 directions for u, a1 and a3 ones for v. u = 0 corresponds to a1 direction and v = 0 correspond to a3 one. So the normal vector is set at the parametrized position :

\[\begin{split}\begin{align} x & = a_{1} \cos(u) \sin(v) \\ y & = a_{2} \sin(u) \sin(v) \\ z & = a_{3} \cos(v) \end{align}\end{split}\]

const double Pi = 3.14159265358979323846;

double u = (double)rand()/(double)(RAND_MAX) % (2*Pi) - (2*Pi);
double v = (double)rand()/(double)(RAND_MAX) % Pi - Pi;
double a1 = (double)rand();
double a2 = (double)rand();
double a3 = (double)rand();
vec normal = normal_ellipsoid(u, v, a1, a2, a3);

Parameters:

u – The angle in the plane defined by the center of the ellipsoid, a1 and a2 directions (double)
v – The angle in the plane defined by the center of the ellipsoid, a1 and a3 directions (double)
a1 – The length of the semi-principal axis a1 (double)
a2 – The length of the semi-principal axis a2 (double)
a3 – The length of the semi-principal axis a3 (double)

Returns:

The normal vector (arma::vec)

double curvature_ellipsoid(const double &u, const double &v, const double &a1, const double &a2, const double &a3)

Provides the curvature of an ellipsoid with semi-principal axes of length a1, a2, a3 at the angle u,v.

Provides the normalized curvature of an ellipsoid with semi-principal axes of length a1, a2, a3. The position of the evaluated curvature is set by angles u and v. These 2 angles correspond to the rotations in the plan defined by the center of the ellipsoid, a1 and a2 directions for u, a1 and a3 ones for v. u = 0 corresponds to a1 direction and v = 0 correspond to a3 one. So the curvature is set at the parametrized position :

\[\begin{split}\begin{align} x & = a_{1} \cos(u) \sin(v) \\ y & = a_{2} \sin(u) \sin(v) \\ z & = a_{3} \cos(v) \end{align}\end{split}\]

and the curvature is

\[ \Xi = \frac{a_1^2\,a_2^2\,a_3^2}{\left( a_1^2\,a_2^2\,\cos^2 v + a_3^2\,\sin^2 v \left( a_2^2\,\cos^2 u + a_1^2\,\sin^2 u \right) \right)^2} \]

const double Pi = 3.14159265358979323846;

double u = (double)rand()/(double)(RAND_MAX) % (2*Pi) - (2*Pi);
double v = (double)rand()/(double)(RAND_MAX) % Pi - Pi;
double a1 = (double)rand();
double a2 = (double)rand();
double a3 = (double)rand();
double curvature = curvature_ellipsoid(u, v, a1, a2, a3);

Parameters:

u – The angle in the plane defined by the center of the ellipsoid, a1 and a2 directions (double)
v – The angle in the plane defined by the center of the ellipsoid, a1 and a3 directions (double)
a1 – The length of the semi-principal axis a1 (double)
a2 – The length of the semi-principal axis a2 (double)
a3 – The length of the semi-principal axis a3 (double)

Returns:

The curvature (double)

arma::vec sigma_int(const arma::vec &sigma_in, const double &u, const double &v, const double &a1, const double &a2, const double &a3)

Provides the normal and tangent components of the traction vector in the normal direction n to an ellipsoid with axes a1, a2, a3 from an input inside stress \( \sigma \). The direction of the normalized vector is set by angles u and v.

Provides the normal and tangent components of a stress vector σin in accordance with the normal direction n to an ellipsoid with axes a1, a2, a3. The normal vector is set at the parametrized position :

\[\begin{split}\begin{align} x & = a_{1} \cos(u) \sin(v) \\ y & = a_{2} \sin(u) \sin(v) \\ z & = a_{3} \cos(v) \end{align}\end{split}\]

const double Pi = 3.14159265358979323846;

vec sigma_in = randu(6);
double u = (double)rand()/(double)(RAND_MAX) % Pi - Pi/2;
double v = (double)rand()/(double)(RAND_MAX) % (2*Pi) - Pi;
double a1 = (double)rand();
double a2 = (double)rand();
double a3 = (double)rand();
vec sigma_i = sigma_int(sigma_in, u, v, a1, a2, a3);

Parameters:

sigma_in – The input stress tensor (arma::vec)
u – The angle in the plane defined by the center of the ellipsoid, a1 and a2 directions (double)
v – The angle in the plane defined by the center of the ellipsoid, a1 and a3 directions (double)
a1 – The length of the semi-principal axis a1 (double)
a2 – The length of the semi-principal axis a2 (double)
a3 – The length of the semi-principal axis a3 (double)

Returns:

The normal and tangent components of the traction vector (arma::vec)

arma::mat p_ikjl(const arma::vec &a)

Provides the Hill interfacial operator according to a normal a (see papers of Siredey and Entemeyer phD dissertation).

The Hill interfacial operator \( \mathbf{P} \) is defined in index notation as:

\[ P_{ikjl} = \frac{1}{2} \left( a_i \, \delta_{kj} \, a_l + a_j \, \delta_{ik} \, a_l + a_i \, \delta_{lj} \, a_k + a_j \, \delta_{il} \, a_k \right) \]

where \( \mathbf{a} \) is the normal vector and \( \delta_{ij} \) is the Kronecker delta.

vec a = randu(3);
mat H = p_ikjl(a);

Parameters:: a – The input normal vector (arma::vec)
Returns:: Hill interfacial operator 6x6 matrix (arma::mat)

arma::mat auto_sym_dyadic(const arma::mat &a)

Provides the dyadic product of a symmetric tensor with itself (auto dyadic product).

This function returns the dyadic (tensor) product of a symmetric tensor with itself:

\[ \mathbf{C} = \mathbf{a} \otimes \mathbf{a}, \quad C_{ijkl} = a_{ij} \, a_{kl} \]

The function returns a 6x6 matrix that corresponds to a 4th order tensor in Voigt notation (such as stiffness matrices).

Parameters:

a – The input symmetric tensor (arma::mat)

Returns:

The 6x6 matrix that represent the dyadic product (arma::mat)

mat a = randu(3,3);
mat sym_a = 0.5*(a + a.t());
mat c = auto_sym_dyadic(sym_a);

arma::mat sym_dyadic(const arma::mat &a, const arma::mat &b)

Provides the dyadic product of two symmetric tensors.

This function returns the dyadic (tensor) product of two symmetric tensors:

\[ \mathbf{C} = \mathbf{a} \otimes \mathbf{b}, \quad C_{ijkl} = a_{ij} \, b_{kl} \]

The function returns a 6x6 matrix that corresponds to a 4th order tensor in Voigt notation (such as stiffness matrices).

Parameters:

a – The first input symmetric tensor (arma::mat)
b – The second input symmetric tensor (arma::mat)

Returns:

The 6x6 matrix that represent the dyadic product (arma::mat)

mat a = randu(3,3);
mat sym_a = 0.5*(a + a.t());
mat b = randu(3,3);
mat sym_b = 0.5*(b + b.t());
mat c = sym_dyadic(sym_a,sym_b);

arma::mat auto_dyadic(const arma::mat &a)

Provides the dyadic product of a tensor with itself (auto dyadic product).

This function returns the dyadic (tensor) product of a tensor with itself:

\[ \mathbf{C} = \mathbf{a} \otimes \mathbf{a}, \quad C_{ijkl} = a_{ij} \, a_{kl} \]

The function returns a 6x6 matrix that corresponds to a 4th order tensor in Voigt notation. Unlike auto_sym_dyadic, this function does not assume the input tensor is symmetric.

Parameters:

a – The input tensor (arma::mat)

Returns:

The 6x6 matrix that represent the dyadic product (arma::mat)

mat a = randu(3,3);
mat c = auto_dyadic(a);

arma::mat dyadic_4vectors_sym(const arma::vec &n_a, const arma::vec &n_b, const std::string &conv)

Provides the dyadic product of four vectors to provide a symmetric 4th order tensor (minor symmetry).

This function returns a 4th order tensor from two vectors with minor symmetry. For convention “aabb”:

\[ C_{ijkl} = (n_a)_i \, (n_a)_j \, (n_b)_k \, (n_b)_l \]

For convention “abab”:

\[ C_{ijkl} = \frac{1}{2} \left( (n_a)_i \, (n_b)_j \, (n_a)_k \, (n_b)_l + (n_a)_i \, (n_b)_j \, (n_b)_k \, (n_a)_l \right) \]

The function returns a 6x6 matrix corresponding to a 4th order tensor in Voigt notation.

Parameters:

n_a – The first input vector (arma::vec)
n_b – The second input vector (arma::vec)
conv – The convention for the dyadic product (“aabb” or “abab”) (std::string)

Returns:

The 6x6 matrix that represent the dyadic product (arma::mat)

vec n_a = randu(3);
vec n_b = randu(3);
mat c = dyadic_4vectors_sym(n_a, n_b, "aabb");

arma::mat dyadic(const arma::mat &a, const arma::mat &b)

Provides the dyadic product of two tensors.

This function returns the dyadic (tensor) product of two tensors:

\[ \mathbf{C} = \mathbf{a} \otimes \mathbf{b}, \quad C_{ijkl} = a_{ij} \, b_{kl} \]

The function returns a 6x6 matrix corresponding to a 4th order tensor in Voigt notation (such as stiffness matrices).

Parameters:

a – The first input tensor (arma::mat)
b – The second input tensor (arma::mat)

Returns:

The 6x6 matrix that represent the dyadic product (arma::mat)

mat a = randu(3,3);
mat b = randu(3,3);
mat c = dyadic(a,b);

arma::mat auto_sym_dyadic_operator(const arma::mat &a)

Provides the symmetric 4th-order dyadic product of a symmetric tensor with itself.

This function returns the symmetric dyadic product \( \mathbf{C} = \mathbf{a} \odot \mathbf{a} \), defined as:

\[ \left(\mathbf{a} \odot \mathbf{a} \right)_{ijkl} = \frac{1}{2} \left( a_{ik} \, a_{jl} + a_{il} \, a_{jk} \right) \]

The function returns a 6x6 matrix corresponding to a 4th order tensor in Voigt notation (possessing both major and minor symmetries).

Parameters:

a – The input symmetric tensor (arma::mat)

Returns:

The 6x6 matrix that represent the symmetric dyadic product (arma::mat)

mat a = randu(3,3);
mat sym_a = 0.5*(a + a.t());
mat c = auto_sym_dyadic_operator(sym_a);

arma::mat sym_dyadic_operator(const arma::mat &a, const arma::mat &b)

Provides the symmetric 4th-order dyadic product of two symmetric tensors.

This function returns the symmetric dyadic product \( \mathbf{C} = \mathbf{a} \odot \mathbf{b} \), defined as:

\[ \left(\mathbf{a} \odot \mathbf{b} \right)_{ijkl} = \frac{1}{2} \left( a_{ik} \, b_{jl} + a_{il} \, b_{jk} \right) \]

The function returns a 6x6 matrix corresponding to a 4th order tensor in Voigt notation (possessing minor symmetry).

Parameters:

a – The first input symmetric tensor (arma::mat)
b – The second input symmetric tensor (arma::mat)

Returns:

The 6x6 matrix that represent the symmetric dyadic product (arma::mat)

mat a = randu(3,3);
mat sym_a = 0.5*(a + a.t());
mat b = randu(3,3);
mat sym_b = 0.5*(b + b.t());
mat c = sym_dyadic_operator(sym_a,sym_b);

arma::mat linearop_eigsym(const arma::vec &b_i, const arma::vec &b_j)

Provides a symmetric 4th-order tensor \( \mathbb{G}^{ij} \), defined such that \( \mathbf{B}_i : \mathbf{D} : \mathbf{B}_j = \mathbb{G}^{ij} : \mathbf{D} \), where \( \mathbf{B}_i \) and \( \mathbf{B}_j \) are eigenprojection tensors.

The 4th-order tensor \( \mathbb{G}^{ij} \) is constructed from eigenvectors \( \mathbf{n}^i \) and \( \mathbf{n}^j \). It is defined such that the triple product can be expressed as a single double contraction:

\[ \mathbf{B}_i : \mathbf{D} : \mathbf{B}_j = \mathbb{G}^{ij} : \mathbf{D} \]

where \( \mathbf{B}_i = \mathbf{n}^i \otimes \mathbf{n}^i \) and \( \mathbf{B}_j = \mathbf{n}^j \otimes \mathbf{n}^j \) are the eigenprojection tensors. Defining \( \mathbf{N} = \mathbf{n}^i \otimes \mathbf{n}^j \), the tensor components are:

\[ \mathbb{G}^{ij}_{klmn} = \frac{1}{2} N_{kl} \left( N_{mn} + N_{nm} \right) \]

This tensor is useful in logarithmic strain computations, particularly for the derivative of the Hencky strain with respect to the Cauchy-Green tensor (see Xiao, Bruhns & Meyers, 1997). The function returns a 6x6 matrix corresponding to a 4th order tensor in Voigt notation.

Parameters:

b_i – The first eigenvector (arma::vec)
b_j – The second eigenvector (arma::vec)

Returns:

The 6x6 matrix that represents the 4th-order tensor \( \mathbb{G}^{ij} \)

contimech.hpp

Overview

API Reference

Overview 

API Reference 