From b382cfebb29b31fad103a87e68a713cc615cfa05 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Sat, 16 Aug 2025 16:56:04 +0200 Subject: [PATCH 01/17] Implemented base eq stress calc functions --- .../__init__.py | 0 src/fatpy/struct_mech/stress.py | 250 ++++++++++++++++++ src/fatpy/struct_mech/test.py | 0 .../transformations.py | 0 src/fatpy/structural_mechanics/eq_stresses.py | 1 - 5 files changed, 250 insertions(+), 1 deletion(-) rename src/fatpy/{structural_mechanics => struct_mech}/__init__.py (100%) create mode 100644 src/fatpy/struct_mech/stress.py create mode 100644 src/fatpy/struct_mech/test.py rename src/fatpy/{structural_mechanics => struct_mech}/transformations.py (100%) delete mode 100644 src/fatpy/structural_mechanics/eq_stresses.py diff --git a/src/fatpy/structural_mechanics/__init__.py b/src/fatpy/struct_mech/__init__.py similarity index 100% rename from src/fatpy/structural_mechanics/__init__.py rename to src/fatpy/struct_mech/__init__.py diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py new file mode 100644 index 0000000..472f8c8 --- /dev/null +++ b/src/fatpy/struct_mech/stress.py @@ -0,0 +1,250 @@ +"""Calculate fundamental stress metrics and invariants. + +These tools provide principal stresses, maximum shear stress, hydrostatic stress, +von Mises equivalent stress, and invariants of the stress tensor. They are +essential for strength, fatigue, and fracture analyses under both uniaxial and +multiaxial loading conditions. + +""" + +import numpy as np +from numpy.typing import NDArray + + +def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + """Calculate the hydrostatic (mean normal) stress for each stress state. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Returns: + Array of shape (n,). Hydrostatic stress for each input state. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: + raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + + return (stress_voigt[:, 0] + stress_voigt[:, 1] + stress_voigt[:, 2]) / 3.0 + + +def _voigt_to_tensor(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + """Convert Voigt stress vectors to symmetric 3x3 tensors. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Returns: + Array of shape (n, 3, 3). Symmetric stress tensors. + """ + n = stress_voigt.shape[0] + tensor = np.zeros((n, 3, 3), dtype=np.float64) + tensor[:, 0, 0] = stress_voigt[:, 0] + tensor[:, 1, 1] = stress_voigt[:, 1] + tensor[:, 2, 2] = stress_voigt[:, 2] + tensor[:, 1, 2] = tensor[:, 2, 1] = stress_voigt[:, 3] + tensor[:, 0, 2] = tensor[:, 2, 0] = stress_voigt[:, 4] + tensor[:, 0, 1] = tensor[:, 1, 0] = stress_voigt[:, 5] + return tensor + + +def calc_principal_stresses_and_directions( + stress_voigt: NDArray[np.float64], +) -> tuple[NDArray[np.float64], NDArray[np.float64]]: + """Compute principal stresses and principal directions for each stress state. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Returns: + Tuple (eigvals, eigvecs): + - eigvals: Array of shape (n, 3). Principal stresses (ascending order). + - eigvecs: Array of shape (n, 3, 3). Principal directions (columns are + eigenvectors). + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: + raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + + tensor = _voigt_to_tensor(stress_voigt) + return np.linalg.eigh(tensor) + + +def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + """Compute principal stresses for each stress state. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Returns: + Array of shape (n, 3). Principal stresses (sorted descending for each row). + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: + raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + + tensor = _voigt_to_tensor(stress_voigt) + eigvals = np.linalg.eigvalsh(tensor) + + return np.sort(eigvals, axis=1)[:, ::-1] + + +def calc_principal_directions(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + """Compute principal directions (eigenvectors) for each stress state. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Returns: + Array of shape (n, 3, 3). Principal directions (columns are eigenvectors). + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: + raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + + tensor = _voigt_to_tensor(stress_voigt) + _, eigvecs = np.linalg.eigh(tensor) + + return eigvecs + + +def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + """Compute the first, second, and third invariants for each stress state. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Returns: + Array of shape (n, 3). Columns are (I1, I2, I3) for each row. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: + raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + + tensor = _voigt_to_tensor(stress_voigt) + i1 = np.trace(tensor, axis1=1, axis2=2) + i2 = 0.5 * (i1**2 - np.trace(np.matmul(tensor, tensor), axis1=1, axis2=2)) + i3 = np.linalg.det(tensor) + + return np.stack((i1, i2, i3), axis=1) + + +def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + """Compute von Mises equivalent stress for each stress state. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Returns: + Array of shape (n,). Von Mises equivalent stress for each row. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: + raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + + sx, sy, sz, syz, sxz, sxy = stress_voigt.T + + return np.sqrt( # type: ignore[no-any-return] + 0.5 + * ( + (sx - sy) ** 2 + + (sy - sz) ** 2 + + (sz - sx) ** 2 + + 6 * (sxy**2 + syz**2 + sxz**2) + ) + ) + + +def calc_signed_von_mises(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + """Compute signed von Mises stress for each stress state. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Returns: + Array of shape (n,). Signed von Mises stress for each row (sign from hydrostatic + part). + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + vm = calc_von_mises_stress(stress_voigt) + sign = np.sign(calc_hydrostatic_stress(stress_voigt)) + return sign * vm + + +def calc_tresca_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + """Compute Tresca (maximum shear) stress for each stress state. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Returns: + Array of shape (n,). Tresca stress for each row. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + principal = calc_principal_stresses(stress_voigt) + return 0.5 * (principal[:, 0] - principal[:, 2]) + + +def calc_signed_tresca(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + """Compute signed Tresca stress for each stress state. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Returns: + Array of shape (n,). Signed Tresca stress for each row (sign from hydrostatic + part). + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + t = calc_tresca_stress(stress_voigt) + sign = np.sign(calc_hydrostatic_stress(stress_voigt)) + return sign * t + + +if __name__ == "__main__": + # Example usage + stress_example = np.array( + [ + [100, 50, 25, 0, 0, 0], + [-50, 50, 0, 0, 0, 0], + [0, 0, 0, 50, 0, 0], + [100, 0, 0, 0, 0, 0], + ] + ) + print("Stress:", stress_example) + # print("Principal Stresses:", calc_principal_stresses(stress_example)) + # print("Hydrostatic Stress:", calc_hydrostatic_stress(stress_example)) + # print("Von Mises Stress:", calc_von_mises_stress(stress_example)) + # print("Tresca Stress:", calc_tresca_stress(stress_example)) + + eigvals, eigvecs = calc_principal_stresses_and_directions(stress_example) + + print("Principal Stresses:", eigvals) + print("Principal Directions:", eigvecs) diff --git a/src/fatpy/struct_mech/test.py b/src/fatpy/struct_mech/test.py new file mode 100644 index 0000000..e69de29 diff --git a/src/fatpy/structural_mechanics/transformations.py b/src/fatpy/struct_mech/transformations.py similarity index 100% rename from src/fatpy/structural_mechanics/transformations.py rename to src/fatpy/struct_mech/transformations.py diff --git a/src/fatpy/structural_mechanics/eq_stresses.py b/src/fatpy/structural_mechanics/eq_stresses.py deleted file mode 100644 index 1adeb21..0000000 --- a/src/fatpy/structural_mechanics/eq_stresses.py +++ /dev/null @@ -1 +0,0 @@ -"""Contains the function to calculate equivalent stresses.""" From 52303ad18e44ad8787c4d8cf38d52da431c25d17 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Sat, 16 Aug 2025 18:17:51 +0200 Subject: [PATCH 02/17] check shape function added, types modified --- .pre-commit-config.yaml | 6 ++-- src/fatpy/struct_mech/stress.py | 56 ++++++++++++++++++++------------- src/fatpy/struct_mech/test.py | 0 3 files changed, 38 insertions(+), 24 deletions(-) delete mode 100644 src/fatpy/struct_mech/test.py diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 128c5b1..6731f47 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,6 +1,6 @@ repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v5.0.0 + rev: v6.0.0 hooks: - id: trailing-whitespace - id: end-of-file-fixer @@ -8,7 +8,7 @@ repos: - id: check-added-large-files - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.11.6 + rev: v0.12.9 hooks: - id: ruff args: @@ -17,6 +17,6 @@ repos: - id: ruff-format - repo: https://github.com/pre-commit/mirrors-mypy - rev: v1.15.0 + rev: v1.17.1 hooks: - id: mypy diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index 472f8c8..1d748e5 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -11,6 +11,20 @@ from numpy.typing import NDArray +def _check_shape(stress_voigt: NDArray[np.float64]) -> None: + """Check the shape of the input stress array. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt + notation. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: + raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + + def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: """Calculate the hydrostatic (mean normal) stress for each stress state. @@ -24,8 +38,7 @@ def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.flo Raises: ValueError: If input is not a 2D array with 6 columns. """ - if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: - raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + _check_shape(stress_voigt) return (stress_voigt[:, 0] + stress_voigt[:, 1] + stress_voigt[:, 2]) / 3.0 @@ -69,11 +82,11 @@ def calc_principal_stresses_and_directions( Raises: ValueError: If input is not a 2D array with 6 columns. """ - if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: - raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + _check_shape(stress_voigt) tensor = _voigt_to_tensor(stress_voigt) - return np.linalg.eigh(tensor) + eigvals, eigvecs = np.linalg.eigh(tensor) + return eigvals, eigvecs def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: @@ -89,8 +102,7 @@ def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.flo Raises: ValueError: If input is not a 2D array with 6 columns. """ - if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: - raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + _check_shape(stress_voigt) tensor = _voigt_to_tensor(stress_voigt) eigvals = np.linalg.eigvalsh(tensor) @@ -111,8 +123,7 @@ def calc_principal_directions(stress_voigt: NDArray[np.float64]) -> NDArray[np.f Raises: ValueError: If input is not a 2D array with 6 columns. """ - if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: - raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + _check_shape(stress_voigt) tensor = _voigt_to_tensor(stress_voigt) _, eigvecs = np.linalg.eigh(tensor) @@ -133,8 +144,7 @@ def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.floa Raises: ValueError: If input is not a 2D array with 6 columns. """ - if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: - raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + _check_shape(stress_voigt) tensor = _voigt_to_tensor(stress_voigt) i1 = np.trace(tensor, axis1=1, axis2=2) @@ -157,12 +167,16 @@ def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float Raises: ValueError: If input is not a 2D array with 6 columns. """ - if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: - raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + _check_shape(stress_voigt) - sx, sy, sz, syz, sxz, sxy = stress_voigt.T + sx = stress_voigt[:, 0] + sy = stress_voigt[:, 1] + sz = stress_voigt[:, 2] + syz = stress_voigt[:, 3] + sxz = stress_voigt[:, 4] + sxy = stress_voigt[:, 5] - return np.sqrt( # type: ignore[no-any-return] + return np.sqrt( 0.5 * ( (sx - sy) ** 2 @@ -239,12 +253,12 @@ def calc_signed_tresca(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64] ] ) print("Stress:", stress_example) - # print("Principal Stresses:", calc_principal_stresses(stress_example)) - # print("Hydrostatic Stress:", calc_hydrostatic_stress(stress_example)) - # print("Von Mises Stress:", calc_von_mises_stress(stress_example)) - # print("Tresca Stress:", calc_tresca_stress(stress_example)) - + print("Principal Stresses:", calc_principal_stresses(stress_example)) + print("Hydrostatic Stress:", calc_hydrostatic_stress(stress_example)) + print("Von Mises Stress:", calc_von_mises_stress(stress_example)) + print("signed Von Mises Stress:", calc_signed_von_mises(stress_example)) + print("Tresca Stress:", calc_tresca_stress(stress_example)) + print("signed Tresca Stress:", calc_signed_tresca(stress_example)) eigvals, eigvecs = calc_principal_stresses_and_directions(stress_example) - print("Principal Stresses:", eigvals) print("Principal Directions:", eigvecs) diff --git a/src/fatpy/struct_mech/test.py b/src/fatpy/struct_mech/test.py deleted file mode 100644 index e69de29..0000000 From 4457ec3d6c15d12527b112fa7a6fbbb8a8f9885f Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Sun, 17 Aug 2025 17:37:19 +0200 Subject: [PATCH 03/17] added principal signed variants, voigt functions moved to utils --- .python-version | 1 + src/fatpy/struct_mech/stress.py | 268 ++++++++++++++++++-------------- src/fatpy/utils/voigt.py | 38 +++++ 3 files changed, 194 insertions(+), 113 deletions(-) create mode 100644 .python-version create mode 100644 src/fatpy/utils/voigt.py diff --git a/.python-version b/.python-version new file mode 100644 index 0000000..e4fba21 --- /dev/null +++ b/.python-version @@ -0,0 +1 @@ +3.12 diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index 1d748e5..2b634ba 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -1,36 +1,30 @@ """Calculate fundamental stress metrics and invariants. -These tools provide principal stresses, maximum shear stress, hydrostatic stress, -von Mises equivalent stress, and invariants of the stress tensor. They are -essential for strength, fatigue, and fracture analyses under both uniaxial and -multiaxial loading conditions. +These functions provide principal stresses, maximum shear (Tresca) stress, +hydrostatic stress, von Mises equivalent stress, and invariants of the stress +tensor. They are essential for strength, fatigue, and fracture analyses under +both uniaxial and multiaxial loading conditions. + +Conventions: +- Principal stresses 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 - -def _check_shape(stress_voigt: NDArray[np.float64]) -> None: - """Check the shape of the input stress array. - - Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. - - Raises: - ValueError: If input is not a 2D array with 6 columns. - """ - if stress_voigt.ndim != 2 or stress_voigt.shape[1] != 6: - raise ValueError("Input must be a n x 6 matrix in Voigt notation.") +from fatpy.utils import voigt def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: """Calculate the hydrostatic (mean normal) stress for each stress state. Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. Returns: Array of shape (n,). Hydrostatic stress for each input state. @@ -38,105 +32,90 @@ def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.flo Raises: ValueError: If input is not a 2D array with 6 columns. """ - _check_shape(stress_voigt) + voigt.check_shape(stress_voigt) return (stress_voigt[:, 0] + stress_voigt[:, 1] + stress_voigt[:, 2]) / 3.0 -def _voigt_to_tensor(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Convert Voigt stress vectors to symmetric 3x3 tensors. - - Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. - - Returns: - Array of shape (n, 3, 3). Symmetric stress tensors. - """ - n = stress_voigt.shape[0] - tensor = np.zeros((n, 3, 3), dtype=np.float64) - tensor[:, 0, 0] = stress_voigt[:, 0] - tensor[:, 1, 1] = stress_voigt[:, 1] - tensor[:, 2, 2] = stress_voigt[:, 2] - tensor[:, 1, 2] = tensor[:, 2, 1] = stress_voigt[:, 3] - tensor[:, 0, 2] = tensor[:, 2, 0] = stress_voigt[:, 4] - tensor[:, 0, 1] = tensor[:, 1, 0] = stress_voigt[:, 5] - return tensor - - def calc_principal_stresses_and_directions( stress_voigt: NDArray[np.float64], ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: - """Compute principal stresses and principal directions for each stress state. + """Calculate principal stresses and principal directions for each state. Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. Returns: Tuple (eigvals, eigvecs): - - eigvals: Array of shape (n, 3). Principal stresses (ascending order). + - eigvals: Array of shape (n, 3). Principal stresses (descending: σ1 ≥ σ2 ≥ σ3) - eigvecs: Array of shape (n, 3, 3). Principal directions (columns are - eigenvectors). + eigenvectors) aligned with eigvals in the same order. Raises: ValueError: If input is not a 2D array with 6 columns. """ - _check_shape(stress_voigt) + voigt.check_shape(stress_voigt) + + tensor = voigt.voigt_to_tensor(stress_voigt) + eigvals, eigvecs = np.linalg.eigh(tensor) # ascending by default + + # Reorder to descending and keep eigenvectors aligned with eigenvalues + idx = np.argsort(eigvals, axis=1)[:, ::-1] # (n,3) + eigvals_sorted = np.take_along_axis(eigvals, idx, axis=1) + eigvecs_sorted = np.take_along_axis(eigvecs, idx[:, None, :], axis=2) - tensor = _voigt_to_tensor(stress_voigt) - eigvals, eigvecs = np.linalg.eigh(tensor) - return eigvals, eigvecs + return eigvals_sorted, eigvecs_sorted def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Compute principal stresses for each stress state. + """Calculate principal stresses for each stress state. Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. Returns: - Array of shape (n, 3). Principal stresses (sorted descending for each row). + Array of shape (n, 3). Principal stresses (descending: σ1 ≥ σ2 ≥ σ3). Raises: ValueError: If input is not a 2D array with 6 columns. """ - _check_shape(stress_voigt) + voigt.check_shape(stress_voigt) - tensor = _voigt_to_tensor(stress_voigt) + tensor = voigt.voigt_to_tensor(stress_voigt) eigvals = np.linalg.eigvalsh(tensor) return np.sort(eigvals, axis=1)[:, ::-1] def calc_principal_directions(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Compute principal directions (eigenvectors) for each stress state. + """Calculate principal directions (eigenvectors) for each stress state. Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. Returns: - Array of shape (n, 3, 3). Principal directions (columns are eigenvectors). + Array of shape (n, 3, 3). Principal directions (columns are eigenvectors) + aligned with descending principal stresses: σ1, σ2, σ3. Raises: ValueError: If input is not a 2D array with 6 columns. """ - _check_shape(stress_voigt) - - tensor = _voigt_to_tensor(stress_voigt) - _, eigvecs = np.linalg.eigh(tensor) + voigt.check_shape(stress_voigt) + # Reuse the ordering logic from the combined function to ensure consistency + _, eigvecs = calc_principal_stresses_and_directions(stress_voigt) return eigvecs def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Compute the first, second, and third invariants for each stress state. + """Calculate the first, second, and third invariants for each stress state. Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. Returns: Array of shape (n, 3). Columns are (I1, I2, I3) for each row. @@ -144,22 +123,24 @@ def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.floa Raises: ValueError: If input is not a 2D array with 6 columns. """ - _check_shape(stress_voigt) + voigt.check_shape(stress_voigt) - tensor = _voigt_to_tensor(stress_voigt) - i1 = np.trace(tensor, axis1=1, axis2=2) - i2 = 0.5 * (i1**2 - np.trace(np.matmul(tensor, tensor), axis1=1, axis2=2)) - i3 = np.linalg.det(tensor) + tensor = voigt.voigt_to_tensor(stress_voigt) + invariant_1 = np.trace(tensor, axis1=1, axis2=2) + invariant_2 = 0.5 * ( + invariant_1**2 - np.trace(np.matmul(tensor, tensor), axis1=1, axis2=2) + ) + invariant_3 = np.linalg.det(tensor) - return np.stack((i1, i2, i3), axis=1) + return np.stack((invariant_1, invariant_2, invariant_3), axis=1) def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Compute von Mises equivalent stress for each stress state. + """Calculate von Mises equivalent stress for each stress state. Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. Returns: Array of shape (n,). Von Mises equivalent stress for each row. @@ -167,7 +148,7 @@ def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float Raises: ValueError: If input is not a 2D array with 6 columns. """ - _check_shape(stress_voigt) + voigt.check_shape(stress_voigt) sx = stress_voigt[:, 0] sy = stress_voigt[:, 1] @@ -187,78 +168,139 @@ def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float ) -def calc_signed_von_mises(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Compute signed von Mises stress for each stress state. +def calc_signed_von_mises_by_hydrostatic( + stress_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + """Calculate signed von Mises stress for each stress state. + + Sign is determined by the hydrostatic stress. Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. Returns: - Array of shape (n,). Signed von Mises stress for each row (sign from hydrostatic - part). + Array of shape (n,). Signed von Mises stress for each row. Raises: ValueError: If input is not a 2D array with 6 columns. """ - vm = calc_von_mises_stress(stress_voigt) + von_mises = calc_von_mises_stress(stress_voigt) sign = np.sign(calc_hydrostatic_stress(stress_voigt)) - return sign * vm + return sign * von_mises + + +def calc_signed_von_mises_by_max_abs_principal( + stress_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + """Calculate signed von Mises stress for each stress state. + + Sign is determined by the maximum absolute principal stress value. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. + + Returns: + Array of shape (n,). Signed von Mises stress for each row. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + von_mises = calc_von_mises_stress(stress_voigt) + principals = calc_principal_stresses(stress_voigt) + + # Find the principal stress with maximum absolute value and get its sign + 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) + + return sign * von_mises def calc_tresca_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Compute Tresca (maximum shear) stress for each stress state. + """Calculate Tresca (maximum shear) stress for each stress state. Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. Returns: - Array of shape (n,). Tresca stress for each row. + Array of shape (n,). Tresca stress for each row. For principal stresses + σ1 ≥ σ2 ≥ σ3, Tresca = (σ1 − σ3)/2. Raises: ValueError: If input is not a 2D array with 6 columns. """ - principal = calc_principal_stresses(stress_voigt) - return 0.5 * (principal[:, 0] - principal[:, 2]) + principals = calc_principal_stresses(stress_voigt) + return 0.5 * (principals[:, 0] - principals[:, 2]) + +def calc_signed_tresca_by_hydrostatic( + stress_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + """Calculate signed Tresca stress for each stress state. -def calc_signed_tresca(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Compute signed Tresca stress for each stress state. + Sign is determined by the hydrostatic stress. Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt - notation. + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. Returns: - Array of shape (n,). Signed Tresca stress for each row (sign from hydrostatic - part). + Array of shape (n,). Signed Tresca stress for each row. Raises: ValueError: If input is not a 2D array with 6 columns. """ - t = calc_tresca_stress(stress_voigt) + tresca = calc_tresca_stress(stress_voigt) sign = np.sign(calc_hydrostatic_stress(stress_voigt)) - return sign * t + return sign * tresca + + +def calc_signed_tresca_by_max_abs_principal( + stress_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + """Calculate signed Tresca stress for each stress state. + + Sign is determined by the maximum absolute principal stress value. + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. + + Returns: + Array of shape (n,). Signed Tresca stress for each row. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + tresca = calc_tresca_stress(stress_voigt) + principals = calc_principal_stresses(stress_voigt) + + # Find the principal stress with maximum absolute value and get its sign + 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) + + return sign * tresca if __name__ == "__main__": # Example usage stress_example = np.array( [ - [100, 50, 25, 0, 0, 0], - [-50, 50, 0, 0, 0, 0], - [0, 0, 0, 50, 0, 0], - [100, 0, 0, 0, 0, 0], + [10, 20, 30, 0, 0, 0], + [50, -50, 0, 0, 0, 10], + [0, 0, 1000, 0, 0, -50], ] ) - print("Stress:", stress_example) - print("Principal Stresses:", calc_principal_stresses(stress_example)) - print("Hydrostatic Stress:", calc_hydrostatic_stress(stress_example)) - print("Von Mises Stress:", calc_von_mises_stress(stress_example)) - print("signed Von Mises Stress:", calc_signed_von_mises(stress_example)) - print("Tresca Stress:", calc_tresca_stress(stress_example)) - print("signed Tresca Stress:", calc_signed_tresca(stress_example)) - eigvals, eigvecs = calc_principal_stresses_and_directions(stress_example) - print("Principal Stresses:", eigvals) - print("Principal Directions:", eigvecs) + # print("Principal Stresses:", calc_principal_stresses(stress_example)) + # print("Hydrostatic Stress:", calc_hydrostatic_stress(stress_example)) + # print("Von Mises Stress:", calc_von_mises_stress(stress_example)) + # print("signed Von Mises Stress:", calc_signed_von_mises(stress_example)) + # print("Tresca Stress:", calc_tresca_stress(stress_example)) + # print("signed Tresca Stress:", calc_signed_tresca(stress_example)) + # eigvals, eigvecs = calc_principal_stresses_and_directions(stress_example) + # print("Principal Stresses:", eigvals) + # print("Principal Directions:", eigvecs) diff --git a/src/fatpy/utils/voigt.py b/src/fatpy/utils/voigt.py new file mode 100644 index 0000000..23cd4cb --- /dev/null +++ b/src/fatpy/utils/voigt.py @@ -0,0 +1,38 @@ +"""Helper functions for handling vectors using Voigt notation.""" + +import numpy as np +from numpy.typing import NDArray + + +def check_shape(voigt_vec: NDArray[np.float64]) -> None: + """Check the shape of the input Voigt array. + + Args: + voigt_vec: Array of shape (n, 6). Each row is a vector in Voigt notation. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + if voigt_vec.ndim != 2 or voigt_vec.shape[1] != 6: + raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + + +def voigt_to_tensor(voigt_vec: NDArray[np.float64]) -> NDArray[np.float64]: + """Convert Voigt vectors to symmetric 3x3 tensors. + + Args: + voigt_vec: Array of shape (n, 6). Each row is a vector in Voigt notation. + (σ11, σ22, σ33, σ23, σ13, σ12) + + Returns: + Array of shape (n, 3, 3). Symmetric tensors. + """ + n = voigt_vec.shape[0] + tensor = np.zeros((n, 3, 3), dtype=np.float64) + tensor[:, 0, 0] = voigt_vec[:, 0] + tensor[:, 1, 1] = voigt_vec[:, 1] + tensor[:, 2, 2] = voigt_vec[:, 2] + tensor[:, 1, 2] = tensor[:, 2, 1] = voigt_vec[:, 3] + tensor[:, 0, 2] = tensor[:, 2, 0] = voigt_vec[:, 4] + tensor[:, 0, 1] = tensor[:, 1, 0] = voigt_vec[:, 5] + return tensor From fd1149f7fe2b0bec0c192deb8830907262f969ae Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Sun, 17 Aug 2025 17:39:39 +0200 Subject: [PATCH 04/17] added missing tabulator in docstring --- src/fatpy/struct_mech/stress.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index 2b634ba..d47f8f9 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -267,7 +267,7 @@ def calc_signed_tresca_by_max_abs_principal( Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + Voigt notation. Returns: Array of shape (n,). Signed Tresca stress for each row. From 7b68691bae7e36a81742d8ffdeed9680d0b5d946 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Tue, 16 Sep 2025 01:11:37 +0200 Subject: [PATCH 05/17] pr docstrings update, wip --- src/fatpy/struct_mech/stress.py | 35 +++++++++++++++++++++++++++++---- 1 file changed, 31 insertions(+), 4 deletions(-) diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index d47f8f9..d3f52a7 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -20,7 +20,13 @@ def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Calculate the hydrostatic (mean normal) stress for each stress state. + r"""Calculate the hydrostatic (mean normal) stress for each stress state. + + ??? abstract "Math Equation" + $$ + \sigma_H = \frac{1}{3} tr(\sigma) = + \frac{1}{3} (\sigma_{xx} + \sigma_{yy} + \sigma_{zz}) + $$ Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in @@ -40,7 +46,7 @@ def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.flo def calc_principal_stresses_and_directions( stress_voigt: NDArray[np.float64], ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: - """Calculate principal stresses and principal directions for each state. + r"""Calculate principal stresses and principal directions for each state. Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in @@ -69,7 +75,17 @@ def calc_principal_stresses_and_directions( def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Calculate principal stresses for each stress state. + r"""Calculate principal stresses for each stress state. + + .. math:: + + Ai(z) = \frac{1}{\\pi \\sqrt{3}} K_{1/3}(t) + + Ai'(z) = -\frac{z}{\\pi \\sqrt{3}} K_{2/3}(t) + + Bi(z) = \\sqrt{\frac{z}{3}} \\left(I_{-1/3}(t) + I_{1/3}(t) \right) + + Bi'(z) = \frac{z}{\\sqrt{3}} \\left(I_{-2/3}(t) + I_{2/3}(t)\right) Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in @@ -111,7 +127,16 @@ def calc_principal_directions(stress_voigt: NDArray[np.float64]) -> NDArray[np.f def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Calculate the first, second, and third invariants for each stress state. + r"""Calculate the first, second, and third invariants for each stress state. + + ??? abstract "Math Equations" + $$ + \begin{align*} + I_1 &= tr(\sigma), \\ + I_2 &= \frac{1}{2}(I_1^{2} - tr(\sigma^{2})), \\ + I_3 &= \det(\sigma) + \end{align*} + $$ Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in @@ -304,3 +329,5 @@ def calc_signed_tresca_by_max_abs_principal( # eigvals, eigvecs = calc_principal_stresses_and_directions(stress_example) # print("Principal Stresses:", eigvals) # print("Principal Directions:", eigvecs) + # print("Principal Stresses:", eigvals) + # print("Principal Directions:", eigvecs) From b22e6aecb70bb28920bc335da637d480472fb570 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Wed, 17 Sep 2025 01:18:59 +0200 Subject: [PATCH 06/17] Added math equations - wip --- src/fatpy/struct_mech/stress.py | 110 ++++++++++++++++++++++++-------- 1 file changed, 85 insertions(+), 25 deletions(-) diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index d3f52a7..f01ebfc 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -22,10 +22,10 @@ def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: r"""Calculate the hydrostatic (mean normal) stress for each stress state. - ??? abstract "Math Equation" + ??? abstract "Math Equations" $$ \sigma_H = \frac{1}{3} tr(\sigma) = - \frac{1}{3} (\sigma_{xx} + \sigma_{yy} + \sigma_{zz}) + \frac{1}{3} (\sigma_{11} + \sigma_{22} + \sigma_{33}) $$ Args: @@ -48,6 +48,12 @@ def calc_principal_stresses_and_directions( ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: r"""Calculate principal stresses and principal directions for each state. + ??? abstract "Math Equations" + Principal stresses and directions are found by solving the eigenvalue problem + for the stress tensor: + + $$ \sigma \mathbf{v} = \lambda \mathbf{v} $$ + Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt notation. @@ -77,15 +83,11 @@ def calc_principal_stresses_and_directions( def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: r"""Calculate principal stresses for each stress state. - .. math:: - - Ai(z) = \frac{1}{\\pi \\sqrt{3}} K_{1/3}(t) - - Ai'(z) = -\frac{z}{\\pi \\sqrt{3}} K_{2/3}(t) - - Bi(z) = \\sqrt{\frac{z}{3}} \\left(I_{-1/3}(t) + I_{1/3}(t) \right) + ??? abstract "Math Equations" + Principal stresses are found by solving the eigenvalue problem + for the stress tensor: - Bi'(z) = \frac{z}{\\sqrt{3}} \\left(I_{-2/3}(t) + I_{2/3}(t)\right) + $$ \sigma \mathbf{v} = \lambda \mathbf{v} $$ Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in @@ -106,7 +108,13 @@ def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.flo def calc_principal_directions(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Calculate principal directions (eigenvectors) for each stress state. + r"""Calculate principal directions (eigenvectors) for each stress state. + + ??? abstract "Math Equations" + Principal directions are found by solving the eigenvalue problem + for the stress tensor: + + $$ \sigma \mathbf{v} = \lambda \mathbf{v} $$ Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in @@ -161,7 +169,16 @@ def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.floa def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Calculate von Mises equivalent stress for each stress state. + r"""Calculate von Mises equivalent stress for each stress state. + + ??? abstract "Math Equations" + $$ + \sigma_{vM} = \tfrac{\sqrt{2}}{2}\sqrt{ + (\sigma_{11}-\sigma_{22})^2 + +(\sigma_{22}-\sigma_{33})^2 + +(\sigma_{33}-\sigma_{11})^2 + + 3(\sigma_{12}^2+\sigma_{23}^2+\sigma_{13}^2)} + $$ Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in @@ -196,10 +213,13 @@ def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float def calc_signed_von_mises_by_hydrostatic( stress_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: - """Calculate signed von Mises stress for each stress state. + r"""Calculate signed von Mises stress for each stress state. Sign is determined by the hydrostatic stress. + ??? abstract "Math Equations" + $$\sigma_{SvM} = sgn(\sigma_H) \cdot \sigma_{vM}$$ + Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt notation. @@ -218,10 +238,13 @@ def calc_signed_von_mises_by_hydrostatic( def calc_signed_von_mises_by_max_abs_principal( stress_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: - """Calculate signed von Mises stress for each stress state. + r"""Calculate signed von Mises stress for each stress state. Sign is determined by the maximum absolute principal stress value. + ??? abstract "Math Equations" + $$\sigma_{SvM} = sgn(\frac{\sigma_{1}+\sigma_{3}}{2}) \cdot \sigma_{vM}$$ + Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt notation. @@ -235,16 +258,45 @@ def calc_signed_von_mises_by_max_abs_principal( von_mises = calc_von_mises_stress(stress_voigt) principals = calc_principal_stresses(stress_voigt) - # Find the principal stress with maximum absolute value and get its sign - 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) + avg_13 = 0.5 * (principals[:, 0] + principals[:, 2]) + sign = np.sign(avg_13).astype(np.float64, copy=False) + + return sign * von_mises + + +def calc_signed_von_mises_by_first_invariant( + stress_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + r"""Calculate signed von Mises stress for each stress state. + + Sign is determined by the first invariant of the stress tensor. + + ??? abstract "Math Equations" + $$\sigma_{SvM} = sgn(tr(\sigma)) \cdot \sigma_{vM}$$ + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. + + Returns: + Array of shape (n,). Signed von Mises stress for each row. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + von_mises = calc_von_mises_stress(stress_voigt) + invariant_1 = stress_voigt[:, 0] + stress_voigt[:, 1] + stress_voigt[:, 2] + + sign = np.sign(invariant_1).astype(np.float64, copy=False) return sign * von_mises def calc_tresca_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - """Calculate Tresca (maximum shear) stress for each stress state. + r"""Calculate Tresca (maximum shear) stress for each stress state. + + ??? abstract "Math Equations" + $$\sigma_{\tau_{max}} = \frac{\sigma_{1} - \sigma_{3}}{2}$$ Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in @@ -264,10 +316,13 @@ def calc_tresca_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64] def calc_signed_tresca_by_hydrostatic( stress_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: - """Calculate signed Tresca stress for each stress state. + r"""Calculate signed Tresca stress for each stress state. Sign is determined by the hydrostatic stress. + ??? abstract "Math Equations" + $$\sigma_{S\tau_{max}} = sgn(\sigma_H) \cdot \sigma_{\tau_{max}}$$ + Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt notation. @@ -286,10 +341,16 @@ def calc_signed_tresca_by_hydrostatic( def calc_signed_tresca_by_max_abs_principal( stress_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: - """Calculate signed Tresca stress for each stress state. + r"""Calculate signed Tresca stress for each stress state. Sign is determined by the maximum absolute principal stress value. + ??? abstract "Math Equations" + $$ + \sigma_{S\tau_{max}} = sgn(\frac{\sigma_{1}+\sigma_{3}}{2}) + \cdot \sigma_{\tau_{max}} + $$ + Args: stress_voigt: Array of shape (n, 6). Each row is a stress vector in Voigt notation. @@ -303,10 +364,8 @@ def calc_signed_tresca_by_max_abs_principal( tresca = calc_tresca_stress(stress_voigt) principals = calc_principal_stresses(stress_voigt) - # Find the principal stress with maximum absolute value and get its sign - 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) + avg_13 = 0.5 * (principals[:, 0] + principals[:, 2]) + sign = np.sign(avg_13).astype(np.float64, copy=False) return sign * tresca @@ -331,3 +390,4 @@ def calc_signed_tresca_by_max_abs_principal( # print("Principal Directions:", eigvecs) # print("Principal Stresses:", eigvals) # print("Principal Directions:", eigvecs) + # print("Principal Directions:", eigvecs) From f0af8cbb8c5c4503eea04c0e2b3e9907148d9cf5 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Wed, 17 Sep 2025 08:51:07 +0200 Subject: [PATCH 07/17] Added stress deviator calculation --- src/fatpy/struct_mech/stress.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index f01ebfc..9f221f3 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -370,6 +370,33 @@ def calc_signed_tresca_by_max_abs_principal( return sign * tresca +def calc_stress_deviator( + stress_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + r"""Calculate stress deviator for each stress state. + + ??? abstract "Math Equations" + $$ \mathbf{s} = \sigma - \frac{1}{3} tr(\sigma) $$ + + Args: + stress_voigt: Array of shape (n, 6). Each row is a stress vector in + Voigt notation. + + Returns: + Array of shape (n,). Signed Tresca stress for each row. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + voigt.check_shape(stress_voigt) + hydrostatic = calc_hydrostatic_stress(stress_voigt) + + deviator = stress_voigt.copy() + deviator[:, :3] -= hydrostatic[:, None] + + return deviator + + if __name__ == "__main__": # Example usage stress_example = np.array( @@ -391,3 +418,5 @@ def calc_signed_tresca_by_max_abs_principal( # print("Principal Stresses:", eigvals) # print("Principal Directions:", eigvecs) # print("Principal Directions:", eigvecs) + # print("Principal Directions:", eigvecs) + # print("Principal Directions:", eigvecs) From f900cda2325d483d68f2a72e4fec26a2dddb4f19 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Tue, 23 Sep 2025 08:11:18 +0200 Subject: [PATCH 08/17] wip - rewriting to n dimensional --- src/fatpy/struct_mech/stress.py | 296 +++++++++++++++++-------------- src/fatpy/utils/voigt.py | 80 +++++++-- tests/struct_mech/test_stress.py | 248 ++++++++++++++++++++++++++ 3 files changed, 473 insertions(+), 151 deletions(-) create mode 100644 tests/struct_mech/test_stress.py diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index 9f221f3..7bfa259 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -6,41 +6,29 @@ both uniaxial and multiaxial loading conditions. Conventions: -- Principal stresses 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 +- Vectors use Voigt notation with shape (n, 6, ...). -from fatpy.utils import voigt + For stress tensors, the six Voigt components are: + (σ_11, σ_22, σ_33, σ_23, σ_13, σ_12) + (σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy) -def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: - r"""Calculate the hydrostatic (mean normal) stress for each stress state. +- Principal stresses are ordered in descending order throughout the module +(σ_1 ≥ σ_2 ≥ σ_3). - ??? abstract "Math Equations" - $$ - \sigma_H = \frac{1}{3} tr(\sigma) = - \frac{1}{3} (\sigma_{11} + \sigma_{22} + \sigma_{33}) - $$ +- Principal directions (eigenvectors) are aligned to this ordering +(columns correspond to σ_1, σ_2, σ_3). - Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. +""" - Returns: - Array of shape (n,). Hydrostatic stress for each input state. +# TODO: N-Dimensional support for trailing dims +# - voigt, tensor dimensions have to be the last axes +# - (..., 6), (..., 3, 3) - Raises: - ValueError: If input is not a 2D array with 6 columns. - """ - voigt.check_shape(stress_voigt) +import numpy as np +from numpy.typing import NDArray - return (stress_voigt[:, 0] + stress_voigt[:, 1] + stress_voigt[:, 2]) / 3.0 +from fatpy.utils import voigt def calc_principal_stresses_and_directions( @@ -55,14 +43,17 @@ def calc_principal_stresses_and_directions( $$ \sigma \mathbf{v} = \lambda \mathbf{v} $$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second axis of size 6 + contains the Voigt stress components. Any trailing dimensions are + preserved and the returned arrays will have those trailing dims. Returns: Tuple (eigvals, eigvecs): - - eigvals: Array of shape (n, 3). Principal stresses (descending: σ1 ≥ σ2 ≥ σ3) - - eigvecs: Array of shape (n, 3, 3). Principal directions (columns are - eigenvectors) aligned with eigvals in the same order. + - eigvals: Array of shape (n, 3, ...). Principal stresses + (descending: σ_1 ≥ σ_2 ≥ σ_3) with trailing dims preserved. + - eigvecs: Array of shape (n, 3, 3, ...). Principal directions (columns are + eigenvectors) aligned with eigvals in the same order. The last two + axes of this array are the 3x3 eigenvector matrix for each input. Raises: ValueError: If input is not a 2D array with 6 columns. @@ -70,12 +61,24 @@ def calc_principal_stresses_and_directions( voigt.check_shape(stress_voigt) tensor = voigt.voigt_to_tensor(stress_voigt) - eigvals, eigvecs = np.linalg.eigh(tensor) # ascending by default - # Reorder to descending and keep eigenvectors aligned with eigenvalues - idx = np.argsort(eigvals, axis=1)[:, ::-1] # (n,3) - eigvals_sorted = np.take_along_axis(eigvals, idx, axis=1) - eigvecs_sorted = np.take_along_axis(eigvecs, idx[:, None, :], axis=2) + # tensor has shape (n, ..., 3, 3). Use eigh over the last two axes. + eigvals, eigvecs = np.linalg.eigh(tensor) # eigvals shape: (n, ..., 3) + + # Sort eigenvalues descending along the last axis and reorder eigenvectors + # accordingly. eigvecs has shape (n, ..., 3, 3) where the last axis is + # the eigenvector corresponding to the eigenvalue at the same position. + sorted_indices = np.argsort(eigvals, axis=-1)[..., ::-1] + eigvals_sorted = np.take_along_axis(eigvals, sorted_indices, axis=-1) + + # For eigvecs we need to take along the last axis (eigenvector index) + # while preserving the matrix axis. Move axes so we can index easily: + # eigvecs currently: (..., 3, 3) where the last axis indexes components + # and the second-last axis indexes eigenvector number. We want to reorder + # the eigenvector-number axis. + eigvecs_sorted = np.take_along_axis( + eigvecs, np.expand_dims(sorted_indices, axis=-2), axis=-1 + ) return eigvals_sorted, eigvecs_sorted @@ -90,11 +93,12 @@ def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.flo $$ \sigma \mathbf{v} = \lambda \mathbf{v} $$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. Returns: - Array of shape (n, 3). Principal stresses (descending: σ1 ≥ σ2 ≥ σ3). + Array of shape (n, 3, ...). Principal stresses (descending: σ1 ≥ σ2 ≥ σ3). Raises: ValueError: If input is not a 2D array with 6 columns. @@ -104,7 +108,8 @@ def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.flo tensor = voigt.voigt_to_tensor(stress_voigt) eigvals = np.linalg.eigvalsh(tensor) - return np.sort(eigvals, axis=1)[:, ::-1] + # eigvals shape: (n, ..., 3). Return sorted descending along last axis. + return np.sort(eigvals, axis=-1)[..., ::-1] def calc_principal_directions(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: @@ -117,11 +122,12 @@ def calc_principal_directions(stress_voigt: NDArray[np.float64]) -> NDArray[np.f $$ \sigma \mathbf{v} = \lambda \mathbf{v} $$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. Returns: - Array of shape (n, 3, 3). Principal directions (columns are eigenvectors) + Array of shape (n, 3, 3, ...). Principal directions (columns are eigenvectors) aligned with descending principal stresses: σ1, σ2, σ3. Raises: @@ -147,11 +153,12 @@ def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.floa $$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. Returns: - Array of shape (n, 3). Columns are (I1, I2, I3) for each row. + Array of shape (n, 3, ...). Columns are (I1, I2, I3) for each entry. Raises: ValueError: If input is not a 2D array with 6 columns. @@ -159,15 +166,73 @@ def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.floa voigt.check_shape(stress_voigt) tensor = voigt.voigt_to_tensor(stress_voigt) - invariant_1 = np.trace(tensor, axis1=1, axis2=2) + invariant_1 = np.trace(tensor, axis1=-2, axis2=-1) invariant_2 = 0.5 * ( - invariant_1**2 - np.trace(np.matmul(tensor, tensor), axis1=1, axis2=2) + 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) + return np.stack((invariant_1, invariant_2, invariant_3), axis=-1) + + +def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: + r"""Calculate the hydrostatic (mean normal) stress for each stress state. + + ??? abstract "Math Equations" + $$ + \sigma_H = \frac{1}{3} tr(\sigma) = + \frac{1}{3} (\sigma_{11} + \sigma_{22} + \sigma_{33}) + $$ + + Args: + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. + + Returns: + Array of shape (n, ...). Hydrostatic stress for each input state. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + voigt.check_shape(stress_voigt) + + # Voigt normal components are at axis=1 indices 0,1,2. Preserve trailing dims. + return ( + stress_voigt[:, 0, ...] + stress_voigt[:, 1, ...] + stress_voigt[:, 2, ...] + ) / 3.0 + +def calc_stress_deviator( + stress_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + r"""Calculate stress deviator for each stress state. + + ??? abstract "Math Equations" + $$ \mathbf{s} = \sigma - \frac{1}{3} tr(\sigma) $$ + + Args: + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. + + Returns: + Array of shape (n, 6, ...). Stress deviator for each entry. + + Raises: + ValueError: If input is not a 2D array with 6 columns. + """ + voigt.check_shape(stress_voigt) + hydrostatic = calc_hydrostatic_stress(stress_voigt) + + deviator = stress_voigt.copy() + # Subtract hydrostatic from the first three Voigt components (axis=1) + deviator[:, :3, ...] = deviator[:, :3, ...] - hydrostatic[..., None] + + return deviator + +# Von Mises functions def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: r"""Calculate von Mises equivalent stress for each stress state. @@ -181,23 +246,24 @@ def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float $$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. Returns: - Array of shape (n,). Von Mises equivalent stress for each row. + Array of shape (n, ...). Von Mises equivalent stress for each entry. Raises: ValueError: If input is not a 2D array with 6 columns. """ voigt.check_shape(stress_voigt) - sx = stress_voigt[:, 0] - sy = stress_voigt[:, 1] - sz = stress_voigt[:, 2] - syz = stress_voigt[:, 3] - sxz = stress_voigt[:, 4] - sxy = stress_voigt[:, 5] + sx = stress_voigt[:, 0, ...] + sy = stress_voigt[:, 1, ...] + sz = stress_voigt[:, 2, ...] + syz = stress_voigt[:, 3, ...] + sxz = stress_voigt[:, 4, ...] + sxy = stress_voigt[:, 5, ...] return np.sqrt( 0.5 @@ -221,11 +287,12 @@ def calc_signed_von_mises_by_hydrostatic( $$\sigma_{SvM} = sgn(\sigma_H) \cdot \sigma_{vM}$$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. Returns: - Array of shape (n,). Signed von Mises stress for each row. + Array of shape (n, ...). Signed von Mises stress for each entry. Raises: ValueError: If input is not a 2D array with 6 columns. @@ -246,11 +313,12 @@ def calc_signed_von_mises_by_max_abs_principal( $$\sigma_{SvM} = sgn(\frac{\sigma_{1}+\sigma_{3}}{2}) \cdot \sigma_{vM}$$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. Returns: - Array of shape (n,). Signed von Mises stress for each row. + Array of shape (n, ...). Signed von Mises stress for each entry. Raises: ValueError: If input is not a 2D array with 6 columns. @@ -258,7 +326,7 @@ def calc_signed_von_mises_by_max_abs_principal( von_mises = calc_von_mises_stress(stress_voigt) principals = calc_principal_stresses(stress_voigt) - avg_13 = 0.5 * (principals[:, 0] + principals[:, 2]) + avg_13 = 0.5 * (principals[..., 0] + principals[..., 2]) sign = np.sign(avg_13).astype(np.float64, copy=False) return sign * von_mises @@ -275,23 +343,27 @@ def calc_signed_von_mises_by_first_invariant( $$\sigma_{SvM} = sgn(tr(\sigma)) \cdot \sigma_{vM}$$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. Returns: - Array of shape (n,). Signed von Mises stress for each row. + Array of shape (n, ...). Signed von Mises stress for each entry. Raises: ValueError: If input is not a 2D array with 6 columns. """ von_mises = calc_von_mises_stress(stress_voigt) - invariant_1 = stress_voigt[:, 0] + stress_voigt[:, 1] + stress_voigt[:, 2] + invariant_1 = ( + stress_voigt[:, 0, ...] + stress_voigt[:, 1, ...] + stress_voigt[:, 2, ...] + ) sign = np.sign(invariant_1).astype(np.float64, copy=False) return sign * von_mises +# Tresca functions def calc_tresca_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: r"""Calculate Tresca (maximum shear) stress for each stress state. @@ -299,18 +371,19 @@ def calc_tresca_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64] $$\sigma_{\tau_{max}} = \frac{\sigma_{1} - \sigma_{3}}{2}$$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. Returns: - Array of shape (n,). Tresca stress for each row. For principal stresses - σ1 ≥ σ2 ≥ σ3, Tresca = (σ1 − σ3)/2. + Array of shape (n, ...). Tresca stress for each entry. For principal + stresses σ1 ≥ σ2 ≥ σ3, Tresca = (σ1 − σ3)/2. Raises: ValueError: If input is not a 2D array with 6 columns. """ principals = calc_principal_stresses(stress_voigt) - return 0.5 * (principals[:, 0] - principals[:, 2]) + return 0.5 * (principals[..., 0] - principals[..., 2]) def calc_signed_tresca_by_hydrostatic( @@ -324,11 +397,12 @@ def calc_signed_tresca_by_hydrostatic( $$\sigma_{S\tau_{max}} = sgn(\sigma_H) \cdot \sigma_{\tau_{max}}$$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. Returns: - Array of shape (n,). Signed Tresca stress for each row. + Array of shape (n, ...). Signed Tresca stress for each entry. Raises: ValueError: If input is not a 2D array with 6 columns. @@ -352,11 +426,12 @@ def calc_signed_tresca_by_max_abs_principal( $$ Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. + stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 + contains the Voigt stress components. Additional trailing dimensions + are supported. Returns: - Array of shape (n,). Signed Tresca stress for each row. + Array of shape (n, ...). Signed Tresca stress for each entry. Raises: ValueError: If input is not a 2D array with 6 columns. @@ -364,59 +439,16 @@ def calc_signed_tresca_by_max_abs_principal( tresca = calc_tresca_stress(stress_voigt) principals = calc_principal_stresses(stress_voigt) - avg_13 = 0.5 * (principals[:, 0] + principals[:, 2]) + avg_13 = 0.5 * (principals[..., 0] + principals[..., 2]) sign = np.sign(avg_13).astype(np.float64, copy=False) return sign * tresca -def calc_stress_deviator( - stress_voigt: NDArray[np.float64], -) -> NDArray[np.float64]: - r"""Calculate stress deviator for each stress state. - - ??? abstract "Math Equations" - $$ \mathbf{s} = \sigma - \frac{1}{3} tr(\sigma) $$ - - Args: - stress_voigt: Array of shape (n, 6). Each row is a stress vector in - Voigt notation. - - Returns: - Array of shape (n,). Signed Tresca stress for each row. - - Raises: - ValueError: If input is not a 2D array with 6 columns. - """ - voigt.check_shape(stress_voigt) - hydrostatic = calc_hydrostatic_stress(stress_voigt) - - deviator = stress_voigt.copy() - deviator[:, :3] -= hydrostatic[:, None] - - return deviator - - -if __name__ == "__main__": - # Example usage - stress_example = np.array( - [ - [10, 20, 30, 0, 0, 0], - [50, -50, 0, 0, 0, 10], - [0, 0, 1000, 0, 0, -50], - ] - ) - # print("Principal Stresses:", calc_principal_stresses(stress_example)) - # print("Hydrostatic Stress:", calc_hydrostatic_stress(stress_example)) - # print("Von Mises Stress:", calc_von_mises_stress(stress_example)) - # print("signed Von Mises Stress:", calc_signed_von_mises(stress_example)) - # print("Tresca Stress:", calc_tresca_stress(stress_example)) - # print("signed Tresca Stress:", calc_signed_tresca(stress_example)) - # eigvals, eigvecs = calc_principal_stresses_and_directions(stress_example) - # print("Principal Stresses:", eigvals) - # print("Principal Directions:", eigvecs) - # print("Principal Stresses:", eigvals) - # print("Principal Directions:", eigvecs) - # print("Principal Directions:", eigvecs) - # print("Principal Directions:", eigvecs) - # print("Principal Directions:", eigvecs) +# The function using np.linalg.eigh(tensor) is compatible with tensors constructed as: +# tensor[:, 0, 0, ...] = vector[:, 0, ...] +# tensor[:, 1, 1, ...] = vector[:, 1, ...] +# tensor[:, 2, 2, ...] = vector[:, 2, ...] +# tensor[:, [1, 2], [2, 1], ...] = vector[:, [3], ...] +# tensor[:, [0, 2], [2, 0], ...] = vector[:, [4], ...] +# tensor[:, [0, 1], [1, 0], ...] = vector[:, [5], ...] diff --git a/src/fatpy/utils/voigt.py b/src/fatpy/utils/voigt.py index 23cd4cb..1e9dbe5 100644 --- a/src/fatpy/utils/voigt.py +++ b/src/fatpy/utils/voigt.py @@ -1,38 +1,80 @@ -"""Helper functions for handling vectors using Voigt notation.""" +"""Helper functions for handling vectors using Voigt notation. + +All functions accept arrays with shape (n, 6, ...) where any trailing +dimensions after the 6-component Voigt axis are allowed. The 3x3 tensor +representation uses the last two axes for matrix indices, i.e. +returned tensors have shape (n, ..., 3, 3). +""" + +# TODO: N-Dimensional support for trailing dims +# - voigt, tensor dimensions have to be the last axes +# - (..., 6), (..., 3, 3) import numpy as np from numpy.typing import NDArray +VOIGT_COMPONENTS_COUNT = 6 +MIN_VECTOR_DIMS = 2 + -def check_shape(voigt_vec: NDArray[np.float64]) -> None: - """Check the shape of the input Voigt array. +def check_shape(vector: NDArray[np.float64]) -> None: + """Validate the Voigt vector shape. Args: - voigt_vec: Array of shape (n, 6). Each row is a vector in Voigt notation. + vector: Array with shape (n, 6, ...) where axis=1 has length 6. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If input does not have at least two dimensions or the + second dimension is not 6. """ - if voigt_vec.ndim != 2 or voigt_vec.shape[1] != 6: - raise ValueError("Input must be a n x 6 matrix in Voigt notation.") + if vector.ndim < MIN_VECTOR_DIMS: + raise ValueError("Input must be at least a 2D array (n, 6, ...).") + if vector.shape[1] != VOIGT_COMPONENTS_COUNT: + raise ValueError( + "Second dimension must correspond to 6 Voigt " + "stress/strain components (n, 6, ...)." + ) -def voigt_to_tensor(voigt_vec: NDArray[np.float64]) -> NDArray[np.float64]: + +def voigt_to_tensor(vector: NDArray[np.float64]) -> NDArray[np.float64]: """Convert Voigt vectors to symmetric 3x3 tensors. Args: - voigt_vec: Array of shape (n, 6). Each row is a vector in Voigt notation. - (σ11, σ22, σ33, σ23, σ13, σ12) + vector: Array of shape (n, 6, ...) where axis=1 contains the Voigt + components (σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy). Returns: - Array of shape (n, 3, 3). Symmetric tensors. + Array with shape (n, ..., 3, 3). The last two axes are the symmetric + 3x3 tensor for each input index. Trailing dimensions from the input + are preserved before the final two matrix axes. """ - n = voigt_vec.shape[0] - tensor = np.zeros((n, 3, 3), dtype=np.float64) - tensor[:, 0, 0] = voigt_vec[:, 0] - tensor[:, 1, 1] = voigt_vec[:, 1] - tensor[:, 2, 2] = voigt_vec[:, 2] - tensor[:, 1, 2] = tensor[:, 2, 1] = voigt_vec[:, 3] - tensor[:, 0, 2] = tensor[:, 2, 0] = voigt_vec[:, 4] - tensor[:, 0, 1] = tensor[:, 1, 0] = voigt_vec[:, 5] + check_shape(vector) + + # Move the Voigt component axis to the end for easy broadcasting, then + # construct tensors with matrix axes as the last two dims. + # Input shape: (n, 6, d1, d2, ...) + # We'll output shape: (n, d1, d2, ..., 3, 3) + leading = vector.shape[0] + trailing = vector.shape[2:] + + out_shape = (leading,) + trailing + (3, 3) + tensor = np.zeros(out_shape, dtype=vector.dtype) + + # Helper to index into input with preserved trailing dims + # Normal components + tensor[(..., 0, 0)] = vector[:, 0].reshape((leading,) + trailing) + tensor[(..., 1, 1)] = vector[:, 1].reshape((leading,) + trailing) + tensor[(..., 2, 2)] = vector[:, 2].reshape((leading,) + trailing) + + # Off-diagonal: Voigt indices 3 -> σ_yz, 4 -> σ_xz, 5 -> σ_xy + tensor[(..., 1, 2)] = vector[:, 3].reshape((leading,) + trailing) + tensor[(..., 2, 1)] = tensor[(..., 1, 2)] + + tensor[(..., 0, 2)] = vector[:, 4].reshape((leading,) + trailing) + tensor[(..., 2, 0)] = tensor[(..., 0, 2)] + + tensor[(..., 0, 1)] = vector[:, 5].reshape((leading,) + trailing) + tensor[(..., 1, 0)] = tensor[(..., 0, 1)] + return tensor diff --git a/tests/struct_mech/test_stress.py b/tests/struct_mech/test_stress.py new file mode 100644 index 0000000..b3bd6f6 --- /dev/null +++ b/tests/struct_mech/test_stress.py @@ -0,0 +1,248 @@ +"""Test functions for stress calculations in structural mechanics. + +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 stresses (descending) + +""" + +import numpy as np +import pytest +from numpy.typing import NDArray + +from fatpy.struct_mech import stress + + +@pytest.fixture +def stress_vector_2d() -> NDArray[np.float64]: + """Fixture providing sample stress vectors in Voigt notation. + + Returns: + NDArray[np.float64]: Sample stress vectors. + """ + # Provide a collection of representative stress states (n,6): + # 1) Uniaxial tension in x (sigma_xx = 100) + # 2) Uniaxial compression in z (sigma_zz = -50) + # 3) Pure hydrostatic (all normal = 30) + # 4) Pure shear (sxy = 40) + # 5) Mixed state + arr = np.array( + [ + [100.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, -50.0, 0.0, 0.0, 0.0], + [30.0, 30.0, 30.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 40.0], + [np.sqrt(2.0), -np.sqrt(2), 0.0, 0.0, 0.0, np.sqrt(2.0)], + ], + dtype=np.float64, + ) + + return arr + + +@pytest.fixture +def principal_stresses_2d() -> NDArray[np.float64]: + """Fixture providing expected principal stresses for the sample stress vectors. + + Returns: + NDArray[np.float64]: Expected principal stresses. + """ + arr = np.array( + [ + [100.0, 0.0, 0.0], + [0.0, 0.0, -50.0], + [30.0, 30.0, 30.0], + [40.0, 0.0, -40.0], + [2.0, 0.0, -2.0], + ], + dtype=np.float64, + ) + + return arr + + +# TODO: implement 3D cases +def stress_vector_3d(stress_vector_2d: NDArray[np.float64]) -> NDArray[np.float64]: + """Fixture providing sample stress vectors in Voigt notation. 3D case. + + This is similar to the 2D case but with additional axis possibly representing time. + + Returns: + NDArray[np.float64]: Sample stress vectors. + """ + # Expand to 3D by adding an additional axis (e.g., time) + return np.repeat(stress_vector_2d[np.newaxis, :, :], repeats=4, axis=0) + + +def test_calc_principal_stresses_and_directions_shape_2d( + stress_vector_2d: NDArray[np.float64], +) -> None: + """Test shape of principal stresses and directions output.""" + principals, directions = stress.calc_principal_stresses_and_directions( + stress_vector_2d + ) + assert principals.shape[0] == stress_vector_2d.shape[0] + assert principals.shape[1] == 3 + + assert directions.shape[0] == stress_vector_2d.shape[0] + assert directions.shape[1] == 3 + assert directions.shape[2] == 3 + + +def test_calc_principal_stresses_and_directions_ordering_2d( + stress_vector_2d: NDArray[np.float64], +) -> None: + """Test ordering of principal stresses (descending).""" + principals, dirctions = stress.calc_principal_stresses_and_directions( + stress_vector_2d + ) + assert np.all(principals[:, 0] >= principals[:, 1]) + assert np.all(principals[:, 1] >= principals[:, 2]) + + # test if directions order matches principal stresses order + # by checking A*v = λ*v for each principal stress/direction pair + for i in range(stress_vector_2d.shape[0]): + sigma = stress_vector_2d[i] + tensor = np.array( + [ + [sigma[0], sigma[5], sigma[4]], + [sigma[5], sigma[1], sigma[3]], + [sigma[4], sigma[3], sigma[2]], + ] + ) + for j in range(3): + eigvec = dirctions[i, :, j] + eigval = principals[i, j] + # Check if A*v = λ*v + Av = tensor @ eigvec + lv = eigval * eigvec + assert np.allclose(Av, lv, atol=1e-12) + + +def test_calc_principal_stresses_and_directions_values_2d( + stress_vector_2d: NDArray[np.float64], + principal_stresses_2d: NDArray[np.float64], +) -> None: + """Test specific known values of principal stresses.""" + principals, _ = stress.calc_principal_stresses_and_directions(stress_vector_2d) + assert np.allclose(principals, principal_stresses_2d, atol=1e-12) + + +def test_calc_principal_stresses_shape_2d( + stress_vector_2d: NDArray[np.float64], +) -> None: + """Test shape of principal stresses output.""" + principals = stress.calc_principal_stresses(stress_vector_2d) + assert principals.shape[0] == stress_vector_2d.shape[0] + assert principals.shape[1] == 3 + + +def test_calc_principal_stresses_ordering_2d( + stress_vector_2d: NDArray[np.float64], +) -> None: + """Test ordering of principal stresses (descending).""" + principals = stress.calc_principal_stresses(stress_vector_2d) + assert np.all(principals[:, 0] >= principals[:, 1]) + assert np.all(principals[:, 1] >= principals[:, 2]) + + +def test_calc_principal_stresses_values_2d( + stress_vector_2d: NDArray[np.float64], + principal_stresses_2d: NDArray[np.float64], +) -> None: + """Test specific known values of principal stresses.""" + principals = stress.calc_principal_stresses(stress_vector_2d) + assert np.allclose(principals, principal_stresses_2d, atol=1e-12) + + +def test_calc_principal_directions_shape_2d( + stress_vector_2d: NDArray[np.float64], +) -> None: + """Test shape of principal directions output.""" + directions = stress.calc_principal_directions(stress_vector_2d) + assert directions.shape[0] == stress_vector_2d.shape[0] + assert directions.shape[1] == 3 + assert directions.shape[2] == 3 + + +def test_calc_principal_directions_ordering_2d( + stress_vector_2d: NDArray[np.float64], +) -> None: + """Test ordering of principal directions matches principal stresses order.""" + principals = stress.calc_principal_stresses(stress_vector_2d) + directions = stress.calc_principal_directions(stress_vector_2d) + # For each stress state, check A*v = λ*v for each direction + for i in range(stress_vector_2d.shape[0]): + sigma = stress_vector_2d[i] + tensor = np.array( + [ + [sigma[0], sigma[5], sigma[4]], + [sigma[5], sigma[1], sigma[3]], + [sigma[4], sigma[3], sigma[2]], + ] + ) + for j in range(3): + eigvec = directions[i, :, j] + eigval = principals[i, j] + Av = tensor @ eigvec + lv = eigval * eigvec + assert np.allclose(Av, lv, atol=1e-12) + + +def test_invariants_and_hydrostatic_deviator( + stress_vector_2d: NDArray[np.float64], +) -> None: + invariants = stress.calc_stress_invariants(stress_vector_2d) + hydro = stress.calc_hydrostatic_stress(stress_vector_2d) + deviator = stress.calc_stress_deviator(stress_vector_2d) + # Shape checks + assert invariants.shape == (stress_vector_2d.shape[0], 3) + assert hydro.shape == (stress_vector_2d.shape[0],) + assert deviator.shape == stress_vector_2d.shape + # For hydrostatic case (index 2), deviator should be zero + assert np.allclose(deviator[2], 0.0, atol=1e-12) + # For hydrostatic case, invariants: I1 = 90, I2 = ? compute from definition + I1 = invariants[2, 0] + assert np.isclose(I1, 90.0) + + +def test_von_mises_and_signed_variants(stress_vector_2d: NDArray[np.float64]) -> None: + vm = stress.calc_von_mises_stress(stress_vector_2d) + svm_h = stress.calc_signed_von_mises_by_hydrostatic(stress_vector_2d) + svm_p = stress.calc_signed_von_mises_by_max_abs_principal(stress_vector_2d) + svm_i1 = stress.calc_signed_von_mises_by_first_invariant(stress_vector_2d) + # Shape + assert vm.shape == (stress_vector_2d.shape[0],) + # Non-negative for von Mises + assert np.all(vm >= 0.0) + # For hydrostatic pure case index 2, von Mises == 0 and signed variants 0 + assert np.isclose(vm[2], 0.0) + assert np.isclose(svm_h[2], 0.0) + assert np.isclose(svm_p[2], 0.0) + assert np.isclose(svm_i1[2], 0.0) + # For pure shear case index 3 (sxy = 40): von Mises = sqrt(3/2)*|shear|*? check + # Compute expected manually from definition used in implementation + sx, sy, sz = stress_vector_2d[3, :3] + sxy = stress_vector_2d[3, 5] + expected_vm3 = np.sqrt( + 0.5 * ((sx - sy) ** 2 + (sy - sz) ** 2 + (sz - sx) ** 2 + 6 * (sxy**2)) + ) + assert np.isclose(vm[3], expected_vm3) + + +def test_tresca_and_signed_variants(stress_vector_2d: NDArray[np.float64]) -> None: + tresca = stress.calc_tresca_stress(stress_vector_2d) + stresca_h = stress.calc_signed_tresca_by_hydrostatic(stress_vector_2d) + stresca_p = stress.calc_signed_tresca_by_max_abs_principal(stress_vector_2d) + # Shape + assert tresca.shape == (stress_vector_2d.shape[0],) + # Hydrostatic case (index 2) tresca == 0 + assert np.isclose(tresca[2], 0.0) + assert np.isclose(stresca_h[2], 0.0) + # For uniaxial tension (index 0), principal stresses should include 100 and two zeros -> tresca = (100 - 0)/2 + assert np.isclose(tresca[0], 50.0) From 7f9f208c55bd2b9397247c7fa4ab4fb19bd0322a Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Tue, 23 Sep 2025 21:57:16 +0200 Subject: [PATCH 09/17] wip - rewriting to n dimensional 2 --- src/fatpy/struct_mech/stress.py | 7 +-- src/fatpy/utils/voigt.py | 84 ++++++++++++++++----------------- 2 files changed, 44 insertions(+), 47 deletions(-) diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index 7bfa259..ad2c159 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -50,14 +50,15 @@ def calc_principal_stresses_and_directions( Returns: Tuple (eigvals, eigvecs): - eigvals: Array of shape (n, 3, ...). Principal stresses - (descending: σ_1 ≥ σ_2 ≥ σ_3) with trailing dims preserved. + (descending: σ_1 ≥ σ_2 ≥ σ_3) with trailing dims preserved. - eigvecs: Array of shape (n, 3, 3, ...). Principal directions (columns are - eigenvectors) aligned with eigvals in the same order. The last two - axes of this array are the 3x3 eigenvector matrix for each input. + eigenvectors) aligned with eigvals in the same order. The last two + axes of this array are the 3x3 eigenvector matrix for each input. Raises: ValueError: If input is not a 2D array with 6 columns. """ + # TODO: check docsrtring agains `MkDocs` formatting voigt.check_shape(stress_voigt) tensor = voigt.voigt_to_tensor(stress_voigt) diff --git a/src/fatpy/utils/voigt.py b/src/fatpy/utils/voigt.py index 1e9dbe5..5aaf622 100644 --- a/src/fatpy/utils/voigt.py +++ b/src/fatpy/utils/voigt.py @@ -1,39 +1,47 @@ """Helper functions for handling vectors using Voigt notation. -All functions accept arrays with shape (n, 6, ...) where any trailing -dimensions after the 6-component Voigt axis are allowed. The 3x3 tensor -representation uses the last two axes for matrix indices, i.e. -returned tensors have shape (n, ..., 3, 3). -""" +Overview: + This module provides utilities for converting between Voigt notation vectors + and their corresponding symmetric 3x3 tensor representations. Voigt notation + is commonly used in continuum mechanics to represent stress and strain tensors + in a compact vector form. + +Supported Shapes: + - Input arrays: (..., 6) + The last axis contains the 6 Voigt components: + (σ_11, σ_22, σ_33, σ_23, σ_13, σ_12) + or equivalently (σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy) + - Output arrays: (..., 3, 3) + The last two axes represent the symmetric 3x3 tensor for each input index. + Any number of leading dimensions is supported and preserved. + +Usage: + Use these functions to convert between Voigt vectors and tensor representations + for stress or strain, ensuring compatibility with numerical routines that expect + either format. -# TODO: N-Dimensional support for trailing dims -# - voigt, tensor dimensions have to be the last axes -# - (..., 6), (..., 3, 3) +""" import numpy as np from numpy.typing import NDArray VOIGT_COMPONENTS_COUNT = 6 -MIN_VECTOR_DIMS = 2 def check_shape(vector: NDArray[np.float64]) -> None: """Validate the Voigt vector shape. Args: - vector: Array with shape (n, 6, ...) where axis=1 has length 6. + vector: Array with shape (..., 6) where the last dimension has length 6. Raises: ValueError: If input does not have at least two dimensions or the second dimension is not 6. """ - if vector.ndim < MIN_VECTOR_DIMS: - raise ValueError("Input must be at least a 2D array (n, 6, ...).") - - if vector.shape[1] != VOIGT_COMPONENTS_COUNT: + if vector.shape[-1] != VOIGT_COMPONENTS_COUNT: raise ValueError( - "Second dimension must correspond to 6 Voigt " - "stress/strain components (n, 6, ...)." + "Last dimension must correspond to 6 Voigt " + "stress/strain components (..., 6)." ) @@ -41,40 +49,28 @@ def voigt_to_tensor(vector: NDArray[np.float64]) -> NDArray[np.float64]: """Convert Voigt vectors to symmetric 3x3 tensors. Args: - vector: Array of shape (n, 6, ...) where axis=1 contains the Voigt - components (σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy). + vector: Array of shape (..., 6) where the last dimension contains the + stress/strain components in order according to Voigt notation: + (σ_11, σ_22, σ_33, σ_23, σ_13, σ_12) + (σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy) Returns: - Array with shape (n, ..., 3, 3). The last two axes are the symmetric - 3x3 tensor for each input index. Trailing dimensions from the input - are preserved before the final two matrix axes. + Array with shape (..., 3, 3). The last two axes are the symmetric + 3x3 tensor for each input index. Trailing dimensions are preserved. """ check_shape(vector) - # Move the Voigt component axis to the end for easy broadcasting, then - # construct tensors with matrix axes as the last two dims. - # Input shape: (n, 6, d1, d2, ...) - # We'll output shape: (n, d1, d2, ..., 3, 3) - leading = vector.shape[0] - trailing = vector.shape[2:] - - out_shape = (leading,) + trailing + (3, 3) - tensor = np.zeros(out_shape, dtype=vector.dtype) + shape = vector.shape[:-1] + (3, 3) + tensor = np.zeros(shape, dtype=vector.dtype) - # Helper to index into input with preserved trailing dims # Normal components - tensor[(..., 0, 0)] = vector[:, 0].reshape((leading,) + trailing) - tensor[(..., 1, 1)] = vector[:, 1].reshape((leading,) + trailing) - tensor[(..., 2, 2)] = vector[:, 2].reshape((leading,) + trailing) - - # Off-diagonal: Voigt indices 3 -> σ_yz, 4 -> σ_xz, 5 -> σ_xy - tensor[(..., 1, 2)] = vector[:, 3].reshape((leading,) + trailing) - tensor[(..., 2, 1)] = tensor[(..., 1, 2)] - - tensor[(..., 0, 2)] = vector[:, 4].reshape((leading,) + trailing) - tensor[(..., 2, 0)] = tensor[(..., 0, 2)] - - tensor[(..., 0, 1)] = vector[:, 5].reshape((leading,) + trailing) - tensor[(..., 1, 0)] = tensor[(..., 0, 1)] + tensor[(..., 0, 0)] = vector[..., 0] # xx + tensor[(..., 1, 1)] = vector[..., 1] # yy + tensor[(..., 2, 2)] = vector[..., 2] # zz + + # Shear components + tensor[(..., [1, 2], [2, 1])] = vector[..., [3]] # yz + tensor[(..., [0, 2], [2, 0])] = vector[..., [4]] # xz + tensor[(..., [0, 1], [1, 0])] = vector[..., [5]] # xy return tensor From a31628527fb58f7a77f57ef9f7920646e2df3c10 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Wed, 24 Sep 2025 02:14:18 +0200 Subject: [PATCH 10/17] utils - voigt tested --- .gitignore | 3 + .pre-commit-config.yaml | 5 +- pyproject.toml | 36 +++++---- tests/utils/test_voigt.py | 161 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 189 insertions(+), 16 deletions(-) create mode 100644 tests/utils/test_voigt.py diff --git a/.gitignore b/.gitignore index b045860..c59c7c7 100644 --- a/.gitignore +++ b/.gitignore @@ -42,3 +42,6 @@ htmlcov/ # Generated by setuptools_scm src/fatpy/_version.py + +# Pre-commit hooks +.cache_pre_commit/ diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6731f47..b77d127 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -8,7 +8,7 @@ repos: - id: check-added-large-files - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.12.9 + rev: v0.13.1 hooks: - id: ruff args: @@ -17,6 +17,7 @@ repos: - id: ruff-format - repo: https://github.com/pre-commit/mirrors-mypy - rev: v1.17.1 + rev: v1.18.2 hooks: - id: mypy + files: ^src/ diff --git a/pyproject.toml b/pyproject.toml index 0558a77..c634b10 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,9 +18,9 @@ classifiers = [ ] dependencies = [ - "matplotlib>=3.10.0", - "numpy>=2.2.2", - "pandas>=2.2.3", + "matplotlib>=3.10.6", + "numpy>=2.3.3", + "pandas>=2.3.2", # "scipy>=1.15.1", ] @@ -33,18 +33,18 @@ dependencies = [ [tool.uv] dev-dependencies = [ - "mypy>=1.15.0", - "ruff>=0.11.3", - "pre-commit>=4.2.0", - "pytest>=8.3.5", - "pytest-cov>=6.0.0", - "pytest-mock>=3.14.0", + "mypy>=1.18.2", + "ruff>=0.13.1", + "pre-commit>=4.3.0", + "pytest>=8.4.2", + "pytest-cov>=7.0.0", + "pytest-mock>=3.15.1", "pytest-mypy>=1.0.1", - "pytest-ruff>=0.4.1", + "pytest-ruff>=0.5", "mkdocs>=1.6.1", - "mkdocs-material[imaging]>=9.5.15", - "mkdocstrings[python]>=0.29.1", - "mkdocs-autorefs>=1.0.0", + "mkdocs-material[imaging]>=9.6.20", + "mkdocstrings[python]>=0.30.1", + "mkdocs-autorefs>=1.4.3", ] [tool.pytest.ini_options] @@ -129,12 +129,16 @@ docstring-code-line-length = "dynamic" [tool.ruff.lint.per-file-ignores] "__init__.py" = ["E402"] +"tests/**" = [ + "D", # pydocstyle + "B", # bugbear + "S", # bandit +] [tool.mypy] python_version = "3.12" files = "src" strict = true -ignore_missing_imports = true check_untyped_defs = true disallow_subclassing_any = true warn_unreachable = true @@ -142,6 +146,10 @@ warn_unused_configs = true show_error_context = true show_column_numbers = true +[[tool.mypy.overrides]] +module = ["tests", "tests.*"] +disallow_untyped_decorators = false + [build-system] requires = ["setuptools>=61.0", "setuptools_scm[toml]>=6.2", "wheel"] build-backend = "setuptools.build_meta" diff --git a/tests/utils/test_voigt.py b/tests/utils/test_voigt.py new file mode 100644 index 0000000..daf7493 --- /dev/null +++ b/tests/utils/test_voigt.py @@ -0,0 +1,161 @@ +import numpy as np +import pytest +from numpy.typing import NDArray + +from fatpy.utils import voigt + + +def get_tensor(array: list[int]) -> NDArray[np.float64]: + """Helper to create a tensor from a nested list.""" + return np.array( + [ + [array[0], array[5], array[4]], + [array[5], array[1], array[3]], + [array[4], array[3], array[2]], + ], + dtype=np.float64, + ) + + +@pytest.fixture +def vector_1d() -> NDArray[np.float64]: + """Fixture providing a 1D Voigt vector. + + Returns: + NDArray[np.float64]: 1D Voigt vector. + """ + return np.arange(1, 7, dtype=np.float64) + + +@pytest.fixture +def vector_2d() -> NDArray[np.float64]: + vector = np.arange(1, 13, dtype=np.float64) + return vector.reshape((2, 6)) + + +@pytest.fixture +def vector_3d() -> NDArray[np.float64]: + vector = np.arange(1, 25, dtype=np.float64) + return vector.reshape((2, 2, 6)) + + +@pytest.fixture +def vector_4d() -> NDArray[np.float64]: + vector = np.arange(1, 49, dtype=np.float64) + return vector.reshape((2, 2, 2, 6)) + + +@pytest.mark.parametrize( + "shape,raises_value_error", + [ + ((6,), False), + ((2, 6), False), + ((2, 2, 6), False), + ((2, 2, 2, 6), False), + ((2, 2, 2, 2, 6), False), + ((2,), True), + ((6, 2), True), + ((6, 2, 2), True), + ((6, 2, 2, 2), True), + ((6, 2, 2, 2, 2), True), + ], +) +def test_check_shape(shape: tuple[int, ...], raises_value_error: bool) -> None: + """Test the shape validation function for Voigt vectors.""" + if not raises_value_error: + voigt.check_shape(np.ones(shape, dtype=np.float64)) + + else: + with pytest.raises(ValueError): + voigt.check_shape(np.ones(shape, dtype=np.float64)) + + +@pytest.mark.parametrize( + "shape", + [ + ((6,)), + ((2, 6)), + ((2, 2, 6)), + ((2, 2, 2, 6)), + ((2, 2, 2, 2, 6)), + ], +) +def test_voigt_to_tensor_shape(shape: tuple[int, ...]) -> None: + """Test the shape validation function for Voigt vectors.""" + tensor = voigt.voigt_to_tensor(np.ones(shape, dtype=np.float64)) + expected_shape = shape[:-1] + (3, 3) + assert tensor.shape == expected_shape + + +def test_voigt_to_tensor_1d(vector_1d: NDArray[np.float64]) -> None: + """Test conversion from 1D Voigt vector to tensor.""" + tensor = voigt.voigt_to_tensor(vector_1d) + expected = get_tensor(list(range(1, 7))) + np.testing.assert_array_equal(tensor, expected) + assert tensor.shape == (3, 3) + + +def test_voigt_to_tensor_2d(vector_2d: NDArray[np.float64]) -> None: + """Test conversion from 2D Voigt vector to tensor.""" + tensor = voigt.voigt_to_tensor(vector_2d) + expected = np.array( + [ + get_tensor(list(range(1, 7))), + get_tensor(list(range(7, 13))), + ], + dtype=np.float64, + ) + np.testing.assert_array_equal(tensor, expected) + assert tensor.shape == (2, 3, 3) + + +def test_voigt_to_tensor_3d(vector_3d: NDArray[np.float64]) -> None: + """Test conversion from 3D Voigt vector to tensor.""" + tensor = voigt.voigt_to_tensor(vector_3d) + expected = np.array( + [ + [ + get_tensor(list(range(1, 7))), + get_tensor(list(range(7, 13))), + ], + [ + get_tensor(list(range(13, 19))), + get_tensor(list(range(19, 25))), + ], + ], + dtype=np.float64, + ) + np.testing.assert_array_equal(tensor, expected) + assert tensor.shape == (2, 2, 3, 3) + + +def test_voigt_to_tensor_4d(vector_4d: NDArray[np.float64]) -> None: + """Test conversion from 4D Voigt vector to tensor.""" + tensor = voigt.voigt_to_tensor(vector_4d) + expected = np.array( + [ + [ + [ + get_tensor(list(range(1, 7))), + get_tensor(list(range(7, 13))), + ], + [ + get_tensor(list(range(13, 19))), + get_tensor(list(range(19, 25))), + ], + ], + [ + [ + get_tensor(list(range(25, 31))), + get_tensor(list(range(31, 37))), + ], + [ + get_tensor(list(range(37, 43))), + get_tensor(list(range(43, 49))), + ], + ], + ], + dtype=np.float64, + ) + np.testing.assert_array_equal(tensor, expected) + assert tensor.shape == (2, 2, 2, 3, 3) From c5d4b66c21245a26af8c95ecd9a01b4ccbf9b2bc Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Wed, 24 Sep 2025 09:19:56 +0200 Subject: [PATCH 11/17] mypy settings adjusted --- .pre-commit-config.yaml | 1 - pyproject.toml | 3 --- 2 files changed, 4 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b77d127..535c7a6 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -20,4 +20,3 @@ repos: rev: v1.18.2 hooks: - id: mypy - files: ^src/ diff --git a/pyproject.toml b/pyproject.toml index c634b10..f3c6ff3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -145,9 +145,6 @@ warn_unreachable = true warn_unused_configs = true show_error_context = true show_column_numbers = true - -[[tool.mypy.overrides]] -module = ["tests", "tests.*"] disallow_untyped_decorators = false [build-system] From cf59567ffe5a9f5634485ceffce8c3e7978f8346 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Sun, 28 Sep 2025 17:27:25 +0200 Subject: [PATCH 12/17] voigt docstrings adjusted --- src/fatpy/utils/voigt.py | 40 ++++++++++++++++++---------------------- 1 file changed, 18 insertions(+), 22 deletions(-) diff --git a/src/fatpy/utils/voigt.py b/src/fatpy/utils/voigt.py index 5aaf622..7954007 100644 --- a/src/fatpy/utils/voigt.py +++ b/src/fatpy/utils/voigt.py @@ -1,24 +1,17 @@ -"""Helper functions for handling vectors using Voigt notation. +"""Tools for working with Voigt notation in continuum mechanics. Overview: - This module provides utilities for converting between Voigt notation vectors - and their corresponding symmetric 3x3 tensor representations. Voigt notation - is commonly used in continuum mechanics to represent stress and strain tensors - in a compact vector form. - -Supported Shapes: - - Input arrays: (..., 6) - The last axis contains the 6 Voigt components: - (σ_11, σ_22, σ_33, σ_23, σ_13, σ_12) - or equivalently (σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy) - - Output arrays: (..., 3, 3) - The last two axes represent the symmetric 3x3 tensor for each input index. - Any number of leading dimensions is supported and preserved. - -Usage: - Use these functions to convert between Voigt vectors and tensor representations - for stress or strain, ensuring compatibility with numerical routines that expect - either format. + This module provides general utilities for handling vectors and tensors + represented in Voigt notation. Voigt notation is widely used to express + symmetric 3x3 tensors, such as stress and strain, in a compact vector form. + +Input Shape Convention: + - Arrays are expected to have shape (..., 6), where the last dimension + contains the six Voigt components: + (σ_11, σ_22, σ_33, σ_23, σ_13, σ_12) + (σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy) + Leading dimensions (the '...' part) are preserved and can be arbitrary, + allowing batch processing of multiple vectors or tensors. """ @@ -35,8 +28,7 @@ def check_shape(vector: NDArray[np.float64]) -> None: vector: Array with shape (..., 6) where the last dimension has length 6. Raises: - ValueError: If input does not have at least two dimensions or the - second dimension is not 6. + ValueError: If the last dimension is not of size 6. """ if vector.shape[-1] != VOIGT_COMPONENTS_COUNT: raise ValueError( @@ -56,10 +48,14 @@ def voigt_to_tensor(vector: NDArray[np.float64]) -> NDArray[np.float64]: Returns: Array with shape (..., 3, 3). The last two axes are the symmetric - 3x3 tensor for each input index. Trailing dimensions are preserved. + 3x3 tensor for each input index. Leading dimensions are preserved. + + Raises: + ValueError: If the last dimension is not of size 6. """ check_shape(vector) + # (1e6, 50, 8, 6) -> (1e6, 50, 8) + (3, 3) -> (1e6, 50, 8, 3, 3) shape = vector.shape[:-1] + (3, 3) tensor = np.zeros(shape, dtype=vector.dtype) From aeec674591ff89468c783ef4206c06584ede37f8 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Sun, 28 Sep 2025 21:02:24 +0200 Subject: [PATCH 13/17] added tensor to voigt --- src/fatpy/utils/voigt.py | 36 ++++++++++++++++++++++++++++++++++++ tests/utils/test_voigt.py | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 68 insertions(+) diff --git a/src/fatpy/utils/voigt.py b/src/fatpy/utils/voigt.py index 7954007..bcd269c 100644 --- a/src/fatpy/utils/voigt.py +++ b/src/fatpy/utils/voigt.py @@ -70,3 +70,39 @@ def voigt_to_tensor(vector: NDArray[np.float64]) -> NDArray[np.float64]: tensor[(..., [0, 1], [1, 0])] = vector[..., [5]] # xy return tensor + + +def tensor_to_voigt(tensor: NDArray[np.float64]) -> NDArray[np.float64]: + """Convert symmetric 3x3 tensors to Voigt vectors. + + Args: + tensor: Array of shape (..., 3, 3) where the last two dimensions + represent symmetric 3x3 tensors. + + Returns: + Array with shape (..., 6). The last dimension contains the stress/strain + components in order according to Voigt notation: + (σ_11, σ_22, σ_33, σ_23, σ_13, σ_12) + (σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy) + Leading dimensions are preserved. + + Raises: + ValueError: If the last two dimensions are not of size (3, 3). + """ + if tensor.shape[-2:] != (3, 3): + raise ValueError("Last two dimensions must correspond to (3, 3) tensors.") + + shape = tensor.shape[:-2] + (VOIGT_COMPONENTS_COUNT,) + array = np.zeros(shape, dtype=tensor.dtype) + + # Normal components + array[..., 0] = tensor[..., 0, 0] # xx + array[..., 1] = tensor[..., 1, 1] # yy + array[..., 2] = tensor[..., 2, 2] # zz + + # Shear components + array[..., 3] = tensor[..., 1, 2] # yz + array[..., 4] = tensor[..., 0, 2] # xz + array[..., 5] = tensor[..., 0, 1] # xy + + return array diff --git a/tests/utils/test_voigt.py b/tests/utils/test_voigt.py index daf7493..0ee5c6b 100644 --- a/tests/utils/test_voigt.py +++ b/tests/utils/test_voigt.py @@ -159,3 +159,35 @@ def test_voigt_to_tensor_4d(vector_4d: NDArray[np.float64]) -> None: ) np.testing.assert_array_equal(tensor, expected) assert tensor.shape == (2, 2, 2, 3, 3) + + +def test_tensor_to_voigt_1d(vector_1d: NDArray[np.float64]) -> None: + """Test conversion from 1D tensor to Voigt vector.""" + tensor = voigt.voigt_to_tensor(vector_1d) + voigt_vector = voigt.tensor_to_voigt(tensor) + np.testing.assert_array_equal(voigt_vector, vector_1d) + assert voigt_vector.shape == (6,) + + +def test_tensor_to_voigt_2d(vector_2d: NDArray[np.float64]) -> None: + """Test conversion from 2D tensor to Voigt vector.""" + tensor = voigt.voigt_to_tensor(vector_2d) + voigt_vector = voigt.tensor_to_voigt(tensor) + np.testing.assert_array_equal(voigt_vector, vector_2d) + assert voigt_vector.shape == (2, 6) + + +def test_tensor_to_voigt_3d(vector_3d: NDArray[np.float64]) -> None: + """Test conversion from 3D tensor to Voigt vector.""" + tensor = voigt.voigt_to_tensor(vector_3d) + voigt_vector = voigt.tensor_to_voigt(tensor) + np.testing.assert_array_equal(voigt_vector, vector_3d) + assert voigt_vector.shape == (2, 2, 6) + + +def test_tensor_to_voigt_4d(vector_4d: NDArray[np.float64]) -> None: + """Test conversion from 4D tensor to Voigt vector.""" + tensor = voigt.voigt_to_tensor(vector_4d) + voigt_vector = voigt.tensor_to_voigt(tensor) + np.testing.assert_array_equal(voigt_vector, vector_4d) + assert voigt_vector.shape == (2, 2, 2, 6) From 32cfc1d448ea0074e65fd81e4abb58216a25aa10 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Mon, 29 Sep 2025 01:24:28 +0200 Subject: [PATCH 14/17] stress functions reworked, tests implemented --- src/fatpy/py.typed | 0 src/fatpy/struct_mech/stress.py | 294 +++++++++++---------- tests/struct_mech/test_stress.py | 435 ++++++++++++++++++++----------- 3 files changed, 423 insertions(+), 306 deletions(-) create mode 100644 src/fatpy/py.typed diff --git a/src/fatpy/py.typed b/src/fatpy/py.typed new file mode 100644 index 0000000..e69de29 diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index ad2c159..c7269eb 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -6,24 +6,21 @@ both uniaxial and multiaxial loading conditions. Conventions: -- Vectors use Voigt notation with shape (n, 6, ...). - - For stress tensors, the six Voigt components are: +- 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 stresses are ordered in descending order throughout the module -(σ_1 ≥ σ_2 ≥ σ_3). + (σ_1 ≥ σ_2 ≥ σ_3). - Principal directions (eigenvectors) are aligned to this ordering -(columns correspond to σ_1, σ_2, σ_3). + (columns correspond to σ_1, σ_2, σ_3). """ -# TODO: N-Dimensional support for trailing dims -# - voigt, tensor dimensions have to be the last axes -# - (..., 6), (..., 3, 3) +# TODO: implement atol, rtol to signed stress functions import numpy as np from numpy.typing import NDArray @@ -32,7 +29,7 @@ def calc_principal_stresses_and_directions( - stress_voigt: NDArray[np.float64], + stress_vector_voigt: NDArray[np.float64], ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: r"""Calculate principal stresses and principal directions for each state. @@ -43,40 +40,26 @@ def calc_principal_stresses_and_directions( $$ \sigma \mathbf{v} = \lambda \mathbf{v} $$ Args: - stress_voigt: Array of shape (n, 6, ...). The second axis of size 6 - contains the Voigt stress components. Any trailing dimensions are - preserved and the returned arrays will have those trailing dims. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: Tuple (eigvals, eigvecs): - - eigvals: Array of shape (n, 3, ...). Principal stresses - (descending: σ_1 ≥ σ_2 ≥ σ_3) with trailing dims preserved. - - eigvecs: Array of shape (n, 3, 3, ...). Principal directions (columns are - eigenvectors) aligned with eigvals in the same order. The last two - axes of this array 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 input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - # TODO: check docsrtring agains `MkDocs` formatting - voigt.check_shape(stress_voigt) - - tensor = voigt.voigt_to_tensor(stress_voigt) - - # tensor has shape (n, ..., 3, 3). Use eigh over the last two axes. - eigvals, eigvecs = np.linalg.eigh(tensor) # eigvals shape: (n, ..., 3) + voigt.check_shape(stress_vector_voigt) - # Sort eigenvalues descending along the last axis and reorder eigenvectors - # accordingly. eigvecs has shape (n, ..., 3, 3) where the last axis is - # the eigenvector corresponding to the eigenvalue at the same position. + tensor = voigt.voigt_to_tensor(stress_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) - - # For eigvecs we need to take along the last axis (eigenvector index) - # while preserving the matrix axis. Move axes so we can index easily: - # eigvecs currently: (..., 3, 3) where the last axis indexes components - # and the second-last axis indexes eigenvector number. We want to reorder - # the eigenvector-number axis. eigvecs_sorted = np.take_along_axis( eigvecs, np.expand_dims(sorted_indices, axis=-2), axis=-1 ) @@ -84,7 +67,9 @@ def calc_principal_stresses_and_directions( return eigvals_sorted, eigvecs_sorted -def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: +def calc_principal_stresses( + stress_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: r"""Calculate principal stresses for each stress state. ??? abstract "Math Equations" @@ -94,26 +79,27 @@ def calc_principal_stresses(stress_voigt: NDArray[np.float64]) -> NDArray[np.flo $$ \sigma \mathbf{v} = \lambda \mathbf{v} $$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, 3, ...). Principal stresses (descending: σ1 ≥ σ2 ≥ σ3). + Array of shape (..., 3). Principal stresses (descending: σ1 ≥ σ2 ≥ σ3). Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(stress_voigt) + voigt.check_shape(stress_vector_voigt) - tensor = voigt.voigt_to_tensor(stress_voigt) + tensor = voigt.voigt_to_tensor(stress_vector_voigt) eigvals = np.linalg.eigvalsh(tensor) # eigvals shape: (n, ..., 3). Return sorted descending along last axis. return np.sort(eigvals, axis=-1)[..., ::-1] -def calc_principal_directions(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: +def calc_principal_directions( + stress_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: r"""Calculate principal directions (eigenvectors) for each stress state. ??? abstract "Math Equations" @@ -123,25 +109,26 @@ def calc_principal_directions(stress_voigt: NDArray[np.float64]) -> NDArray[np.f $$ \sigma \mathbf{v} = \lambda \mathbf{v} $$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, 3, 3, ...). Principal directions (columns are eigenvectors) + Array of shape (..., 3, 3). Principal directions (columns are eigenvectors) aligned with descending principal stresses: σ1, σ2, σ3. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(stress_voigt) + voigt.check_shape(stress_vector_voigt) # Reuse the ordering logic from the combined function to ensure consistency - _, eigvecs = calc_principal_stresses_and_directions(stress_voigt) + _, eigvecs = calc_principal_stresses_and_directions(stress_vector_voigt) return eigvecs -def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: +def calc_stress_invariants( + stress_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: r"""Calculate the first, second, and third invariants for each stress state. ??? abstract "Math Equations" @@ -154,19 +141,19 @@ def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.floa $$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, 3, ...). Columns are (I1, I2, I3) for each entry. + Array of shape (..., 3). The last dimension contains (I1, I2, I3) for + each entry. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(stress_voigt) + voigt.check_shape(stress_vector_voigt) - tensor = voigt.voigt_to_tensor(stress_voigt) + tensor = voigt.voigt_to_tensor(stress_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) @@ -176,7 +163,9 @@ def calc_stress_invariants(stress_voigt: NDArray[np.float64]) -> NDArray[np.floa return np.stack((invariant_1, invariant_2, invariant_3), axis=-1) -def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: +def calc_hydrostatic_stress( + stress_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: r"""Calculate the hydrostatic (mean normal) stress for each stress state. ??? abstract "Math Equations" @@ -186,26 +175,28 @@ def calc_hydrostatic_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.flo $$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, ...). Hydrostatic stress for each input state. + Array of shape (...). Hydrostatic stress for each input state. + Tensor rank is reduced by one. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(stress_voigt) + voigt.check_shape(stress_vector_voigt) - # Voigt normal components are at axis=1 indices 0,1,2. Preserve trailing dims. + # Voigt normal components are at last axis indices 0,1,2 return ( - stress_voigt[:, 0, ...] + stress_voigt[:, 1, ...] + stress_voigt[:, 2, ...] + stress_vector_voigt[..., 0] + + stress_vector_voigt[..., 1] + + stress_vector_voigt[..., 2] ) / 3.0 def calc_stress_deviator( - stress_voigt: NDArray[np.float64], + stress_vector_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: r"""Calculate stress deviator for each stress state. @@ -213,28 +204,29 @@ def calc_stress_deviator( $$ \mathbf{s} = \sigma - \frac{1}{3} tr(\sigma) $$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, 6, ...). Stress deviator for each entry. + Array of shape (..., 6). Stress deviator for each entry. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(stress_voigt) - hydrostatic = calc_hydrostatic_stress(stress_voigt) + voigt.check_shape(stress_vector_voigt) + hydrostatic = calc_hydrostatic_stress(stress_vector_voigt) - deviator = stress_voigt.copy() - # Subtract hydrostatic from the first three Voigt components (axis=1) - deviator[:, :3, ...] = deviator[:, :3, ...] - hydrostatic[..., None] + deviator = stress_vector_voigt.copy() + # Subtract hydrostatic from the first three Voigt components (last axis) + deviator[..., :3] = deviator[..., :3] - hydrostatic[..., None] return deviator # Von Mises functions -def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: +def calc_von_mises_stress( + stress_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: r"""Calculate von Mises equivalent stress for each stress state. ??? abstract "Math Equations" @@ -247,25 +239,26 @@ def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float $$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, ...). Von Mises equivalent stress for each entry. + Array of shape (...). Von Mises equivalent stress for each entry. + Tensor rank is reduced by one. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - voigt.check_shape(stress_voigt) + voigt.check_shape(stress_vector_voigt) - sx = stress_voigt[:, 0, ...] - sy = stress_voigt[:, 1, ...] - sz = stress_voigt[:, 2, ...] - syz = stress_voigt[:, 3, ...] - sxz = stress_voigt[:, 4, ...] - sxy = stress_voigt[:, 5, ...] + sx = stress_vector_voigt[..., 0] + sy = stress_vector_voigt[..., 1] + sz = stress_vector_voigt[..., 2] + syz = stress_vector_voigt[..., 3] + sxz = stress_vector_voigt[..., 4] + sxy = stress_vector_voigt[..., 5] + # Von Mises formula expanded to simplify computation return np.sqrt( 0.5 * ( @@ -278,7 +271,7 @@ def calc_von_mises_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float def calc_signed_von_mises_by_hydrostatic( - stress_voigt: NDArray[np.float64], + stress_vector_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: r"""Calculate signed von Mises stress for each stress state. @@ -288,18 +281,22 @@ def calc_signed_von_mises_by_hydrostatic( $$\sigma_{SvM} = sgn(\sigma_H) \cdot \sigma_{vM}$$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, ...). Signed von Mises stress for each entry. + Array of shape (...). Signed von Mises stress for each entry. + Tensor rank is reduced by one. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - von_mises = calc_von_mises_stress(stress_voigt) - sign = np.sign(calc_hydrostatic_stress(stress_voigt)) + von_mises = calc_von_mises_stress(stress_vector_voigt) + hydrostatic_stress = calc_hydrostatic_stress(stress_vector_voigt) + + sign = np.sign(hydrostatic_stress).astype(np.float64, copy=False) + sign[np.isclose(hydrostatic_stress, 0)] = 1.0 + return sign * von_mises @@ -308,33 +305,35 @@ def calc_signed_von_mises_by_max_abs_principal( ) -> NDArray[np.float64]: r"""Calculate signed von Mises stress for each stress state. - Sign is determined by the maximum absolute principal stress value. + Sign is determined by average of the maximum and minimum principal stresses. + In case the maximum absolute principal stress is zero, the sign is set to +1. + In case the maximum absolute principal stress is equal to negative value of + the minimum principal stress, the sign is set to +1 as well. ??? abstract "Math Equations" - $$\sigma_{SvM} = sgn(\frac{\sigma_{1}+\sigma_{3}}{2}) \cdot \sigma_{vM}$$ + $$\sigma_{SvM} = sgn(\sigma_{max,abs}) \cdot \sigma_{vM}$$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, ...). Signed von Mises stress for each entry. + Array of shape (...). Signed von Mises stress for each entry. + Tensor rank is reduced by one. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ von_mises = calc_von_mises_stress(stress_voigt) principals = calc_principal_stresses(stress_voigt) - avg_13 = 0.5 * (principals[..., 0] + principals[..., 2]) sign = np.sign(avg_13).astype(np.float64, copy=False) - + sign[np.isclose(avg_13, 0)] = 1.0 return sign * von_mises def calc_signed_von_mises_by_first_invariant( - stress_voigt: NDArray[np.float64], + stress_vector_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: r"""Calculate signed von Mises stress for each stress state. @@ -344,51 +343,54 @@ def calc_signed_von_mises_by_first_invariant( $$\sigma_{SvM} = sgn(tr(\sigma)) \cdot \sigma_{vM}$$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, ...). Signed von Mises stress for each entry. + Array of shape (...). Signed von Mises stress for each entry. + Tensor rank is reduced by one. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - von_mises = calc_von_mises_stress(stress_voigt) + von_mises = calc_von_mises_stress(stress_vector_voigt) invariant_1 = ( - stress_voigt[:, 0, ...] + stress_voigt[:, 1, ...] + stress_voigt[:, 2, ...] + stress_vector_voigt[..., 0] + + stress_vector_voigt[..., 1] + + stress_vector_voigt[..., 2] ) sign = np.sign(invariant_1).astype(np.float64, copy=False) + sign[np.isclose(invariant_1, 0)] = 1.0 return sign * von_mises # Tresca functions -def calc_tresca_stress(stress_voigt: NDArray[np.float64]) -> NDArray[np.float64]: +def calc_tresca_stress(stress_vector_voigt: NDArray[np.float64]) -> NDArray[np.float64]: r"""Calculate Tresca (maximum shear) stress for each stress state. ??? abstract "Math Equations" $$\sigma_{\tau_{max}} = \frac{\sigma_{1} - \sigma_{3}}{2}$$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, ...). Tresca stress for each entry. For principal + Array of shape (...). Tresca stress for each entry. For principal stresses σ1 ≥ σ2 ≥ σ3, Tresca = (σ1 − σ3)/2. + Tensor rank is reduced by one. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - principals = calc_principal_stresses(stress_voigt) + principals = calc_principal_stresses(stress_vector_voigt) return 0.5 * (principals[..., 0] - principals[..., 2]) def calc_signed_tresca_by_hydrostatic( - stress_voigt: NDArray[np.float64], + stress_vector_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: r"""Calculate signed Tresca stress for each stress state. @@ -398,23 +400,27 @@ def calc_signed_tresca_by_hydrostatic( $$\sigma_{S\tau_{max}} = sgn(\sigma_H) \cdot \sigma_{\tau_{max}}$$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, ...). Signed Tresca stress for each entry. + Array of shape (...). Signed Tresca stress for each entry. + Tensor rank is reduced by one. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - tresca = calc_tresca_stress(stress_voigt) - sign = np.sign(calc_hydrostatic_stress(stress_voigt)) + tresca = calc_tresca_stress(stress_vector_voigt) + hydrostatic_stress = calc_hydrostatic_stress(stress_vector_voigt) + + sign = np.sign(hydrostatic_stress) + sign[np.isclose(hydrostatic_stress, 0)] = 1.0 + return sign * tresca def calc_signed_tresca_by_max_abs_principal( - stress_voigt: NDArray[np.float64], + stress_vector_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: r"""Calculate signed Tresca stress for each stress state. @@ -427,29 +433,21 @@ def calc_signed_tresca_by_max_abs_principal( $$ Args: - stress_voigt: Array of shape (n, 6, ...). The second dimension of size 6 - contains the Voigt stress components. Additional trailing dimensions - are supported. + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. Returns: - Array of shape (n, ...). Signed Tresca stress for each entry. + Array of shape (...). Signed Tresca stress for each entry. + Tensor rank is reduced by one. Raises: - ValueError: If input is not a 2D array with 6 columns. + ValueError: If the last dimension is not of size 6. """ - tresca = calc_tresca_stress(stress_voigt) - principals = calc_principal_stresses(stress_voigt) + tresca = calc_tresca_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) + sign[np.isclose(avg_13, 0)] = 1.0 return sign * tresca - - -# The function using np.linalg.eigh(tensor) is compatible with tensors constructed as: -# tensor[:, 0, 0, ...] = vector[:, 0, ...] -# tensor[:, 1, 1, ...] = vector[:, 1, ...] -# tensor[:, 2, 2, ...] = vector[:, 2, ...] -# tensor[:, [1, 2], [2, 1], ...] = vector[:, [3], ...] -# tensor[:, [0, 2], [2, 0], ...] = vector[:, [4], ...] -# tensor[:, [0, 1], [1, 0], ...] = vector[:, [5], ...] diff --git a/tests/struct_mech/test_stress.py b/tests/struct_mech/test_stress.py index b3bd6f6..f24d394 100644 --- a/tests/struct_mech/test_stress.py +++ b/tests/struct_mech/test_stress.py @@ -1,5 +1,8 @@ """Test functions for stress 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 @@ -16,28 +19,29 @@ from numpy.typing import NDArray from fatpy.struct_mech import stress +from fatpy.utils import voigt @pytest.fixture -def stress_vector_2d() -> NDArray[np.float64]: - """Fixture providing sample stress vectors in Voigt notation. +def stress_vector_sample() -> NDArray[np.float64]: + """Fixture providing sample stress vectors in Voigt notation (3D: shape (2, 3, 6)). Returns: NDArray[np.float64]: Sample stress vectors. """ - # Provide a collection of representative stress states (n,6): - # 1) Uniaxial tension in x (sigma_xx = 100) - # 2) Uniaxial compression in z (sigma_zz = -50) - # 3) Pure hydrostatic (all normal = 30) - # 4) Pure shear (sxy = 40) - # 5) Mixed state + # Two sets, each with three representative stress states (3,6): arr = np.array( [ - [100.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [0.0, 0.0, -50.0, 0.0, 0.0, 0.0], - [30.0, 30.0, 30.0, 0.0, 0.0, 0.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 40.0], - [np.sqrt(2.0), -np.sqrt(2), 0.0, 0.0, 0.0, np.sqrt(2.0)], + [ + [100, 0, 0, 0, 0, 0], # Uniaxial tension in x + [0, 0, -50, 0, 0, 0], # Uniaxial compression in z + [30, 30, 30, 0, 0, 0], # Pure hydrostatic + ], + [ + [0, 0, 0, 0, 0, 40], # Pure shear + [np.sqrt(2), -np.sqrt(2), 0, 0, 0, np.sqrt(2)], # Mixed state + [14, 0, 6, 0, 3, 0], # Another mixed state + ], ], dtype=np.float64, ) @@ -46,7 +50,7 @@ def stress_vector_2d() -> NDArray[np.float64]: @pytest.fixture -def principal_stresses_2d() -> NDArray[np.float64]: +def principal_stresses_sample() -> NDArray[np.float64]: """Fixture providing expected principal stresses for the sample stress vectors. Returns: @@ -54,11 +58,16 @@ def principal_stresses_2d() -> NDArray[np.float64]: """ arr = np.array( [ - [100.0, 0.0, 0.0], - [0.0, 0.0, -50.0], - [30.0, 30.0, 30.0], - [40.0, 0.0, -40.0], - [2.0, 0.0, -2.0], + [ + [100, 0, 0], + [0, 0, -50], + [30, 30, 30], + ], + [ + [40, 0, -40], + [2, 0, -2], + [15, 5, 0], + ], ], dtype=np.float64, ) @@ -66,183 +75,293 @@ def principal_stresses_2d() -> NDArray[np.float64]: return arr -# TODO: implement 3D cases -def stress_vector_3d(stress_vector_2d: NDArray[np.float64]) -> NDArray[np.float64]: - """Fixture providing sample stress vectors in Voigt notation. 3D case. - - This is similar to the 2D case but with additional axis possibly representing time. +@pytest.fixture +def hydrostatic_stress_sample() -> NDArray[np.float64]: + """Fixture providing expected hydrostatic stresses for the sample stress vectors. Returns: - NDArray[np.float64]: Sample stress vectors. + NDArray[np.float64]: Expected hydrostatic stresses. """ - # Expand to 3D by adding an additional axis (e.g., time) - return np.repeat(stress_vector_2d[np.newaxis, :, :], repeats=4, axis=0) + arr = np.array( + [ + [100 / 3, -50 / 3, 30], + [0, 0, 20 / 3], + ], + dtype=np.float64, + ) + return arr -def test_calc_principal_stresses_and_directions_shape_2d( - stress_vector_2d: NDArray[np.float64], + +def test_calc_principal_stresses_and_directions_shape( + stress_vector_sample: NDArray[np.float64], ) -> None: """Test shape of principal stresses and directions output.""" principals, directions = stress.calc_principal_stresses_and_directions( - stress_vector_2d + stress_vector_sample ) - assert principals.shape[0] == stress_vector_2d.shape[0] - assert principals.shape[1] == 3 - - assert directions.shape[0] == stress_vector_2d.shape[0] - assert directions.shape[1] == 3 - assert directions.shape[2] == 3 + assert principals.shape[:-1] == stress_vector_sample.shape[:-1] + assert principals.shape[-1] == 3 -def test_calc_principal_stresses_and_directions_ordering_2d( - stress_vector_2d: NDArray[np.float64], +def test_calc_principal_stresses_and_directions_ordering( + stress_vector_sample: NDArray[np.float64], ) -> None: """Test ordering of principal stresses (descending).""" - principals, dirctions = stress.calc_principal_stresses_and_directions( - stress_vector_2d + principals, directions = stress.calc_principal_stresses_and_directions( + stress_vector_sample ) - assert np.all(principals[:, 0] >= principals[:, 1]) - assert np.all(principals[:, 1] >= principals[:, 2]) + assert np.all(principals[..., 0] >= principals[..., 1]) + assert np.all(principals[..., 1] >= principals[..., 2]) # test if directions order matches principal stresses order # by checking A*v = λ*v for each principal stress/direction pair - for i in range(stress_vector_2d.shape[0]): - sigma = stress_vector_2d[i] - tensor = np.array( - [ - [sigma[0], sigma[5], sigma[4]], - [sigma[5], sigma[1], sigma[3]], - [sigma[4], sigma[3], sigma[2]], - ] - ) - for j in range(3): - eigvec = dirctions[i, :, j] - eigval = principals[i, j] - # Check if A*v = λ*v - Av = tensor @ eigvec - lv = eigval * eigvec + for idx in np.ndindex(principals.shape[:-1]): + stress_tensor = voigt.voigt_to_tensor(stress_vector_sample[idx]) + for i in range(3): + principal_stress = principals[idx + (i,)] + direction = directions[idx + (slice(None), i)] + Av = stress_tensor @ direction + lv = principal_stress * direction assert np.allclose(Av, lv, atol=1e-12) -def test_calc_principal_stresses_and_directions_values_2d( - stress_vector_2d: NDArray[np.float64], - principal_stresses_2d: NDArray[np.float64], +def test_calc_principal_stresses_and_directions( + stress_vector_sample: NDArray[np.float64], + principal_stresses_sample: NDArray[np.float64], ) -> None: - """Test specific known values of principal stresses.""" - principals, _ = stress.calc_principal_stresses_and_directions(stress_vector_2d) - assert np.allclose(principals, principal_stresses_2d, atol=1e-12) + """Test shape of principal stresses and directions output.""" + principals, directions = stress.calc_principal_stresses_and_directions( + stress_vector_sample + ) + assert np.allclose(principals, principal_stresses_sample, atol=1e-12) -def test_calc_principal_stresses_shape_2d( - stress_vector_2d: NDArray[np.float64], + +def test_calc_principal_stresses( + stress_vector_sample: NDArray[np.float64], + principal_stresses_sample: NDArray[np.float64], ) -> None: """Test shape of principal stresses output.""" - principals = stress.calc_principal_stresses(stress_vector_2d) - assert principals.shape[0] == stress_vector_2d.shape[0] - assert principals.shape[1] == 3 + principals = stress.calc_principal_stresses(stress_vector_sample) + + assert principals.shape[:-1] == stress_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_stresses_sample, atol=1e-12) -def test_calc_principal_stresses_ordering_2d( - stress_vector_2d: NDArray[np.float64], + +def test_calc_stress_invariants( + stress_vector_sample: NDArray[np.float64], ) -> None: - """Test ordering of principal stresses (descending).""" - principals = stress.calc_principal_stresses(stress_vector_2d) - assert np.all(principals[:, 0] >= principals[:, 1]) - assert np.all(principals[:, 1] >= principals[:, 2]) + invariants = stress.calc_stress_invariants(stress_vector_sample) + + assert invariants.shape == stress_vector_sample.shape[:-1] + (3,) + for idx in np.ndindex(invariants.shape[:-1]): + stress_tensor = voigt.voigt_to_tensor(stress_vector_sample[idx]) -def test_calc_principal_stresses_values_2d( - stress_vector_2d: NDArray[np.float64], - principal_stresses_2d: NDArray[np.float64], + I1 = np.trace(stress_tensor) + I2 = 0.5 * (I1**2 - np.trace(stress_tensor @ stress_tensor)) + I3 = np.linalg.det(stress_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_hydrostatic_stress( + stress_vector_sample: NDArray[np.float64], + hydrostatic_stress_sample: NDArray[np.float64], ) -> None: - """Test specific known values of principal stresses.""" - principals = stress.calc_principal_stresses(stress_vector_2d) - assert np.allclose(principals, principal_stresses_2d, atol=1e-12) + hydrostatic = stress.calc_hydrostatic_stress(stress_vector_sample) + + assert hydrostatic.shape == stress_vector_sample.shape[:-1] + + assert np.allclose(hydrostatic, hydrostatic_stress_sample, atol=1e-12) -def test_calc_principal_directions_shape_2d( - stress_vector_2d: NDArray[np.float64], +def test_calc_stress_deviator( + stress_vector_sample: NDArray[np.float64], ) -> None: - """Test shape of principal directions output.""" - directions = stress.calc_principal_directions(stress_vector_2d) - assert directions.shape[0] == stress_vector_2d.shape[0] - assert directions.shape[1] == 3 - assert directions.shape[2] == 3 + deviator = stress.calc_stress_deviator(stress_vector_sample) + assert deviator.shape == stress_vector_sample.shape -def test_calc_principal_directions_ordering_2d( - stress_vector_2d: NDArray[np.float64], + for idx in np.ndindex(deviator.shape[:-1]): + stress_tensor = voigt.voigt_to_tensor(stress_vector_sample[idx]) + hydrostatic = np.trace(stress_tensor) / 3.0 + hydrostatic_tensor = np.eye(3) * hydrostatic + deviator_tensor = stress_tensor - hydrostatic_tensor + deviator_voigt = voigt.tensor_to_voigt(deviator_tensor) + assert np.allclose(deviator[idx], deviator_voigt, atol=1e-12) + + +def test_calc_von_mises( + stress_vector_sample: NDArray[np.float64], ) -> None: - """Test ordering of principal directions matches principal stresses order.""" - principals = stress.calc_principal_stresses(stress_vector_2d) - directions = stress.calc_principal_directions(stress_vector_2d) - # For each stress state, check A*v = λ*v for each direction - for i in range(stress_vector_2d.shape[0]): - sigma = stress_vector_2d[i] - tensor = np.array( - [ - [sigma[0], sigma[5], sigma[4]], - [sigma[5], sigma[1], sigma[3]], - [sigma[4], sigma[3], sigma[2]], - ] - ) - for j in range(3): - eigvec = directions[i, :, j] - eigval = principals[i, j] - Av = tensor @ eigvec - lv = eigval * eigvec - assert np.allclose(Av, lv, atol=1e-12) + """Test calculation of von Mises stress. + + Uses the definition based on the stress deviator tensor. + + Args: + stress_vector_sample (NDArray[np.float64]): Sample stress vectors. + """ + von_mises_stress = stress.calc_von_mises_stress(stress_vector_sample) + + assert von_mises_stress.shape == stress_vector_sample.shape[:-1] + + for idx in np.ndindex(von_mises_stress.shape): + stress_tensor = voigt.voigt_to_tensor(stress_vector_sample[idx]) + s = stress_tensor - np.eye(3) * (np.trace(stress_tensor) / 3.0) + vm = np.sqrt(1.5 * np.sum(s**2)) + assert np.isclose(von_mises_stress[idx], vm, atol=1e-12) + +def test_calc_signed_von_mises_by_hydrostatic( + stress_vector_sample: NDArray[np.float64], +) -> None: + """Test calculation of signed von Mises stress. + + Uses the definition based on the hydrostatic stress. + + Args: + stress_vector_sample (NDArray[np.float64]): Sample stress vectors. + """ + signed_von_mises = stress.calc_signed_von_mises_by_hydrostatic(stress_vector_sample) + + assert signed_von_mises.shape == stress_vector_sample.shape[:-1] + + for idx in np.ndindex(signed_von_mises.shape): + stress_tensor = voigt.voigt_to_tensor(stress_vector_sample[idx]) + s = stress_tensor - np.eye(3) * (np.trace(stress_tensor) / 3.0) + vm = np.sqrt(1.5 * np.sum(s**2)) + hydrostatic = np.trace(stress_tensor) / 3.0 + sign = np.sign(hydrostatic) + sign = 1.0 if np.isclose(hydrostatic, 0.0) else sign + assert np.isclose(signed_von_mises[idx], sign * vm, atol=1e-12) -def test_invariants_and_hydrostatic_deviator( - stress_vector_2d: NDArray[np.float64], + +def test_calc_signed_von_mises_by_max_abs_principal( + stress_vector_sample: NDArray[np.float64], ) -> None: - invariants = stress.calc_stress_invariants(stress_vector_2d) - hydro = stress.calc_hydrostatic_stress(stress_vector_2d) - deviator = stress.calc_stress_deviator(stress_vector_2d) - # Shape checks - assert invariants.shape == (stress_vector_2d.shape[0], 3) - assert hydro.shape == (stress_vector_2d.shape[0],) - assert deviator.shape == stress_vector_2d.shape - # For hydrostatic case (index 2), deviator should be zero - assert np.allclose(deviator[2], 0.0, atol=1e-12) - # For hydrostatic case, invariants: I1 = 90, I2 = ? compute from definition - I1 = invariants[2, 0] - assert np.isclose(I1, 90.0) - - -def test_von_mises_and_signed_variants(stress_vector_2d: NDArray[np.float64]) -> None: - vm = stress.calc_von_mises_stress(stress_vector_2d) - svm_h = stress.calc_signed_von_mises_by_hydrostatic(stress_vector_2d) - svm_p = stress.calc_signed_von_mises_by_max_abs_principal(stress_vector_2d) - svm_i1 = stress.calc_signed_von_mises_by_first_invariant(stress_vector_2d) - # Shape - assert vm.shape == (stress_vector_2d.shape[0],) - # Non-negative for von Mises - assert np.all(vm >= 0.0) - # For hydrostatic pure case index 2, von Mises == 0 and signed variants 0 - assert np.isclose(vm[2], 0.0) - assert np.isclose(svm_h[2], 0.0) - assert np.isclose(svm_p[2], 0.0) - assert np.isclose(svm_i1[2], 0.0) - # For pure shear case index 3 (sxy = 40): von Mises = sqrt(3/2)*|shear|*? check - # Compute expected manually from definition used in implementation - sx, sy, sz = stress_vector_2d[3, :3] - sxy = stress_vector_2d[3, 5] - expected_vm3 = np.sqrt( - 0.5 * ((sx - sy) ** 2 + (sy - sz) ** 2 + (sz - sx) ** 2 + 6 * (sxy**2)) + """Test calculation of signed von Mises stress. + + Uses the definition based on the principal stress with the largest absolute value. + + Args: + stress_vector_sample (NDArray[np.float64]): Sample stress vectors. + """ + signed_von_mises = stress.calc_signed_von_mises_by_max_abs_principal( + stress_vector_sample ) - assert np.isclose(vm[3], expected_vm3) - - -def test_tresca_and_signed_variants(stress_vector_2d: NDArray[np.float64]) -> None: - tresca = stress.calc_tresca_stress(stress_vector_2d) - stresca_h = stress.calc_signed_tresca_by_hydrostatic(stress_vector_2d) - stresca_p = stress.calc_signed_tresca_by_max_abs_principal(stress_vector_2d) - # Shape - assert tresca.shape == (stress_vector_2d.shape[0],) - # Hydrostatic case (index 2) tresca == 0 - assert np.isclose(tresca[2], 0.0) - assert np.isclose(stresca_h[2], 0.0) - # For uniaxial tension (index 0), principal stresses should include 100 and two zeros -> tresca = (100 - 0)/2 - assert np.isclose(tresca[0], 50.0) + + assert signed_von_mises.shape == stress_vector_sample.shape[:-1] + + for idx in np.ndindex(signed_von_mises.shape): + stress_tensor = voigt.voigt_to_tensor(stress_vector_sample[idx]) + s = stress_tensor - np.eye(3) * (np.trace(stress_tensor) / 3.0) + vm = np.sqrt(1.5 * np.sum(s**2)) + principals = np.linalg.eigvalsh(stress_tensor) + avg13 = 0.5 * (principals[0] + principals[2]) + sign = np.sign(avg13) + sign = 1.0 if np.isclose(avg13, 0.0) else sign + svm = sign * vm + assert np.isclose(signed_von_mises[idx], svm, atol=1e-12) + + +def test_calc_signed_von_mises_by_first_invariant( + stress_vector_sample: NDArray[np.float64], +) -> None: + """Test calculation of signed von Mises stress. + + Uses the definition based on the first stress invariant. + + Args: + stress_vector_sample (NDArray[np.float64]): Sample stress vectors. + """ + signed_von_mises = stress.calc_signed_von_mises_by_first_invariant( + stress_vector_sample + ) + + assert signed_von_mises.shape == stress_vector_sample.shape[:-1] + + for idx in np.ndindex(signed_von_mises.shape): + stress_tensor = voigt.voigt_to_tensor(stress_vector_sample[idx]) + s = stress_tensor - np.eye(3) * (np.trace(stress_tensor) / 3.0) + vm = np.sqrt(1.5 * np.sum(s**2)) + I1 = np.trace(stress_tensor) + sign = np.sign(I1) + sign = 1.0 if np.isclose(I1, 0.0) else sign + svm = sign * vm + assert np.isclose(signed_von_mises[idx], svm, atol=1e-12) + + +def test_calc_tresca( + stress_vector_sample: NDArray[np.float64], +) -> None: + """Test calculation of Tresca stress. + + Args: + stress_vector_sample (NDArray[np.float64]): Sample stress vectors. + """ + tresca_stress = stress.calc_tresca_stress(stress_vector_sample) + + assert tresca_stress.shape == stress_vector_sample.shape[:-1] + + for idx in np.ndindex(tresca_stress.shape): + stress_tensor = voigt.voigt_to_tensor(stress_vector_sample[idx]) + principals = np.linalg.eigvalsh(stress_tensor) + tresca = (principals[2] - principals[0]) / 2 + assert np.isclose(tresca_stress[idx], tresca, atol=1e-12) + + +def test_calc_signed_tresca_by_hydrostatic( + stress_vector_sample: NDArray[np.float64], +) -> None: + """Test calculation of signed Tresca stress. + + Uses the definition based on the hydrostatic stress. + + Args: + stress_vector_sample (NDArray[np.float64]): Sample stress vectors. + """ + signed_tresca = stress.calc_signed_tresca_by_hydrostatic(stress_vector_sample) + + assert signed_tresca.shape == stress_vector_sample.shape[:-1] + + for idx in np.ndindex(signed_tresca.shape): + stress_tensor = voigt.voigt_to_tensor(stress_vector_sample[idx]) + principals = np.linalg.eigvalsh(stress_tensor) + tresca = (principals[2] - principals[0]) / 2 + hydrostatic = np.trace(stress_tensor) / 3.0 + sign = np.sign(hydrostatic) + sign = 1.0 if np.isclose(hydrostatic, 0.0) else sign + assert np.isclose(signed_tresca[idx], sign * tresca, atol=1e-12) + + +def test_calc_signed_tresca_by_max_abs_principal( + stress_vector_sample: NDArray[np.float64], +) -> None: + """Test calculation of signed Tresca stress. + + Uses the definition based on the principal stress with the largest absolute value. + + Args: + stress_vector_sample (NDArray[np.float64]): Sample stress vectors. + """ + signed_tresca = stress.calc_signed_tresca_by_max_abs_principal(stress_vector_sample) + + assert signed_tresca.shape == stress_vector_sample.shape[:-1] + + for idx in np.ndindex(signed_tresca.shape): + stress_tensor = voigt.voigt_to_tensor(stress_vector_sample[idx]) + principals = np.linalg.eigvalsh(stress_tensor) + tresca = (principals[2] - principals[0]) / 2 + avg13 = 0.5 * (principals[0] + principals[2]) + sign = np.sign(avg13) + sign = 1.0 if np.isclose(avg13, 0.0) else sign + assert np.isclose(signed_tresca[idx], sign * tresca, atol=1e-12) From 8097e6205afef61e0a5542322a79fa2716012166 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Mon, 29 Sep 2025 01:30:02 +0200 Subject: [PATCH 15/17] Fixed docstring spacing --- src/fatpy/struct_mech/stress.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index c7269eb..f8f7ce0 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -46,10 +46,10 @@ 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. + (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. + 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. From 58e1a5f78c551b6b2b3dc9aa397d70c7f70ae56a Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Sat, 4 Oct 2025 20:37:12 +0200 Subject: [PATCH 16/17] Fixed bug in signed variants with array of shape (6,) --- src/fatpy/struct_mech/stress.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index f8f7ce0..ee3eb06 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -295,7 +295,7 @@ def calc_signed_von_mises_by_hydrostatic( hydrostatic_stress = calc_hydrostatic_stress(stress_vector_voigt) sign = np.sign(hydrostatic_stress).astype(np.float64, copy=False) - sign[np.isclose(hydrostatic_stress, 0)] = 1.0 + sign = np.where(np.isclose(hydrostatic_stress, 0), 1.0, sign) return sign * von_mises @@ -326,9 +326,11 @@ def calc_signed_von_mises_by_max_abs_principal( """ von_mises = calc_von_mises_stress(stress_voigt) principals = calc_principal_stresses(stress_voigt) + avg_13 = 0.5 * (principals[..., 0] + principals[..., 2]) sign = np.sign(avg_13).astype(np.float64, copy=False) - sign[np.isclose(avg_13, 0)] = 1.0 + sign = np.where(np.isclose(avg_13, 0), 1.0, sign) + return sign * von_mises @@ -361,7 +363,7 @@ def calc_signed_von_mises_by_first_invariant( ) sign = np.sign(invariant_1).astype(np.float64, copy=False) - sign[np.isclose(invariant_1, 0)] = 1.0 + sign = np.where(np.isclose(invariant_1, 0), 1.0, sign) return sign * von_mises @@ -414,7 +416,7 @@ def calc_signed_tresca_by_hydrostatic( hydrostatic_stress = calc_hydrostatic_stress(stress_vector_voigt) sign = np.sign(hydrostatic_stress) - sign[np.isclose(hydrostatic_stress, 0)] = 1.0 + sign = np.where(np.isclose(hydrostatic_stress, 0), 1.0, sign) return sign * tresca @@ -448,6 +450,6 @@ def calc_signed_tresca_by_max_abs_principal( avg_13 = 0.5 * (principals[..., 0] + principals[..., 2]) sign = np.sign(avg_13).astype(np.float64, copy=False) - sign[np.isclose(avg_13, 0)] = 1.0 + sign = np.where(np.isclose(avg_13, 0), 1.0, sign) return sign * tresca From 0abd1ff30bfd2fd8abb122404b88db4f4b73f642 Mon Sep 17 00:00:00 2001 From: Vybornak2 Date: Sun, 5 Oct 2025 23:08:04 +0200 Subject: [PATCH 17/17] Added atol, rtol, modified singed variants docstrings, removed principal_directions --- src/fatpy/struct_mech/stress.py | 153 ++++++++++++++++++++++---------- 1 file changed, 108 insertions(+), 45 deletions(-) diff --git a/src/fatpy/struct_mech/stress.py b/src/fatpy/struct_mech/stress.py index ee3eb06..5f5f702 100644 --- a/src/fatpy/struct_mech/stress.py +++ b/src/fatpy/struct_mech/stress.py @@ -20,8 +20,6 @@ """ -# TODO: implement atol, rtol to signed stress functions - import numpy as np from numpy.typing import NDArray @@ -97,35 +95,6 @@ def calc_principal_stresses( return np.sort(eigvals, axis=-1)[..., ::-1] -def calc_principal_directions( - stress_vector_voigt: NDArray[np.float64], -) -> NDArray[np.float64]: - r"""Calculate principal directions (eigenvectors) for each stress state. - - ??? abstract "Math Equations" - Principal directions are found by solving the eigenvalue problem - for the stress tensor: - - $$ \sigma \mathbf{v} = \lambda \mathbf{v} $$ - - Args: - stress_vector_voigt: Array of shape (..., 6). The last dimension contains the - Voigt stress components. Leading dimensions are preserved. - - Returns: - Array of shape (..., 3, 3). Principal directions (columns are eigenvectors) - aligned with descending principal stresses: σ1, σ2, σ3. - - Raises: - ValueError: If the last dimension is not of size 6. - """ - voigt.check_shape(stress_vector_voigt) - - # Reuse the ordering logic from the combined function to ensure consistency - _, eigvecs = calc_principal_stresses_and_directions(stress_vector_voigt) - return eigvecs - - def calc_stress_invariants( stress_vector_voigt: NDArray[np.float64], ) -> NDArray[np.float64]: @@ -272,17 +241,37 @@ def calc_von_mises_stress( def calc_signed_von_mises_by_hydrostatic( stress_vector_voigt: NDArray[np.float64], + rtol: float = 1e-5, + atol: float = 1e-8, ) -> NDArray[np.float64]: r"""Calculate signed von Mises stress for each stress state. 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) + + 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" - $$\sigma_{SvM} = sgn(\sigma_H) \cdot \sigma_{vM}$$ + $$ + \sigma_{SvM} = \begin{cases} + +\sigma_{vM} & \text{if } \sigma_H \geq 0 \\ + -\sigma_{vM} & \text{if } \sigma_H < 0 + \end{cases} + $$ Args: stress_vector_voigt: Array of shape (..., 6). The last dimension contains the Voigt stress components. Leading dimensions are preserved. + rtol: Relative tolerance for comparing hydrostatic stress to zero. + Default is 1e-5. + atol: Absolute tolerance for comparing hydrostatic stress to zero. + Default is 1e-8. Returns: Array of shape (...). Signed von Mises stress for each entry. @@ -295,27 +284,44 @@ def calc_signed_von_mises_by_hydrostatic( hydrostatic_stress = calc_hydrostatic_stress(stress_vector_voigt) sign = np.sign(hydrostatic_stress).astype(np.float64, copy=False) - sign = np.where(np.isclose(hydrostatic_stress, 0), 1.0, sign) + sign = np.where(np.isclose(hydrostatic_stress, 0, rtol=rtol, atol=atol), 1.0, sign) return sign * von_mises def calc_signed_von_mises_by_max_abs_principal( stress_voigt: NDArray[np.float64], + rtol: float = 1e-5, + atol: float = 1e-8, ) -> NDArray[np.float64]: r"""Calculate signed von Mises stress for each stress state. Sign is determined by average of the maximum and minimum principal stresses. - In case the maximum absolute principal stress is zero, the sign is set to +1. - In case the maximum absolute principal stress is equal to negative value of - the minimum principal stress, the sign is set to +1 as well. + + 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" - $$\sigma_{SvM} = sgn(\sigma_{max,abs}) \cdot \sigma_{vM}$$ + $$ + \sigma_{SvM} = \begin{cases} + +\sigma_{vM} & \text{if } \frac{\sigma_1 + \sigma_3}{2} \geq 0 \\ + -\sigma_{vM} & \text{if } \frac{\sigma_1 + \sigma_3}{2} < 0 + \end{cases} + $$ Args: stress_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. + atol: Absolute tolerance for comparing the average of maximum and minimum + principal stresses to zero. Default is 1e-8. Returns: Array of shape (...). Signed von Mises stress for each entry. @@ -329,24 +335,44 @@ def calc_signed_von_mises_by_max_abs_principal( 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), 1.0, sign) + sign = np.where(np.isclose(avg_13, 0, rtol=rtol, atol=atol), 1.0, sign) return sign * von_mises def calc_signed_von_mises_by_first_invariant( stress_vector_voigt: NDArray[np.float64], + rtol: float = 1e-5, + atol: float = 1e-8, ) -> NDArray[np.float64]: r"""Calculate signed von Mises stress for each stress state. 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) + + 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" - $$\sigma_{SvM} = sgn(tr(\sigma)) \cdot \sigma_{vM}$$ + $$ + \sigma_{SvM} = \begin{cases} + +\sigma_{vM} & \text{if } tr(\sigma) \geq 0 \\ + -\sigma_{vM} & \text{if } tr(\sigma) < 0 + \end{cases} + $$ Args: 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 first invariant to zero. + Default is 1e-5. + atol: Absolute tolerance for comparing the first invariant to zero. + Default is 1e-8. Returns: Array of shape (...). Signed von Mises stress for each entry. @@ -363,7 +389,7 @@ def calc_signed_von_mises_by_first_invariant( ) sign = np.sign(invariant_1).astype(np.float64, copy=False) - sign = np.where(np.isclose(invariant_1, 0), 1.0, sign) + sign = np.where(np.isclose(invariant_1, 0, rtol=rtol, atol=atol), 1.0, sign) return sign * von_mises @@ -393,17 +419,37 @@ def calc_tresca_stress(stress_vector_voigt: NDArray[np.float64]) -> NDArray[np.f def calc_signed_tresca_by_hydrostatic( stress_vector_voigt: NDArray[np.float64], + rtol: float = 1e-5, + atol: float = 1e-8, ) -> NDArray[np.float64]: r"""Calculate signed Tresca stress for each stress state. 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) + + 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" - $$\sigma_{S\tau_{max}} = sgn(\sigma_H) \cdot \sigma_{\tau_{max}}$$ + $$ + \sigma_{S\tau_{max}} = \begin{cases} + +\sigma_{\tau_{max}} & \text{if } \sigma_H \geq 0 \\ + -\sigma_{\tau_{max}} & \text{if } \sigma_H < 0 + \end{cases} + $$ Args: stress_vector_voigt: Array of shape (..., 6). The last dimension contains the Voigt stress components. Leading dimensions are preserved. + rtol: Relative tolerance for comparing hydrostatic stress to zero. + Default is 1e-5. + atol: Absolute tolerance for comparing hydrostatic stress to zero. + Default is 1e-8. Returns: Array of shape (...). Signed Tresca stress for each entry. @@ -416,27 +462,44 @@ def calc_signed_tresca_by_hydrostatic( hydrostatic_stress = calc_hydrostatic_stress(stress_vector_voigt) sign = np.sign(hydrostatic_stress) - sign = np.where(np.isclose(hydrostatic_stress, 0), 1.0, sign) + sign = np.where(np.isclose(hydrostatic_stress, 0, rtol=rtol, atol=atol), 1.0, sign) return sign * tresca def calc_signed_tresca_by_max_abs_principal( stress_vector_voigt: NDArray[np.float64], + rtol: float = 1e-5, + atol: float = 1e-8, ) -> NDArray[np.float64]: r"""Calculate signed Tresca stress for each stress state. 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) + + 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" $$ - \sigma_{S\tau_{max}} = sgn(\frac{\sigma_{1}+\sigma_{3}}{2}) - \cdot \sigma_{\tau_{max}} + \sigma_{S\tau_{max}} = \begin{cases} + +\sigma_{\tau_{max}} & \text{if } \frac{\sigma_1 + \sigma_3}{2} \geq 0 \\ + -\sigma_{\tau_{max}} & \text{if } \frac{\sigma_1 + \sigma_3}{2} < 0 + \end{cases} $$ Args: 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. + atol: Absolute tolerance for comparing the average of maximum and minimum + principal stresses to zero. Default is 1e-8. Returns: Array of shape (...). Signed Tresca stress for each entry. @@ -450,6 +513,6 @@ def calc_signed_tresca_by_max_abs_principal( 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), 1.0, sign) + sign = np.where(np.isclose(avg_13, 0, rtol=rtol, atol=atol), 1.0, sign) return sign * tresca