Tensile Strength Evaluation as part of Steady State Analysis (#34)

This pull request refactors and improves the stress and steady-state
evaluation logic in the `weac.analysis` package. It introduces new
result data classes and fixes plotting issues related to stress
calculations and visualization.

### Refactoring and Data Model Improvements
* Introduced `MaximalStressResult` and `SteadyStateResult` data classes,
replacing the old `SSERRResult` for clearer and more comprehensive
result encapsulation. The new classes include normalized stress fields
and system state at touchdown.
(`[[1]](diffhunk://#diff-cd31ff029b0e336203154640e70d582b1554ca04f41c94048eaec4dee7daab6dL98-R129)`,
`[[2]](diffhunk://#diff-cd31ff029b0e336203154640e70d582b1554ca04f41c94048eaec4dee7daab6dL110-R153)`,
`[[3]](diffhunk://#diff-1b897e0700621752c15b1115da95c4b2a2d61aa57dde50d6b780bf490269b16cL11-R11)`,
`[[4]](diffhunk://#diff-1b897e0700621752c15b1115da95c4b2a2d61aa57dde50d6b780bf490269b16cL21-R21)`)
* Updated method names and references from `SSERR` to `SteadyState` for
consistency and clarity, including renaming the evaluation method and
result types.
(`[[1]](diffhunk://#diff-cd31ff029b0e336203154640e70d582b1554ca04f41c94048eaec4dee7daab6dL644-R684)`,
`[[2]](diffhunk://#diff-cd31ff029b0e336203154640e70d582b1554ca04f41c94048eaec4dee7daab6dL691-R734)`,
`[[3]](diffhunk://#diff-1b897e0700621752c15b1115da95c4b2a2d61aa57dde50d6b780bf490269b16cL11-R11)`,
`[[4]](diffhunk://#diff-1b897e0700621752c15b1115da95c4b2a2d61aa57dde50d6b780bf490269b16cL21-R21)`)

### Stress Calculation Enhancements
* Refactored the `Sxx` method in `analyzer.py` to support normalization
of axial normal stress to tensile strength, improving numerical
stability and result interpretability.
(`[[1]](diffhunk://#diff-8468e03357e32957d46ec295646154a4a64e750fa249655a80c6620d9ad32440L228-R228)`,
`[[2]](diffhunk://#diff-8468e03357e32957d46ec295646154a4a64e750fa249655a80c6620d9ad32440R242-R243)`,
`[[3]](diffhunk://#diff-8468e03357e32957d46ec295646154a4a64e750fa249655a80c6620d9ad32440L261-R269)`,
`[[4]](diffhunk://#diff-8468e03357e32957d46ec295646154a4a64e750fa249655a80c6620d9ad32440L277-R291)`)
* Improved normalization logic in `principal_stress_slab` to correctly
convert tensile strength units and normalize principal stresses.
(`[[1]](diffhunk://#diff-8468e03357e32957d46ec295646154a4a64e750fa249655a80c6620d9ad32440R451-R452)`,
`[[2]](diffhunk://#diff-8468e03357e32957d46ec295646154a4a64e750fa249655a80c6620d9ad32440L463-R478)`)
* Added `_calculate_maximal_stresses` helper method to centralize
maximal stress computation for steady-state evaluation.
(`[src/weac/analysis/criteria_evaluator.pyR1211-R1241](diffhunk://#diff-cd31ff029b0e336203154640e70d582b1554ca04f41c94048eaec4dee7daab6dR1211-R1241)`)

### Plotting and Visualization Fixes
* Updated `plot_deformed` in `plotter.py` to robustly handle NaN values
in weak-layer coordinates, ensuring correct plotting and preventing
errors when visualizing stress/displacement fields.
(`[[1]](diffhunk://#diff-737cdaa926ddb7c8e3803d8568adddbe44745a23baca2d6bd39fc27e8a41a1dcL888-R934)`,
`[[2]](diffhunk://#diff-737cdaa926ddb7c8e3803d8568adddbe44745a23baca2d6bd39fc27e8a41a1dcL971-R978)`)

These changes collectively enhance the clarity of stress analysis and
steady-state evaluation in the codebase.

<!-- This is an auto-generated comment: release notes by coderabbit.ai
-->
## Summary by CodeRabbit

* **New Features**
* Per-layer stress normalization option and a new deformation
visualization routine with selectable fields and normalization.

* **Bug Fixes**
* Robust handling and masking of NaN/invalid weak-layer coordinates in
deformation plotting.

* **Refactor**
* Steady‑state evaluation payload reorganized and renamed: includes
energy-release metric, detailed maximal-stress metrics, and system info.

* **Tests**
* Added tests for stress unit conversion and normalization consistency.

* **Chores**
  * Gitignore reinstates .weac-reference and adds a dev/ entry.

<sub>✏️ Tip: You can customize this high-level summary in your review
settings.</sub>
<!-- end of auto-generated comment: release notes by coderabbit.ai -->

Author: pillowbeast

.gitignore

@@ -56,4 +56,7 @@ scratch/
 temp*
 old*
 
-.weac-reference/
+.weac-reference/
+
+# Folder for development and local testing
+dev/

demo/demo.ipynb

@@ -8,7 +8,7 @@
     CoupledCriterionResult,
     CriteriaEvaluator,
     FindMinimumForceResult,
-    SSERRResult,
+    SteadyStateResult,
 )
 from .plotter import Plotter
 
