From 6438e0ba4df07aa15da29f976bd56e819144961a Mon Sep 17 00:00:00 2001 From: Tomas Karas Date: Wed, 15 Oct 2025 15:58:54 +0200 Subject: [PATCH 1/9] Strain funkcionality --- src/fatpy/struct_mech/strain.py | 297 ++++++++++++++++++++++++++++++++ 1 file changed, 297 insertions(+) create mode 100644 src/fatpy/struct_mech/strain.py diff --git a/src/fatpy/struct_mech/strain.py b/src/fatpy/struct_mech/strain.py new file mode 100644 index 0000000..f99a3ed --- /dev/null +++ b/src/fatpy/struct_mech/strain.py @@ -0,0 +1,297 @@ +"""Calculate fundamental strain metrics and invariants. + +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_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_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_voigt) + + tensor = voigt.voigt_to_tensor(strain_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_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_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_voigt) + + tensor = voigt.voigt_to_tensor(strain_voigt) + eigvals = np.linalg.eigvalsh(tensor) + + return np.sort(eigvals, axis=-1)[..., ::-1] + + +def calc_strain_invariants(strain_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_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_voigt) + + tensor = voigt.voigt_to_tensor(strain_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_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + r"""Calculate the volumetric (mean normal / hydrostatic) 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_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_voigt) + + return (strain_voigt[..., 0] + strain_voigt[..., 1] + strain_voigt[..., 2]) / 3.0 + + +def calc_deviatoric_strain(strain_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_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_voigt) + + volumetric = calc_volumetric_strain(strain_voigt) + deviatoric = strain_voigt.copy() + deviatoric[..., 0:3] = deviatoric[..., 0:3] - volumetric[..., None] + + return deviatoric + + +# Von Mises functions +def calc_von_mises_strain_from_principals( + strain_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + r"""Calculate von Mises equivalent strain for each strain state using principal strains. + + ??? abstract "Math Equations" + $$ \varepsilon_{vM} = \frac{\sqrt{2}}{3} \sqrt{ + (\varepsilon_1-\varepsilon_2)^2 + (\varepsilon_2-\varepsilon_3)^2 + + (\varepsilon_3-\varepsilon_1)^2} $$ + + Args: + strain_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_voigt) + + principals = calc_principal_strains(strain_voigt) + e1 = principals[..., 0] + e2 = principals[..., 1] + e3 = principals[..., 2] + + return np.sqrt((2 / 9.0) * ((e1 - e2) ** 2 + (e2 - e3) ** 2 + (e3 - e1) ** 2)) + + +# ? Definition https://www.sciencedirect.com/topics/engineering/equivalent-strain +def calc_von_mises_strain_voigt( + strain_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_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_voigt) + + e11 = strain_voigt[..., 0] + e22 = strain_voigt[..., 1] + e33 = strain_voigt[..., 2] + e23 = strain_voigt[..., 3] # epsilon_23 + e13 = strain_voigt[..., 4] # epsilon_13 + e12 = strain_voigt[..., 5] # epsilon_12 + + term_norm = (e11 - e22) ** 2 + (e22 - e33) ** 2 + (e33 - e11) ** 2 + term_shear = 6.0 * (e12**2 + e23**2 + e13**2) + + return np.sqrt((2 / 9.0) * term_norm + term_shear) + + +def calc_signed_von_mises_by_max_abs_principal( + strain_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 the principal strain with the largest absolute value. + + ??? note "Sign Convention" + The sign assignment follows these rules: + + - **Positive (+)**: When the max absolute principal strain > 0 (tension dominant) + - **Negative (-)**: When the max absolute principal strain < 0 (compression dominant) + - **Positive (+)**: When max absolute principal strain ≈ 0 (within tolerance, default fallback) + + Args: + strain_voigt: Array of shape (..., 6). The last dimension contains the + Voigt strain components. Leading dimensions are preserved. + rtol: Relative tolerance for comparing the maximum absolute principal strain to zero. + Default is 1e-5. + atol: Absolute tolerance for comparing the maximum absolute 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_voigt) + + von_mises = calc_von_mises_strain_voigt(strain_voigt) + principals = calc_principal_strains(strain_voigt) + + indices = np.argmax(np.abs(principals), axis=-1) + max_abs_values = np.take_along_axis(principals, indices[..., None], axis=-1)[..., 0] + sign = np.sign(max_abs_values).astype(np.float64, copy=False) + sign = np.where(np.isclose(max_abs_values, 0, rtol=rtol, atol=atol), 1.0, sign) + + return sign * von_mises From e426e832f7faf561c20fd13cc2ebfe9410677361 Mon Sep 17 00:00:00 2001 From: Tomas Karas Date: Wed, 15 Oct 2025 16:52:39 +0200 Subject: [PATCH 2/9] updated docs strings --- src/fatpy/struct_mech/strain.py | 37 ++++++++++-------- src/fatpy/struct_mech/stress.py | 69 +++++++++++++++++++++------------ 2 files changed, 65 insertions(+), 41 deletions(-) diff --git a/src/fatpy/struct_mech/strain.py b/src/fatpy/struct_mech/strain.py index f99a3ed..467981c 100644 --- a/src/fatpy/struct_mech/strain.py +++ b/src/fatpy/struct_mech/strain.py @@ -40,11 +40,11 @@ def calc_principal_strains_and_directions( 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. + 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. @@ -107,7 +107,7 @@ def calc_strain_invariants(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa Returns: Array of shape (..., 3). The last dimension contains (I1, I2, I3) for - each entry. + each entry. Raises: ValueError: If the last dimension is not of size 6. @@ -125,7 +125,7 @@ def calc_strain_invariants(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa def calc_volumetric_strain(strain_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - r"""Calculate the volumetric (mean normal / hydrostatic) strain for each strain state. + r"""Calculate the volumetric (mean normal)strain for each strain state. ??? abstract "Math Equations" $$ \varepsilon_{vol} = \frac{1}{3} \, tr(\varepsilon) = @@ -138,7 +138,7 @@ def calc_volumetric_strain(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa Returns: Array of shape (...). Volumetric (mean normal) strain for each input state. - Tensor rank is reduced by one. + Tensor rank is reduced by one. Raises: ValueError: If the last dimension is not of size 6. @@ -197,7 +197,7 @@ def calc_von_mises_strain_from_principals( Returns: Array of shape (...). Von Mises equivalent strain for each entry. - Tensor rank is reduced by one. + Tensor rank is reduced by one. Raises: ValueError: If the last dimension is not of size 6. @@ -233,7 +233,7 @@ def calc_von_mises_strain_voigt( Returns: Array of shape (...). Von Mises equivalent strain for each entry. - Tensor rank is reduced by one. + Tensor rank is reduced by one. Raises: ValueError: If the last dimension is not of size 6. @@ -265,21 +265,26 @@ def calc_signed_von_mises_by_max_abs_principal( ??? note "Sign Convention" The sign assignment follows these rules: - - **Positive (+)**: When the max absolute principal strain > 0 (tension dominant) - - **Negative (-)**: When the max absolute principal strain < 0 (compression dominant) - - **Positive (+)**: When max absolute principal strain ≈ 0 (within tolerance, default fallback) + - **Positive (+)**: When the max absolute principal strain > 0 (tension + dominant) + - **Negative (-)**: When the max absolute principal strain < 0 (compression + dominant) + - **Positive (+)**: When max absolute principal strain ≈ 0 (within tolerance, + default fallback) Args: strain_voigt: Array of shape (..., 6). The last dimension contains the Voigt strain components. Leading dimensions are preserved. - rtol: Relative tolerance for comparing the maximum absolute principal strain to zero. + rtol: Relative tolerance for comparing the maximum absolute principal strain + to zero. Default is 1e-5. - atol: Absolute tolerance for comparing the maximum absolute principal strain to zero. + atol: Absolute tolerance for comparing the maximum absolute 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. + Tensor rank is reduced by one. Raises: ValueError: If the last dimension is not of size 6. diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index 5f5f702..d32ccc1 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 @@ -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 @@ -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 From c641f82eeb08e843bf04097f0dc31a4d5d18dda5 Mon Sep 17 00:00:00 2001 From: Tomas Karas Date: Wed, 15 Oct 2025 17:12:59 +0200 Subject: [PATCH 3/9] docs fix --- src/fatpy/struct_mech/strain.py | 2 +- src/fatpy/struct_mech/stress.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/fatpy/struct_mech/strain.py b/src/fatpy/struct_mech/strain.py index 467981c..73b2850 100644 --- a/src/fatpy/struct_mech/strain.py +++ b/src/fatpy/struct_mech/strain.py @@ -184,7 +184,7 @@ def calc_deviatoric_strain(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa def calc_von_mises_strain_from_principals( strain_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: - r"""Calculate von Mises equivalent strain for each strain state using principal strains. + r"""Calculate von Mises equivalent strain for each strain state using principals. ??? abstract "Math Equations" $$ \varepsilon_{vM} = \frac{\sqrt{2}}{3} \sqrt{ diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index d32ccc1..13dba85 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -174,8 +174,8 @@ def calc_stress_deviator( $$ \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. + 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) $$ From 9400f68d85108569b91d3149b96f66964d51b64d Mon Sep 17 00:00:00 2001 From: Tomas Karas Date: Tue, 21 Oct 2025 12:07:28 +0200 Subject: [PATCH 4/9] Added tests --- src/fatpy/struct_mech/strain.py | 2 +- tests/struct_mech/test_strain.py | 255 +++++++++++++++++++++++++++++++ 2 files changed, 256 insertions(+), 1 deletion(-) create mode 100644 tests/struct_mech/test_strain.py diff --git a/src/fatpy/struct_mech/strain.py b/src/fatpy/struct_mech/strain.py index 73b2850..9c16214 100644 --- a/src/fatpy/struct_mech/strain.py +++ b/src/fatpy/struct_mech/strain.py @@ -250,7 +250,7 @@ def calc_von_mises_strain_voigt( term_norm = (e11 - e22) ** 2 + (e22 - e33) ** 2 + (e33 - e11) ** 2 term_shear = 6.0 * (e12**2 + e23**2 + e13**2) - return np.sqrt((2 / 9.0) * term_norm + term_shear) + return np.sqrt((2.0 / 9.0) * (term_norm + term_shear)) def calc_signed_von_mises_by_max_abs_principal( diff --git a/tests/struct_mech/test_strain.py b/tests/struct_mech/test_strain.py new file mode 100644 index 0000000..717bd8b --- /dev/null +++ b/tests/struct_mech/test_strain.py @@ -0,0 +1,255 @@ +"""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: + 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: + 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: + 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: + 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_from_principals( + strain_vector_sample: NDArray[np.float64], +) -> None: + """Test calculation of von Mises strain from principals. + + Uses the definition based on the strain deviator tensor. + + Args: + strain_vector_sample (NDArray[np.float64]): Sample strain vectors. + """ + von_mises = strain.calc_von_mises_strain_from_principals(strain_vector_sample) + + assert von_mises.shape == strain_vector_sample.shape[:-1] + + for idx in np.ndindex(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)) + assert np.isclose(von_mises[idx], vm, atol=1e-12) + + +def test_calc_von_mises_voigt( + 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_voigt(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: + 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) + indices = np.argmax(np.abs(principals)) + max_abs = principals[indices] + sign = np.sign(max_abs) + sign = np.where(np.isclose(max_abs, 0.0), 1.0, sign) + assert np.isclose(signed_von_mises[idx], sign * vm, atol=1e-12) From 456acfbad4c8c08300fe91d00891de5b25ff7379 Mon Sep 17 00:00:00 2001 From: Tomas Karas Date: Wed, 22 Oct 2025 17:03:54 +0200 Subject: [PATCH 5/9] sign_von_mises fix --- src/fatpy/struct_mech/strain.py | 25 ++++++++++++++++++------- tests/struct_mech/test_strain.py | 18 ++++++++++++------ 2 files changed, 30 insertions(+), 13 deletions(-) diff --git a/src/fatpy/struct_mech/strain.py b/src/fatpy/struct_mech/strain.py index 9c16214..c6b3ec6 100644 --- a/src/fatpy/struct_mech/strain.py +++ b/src/fatpy/struct_mech/strain.py @@ -213,7 +213,7 @@ def calc_von_mises_strain_from_principals( # ? Definition https://www.sciencedirect.com/topics/engineering/equivalent-strain -def calc_von_mises_strain_voigt( +def calc_von_mises_strain( strain_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: r"""Von Mises equivalent strain computed directly from Voigt components. @@ -260,7 +260,7 @@ def calc_signed_von_mises_by_max_abs_principal( ) -> NDArray[np.float64]: r"""Calculate signed von Mises equivalent strain for each strain state. - Sign is determined by the principal strain with the largest absolute value. + Sign is determined by average of the maximum and minimum principal strains. ??? note "Sign Convention" The sign assignment follows these rules: @@ -272,6 +272,18 @@ def calc_signed_von_mises_by_max_abs_principal( - **Positive (+)**: When max absolute principal strain ≈ 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_voigt: Array of shape (..., 6). The last dimension contains the Voigt strain components. Leading dimensions are preserved. @@ -291,12 +303,11 @@ def calc_signed_von_mises_by_max_abs_principal( """ voigt.check_shape(strain_voigt) - von_mises = calc_von_mises_strain_voigt(strain_voigt) + von_mises = calc_von_mises_strain(strain_voigt) principals = calc_principal_strains(strain_voigt) - indices = np.argmax(np.abs(principals), axis=-1) - max_abs_values = np.take_along_axis(principals, indices[..., None], axis=-1)[..., 0] - sign = np.sign(max_abs_values).astype(np.float64, copy=False) - sign = np.where(np.isclose(max_abs_values, 0, rtol=rtol, atol=atol), 1.0, sign) + 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/tests/struct_mech/test_strain.py b/tests/struct_mech/test_strain.py index 717bd8b..8340aa2 100644 --- a/tests/struct_mech/test_strain.py +++ b/tests/struct_mech/test_strain.py @@ -133,6 +133,7 @@ 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] @@ -147,6 +148,7 @@ def test_calc_principal_strains( 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,) @@ -167,6 +169,7 @@ 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] @@ -177,6 +180,7 @@ def test_calc_volumetric_strain( 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 @@ -211,7 +215,7 @@ def test_calc_von_mises_from_principals( assert np.isclose(von_mises[idx], vm, atol=1e-12) -def test_calc_von_mises_voigt( +def test_calc_von_mises( strain_vector_sample: NDArray[np.float64], ) -> None: """Test calculation of von Mises strain from Voigt components. @@ -221,7 +225,7 @@ def test_calc_von_mises_voigt( Args: strain_vector_sample (NDArray[np.float64]): Sample strain vectors. """ - von_mises = strain.calc_von_mises_strain_voigt(strain_vector_sample) + von_mises = strain.calc_von_mises_strain(strain_vector_sample) assert von_mises.shape == strain_vector_sample.shape[:-1] @@ -237,6 +241,9 @@ def test_calc_von_mises_voigt( 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 ) @@ -248,8 +255,7 @@ def test_calc_signed_von_mises_by_max_abs_principal( 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) - indices = np.argmax(np.abs(principals)) - max_abs = principals[indices] - sign = np.sign(max_abs) - sign = np.where(np.isclose(max_abs, 0.0), 1.0, sign) + 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) From 800bd93620188aa8d70403a3a8f07f131cbed265 Mon Sep 17 00:00:00 2001 From: Tomas Karas Date: Tue, 28 Oct 2025 14:45:10 +0100 Subject: [PATCH 6/9] fix based on review --- src/fatpy/struct_mech/strain.py | 122 ++++++++++++++++++-------------- src/fatpy/struct_mech/stress.py | 12 ++-- 2 files changed, 73 insertions(+), 61 deletions(-) diff --git a/src/fatpy/struct_mech/strain.py b/src/fatpy/struct_mech/strain.py index c6b3ec6..56551cb 100644 --- a/src/fatpy/struct_mech/strain.py +++ b/src/fatpy/struct_mech/strain.py @@ -24,7 +24,7 @@ def calc_principal_strains_and_directions( - strain_voigt: NDArray[np.float64], + strain_vector_voigt: NDArray[np.float64], ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: r"""Calculate principal strains and principal directions for each state. @@ -35,23 +35,23 @@ def calc_principal_strains_and_directions( $$ \varepsilon \mathbf{v} = \lambda \mathbf{v} $$ Args: - strain_voigt: Array of shape (..., 6). The last dimension contains the + 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 + - 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 + - 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_voigt) + voigt.check_shape(strain_vector_voigt) - tensor = voigt.voigt_to_tensor(strain_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) @@ -62,7 +62,9 @@ def calc_principal_strains_and_directions( return eigvals_sorted, eigvecs_sorted -def calc_principal_strains(strain_voigt: NDArray[np.float64]) -> NDArray[np.float64]: +def calc_principal_strains( + strain_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: r"""Calculate principal strains for each strain state. ??? abstract "Math Equations" @@ -72,7 +74,7 @@ def calc_principal_strains(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa $$ \varepsilon \mathbf{v} = \lambda \mathbf{v} $$ Args: - strain_voigt: Array of shape (..., 6). The last dimension contains the + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the Voigt strain components. Leading dimensions are preserved. Returns: @@ -81,15 +83,17 @@ def calc_principal_strains(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa Raises: ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(strain_voigt) + voigt.check_shape(strain_vector_voigt) - tensor = voigt.voigt_to_tensor(strain_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_voigt: NDArray[np.float64]) -> NDArray[np.float64]: +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" @@ -102,7 +106,7 @@ def calc_strain_invariants(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa $$ Args: - strain_voigt: Array of shape (..., 6). The last dimension contains the + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the Voigt strain components. Leading dimensions are preserved. Returns: @@ -112,9 +116,9 @@ def calc_strain_invariants(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa Raises: ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(strain_voigt) + voigt.check_shape(strain_vector_voigt) - tensor = voigt.voigt_to_tensor(strain_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) @@ -124,7 +128,9 @@ def calc_strain_invariants(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa return np.stack((invariant_1, invariant_2, invariant_3), axis=-1) -def calc_volumetric_strain(strain_voigt: NDArray[np.float64]) -> NDArray[np.float64]: +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" @@ -133,7 +139,7 @@ def calc_volumetric_strain(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa $$ Args: - strain_voigt: Array of shape (..., 6). The last dimension contains the + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the Voigt strain components. Leading dimensions are preserved. Returns: @@ -143,12 +149,18 @@ def calc_volumetric_strain(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa Raises: ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(strain_voigt) + voigt.check_shape(strain_vector_voigt) - return (strain_voigt[..., 0] + strain_voigt[..., 1] + strain_voigt[..., 2]) / 3.0 + return ( + strain_vector_voigt[..., 0] + + strain_vector_voigt[..., 1] + + strain_vector_voigt[..., 2] + ) / 3.0 -def calc_deviatoric_strain(strain_voigt: NDArray[np.float64]) -> NDArray[np.float64]: +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" @@ -162,7 +174,7 @@ def calc_deviatoric_strain(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa $$ \varepsilon_{dev} = \varepsilon - \frac{1}{3} tr(\varepsilon) $$ Args: - strain_voigt: Array of shape (..., 6). The last dimension contains the + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the Voigt strain components. Leading dimensions are preserved. Returns: @@ -171,18 +183,18 @@ def calc_deviatoric_strain(strain_voigt: NDArray[np.float64]) -> NDArray[np.floa Raises: ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(strain_voigt) + voigt.check_shape(strain_vector_voigt) - volumetric = calc_volumetric_strain(strain_voigt) - deviatoric = strain_voigt.copy() - deviatoric[..., 0:3] = deviatoric[..., 0:3] - volumetric[..., None] + 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_from_principals( - strain_voigt: NDArray[np.float64], + strain_vector_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: r"""Calculate von Mises equivalent strain for each strain state using principals. @@ -192,7 +204,7 @@ def calc_von_mises_strain_from_principals( (\varepsilon_3-\varepsilon_1)^2} $$ Args: - strain_voigt: Array of shape (..., 6). The last dimension contains the + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the Voigt strain components. Leading dimensions are preserved. Returns: @@ -202,9 +214,9 @@ def calc_von_mises_strain_from_principals( Raises: ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(strain_voigt) + voigt.check_shape(strain_vector_voigt) - principals = calc_principal_strains(strain_voigt) + principals = calc_principal_strains(strain_vector_voigt) e1 = principals[..., 0] e2 = principals[..., 1] e3 = principals[..., 2] @@ -212,9 +224,8 @@ def calc_von_mises_strain_from_principals( return np.sqrt((2 / 9.0) * ((e1 - e2) ** 2 + (e2 - e3) ** 2 + (e3 - e1) ** 2)) -# ? Definition https://www.sciencedirect.com/topics/engineering/equivalent-strain def calc_von_mises_strain( - strain_voigt: NDArray[np.float64], + strain_vector_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: r"""Von Mises equivalent strain computed directly from Voigt components. @@ -228,7 +239,7 @@ def calc_von_mises_strain( $$ Args: - strain_voigt: Array of shape (..., 6). The last dimension contains the + strain_vector_voigt: Array of shape (..., 6). The last dimension contains the Voigt strain components. Leading dimensions are preserved. Returns: @@ -238,23 +249,27 @@ def calc_von_mises_strain( Raises: ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(strain_voigt) - - e11 = strain_voigt[..., 0] - e22 = strain_voigt[..., 1] - e33 = strain_voigt[..., 2] - e23 = strain_voigt[..., 3] # epsilon_23 - e13 = strain_voigt[..., 4] # epsilon_13 - e12 = strain_voigt[..., 5] # epsilon_12 - - term_norm = (e11 - e22) ** 2 + (e22 - e33) ** 2 + (e33 - e11) ** 2 - term_shear = 6.0 * (e12**2 + e23**2 + e13**2) - - return np.sqrt((2.0 / 9.0) * (term_norm + term_shear)) + 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_voigt: NDArray[np.float64], + strain_vector_voigt: NDArray[np.float64], rtol: float = 1e-5, atol: float = 1e-8, ) -> NDArray[np.float64]: @@ -265,12 +280,9 @@ def calc_signed_von_mises_by_max_abs_principal( ??? note "Sign Convention" The sign assignment follows these rules: - - **Positive (+)**: When the max absolute principal strain > 0 (tension - dominant) - - **Negative (-)**: When the max absolute principal strain < 0 (compression - dominant) - - **Positive (+)**: When max absolute principal strain ≈ 0 (within tolerance, - default fallback) + - **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 @@ -285,7 +297,7 @@ def calc_signed_von_mises_by_max_abs_principal( $$ Args: - strain_voigt: Array of shape (..., 6). The last dimension contains the + 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 maximum absolute principal strain to zero. @@ -301,10 +313,10 @@ def calc_signed_von_mises_by_max_abs_principal( Raises: ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(strain_voigt) + voigt.check_shape(strain_vector_voigt) - von_mises = calc_von_mises_strain(strain_voigt) - principals = calc_principal_strains(strain_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) diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index 13dba85..6d0df81 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -43,9 +43,9 @@ def calc_principal_stresses_and_directions( Returns: Tuple (eigvals, eigvecs): - eigvals: Array of shape (..., 3). Principal stresses + - 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 + - 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. @@ -300,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]: @@ -328,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. @@ -342,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) From 745e5c2439ea30c82e4c811a8ce64b20a31efbba Mon Sep 17 00:00:00 2001 From: Tomas Karas Date: Tue, 28 Oct 2025 14:56:14 +0100 Subject: [PATCH 7/9] von mises func moved --- src/fatpy/struct_mech/strain.py | 31 ------------------------------- tests/struct_mech/test_strain.py | 21 --------------------- 2 files changed, 52 deletions(-) diff --git a/src/fatpy/struct_mech/strain.py b/src/fatpy/struct_mech/strain.py index 56551cb..a12db1c 100644 --- a/src/fatpy/struct_mech/strain.py +++ b/src/fatpy/struct_mech/strain.py @@ -193,37 +193,6 @@ def calc_deviatoric_strain( # Von Mises functions -def calc_von_mises_strain_from_principals( - strain_vector_voigt: NDArray[np.float64], -) -> NDArray[np.float64]: - r"""Calculate von Mises equivalent strain for each strain state using principals. - - ??? abstract "Math Equations" - $$ \varepsilon_{vM} = \frac{\sqrt{2}}{3} \sqrt{ - (\varepsilon_1-\varepsilon_2)^2 + (\varepsilon_2-\varepsilon_3)^2 + - (\varepsilon_3-\varepsilon_1)^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) - - principals = calc_principal_strains(strain_vector_voigt) - e1 = principals[..., 0] - e2 = principals[..., 1] - e3 = principals[..., 2] - - return np.sqrt((2 / 9.0) * ((e1 - e2) ** 2 + (e2 - e3) ** 2 + (e3 - e1) ** 2)) - - def calc_von_mises_strain( strain_vector_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: diff --git a/tests/struct_mech/test_strain.py b/tests/struct_mech/test_strain.py index 8340aa2..b0417a0 100644 --- a/tests/struct_mech/test_strain.py +++ b/tests/struct_mech/test_strain.py @@ -194,27 +194,6 @@ def test_calc_deviatoric_strain( assert np.allclose(deviator[idx], deviator_voigt, atol=1e-12) -def test_calc_von_mises_from_principals( - strain_vector_sample: NDArray[np.float64], -) -> None: - """Test calculation of von Mises strain from principals. - - Uses the definition based on the strain deviator tensor. - - Args: - strain_vector_sample (NDArray[np.float64]): Sample strain vectors. - """ - von_mises = strain.calc_von_mises_strain_from_principals(strain_vector_sample) - - assert von_mises.shape == strain_vector_sample.shape[:-1] - - for idx in np.ndindex(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)) - assert np.isclose(von_mises[idx], vm, atol=1e-12) - - def test_calc_von_mises( strain_vector_sample: NDArray[np.float64], ) -> None: From 837471c1244acae42caef1be6f055981b45e315d Mon Sep 17 00:00:00 2001 From: Tomas Karas Date: Wed, 29 Oct 2025 15:25:52 +0100 Subject: [PATCH 8/9] docs fix --- src/fatpy/struct_mech/strain.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/fatpy/struct_mech/strain.py b/src/fatpy/struct_mech/strain.py index a12db1c..0cb2d3c 100644 --- a/src/fatpy/struct_mech/strain.py +++ b/src/fatpy/struct_mech/strain.py @@ -268,11 +268,11 @@ def calc_signed_von_mises_by_max_abs_principal( 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 maximum absolute principal strain - to zero. + 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 maximum absolute principal strain - to zero. + atol: Absolute tolerance for comparing the average of maximum and minimum + principal strain to zero. Default is 1e-8. Returns: From 355b222b1a9e261c7be730c4dcb2a93c9f3a9e49 Mon Sep 17 00:00:00 2001 From: Tomas Karas Date: Wed, 29 Oct 2025 16:04:33 +0100 Subject: [PATCH 9/9] docs format update --- src/fatpy/struct_mech/strain.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/fatpy/struct_mech/strain.py b/src/fatpy/struct_mech/strain.py index 0cb2d3c..cf4dac2 100644 --- a/src/fatpy/struct_mech/strain.py +++ b/src/fatpy/struct_mech/strain.py @@ -269,11 +269,9 @@ def calc_signed_von_mises_by_max_abs_principal( 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. + 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. + principal strain to zero. Default is 1e-8. Returns: Array of shape (...). Signed von Mises equivalent strain for each entry.