pydrex.diagnostics

  1r"""> PyDRex: Methods to calculate texture and strain diagnostics.
  2
  3.. note::
  4    Calculations expect orientation matrices $a$ to represent passive
  5    (i.e. alias) rotations, which are defined in terms of the extrinsic ZXZ
  6    euler angles $ϕ, θ, φ$ as
  7    $$
  8    a = \begin{bmatrix}
  9            \cosφ\cosϕ - \cosθ\sinϕ\sinφ & \cosθ\cosϕ\sinφ + \cosφ\sinϕ & \sinφ\sinθ \cr
 10           -\sinφ\cosϕ - \cosθ\sinϕ\cosφ & \cosθ\cosϕ\cosφ - \sinφ\sinϕ & \cosφ\sinθ \cr
 11            \sinθ\sinϕ & -\sinθ\cosϕ & \cosθ
 12        \end{bmatrix}
 13    $$
 14    such that a[i, j] gives the direction cosine of the angle between the i-th
 15    grain axis and the j-th external axis (in the global Eulerian frame).
 16
 17"""
 18
 19import functools as ft
 20
 21import numba as nb
 22import numpy as np
 23import scipy.linalg as la
 24
 25from pydrex import geometry as _geo
 26from pydrex import logger as _log
 27from pydrex import stats as _stats
 28from pydrex import tensors as _tensors
 29from pydrex import utils as _utils
 30
 31Pool, HAS_RAY = _utils.import_proc_pool()
 32if HAS_RAY:
 33    import ray
 34
 35    from pydrex import distributed as _dstr
 36
 37
 38def elasticity_components(voigt_matrices: np.ndarray):
 39    """Calculate elasticity decompositions for the given elasticity tensors.
 40
 41    - `voigt_matrices` — the Nx6x6 Voigt matrix representations of the averaged
 42      elasticity tensors for a series of N polycrystal textures
 43
 44    Returns a dictionary with the following data series:
 45    - `'bulk_modulus'` — the computed bulk modulus for each Voigt matrix C
 46    - `'shear_modulus'` — the computed shear modulus for each Voigt matrix C
 47    - `'percent_anisotropy'` — for each Voigt matrix C, the total percentage of the
 48      norm of the elastic tensor ||C|| that deviates from the norm of the "closest"
 49      isotropic elasticity tensor
 50    - `'percent_hexagonal'` — for each C, the percentage of ||C|| attributed to
 51      a hexagonally symmetric bulk medium
 52    - `'percent_tetragonal'` — for each C, the percentage of ||C|| attributed to
 53      a tetragonally symmetric bulk medium
 54    - `'percent_orthorhombic'` — for each C, the percentage of ||C|| attributed to
 55      an orthorhombically symmetric bulk medium
 56    - `'percent_monoclinic'` — for each C, the percentage of ||C|| attributed to
 57      a monoclinically symmetric bulk medium
 58    - `'percent_triclinic'` — for each C, the percentage of ||C|| attributed to
 59      a triclinically "symmetric" bulk medium (no mirror planes)
 60    - `'hexagonal_axis'` — for each C, the axis of hexagonal symmetry for the "closest"
 61      hexagonally symmetric approximation to C, a.k.a. the "transverse isotropy" axis
 62
 63    .. note::
 64        Only 5 symmetry classes are relevant for elasticity decomposition,
 65        compared to the usual 6 used to describe crystal families.
 66        Crystals with cubic symmetry contribute to the isotropic elasticity tensor,
 67        because the lattice spacing is identical in all orthogonal directions.
 68        Note also that the trigonal crystal *system* is not a crystal family
 69        (it belongs to the hexagonal family).
 70
 71    """
 72    n_matrices = len(voigt_matrices)
 73    out = {
 74        "bulk_modulus": np.empty(n_matrices),
 75        "shear_modulus": np.empty(n_matrices),
 76        "percent_anisotropy": np.empty(n_matrices),
 77        "percent_hexagonal": np.empty(n_matrices),
 78        "percent_tetragonal": np.empty(n_matrices),
 79        "percent_orthorhombic": np.empty(n_matrices),
 80        "percent_monoclinic": np.empty(n_matrices),
 81        "percent_triclinic": np.empty(n_matrices),
 82        "hexagonal_axis": np.empty((n_matrices, 3)),
 83    }
 84    for m, matrix in enumerate(voigt_matrices):
 85        voigt_matrix = _tensors.upper_tri_to_symmetric(matrix)
 86        stiffness_dilat, stiffness_deviat = _tensors.voigt_decompose(voigt_matrix)
 87        K = np.trace(stiffness_dilat) / 9  # Bulk modulus
 88        G = (np.trace(stiffness_deviat) - 3 * K) / 10  # Shear modulus
 89        out["bulk_modulus"][m] = K
 90        out["shear_modulus"][m] = G
 91
 92        # Appendix A5, Browaeys & Chevrot, 2004.
 93        # The isotropic stiffness vector is independent of the coordinate system.
 94        isotropic_vector = np.hstack(
 95            (
 96                np.repeat(K + 4 * G / 3, 3),
 97                np.repeat(np.sqrt(2) * (K - 2 * G / 3), 3),
 98                np.repeat(2 * G, 3),
 99                np.repeat(0, 12),
100            )
101        )
102        voigt_vector = _tensors.voigt_matrix_to_vector(voigt_matrix)
103        out["percent_anisotropy"][m] = (
104            la.norm(voigt_vector - isotropic_vector) / la.norm(voigt_vector) * 100
105        )
106
107        # Search for SCCS axes (orthogonal axes defined intrinsically by
108        # eigenstrains of the elastic tensor).
109        unpermuted_SCCS = np.empty((3, 3))
110        eigv_dij = la.eigh(stiffness_dilat)[1]
111        eigv_vij = la.eigh(stiffness_deviat)[1]
112        # Averaging of eigenvectors to get the symmetry axes.
113        for i in range(3):
114            index_vij = 0  # [i(+1?)] of v_ij = eigvect that is nearest to the d_ij one.
115            angle = 10  # Initialise to any number > 2π radians.
116            for j in range(3):
117                # Calculate angle between a pair of eigenvectors.
118                # One eigenvector is taken from v_ij and one from d_ij.
119                # Do not distinguish between vectors in opposite directions.
120                dot_eigvects = np.dot(eigv_dij[:, i], eigv_vij[:, j])
121                angle_eigvects = smallest_angle(eigv_dij[:, i], eigv_vij[:, j])
122                if angle_eigvects < angle:
123                    angle = angle_eigvects
124                    # NOTE: Differences between ASPECT And original implementation:
125                    # - In C++, std::copysign(1.0, 0.0) returns -1.
126                    # - In Fortran 90, sign(1.0, 0.0) returns 1.
127                    # - ASPECT implementation uses 'j' instead of 'j+1', doesn't seem to
128                    # make a difference using the equivalent change for vec_SCSS = ...
129                    index_vij = np.sign(dot_eigvects) * j if dot_eigvects != 0 else j
130                    # index_vij = (
131                    #     np.sign(dot_eigvects) * (j + 1)
132                    #     if dot_eigvects != 0
133                    #     else (j + 1)
134                    # )
135
136            # Add/subtract the nearest v_ij eigenvector multiplied by the signed index,
137            # which effectively rotates the d_ij eigenvector, and then normalise it.
138            vec_SCCS = (
139                eigv_dij[:, i] + index_vij * eigv_vij[:, int(abs(index_vij))]
140            ) / 2
141            # vec_SCSS = (
142            #     eigv_dij[:, i] + index_vij * eigv_vij[:, int(abs(index_vij) - 1)]
143            # ) / 2
144            vec_SCCS /= la.norm(vec_SCCS)
145            unpermuted_SCCS[:, i] = vec_SCCS
146
147        # Determine SCCS permutation that gives best hexagonal approximation.
148        # This is achieved by minimising the "distance" between the voigt vector and its
149        # projection onto the hexagonal symmetry subspace.
150        # The elastic tensor is rotated into this system for decomposition.
151        elastic_tensor = _tensors.voigt_to_elastic_tensor(voigt_matrix)
152        distance = la.norm(voigt_vector)
153        for i in range(3):
154            permuted_SCCS = unpermuted_SCCS[:, [(i + j) % 3 for j in range(3)]]
155            rotated_voigt_matrix = _tensors.elastic_tensor_to_voigt(
156                _tensors.rotate(elastic_tensor, permuted_SCCS.transpose())
157            )
158            rotated_voigt_vector = _tensors.voigt_matrix_to_vector(rotated_voigt_matrix)
159            mono_and_higher_vector = _tensors.mono_project(rotated_voigt_vector)
160            tric_vector = rotated_voigt_vector - mono_and_higher_vector
161            ortho_and_higher_vector = _tensors.ortho_project(mono_and_higher_vector)
162            tetr_and_higher_vector = _tensors.tetr_project(ortho_and_higher_vector)
163            hex_and_higher_vector = _tensors.hex_project(tetr_and_higher_vector)
164            mono_vector = mono_and_higher_vector - ortho_and_higher_vector
165            ortho_vector = ortho_and_higher_vector - tetr_and_higher_vector
166            tetr_vector = tetr_and_higher_vector - hex_and_higher_vector
167            hex_vector = hex_and_higher_vector - isotropic_vector
168
169            δ = la.norm(rotated_voigt_vector - hex_and_higher_vector)
170            if δ < distance:
171                distance = δ
172
173                percent = 100 / la.norm(voigt_vector)
174                # print("\n", _tensors.voigt_vector_to_matrix(hex_vector))
175                out["percent_triclinic"][m] = la.norm(tric_vector) * percent
176                out["percent_monoclinic"][m] = la.norm(mono_vector) * percent
177                out["percent_orthorhombic"][m] = la.norm(ortho_vector) * percent
178                out["percent_tetragonal"][m] = la.norm(tetr_vector) * percent
179                out["percent_hexagonal"][m] = la.norm(hex_vector) * percent
180                # Last SCCS axis is always the hexagonal symmetry axis.
181                out["hexagonal_axis"][m, ...] = permuted_SCCS[:, 2]
182    return out
183
184
185def bingham_average(orientations: np.ndarray, axis: str = "a"):
186    """Compute Bingham average of orientation matrices.
187
188    Returns the antipodally symmetric average orientation
189    of the given crystallographic `axis`, or the a-axis by default.
190    Valid axis specifiers are "a" for [100], "b" for [010] and "c" for [001].
191
192    See also: [Watson (1966)](https://doi.org/10.1086%2F627211),
193    [Mardia & Jupp, “Directional Statistics”](https://doi.org/10.1002/9780470316979).
194
195    """
196    match axis:
197        case "a":
198            row = 0
199        case "b":
200            row = 1
201        case "c":
202            row = 2
203        case _:
204            raise ValueError(f"axis must be 'a', 'b', or 'c', not {axis}")
205
206    # https://courses.eas.ualberta.ca/eas421/lecturepages/orientation.html
207    # Eigenvector corresponding to largest eigenvalue is the mean direction.
208    # SciPy returns eigenvalues in ascending order (same order for vectors).
209    # SciPy uses lower triangular entries by default, no need for all components.
210    mean_vector = la.eigh(_stats._scatter_matrix(np.asarray(orientations), row))[1][
211        :, -1
212    ]
213    return mean_vector / la.norm(mean_vector)
214
215
216def finite_strain(deformation_gradient: np.ndarray, driver="ev"):
217    """Extract measures of finite strain from the deformation gradient.
218
219    Extracts the largest principal strain value and the vector defining the axis of
220    maximum extension (longest semiaxis of the finite strain ellipsoid) from the 3x3
221    deformation gradient tensor.
222
223    """
224    # Get eigenvalues and eigenvectors of the left stretch (Cauchy-Green) tensor.
225    B_λ, B_v = la.eigh(
226        deformation_gradient @ deformation_gradient.transpose(),
227        driver=driver,
228    )
229    # Stretch ratio is sqrt(λ) - 1, the (-1) is so that undeformed strain is 0 not 1.
230    return np.sqrt(B_λ[-1]) - 1, B_v[:, -1]
231
232
233def symmetry_pgr(orientations: np.ndarray, axis="a"):
234    r"""Compute texture symmetry eigenvalue diagnostics from grain orientation matrices.
235
236    Compute Point, Girdle and Random symmetry diagnostics
237    for ternary texture classification.
238    Returns the tuple (P, G, R) where
239    $$
240    \begin{align\*}
241    P &= (λ_{1} - λ_{2}) / N \cr
242    G &= 2 (λ_{2} - λ_{3}) / N \cr
243    R &= 3 λ_{3} / N
244    \end{align\*}
245    $$
246    with $N$ the sum of the eigenvalues $λ_{1} ≥ λ_{2} ≥ λ_{3}$
247    of the scatter (inertia) matrix.
248
249    See e.g. [Vollmer (1990)](https://doi.org/10.1130/0016-7606(1990)102%3C0786:AAOEMT%3E2.3.CO;2).
250
251    """
252    match axis:
253        case "a":
254            row = 0
255        case "b":
256            row = 1
257        case "c":
258            row = 2
259        case _:
260            raise ValueError(f"axis must be 'a', 'b', or 'c', not {axis}")
261
262    scatter = _stats._scatter_matrix(orientations, row)
263    # SciPy uses lower triangular entries by default, no need for all components.
264    eigvals_descending = la.eigvalsh(scatter)[::-1]
265    sum_eigvals = np.sum(eigvals_descending)
266    return (
267        (eigvals_descending[0] - eigvals_descending[1]) / sum_eigvals,
268        2 * (eigvals_descending[1] - eigvals_descending[2]) / sum_eigvals,
269        3 * eigvals_descending[2] / sum_eigvals,
270    )
271
272
273def misorientation_indices(
274    orientation_stack: np.ndarray,
275    system: _geo.LatticeSystem,
276    bins: int | None = None,
277    ncpus: int | None = None,
278    pool=None,
279):
280    """Calculate M-indices for a series of polycrystal textures.
281
282    Calculate M-index using `misorientation_index` for a series of texture snapshots.
283    The `orientation_stack` is a NxMx3x3 array of orientations where N is the number of
284    texture snapshots and M is the number of grains.
285
286    Uses either Ray or the Python multiprocessing library to calculate texture indices
287    for multiple snapshots simultaneously. The arguments `ncpus` and `pool` are only
288    relevant to the latter option: if `ncpus` is `None` the number of CPU cores to use
289    is chosen automatically based on the maximum number available to the Python
290    interpreter, otherwise the specified number of cores is requested. Alternatively, an
291    existing instance of `multiprocessing.Pool` can be provided.
292
293    If Ray is installed, it will be automatically preferred. In this case, the number of
294    parallel workers should be set upon initialisation of the Ray cluster (which can be
295    distributed over the network).
296
297    See `misorientation_index` for documentation of the remaining arguments.
298
299    """
300    if ncpus is not None and pool is not None:
301        _log.warning("ignoring `ncpus` argument because a Pool was provided")
302    m_indices = np.empty(len(orientation_stack))
303    _run = ft.partial(
304        misorientation_index,
305        system=system,
306        bins=bins,
307    )
308    if pool is None:
309        if ncpus is None:
310            ncpus = _utils.default_ncpus()
311        with Pool(processes=ncpus) as pool:
312            for i, out in enumerate(pool.imap(_run, orientation_stack)):
313                m_indices[i] = out
314    else:
315        if HAS_RAY:
316            m_indices = np.array(
317                ray.get(
318                    [
319                        _dstr.misorientation_index.remote(
320                            ray.put(a), system=system, bins=bins
321                        )
322                        for a in orientation_stack
323                    ]
324                )
325            )
326        else:
327            for i, out in enumerate(pool.imap(_run, orientation_stack)):
328                m_indices[i] = out
329    return m_indices
330
331
332def misorientation_index(
333    orientations: np.ndarray, system: _geo.LatticeSystem, bins: int | None = None
334):
335    r"""Calculate M-index for polycrystal orientations.
336
337    The `bins` argument is passed to `numpy.histogram`.
338    If left as `None`, 1° bins will be used as recommended by the reference paper.
339    The `symmetry` argument specifies the lattice system which determines intrinsic
340    symmetry degeneracies and the maximum allowable misorientation angle.
341    See `_geo.LatticeSystem` for supported systems.
342
343    .. warning::
344        This method must be able to allocate an array of shape
345        $ \frac{N!}{2(N-2)!}× M^{2} $
346        for N the length of `orientations` and M the number of symmetry operations for
347        the given `system`.
348
349    See [Skemer et al. (2005)](https://doi.org/10.1016/j.tecto.2005.08.023).
350
351    """
352    θmax = _stats._max_misorientation(system)
353    misorientations_count, bin_edges = _stats.misorientation_hist(
354        orientations, system, bins
355    )
356    # Compute counts of theoretical misorientation for an isotropic aggregate,
357    # using the same bins (Skemer 2005 recommend 1° bins).
358    misorientations_theory = np.array(
359        [
360            _stats.misorientations_random(bin_edges[i], bin_edges[i + 1], system)
361            for i in range(len(misorientations_count))
362        ]
363    )
364    # Equation 2 in Skemer 2005.
365    return (θmax / (2 * len(misorientations_count))) * np.sum(
366        np.abs(misorientations_theory - misorientations_count)
367    )
368
369
370def coaxial_index(orientations: np.ndarray, axis1: str = "b", axis2: str = "a"):
371    r"""Calculate coaxial “BA” index for a combination of two crystal axes.
372
373    The BA index of [Mainprice et al. (2015)](https://doi.org/10.1144/SP409.8)
374    is derived from the eigenvalue `symmetry` diagnostics and measures point vs girdle
375    symmetry in the aggregate. $BA ∈ [0, 1]$ where $BA = 0$ corresponds to a perfect
376    axial girdle texture and $BA = 1$ represents a point symmetry texture assuming that
377    the random component $R$ is negligible. May be less susceptible to random
378    fluctuations compared to the raw eigenvalue diagnostics.
379
380    """
381    P1, G1, _ = symmetry_pgr(orientations, axis=axis1)
382    P2, G2, _ = symmetry_pgr(orientations, axis=axis2)
383    return 0.5 * (2 - (P1 / (G1 + P1)) - (G2 / (G2 + P2)))
384
385
386@nb.njit(fastmath=True)
387def smallest_angle(
388    vector: np.ndarray, axis: np.ndarray, plane: np.ndarray | None = None
389):
390    """Get smallest angle between a unit `vector` and a bidirectional `axis`.
391
392    The axis is specified using either of its two parallel unit vectors.
393    Optionally project the vector onto the `plane` (given by its unit normal)
394    before calculating the angle.
395
396    Examples:
397
398    >>> import numpy as np
399    >>> def f64(x): return np.asarray(x, dtype=np.float64)
400    >>> int(smallest_angle(f64([1e0, 0e0, 0e0]), f64([1e0, 0e0, 0e0])))
401    0
402    >>> int(smallest_angle(f64([1e0, 0e0, 0e0]), f64([0e0, 1e0, 0e0])))
403    90
404    >>> int(smallest_angle(f64([1e0, 0e0, 0e0]), f64([0e0, -1e0, 0e0])))
405    90
406    >>> int(smallest_angle(f64([1e0, 0e0, 0e0]), f64([np.sqrt(2), np.sqrt(2), 0e0])))
407    45
408    >>> int(smallest_angle(f64([1e0, 0e0, 0e0]), f64([-np.sqrt(2), np.sqrt(2), 0e0])))
409    45
410
411    """
412    if plane is not None:
413        _vector = vector - plane * np.dot(vector, plane)
414    else:
415        _vector = vector
416    angle = np.rad2deg(
417        np.arccos(
418            np.clip(
419                np.asarray(  # https://github.com/numba/numba/issues/3469
420                    np.dot(_vector, axis)
421                    / (np.linalg.norm(_vector) * np.linalg.norm(axis))
422                ),
423                -1,
424                1,
425            )
426        )
427    )
428    if angle > 90:
429        return 180 - angle
430    return angle