@@ -18,6 +18,6 @@
     "CoupledCriterionHistory",
     "CoupledCriterionResult",
     "FindMinimumForceResult",
-    "SSERRResult",
+    "SteadyStateResult",
     "Plotter",
 ]

src/weac/analysis/__init__.py

@@ -214,7 +214,7 @@ def get_zmesh(self, dz=2):
             ],  # Convert to t/mm^3
             "tensile_strength": [
                 layer.tensile_strength for layer in self.sm.slab.layers
-            ],
+            ],  # in kPa
         }
 
         # Repeat properties for each grid point in the layer
@@ -225,7 +225,7 @@ def get_zmesh(self, dz=2):
         return si
 
     @track_analyzer_call
-    def Sxx(self, Z, phi, dz=2, unit="kPa"):
+    def Sxx(self, Z, phi, dz=2, unit="kPa", normalize: bool = False):
         """
         Compute axial normal stress in slab layers.
 
@@ -239,6 +239,10 @@ def Sxx(self, Z, phi, dz=2, unit="kPa"):
             Element size along z-axis (mm). Default is 2 mm.
         unit : {'kPa', 'MPa'}, optional
             Desired output unit. Default is 'kPa'.
+        normalize : bool, optional
+            Toggle normalization. If True, normalize stress values to the tensile strength of each layer (dimensionless).
+            When normalized, the `unit` parameter is ignored and values are returned as ratios.
+            Default is False.
 
         Returns
         -------
@@ -258,13 +262,13 @@ def Sxx(self, Z, phi, dz=2, unit="kPa"):
         m = Z.shape[1]
 
         # Initialize axial normal stress Sxx
-        Sxx = np.zeros(shape=[n, m])
+        Sxx_MPa = np.zeros(shape=[n, m])
 
         # Compute axial normal stress Sxx at grid points in MPa
         for i, z in enumerate(zi):
-            E = zmesh["E"][i]
+            E_MPa = zmesh["E"][i]
             nu = zmesh["nu"][i]
-            Sxx[i, :] = E / (1 - nu**2) * self.sm.fq.du_dx(Z, z)
+            Sxx_MPa[i, :] = E_MPa / (1 - nu**2) * self.sm.fq.du_dx(Z, z)
 
         # Calculate weight load at grid points and superimpose on stress field
         qt = -rho * G_MM_S2 * np.sin(np.deg2rad(phi))
@@ -274,14 +278,22 @@ def Sxx(self, Z, phi, dz=2, unit="kPa"):
         # Sxx[-1, :] += qt[-1] * (zi[-1] - zi[-2])
         # New Implementation: Changed for numerical stability
         dz = np.diff(zi)
-        Sxx[:-1, :] += qt[:-1, np.newaxis] * dz[:, np.newaxis]
-        Sxx[-1, :] += qt[-1] * dz[-1]
+        Sxx_MPa[:-1, :] += qt[:-1, np.newaxis] * dz[:, np.newaxis]
+        Sxx_MPa[-1, :] += qt[-1] * dz[-1]
+
+        # Normalize tensile stresses to tensile strength
+        if normalize:
+            tensile_strength_kPa = zmesh["tensile_strength"]
+            tensile_strength_MPa = tensile_strength_kPa / 1e3
+            # Normalize axial normal stress to layers' tensile strength
+            normalized_Sxx = Sxx_MPa / tensile_strength_MPa[:, None]
+            return normalized_Sxx
 
         # Return axial normal stress in specified unit
-        return convert[unit] * Sxx
+        return convert[unit] * Sxx_MPa
 
     @track_analyzer_call
-    def Txz(self, Z, phi, dz=2, unit="kPa"):
+    def Txz(self, Z, phi, dz=2, unit="kPa", normalize: bool = False):
         """
         Compute shear stress in slab layers.
 
