diff --git a/src/fatpy/struct_mech/strain.py b/src/fatpy/struct_mech/strain.py index 851b1b4..cf4dac2 100644 --- a/src/fatpy/struct_mech/strain.py +++ b/src/fatpy/struct_mech/strain.py @@ -1,4 +1,292 @@ -"""Strain analysis methods of structural mechanics. +"""Calculate fundamental strain metrics and invariants. -A group of methods evaluating equivalent strain based on full-tensor input. +These functions provide principal strains, hydrostatic strain, von Mises equivalent +strain, and invariants of the strain tensor. They are essential for strength, fatigue, +and fracture analyses under both uniaxial and multiaxial loading conditions. + +Conventions: +- Vectors use Voigt notation with shape (..., 6), where the last dimension + contains the six Voigt components and leading dimensions are preserved: + + (ε_11, ε_22, ε_33, ε_23, ε_13, ε_12) + (ε_xx, ε_yy, ε_zz, ε_yz, ε_xz, ε_xy) + +- Principal strains are ordered in descending order throughout the module: + ε1 ≥ ε2 ≥ ε3. +- Principal directions (eigenvectors) are aligned to this ordering + (columns correspond to ε1, ε2, ε3). """ + +import numpy as np +from numpy.typing import NDArray + +from fatpy.utils import voigt + + +def calc_principal_strains_and_directions( + strain_vector_voigt: NDArray[np.float64], +) -> tuple[NDArray[np.float64], NDArray[np.float64]]: + r"""Calculate principal strains and principal directions for each state. + + ??? abstract "Math Equations" + Principal strains and directions are found by solving the eigenvalue problem + for the strain tensor: + + $$ \varepsilon \mathbf{v} = \lambda \mathbf{v} $$ + + Args: + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt strain components. Leading dimensions are preserved. + + Returns: + Tuple (eigvals, eigvecs): + - eigvals: Array of shape (..., 3). Principal strains + (descending: ε_1 ≥ ε_2 ≥ ε_3) with leading dimensions preserved. + - eigvecs: Array of shape (..., 3, 3). Principal directions (columns are + eigenvectors) aligned with eigvals in the same order. The last two + dimensions are the 3x3 eigenvector matrix for each input. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(strain_vector_voigt) + + tensor = voigt.voigt_to_tensor(strain_vector_voigt) + eigvals, eigvecs = np.linalg.eigh(tensor) + sorted_indices = np.argsort(eigvals, axis=-1)[..., ::-1] + eigvals_sorted = np.take_along_axis(eigvals, sorted_indices, axis=-1) + eigvecs_sorted = np.take_along_axis( + eigvecs, np.expand_dims(sorted_indices, axis=-2), axis=-1 + ) + + return eigvals_sorted, eigvecs_sorted + + +def calc_principal_strains( + strain_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + r"""Calculate principal strains for each strain state. + + ??? abstract "Math Equations" + Principal strains are found by solving the eigenvalue problem + for the strain tensor: + + $$ \varepsilon \mathbf{v} = \lambda \mathbf{v} $$ + + Args: + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt strain components. Leading dimensions are preserved. + + Returns: + Array of shape (..., 3). Principal strains (descending: ε1 ≥ ε2 ≥ ε3). + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(strain_vector_voigt) + + tensor = voigt.voigt_to_tensor(strain_vector_voigt) + eigvals = np.linalg.eigvalsh(tensor) + + return np.sort(eigvals, axis=-1)[..., ::-1] + + +def calc_strain_invariants( + strain_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + r"""Calculate the first, second, and third invariants for each strain state. + + ??? abstract "Math Equations" + $$ + \begin{align*} + I_1 &= tr(\varepsilon), \\ + I_2 &= \tfrac{1}{2}\big(I_1^{2} - tr(\varepsilon^{2})\big), \\ + I_3 &= \det(\varepsilon) + \end{align*} + $$ + + Args: + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt strain components. Leading dimensions are preserved. + + Returns: + Array of shape (..., 3). The last dimension contains (I1, I2, I3) for + each entry. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(strain_vector_voigt) + + tensor = voigt.voigt_to_tensor(strain_vector_voigt) + invariant_1 = np.trace(tensor, axis1=-2, axis2=-1) + invariant_2 = 0.5 * ( + invariant_1**2 - np.trace(np.matmul(tensor, tensor), axis1=-2, axis2=-1) + ) + invariant_3 = np.linalg.det(tensor) + + return np.stack((invariant_1, invariant_2, invariant_3), axis=-1) + + +def calc_volumetric_strain( + strain_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + r"""Calculate the volumetric (mean normal)strain for each strain state. + + ??? abstract "Math Equations" + $$ \varepsilon_{vol} = \frac{1}{3} \, tr(\varepsilon) = + \frac{1}{3}(\varepsilon_{11} + \varepsilon_{22} + \varepsilon_{33}) + $$ + + Args: + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt strain components. Leading dimensions are preserved. + + Returns: + Array of shape (...). Volumetric (mean normal) strain for each input state. + Tensor rank is reduced by one. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(strain_vector_voigt) + + return ( + strain_vector_voigt[..., 0] + + strain_vector_voigt[..., 1] + + strain_vector_voigt[..., 2] + ) / 3.0 + + +def calc_deviatoric_strain( + strain_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + r"""Calculate the deviatoric strain for each strain state. + + ??? abstract "Math Equations" + The strain tensor decomposes as: + + $$ \varepsilon = \varepsilon_{dev} + \varepsilon_{vol} \mathbf{I} $$ + + where the deviatoric part is traceless and obtained by subtracting the + volumetric part from the normal components. + + $$ \varepsilon_{dev} = \varepsilon - \frac{1}{3} tr(\varepsilon) $$ + + Args: + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt strain components. Leading dimensions are preserved. + + Returns: + Array of shape (..., 6). Deviatoric strain for each input state. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(strain_vector_voigt) + + volumetric = calc_volumetric_strain(strain_vector_voigt) + deviatoric = strain_vector_voigt.copy() + deviatoric[..., :3] = deviatoric[..., :3] - volumetric[..., None] + + return deviatoric + + +# Von Mises functions +def calc_von_mises_strain( + strain_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + r"""Von Mises equivalent strain computed directly from Voigt components. + + ??? abstract "Math Equations" + $$ + \varepsilon_{vM} = \tfrac{\sqrt{2}}{3}\sqrt{ + (\varepsilon_{11}-\varepsilon_{22})^2 + +(\varepsilon_{22}-\varepsilon_{33})^2 + +(\varepsilon_{33}-\varepsilon_{11})^2 + + 6(\varepsilon_{12}^2+\varepsilon_{23}^2+\varepsilon_{13}^2)} + $$ + + Args: + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt strain components. Leading dimensions are preserved. + + Returns: + Array of shape (...). Von Mises equivalent strain for each entry. + Tensor rank is reduced by one. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(strain_vector_voigt) + + e11 = strain_vector_voigt[..., 0] + e22 = strain_vector_voigt[..., 1] + e33 = strain_vector_voigt[..., 2] + e23 = strain_vector_voigt[..., 3] # epsilon_23 + e13 = strain_vector_voigt[..., 4] # epsilon_13 + e12 = strain_vector_voigt[..., 5] # epsilon_12 + return np.sqrt( + (2.0 / 9.0) + * ( + (e11 - e22) ** 2 + + (e22 - e33) ** 2 + + (e33 - e11) ** 2 + + 6.0 * (e12**2 + e23**2 + e13**2) + ) + ) + + +def calc_signed_von_mises_by_max_abs_principal( + strain_vector_voigt: NDArray[np.float64], + rtol: float = 1e-5, + atol: float = 1e-8, +) -> NDArray[np.float64]: + r"""Calculate signed von Mises equivalent strain for each strain state. + + Sign is determined by average of the maximum and minimum principal strains. + + ??? note "Sign Convention" + The sign assignment follows these rules: + + - **Positive (+)**: When (ε₁ + ε₃)/2 > 0 (tension dominant) + - **Negative (-)**: When (ε₁ + ε₃)/2 < 0 (compression dominant) + - **Positive (+)**: When (ε₁ + ε₃)/2 ≈ 0 (within tolerance, default fallback) + + Tolerance parameters ensure numerical stability in edge cases where the + determining value is very close to zero, preventing erratic sign changes + that could occur due to floating-point precision limitations. + + ??? abstract "Math Equations" + $$ + \varepsilon_{SvM} = \begin{cases} + +\varepsilon_{vM} & \text{if } \frac{\varepsilon_1 + \varepsilon_3}{2} \geq 0 \\ + -\varepsilon_{vM} & \text{if } \frac{\varepsilon_1 + \varepsilon_3}{2} < 0 + \end{cases} + $$ + + Args: + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt strain components. Leading dimensions are preserved. + rtol: Relative tolerance for comparing the average of maximum and minimum + principal strain to zero. Default is 1e-5. + atol: Absolute tolerance for comparing the average of maximum and minimum + principal strain to zero. Default is 1e-8. + + Returns: + Array of shape (...). Signed von Mises equivalent strain for each entry. + Tensor rank is reduced by one. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(strain_vector_voigt) + + von_mises = calc_von_mises_strain(strain_vector_voigt) + principals = calc_principal_strains(strain_vector_voigt) + + avg_13 = 0.5 * (principals[..., 0] + principals[..., 2]) + sign = np.sign(avg_13).astype(np.float64, copy=False) + sign = np.where(np.isclose(avg_13, 0, rtol=rtol, atol=atol), 1.0, sign) + + return sign * von_mises diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index 5f5f702..6d0df81 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -43,11 +43,11 @@ def calc_principal_stresses_and_directions( Returns: Tuple (eigvals, eigvecs): - - eigvals: Array of shape (..., 3). Principal stresses - (descending: σ_1 ≥ σ_2 ≥ σ_3) with leading dimensions preserved. - - eigvecs: Array of shape (..., 3, 3). Principal directions (columns are - eigenvectors) aligned with eigvals in the same order. The last two - dimensions are the 3x3 eigenvector matrix for each input. + - eigvals: Array of shape (..., 3). Principal stresses + (descending: σ_1 ≥ σ_2 ≥ σ_3) with leading dimensions preserved. + - eigvecs: Array of shape (..., 3, 3). Principal directions (columns are + eigenvectors) aligned with eigvals in the same order. The last two + dimensions are the 3x3 eigenvector matrix for each input. Raises: ValueError: If the last dimension is not of size 6. @@ -170,6 +170,13 @@ def calc_stress_deviator( r"""Calculate stress deviator for each stress state. ??? abstract "Math Equations" + The stress tensor decomposes as: + + $$ \sigma = \mathbf{s} + \sigma_H \mathbf{I} $$ + + where the deviatoric part $\mathbf{s}$ is traceless and obtained by subtracting + the hydrostatic part from the normal components. + $$ \mathbf{s} = \sigma - \frac{1}{3} tr(\sigma) $$ Args: @@ -248,10 +255,13 @@ def calc_signed_von_mises_by_hydrostatic( Sign is determined by the hydrostatic stress. - The sign assignment follows these rules: - - **Positive (+)**: When hydrostatic stress > 0 (tensile dominant state) - - **Negative (-)**: When hydrostatic stress < 0 (compressive dominant state) - - **Positive (+)**: When hydrostatic stress ≈ 0 (within tolerance, default fallback) + ??? note "Sign Convention" + The sign assignment follows these rules: + + - **Positive (+)**: When hydrostatic stress > 0 (tensile dominant state) + - **Negative (-)**: When hydrostatic stress < 0 (compressive dominant state) + - **Positive (+)**: When hydrostatic stress ≈ 0 (within tolerance, default + fallback) Tolerance parameters ensure numerical stability in edge cases where the determining value is very close to zero, preventing erratic sign changes @@ -290,7 +300,7 @@ def calc_signed_von_mises_by_hydrostatic( def calc_signed_von_mises_by_max_abs_principal( - stress_voigt: NDArray[np.float64], + stress_vector_voigt: NDArray[np.float64], rtol: float = 1e-5, atol: float = 1e-8, ) -> NDArray[np.float64]: @@ -298,10 +308,12 @@ def calc_signed_von_mises_by_max_abs_principal( Sign is determined by average of the maximum and minimum principal stresses. - The sign assignment follows these rules: - - **Positive (+)**: When (σ₁ + σ₃)/2 > 0 (tension dominant) - - **Negative (-)**: When (σ₁ + σ₃)/2 < 0 (compression dominant) - - **Positive (+)**: When (σ₁ + σ₃)/2 ≈ 0 (within tolerance, default fallback) + ??? note "Sign Convention" + The sign assignment follows these rules: + + - **Positive (+)**: When (σ₁ + σ₃)/2 > 0 (tension dominant) + - **Negative (-)**: When (σ₁ + σ₃)/2 < 0 (compression dominant) + - **Positive (+)**: When (σ₁ + σ₃)/2 ≈ 0 (within tolerance, default fallback) Tolerance parameters ensure numerical stability in edge cases where the determining value is very close to zero, preventing erratic sign changes @@ -316,7 +328,7 @@ def calc_signed_von_mises_by_max_abs_principal( $$ Args: - stress_voigt: Array of shape (..., 6). The last dimension contains the + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the Voigt stress components. Leading dimensions are preserved. rtol: Relative tolerance for comparing the average of maximum and minimum principal stresses to zero. Default is 1e-5. @@ -330,8 +342,8 @@ def calc_signed_von_mises_by_max_abs_principal( Raises: ValueError: If the last dimension is not of size 6. """ - von_mises = calc_von_mises_stress(stress_voigt) - principals = calc_principal_stresses(stress_voigt) + von_mises = calc_von_mises_stress(stress_vector_voigt) + principals = calc_principal_stresses(stress_vector_voigt) avg_13 = 0.5 * (principals[..., 0] + principals[..., 2]) sign = np.sign(avg_13).astype(np.float64, copy=False) @@ -349,10 +361,12 @@ def calc_signed_von_mises_by_first_invariant( Sign is determined by the first invariant of the stress tensor. - The sign assignment follows these rules: - - **Positive (+)**: When tr(σ) > 0 (tensile volumetric stress) - - **Negative (-)**: When tr(σ) < 0 (compressive volumetric stress) - - **Positive (+)**: When tr(σ) ≈ 0 (within tolerance, default fallback) + ??? note "Sign Convention" + The sign assignment follows these rules: + + - **Positive (+)**: When tr(σ) > 0 (tensile volumetric stress) + - **Negative (-)**: When tr(σ) < 0 (compressive volumetric stress) + - **Positive (+)**: When tr(σ) ≈ 0 (within tolerance, default fallback) Tolerance parameters ensure numerical stability in edge cases where the determining value is very close to zero, preventing erratic sign changes @@ -426,10 +440,13 @@ def calc_signed_tresca_by_hydrostatic( Sign is determined by the hydrostatic stress. - The sign assignment follows these rules: - - **Positive (+)**: When hydrostatic stress > 0 (tensile dominant state) - - **Negative (-)**: When hydrostatic stress < 0 (compressive dominant state) - - **Positive (+)**: When hydrostatic stress ≈ 0 (within tolerance, default fallback) + ??? note "Sign Convention" + The sign assignment follows these rules: + + - **Positive (+)**: When hydrostatic stress > 0 (tensile dominant state) + - **Negative (-)**: When hydrostatic stress < 0 (compressive dominant state) + - **Positive (+)**: When hydrostatic stress ≈ 0 (within tolerance, default + fallback) Tolerance parameters ensure numerical stability in edge cases where the determining value is very close to zero, preventing erratic sign changes @@ -476,10 +493,12 @@ def calc_signed_tresca_by_max_abs_principal( Sign is determined by the maximum absolute principal stress value. - The sign assignment follows these rules: - - **Positive (+)**: When (σ₁ + σ₃)/2 > 0 (tension dominant) - - **Negative (-)**: When (σ₁ + σ₃)/2 < 0 (compression dominant) - - **Positive (+)**: When (σ₁ + σ₃)/2 ≈ 0 (within tolerance, default fallback) + ??? note "Sign Convention" + The sign assignment follows these rules: + + - **Positive (+)**: When (σ₁ + σ₃)/2 > 0 (tension dominant) + - **Negative (-)**: When (σ₁ + σ₃)/2 < 0 (compression dominant) + - **Positive (+)**: When (σ₁ + σ₃)/2 ≈ 0 (within tolerance, default fallback) Tolerance parameters ensure numerical stability in edge cases where the determining value is very close to zero, preventing erratic sign changes diff --git a/tests/struct_mech/test_strain.py b/tests/struct_mech/test_strain.py new file mode 100644 index 0000000..b0417a0 --- /dev/null +++ b/tests/struct_mech/test_strain.py @@ -0,0 +1,240 @@ +"""Test functions for strain calculations in structural mechanics. + +Voigt notation: + Voigt notation tests are not repeated here as they are covered in related tests. + +Eigenvalues and eigenvectors: + + Only features related to eigenvalue/vector modification from `np.linalg.eig` are + tested as testing of eigenvalue/vector calculation itself is out of scope. + Following features are tested: + + 1. Shape of outputs + 2. Ordering of principal strains (descending) + +""" + +import numpy as np +import pytest +from numpy.typing import NDArray + +from fatpy.struct_mech import strain +from fatpy.utils import voigt + + +@pytest.fixture +def strain_vector_sample() -> NDArray[np.float64]: + """Fixture providing sample strain vectors in Voigt notation (3D: shape (2, 3, 6)). + + Returns: + NDArray[np.float64]: Sample strain vectors. + """ + # Two sets, each with three representative strain states (3,6): + arr = np.array( + [ + [ + [0.01, 0.0, 0.0, 0.0, 0.0, 0.0], # Uniaxial strain in x + [0.0, 0.0, -0.005, 0.0, 0.0, 0.0], # Uniaxial compression in z + [0.02, 0.02, 0.02, 0.0, 0.0, 0.0], # Pure volumetric + ], + [ + [0.0, 0.0, 0.0, 0.0, 0.0, 0.03], # Pure shear + [0.008, 0.008, 0.0, 0.0, 0.0, 0.006], # Mixed state + [0.014, 0.0, 0.006, 0.0, 0.003, 0.0], # Another mixed + ], + ], + dtype=np.float64, + ) + + return arr + + +@pytest.fixture +def principal_strains_sample() -> NDArray[np.float64]: + """Fixture providing expected principal strains for the sample strain vectors. + + Returns: + NDArray[np.float64]: Expected principal strains. + """ + arr = np.array( + [ + [ + [0.01, 0.0, 0.0], + [0.0, 0.0, -0.005], + [0.02, 0.02, 0.02], + ], + [ + [0.03, 0.0, -0.03], + [0.014, 0.002, 0.0], + [0.015, 0.005, 0.0], + ], + ], + dtype=np.float64, + ) + + return arr + + +@pytest.fixture +def volumetric_strain_sample() -> NDArray[np.float64]: + """Fixture providing expected volumetric strains for the sample strain vectors. + + Returns: + NDArray[np.float64]: Expected volumetric strains. + """ + arr = np.array( + [ + [(0.01 + 0.0 + 0.0) / 3.0, (0.0 + 0.0 + -0.005) / 3.0, 0.02], + [0.0, (0.008 + 0.008 + 0.0) / 3.0, (0.014 + 0.0 + 0.006) / 3.0], + ], + dtype=np.float64, + ) + + return arr + + +def test_calc_principal_strains_and_directions_shape( + strain_vector_sample: NDArray[np.float64], +) -> None: + """Test shape of principal strains and directions output.""" + principals, directions = strain.calc_principal_strains_and_directions( + strain_vector_sample + ) + assert principals.shape[:-1] == strain_vector_sample.shape[:-1] + assert principals.shape[-1] == 3 + + assert directions.shape[:-2] == strain_vector_sample.shape[:-1] + assert directions.shape[-2:] == (3, 3) # 3x3 eigenvector matrix + + +def test_calc_principal_strains_and_directions_ordering( + strain_vector_sample: NDArray[np.float64], +) -> None: + """Test ordering of principal strains (descending).""" + principals, directions = strain.calc_principal_strains_and_directions( + strain_vector_sample + ) + assert np.all(principals[..., 0] >= principals[..., 1]) + assert np.all(principals[..., 1] >= principals[..., 2]) + + # test if directions order matches principal strains order + # by checking A*v = λ*v for each principal strain/direction pair + for idx in np.ndindex(principals.shape[:-1]): + strain_tensor = voigt.voigt_to_tensor(strain_vector_sample[idx]) + for i in range(3): + principal_strain = principals[idx + (i,)] + direction = directions[idx + (slice(None), i)] + Av = strain_tensor @ direction + lv = principal_strain * direction + assert np.allclose(Av, lv, atol=1e-12) + + +def test_calc_principal_strains( + strain_vector_sample: NDArray[np.float64], + principal_strains_sample: NDArray[np.float64], +) -> None: + """Test shape and values of principal strains calculation.""" + principals = strain.calc_principal_strains(strain_vector_sample) + + assert principals.shape[:-1] == strain_vector_sample.shape[:-1] + assert principals.shape[-1] == 3 + + assert np.all(principals[..., 0] >= principals[..., 1]) + assert np.all(principals[..., 1] >= principals[..., 2]) + + assert np.allclose(principals, principal_strains_sample, atol=1e-12) + + +def test_calc_strain_invariants( + strain_vector_sample: NDArray[np.float64], +) -> None: + """Test calculation of strain invariants.""" + invariants = strain.calc_strain_invariants(strain_vector_sample) + + assert invariants.shape == strain_vector_sample.shape[:-1] + (3,) + + for idx in np.ndindex(invariants.shape[:-1]): + strain_tensor = voigt.voigt_to_tensor(strain_vector_sample[idx]) + + I1 = np.trace(strain_tensor) + I2 = 0.5 * (I1**2 - np.trace(strain_tensor @ strain_tensor)) + I3 = np.linalg.det(strain_tensor) + + assert np.isclose(invariants[idx + (0,)], I1, atol=1e-12) + assert np.isclose(invariants[idx + (1,)], I2, atol=1e-12) + assert np.isclose(invariants[idx + (2,)], I3, atol=1e-12) + + +def test_calc_volumetric_strain( + strain_vector_sample: NDArray[np.float64], + volumetric_strain_sample: NDArray[np.float64], +) -> None: + """Test calculation of volumetric strain.""" + volumetric = strain.calc_volumetric_strain(strain_vector_sample) + + assert volumetric.shape == strain_vector_sample.shape[:-1] + + assert np.allclose(volumetric, volumetric_strain_sample, atol=1e-12) + + +def test_calc_deviatoric_strain( + strain_vector_sample: NDArray[np.float64], +) -> None: + """Test calculation of deviatoric strain.""" + deviator = strain.calc_deviatoric_strain(strain_vector_sample) + + assert deviator.shape == strain_vector_sample.shape + + for idx in np.ndindex(deviator.shape[:-1]): + strain_tensor = voigt.voigt_to_tensor(strain_vector_sample[idx]) + hydro = np.trace(strain_tensor) / 3.0 + hydro_tensor = np.eye(3) * hydro + deviator_tensor = strain_tensor - hydro_tensor + deviator_voigt = voigt.tensor_to_voigt(deviator_tensor) + assert np.allclose(deviator[idx], deviator_voigt, atol=1e-12) + + +def test_calc_von_mises( + strain_vector_sample: NDArray[np.float64], +) -> None: + """Test calculation of von Mises strain from Voigt components. + + Uses the principals method for verification. + + Args: + strain_vector_sample (NDArray[np.float64]): Sample strain vectors. + """ + von_mises = strain.calc_von_mises_strain(strain_vector_sample) + + assert von_mises.shape == strain_vector_sample.shape[:-1] + + strain_tensor = voigt.voigt_to_tensor(strain_vector_sample) + principals = np.linalg.eigvalsh(strain_tensor) + e1 = principals[..., 0] + e2 = principals[..., 1] + e3 = principals[..., 2] + vm = np.sqrt((2.0 / 9.0) * ((e1 - e2) ** 2 + (e2 - e3) ** 2 + (e3 - e1) ** 2)) + assert np.allclose(von_mises, vm, atol=1e-12) + + +def test_calc_signed_von_mises_by_max_abs_principal( + strain_vector_sample: NDArray[np.float64], +) -> None: + """Test calculation of signed von Mises strain. + Utilizes average of maximum and minimum principal strains for sign determination. + """ + signed_von_mises = strain.calc_signed_von_mises_by_max_abs_principal( + strain_vector_sample + ) + + assert signed_von_mises.shape == strain_vector_sample.shape[:-1] + + for idx in np.ndindex(signed_von_mises.shape): + strain_tensor = voigt.voigt_to_tensor(strain_vector_sample[idx]) + s = strain_tensor - np.eye(3) * (np.trace(strain_tensor) / 3.0) + vm = np.sqrt((2.0 / 3.0) * np.sum(s**2)) + principals = np.linalg.eigvalsh(strain_tensor) + avg13 = 0.5 * (principals[0] + principals[2]) + sign = np.sign(avg13) + sign = np.where(np.isclose(avg13, 0.0), 1.0, sign) + assert np.isclose(signed_von_mises[idx], sign * vm, atol=1e-12)