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 128c5b1..535c7a6 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.13.1 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.18.2 hooks: - id: mypy 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/pyproject.toml b/pyproject.toml index c4d05ee..beecc56 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,18 +129,23 @@ 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 warn_unused_configs = true show_error_context = true show_column_numbers = true +disallow_untyped_decorators = false [build-system] requires = ["setuptools>=61.0", "setuptools_scm[toml]>=6.2", "wheel"] diff --git a/src/fatpy/py.typed b/src/fatpy/py.typed new file mode 100644 index 0000000..e69de29 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..5f5f702 --- /dev/null +++ b/src/fatpy/struct_mech/stress.py @@ -0,0 +1,518 @@ +"""Calculate fundamental stress metrics and invariants. + +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: +- 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). + +- Principal directions (eigenvectors) are aligned to this ordering + (columns correspond to σ_1, σ_2, σ_3). + +""" + +import numpy as np +from numpy.typing import NDArray + +from fatpy.utils import voigt + + +def calc_principal_stresses_and_directions( + stress_vector_voigt: NDArray[np.float64], +) -> 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_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 (..., 3). Principal stresses + (descending: σ_1 ≥ σ_2 ≥ σ_3) with leading dimensions preserved. + - eigvecs: Array of shape (..., 3, 3). Principal directions (columns are + eigenvectors) aligned with eigvals in the same order. The last two + dimensions are the 3x3 eigenvector matrix for each input. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(stress_vector_voigt) + + 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) + eigvecs_sorted = np.take_along_axis( + eigvecs, np.expand_dims(sorted_indices, axis=-2), axis=-1 + ) + + return eigvals_sorted, eigvecs_sorted + + +def calc_principal_stresses( + stress_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + r"""Calculate principal stresses for each stress state. + + ??? abstract "Math Equations" + Principal stresses 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). Principal stresses (descending: σ1 ≥ σ2 ≥ σ3). + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(stress_vector_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_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" + $$ + \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_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. + + Returns: + Array of shape (..., 3). The last dimension contains (I1, I2, I3) for + each entry. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(stress_vector_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) + ) + invariant_3 = np.linalg.det(tensor) + + return np.stack((invariant_1, invariant_2, invariant_3), axis=-1) + + +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" + $$ + \sigma_H = \frac{1}{3} tr(\sigma) = + \frac{1}{3} (\sigma_{11} + \sigma_{22} + \sigma_{33}) + $$ + + Args: + stress_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. + + Returns: + Array of shape (...). Hydrostatic stress for each input state. + Tensor rank is reduced by one. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(stress_vector_voigt) + + # Voigt normal components are at last axis indices 0,1,2 + return ( + stress_vector_voigt[..., 0] + + stress_vector_voigt[..., 1] + + stress_vector_voigt[..., 2] + ) / 3.0 + + +def calc_stress_deviator( + stress_vector_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_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. + + Returns: + Array of shape (..., 6). Stress deviator for each entry. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(stress_vector_voigt) + hydrostatic = calc_hydrostatic_stress(stress_vector_voigt) + + 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_vector_voigt: NDArray[np.float64], +) -> NDArray[np.float64]: + 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_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. + + Returns: + Array of shape (...). Von Mises equivalent stress for each entry. + Tensor rank is reduced by one. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + voigt.check_shape(stress_vector_voigt) + + 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 + * ( + (sx - sy) ** 2 + + (sy - sz) ** 2 + + (sz - sx) ** 2 + + 6 * (sxy**2 + syz**2 + sxz**2) + ) + ) + + +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} = \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. + Tensor rank is reduced by one. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + 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.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. + + 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} = \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. + Tensor rank is reduced by one. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + von_mises = calc_von_mises_stress(stress_voigt) + principals = calc_principal_stresses(stress_voigt) + + avg_13 = 0.5 * (principals[..., 0] + principals[..., 2]) + sign = np.sign(avg_13).astype(np.float64, copy=False) + sign = np.where(np.isclose(avg_13, 0, rtol=rtol, atol=atol), 1.0, sign) + + return sign * von_mises + + +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} = \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. + Tensor rank is reduced by one. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + von_mises = calc_von_mises_stress(stress_vector_voigt) + invariant_1 = ( + stress_vector_voigt[..., 0] + + stress_vector_voigt[..., 1] + + stress_vector_voigt[..., 2] + ) + + sign = np.sign(invariant_1).astype(np.float64, copy=False) + sign = np.where(np.isclose(invariant_1, 0, rtol=rtol, atol=atol), 1.0, sign) + + return sign * von_mises + + +# Tresca functions +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_vector_voigt: Array of shape (..., 6). The last dimension contains the + Voigt stress components. Leading dimensions are preserved. + + Returns: + 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 the last dimension is not of size 6. + """ + principals = calc_principal_stresses(stress_vector_voigt) + return 0.5 * (principals[..., 0] - principals[..., 2]) + + +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}} = \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. + Tensor rank is reduced by one. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + tresca = calc_tresca_stress(stress_vector_voigt) + hydrostatic_stress = calc_hydrostatic_stress(stress_vector_voigt) + + sign = np.sign(hydrostatic_stress) + 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}} = \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. + Tensor rank is reduced by one. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + 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.where(np.isclose(avg_13, 0, rtol=rtol, atol=atol), 1.0, sign) + + return sign * tresca 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.""" diff --git a/src/fatpy/utils/voigt.py b/src/fatpy/utils/voigt.py new file mode 100644 index 0000000..bcd269c --- /dev/null +++ b/src/fatpy/utils/voigt.py @@ -0,0 +1,108 @@ +"""Tools for working with Voigt notation in continuum mechanics. + +Overview: + 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. + +""" + +import numpy as np +from numpy.typing import NDArray + +VOIGT_COMPONENTS_COUNT = 6 + + +def check_shape(vector: NDArray[np.float64]) -> None: + """Validate the Voigt vector shape. + + Args: + vector: Array with shape (..., 6) where the last dimension has length 6. + + Raises: + ValueError: If the last dimension is not of size 6. + """ + if vector.shape[-1] != VOIGT_COMPONENTS_COUNT: + raise ValueError( + "Last dimension must correspond to 6 Voigt " + "stress/strain components (..., 6)." + ) + + +def voigt_to_tensor(vector: NDArray[np.float64]) -> NDArray[np.float64]: + """Convert Voigt vectors to symmetric 3x3 tensors. + + Args: + 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 (..., 3, 3). The last two axes are the symmetric + 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) + + # Normal components + 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 + + +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/struct_mech/test_stress.py b/tests/struct_mech/test_stress.py new file mode 100644 index 0000000..f24d394 --- /dev/null +++ b/tests/struct_mech/test_stress.py @@ -0,0 +1,367 @@ +"""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 + 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 +from fatpy.utils import voigt + + +@pytest.fixture +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. + """ + # Two sets, each with three representative stress states (3,6): + arr = np.array( + [ + [ + [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, + ) + + return arr + + +@pytest.fixture +def principal_stresses_sample() -> 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, -50], + [30, 30, 30], + ], + [ + [40, 0, -40], + [2, 0, -2], + [15, 5, 0], + ], + ], + dtype=np.float64, + ) + + return arr + + +@pytest.fixture +def hydrostatic_stress_sample() -> NDArray[np.float64]: + """Fixture providing expected hydrostatic stresses for the sample stress vectors. + + Returns: + NDArray[np.float64]: Expected hydrostatic stresses. + """ + 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( + 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_sample + ) + assert principals.shape[:-1] == stress_vector_sample.shape[:-1] + assert principals.shape[-1] == 3 + + +def test_calc_principal_stresses_and_directions_ordering( + stress_vector_sample: NDArray[np.float64], +) -> None: + """Test ordering of principal stresses (descending).""" + 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]) + + # test if directions order matches principal stresses order + # by checking A*v = λ*v for each principal stress/direction pair + 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( + stress_vector_sample: NDArray[np.float64], + principal_stresses_sample: NDArray[np.float64], +) -> None: + """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( + 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_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_stress_invariants( + stress_vector_sample: NDArray[np.float64], +) -> None: + 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]) + + 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: + 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_stress_deviator( + stress_vector_sample: NDArray[np.float64], +) -> None: + deviator = stress.calc_stress_deviator(stress_vector_sample) + + assert deviator.shape == stress_vector_sample.shape + + 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 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_calc_signed_von_mises_by_max_abs_principal( + stress_vector_sample: NDArray[np.float64], +) -> None: + """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 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) diff --git a/tests/utils/test_voigt.py b/tests/utils/test_voigt.py new file mode 100644 index 0000000..0ee5c6b --- /dev/null +++ b/tests/utils/test_voigt.py @@ -0,0 +1,193 @@ +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) + + +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)