@@ -295,6 +307,9 @@ def Txz(self, Z, phi, dz=2, unit="kPa"):
             Element size along z-axis (mm). Default is 2 mm.
         unit : {'kPa', 'MPa'}, optional
             Desired output unit. Default is 'kPa'.
+        normalize : bool, optional
+            Toggle normalization.  If True, normalize shear stress values to the tensile strength of each layer (dimensionless).
+            When normalized, the `unit` parameter is ignored and values are returned as ratios. Default is False.
 
         Returns
         -------
@@ -332,14 +347,22 @@ def Txz(self, Z, phi, dz=2, unit="kPa"):
 
         # Integrate -dsxx_dx along z and add cumulative weight load
         # to obtain shear stress Txz in MPa
-        Txz = cumulative_trapezoid(dsxx_dx, zi, axis=0, initial=0)
-        Txz += cumulative_trapezoid(qt, zi, initial=0)[:, None]
+        Txz_MPa = cumulative_trapezoid(dsxx_dx, zi, axis=0, initial=0)
+        Txz_MPa += cumulative_trapezoid(qt, zi, initial=0)[:, None]
+
+        # Normalize shear stresses to tensile strength
+        if normalize:
+            tensile_strength_kPa = zmesh["tensile_strength"]
+            tensile_strength_MPa = tensile_strength_kPa / 1e3
+            # Normalize shear stress to layers' tensile strength
+            normalized_Txz = Txz_MPa / tensile_strength_MPa[:, None]
+            return normalized_Txz
 
         # Return shear stress Txz in specified unit
-        return convert[unit] * Txz
+        return convert[unit] * Txz_MPa
 
     @track_analyzer_call
