Elasticipy.tensors.elasticity

class Elasticipy.tensors.elasticity.ComplianceTensor(C, check_positive_definite=True, mapping=<Elasticipy.tensors.mapping.VoigtMapping object>, **kwargs)[source]

Bases: StiffnessTensor

Class for manipulating compliance tensors

Construct a compliance tensor or an array of commpliance tensors.

The compliance tensor can be constructed from a (6,6) matrix or slices of (6,6) matrices. These matrices must be symmetric. An error is thrown if this matrix in not definite positive (except if check_positive_definite==False, see below). The input argument can also be the full tensor (array of shape (…,3,3,3,3)).

Parameters:

M (list or np.ndarray) – (…,6,6) matrix corresponding to the compliance tensor, written using the Voigt notation, or array of shape (…,3,3,3,3).
phase_name (str, optional) – Phase name to display
check_positive_definite (bool, optional) – Whether to check if the input matrix is positive definite or not. True by default.
check_symmetry (bool, optional) – Whether to check or not that the input matrix is symmetric.
force_symmetry (bool, optional) – If true, the major symmetry of the tensor is forced
mapping (str or MappingConvention) – mapping convention to use. Default is VoigtMapping.

Notes

The units used when building the comliance tensor are up to the user (/GPa, /MPa, /psi etc.). Therefore, the results you will get when performing operations (Young’s modulus, “product” with strain tensor etc.) will be consistent with these units. For instance, if the compliance tensor is defined in /GPa, the computed stress will be given in GPa.

Examples

Create a compliance tensor for cubic symmetry:

>>> matrix = [[200, 40,  40,  0,  0,  0 ],
...           [40,  200, 40,  0,  0,  0 ],
...           [40,  40,  200, 0,  0,  0 ],
...           [0,   0,   0,   20, 0,  0 ],
...           [0,   0,   0,   0,  20, 0 ],
...           [0,   0,   0,   0,   0, 20]]
>>> from Elasticipy.tensors.elasticity import ComplianceTensor
>>> S = ComplianceTensor(matrix)
>>> print(S)
Compliance tensor (in Voigt mapping):
[[200.  40.  40.   0.   0.   0.]
 [ 40. 200.  40.   0.   0.   0.]
 [ 40.  40. 200.   0.   0.   0.]
 [  0.   0.   0.  20.   0.   0.]
 [  0.   0.   0.   0.  20.   0.]
 [  0.   0.   0.   0.   0.  20.]]

Create a stiffness tensor from full (3,3,3,3) array:

>>> S_full = S.full_tensor # (3,3,3,3) numpy array
>>> print((type(S_full), S_full.shape))
(<class 'numpy.ndarray'>, (3, 3, 3, 3))

>>> ComplianceTensor(S_full)
Compliance tensor (in Voigt mapping):
[[200.  40.  40.   0.   0.   0.]
 [ 40. 200.  40.   0.   0.   0.]
 [ 40.  40. 200.   0.   0.   0.]
 [  0.   0.   0.  20.   0.   0.]
 [  0.   0.   0.   0.  20.   0.]
 [  0.   0.   0.   0.   0.  20.]]

Create an array of compliance tensors:

First, we create two slices of (6,6) matrices, corresponding to two compliance values:

>>> slices = [[[200, 40,  40,  0,  0,  0 ],
...            [40,  200, 40,  0,  0,  0 ],
...            [40,  40,  200, 0,  0,  0 ],
...            [0,   0,   0,   20, 0,  0 ],
...            [0,   0,   0,   0,  20, 0 ],
...            [0,   0,   0,   0,   0, 20]],
...           [[250, 80,  80,  0,  0,  0 ],
...            [80,  250, 80,  0,  0,  0 ],
...            [80,  80,  250, 0,  0,  0 ],
...            [0,   0,   0,   40, 0,  0 ],
...            [0,   0,   0,   0,  40, 0 ],
...            [0,   0,   0,   0,   0, 40]]]

Then, one can create an array of compliance tensors:

>>> S_array=ComplianceTensor(slices)
>>> print(S_array)
Compliance tensor (in Voigt mapping):
[[[200.  40.  40.   0.   0.   0.]
  [ 40. 200.  40.   0.   0.   0.]
  [ 40.  40. 200.   0.   0.   0.]
  [  0.   0.   0.  20.   0.   0.]
  [  0.   0.   0.   0.  20.   0.]
  [  0.   0.   0.   0.   0.  20.]]

 [[250.  80.  80.   0.   0.   0.]
  [ 80. 250.  80.   0.   0.   0.]
  [ 80.  80. 250.   0.   0.   0.]
  [  0.   0.   0.  40.   0.   0.]
  [  0.   0.   0.   0.  40.   0.]
  [  0.   0.   0.   0.   0.  40.]]]

This array can be subindexed. E.g.:

>>> S_array[0]
Compliance tensor (in Voigt mapping):
[[200.  40.  40.   0.   0.   0.]
 [ 40. 200.  40.   0.   0.   0.]
 [ 40.  40. 200.   0.   0.   0.]
 [  0.   0.   0.  20.   0.   0.]
 [  0.   0.   0.   0.  20.   0.]
 [  0.   0.   0.   0.   0.  20.]]

Christoffel_tensor(u)[source]

Create the Christoffel tensor along a given direction, or set or directions.

Parameters:: u (list or np.ndarray) – 3D direction(s) to compute the Christoffel tensor along with
Returns:: Christoffel tensor(s). if u is a list of directions, Gamma[i] is the Christoffel tensor for direction u[i].
Return type:: SymmetricSecondOrderTensor

See also

wave_velocity: computes the p- and s-wave velocities.

Notes

For a given stiffness tensor C and a given unit vector u, the Christoffel tensor is defined as [Jaeken] :

\[M_{ij} = C_{iklj}.u_k.u_l\]

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77) # Cubic Cu
>>> C.Christoffel_tensor([0, 0, 1])
Symmetric second-order tensor
[[ 77.   0.   0.]
 [  0.  77.   0.]
 [  0.   0. 186.]]

Hill_average(axis=None)[source]

Compute the Voigt-Reuss-Hill average (mean of Voigt and Reuss averages for stiffness tensors).

If the object is a single tensor, the returned tensor corresponds to the average over an infinite set of random rotations, resulting in an isotropic behavior. Otherwise, the mean is computed over the given axis.

Parameters:: axis (int, optional) – If provided, axis to compute the average along with. If none, the average is computed on the flattened array
Returns:: Voigt-Reuss-Hill average of tensor
Return type:: StiffnessTensor

See also

Voigt_average: compute the Voigt average
Reuss_average: compute the Reuss average
average: generic function for calling either the Voigt, Reuss or Hill average

Examples

Let us consider the stiffness of pure copper:

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)

If we assume that an aggregate is composed of an infinite set grains whose orientations are uniformly distributed, the Voigt average is:

>>> C.Hill_average()
Stiffness tensor (in Voigt mapping):
[[217.83103448 118.08448276 118.08448276   0.           0.
    0.        ]
 [118.08448276 217.83103448 118.08448276   0.           0.
    0.        ]
 [118.08448276 118.08448276 217.83103448   0.           0.
    0.        ]
 [  0.           0.           0.          49.87327586   0.
    0.        ]
 [  0.           0.           0.           0.          49.87327586
    0.        ]
 [  0.           0.           0.           0.           0.
   49.87327586]]

Now, we consider an aggregate of grains whose orientations follow a pure fibre texture along the X axis. First, generate a (large) random set of rotation corresponding to this texture:

>>> from scipy.spatial.transform import Rotation
>>> import numpy as np
>>> g = Rotation.from_euler('X', np.linspace(0, 90, 1000), degrees=True)    # 1000 rotations around X

Now apply rotations to compute the array of stiffness tensors:

>>> C_rotated = C * g
>>> C_rotated
Stiffness tensor array of shape (1000,)

Finally, one can check that the Voigt average, computed from the rotated stiffness tensors is:

>>> C_rotated.Hill_average()
Stiffness tensor (in Voigt mapping):
[[ 1.86000000e+02  1.34000000e+02  1.34000000e+02 -8.62027175e-16
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  2.05164524e+02  1.14835476e+02 -5.48889226e-16
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  1.14835476e+02  2.05164524e+02 -1.44776034e-15
   0.00000000e+00  0.00000000e+00]
 [-8.09683464e-17  8.46209730e-17 -3.98365128e-16  4.52092721e+01
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.70000000e+01 -7.54674101e-17]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.54674101e-17  7.70000000e+01]]

property Poisson_ratio[source]

Directional Poisson’s ratio

Returns:: Poisson’s ratio
Return type:: HyperSphericalFunction

Notes

If the material undergoes tensile strain \(\varepsilon_{ii}\) along the i-th direction, the Poisson ratios are defined as:

\[\nu_{ij}=-\frac{\partial \varepsilon_{jj}}{\partial \varepsilon_{ii}}\]

where \(\varepsilon_{jj}\) denotes the (compressive) longitudinal strain along the j-th direction.

Reuss_average(axis=None)[source]

Compute the Reuss average (from the mean of compliance tensors).

If the object is a single tensor, the returned tensor corresponds to the average over an infinite set of random rotations, resulting in an isotropic behavior. Otherwise, the mean is computed over the given axis.

Parameters:: axis (int, optional) – If provided, axis to compute the average along with. If none, the average is computed on the flattened array
Returns:: Reuss average of stiffness tensor
Return type:: StiffnessTensor

See also

Voigt_average: compute the Voigt average
Hill_average: compute the Voigt-Reuss-Hill average
average: generic function for calling either the Voigt, Reuss or Hill average

Examples

Let us consider the stiffness of pure copper:

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)

If we assume that an aggregate is composed of an infinite set grains whose orientations are uniformly distributed, the Voigt average is:

>>> C.Reuss_average()
Stiffness tensor (in Voigt mapping):
[[208.86206897 122.56896552 122.56896552   0.           0.
    0.        ]
 [122.56896552 208.86206897 122.56896552   0.           0.
    0.        ]
 [122.56896552 122.56896552 208.86206897   0.           0.
    0.        ]
 [  0.           0.           0.          43.14655172   0.
    0.        ]
 [  0.           0.           0.           0.          43.14655172
    0.        ]
 [  0.           0.           0.           0.           0.
   43.14655172]]

Now, we consider an aggregate of grains whose orientations follow a pure fibre texture along the X axis. First, generate a (large) random set of rotation corresponding to this texture:

>>> from scipy.spatial.transform import Rotation
>>> import numpy as np
>>> g = Rotation.from_euler('X', np.linspace(0, 90, 1000), degrees=True)    # 1000 rotations around X

Now apply rotations to compute the array of stiffness tensors:

>>> C_rotated = C * g
>>> C_rotated
Stiffness tensor array of shape (1000,)

Finally, one can check that the Voigt average, computed from the rotated stiffness tensors is:

>>> C_rotated.Reuss_average()
Stiffness tensor (in Voigt mapping):
[[ 1.86000000e+02  1.34000000e+02  1.34000000e+02 -1.68008952e-15
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  1.98854548e+02  1.21145452e+02 -1.09777845e-15
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  1.21145452e+02  1.98854548e+02 -2.89552069e-15
   0.00000000e+00  0.00000000e+00]
 [-1.61936693e-16  1.69241946e-16 -7.96730255e-16  3.88930441e+01
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.70000000e+01 -7.54674101e-17]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.54674101e-17  7.70000000e+01]]

Voigt_average(axis=None)[source]

Compute the Voigt average (from the mean of stiffness tensors).

If the object is a single tensor, the returned tensor corresponds to the average over an infinite set of random rotations, resulting in an isotropic behavior. Otherwise, the mean is computed over the given axis.

Parameters:: axis (int, optional) – If provided, the average is computed along this axis. Otherwise, the mean is computed on the flattened array.
Returns:: Voigt average of stiffness tensor
Return type:: StiffnessTensor

See also

Reuss_average: compute the Reuss average
Hill_average: compute the Voigt-Reuss-Hill average
average: generic function for calling either the Voigt, Reuss or Hill average

Examples

Let us consider the stiffness of pure copper:

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)

If we assume that an aggregate is composed of an infinite set grains whose orientations are uniformly distributed, the Voigt average is:

>>> C.Voigt_average()
Stiffness tensor (in Voigt mapping):
[[226.8 113.6 113.6   0.    0.    0. ]
 [113.6 226.8 113.6   0.    0.    0. ]
 [113.6 113.6 226.8   0.    0.    0. ]
 [  0.    0.    0.   56.6   0.    0. ]
 [  0.    0.    0.    0.   56.6   0. ]
 [  0.    0.    0.    0.    0.   56.6]]

Now, we consider an aggregate of grains whose orientations follow a pure fibre texture along the X axis. First, generate a (large) random set of rotation corresponding to this texture:

>>> from scipy.spatial.transform import Rotation
>>> import numpy as np
>>> g = Rotation.from_euler('X', np.linspace(0, 90, 1000), degrees=True)    # 1000 rotations around X

Now apply rotations to compute the array of stiffness tensors:

>>> C_rotated = C * g
>>> C_rotated
Stiffness tensor array of shape (1000,)

Finally, one can check that the Voigt average, computed from the rotated stiffness tensors is:

>>> C_rotated.Voigt_average()
Stiffness tensor (in Voigt mapping):
[[ 1.86000000e+02  1.34000000e+02  1.34000000e+02 -4.39648318e-17
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  2.11474500e+02  1.08525500e+02  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  1.08525500e+02  2.11474500e+02  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  5.15255000e+01
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.70000000e+01 -7.54674101e-17]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.54674101e-17  7.70000000e+01]]

property Young_modulus[source]

Directional Young’s modulus

Returns:: Young’s modulus
Return type:: SphericalFunction

Zener_ratio(tol=0.0001)[source]

Compute the Zener ratio (Z). Only valid for cubic symmetry.

This function first checks that the tensor has cubic symmetry within a given tolerance. If not, an error is raised.

Parameters:: tol (float, optional) – Tolerance to consider that the material has cubic symmetry
Returns:: Zener ratio (NaN is the symmetry is not cubic)
Return type:: float

Notes

If the tensor is written in canonical base with Voigt mapping, the Zener ratio is defined as:

\[Z=\frac{ 2C_{44} }{C_{11} - C_{12}}\]

The present implementation takes advantage of eigenstiffness to compute the Zener ratio in any base, i.e. even if the tensor is not given in canonical base (e.g. if rotated).

See also

universal_anisotropy: compute the universal anisotropy factor
eig_stiffnesses: eigenstiffnesses of the tensor
is_cubic: check whether the tensor has cubic symmetry or not

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=200, C12=40, C44=20)
>>> C.Zener_ratio()
0.25

which obvisouly corresponds to 2.C44/(C11-C12).

Now, rotate the tensor and see how it looks like:

>>> from scipy.spatial.transform import Rotation
>>> g = Rotation.from_euler('Z', 30, degrees=True)
>>> C_rot = C*g
>>> C_rot
Stiffness tensor (in Voigt mapping):
[[155.          85.          40.           0.           0.
   25.98076211]
 [ 85.         155.          40.           0.           0.
  -25.98076211]
 [ 40.          40.         200.           0.           0.
    0.        ]
 [  0.           0.           0.          20.           0.
    0.        ]
 [  0.           0.           0.           0.          20.
    0.        ]
 [ 25.98076211 -25.98076211   0.           0.           0.
   65.        ]]

Still, we have >>> C_rot.Zener_ratio() 0.24999999999999983

average(method, axis=None)[source]

Compute either the Voigt, Reuss, or Hill average of the stiffness tensor.

This function is just a shortcut for Voigt_average(), Reuss_average(), or Hill_average() and Hill_average().

Parameters:

axis (int, optional) – If provided, axis to compute the average along with. If none, the average is computed on the flattened array
method (str {'Voigt', 'Reuss', 'Hill'}) – Method to use to compute the average.

Return type:

StiffnessTensor

See also

Voigt_average: compute the Voigt average
Reuss_average: compute the Reuss average
Hill_average: compute the Voigt-Reuss-Hill average

property bulk_modulus[source]

Compute the bulk modulus of the material

Returns:: Bulk modulus
Return type:: float or numpy.ndarray

See also

linear_compressibility: directional linear compressibility

copy()[source]: Create a copy of the tensor

