pydrex.tensors

  1"""> PyDRex: Tensor operation functions and helpers.
  2
  3For Voigt notation, the symmetric 6x6 matrix representation is used,
  4which assumes that the fourth order tensor being represented as such is also symmetric.
  5The vectorial notation uses 21 components which are the independent components of the
  6symmetric 6x6 matrix.
  7
  8"""
  9
 10import numba as nb
 11import numpy as np
 12
 13
 14@nb.njit(fastmath=True)
 15def polar_decompose(matrix, left=True):
 16    """Compute polar decomposition of M as either M = RU (right) or M = VP (left)."""
 17    U, S, Vh = np.linalg.svd(matrix)
 18    if left:
 19        return U @ Vh, U @ (np.diag(S) @ U.transpose())
 20    U_matrix = Vh.transpose() @ (np.diag(S) @ Vh)
 21    return matrix @ np.linalg.inv(U_matrix), U_matrix
 22
 23
 24@nb.njit(fastmath=True)
 25def invariants_second_order(tensor):
 26    """Calculate invariants of a second order tensor."""
 27    return (
 28        np.trace(tensor),
 29        tensor[0, 0] * tensor[1, 1]
 30        + tensor[1, 1] * tensor[2, 2]
 31        + tensor[2, 2] * tensor[0, 0]
 32        - tensor[0, 1] * tensor[1, 0]
 33        - tensor[1, 2] * tensor[2, 1]
 34        - tensor[2, 0] * tensor[0, 2],
 35        np.linalg.det(tensor),
 36    )
 37
 38
 39@nb.njit(fastmath=True)
 40def voigt_decompose(matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
 41    """Decompose elastic tensor (as 6x6 Voigt matrix) into distinct contractions.
 42
 43    Return the only two independent contractions of the elastic tensor given as a 6x6
 44    Voigt `matrix`. With reference to the full 4-th order elastic tensor, the
 45    contractions are defined as:
 46    - $d_{ij} = C_{ijkk}$ (dilatational stiffness tensor)
 47    - $v_{ij} = C_{ijkj}$ (deviatoric stiffness tensor)
 48
 49    Any vector which is an eigenvector of both $d_{ij}$ and $v_{ij}$ is always normal to
 50    a symmetry plane of the elastic medium.
 51
 52    See Equations 3.4 & 3.5 in [Browaeys & Chevrot (2004)](https://doi.org/10.1111/j.1365-246X.2004.02415.x).
 53
 54    """
 55    # 1. Compose dᵢⱼ = Cᵢⱼₖₖ (dilatational stiffness tensor)
 56    # Eq. 3.4 in Browaeys & Chevrot, 2004.
 57    stiffness_dilat = np.empty((3, 3))
 58    for i in range(3):
 59        stiffness_dilat[i, i] = matrix[:3, i].sum()
 60    stiffness_dilat[0, 1] = stiffness_dilat[1, 0] = matrix[:3, 5].sum()
 61    stiffness_dilat[0, 2] = stiffness_dilat[2, 0] = matrix[:3, 4].sum()
 62    stiffness_dilat[1, 2] = stiffness_dilat[2, 1] = matrix[:3, 3].sum()
 63    # 2. Compose vᵢⱼ = Cᵢⱼₖⱼ (deviatoric stiffness tensor)
 64    # Eq. 3.5, Browaeys & Chevrot, 2004.
 65    stiffness_deviat = np.empty((3, 3))
 66    stiffness_deviat[0, 0] = matrix[0, 0] + matrix[4, 4] + matrix[5, 5]
 67    stiffness_deviat[1, 1] = matrix[1, 1] + matrix[3, 3] + matrix[5, 5]
 68    stiffness_deviat[2, 2] = matrix[2, 2] + matrix[3, 3] + matrix[4, 4]
 69    stiffness_deviat[0, 1] = matrix[0, 5] + matrix[1, 5] + matrix[3, 4]
 70    stiffness_deviat[0, 2] = matrix[0, 4] + matrix[2, 4] + matrix[3, 5]
 71    stiffness_deviat[1, 2] = matrix[1, 3] + matrix[2, 3] + matrix[4, 5]
 72    stiffness_deviat = upper_tri_to_symmetric(stiffness_deviat)
 73    return stiffness_dilat, stiffness_deviat
 74
 75
 76@nb.njit(fastmath=True)
 77def mono_project(voigt_vector):
 78    """Project 21-component `voigt_vector` onto monoclinic symmetry subspace.
 79
 80    Monoclinic symmetry is characterised by 13 independent elasticity components.
 81
 82    See [Browaeys & Chevrot (2004)](https://doi.org/10.1111/j.1365-246X.2004.02415.x).
 83
 84    """
 85    out = voigt_vector.copy()
 86    inds = (9, 10, 12, 13, 15, 16, 18, 19)
 87    for i in inds:
 88        out[i] = 0
 89    return out
 90
 91
 92@nb.njit(fastmath=True)
 93def ortho_project(voigt_vector):
 94    """Project 21-component `voigt_vector` onto orthorhombic symmetry subspace.
 95
 96    Orthorhombic symmetry is characterised by 9 independent elasticity components.
 97
 98    See [Browaeys & Chevrot (2004)](https://doi.org/10.1111/j.1365-246X.2004.02415.x).
 99
100    """
101    out = voigt_vector.copy()
102    out[9:] = 0
103    return out
104
105
106@nb.njit(fastmath=True)
107def tetr_project(voigt_vector):
108    """Project 21-component `voigt_vector` onto tetragonal symmetry subspace.
109
110    Tetragonal symmetry is characterised by 6 independent elasticity components.
111
112    See [Browaeys & Chevrot (2004)](https://doi.org/10.1111/j.1365-246X.2004.02415.x).
113
114    """
115    out = ortho_project(voigt_vector)
116    for i, j in ((0, 1), (3, 4), (6, 7)):
117        for k in range(2):
118            out[i + k] = 0.5 * (voigt_vector[i] + voigt_vector[j])
119    return out
120
121
122@nb.njit(fastmath=True)
123def hex_project(voigt_vector):
124    """Project 21-component `voigt_vector` onto hexagonal symmetry subspace.
125
126    Hexagonal symmetry (a.k.a. transverse isotropy) is characterised by 5 independent
127    elasticity components.
128
129    See [Browaeys & Chevrot (2004)](https://doi.org/10.1111/j.1365-246X.2004.02415.x).
130
131    """
132    x = voigt_vector
133    out = np.zeros(21)
134    out[0] = out[1] = 3 / 8 * (x[0] + x[1]) + x[5] / 4 / np.sqrt(2) + x[8] / 4
135    out[2] = x[2]
136    out[3] = out[4] = (x[3] + x[4]) / 2
137    out[5] = (x[0] + x[1]) / 4 / np.sqrt(2) + 3 / 4 * x[5] - x[8] / 2 / np.sqrt(2)
138    out[6] = out[7] = (x[6] + x[7]) / 2
139    out[8] = (x[0] + x[1]) / 4 - x[5] / 2 / np.sqrt(2) + x[8] / 2
140    return out
141
142
143@nb.njit(fastmath=True)
144def upper_tri_to_symmetric(arr):
145    """Create symmetric array using upper triangle of input array.
146
147    Examples:
148
149    >>> import numpy as np
150    >>> upper_tri_to_symmetric(np.array([
151    ...         [ 1.,  2.,  3.,  4.],
152    ...         [ 0.,  5.,  6.,  7.],
153    ...         [ 0.,  0.,  8.,  9.],
154    ...         [ 9.,  0.,  0., 10.]
155    ... ]))
156    array([[ 1.,  2.,  3.,  4.],
157           [ 2.,  5.,  6.,  7.],
158           [ 3.,  6.,  8.,  9.],
159           [ 4.,  7.,  9., 10.]])
160
161    """
162    # <https://stackoverflow.com/questions/58718365/fast-way-to-convert-upper-triangular-matrix-into-symmetric-matrix>
163    upper_tri = np.triu(arr)
164    return np.where(upper_tri, upper_tri, upper_tri.transpose())
165
166
167@nb.njit(fastmath=True)
168def voigt_to_elastic_tensor(matrix):
169    """Create 4-th order elastic tensor from an equivalent Voigt matrix.
170
171    See also: `elastic_tensor_to_voigt`.
172
173    """
174    tensor = np.empty((3, 3, 3, 3))
175    for p in range(3):
176        for q in range(3):
177            delta_pq = 1 if p == q else 0
178            i = (p + 1) * delta_pq + (1 - delta_pq) * (7 - p - q) - 1
179            for r in range(3):
180                for s in range(3):
181                    delta_rs = 1 if r == s else 0
182                    j = (r + 1) * delta_rs + (1 - delta_rs) * (7 - r - s) - 1
183                    tensor[p, q, r, s] = matrix[i, j]
184    return tensor
185
186
187@nb.njit(fastmath=True)
188def elastic_tensor_to_voigt(tensor):
189    """Create a 6x6 Voigt matrix from an equivalent 4-th order elastic tensor."""
190    matrix = np.zeros((6, 6))
191    matrix_indices = np.zeros((6, 6))
192    for p in range(3):
193        for q in range(3):
194            delta_pq = 1 if p == q else 0
195            i = (p + 1) * delta_pq + (1 - delta_pq) * (7 - p - q) - 1
196            for r in range(3):
197                for s in range(3):
198                    delta_rs = 1 if r == s else 0
199                    j = (r + 1) * delta_rs + (1 - delta_rs) * (7 - r - s) - 1
200                    matrix[i, j] += tensor[p, q, r, s]
201                    matrix_indices[i, j] += 1
202    matrix /= matrix_indices
203    return (matrix + matrix.transpose()) / 2
204
205
206@nb.njit(fastmath=True)
207def voigt_matrix_to_vector(matrix):
208    """Create the 21-component Voigt vector equivalent to the 6x6 Voigt matrix."""
209    vector = np.zeros(21)
210    for i in range(3):
211        vector[i] = matrix[i, i]
212        vector[i + 3] = np.sqrt(2) * matrix[(i + 1) % 3, (i + 2) % 3]
213        vector[i + 6] = 2 * matrix[i + 3, i + 3]
214        vector[i + 9] = 2 * matrix[i, i + 3]
215        vector[i + 12] = 2 * matrix[(i + 2) % 3, i + 3]
216        vector[i + 15] = 2 * matrix[(i + 1) % 3, i + 3]
217        vector[i + 18] = 2 * np.sqrt(2) * matrix[(i + 1) % 3 + 3, (i + 2) % 3 + 3]
218    return vector
219
220
221@nb.njit(fastmath=True)
222def voigt_vector_to_matrix(vector):
223    """Create the 6x6 matrix representation of the 21-component Voigt vector.
224
225    See also: `voigt_matrix_to_vector`.
226
227    """
228    matrix = np.zeros((6, 6))
229    for i in range(3):
230        matrix[i, i] = vector[i]
231        matrix[i + 3, i + 3] = 0.5 * vector[i + 6]
232
233    matrix[1, 2] = 1 / np.sqrt(2) * vector[3]
234    matrix[0, 2] = 1 / np.sqrt(2) * vector[4]
235    matrix[0, 1] = 1 / np.sqrt(2) * vector[5]
236
237    matrix[0, 3] = 0.5 * vector[9]
238    matrix[1, 4] = 0.5 * vector[10]
239    matrix[2, 5] = 0.5 * vector[11]
240    matrix[2, 3] = 0.5 * vector[12]
241
242    matrix[0, 4] = 0.5 * vector[13]
243    matrix[1, 5] = 0.5 * vector[14]
244    matrix[1, 3] = 0.5 * vector[15]
245
246    matrix[2][4] = 0.5 * vector[16]
247    matrix[0][5] = 0.5 * vector[17]
248    matrix[4][5] = 0.5 * 1 / np.sqrt(2) * vector[18]
249    matrix[3][5] = 0.5 * 1 / np.sqrt(2) * vector[19]
250    matrix[3][4] = 0.5 * 1 / np.sqrt(2) * vector[20]
251    return upper_tri_to_symmetric(matrix)
252
253
254@nb.njit(fastmath=True)
255def rotate(tensor, rotation):
256    """Rotate 4-th order tensor using a 3x3 rotation matrix."""
257    rotated_tensor = np.zeros((3, 3, 3, 3))
258    for i in range(3):
259        for j in range(3):
260            for k in range(3):
261                for L in range(3):
262                    for a in range(3):
263                        for b in range(3):
264                            for c in range(3):
265                                for d in range(3):
266                                    rotated_tensor[i, j, k, L] += (
267                                        rotation[i, a]
268                                        * rotation[j, b]
269                                        * rotation[k, c]
270                                        * rotation[L, d]
271                                        * tensor[a, b, c, d]
272                                    )
273    return rotated_tensor