-    def Szz(self, Z, phi, dz=2, unit="kPa"):
+    def Szz(self, Z, phi, dz=2, unit="kPa", normalize: bool = False):
         """
         Compute transverse normal stress in slab layers.
 
@@ -353,6 +376,10 @@ def Szz(self, Z, phi, dz=2, unit="kPa"):
             Element size along z-axis (mm). Default is 2 mm.
         unit : {'kPa', 'MPa'}, optional
             Desired output unit. Default is 'kPa'.
+        normalize : bool, optional
+            Toggle normalization. If True, normalize stress values to the tensile strength of each layer (dimensionless).
+            When normalized, the `unit` parameter is ignored and values are returned as ratios.
+            Default is False.
 
         Returns
         -------
@@ -392,11 +419,19 @@ def Szz(self, Z, phi, dz=2, unit="kPa"):
         # Integrate dsxx_dxdx twice along z to obtain transverse
         # normal stress Szz in MPa
         integrand = cumulative_trapezoid(dsxx_dxdx, zi, axis=0, initial=0)
-        Szz = cumulative_trapezoid(integrand, zi, axis=0, initial=0)
-        Szz += cumulative_trapezoid(-qn, zi, initial=0)[:, None]
+        Szz_MPa = cumulative_trapezoid(integrand, zi, axis=0, initial=0)
+        Szz_MPa += cumulative_trapezoid(-qn, zi, initial=0)[:, None]
 
-        # Return shear stress txz in specified unit
-        return convert[unit] * Szz
+        # Normalize tensile stresses to tensile strength
+        if normalize:
+            tensile_strength_kPa = zmesh["tensile_strength"]
+            tensile_strength_MPa = tensile_strength_kPa / 1e3
+            # Normalize transverse normal stress to layers' tensile strength
+            normalized_Szz = Szz_MPa / tensile_strength_MPa[:, None]
+            return normalized_Szz
+
+        # Return transverse normal stress Szz in specified unit
+        return convert[unit] * Szz_MPa
 
     @track_analyzer_call
     def principal_stress_slab(
@@ -438,6 +473,8 @@ def principal_stress_slab(
             'min', or if normalization of compressive principal stress
             is requested.
         """
+        convert = {"kPa": 1e3, "MPa": 1}
+
         # Raise error if specified component is not available
         if val not in ["min", "max"]:
             raise ValueError(f"Component {val} not defined.")