classmethod cubic(*, S11=0.0, S12=0.0, S44=0.0, phase_name=None)[source]

Create a fourth-order tensor from cubic symmetry.

Parameters:

S11 (float)
S12 (float)
S44 (float)
phase_name (str, optional) – Phase name to display

Return type:

StiffnessTensor

See also

hexagonal: create a tensor from hexagonal symmetry
orthorhombic: create a tensor from orthorhombic symmetry

ddot(other, mode='pair')[source]

Perform tensor product contracted twice (“:”) between two fourth-order tensors

Parameters:

other (FourthOrderTensor or SecondOrderTensor) – Right-hand side of “:” symbol
mode (str, optional) – If mode==”pair”, the tensors must be broadcastable, and the tensor product are performed on the last axes. If mode==”cross”, all cross-combinations are considered.

Returns:

Result from double-contraction

Return type:

FourthOrderTensor

Examples

First, let consider two random arrays of Fourth-order tensors:

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> T1 = FourthOrderTensor.rand(shape=(2,3))
>>> T2 = FourthOrderTensor.rand(shape=3)
>>> T1T2_pair = T1.ddot(T2)
>>> T1T2_pair
4th-order tensor array of shape (2, 3)

whereas:

>>> T1T2_cross = T1.ddot(T2, mode='cross')
>>> T1T2_cross
4th-order tensor array of shape (2, 3, 3)

The command above is equivalent (but way faster) to:

>>> T1T2_cross_loop = FourthOrderTensor.zeros(shape=(2,3,3))
>>> for i in range(2):
...     for j in range(3):
...         for k in range(3):
...             T1T2_cross_loop[i,j,k] = T1[i,j].ddot(T2[k])

One can check that the results are consistent with:

>>> T1T2_cross_loop == T1T2_cross
array([[[ True,  True,  True],
        [ True,  True,  True],
        [ True,  True,  True]],

       [[ True,  True,  True],
        [ True,  True,  True],
        [ True,  True,  True]]])

deviatoric_part()[source]

Return the deviatoric part of the tensor

Returns:: Deviatoric part of the tensor
Return type:: FourthOrderTensor

See also

identity_tensor: return the identity tensor
spherical_part: return the spherical part of the tensor

eig()[source]

Compute the eigencompliances and the eigenstresses.

Solve the eigenvalue problem from the Kelvin matrix of the compliance tensor (see Notes).

Returns:

numpy.ndarray – Array of 6 eigencompliances (eigenvalues of the stiffness matrix)
numpy.ndarray – (6,6) array of eigenstresses (eigenvectors of the stiffness matrix)

See also

Kelvin: returns the stiffness components as a (6,6) matrix, according to the Kelvin mapping convention.
eig_compliances: returns the eigencompliances only
eig_stresses: returns the eigenstresses only

Notes

The definition for eigencompliances and the eigenstresses are introduced in [Helbig].

property eig_compliances[source]

Compute the eigencompliances given by the Kelvin’s matrix for stiffness.

Returns:: 6 eigenvalues of the Kelvin’s compliance matrix, in ascending order
Return type:: numpy.ndarray

See also

eig: returns the eigencompliances and the eigenstresses
eig_strains: returns the eigenstresses only

property eig_stiffnesses[source]

Compute the eigenstiffnesses from the Kelvin’s matrix of compliance

Returns:: inverses of 6 eigenvalues of the Kelvin’s compliance matrix, in descending order
Return type:: numpy.ndarray

See also

eig_compliances: compute the eigencompliances from the Kelvin’s matrix of compliance

eig_stiffnesses_multiplicity(tol=0.0001)[source]

Compute the eigenstiffnesses, then returns the multiplicity of each eigenstiffness.

Given an absolute tolerance, duplicates in eigenstiffnesses are considered to compute the multiplicity of each value.

Parameters:

tol (float, optional) – Absolute tolerance to assume that two distinct eigenstiffnesses are the same

Returns:

numpy.ndarray – Unique values of eigenstiffnesses, sorted by increasing multiplicity
numpy.ndarray – Multiplicity of each unique eigenstiffness, sorted by ascending value

See also

eig_stiffnesses: compute the eigenstiffnesses

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)
>>> C.eig_stiffnesses
array([ 52.,  52., 154., 154., 154., 454.])
>>> C.eig_stiffnesses_multiplicity()
(array([ 52., 154., 454.]), array([2, 3, 1]))

property eig_strains[source]

Compute the eigenstrains from the Kelvin’s matrix for stiffness

Returns:: (6,6) matrix of eigenstrains, sorted by ascending order of eigenstiffnesses.
Return type:: numpy.ndarray

See also

eig: returns both the eigenvalues and the eigenvectors of the Kelvin matrix

property eig_stresses[source]

Compute the eigenstresses from the Kelvin’s matrix for stiffness

Returns:: (6,6) matrix of eigenstresses, sorted by ascending order of eigencompliances.
Return type:: numpy.ndarray

See also

eig: returns both the eigencompliances and the eigenstresses

classmethod eye(shape=(), **kwargs)[source]

Create a 4th-order identity tensor.

See notes for definition.

Parameters:

shape (int or tuple, optional) – Shape of the tensor to create
mapping (Kelvin mapping, optional) – Mapping convention to use. Must be either Kelvin or Voigt.

Returns:

Identity tensor

Return type:

FourthOrderTensor or SymmetricFourthOrderTensor

Notes

The Fourth-order identity tensor is defined as:

\[I_{ijkl} = \frac12\left( \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}\right)\]

Examples

Create a (single) identity tensor:

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> I = FourthOrderTensor.eye()
>>> print(I)
4th-order tensor (in Kelvin mapping):
[[1. 0. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 1.]]

Alternatively, one can use another mapping convention, e.g. Voigt:

>>> from Elasticipy.tensors.mapping import VoigtMapping
>>> Iv = FourthOrderTensor.eye(mapping=VoigtMapping())
>>> print(Iv)
4th-order tensor (in Voigt mapping):
[[1.  0.  0.  0.  0.  0. ]
 [0.  1.  0.  0.  0.  0. ]
 [0.  0.  1.  0.  0.  0. ]
 [0.  0.  0.  0.5 0.  0. ]
 [0.  0.  0.  0.  0.5 0. ]
 [0.  0.  0.  0.  0.  0.5]]

Still, we have:

>>> I == Iv
True

as they correspond to the same tensor, but expressed as a matrix with different mapping conventions. Indeed, one can check that:

>>> import numpy as np
>>> np.array_equal(I.full_tensor, Iv.full_tensor)
True

flatten()[source]

Flatten the tensor

If the tensor array is of shape (m,n,o…,r), the flattened array will be of shape (m*n*o*…*r,).

Returns:: Flattened tensor
Return type:: SymmetricFourthOrderTensor

Examples

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> T = FourthOrderTensor.rand(shape=(5,6))
>>> T
4th-order tensor array of shape (5, 6)
>>> T.flatten()
4th-order tensor array of shape (30,)

classmethod fromCrystalSymmetry(symmetry='Triclinic', point_group=None, diad='y', phase_name=None, prefix=None, **kwargs)[source]

Create a fourth-order tensor from limited number of components, taking advantage of crystallographic symmetries

Parameters:

symmetry (str, default Triclinic) – Name of the crystallographic symmetry
point_group (str) – Point group of the considered crystal. Only used (and mandatory) for tetragonal and trigonal symmetries.
diad (str {'x', 'y'}, default 'x') – Alignment convention. Sets whether x||a or y||b. Only used for monoclinic symmetry.
phase_name (str, default None) – Name to use when printing the tensor
prefix (str, default None) – Define the prefix to use when providing the components. By default, it is ‘C’ for stiffness tensors, ‘S’ for compliance.
kwargs – Keywords describing all the necessary components, depending on the crystal’s symmetry and the type of tensor. For Stiffness, they should be named as ‘Cij’ (e.g. C11=…, C12=…). For Comliance, they should be named as ‘Sij’ (e.g. S11=…, S12=…). See examples below. The behaviour can be overriten with the prefix option (see above)

Return type:

FourthOrderTensor

See also

isotropic: creates an isotropic stiffness tensor from two paremeters (e.g. E and v).

Notes

The relationships between the tensor’s components depend on the crystallogrpahic symmetry [Nye].

References

[Nye]

Nye, J. F. Physical Properties of Crystals. London: Oxford University Press, 1959.

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor

>>> StiffnessTensor.fromCrystalSymmetry(symmetry='monoclinic', diad='y', phase_name='TiNi',
...                                     C11=231, C12=127, C13=104,
...                                     C22=240, C23=131, C33=175,
...                                     C44=81, C55=11, C66=85,
...                                     C15=-18, C25=1, C35=-3, C46=3)
Stiffness tensor (in Voigt mapping):
[[231. 127. 104.   0. -18.   0.]
 [127. 240. 131.   0.   1.   0.]
 [104. 131. 175.   0.  -3.   0.]
 [  0.   0.   0.  81.   0.   3.]
 [-18.   1.  -3.   0.  11.   0.]
 [  0.   0.   0.   3.   0.  85.]]
Phase: TiNi

>>> from Elasticipy.tensors.elasticity import ComplianceTensor

>>> ComplianceTensor.fromCrystalSymmetry(symmetry='monoclinic', diad='y', phase_name='TiNi',
...                                      S11=8, S12=-3, S13=-2,
...                                      S22=8, S23=-5, S33=10,
...                                      S44=12, S55=116, S66=12,
...                                      S15=14, S25=-8, S35=0, S46=0)
Compliance tensor (in Voigt mapping):
[[  8.  -3.  -2.   0.  14.   0.]
 [ -3.   8.  -5.   0.  -8.   0.]
 [ -2.  -5.  10.   0.   0.   0.]
 [  0.   0.   0.  12.   0.   0.]
 [ 14.  -8.   0.   0. 116.   0.]
 [  0.   0.   0.   0.   0.  12.]]
Phase: TiNi

classmethod from_Kelvin(matrix, **kwargs)[source]

Create a tensor from the (6,6) matrix following the Kelvin(-Mandel) mapping convention

Parameters:

matrix (list or numpy.ndarray) – (6,6) matrix of components
kwargs – keyword arguments passed to the constructor

Return type:

StiffnessTensor

See also

to_Kelvin: return the components as a (6,6) matrix following the Kelvin convention

classmethod from_MP(ids, api_key=None)[source]

Import stiffness tensor(s) from the Materials Project API, given their material ids.

You need to register to https://materialsproject.org first to get an API key. This key can be explicitly passed as an argument (see below), or provided as an environment variable named MP_API_KEY.

Parameters:

ids (str or list of str) – ID(s) of the material to import (e.g. “mp-1048”)
api_key (str, optional) – API key to the Materials Project API. If not provided, it should be available as the API_KEY environment variable.

Returns:

If one of the requested material ids was not found, the corresponding value in the list will be None.

Return type:

list of StiffnessTensor

classmethod from_txt_file(filename)[source]

Load the tensor from a text file.

The two first lines can have data about phase name and symmetry, but this is not mandatory.

Parameters:: filename (str) – Filename to load the tensor from.
Returns:: The reconstructed tensor read from the file.
Return type:: SymmetricFourthOrderTensor

See also

save_to_txt: create a tensor from text file

property full_tensor[source]

Returns the full (unvoigted) tensor as a (3, 3, 3, 3) or (…, 3, 3, 3, 3) array

Returns:: Full tensor (4-index notation)
Return type:: np.ndarray

Examples

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> I = FourthOrderTensor.eye() # 4th order identity tensor
>>> print(I)
4th-order tensor (in Kelvin mapping):
[[1. 0. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 1.]]

>>> I_full = I.full_tensor
>>> type(I_full)
<class 'numpy.ndarray'>
>>> I_full.shape
(3, 3, 3, 3)

When working on tensor arrays, the shape of the resulting numpy array will change accordlingly. E.g.:

>>> I_array = FourthOrderTensor.eye(shape=(5,6)) # Array of 4th order identity tensor
>>> I_array.full_tensor.shape
(5, 6, 3, 3, 3, 3)

classmethod hexagonal(*, S11=0.0, S12=0.0, S13=0.0, S33=0.0, S44=0.0, phase_name=None)[source]

Create a fourth-order tensor from hexagonal symmetry.

Parameters:

S11 (float) – Components of the tensor, using the Voigt notation
S12 (float) – Components of the tensor, using the Voigt notation
S13 (float) – Components of the tensor, using the Voigt notation
S33 (float) – Components of the tensor, using the Voigt notation
S44 (float) – Components of the tensor, using the Voigt notation
phase_name (str, optional) – Phase name to display

Return type:

FourthOrderTensor

See also

transverse_isotropic: creates a transverse-isotropic tensor from engineering parameters
cubic: create a tensor from cubic symmetry
tetragonal: create a tensor from tetragonal symmetry

classmethod identity(**kwargs)[source]

Construct the Fourth-order identity tensor.

This is actually an alias for eye().

Parameters:: kwargs – Keyword arguments passed to the Fourth-order tensor constructor.
Return type:: FourthOrderTensor

See also

eye: Fourth-order identity tensor

classmethod identity_deviatoric_part(**kwargs)[source]

Return the deviatoric part of the identity tensor.

See notes for the mathematical definition.

Parameters:: kwargs – keyword arguments passed to eye constructor
Return type:: FourthOrderTensor or SymmetricTensor

See also

identity_tensor: return the identity tensor
identity_spherical_part: return the spherical part of the identity tensor

Notes

The deviatoric part of the identity tensor is defined as:

\[K = I - J\]

where \(I\) and \(J\) denote the identity and the deviatoric part of the identity tensor, respectively.

Examples

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> K = FourthOrderTensor.identity_deviatoric_part()
>>> print(K)
4th-order tensor (in Kelvin mapping):
[[ 0.66666667 -0.33333333 -0.33333333  0.          0.          0.        ]
 [-0.33333333  0.66666667 -0.33333333  0.          0.          0.        ]
 [-0.33333333 -0.33333333  0.66666667  0.          0.          0.        ]
 [ 0.          0.          0.          1.          0.          0.        ]
 [ 0.          0.          0.          0.          1.          0.        ]
 [ 0.          0.          0.          0.          0.          1.        ]]

One can check that K has zero spherical part:

>>> print(K.spherical_part())
4th-order tensor (in Kelvin mapping):
[[2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]]

classmethod identity_spherical_part(shape=(), **kwargs)[source]

Return the spherical part of the identity tensor.

See Notes for mathematical definition.

Parameters:

shape (tuple of int, optional) – Shape of the tensor to create
kwargs – Keyword arguments passed to the Fourth-order tensor constructor.

Return type:

FourthOrderTensor

See also

identity_tensor: return the identity tensor
identity_deviatoric_part: return the deviatoric part of the identity tensor

Notes

The spherical part of the identity tensor is defined as:

\[J_{ijkl} = \frac13 \delta_{ij}\delta_{kl}\]

Examples

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> J = FourthOrderTensor.identity_spherical_part()
>>> print(J)
4th-order tensor (in Kelvin mapping):
[[0.33333333 0.33333333 0.33333333 0.         0.         0.        ]
 [0.33333333 0.33333333 0.33333333 0.         0.         0.        ]
 [0.33333333 0.33333333 0.33333333 0.         0.         0.        ]
 [0.         0.         0.         0.         0.         0.        ]
 [0.         0.         0.         0.         0.         0.        ]
 [0.         0.         0.         0.         0.         0.        ]]

On can check that J has zero deviatoric part:

>>> J.deviatoric_part()
4th-order tensor (in Kelvin mapping):
[[2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]]

infinite_random_average()[source]

Compute the average of the tensor, assuming that an infinite number of random orientations is applied.

Returns:: Average tensor or tensor array. The returned array will be of the same shape as the input object.
Return type:: FourthOrderTensor

inv()[source]

Compute the reciprocal stiffness tensor

Returns:: Reciprocal tensor
Return type:: StiffnessTensor

is_cubic(tol=0.01)[source]

Check that the tensor corresponds to cubic symmetry, within a given tolerance.

The method relies on the multiplicity of eigenstiffnesses.

Parameters:: tol (float, optional) – Absolute tolerance to consider multiplicity of eigenstiffnesses
Returns:: If the tensor is single, the returned value is boolean. If the object is a tensor array, the returned value is an array of bools, the same shape as the tensor array.
Return type:: bool or numpy.ndarray

