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::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} \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 + a_2^2\,\cos^2 v + a_1^2\,\sin u \right)^2 } \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();
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} \)

mat B = L_Cauchy_Green(F1);

vec bi = zeros(3);
mat Bi;
eig_sym(bi, Bi, B);
mat BBBB = zeros(6,6);

double f_z = 0.;
for (unsigned int i=0; i<3; i++) {
    for (unsigned int j=0; j<3; j++) {
        if ((i!=j)&&(fabs(bi(i)-bi(j))>simcoon::iota)) {
            f_z = (1.+(bi(i)/bi(j)))/(1.-(bi(i)/bi(j)))+2./log(bi(i)/bi(j));
            BBBB = BBBB + f_z*B_klmn(Bi.col(i),Bi.col(j));
        }
    }
}

contimech.hpp

Overview

API Reference

Overview 

API Reference 