"""Strain → stress + stress-invariant helpers (T1-P5).
Isotropic linear-elastic stress recovery:
* :func:`stress_from_strain` — contracts an engineering-shear Voigt
strain ``(..., 6)`` with an elasticity matrix ``C (6, 6)`` to produce
Voigt stress ``(..., 6)``. The leading shape is preserved, so the
caller can pass a single point ``(6,)``, one element ``(n_nodes, 6)``,
or a stack ``(n_elems, n_nodes, 6)``.
* :func:`stress_invariants` — derives the scalar fields engineers
actually plot (von Mises, principal σ₁/σ₂/σ₃, max-shear, hydrostatic,
stress intensity) from Voigt stress.
Voigt convention matches the strain kernels and
:meth:`femorph_solver.result._material.MaterialTable.eval_c_isotropic`:
``[σxx, σyy, σzz, τxy, τyz, τxz]`` with engineering shears.
References
----------
* Hooke's law in Voigt form (σ = C·ε) for isotropic linear
elasticity: Cook, Malkus, Plesha, Witt, *Concepts and
Applications of Finite Element Analysis*, 4th ed., Wiley, 2002,
§3.7 (stress--strain relation, engineering-shear convention).
* Stress invariants (principal stresses, von Mises equivalent
``σ_vm``, stress intensity ``σ_int``, hydrostatic ``σ_h``,
max-shear ``τ_max``): Cook et al. §7.5. Principal-stress
cubic from the characteristic equation of the 3×3 stress
tensor; von Mises ``σ_vm = √(1.5 s:s)`` with deviatoric
``s = σ − σ_h·I``. Continuum-mechanics reference:
Timoshenko, S.P. and Goodier, J.N., *Theory of Elasticity*,
3rd ed., McGraw-Hill, 1970, §4.
"""
from __future__ import annotations
import numpy as np
__all__ = [
"stress_from_strain",
"stress_invariants",
]
def stress_from_strain(strain: np.ndarray, c_matrix: np.ndarray) -> np.ndarray:
"""Contract ``C (6, 6)`` with a Voigt strain to produce Voigt stress.
Parameters
----------
strain : numpy.ndarray
Voigt strain, shape ``(..., 6)``. Any leading shape is fine;
the contraction runs over the trailing axis.
c_matrix : numpy.ndarray
``(6, 6)`` elasticity matrix in Voigt order. Typically
produced by
:meth:`femorph_solver.result._material.MaterialTable.eval_c_isotropic`.
Returns
-------
numpy.ndarray
Voigt stress with the same shape as ``strain``.
"""
strain = np.asarray(strain, dtype=np.float64)
c_matrix = np.asarray(c_matrix, dtype=np.float64)
if strain.shape[-1] != 6:
raise ValueError(f"strain last axis must be 6, got shape {strain.shape}")
if c_matrix.shape != (6, 6):
raise ValueError(f"c_matrix must be (6, 6), got {c_matrix.shape}")
return np.einsum("ij,...j->...i", c_matrix, strain)
[docs]
def stress_invariants(stress: np.ndarray) -> dict[str, np.ndarray]:
"""Scalar invariants of a Voigt stress field.
Parameters
----------
stress : numpy.ndarray
``(..., 6)`` Voigt stress ``[σxx, σyy, σzz, τxy, τyz, τxz]``.
Returns
-------
dict[str, numpy.ndarray]
Each value has the same leading shape as ``stress`` minus the
trailing ``6`` axis. Keys:
* ``"von_mises"`` — :math:`\\sigma_{vm} = \\sqrt{\\tfrac{3}{2}\\,s_{ij}s_{ij}}`.
* ``"s1"``, ``"s2"``, ``"s3"`` — principal stresses, sorted
descending (``s1 >= s2 >= s3``).
* ``"max_shear"`` — :math:`(s_1 - s_3)/2`.
* ``"hydrostatic"`` — :math:`(\\sigma_{xx}+\\sigma_{yy}+\\sigma_{zz})/3`.
* ``"stress_intensity"`` — :math:`s_1 - s_3`.
"""
stress = np.asarray(stress, dtype=np.float64)
if stress.shape[-1] != 6:
raise ValueError(f"stress last axis must be 6, got shape {stress.shape}")
sxx = stress[..., 0]
syy = stress[..., 1]
szz = stress[..., 2]
txy = stress[..., 3]
tyz = stress[..., 4]
txz = stress[..., 5]
hydrostatic = (sxx + syy + szz) / 3.0
# von Mises via the deviatoric-double-dot form — numerically
# identical to the (σxx−σyy)²+... expansion but cleaner.
von_mises = np.sqrt(
0.5
* (
(sxx - syy) ** 2
+ (syy - szz) ** 2
+ (szz - sxx) ** 2
+ 6.0 * (txy**2 + tyz**2 + txz**2)
)
)
# Principal stresses via eigendecomposition of the 3×3 Cauchy
# tensor. np.linalg.eigvalsh returns ascending; reverse to the
# canonical s1 >= s2 >= s3 convention.
leading_shape = stress.shape[:-1]
flat = stress.reshape(-1, 6)
tensors = np.empty((flat.shape[0], 3, 3), dtype=np.float64)
tensors[:, 0, 0] = flat[:, 0]
tensors[:, 1, 1] = flat[:, 1]
tensors[:, 2, 2] = flat[:, 2]
tensors[:, 0, 1] = tensors[:, 1, 0] = flat[:, 3]
tensors[:, 1, 2] = tensors[:, 2, 1] = flat[:, 4]
tensors[:, 0, 2] = tensors[:, 2, 0] = flat[:, 5]
eigs = np.linalg.eigvalsh(tensors) # ascending
s3 = eigs[:, 0].reshape(leading_shape)
s2 = eigs[:, 1].reshape(leading_shape)
s1 = eigs[:, 2].reshape(leading_shape)
return {
"von_mises": von_mises,
"s1": s1,
"s2": s2,
"s3": s3,
"max_shear": 0.5 * (s1 - s3),
"hydrostatic": hydrostatic,
"stress_intensity": s1 - s3,
}