See also

is_isotropic: check if the stiffness tensor is isotropic
is_tetragonal: check if the stiffness tensor has tetragonal symmetry
eig_stiffnesses: compute eigenstiffnesses

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> from scipy.spatial.transform import Rotation
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)
>>> C_rotated = C * Rotation.random(random_state=123)
>>> C_rotated
Stiffness tensor (in Voigt mapping):
[[237.71171578  96.41409344 119.87419078   8.1901353   -3.63846312
  -20.34233446]
 [ 96.41409344 250.74909842 106.83680814   9.33462785  -6.52548033
    0.99714278]
 [119.87419078 106.83680814 227.28900108 -17.52476315  10.16394345
   19.34519167]
 [  8.1901353    9.33462785 -17.52476315  49.83680814  19.34519167
   -6.52548033]
 [ -3.63846312  -6.52548033  10.16394345  19.34519167  62.87419078
    8.1901353 ]
 [-20.34233446   0.99714278  19.34519167  -6.52548033   8.1901353
   39.41409344]]

Once rotated, it is not clear if the stiffness tensors has cubic symmetry on sight. Yet:

>>> C_rotated.is_cubic()
True

is_isotropic(tol=0.01)[source]

Check that the tensor corresponds to isotropic symmetry, within a given tolerance.

The method relies on the multiplicity of eigenstiffnesses

Parameters:: tol (float) – Absolute tolerance to consider multiplicity of eigenstiffnesses
Returns:: If the tensor is single, the returned value is boolean. If the object is a tensor array, the returned value is an array of bools, the same shape as the tensor array.
Return type:: bool or numpy.ndarray

See also

is_cubic: check if the stiffness tensor has cubic symmetry
is_tetragonal: check if the stiffness tensor has tetragonal symmetry
eig_stiffnesses: compute eigenstiffnesses

is_tetragonal(tol=0.01)[source]

Check that the tensor corresponds to tetragonal symmetry, within a given tolerance.

The method relies on the multiplicity of eigenstiffnesses.

Parameters:: tol (float) – Absolute tolerance to consider multiplicity of eigenstiffnesses
Returns:: If the tensor is single, the returned value is boolean. If the object is a tensor array, the returned value is an array of bools, the same shape as the tensor array.
Return type:: bool or numpy.ndarray

See also

is_isotropic: check if the stiffness tensor is isotropic
is_cubic: check if the stiffness tensor has cubic symmetry
eig_stiffnesses: compute eigenstiffnesses

classmethod isotropic(E=None, nu=None, G=None, lame1=None, lame2=None, K=None, phase_name=None)[source]

Create an isotropic stiffness tensor.

The stiffness tensor must be constructed from exactly two elasticity coefficients, namely: E, nu, G, lame1, or lame2. Note that lame2 is just an alias for G. Each of these coefficient can be a list; in this case, the returned object is a tensor array (see examples).

Parameters:

E (float or list, None) – Young modulus
nu (float or list, None) – Poisson ratio
G (float or list, None) – Shear modulus
lame1 (float or list, None) – First Lamé coefficient
lame2 (float or list, None) – Second Lamé coefficient (alias for G)
K (float or list, None) – Bulk modulus
phase_name (str, None) – Name to print

Return type:

Corresponding isotropic stiffness tensor

See also

transverse_isotropic: create a transverse-isotropic tensor

Notes

The units you use when passing the elastic moduli must be consistent with that of the stress tensor. For instance, if you expect to work in MPa, the Young’s modulus and the Lamé’s coefficient must be given in MPa as well.

Examples

On can check that the shear modulus for steel is around 82 GPa:

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C=StiffnessTensor.isotropic(E=210e3, nu=0.28)
>>> C.shear_modulus
Hyperspherical function
Min=82031.24999999991, Max=82031.25000000006

Similarly, the same tensor can be defined from another pair of parameters, e.g. Young and shear moduli:

>>> C=StiffnessTensor.isotropic(E=210e3, G=82031)
>>> C.Young_modulus
Spherical function
Min=210000.0000000001, Max=210000.00000000032

One can build a tensor array by providing a list of values for the input arguments, instead of floats. For instance:

>>> C = StiffnessTensor.isotropic(E=(210, 70), nu=(0.28, 0.35)) # Elastic moduli for steel and aluminium
>>> C.shape
(2,)

We can easily check that the shear moduli for steel and aluminium are:

>>> C[0].shear_modulus
Hyperspherical function
Min=82.0312499999999, Max=82.03125000000006

>>> C[1].shear_modulus
Hyperspherical function
Min=25.925925925925895, Max=25.925925925925952

property lame1[source]

Compute the first Lamé’s parameter (only for isotropic materials).

If the stiffness/compliance tensor is not isotropic, NaN is returned.

Returns:: First Lamé’s parameter
Return type:: float

See also

lame2: second Lamé’s parameter

property lame2[source]

Compute the second Lamé’s parameter (only for isotropic materials).

If the stiffness/compliance tensor is not isotropic, NaN is returned.

Returns:: Second Lamé’s parameter
Return type:: float

See also

lame1: first Lamé’s parameter

property linear_compressibility[source]

Compute the directional linear compressibility.

Returns:: Directional linear compressibility
Return type:: SphericalFunction

See also

bulk_modulus: bulk modulus of the material

linear_invariants()[source]

Compute the linear invariants of the tensor, or tensor array.

If the object is a tensor array, the linear invariants are returned as arrays of each invariant. See notes for the actual definitions.

Returns:

A1 (float or np.ndarray) – First linear invariant
A2 (float or np.ndarray) – Second linear invariant

See also

quadratic_invariants: compute the quadratic invariants of a fourth-order tensor

Notes

The linear invariants are:

\[ \begin{align}\begin{aligned}A_1=C_{ijij}\\A_2=C_{iijj}\end{aligned}\end{align} \]

matrix(mapping_convention=None)[source]

Returns the components of the tensor as a matrix.

Parameters:: mapping_convention (VoigtMapping, optional) – Mapping convention to use for the returned matrix. If not provided, that of the tensor is used.
Returns:: Components of the tensor as a matrix
Return type:: numpy.ndarray

Examples

Create an identity 4th-order tensor:

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> t = FourthOrderTensor.eye()

Its matrix with respect to Kelvin mapping is:

>>> t.matrix()
array([[1., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 1.]])

whereas, when using the Voigt mapping, we have:

>>> from Elasticipy.tensors.mapping import VoigtMapping
>>> t.matrix(mapping_convention=VoigtMapping())
array([[1. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 1. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 1. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0.5, 0. , 0. ],
       [0. , 0. , 0. , 0. , 0.5, 0. ],
       [0. , 0. , 0. , 0. , 0. , 0.5]])

For stiffness tensors, the default mapping convention is Voigt, so that:

>>> from Elasticipy.tensors.elasticity import StiffnessTensor, ComplianceTensor
>>> StiffnessTensor.eye().matrix()
array([[1. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 1. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 1. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0.5, 0. , 0. ],
       [0. , 0. , 0. , 0. , 0.5, 0. ],
       [0. , 0. , 0. , 0. , 0. , 0.5]])

whereas for compliance tensor, the default mapping convention gives:

>>> ComplianceTensor.eye().matrix()
array([[1., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 2., 0., 0.],
       [0., 0., 0., 0., 2., 0.],
       [0., 0., 0., 0., 0., 2.]])

mean(axis=None)[source]

Compute the mean value of the tensor T

Parameters:: axis (int or list of int or tuple of int, optional) – axis along which to compute the mean. If None, the mean is computed on the flattened tensor
Returns:: If no axis is given, the result will be of shape (3,3,3,3).
Return type:: numpy.ndarray

Examples

Create a random tensor array of shape (5,6):

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> T = FourthOrderTensor.rand(shape=(5,6))
>>> Overall_mean = T.mean()
>>> Overall_mean.shape
()
>>> Overall_mean
4th-order tensor (in Kelvin mapping):
[[0.514295   0.52259217 0.42899181 0.77148692 0.64073221 0.73211491]
 [0.49422678 0.43718365 0.40786118 0.8170971  0.68435571 0.67262655]
 [0.48753674 0.51142541 0.44650454 0.76310921 0.67724973 0.69430165]
 [0.53946846 0.75101474 0.73578098 1.04338905 1.21598419 0.99489014]
 [0.75354555 0.61193555 0.82341479 1.11197826 0.89183143 1.20986243]
 [0.66078807 0.70126535 0.63719147 0.87567139 1.05671229 1.03004098]]

>>> axis_0_mean = T.mean(axis=0)
>>> axis_0_mean.shape
(6,)
>>> axis_1_mean = T.mean(axis=1)
>>> axis_1_mean.shape
(5,)

classmethod monoclinic(*, S11=0.0, S12=0.0, S13=0.0, S22=0.0, S23=0.0, S33=0.0, S44=0.0, S55=0.0, S66=0.0, S15=None, S25=None, S35=None, S46=None, S16=None, S26=None, S36=None, S45=None, phase_name=None)[source]

Create a fourth-order tensor from monoclinic symmetry. It automatically detects whether the components are given according to the Y or Z diad, depending on the input arguments.

For Diad || y, S15, S25, S35 and S46 must be provided. For Diad || z, S16, S26, S36 and S45 must be provided.

Parameters:

S11 (float) – Components of the tensor, using the Voigt notation
S12 (float) – Components of the tensor, using the Voigt notation
S13 (float) – Components of the tensor, using the Voigt notation
S22 (float) – Components of the tensor, using the Voigt notation
S23 (float) – Components of the tensor, using the Voigt notation
S33 (float) – Components of the tensor, using the Voigt notation
S44 (float) – Components of the tensor, using the Voigt notation
S55 (float) – Components of the tensor, using the Voigt notation
S66 (float) – Components of the tensor, using the Voigt notation
S15 (float, optional) – S15 component of the tensor (if Diad || y)
S25 (float, optional) – S25 component of the tensor (if Diad || y)
S35 (float, optional) – S35 component of the tensor (if Diad || y)
S46 (float, optional) – S46 component of the tensor (if Diad || y)
S16 (float, optional) – S16 component of the tensor (if Diad || z)
S26 (float, optional) – S26 component of the tensor (if Diad || z)
S36 (float, optional) – S36 component of the tensor (if Diad || z)
S45 (float, optional) – S45 component of the tensor (if Diad || z)
phase_name (str, optional) – Name to display

Return type:

FourthOrderTensor

See also

triclinic: create a tensor from triclinic symmetry
orthorhombic: create a tensor from orthorhombic symmetry

property ndim[source]

Returns the dimensionality of the tensor (number of dimensions in the orientation array)

Returns:: Number of dimensions
Return type:: int

classmethod ones(shape=None, **kwargs)[source]

Create a 4th-order tensor full of ones.

Parameters:

shape (int or tuple, optional) – Shape of the tensor to create
kwargs – keyword arguments passed to the constructor

Return type:

FourthOrderTensor

Examples

>>> tensor_of_ones = FourthOrderTensor.ones()
>>> tensor_of_ones
4th-order tensor (in Kelvin mapping):
[[1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]]

At first sight, the tensor may appear not full of ones at all, but the representation above uses the Kelvin mapping convention. Indeed, one can check that the full tensor is actually full of ones. E.g.:

>>> tensor_of_ones.full_tensor[0,1,0,2]
1.0

Alternatively, the Voigt mapping convention may help figuring it out:

>>> from Elasticipy.tensors.mapping import VoigtMapping
>>> tensor_of_ones_voigt = FourthOrderTensor.ones(mapping=VoigtMapping())
>>> tensor_of_ones_voigt
4th-order tensor (in Voigt mapping):
[[1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]]

although both tensors are actually the same:

>>> tensor_of_ones == tensor_of_ones_voigt
True

classmethod orthorhombic(*, S11=0.0, S12=0.0, S13=0.0, S22=0.0, S23=0.0, S33=0.0, S44=0.0, S55=0.0, S66=0.0, phase_name=None)[source]

Create a fourth-order tensor from orthorhombic symmetry.

Parameters:

S11 (float) – Components of the tensor, using the Voigt notation
S12 (float) – Components of the tensor, using the Voigt notation
S13 (float) – Components of the tensor, using the Voigt notation
S22 (float) – Components of the tensor, using the Voigt notation
S23 (float) – Components of the tensor, using the Voigt notation
S33 (float) – Components of the tensor, using the Voigt notation
S44 (float) – Components of the tensor, using the Voigt notation
S55 (float) – Components of the tensor, using the Voigt notation
S66 (float) – Components of the tensor, using the Voigt notation
phase_name (str, optional) – Phase name to display

Return type:

FourthOrderTensor

See also

monoclinic: create a tensor from monoclinic symmetry
orthorhombic: create a tensor from orthorhombic symmetry

classmethod orthotropic(*, Ex, Ey, Ez, Gxy, Gxz, Gyz, nu_yx=None, nu_zx=None, nu_zy=None, nu_xy=None, nu_xz=None, nu_yz=None, **kwargs)[source]

Create a stiffness tensor corresponding to orthotropic symmetry, given the engineering constants.

Exactly three Poisson ratios must be provided. See Notes for details.

Parameters:

Ex (float) – Young modulus along the x axis
Ey (float) – Young modulus along the y axis
Ez (float) – Young modulus along the z axis
Gxy (float) – Shear modulus in the x-y plane
Gxz (float) – Shear modulus in the x-z plane
Gyz (float) – Shear modulus in the y-z plane
nu_xy (float, optional) – Poisson ratio along x and y axes. Either nu_xy or nu_yx must be provided, not both.
nu_yx (float, optional) – Poisson ratio along x and y axes. Either nu_xy or nu_yx must be provided, not both.
nu_xz (float, optional) – Poisson ratio along x and z axes. Either nu_xz or nu_zx must be provided, not both.
nu_zx (float, optional) – Poisson ratio along x and z axes. Either nu_xz or nu_zx must be provided, not both.
nu_yz (float, optional) – Poisson ratio along y and z axes. Either nu_yz or nu_zy must be provided, not both.
nu_zy (float, optional) – Poisson ratio along y and z axes. Either nu_yz or nu_zy must be provided, not both.
kwargs (dict, optional) – Keyword arguments to pass to the StiffnessTensor constructor

Return type:

StiffnessTensor

See also

transverse_isotropic: create a stiffness tensor for transverse-isotropic symmetry

Notes

If the material undergoes tensile strain \(\varepsilon_{ii}\) along the i-th direction, the Poisson ratios are defined as:

\[\nu_{ij}=-\frac{\partial \varepsilon_{jj}}{\partial \varepsilon_{ii}}\]

where \(\varepsilon_{jj}\) denotes the (compressive) longitudinal strain along the j-th direction. If \(E_x\) and \(E_y\) are the Young moduli along x and y, we have:

\[\frac{\nu_{xy}}{E_x} = \frac{\nu_{yx}}{E_y}\]

quadratic_invariants()[source]

Compute the quadratic invariants of the tensor, or tensor array.

If the object is a tensor array, the returned values are arrays of each invariant. See notes for definitions.

Returns:: B1, B2, B3, B4, B5
Return type:: float or np.ndarray

See also

linear_invariants: compute the linear invariants of a Fourth-order tensor

Notes

The quadratic invariants are defined as [Norris]:

\[ \begin{align}\begin{aligned}B_1 = C_{ijkl}C_{ijkl}\\B_2 = C_{iikl}C_{jjkl}\\B_3 = C_{iikl}C_{jkjl}\\B_4 = C_{kiil}C_{kjjl}\\B_5 = C_{ijkl}C_{ikjl}\end{aligned}\end{align} \]

References

[Norris]

Norris, A. N. (22 May 2007). “Quadratic invariants of elastic moduli”. The Quarterly Journal of Mechanics and Applied Mathematics. 60 (3): 367–389. doi:10.1093/qjmam/hbm007

classmethod rand(shape=None, **kwargs)[source]

Populate a Fourth-order tensor with random values in half-open interval [0.0, 1.0).

Parameters:

shape (tuple or int, optional) – Set the shape of the tensor array. If None, the returned tensor will be single.
kwargs – Keyword arguments passed to the Fourth-order tensor constructor.

Returns:

Fourth-order tensor

Return type:

FourthOrderTensor

rotate(rotation)[source]