@@ -460,9 +497,10 @@ def principal_stress_slab(
         # Normalize tensile stresses to tensile strength
         if normalize and val == "max":
             zmesh = self.get_zmesh(dz=dz)
-            tensile_strength = zmesh["tensile_strength"]
+            tensile_strength_kPa = zmesh["tensile_strength"]
+            tensile_strength_converted = tensile_strength_kPa / 1e3 * convert[unit]
             # Normalize maximum principal stress to layers' tensile strength
-            normalized_Ps = Ps / tensile_strength[:, None]
+            normalized_Ps = Ps / tensile_strength_converted[:, None]
             return normalized_Ps
 
         # Return absolute principal stresses

src/weac/analysis/analyzer.py

@@ -95,9 +95,38 @@ class CoupledCriterionResult:
 
 
 @dataclass
-class SSERRResult:
+class MaximalStressResult:
     """
-    Holds the results of the SSERR evaluation.
+    Holds the results of the maximal stress evaluation.
+
+    Attributes:
+    -----------
+    principal_stress_kPa: np.ndarray
+        The principal stress in kPa.
+    Sxx_kPa: np.ndarray
+        The axial normal stress in kPa.
+    principal_stress_norm: np.ndarray
+        The normalized principal stress to the tensile strength of the layers.
+    Sxx_norm: np.ndarray
+        The normalized axial normal stress to the tensile strength of the layers.
+    max_principal_stress_norm: float
+        The normalized maximum principal stress to the tensile strength of the layers.
+    max_Sxx_norm: float
+        The normalized maximum axial normal stress to the tensile strength of the layers.
+    """
+
+    principal_stress_kPa: np.ndarray
+    Sxx_kPa: np.ndarray
+    principal_stress_norm: np.ndarray
+    Sxx_norm: np.ndarray
+    max_principal_stress_norm: float
+    max_Sxx_norm: float
+
+
+@dataclass
+class SteadyStateResult:
+    """
+    Holds the results of the Steady State evaluation.
 
     Attributes:
     -----------
@@ -107,15 +136,21 @@ class SSERRResult:
         The message of the evaluation.
     touchdown_distance : float
         The touchdown distance.
-    SSERR : float
-        The Steady-State Energy Release Rate calculated with the
-        touchdown distance from G_I and G_II.
+    energy_release_rate : float
+        The steady-state energy release rate calculated with the
+        touchdown distance from the differential energy release rate.
+    maximal_stress_result: MaximalStressResult
+        The maximal stresses in the system at the touchdown distance.
+    system: SystemModel
+        The modified system model used for the steady state evaluation.
     """
 
     converged: bool
     message: str
     touchdown_distance: float
-    SSERR: float
+    energy_release_rate: float
+    maximal_stress_result: MaximalStressResult
+    system: SystemModel
 
 
 @dataclass
@@ -641,12 +676,12 @@ def evaluate_coupled_criterion(
             _recursion_depth=_recursion_depth + 1,
         )
 
-    def evaluate_SSERR(
+    def evaluate_SteadyState(
         self,
         system: SystemModel,
         vertical: bool = False,
         print_call_stats: bool = False,
-    ) -> SSERRResult:
+    ) -> SteadyStateResult:
         """
         Evaluates the Touchdown Distance in the Steady State and the Steady State
         Energy Release Rate.
@@ -688,12 +723,19 @@ def evaluate_SSERR(
         system_copy.update_scenario(segments=segments, scenario_config=scenario_config)
         touchdown_distance = system_copy.slab_touchdown.touchdown_distance
         analyzer = Analyzer(system_copy, printing_enabled=print_call_stats)
-        G, _, _ = analyzer.differential_ERR(unit="J/m^2")
-        return SSERRResult(
+        energy_release_rate, _, _ = analyzer.differential_ERR(unit="J/m^2")
+        maximal_stress_result = self._calculate_maximal_stresses(
+            system_copy, print_call_stats=print_call_stats
+        )
+        if print_call_stats:
+            analyzer.print_call_stats(message="evaluate_SteadyState Call Statistics")
+        return SteadyStateResult(
             converged=True,
-            message="SSERR evaluation successful.",
+            message="Steady State evaluation successful.",
             touchdown_distance=touchdown_distance,
-            SSERR=G,
+            energy_release_rate=energy_release_rate,
+            maximal_stress_result=maximal_stress_result,
+            system=system_copy,
         )
 
     def find_minimum_force(
@@ -1170,3 +1212,51 @@ def _fracture_toughness_exceedance(
 
         # Return the difference from the target
         return g_delta_diff - target
+
+    def _calculate_maximal_stresses(
+        self,
+        system: SystemModel,
+        print_call_stats: bool = False,
+    ) -> MaximalStressResult:
+        """
+        Calculate the maximal stresses in the system.
+
+        Parameters
+        ----------
+        system : SystemModel
+            The system model to analyze.
+        print_call_stats : bool, optional
+            Whether to print analyzer call statistics. Default is False.
+
+        Returns
+        -------
+        MaximalStressResult
+            Object containing both absolute (in kPa) and normalized stress fields,
+            along with maximum normalized stress values.
+        """
+        analyzer = Analyzer(system, printing_enabled=print_call_stats)
+        _, Z, _ = analyzer.rasterize_solution(num=4000, mode="cracked")
+        Sxx_kPa = analyzer.Sxx(Z=Z, phi=system.scenario.phi, dz=5, unit="kPa")
+        principal_stress_kPa = analyzer.principal_stress_slab(
+            Z=Z, phi=system.scenario.phi, dz=5, unit="kPa"
+        )
+        Sxx_norm = analyzer.Sxx(
+            Z=Z, phi=system.scenario.phi, dz=5, unit="kPa", normalize=True
+        )
+        principal_stress_norm = analyzer.principal_stress_slab(
+            Z=Z, phi=system.scenario.phi, dz=5, unit="kPa", normalize=True
+        )
+        max_principal_stress_norm = np.max(principal_stress_norm)
+        max_Sxx_norm = np.max(Sxx_norm)
+        if print_call_stats:
+            analyzer.print_call_stats(
+                message="_calculate_maximal_stresses Call Statistics"
+            )
+        return MaximalStressResult(
+            principal_stress_kPa=principal_stress_kPa,
+            Sxx_kPa=Sxx_kPa,
+            principal_stress_norm=principal_stress_norm,
+            Sxx_norm=Sxx_norm,
+            max_principal_stress_norm=max_principal_stress_norm,
+            max_Sxx_norm=max_Sxx_norm,
+        )