Apply a single rotation to a tensor, and return its component into the rotated frame.

Parameters:: rotation (Rotation or orix.quaternion.rotation.Rotation) – Rotation to apply
Returns:: Rotated tensor
Return type:: SymmetricFourthOrderTensor

Examples

Let start from a given tensor, (say ones):

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> T = FourthOrderTensor.ones()
>>> T
4th-order tensor (in Kelvin mapping):
[[1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]]

Define a rotation. E.g.:

>>> from scipy.spatial.transform import Rotation
>>> g = Rotation.from_euler('X', 90, degrees=True)

Then , apply rotation:

>>> Trotated = T.rotate(g)
>>> Trotated
4th-order tensor (in Kelvin mapping):
[[ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [-1.41421356 -1.41421356 -1.41421356  2.         -2.          2.        ]
 [ 1.41421356  1.41421356  1.41421356 -2.          2.         -2.        ]
 [-1.41421356 -1.41421356 -1.41421356  2.         -2.          2.        ]]

Actually, a more simple syntax is:

>>> T * g
4th-order tensor (in Kelvin mapping):
[[ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [-1.41421356 -1.41421356 -1.41421356  2.         -2.          2.        ]
 [ 1.41421356  1.41421356  1.41421356 -2.          2.         -2.        ]
 [-1.41421356 -1.41421356 -1.41421356  2.         -2.          2.        ]]

Obviously, the original tensor can be retrieved by applying the reverse rotation:

>>> Trotated * g.inv()
4th-order tensor (in Kelvin mapping):
[[1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]]

If g is composed of multiple rotations, this will result in a tensor array, corresponding to each rotation:

>>> import numpy as np
>>> theta = np.linspace(0, 90, 100)
>>> g = Rotation.from_euler('X', theta, degrees=True)
>>> Trotated = T * g
>>> Trotated
4th-order tensor array of shape (100,)

save_to_txt(filename, matrix_only=False)[source]

Save the tensor to a text file.

Parameters:

filename (str) – Filename to save the tensor to.
matrix_only (bool, False) – If true, only the components of the stiffness tensor is saved (no data about phase nor symmetry)

See also

from_txt_file: create a tensor from text file

property shape[source]

Return the shape of the tensor array

Returns:: Shape of the tensor array
Return type:: tuple

property shear_modulus[source]

Directional shear modulus

Returns:: Shear modulus
Return type:: HyperSphericalFunction

spherical_part()[source]

Return the spherical part of the tensor

Returns:: Spherical part of the tensor
Return type:: FourthOrderTensor

See also

identity_tensor: return the identity tensor
deviatoric_part: return the deviatoric part of the tensor

classmethod tetragonal(*, S11=0.0, S12=0.0, S13=0.0, S33=0.0, S44=0.0, S16=0.0, S66=0.0, phase_name=None)[source]

Create a fourth-order tensor from tetragonal symmetry.

Parameters:

S11 (float) – Components of the tensor, using the Voigt notation
S12 (float) – Components of the tensor, using the Voigt notation
S13 (float) – Components of the tensor, using the Voigt notation
S33 (float) – Components of the tensor, using the Voigt notation
S44 (float) – Components of the tensor, using the Voigt notation
S66 (float) – Components of the tensor, using the Voigt notation
S16 (float, optional) – S16 component in Voigt notation (for point groups 4, -4 and 4/m only)
phase_name (str, optional) – Phase name to display

Return type:

FourthOrderTensor

See also

trigonal: create a tensor from trigonal symmetry
orthorhombic: create a tensor from orthorhombic symmetry

to_Kelvin()[source]

Returns all the tensor components using the Kelvin(-Mandel) mapping convention.

Returns:: (6,6) matrix, according to the Kelvin mapping
Return type:: numpy.ndarray

See also

eig: returns the eigenvalues and the eigenvectors of the Kelvin’s matrix
from_Kelvin: Construct a fourth-order tensor from its (6,6) Kelvin matrix

Notes

This mapping convention is defined as follows [Helbig]:

\[\begin{split}C_K = \begin{bmatrix} C_{11} & C_{12} & C_{13} & \sqrt{2}C_{14} & \sqrt{2}C_{15} & \sqrt{2}C_{16}\\ C_{12} & C_{22} & C_{23} & \sqrt{2}C_{24} & \sqrt{2}C_{25} & \sqrt{2}C_{26}\\ C_{13} & C_{23} & C_{33} & \sqrt{2}C_{34} & \sqrt{2}C_{35} & \sqrt{2}C_{36}\\ \sqrt{2}C_{14} & \sqrt{2}C_{24} & \sqrt{2}C_{34} & 2C_{44} & 2C_{45} & 2C_{46}\\ \sqrt{2}C_{15} & \sqrt{2}C_{25} & \sqrt{2}C_{35} & 2C_{45} & 2C_{55} & 2C_{56}\\ \sqrt{2}C_{16} & \sqrt{2}C_{26} & \sqrt{2}C_{36} & 2C_{46} & 2C_{56} & 2C_{66}\\ \end{bmatrix}\end{split}\]

References

[Helbig]

Helbig, K. (2013). What Kelvin might have written about Elasticity. Geophysical Prospecting, 61(1), 1-20. doi: 10.1111/j.1365-2478.2011.01049.x

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=200, C12=40, C44=20)
>>> print(C)
Stiffness tensor (in Voigt mapping):
[[200.  40.  40.   0.   0.   0.]
 [ 40. 200.  40.   0.   0.   0.]
 [ 40.  40. 200.   0.   0.   0.]
 [  0.   0.   0.  20.   0.   0.]
 [  0.   0.   0.   0.  20.   0.]
 [  0.   0.   0.   0.   0.  20.]]
>>> C.to_Kelvin()
array([[200.,  40.,  40.,   0.,   0.,   0.],
       [ 40., 200.,  40.,   0.,   0.,   0.],
       [ 40.,  40., 200.,   0.,   0.,   0.],
       [  0.,   0.,   0.,  40.,   0.,   0.],
       [  0.,   0.,   0.,   0.,  40.,   0.],
       [  0.,   0.,   0.,   0.,   0.,  40.]])

to_pymatgen()[source]

Convert the compliance tensor (from Elasticipy) to Python Materials Genomics (Pymatgen) format.

Returns:: Compliance tensor for pymatgen
Return type:: ComplianceTensor

transpose_array()[source]

Transpose the orientations of the tensor array

Returns:: The same tensor, but with transposed axes
Return type:: FourthOrderTensor

Examples

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> A = FourthOrderTensor.rand(shape=(3,4))
>>> A.transpose_array()
4th-order tensor array of shape (4, 3)

classmethod transverse_isotropic(*args, **kwargs)[source]

Create a stiffness tensor corresponding to the transversely isotropic symmetry with respect to Z axis, given the engineering constants.

Exactly two Poisson ratios must be provided (nu_xy or nu_yx, and nu_xz or nu_zx). See Notes for details.

Parameters:

Ex (float) – Young modulus along the x axis
Ez (float) – Young modulus along the y axis
Gxz (float) – Shear modulus in the x-z plane
nu_xy (float, optional) – Poisson ratio along x and y. Either nu_xy or nu_yx must be provided, not both.
nu_yx (float, optional) – Poisson ratio along x and y. Either nu_xy or nu_yx must be provided, not both.
nu_xz (float, optional) – Poisson ratio along x and z. Either nu_xz or nu_zx must be provided, not both.
nu_zx (float, optional) – Poisson ratio along x and z. Either nu_xz or nu_zx must be provided, not both.
kwargs (dict) – Keyword arguments to pass to the StiffnessTensor constructor

Return type:

StiffnessTensor

See also

orthotropic: create a stiffness tensor for orthotropic symmetry

Notes

If the material undergoes tensile strain \(\varepsilon_{ii}\) along the i-th direction, the Poisson ratios are defined as:

\[\nu_{ij}=-\frac{\partial \varepsilon_{jj}}{\partial \varepsilon_{ii}}\]

where \(\varepsilon_{jj}\) denotes the (compressive) longitudinal strain along the j-th direction. If \(E_x\) and \(E_y\) are the Young moduli along x and y, we have:

\[\frac{\nu_{xy}}{E_x} = \frac{\nu_{yx}}{E_y}\]

classmethod triclinic(S11=0.0, S12=0.0, S13=0.0, S14=0.0, S15=0.0, S16=0.0, S22=0.0, S23=0.0, C24=0.0, S25=0.0, S26=0.0, S33=0.0, C34=0.0, S35=0.0, S36=0.0, S44=0.0, S45=0.0, S46=0.0, S55=0.0, C56=0.0, S66=0.0, phase_name=None)[source]

Parameters:

S11 (float) – Components of the tensor
S12 (float) – Components of the tensor
S13 (float) – Components of the tensor
S14 (float) – Components of the tensor
S15 (float) – Components of the tensor
S16 (float) – Components of the tensor
S22 (float) – Components of the tensor
S23 (float) – Components of the tensor
C24 (float) – Components of the tensor
S25 (float) – Components of the tensor
S26 (float) – Components of the tensor
S33 (float) – Components of the tensor
C34 (float) – Components of the tensor
S35 (float) – Components of the tensor
S36 (float) – Components of the tensor
S44 (float) – Components of the tensor
S45 (float) – Components of the tensor
S46 (float) – Components of the tensor
S55 (float) – Components of the tensor
C56 (float) – Components of the tensor
S66 (float) – Components of the tensor
phase_name (str, optional) – Name to display

Return type:

FourthOrderTensor

See also

monoclinic: create a tensor from monoclinic symmetry
orthorhombic: create a tensor from orthorhombic symmetry

classmethod trigonal(*, S11=0.0, S12=0.0, S13=0.0, S14=0.0, S33=0.0, S44=0.0, S15=0.0, phase_name=None)[source]

Create a fourth-order tensor from trigonal symmetry.

Parameters:

S11 (float) – Components of the tensor, using the Voigt notation
S12 (float) – Components of the tensor, using the Voigt notation
S13 (float) – Components of the tensor, using the Voigt notation
S14 (float) – Components of the tensor, using the Voigt notation
S33 (float) – Components of the tensor, using the Voigt notation
S44 (float) – Components of the tensor, using the Voigt notation
S15 (float, optional) – S15 component of the tensor, only used for point groups 3 and -3.
phase_name (str, optional) – Phase name to display

Return type:

FourthOrderTensor

See also

tetragonal: create a tensor from tetragonal symmetry
orthorhombic: create a tensor from orthorhombic symmetry

property universal_anisotropy[source]

Compute the universal anisotropy factor.

It is actually an alias for inv().universal_anisotropy.

Returns:: Universal anisotropy factor
Return type:: float

wave_velocity(rho)[source]

Compute the wave velocities, given the mass density.

Parameters:: rho (float) – mass density. Its unit must be consistent with that of the stiffness tensor. See notes for hints.

See also

Christoffel_tensor: Computes the Christoffel tensor along a given direction

Returns:

c_p (SphericalFunction) – Velocity of the primary (compressive) wave
c_s1 (SphericalFunction) – Velocity of the fast secondary (shear) wave
c_s2 (SphericalFunction) – Velocity of the slow secondary (shear) wave

Notes

The estimation of the wave velocities is made by finding the eigenvalues of the Christoffel tensor [Jaeken].

One should double-check the units. The table below provides hints about the unit you get, depending on the units you use for stiffness and the mass density:

Stiffness	Mass density	Velocities	Notes
Pa (N/m²)	kg/m³	m/s	SI units
GPa (10⁹ Pa)	kg/dm³	km/s	Conversion factor
GPa (10³ N/mm²)	kg/mm³	m/s	Consistent units
MPa (10⁶ Pa)	kg/m³	km/s	Conversion factor
MPa (10³ N/mm²)	g/mm³	m/s	Consistent units

References

[Jaeken]

J. W. Jaeken, S. Cottenier, Solving the Christoffel equation: Phase and group velocities, Computer Physics Communications (207), 2016, https://doi.org/10.1016/j.cpc.2016.06.014.

Examples

We will investigate the wave velocities in cubic Cu. First, define the stiffness tensor:

>>> C = StiffnessTensor.cubic(C11=186e3, C12=134e3, C44=77e3) # Stiffness in MPa
>>> rho = 8960 # mass density in kg/m³
>>> c_p, c_s1, c_s2 = C.wave_velocity(rho)

The range of velocities of the (fast) primary wave range can be printed with:

>>> print(c_p)
Spherical function
Min=4.556196722204669, Max=5.324304112812698

As the stiffness is given in MPa and the mass desity is given in kg/m³, the velocities are returned in km/s (see Notes). If one wants to know the direction for min/max values:

>>> val_min, u_min = c_p.min()
>>> print(u_min)
[[0. 0. 1.]]

>>> val_max, u_max = c_p.max()
>>> print(u_max)
[[0.57735022 0.57735033 0.57735026]]

The secondary wave velocities are given by cs_1 and cs_2:

>>> c_s1
Spherical function
Min=2.1906864565414192, Max=2.9315098498896446

>>> c_s2
Spherical function
Min=1.7034628596749235, Max=2.9315098498896437

classmethod weighted_average(*args)[source]

Compute the weighted average of a list of stiffness tensors, with respect to a given method (Voigt, Reuss or Hill).

Parameters:

Cs (list of StiffnessTensor or list of ComplianceTensor or tuple of StiffnessTensor or tuple of ComplianceTensor) – Series of tensors to compute the average from
volume_fractions (iterable of floats) – Volume fractions of each phase
method (str, {'Voigt', 'Reuss', 'Hill'}) – Method to use. It can be ‘Voigt’, ‘Reuss’, or ‘Hill’.

Returns:

Average tensor

Return type:

StiffnessTensor

classmethod zeros(shape=(), **kwargs)[source]

Create a fourth-order tensor populated with zeros

Parameters:

shape (int or tuple, optional) – Shape of the tensor to create
kwargs – Keyword arguments passed to the FourthOrderTensor constructor

Return type:

FourthOrderTensor

Examples

The single-valued null 4th order tensor is just:

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> FourthOrderTensor.zeros()
4th-order tensor (in Kelvin mapping):
[[0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]]

One can also create an array of such tensors:

>>> zeros_tensor = FourthOrderTensor.zeros(shape=3)

and check that it populated with zeros:

>>> zeros_tensor == 0.
array([ True,  True,  True])

class Elasticipy.tensors.elasticity.StiffnessTensor(M, check_positive_definite=True, phase_name=None, mapping=<Elasticipy.tensors.mapping.VoigtMapping object>, **kwargs)[source]

Bases: SymmetricFourthOrderTensor

Class for manipulating fourth-order stiffness tensors.

Construct a stiffness tensor or an array of stiffness tensors.

The stiffness tensor can be constructed from a (6,6) matrix or slices of (6,6) matrices. These matrices must be symmetric. An error is thrown if this matrix in not definite positive (except if check_positive_definite==False, see below). The input argument can also be the full tensor (array of shape (…,3,3,3,3)).

Parameters:

M (list or np.ndarray) – (…,6,6) matrix corresponding to the stiffness tensor, written using the Voigt notation, or array of shape (…,3,3,3,3).
phase_name (str, optional) – Phase name to display
check_positive_definite (bool, optional) – Whether to check if the input matrix is positive definite or not. True by default.
check_symmetry (bool, optional) – Whether to check or not that the input matrix is symmetric.
force_symmetry (bool, optional) – If true, the major symmetry of the tensor is forced
mapping (str or MappingConvention) – mapping convention to use. Default is VoigtMapping.

Notes

The units used when building the stiffness tensor are up to the user (GPa, MPa, psi etc.). Therefor, the results you will get when performing operations (Young’s modulus, “product” with strain tensor etc.) will be consistent with these units. For instance, if the stiffness tensor is defined in GPa, the computed stress will be given in GPa as well.

Examples

Create a stiffness tensor for cubic symmetry:

>>> matrix = [[200, 40,  40,  0,  0,  0 ],
...           [40,  200, 40,  0,  0,  0 ],
...           [40,  40,  200, 0,  0,  0 ],
...           [0,   0,   0,   20, 0,  0 ],
...           [0,   0,   0,   0,  20, 0 ],
...           [0,   0,   0,   0,   0, 20]]
>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor(matrix)
>>> print(C)
Stiffness tensor (in Voigt mapping):
[[200.  40.  40.   0.   0.   0.]
 [ 40. 200.  40.   0.   0.   0.]
 [ 40.  40. 200.   0.   0.   0.]
 [  0.   0.   0.  20.   0.   0.]
 [  0.   0.   0.   0.  20.   0.]
 [  0.   0.   0.   0.   0.  20.]]

Create a stiffness tensor from full (3,3,3,3) array:

>>> C_full = C.full_tensor # (3,3,3,3) numpy array
>>> print((type(C_full), C_full.shape))
(<class 'numpy.ndarray'>, (3, 3, 3, 3))

>>> StiffnessTensor(C_full)
Stiffness tensor (in Voigt mapping):
[[200.  40.  40.   0.   0.   0.]
 [ 40. 200.  40.   0.   0.   0.]
 [ 40.  40. 200.   0.   0.   0.]
 [  0.   0.   0.  20.   0.   0.]
 [  0.   0.   0.   0.  20.   0.]
 [  0.   0.   0.   0.   0.  20.]]

Create an array of stiffness tensors:

First, we create two slices of (6,6) matrices, corresponding to two stiffness values:

>>> slices = [[[200, 40,  40,  0,  0,  0 ],
...            [40,  200, 40,  0,  0,  0 ],
...            [40,  40,  200, 0,  0,  0 ],
...            [0,   0,   0,   20, 0,  0 ],
...            [0,   0,   0,   0,  20, 0 ],
...            [0,   0,   0,   0,   0, 20]],
...           [[250, 80,  80,  0,  0,  0 ],
...            [80,  250, 80,  0,  0,  0 ],
...            [80,  80,  250, 0,  0,  0 ],
...            [0,   0,   0,   40, 0,  0 ],
...            [0,   0,   0,   0,  40, 0 ],
...            [0,   0,   0,   0,   0, 40]]]

Then, one can create an array of stiffness tensors:

>>> C_array=StiffnessTensor(slices)
>>> print(C_array)
Stiffness tensor (in Voigt mapping):
[[[200.  40.  40.   0.   0.   0.]
  [ 40. 200.  40.   0.   0.   0.]
  [ 40.  40. 200.   0.   0.   0.]
  [  0.   0.   0.  20.   0.   0.]
  [  0.   0.   0.   0.  20.   0.]
  [  0.   0.   0.   0.   0.  20.]]

 [[250.  80.  80.   0.   0.   0.]
  [ 80. 250.  80.   0.   0.   0.]
  [ 80.  80. 250.   0.   0.   0.]
  [  0.   0.   0.  40.   0.   0.]
  [  0.   0.   0.   0.  40.   0.]
  [  0.   0.   0.   0.   0.  40.]]]

This array can be subindexed. E.g.:

>>> C_array[0]
Stiffness tensor (in Voigt mapping):
[[200.  40.  40.   0.   0.   0.]
 [ 40. 200.  40.   0.   0.   0.]
 [ 40.  40. 200.   0.   0.   0.]
 [  0.   0.   0.  20.   0.   0.]
 [  0.   0.   0.   0.  20.   0.]
 [  0.   0.   0.   0.   0.  20.]]

Christoffel_tensor(u)[source]

Create the Christoffel tensor along a given direction, or set or directions.

Parameters:: u (list or np.ndarray) – 3D direction(s) to compute the Christoffel tensor along with
Returns:: Christoffel tensor(s). if u is a list of directions, Gamma[i] is the Christoffel tensor for direction u[i].
Return type:: SymmetricSecondOrderTensor

See also

wave_velocity: computes the p- and s-wave velocities.

Notes

For a given stiffness tensor C and a given unit vector u, the Christoffel tensor is defined as [Jaeken] :

\[M_{ij} = C_{iklj}.u_k.u_l\]

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77) # Cubic Cu
>>> C.Christoffel_tensor([0, 0, 1])
Symmetric second-order tensor
[[ 77.   0.   0.]
 [  0.  77.   0.]
 [  0.   0. 186.]]

Hill_average(axis=None)[source]

Compute the Voigt-Reuss-Hill average (mean of Voigt and Reuss averages for stiffness tensors).

If the object is a single tensor, the returned tensor corresponds to the average over an infinite set of random rotations, resulting in an isotropic behavior. Otherwise, the mean is computed over the given axis.

Parameters:: axis (int, optional) – If provided, axis to compute the average along with. If none, the average is computed on the flattened array
Returns:: Voigt-Reuss-Hill average of tensor
Return type:: StiffnessTensor

See also

Voigt_average: compute the Voigt average
Reuss_average: compute the Reuss average
average: generic function for calling either the Voigt, Reuss or Hill average

Examples

Let us consider the stiffness of pure copper:

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)

If we assume that an aggregate is composed of an infinite set grains whose orientations are uniformly distributed, the Voigt average is:

>>> C.Hill_average()
Stiffness tensor (in Voigt mapping):
[[217.83103448 118.08448276 118.08448276   0.           0.
    0.        ]
 [118.08448276 217.83103448 118.08448276   0.           0.
    0.        ]
 [118.08448276 118.08448276 217.83103448   0.           0.
    0.        ]
 [  0.           0.           0.          49.87327586   0.
    0.        ]
 [  0.           0.           0.           0.          49.87327586
    0.        ]
 [  0.           0.           0.           0.           0.
   49.87327586]]

Now, we consider an aggregate of grains whose orientations follow a pure fibre texture along the X axis. First, generate a (large) random set of rotation corresponding to this texture:

>>> from scipy.spatial.transform import Rotation
>>> import numpy as np
>>> g = Rotation.from_euler('X', np.linspace(0, 90, 1000), degrees=True)    # 1000 rotations around X

Now apply rotations to compute the array of stiffness tensors:

>>> C_rotated = C * g
>>> C_rotated
Stiffness tensor array of shape (1000,)

Finally, one can check that the Voigt average, computed from the rotated stiffness tensors is:

>>> C_rotated.Hill_average()
Stiffness tensor (in Voigt mapping):
[[ 1.86000000e+02  1.34000000e+02  1.34000000e+02 -8.62027175e-16
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  2.05164524e+02  1.14835476e+02 -5.48889226e-16
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  1.14835476e+02  2.05164524e+02 -1.44776034e-15
   0.00000000e+00  0.00000000e+00]
 [-8.09683464e-17  8.46209730e-17 -3.98365128e-16  4.52092721e+01
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.70000000e+01 -7.54674101e-17]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.54674101e-17  7.70000000e+01]]

property Poisson_ratio[source]

Directional Poisson’s ratio

Returns:: Poisson’s ratio
Return type:: HyperSphericalFunction

Notes

If the material undergoes tensile strain \(\varepsilon_{ii}\) along the i-th direction, the Poisson ratios are defined as:

\[\nu_{ij}=-\frac{\partial \varepsilon_{jj}}{\partial \varepsilon_{ii}}\]

where \(\varepsilon_{jj}\) denotes the (compressive) longitudinal strain along the j-th direction.

Reuss_average(axis=None)[source]

Compute the Reuss average (from the mean of compliance tensors).

If the object is a single tensor, the returned tensor corresponds to the average over an infinite set of random rotations, resulting in an isotropic behavior. Otherwise, the mean is computed over the given axis.

Parameters:: axis (int, optional) – If provided, axis to compute the average along with. If none, the average is computed on the flattened array
Returns:: Reuss average of stiffness tensor
Return type:: StiffnessTensor

See also

Voigt_average: compute the Voigt average
Hill_average: compute the Voigt-Reuss-Hill average
average: generic function for calling either the Voigt, Reuss or Hill average

Examples

Let us consider the stiffness of pure copper:

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)

If we assume that an aggregate is composed of an infinite set grains whose orientations are uniformly distributed, the Voigt average is:

>>> C.Reuss_average()
Stiffness tensor (in Voigt mapping):
[[208.86206897 122.56896552 122.56896552   0.           0.
    0.        ]
 [122.56896552 208.86206897 122.56896552   0.           0.
    0.        ]
 [122.56896552 122.56896552 208.86206897   0.           0.
    0.        ]
 [  0.           0.           0.          43.14655172   0.
    0.        ]
 [  0.           0.           0.           0.          43.14655172
    0.        ]
 [  0.           0.           0.           0.           0.
   43.14655172]]

Now, we consider an aggregate of grains whose orientations follow a pure fibre texture along the X axis. First, generate a (large) random set of rotation corresponding to this texture:

>>> from scipy.spatial.transform import Rotation
>>> import numpy as np
>>> g = Rotation.from_euler('X', np.linspace(0, 90, 1000), degrees=True)    # 1000 rotations around X

Now apply rotations to compute the array of stiffness tensors:

>>> C_rotated = C * g
>>> C_rotated
Stiffness tensor array of shape (1000,)

Finally, one can check that the Voigt average, computed from the rotated stiffness tensors is:

>>> C_rotated.Reuss_average()
Stiffness tensor (in Voigt mapping):
[[ 1.86000000e+02  1.34000000e+02  1.34000000e+02 -1.68008952e-15
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  1.98854548e+02  1.21145452e+02 -1.09777845e-15
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  1.21145452e+02  1.98854548e+02 -2.89552069e-15
   0.00000000e+00  0.00000000e+00]
 [-1.61936693e-16  1.69241946e-16 -7.96730255e-16  3.88930441e+01
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.70000000e+01 -7.54674101e-17]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.54674101e-17  7.70000000e+01]]

Voigt_average(axis=None)[source]

Compute the Voigt average (from the mean of stiffness tensors).

If the object is a single tensor, the returned tensor corresponds to the average over an infinite set of random rotations, resulting in an isotropic behavior. Otherwise, the mean is computed over the given axis.

Parameters:: axis (int, optional) – If provided, the average is computed along this axis. Otherwise, the mean is computed on the flattened array.
Returns:: Voigt average of stiffness tensor
Return type:: StiffnessTensor

See also

Reuss_average: compute the Reuss average
Hill_average: compute the Voigt-Reuss-Hill average
average: generic function for calling either the Voigt, Reuss or Hill average

Examples

Let us consider the stiffness of pure copper:

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)

If we assume that an aggregate is composed of an infinite set grains whose orientations are uniformly distributed, the Voigt average is:

>>> C.Voigt_average()
Stiffness tensor (in Voigt mapping):
[[226.8 113.6 113.6   0.    0.    0. ]
 [113.6 226.8 113.6   0.    0.    0. ]
 [113.6 113.6 226.8   0.    0.    0. ]
 [  0.    0.    0.   56.6   0.    0. ]
 [  0.    0.    0.    0.   56.6   0. ]
 [  0.    0.    0.    0.    0.   56.6]]

Now, we consider an aggregate of grains whose orientations follow a pure fibre texture along the X axis. First, generate a (large) random set of rotation corresponding to this texture:

>>> from scipy.spatial.transform import Rotation
>>> import numpy as np
>>> g = Rotation.from_euler('X', np.linspace(0, 90, 1000), degrees=True)    # 1000 rotations around X

Now apply rotations to compute the array of stiffness tensors:

>>> C_rotated = C * g
>>> C_rotated
Stiffness tensor array of shape (1000,)

Finally, one can check that the Voigt average, computed from the rotated stiffness tensors is:

>>> C_rotated.Voigt_average()
Stiffness tensor (in Voigt mapping):
[[ 1.86000000e+02  1.34000000e+02  1.34000000e+02 -4.39648318e-17
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  2.11474500e+02  1.08525500e+02  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 1.34000000e+02  1.08525500e+02  2.11474500e+02  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  5.15255000e+01
   0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.70000000e+01 -7.54674101e-17]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   7.54674101e-17  7.70000000e+01]]

property Young_modulus[source]

Directional Young’s modulus

Returns:: Young’s modulus
Return type:: SphericalFunction

Zener_ratio(tol=0.0001)[source]

Compute the Zener ratio (Z). Only valid for cubic symmetry.

This function first checks that the tensor has cubic symmetry within a given tolerance. If not, an error is raised.

Parameters:: tol (float, optional) – Tolerance to consider that the material has cubic symmetry
Returns:: Zener ratio (NaN is the symmetry is not cubic)
Return type:: float

Notes

If the tensor is written in canonical base with Voigt mapping, the Zener ratio is defined as:

\[Z=\frac{ 2C_{44} }{C_{11} - C_{12}}\]

The present implementation takes advantage of eigenstiffness to compute the Zener ratio in any base, i.e. even if the tensor is not given in canonical base (e.g. if rotated).

See also

universal_anisotropy: compute the universal anisotropy factor
eig_stiffnesses: eigenstiffnesses of the tensor
is_cubic: check whether the tensor has cubic symmetry or not

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=200, C12=40, C44=20)
>>> C.Zener_ratio()
0.25

which obvisouly corresponds to 2.C44/(C11-C12).

Now, rotate the tensor and see how it looks like:

>>> from scipy.spatial.transform import Rotation
>>> g = Rotation.from_euler('Z', 30, degrees=True)
>>> C_rot = C*g
>>> C_rot
Stiffness tensor (in Voigt mapping):
[[155.          85.          40.           0.           0.
   25.98076211]
 [ 85.         155.          40.           0.           0.
  -25.98076211]
 [ 40.          40.         200.           0.           0.
    0.        ]
 [  0.           0.           0.          20.           0.
    0.        ]
 [  0.           0.           0.           0.          20.
    0.        ]
 [ 25.98076211 -25.98076211   0.           0.           0.
   65.        ]]

Still, we have >>> C_rot.Zener_ratio() 0.24999999999999983

average(method, axis=None)[source]

Compute either the Voigt, Reuss, or Hill average of the stiffness tensor.

This function is just a shortcut for Voigt_average(), Reuss_average(), or Hill_average() and Hill_average().

Parameters:

axis (int, optional) – If provided, axis to compute the average along with. If none, the average is computed on the flattened array
method (str {'Voigt', 'Reuss', 'Hill'}) – Method to use to compute the average.

Return type:

StiffnessTensor

See also

Voigt_average: compute the Voigt average
Reuss_average: compute the Reuss average
Hill_average: compute the Voigt-Reuss-Hill average

property bulk_modulus[source]

Compute the bulk modulus of the material

Returns:: Bulk modulus
Return type:: float or numpy.ndarray

See also

linear_compressibility: directional linear compressibility

copy()[source]: Create a copy of the tensor

classmethod cubic(*, C11=0.0, C12=0.0, C44=0.0, phase_name=None)[source]

Create a fourth-order tensor from cubic symmetry.

Parameters:

C11 (float)
C12 (float)
C44 (float)
phase_name (str, optional) – Phase name to display

Return type:

StiffnessTensor

See also

hexagonal: create a tensor from hexagonal symmetry
orthorhombic: create a tensor from orthorhombic symmetry

ddot(other, mode='pair')[source]

Perform tensor product contracted twice (“:”) between two fourth-order tensors

Parameters:

other (FourthOrderTensor or SecondOrderTensor) – Right-hand side of “:” symbol
mode (str, optional) – If mode==”pair”, the tensors must be broadcastable, and the tensor product are performed on the last axes. If mode==”cross”, all cross-combinations are considered.

Returns:

Result from double-contraction

Return type:

FourthOrderTensor

Examples

First, let consider two random arrays of Fourth-order tensors:

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> T1 = FourthOrderTensor.rand(shape=(2,3))
>>> T2 = FourthOrderTensor.rand(shape=3)
>>> T1T2_pair = T1.ddot(T2)
>>> T1T2_pair
4th-order tensor array of shape (2, 3)

whereas:

>>> T1T2_cross = T1.ddot(T2, mode='cross')
>>> T1T2_cross
4th-order tensor array of shape (2, 3, 3)

The command above is equivalent (but way faster) to:

>>> T1T2_cross_loop = FourthOrderTensor.zeros(shape=(2,3,3))
>>> for i in range(2):
...     for j in range(3):
...         for k in range(3):
...             T1T2_cross_loop[i,j,k] = T1[i,j].ddot(T2[k])

One can check that the results are consistent with:

>>> T1T2_cross_loop == T1T2_cross
array([[[ True,  True,  True],
        [ True,  True,  True],
        [ True,  True,  True]],

       [[ True,  True,  True],
        [ True,  True,  True],
        [ True,  True,  True]]])

deviatoric_part()[source]

Return the deviatoric part of the tensor

Returns:: Deviatoric part of the tensor
Return type:: FourthOrderTensor

See also

identity_tensor: return the identity tensor
spherical_part: return the spherical part of the tensor

eig()[source]

Compute the eigenstiffnesses and the eigenstrains.

Solve the eigenvalue problem from the Kelvin matrix of the stiffness tensor (see Notes).

Returns:

numpy.ndarray – Array of 6 eigenstiffnesses (eigenvalues of the stiffness matrix)
numpy.ndarray – (6,6) array of eigenstrains (eigenvectors of the stiffness matrix)

See also

to_Kelvin: returns the stiffness components as a (6,6) matrix, according to the Kelvin mapping convention.
eig_stiffnesses: returns the eigenstiffnesses only
eig_strains: returns the eigenstrains only

Notes

The definition for eigenstiffnesses and the eigenstrains are introduced in [Helbig].

property eig_compliances[source]

Compute the eigencompliances from the Kelvin’s matrix of stiffness

Returns:: Inverses of the 6 eigenvalues of the Kelvin’s stiffness matrix, in descending order
Return type:: numpy.ndarray

See also

eig_stiffnesses: compute the eigenstiffnesses from the Kelvin’s matrix of stiffness

property eig_stiffnesses[source]

Compute the eigenstiffnesses given by the Kelvin’s matrix for stiffness.

Returns:: 6 eigenvalues of the Kelvin’s stiffness matrix, in ascending order
Return type:: numpy.ndarray

See also

eig: returns the eigenstiffnesses and the eigenstrains
eig_strains: returns the eigenstrains only
eig_stiffnesses_multiplicity: returns the unique values of eigenstiffnesses with multiplicity

Notes

The eigenstiffnesses are defined in [Helbig].

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=200, C12=40, C44=20)
>>> C.eig_stiffnesses
array([ 40.,  40.,  40., 160., 160., 280.])

These values actually correspond to 2*C44 (with multiplicity = 3), C11-C12 (with multiplicity = 2) and C11+2C12 (no multiplicity); see [Helbig].

eig_stiffnesses_multiplicity(tol=0.0001)[source]

Compute the eigenstiffnesses, then returns the multiplicity of each eigenstiffness.

Given an absolute tolerance, duplicates in eigenstiffnesses are considered to compute the multiplicity of each value.

Parameters:

tol (float, optional) – Absolute tolerance to assume that two distinct eigenstiffnesses are the same

Returns:

numpy.ndarray – Unique values of eigenstiffnesses, sorted by increasing multiplicity
numpy.ndarray – Multiplicity of each unique eigenstiffness, sorted by ascending value

See also

eig_stiffnesses: compute the eigenstiffnesses

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)
>>> C.eig_stiffnesses
array([ 52.,  52., 154., 154., 154., 454.])
>>> C.eig_stiffnesses_multiplicity()
(array([ 52., 154., 454.]), array([2, 3, 1]))

property eig_strains[source]

Compute the eigenstrains from the Kelvin’s matrix for stiffness

Returns:: (6,6) matrix of eigenstrains, sorted by ascending order of eigenstiffnesses.
Return type:: numpy.ndarray

See also

eig: returns both the eigenvalues and the eigenvectors of the Kelvin matrix

classmethod eye(shape=(), **kwargs)[source]

Create a 4th-order identity tensor.

See notes for definition.

Parameters:

shape (int or tuple, optional) – Shape of the tensor to create
mapping (Kelvin mapping, optional) – Mapping convention to use. Must be either Kelvin or Voigt.

Returns:

Identity tensor

Return type:

FourthOrderTensor or SymmetricFourthOrderTensor

Notes

The Fourth-order identity tensor is defined as:

\[I_{ijkl} = \frac12\left( \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}\right)\]

Examples

Create a (single) identity tensor:

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> I = FourthOrderTensor.eye()
>>> print(I)
4th-order tensor (in Kelvin mapping):
[[1. 0. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 1.]]

Alternatively, one can use another mapping convention, e.g. Voigt:

>>> from Elasticipy.tensors.mapping import VoigtMapping
>>> Iv = FourthOrderTensor.eye(mapping=VoigtMapping())
>>> print(Iv)
4th-order tensor (in Voigt mapping):
[[1.  0.  0.  0.  0.  0. ]
 [0.  1.  0.  0.  0.  0. ]
 [0.  0.  1.  0.  0.  0. ]
 [0.  0.  0.  0.5 0.  0. ]
 [0.  0.  0.  0.  0.5 0. ]
 [0.  0.  0.  0.  0.  0.5]]

Still, we have:

>>> I == Iv
True

as they correspond to the same tensor, but expressed as a matrix with different mapping conventions. Indeed, one can check that:

>>> import numpy as np
>>> np.array_equal(I.full_tensor, Iv.full_tensor)
True

flatten()[source]

Flatten the tensor

If the tensor array is of shape (m,n,o…,r), the flattened array will be of shape (m*n*o*…*r,).

Returns:: Flattened tensor
Return type:: SymmetricFourthOrderTensor

Examples

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> T = FourthOrderTensor.rand(shape=(5,6))
>>> T
4th-order tensor array of shape (5, 6)
>>> T.flatten()
4th-order tensor array of shape (30,)

classmethod fromCrystalSymmetry(symmetry='Triclinic', point_group=None, diad='y', phase_name=None, prefix=None, **kwargs)[source]

Create a fourth-order tensor from limited number of components, taking advantage of crystallographic symmetries

Parameters:

symmetry (str, default Triclinic) – Name of the crystallographic symmetry
point_group (str) – Point group of the considered crystal. Only used (and mandatory) for tetragonal and trigonal symmetries.
diad (str {'x', 'y'}, default 'x') – Alignment convention. Sets whether x||a or y||b. Only used for monoclinic symmetry.
phase_name (str, default None) – Name to use when printing the tensor
prefix (str, default None) – Define the prefix to use when providing the components. By default, it is ‘C’ for stiffness tensors, ‘S’ for compliance.
kwargs – Keywords describing all the necessary components, depending on the crystal’s symmetry and the type of tensor. For Stiffness, they should be named as ‘Cij’ (e.g. C11=…, C12=…). For Comliance, they should be named as ‘Sij’ (e.g. S11=…, S12=…). See examples below. The behaviour can be overriten with the prefix option (see above)

Return type:

FourthOrderTensor

See also

isotropic: creates an isotropic stiffness tensor from two paremeters (e.g. E and v).

Notes

The relationships between the tensor’s components depend on the crystallogrpahic symmetry [Nye].

References

[Nye]

Nye, J. F. Physical Properties of Crystals. London: Oxford University Press, 1959.

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor

>>> StiffnessTensor.fromCrystalSymmetry(symmetry='monoclinic', diad='y', phase_name='TiNi',
...                                     C11=231, C12=127, C13=104,
...                                     C22=240, C23=131, C33=175,
...                                     C44=81, C55=11, C66=85,
...                                     C15=-18, C25=1, C35=-3, C46=3)
Stiffness tensor (in Voigt mapping):
[[231. 127. 104.   0. -18.   0.]
 [127. 240. 131.   0.   1.   0.]
 [104. 131. 175.   0.  -3.   0.]
 [  0.   0.   0.  81.   0.   3.]
 [-18.   1.  -3.   0.  11.   0.]
 [  0.   0.   0.   3.   0.  85.]]
Phase: TiNi

>>> from Elasticipy.tensors.elasticity import ComplianceTensor

>>> ComplianceTensor.fromCrystalSymmetry(symmetry='monoclinic', diad='y', phase_name='TiNi',
...                                      S11=8, S12=-3, S13=-2,
...                                      S22=8, S23=-5, S33=10,
...                                      S44=12, S55=116, S66=12,
...                                      S15=14, S25=-8, S35=0, S46=0)
Compliance tensor (in Voigt mapping):
[[  8.  -3.  -2.   0.  14.   0.]
 [ -3.   8.  -5.   0.  -8.   0.]
 [ -2.  -5.  10.   0.   0.   0.]
 [  0.   0.   0.  12.   0.   0.]
 [ 14.  -8.   0.   0. 116.   0.]
 [  0.   0.   0.   0.   0.  12.]]
Phase: TiNi

classmethod from_Kelvin(matrix, **kwargs)[source]

Create a tensor from the (6,6) matrix following the Kelvin(-Mandel) mapping convention

Parameters:

matrix (list or numpy.ndarray) – (6,6) matrix of components
kwargs – keyword arguments passed to the constructor

Return type:

StiffnessTensor

See also

to_Kelvin: return the components as a (6,6) matrix following the Kelvin convention

classmethod from_MP(ids, api_key=None)[source]

Import stiffness tensor(s) from the Materials Project API, given their material ids.

You need to register to https://materialsproject.org first to get an API key. This key can be explicitly passed as an argument (see below), or provided as an environment variable named MP_API_KEY.

Parameters:

ids (str or list of str) – ID(s) of the material to import (e.g. “mp-1048”)
api_key (str, optional) – API key to the Materials Project API. If not provided, it should be available as the API_KEY environment variable.

Returns:

If one of the requested material ids was not found, the corresponding value in the list will be None.

Return type:

list of StiffnessTensor

classmethod from_txt_file(filename)[source]

Load the tensor from a text file.

The two first lines can have data about phase name and symmetry, but this is not mandatory.

Parameters:: filename (str) – Filename to load the tensor from.
Returns:: The reconstructed tensor read from the file.
Return type:: SymmetricFourthOrderTensor

See also

save_to_txt: create a tensor from text file

property full_tensor[source]

Returns the full (unvoigted) tensor as a (3, 3, 3, 3) or (…, 3, 3, 3, 3) array

Returns:: Full tensor (4-index notation)
Return type:: np.ndarray

Examples

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> I = FourthOrderTensor.eye() # 4th order identity tensor
>>> print(I)
4th-order tensor (in Kelvin mapping):
[[1. 0. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 1.]]

>>> I_full = I.full_tensor
>>> type(I_full)
<class 'numpy.ndarray'>
>>> I_full.shape
(3, 3, 3, 3)

When working on tensor arrays, the shape of the resulting numpy array will change accordlingly. E.g.:

>>> I_array = FourthOrderTensor.eye(shape=(5,6)) # Array of 4th order identity tensor
>>> I_array.full_tensor.shape
(5, 6, 3, 3, 3, 3)

classmethod hexagonal(*, C11=0.0, C12=0.0, C13=0.0, C33=0.0, C44=0.0, phase_name=None)[source]

Create a fourth-order tensor from hexagonal symmetry.

Parameters:

C11 (float) – Components of the tensor, using the Voigt notation
C12 (float) – Components of the tensor, using the Voigt notation
C13 (float) – Components of the tensor, using the Voigt notation
C33 (float) – Components of the tensor, using the Voigt notation
C44 (float) – Components of the tensor, using the Voigt notation
phase_name (str, optional) – Phase name to display

Return type:

FourthOrderTensor

See also

transverse_isotropic: creates a transverse-isotropic tensor from engineering parameters
cubic: create a tensor from cubic symmetry
tetragonal: create a tensor from tetragonal symmetry

classmethod identity(**kwargs)[source]

Construct the Fourth-order identity tensor.

This is actually an alias for eye().

Parameters:: kwargs – Keyword arguments passed to the Fourth-order tensor constructor.
Return type:: FourthOrderTensor

See also

eye: Fourth-order identity tensor

classmethod identity_deviatoric_part(**kwargs)[source]

Return the deviatoric part of the identity tensor.

See notes for the mathematical definition.

Parameters:: kwargs – keyword arguments passed to eye constructor
Return type:: FourthOrderTensor or SymmetricTensor

See also

identity_tensor: return the identity tensor
identity_spherical_part: return the spherical part of the identity tensor

Notes

The deviatoric part of the identity tensor is defined as:

\[K = I - J\]

where \(I\) and \(J\) denote the identity and the deviatoric part of the identity tensor, respectively.

Examples

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> K = FourthOrderTensor.identity_deviatoric_part()
>>> print(K)
4th-order tensor (in Kelvin mapping):
[[ 0.66666667 -0.33333333 -0.33333333  0.          0.          0.        ]
 [-0.33333333  0.66666667 -0.33333333  0.          0.          0.        ]
 [-0.33333333 -0.33333333  0.66666667  0.          0.          0.        ]
 [ 0.          0.          0.          1.          0.          0.        ]
 [ 0.          0.          0.          0.          1.          0.        ]
 [ 0.          0.          0.          0.          0.          1.        ]]

One can check that K has zero spherical part:

>>> print(K.spherical_part())
4th-order tensor (in Kelvin mapping):
[[2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]]

classmethod identity_spherical_part(shape=(), **kwargs)[source]

Return the spherical part of the identity tensor.

See Notes for mathematical definition.

Parameters:

shape (tuple of int, optional) – Shape of the tensor to create
kwargs – Keyword arguments passed to the Fourth-order tensor constructor.

Return type:

FourthOrderTensor

See also

identity_tensor: return the identity tensor
identity_deviatoric_part: return the deviatoric part of the identity tensor

Notes

The spherical part of the identity tensor is defined as:

\[J_{ijkl} = \frac13 \delta_{ij}\delta_{kl}\]

Examples

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> J = FourthOrderTensor.identity_spherical_part()
>>> print(J)
4th-order tensor (in Kelvin mapping):
[[0.33333333 0.33333333 0.33333333 0.         0.         0.        ]
 [0.33333333 0.33333333 0.33333333 0.         0.         0.        ]
 [0.33333333 0.33333333 0.33333333 0.         0.         0.        ]
 [0.         0.         0.         0.         0.         0.        ]
 [0.         0.         0.         0.         0.         0.        ]
 [0.         0.         0.         0.         0.         0.        ]]

On can check that J has zero deviatoric part:

>>> J.deviatoric_part()
4th-order tensor (in Kelvin mapping):
[[2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [2.77555756e-17 2.77555756e-17 2.77555756e-17 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00]]

infinite_random_average()[source]

Compute the average of the tensor, assuming that an infinite number of random orientations is applied.

Returns:: Average tensor or tensor array. The returned array will be of the same shape as the input object.
Return type:: FourthOrderTensor

inv()[source]

Compute the reciprocal compliance tensor

Returns:: Reciprocal tensor
Return type:: ComplianceTensor

Examples

Create a stiffness tensor for cubic symmetry:

>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)
>>> print(C)
Stiffness tensor (in Voigt mapping):
[[186. 134. 134.   0.   0.   0.]
 [134. 186. 134.   0.   0.   0.]
 [134. 134. 186.   0.   0.   0.]
 [  0.   0.   0.  77.   0.   0.]
 [  0.   0.   0.   0.  77.   0.]
 [  0.   0.   0.   0.   0.  77.]]

The reciprocal (compliance) tensor is given by:

>>> S = C.inv()
>>> print(S)
Compliance tensor (in Voigt mapping):
[[ 0.01355473 -0.00567604 -0.00567604  0.          0.          0.        ]
 [-0.00567604  0.01355473 -0.00567604  0.          0.          0.        ]
 [-0.00567604 -0.00567604  0.01355473  0.          0.          0.        ]
 [ 0.          0.          0.          0.01298701  0.          0.        ]
 [ 0.          0.          0.          0.          0.01298701  0.        ]
 [ 0.          0.          0.          0.          0.          0.01298701]]

This also works for tensor arrays. In this case, the returned value will be a compliance tensor array of the same shape as the tensor array. For instance:

>>> slices = [[[200, 40,  40,  0,  0,  0 ],
...            [40,  200, 40,  0,  0,  0 ],
...            [40,  40,  200, 0,  0,  0 ],
...            [0,   0,   0,   20, 0,  0 ],
...            [0,   0,   0,   0,  20, 0 ],
...            [0,   0,   0,   0,   0, 20]],
...           [[250, 80,  80,  0,  0,  0 ],
...            [80,  250, 80,  0,  0,  0 ],
...            [80,  80,  250, 0,  0,  0 ],
...            [0,   0,   0,   40, 0,  0 ],
...            [0,   0,   0,   0,  40, 0 ],
...            [0,   0,   0,   0,   0, 40]]]
>>> C = StiffnessTensor(slices)
>>> S = C.inv()
>>> print(S)
Compliance tensor (in Voigt mapping):
[[[ 0.00535714 -0.00089286 -0.00089286  0.          0.
    0.        ]
  [-0.00089286  0.00535714 -0.00089286  0.          0.
    0.        ]
  [-0.00089286 -0.00089286  0.00535714  0.          0.
    0.        ]
  [ 0.          0.          0.          0.05        0.
    0.        ]
  [ 0.          0.          0.          0.          0.05
    0.        ]
  [ 0.          0.          0.          0.          0.
    0.05      ]]

 [[ 0.00473458 -0.00114778 -0.00114778  0.          0.
    0.        ]
  [-0.00114778  0.00473458 -0.00114778  0.          0.
    0.        ]
  [-0.00114778 -0.00114778  0.00473458  0.          0.
    0.        ]
  [ 0.          0.          0.          0.025       0.
    0.        ]
  [ 0.          0.          0.          0.          0.025
    0.        ]
  [ 0.          0.          0.          0.          0.
    0.025     ]]]

is_cubic(tol=0.01)[source]

Check that the tensor corresponds to cubic symmetry, within a given tolerance.

The method relies on the multiplicity of eigenstiffnesses.

Parameters:: tol (float, optional) – Absolute tolerance to consider multiplicity of eigenstiffnesses
Returns:: If the tensor is single, the returned value is boolean. If the object is a tensor array, the returned value is an array of bools, the same shape as the tensor array.
Return type:: bool or numpy.ndarray

See also

is_isotropic: check if the stiffness tensor is isotropic
is_tetragonal: check if the stiffness tensor has tetragonal symmetry
eig_stiffnesses: compute eigenstiffnesses

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> from scipy.spatial.transform import Rotation
>>> C = StiffnessTensor.cubic(C11=186, C12=134, C44=77)
>>> C_rotated = C * Rotation.random(random_state=123)
>>> C_rotated
Stiffness tensor (in Voigt mapping):
[[237.71171578  96.41409344 119.87419078   8.1901353   -3.63846312
  -20.34233446]
 [ 96.41409344 250.74909842 106.83680814   9.33462785  -6.52548033
    0.99714278]
 [119.87419078 106.83680814 227.28900108 -17.52476315  10.16394345
   19.34519167]
 [  8.1901353    9.33462785 -17.52476315  49.83680814  19.34519167
   -6.52548033]
 [ -3.63846312  -6.52548033  10.16394345  19.34519167  62.87419078
    8.1901353 ]
 [-20.34233446   0.99714278  19.34519167  -6.52548033   8.1901353
   39.41409344]]

Once rotated, it is not clear if the stiffness tensors has cubic symmetry on sight. Yet:

>>> C_rotated.is_cubic()
True

is_isotropic(tol=0.01)[source]

Check that the tensor corresponds to isotropic symmetry, within a given tolerance.

The method relies on the multiplicity of eigenstiffnesses

Parameters:: tol (float) – Absolute tolerance to consider multiplicity of eigenstiffnesses
Returns:: If the tensor is single, the returned value is boolean. If the object is a tensor array, the returned value is an array of bools, the same shape as the tensor array.
Return type:: bool or numpy.ndarray

See also

is_cubic: check if the stiffness tensor has cubic symmetry
is_tetragonal: check if the stiffness tensor has tetragonal symmetry
eig_stiffnesses: compute eigenstiffnesses

is_tetragonal(tol=0.01)[source]

Check that the tensor corresponds to tetragonal symmetry, within a given tolerance.

The method relies on the multiplicity of eigenstiffnesses.

Parameters:: tol (float) – Absolute tolerance to consider multiplicity of eigenstiffnesses
Returns:: If the tensor is single, the returned value is boolean. If the object is a tensor array, the returned value is an array of bools, the same shape as the tensor array.
Return type:: bool or numpy.ndarray

See also

is_isotropic: check if the stiffness tensor is isotropic
is_cubic: check if the stiffness tensor has cubic symmetry
eig_stiffnesses: compute eigenstiffnesses

classmethod isotropic(E=None, nu=None, G=None, lame1=None, lame2=None, K=None, phase_name=None)[source]

Create an isotropic stiffness tensor.

The stiffness tensor must be constructed from exactly two elasticity coefficients, namely: E, nu, G, lame1, or lame2. Note that lame2 is just an alias for G. Each of these coefficient can be a list; in this case, the returned object is a tensor array (see examples).

Parameters:

E (float or list, None) – Young modulus
nu (float or list, None) – Poisson ratio
G (float or list, None) – Shear modulus
lame1 (float or list, None) – First Lamé coefficient
lame2 (float or list, None) – Second Lamé coefficient (alias for G)
K (float or list, None) – Bulk modulus
phase_name (str, None) – Name to print

Return type:

Corresponding isotropic stiffness tensor

See also

transverse_isotropic: create a transverse-isotropic tensor

Notes

The units you use when passing the elastic moduli must be consistent with that of the stress tensor. For instance, if you expect to work in MPa, the Young’s modulus and the Lamé’s coefficient must be given in MPa as well.

Examples

On can check that the shear modulus for steel is around 82 GPa:

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C=StiffnessTensor.isotropic(E=210e3, nu=0.28)
>>> C.shear_modulus
Hyperspherical function
Min=82031.24999999991, Max=82031.25000000006

Similarly, the same tensor can be defined from another pair of parameters, e.g. Young and shear moduli:

>>> C=StiffnessTensor.isotropic(E=210e3, G=82031)
>>> C.Young_modulus
Spherical function
Min=210000.0000000001, Max=210000.00000000032

One can build a tensor array by providing a list of values for the input arguments, instead of floats. For instance:

>>> C = StiffnessTensor.isotropic(E=(210, 70), nu=(0.28, 0.35)) # Elastic moduli for steel and aluminium
>>> C.shape
(2,)

We can easily check that the shear moduli for steel and aluminium are:

>>> C[0].shear_modulus
Hyperspherical function
Min=82.0312499999999, Max=82.03125000000006

>>> C[1].shear_modulus
Hyperspherical function
Min=25.925925925925895, Max=25.925925925925952

property lame1[source]

Compute the first Lamé’s parameter (only for isotropic materials).

If the stiffness/compliance tensor is not isotropic, NaN is returned.

Returns:: First Lamé’s parameter
Return type:: float

See also

lame2: second Lamé’s parameter

property lame2[source]

Compute the second Lamé’s parameter (only for isotropic materials).

If the stiffness/compliance tensor is not isotropic, NaN is returned.

Returns:: Second Lamé’s parameter
Return type:: float

See also

lame1: first Lamé’s parameter

property linear_compressibility[source]

Compute the directional linear compressibility.

Returns:: Directional linear compressibility
Return type:: SphericalFunction

See also

bulk_modulus: bulk modulus of the material

linear_invariants()[source]

Compute the linear invariants of the tensor, or tensor array.

If the object is a tensor array, the linear invariants are returned as arrays of each invariant. See notes for the actual definitions.

Returns:

A1 (float or np.ndarray) – First linear invariant
A2 (float or np.ndarray) – Second linear invariant

See also

quadratic_invariants: compute the quadratic invariants of a fourth-order tensor

Notes

The linear invariants are:

\[ \begin{align}\begin{aligned}A_1=C_{ijij}\\A_2=C_{iijj}\end{aligned}\end{align} \]

matrix(mapping_convention=None)[source]

Returns the components of the tensor as a matrix.

Parameters:: mapping_convention (VoigtMapping, optional) – Mapping convention to use for the returned matrix. If not provided, that of the tensor is used.
Returns:: Components of the tensor as a matrix
Return type:: numpy.ndarray

Examples

Create an identity 4th-order tensor:

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> t = FourthOrderTensor.eye()

Its matrix with respect to Kelvin mapping is:

>>> t.matrix()
array([[1., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 1.]])

whereas, when using the Voigt mapping, we have:

>>> from Elasticipy.tensors.mapping import VoigtMapping
>>> t.matrix(mapping_convention=VoigtMapping())
array([[1. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 1. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 1. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0.5, 0. , 0. ],
       [0. , 0. , 0. , 0. , 0.5, 0. ],
       [0. , 0. , 0. , 0. , 0. , 0.5]])

For stiffness tensors, the default mapping convention is Voigt, so that:

>>> from Elasticipy.tensors.elasticity import StiffnessTensor, ComplianceTensor
>>> StiffnessTensor.eye().matrix()
array([[1. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 1. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 1. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0.5, 0. , 0. ],
       [0. , 0. , 0. , 0. , 0.5, 0. ],
       [0. , 0. , 0. , 0. , 0. , 0.5]])

whereas for compliance tensor, the default mapping convention gives:

>>> ComplianceTensor.eye().matrix()
array([[1., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 2., 0., 0.],
       [0., 0., 0., 0., 2., 0.],
       [0., 0., 0., 0., 0., 2.]])

mean(axis=None)[source]

Compute the mean value of the tensor T

Parameters:: axis (int or list of int or tuple of int, optional) – axis along which to compute the mean. If None, the mean is computed on the flattened tensor
Returns:: If no axis is given, the result will be of shape (3,3,3,3).
Return type:: numpy.ndarray

Examples

Create a random tensor array of shape (5,6):

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> T = FourthOrderTensor.rand(shape=(5,6))
>>> Overall_mean = T.mean()
>>> Overall_mean.shape
()
>>> Overall_mean
4th-order tensor (in Kelvin mapping):
[[0.514295   0.52259217 0.42899181 0.77148692 0.64073221 0.73211491]
 [0.49422678 0.43718365 0.40786118 0.8170971  0.68435571 0.67262655]
 [0.48753674 0.51142541 0.44650454 0.76310921 0.67724973 0.69430165]
 [0.53946846 0.75101474 0.73578098 1.04338905 1.21598419 0.99489014]
 [0.75354555 0.61193555 0.82341479 1.11197826 0.89183143 1.20986243]
 [0.66078807 0.70126535 0.63719147 0.87567139 1.05671229 1.03004098]]

>>> axis_0_mean = T.mean(axis=0)
>>> axis_0_mean.shape
(6,)
>>> axis_1_mean = T.mean(axis=1)
>>> axis_1_mean.shape
(5,)

classmethod monoclinic(*, C11=0.0, C12=0.0, C13=0.0, C22=0.0, C23=0.0, C33=0.0, C44=0.0, C55=0.0, C66=0.0, C15=None, C25=None, C35=None, C46=None, C16=None, C26=None, C36=None, C45=None, phase_name=None)[source]

Create a fourth-order tensor from monoclinic symmetry. It automatically detects whether the components are given according to the Y or Z diad, depending on the input arguments.

For Diad || y, C15, C25, C35 and C46 must be provided. For Diad || z, C16, C26, C36 and C45 must be provided.

Parameters:

C11 (float) – Components of the tensor, using the Voigt notation
C12 (float) – Components of the tensor, using the Voigt notation
C13 (float) – Components of the tensor, using the Voigt notation
C22 (float) – Components of the tensor, using the Voigt notation
C23 (float) – Components of the tensor, using the Voigt notation
C33 (float) – Components of the tensor, using the Voigt notation
C44 (float) – Components of the tensor, using the Voigt notation
C55 (float) – Components of the tensor, using the Voigt notation
C66 (float) – Components of the tensor, using the Voigt notation
C15 (float, optional) – C15 component of the tensor (if Diad || y)
C25 (float, optional) – C25 component of the tensor (if Diad || y)
C35 (float, optional) – C35 component of the tensor (if Diad || y)
C46 (float, optional) – C46 component of the tensor (if Diad || y)
C16 (float, optional) – C16 component of the tensor (if Diad || z)
C26 (float, optional) – C26 component of the tensor (if Diad || z)
C36 (float, optional) – C36 component of the tensor (if Diad || z)
C45 (float, optional) – C45 component of the tensor (if Diad || z)
phase_name (str, optional) – Name to display

Return type:

FourthOrderTensor

See also

triclinic: create a tensor from triclinic symmetry
orthorhombic: create a tensor from orthorhombic symmetry

property ndim[source]

Returns the dimensionality of the tensor (number of dimensions in the orientation array)

Returns:: Number of dimensions
Return type:: int

classmethod ones(shape=None, **kwargs)[source]

Create a 4th-order tensor full of ones.

Parameters:

shape (int or tuple, optional) – Shape of the tensor to create
kwargs – keyword arguments passed to the constructor

Return type:

FourthOrderTensor

Examples

>>> tensor_of_ones = FourthOrderTensor.ones()
>>> tensor_of_ones
4th-order tensor (in Kelvin mapping):
[[1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]]

At first sight, the tensor may appear not full of ones at all, but the representation above uses the Kelvin mapping convention. Indeed, one can check that the full tensor is actually full of ones. E.g.:

>>> tensor_of_ones.full_tensor[0,1,0,2]
1.0

Alternatively, the Voigt mapping convention may help figuring it out:

>>> from Elasticipy.tensors.mapping import VoigtMapping
>>> tensor_of_ones_voigt = FourthOrderTensor.ones(mapping=VoigtMapping())
>>> tensor_of_ones_voigt
4th-order tensor (in Voigt mapping):
[[1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]]

although both tensors are actually the same:

>>> tensor_of_ones == tensor_of_ones_voigt
True

classmethod orthorhombic(*, C11=0.0, C12=0.0, C13=0.0, C22=0.0, C23=0.0, C33=0.0, C44=0.0, C55=0.0, C66=0.0, phase_name=None)[source]

Create a fourth-order tensor from orthorhombic symmetry.

Parameters:

C11 (float) – Components of the tensor, using the Voigt notation
C12 (float) – Components of the tensor, using the Voigt notation
C13 (float) – Components of the tensor, using the Voigt notation
C22 (float) – Components of the tensor, using the Voigt notation
C23 (float) – Components of the tensor, using the Voigt notation
C33 (float) – Components of the tensor, using the Voigt notation
C44 (float) – Components of the tensor, using the Voigt notation
C55 (float) – Components of the tensor, using the Voigt notation
C66 (float) – Components of the tensor, using the Voigt notation
phase_name (str, optional) – Phase name to display

Return type:

FourthOrderTensor

See also

monoclinic: create a tensor from monoclinic symmetry
orthorhombic: create a tensor from orthorhombic symmetry

classmethod orthotropic(*, Ex, Ey, Ez, Gxy, Gxz, Gyz, nu_yx=None, nu_zx=None, nu_zy=None, nu_xy=None, nu_xz=None, nu_yz=None, **kwargs)[source]

Create a stiffness tensor corresponding to orthotropic symmetry, given the engineering constants.

Exactly three Poisson ratios must be provided. See Notes for details.

Parameters:

Ex (float) – Young modulus along the x axis
Ey (float) – Young modulus along the y axis
Ez (float) – Young modulus along the z axis
Gxy (float) – Shear modulus in the x-y plane
Gxz (float) – Shear modulus in the x-z plane
Gyz (float) – Shear modulus in the y-z plane
nu_xy (float, optional) – Poisson ratio along x and y axes. Either nu_xy or nu_yx must be provided, not both.
nu_yx (float, optional) – Poisson ratio along x and y axes. Either nu_xy or nu_yx must be provided, not both.
nu_xz (float, optional) – Poisson ratio along x and z axes. Either nu_xz or nu_zx must be provided, not both.
nu_zx (float, optional) – Poisson ratio along x and z axes. Either nu_xz or nu_zx must be provided, not both.
nu_yz (float, optional) – Poisson ratio along y and z axes. Either nu_yz or nu_zy must be provided, not both.
nu_zy (float, optional) – Poisson ratio along y and z axes. Either nu_yz or nu_zy must be provided, not both.
kwargs (dict, optional) – Keyword arguments to pass to the StiffnessTensor constructor

Return type:

StiffnessTensor

See also

transverse_isotropic: create a stiffness tensor for transverse-isotropic symmetry

Notes

If the material undergoes tensile strain \(\varepsilon_{ii}\) along the i-th direction, the Poisson ratios are defined as:

\[\nu_{ij}=-\frac{\partial \varepsilon_{jj}}{\partial \varepsilon_{ii}}\]

where \(\varepsilon_{jj}\) denotes the (compressive) longitudinal strain along the j-th direction. If \(E_x\) and \(E_y\) are the Young moduli along x and y, we have:

\[\frac{\nu_{xy}}{E_x} = \frac{\nu_{yx}}{E_y}\]

quadratic_invariants()[source]

Compute the quadratic invariants of the tensor, or tensor array.

If the object is a tensor array, the returned values are arrays of each invariant. See notes for definitions.

Returns:: B1, B2, B3, B4, B5
Return type:: float or np.ndarray

See also

linear_invariants: compute the linear invariants of a Fourth-order tensor

Notes

The quadratic invariants are defined as [Norris]:

\[ \begin{align}\begin{aligned}B_1 = C_{ijkl}C_{ijkl}\\B_2 = C_{iikl}C_{jjkl}\\B_3 = C_{iikl}C_{jkjl}\\B_4 = C_{kiil}C_{kjjl}\\B_5 = C_{ijkl}C_{ikjl}\end{aligned}\end{align} \]

References

[Norris]

Norris, A. N. (22 May 2007). “Quadratic invariants of elastic moduli”. The Quarterly Journal of Mechanics and Applied Mathematics. 60 (3): 367–389. doi:10.1093/qjmam/hbm007

classmethod rand(shape=None, **kwargs)[source]

Populate a Fourth-order tensor with random values in half-open interval [0.0, 1.0).

Parameters:

shape (tuple or int, optional) – Set the shape of the tensor array. If None, the returned tensor will be single.
kwargs – Keyword arguments passed to the Fourth-order tensor constructor.

Returns:

Fourth-order tensor

Return type:

FourthOrderTensor

rotate(rotation)[source]

Apply a single rotation to a tensor, and return its component into the rotated frame.

Parameters:: rotation (Rotation or orix.quaternion.rotation.Rotation) – Rotation to apply
Returns:: Rotated tensor
Return type:: SymmetricFourthOrderTensor

Examples

Let start from a given tensor, (say ones):

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> T = FourthOrderTensor.ones()
>>> T
4th-order tensor (in Kelvin mapping):
[[1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]]

Define a rotation. E.g.:

>>> from scipy.spatial.transform import Rotation
>>> g = Rotation.from_euler('X', 90, degrees=True)

Then , apply rotation:

>>> Trotated = T.rotate(g)
>>> Trotated
4th-order tensor (in Kelvin mapping):
[[ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [-1.41421356 -1.41421356 -1.41421356  2.         -2.          2.        ]
 [ 1.41421356  1.41421356  1.41421356 -2.          2.         -2.        ]
 [-1.41421356 -1.41421356 -1.41421356  2.         -2.          2.        ]]

Actually, a more simple syntax is:

>>> T * g
4th-order tensor (in Kelvin mapping):
[[ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [ 1.          1.          1.         -1.41421356  1.41421356 -1.41421356]
 [-1.41421356 -1.41421356 -1.41421356  2.         -2.          2.        ]
 [ 1.41421356  1.41421356  1.41421356 -2.          2.         -2.        ]
 [-1.41421356 -1.41421356 -1.41421356  2.         -2.          2.        ]]

Obviously, the original tensor can be retrieved by applying the reverse rotation:

>>> Trotated * g.inv()
4th-order tensor (in Kelvin mapping):
[[1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.         1.         1.         1.41421356 1.41421356 1.41421356]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]
 [1.41421356 1.41421356 1.41421356 2.         2.         2.        ]]

If g is composed of multiple rotations, this will result in a tensor array, corresponding to each rotation:

>>> import numpy as np
>>> theta = np.linspace(0, 90, 100)
>>> g = Rotation.from_euler('X', theta, degrees=True)
>>> Trotated = T * g
>>> Trotated
4th-order tensor array of shape (100,)

save_to_txt(filename, matrix_only=False)[source]

Save the tensor to a text file.

Parameters:

filename (str) – Filename to save the tensor to.
matrix_only (bool, False) – If true, only the components of the stiffness tensor is saved (no data about phase nor symmetry)

See also

from_txt_file: create a tensor from text file

property shape[source]

Return the shape of the tensor array

Returns:: Shape of the tensor array
Return type:: tuple

property shear_modulus[source]

Directional shear modulus

Returns:: Shear modulus
Return type:: HyperSphericalFunction

spherical_part()[source]

Return the spherical part of the tensor

Returns:: Spherical part of the tensor
Return type:: FourthOrderTensor

See also

identity_tensor: return the identity tensor
deviatoric_part: return the deviatoric part of the tensor

classmethod tetragonal(*, C11=0.0, C12=0.0, C13=0.0, C33=0.0, C44=0.0, C16=0.0, C66=0.0, phase_name=None)[source]

Create a fourth-order tensor from tetragonal symmetry.

Parameters:

C11 (float) – Components of the tensor, using the Voigt notation
C12 (float) – Components of the tensor, using the Voigt notation
C13 (float) – Components of the tensor, using the Voigt notation
C33 (float) – Components of the tensor, using the Voigt notation
C44 (float) – Components of the tensor, using the Voigt notation
C66 (float) – Components of the tensor, using the Voigt notation
C16 (float, optional) – C16 component in Voigt notation (for point groups 4, -4 and 4/m only)
phase_name (str, optional) – Phase name to display

Return type:

FourthOrderTensor

See also

trigonal: create a tensor from trigonal symmetry
orthorhombic: create a tensor from orthorhombic symmetry

to_Kelvin()[source]

Returns all the tensor components using the Kelvin(-Mandel) mapping convention.

Returns:: (6,6) matrix, according to the Kelvin mapping
Return type:: numpy.ndarray

See also

eig: returns the eigenvalues and the eigenvectors of the Kelvin’s matrix
from_Kelvin: Construct a fourth-order tensor from its (6,6) Kelvin matrix

Notes

This mapping convention is defined as follows [Helbig]:

\[\begin{split}C_K = \begin{bmatrix} C_{11} & C_{12} & C_{13} & \sqrt{2}C_{14} & \sqrt{2}C_{15} & \sqrt{2}C_{16}\\ C_{12} & C_{22} & C_{23} & \sqrt{2}C_{24} & \sqrt{2}C_{25} & \sqrt{2}C_{26}\\ C_{13} & C_{23} & C_{33} & \sqrt{2}C_{34} & \sqrt{2}C_{35} & \sqrt{2}C_{36}\\ \sqrt{2}C_{14} & \sqrt{2}C_{24} & \sqrt{2}C_{34} & 2C_{44} & 2C_{45} & 2C_{46}\\ \sqrt{2}C_{15} & \sqrt{2}C_{25} & \sqrt{2}C_{35} & 2C_{45} & 2C_{55} & 2C_{56}\\ \sqrt{2}C_{16} & \sqrt{2}C_{26} & \sqrt{2}C_{36} & 2C_{46} & 2C_{56} & 2C_{66}\\ \end{bmatrix}\end{split}\]

References

[Helbig]

Helbig, K. (2013). What Kelvin might have written about Elasticity. Geophysical Prospecting, 61(1), 1-20. doi: 10.1111/j.1365-2478.2011.01049.x

Examples

>>> from Elasticipy.tensors.elasticity import StiffnessTensor
>>> C = StiffnessTensor.cubic(C11=200, C12=40, C44=20)
>>> print(C)
Stiffness tensor (in Voigt mapping):
[[200.  40.  40.   0.   0.   0.]
 [ 40. 200.  40.   0.   0.   0.]
 [ 40.  40. 200.   0.   0.   0.]
 [  0.   0.   0.  20.   0.   0.]
 [  0.   0.   0.   0.  20.   0.]
 [  0.   0.   0.   0.   0.  20.]]
>>> C.to_Kelvin()
array([[200.,  40.,  40.,   0.,   0.,   0.],
       [ 40., 200.,  40.,   0.,   0.,   0.],
       [ 40.,  40., 200.,   0.,   0.,   0.],
       [  0.,   0.,   0.,  40.,   0.,   0.],
       [  0.,   0.,   0.,   0.,  40.,   0.],
       [  0.,   0.,   0.,   0.,   0.,  40.]])

to_pymatgen()[source]

Convert the stiffness tensor (from Elasticipy) to Python Materials Genomics (Pymatgen) format.

Returns:: Stiffness tensor for pymatgen
Return type:: pymatgen.analysis.elasticity.elastic.ElasticTensor

transpose_array()[source]

Transpose the orientations of the tensor array

Returns:: The same tensor, but with transposed axes
Return type:: FourthOrderTensor

Examples

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> A = FourthOrderTensor.rand(shape=(3,4))
>>> A.transpose_array()
4th-order tensor array of shape (4, 3)

classmethod transverse_isotropic(*, Ex, Ez, Gxz, nu_yx=None, nu_xy=None, nu_zx=None, nu_xz=None, **kwargs)[source]

Create a stiffness tensor corresponding to the transversely isotropic symmetry with respect to Z axis, given the engineering constants.

Exactly two Poisson ratios must be provided (nu_xy or nu_yx, and nu_xz or nu_zx). See Notes for details.

Parameters:

Ex (float) – Young modulus along the x axis
Ez (float) – Young modulus along the y axis
Gxz (float) – Shear modulus in the x-z plane
nu_xy (float, optional) – Poisson ratio along x and y. Either nu_xy or nu_yx must be provided, not both.
nu_yx (float, optional) – Poisson ratio along x and y. Either nu_xy or nu_yx must be provided, not both.
nu_xz (float, optional) – Poisson ratio along x and z. Either nu_xz or nu_zx must be provided, not both.
nu_zx (float, optional) – Poisson ratio along x and z. Either nu_xz or nu_zx must be provided, not both.
kwargs (dict) – Keyword arguments to pass to the StiffnessTensor constructor

Return type:

StiffnessTensor

See also

orthotropic: create a stiffness tensor for orthotropic symmetry

Notes

If the material undergoes tensile strain \(\varepsilon_{ii}\) along the i-th direction, the Poisson ratios are defined as:

\[\nu_{ij}=-\frac{\partial \varepsilon_{jj}}{\partial \varepsilon_{ii}}\]

where \(\varepsilon_{jj}\) denotes the (compressive) longitudinal strain along the j-th direction. If \(E_x\) and \(E_y\) are the Young moduli along x and y, we have:

\[\frac{\nu_{xy}}{E_x} = \frac{\nu_{yx}}{E_y}\]

classmethod triclinic(C11=0.0, C12=0.0, C13=0.0, C14=0.0, C15=0.0, C16=0.0, C22=0.0, C23=0.0, C24=0.0, C25=0.0, C26=0.0, C33=0.0, C34=0.0, C35=0.0, C36=0.0, C44=0.0, C45=0.0, C46=0.0, C55=0.0, C56=0.0, C66=0.0, phase_name=None)[source]

Parameters:

C11 (float) – Components of the tensor
C12 (float) – Components of the tensor
C13 (float) – Components of the tensor
C14 (float) – Components of the tensor
C15 (float) – Components of the tensor
C16 (float) – Components of the tensor
C22 (float) – Components of the tensor
C23 (float) – Components of the tensor
C24 (float) – Components of the tensor
C25 (float) – Components of the tensor
C26 (float) – Components of the tensor
C33 (float) – Components of the tensor
C34 (float) – Components of the tensor
C35 (float) – Components of the tensor
C36 (float) – Components of the tensor
C44 (float) – Components of the tensor
C45 (float) – Components of the tensor
C46 (float) – Components of the tensor
C55 (float) – Components of the tensor
C56 (float) – Components of the tensor
C66 (float) – Components of the tensor
phase_name (str, optional) – Name to display

Return type:

FourthOrderTensor

See also

monoclinic: create a tensor from monoclinic symmetry
orthorhombic: create a tensor from orthorhombic symmetry

classmethod trigonal(*, C11=0.0, C12=0.0, C13=0.0, C14=0.0, C33=0.0, C44=0.0, C15=0.0, phase_name=None)[source]

Create a fourth-order tensor from trigonal symmetry.

Parameters:

C11 (float) – Components of the tensor, using the Voigt notation
C12 (float) – Components of the tensor, using the Voigt notation
C13 (float) – Components of the tensor, using the Voigt notation
C14 (float) – Components of the tensor, using the Voigt notation
C33 (float) – Components of the tensor, using the Voigt notation
C44 (float) – Components of the tensor, using the Voigt notation
C15 (float, optional) – C15 component of the tensor, only used for point groups 3 and -3.
phase_name (str, optional) – Phase name to display

Return type:

FourthOrderTensor

See also

tetragonal: create a tensor from tetragonal symmetry
orthorhombic: create a tensor from orthorhombic symmetry

property universal_anisotropy[source]

Compute the universal anisotropy factor.

The larger the value, the more likely the material will behave in an anisotropic way.

Returns:: The universal anisotropy factor.
Return type:: float

Notes

The universal anisotropy factor is defined as [3]:

\[5\frac{G_v}{G_r} + \frac{K_v}{K_r} - 6\]

References

wave_velocity(rho)[source]

Compute the wave velocities, given the mass density.

Parameters:: rho (float) – mass density. Its unit must be consistent with that of the stiffness tensor. See notes for hints.

See also

Christoffel_tensor: Computes the Christoffel tensor along a given direction

Returns:

c_p (SphericalFunction) – Velocity of the primary (compressive) wave
c_s1 (SphericalFunction) – Velocity of the fast secondary (shear) wave
c_s2 (SphericalFunction) – Velocity of the slow secondary (shear) wave

Notes

The estimation of the wave velocities is made by finding the eigenvalues of the Christoffel tensor [Jaeken].

One should double-check the units. The table below provides hints about the unit you get, depending on the units you use for stiffness and the mass density:

Stiffness	Mass density	Velocities	Notes
Pa (N/m²)	kg/m³	m/s	SI units
GPa (10⁹ Pa)	kg/dm³	km/s	Conversion factor
GPa (10³ N/mm²)	kg/mm³	m/s	Consistent units
MPa (10⁶ Pa)	kg/m³	km/s	Conversion factor
MPa (10³ N/mm²)	g/mm³	m/s	Consistent units

References

[Jaeken]

J. W. Jaeken, S. Cottenier, Solving the Christoffel equation: Phase and group velocities, Computer Physics Communications (207), 2016, https://doi.org/10.1016/j.cpc.2016.06.014.

Examples

We will investigate the wave velocities in cubic Cu. First, define the stiffness tensor:

>>> C = StiffnessTensor.cubic(C11=186e3, C12=134e3, C44=77e3) # Stiffness in MPa
>>> rho = 8960 # mass density in kg/m³
>>> c_p, c_s1, c_s2 = C.wave_velocity(rho)

The range of velocities of the (fast) primary wave range can be printed with:

>>> print(c_p)
Spherical function
Min=4.556196722204669, Max=5.324304112812698

As the stiffness is given in MPa and the mass desity is given in kg/m³, the velocities are returned in km/s (see Notes). If one wants to know the direction for min/max values:

>>> val_min, u_min = c_p.min()
>>> print(u_min)
[[0. 0. 1.]]

>>> val_max, u_max = c_p.max()
>>> print(u_max)
[[0.57735022 0.57735033 0.57735026]]

The secondary wave velocities are given by cs_1 and cs_2:

>>> c_s1
Spherical function
Min=2.1906864565414192, Max=2.9315098498896446

>>> c_s2
Spherical function
Min=1.7034628596749235, Max=2.9315098498896437

classmethod weighted_average(Cs, volume_fractions, method)[source]

Compute the weighted average of a list of stiffness tensors, with respect to a given method (Voigt, Reuss or Hill).

Parameters:

Cs (list of StiffnessTensor or list of ComplianceTensor or tuple of StiffnessTensor or tuple of ComplianceTensor) – Series of tensors to compute the average from
volume_fractions (iterable of floats) – Volume fractions of each phase
method (str, {'Voigt', 'Reuss', 'Hill'}) – Method to use. It can be ‘Voigt’, ‘Reuss’, or ‘Hill’.

Returns:

Average tensor

Return type:

StiffnessTensor

classmethod zeros(shape=(), **kwargs)[source]

Create a fourth-order tensor populated with zeros

Parameters:

shape (int or tuple, optional) – Shape of the tensor to create
kwargs – Keyword arguments passed to the FourthOrderTensor constructor

Return type:

FourthOrderTensor

Examples

The single-valued null 4th order tensor is just:

>>> from Elasticipy.tensors.fourth_order import FourthOrderTensor
>>> FourthOrderTensor.zeros()
4th-order tensor (in Kelvin mapping):
[[0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]]

One can also create an array of such tensors:

>>> zeros_tensor = FourthOrderTensor.zeros(shape=3)

and check that it populated with zeros:

>>> zeros_tensor == 0.
array([ True,  True,  True])