Minor/fixes

TODO.md

@@ -2,12 +2,15 @@
 
 ## Major
 
+- [ ] Layer & Slab using pipelines from Mary-Kate + Attributes (value, calculated [bool], pipeline, uncertainty)
+- [ ] Uncertainties propagation
 - [ ] Use Classes for Boundary Types
 - [ ] Automatically figure out type of system
 - [ ] Automatically set boundary conditions based on system
 
 ## Minor
 
+- [ ] Swap to Pytest from Unittest
 - [ ] resolve fracture criterion also when lower than strength criterion
 - [ ] Florian CriterionEvaluator: clarify and fix damping behavior (find_minimum_force / evaluate_coupled_criterion)
   - Expected behavior

src/weac/analysis/criteria_evaluator.py

@@ -721,6 +721,12 @@ def evaluate_SteadyState(
         -----------
         system: SystemModel
             The system model.
+        mode: TouchdownMode, optional
+            Touchdown evaluation mode. The three supported modes are:
+            ``"C_in_contact"`` for a cut distance of ``2 * l_BC``,
+            ``"B_point_contact"`` for a cut distance just shorter than ``l_BC``,
+            and ``"A_free_hanging"`` for a cut distance just shorter than ``l_AB``.
+            Defaults to ``"C_in_contact"``.
         vertical: bool, optional
             Whether to evaluate the system in a vertical configuration.
             Defaults to False.
@@ -1301,15 +1307,30 @@ def _calculate_maximal_stresses(
         )
         max_principal_stress_norm = np.max(principal_stress_norm)
         max_Sxx_norm = np.max(Sxx_norm)
-        # Count height levels as failure-prone when tensile stress exceeds the layer
-        # strength or the slab density is below the configured weak-snow threshold.
-        # zmesh rho is t/mm^3, layer rho is kg/m^3
+        # zmesh rows are ordered from slab top to bottom. Low-density levels only
+        # fail through downward growth from tensile failures above; when they do,
+        # they are excluded from the slab tensile criterion denominator.
+        # zmesh rho is t/mm^3, layer rho is kg/m^3.
         zmesh = analyzer.get_zmesh(dz=1)
         rho_kg_m3 = zmesh["rho"] * 1e12
         tensile_exceeds = np.max(Sxx_norm, axis=1) > 1
         low_density = rho_kg_m3 <= self.criteria_config.low_density_threshold_kg_m3
-        height_level_prone_to_fail = tensile_exceeds | low_density
-        slab_tensile_criterion = np.mean(height_level_prone_to_fail)
+        height_level_prone_to_fail = np.zeros_like(tensile_exceeds, dtype=bool)
+        all_above_prone_to_fail = True
+        for index, is_low_density in enumerate(low_density):
+            if is_low_density:
+                height_level_prone_to_fail[index] = all_above_prone_to_fail
+            else:
+                height_level_prone_to_fail[index] = tensile_exceeds[index]
+            all_above_prone_to_fail &= height_level_prone_to_fail[index]
+
+        low_density_prone_to_fail = low_density & height_level_prone_to_fail
+        load_bearing_levels = ~low_density_prone_to_fail
+        slab_tensile_criterion = (
+            np.mean(height_level_prone_to_fail[load_bearing_levels])
+            if np.any(load_bearing_levels)
+            else 1.0
+        )
         if print_call_stats:
             analyzer.print_call_stats(
                 message="_calculate_maximal_stresses Call Statistics"

src/weac/components/criteria_config.py

@@ -50,7 +50,8 @@ class CriteriaConfig(BaseModel):
         Order of magnitude for stress envelope. Default is 1.0.
     low_density_threshold_kg_m3 : float
         Slab density threshold in kg/m^3 below which a layer is treated as weak snow
-        and counted as prone to tensile failure in the slab tensile criterion.
+        and excluded from the slab tensile criterion percentage when broken through
+        directional growth from above.
     """
 
     fn: float = Field(
@@ -94,6 +95,6 @@ class CriteriaConfig(BaseModel):
         gt=0,
         description=(
             "Slab density threshold in kg/m^3 below which a layer is treated as weak "
-            "snow in the slab tensile criterion"
+            "snow and conditionally excluded from the slab tensile criterion percentage"
         ),
     )

src/weac/core/field_quantities.py

@@ -107,13 +107,14 @@ def sig(self, Z: np.ndarray, unit: StressUnit = "MPa") -> float | np.ndarray:
         return -self._unit_factor(unit) * self.es.weak_layer.kn * self.w(Z)
 
     def tau(self, Z: np.ndarray, unit: StressUnit = "MPa") -> float | np.ndarray:
-        """Weak-layer shear stress `tau = -kt * (w' * h/2 - u(h=H/2))`"""
+        """Weak-layer shear stress `tau = kt * h * (w' / 2 - u(h=H/2) / h)`"""
         return (
-            -self._unit_factor(unit)
+            self._unit_factor(unit)
             * self.es.weak_layer.kt
+            * self.es.weak_layer.h
             * (
-                self.dw_dx(Z) * self.es.weak_layer.h / 2
-                - self.u(Z, h0=self.es.slab.H / 2)
+                self.dw_dx(Z) / 2
+                - self.u(Z, h0=self.es.slab.H / 2) / self.es.weak_layer.h
             )
         )
 
@@ -122,7 +123,7 @@ def eps(self, Z: np.ndarray) -> float | np.ndarray:
         return -self.w(Z) / self.es.weak_layer.h
 
     def gamma(self, Z: np.ndarray) -> float | np.ndarray:
-        """Weak-layer shear strain `gamma = (w' * h/2 - u(h=H/2)) / h`"""
+        """Weak-layer shear strain `gamma = w' / 2 - u(h=H/2) / h`"""
         return (
             self.dw_dx(Z) / 2 - self.u(Z, h0=self.es.slab.H / 2) / self.es.weak_layer.h
         )

tests/analysis/test_analyzer.py

@@ -156,6 +156,37 @@ def test_energy_release_rates_shapes(self):
         self.assertEqual(Gdif.shape, (4,))
         self.assertTrue(np.isfinite(Gdif).all())
 
+    def test_energy_release_rate_integrands_non_negative(self):
+        """Test that ERR integrands are non-negative for matching stress/strain."""
+        slope_angle = 20.0
+        system = SystemModel(
+            model_input=ModelInput(
+                scenario_config=ScenarioConfig(phi=slope_angle, system_type="skier"),
+                layers=[Layer()],
+                weak_layer=WeakLayer(),
+                segments=[Segment(), Segment()],
+            ),
+            config=Config(),
+        )
+        analyzer = Analyzer(system_model=system, printing_enabled=False)
+
+        z_uncracked = np.array([[0.0], [0.0], [1.0], [0.2], [0.0], [0.0]])
+
+        def constant_solution(x):
+            return np.repeat(z_uncracked, np.size(np.atleast_1d(x)), axis=1)
+
+        mode_i = analyzer._integrand_GI(  # pylint: disable=protected-access
+            np.array([0.0, 1.0]), constant_solution, constant_solution
+        )
+        mode_ii = analyzer._integrand_GII(  # pylint: disable=protected-access
+            np.array([0.0, 1.0]), constant_solution, constant_solution
+        )
+
+        self.assertTrue(np.all(mode_i >= 0), "Mode I integrand should be non-negative")
+        self.assertTrue(
+            np.all(mode_ii >= 0), "Mode II integrand should be non-negative"
+        )
+
     def test_internal_and_external_potentials_pst(self):
         """Test internal and external potentials for PST."""
         # Ensure PST-specific methods run

tests/analysis/test_criteria_evaluator.py

@@ -68,19 +68,19 @@ def test_stress_envelope_adam_unpublished(self):
         self.assertGreater(result[0], 0)
 
     @patch("weac.analysis.criteria_evaluator.Analyzer")
-    def test_calculate_maximal_stresses_uses_configured_low_density_threshold(
+    def test_calculate_maximal_stresses_applies_directional_low_density_exclusion(
         self, mock_analyzer_cls
-    ):
-        """Test that the slab tensile criterion uses the configured density cutoff."""
-        sxx_kpa = np.zeros((3, 1))
-        principal_stress_kPa = np.zeros((3, 1))
-        sxx_norm = np.full((3, 1), 0.5)
-        principal_stress_norm = np.full((3, 1), 0.5)
+    ):  # pylint: disable=protected-access
+        """Test that weak snow is excluded only after downward failure growth."""
+        sxx_kpa = np.zeros((4, 1))
+        principal_stress_kPa = np.zeros((4, 1))
+        sxx_norm = np.array([[1.5], [0.5], [0.5], [0.5]])
+        principal_stress_norm = np.full((4, 1), 0.5)
 
         mock_analyzer = mock_analyzer_cls.return_value
         mock_analyzer.rasterize_solution.return_value = (
             None,
-            np.array([0, 1, 2]),
+            np.array([0, 1, 2, 3]),
             None,
         )
         mock_analyzer.Sxx.side_effect = (
@@ -92,20 +92,22 @@ def test_calculate_maximal_stresses_uses_configured_low_density_threshold(
             )
         )
         mock_analyzer.get_zmesh.return_value = {
-            "rho": np.array([90.0, 110.0, 130.0]) * 1e-12
+            "rho": np.array([130.0, 90.0, 90.0, 130.0]) * 1e-12
         }
         system = SimpleNamespace(scenario=SimpleNamespace(phi=30.0))
 
         # Access the helper directly so the test can isolate the density-threshold logic.
-        default_result = CriteriaEvaluator(
+        top_broken_result = CriteriaEvaluator(
             CriteriaConfig()
         )._calculate_maximal_stresses(system=system)  # pylint: disable=protected-access
-        tuned_result = CriteriaEvaluator(
-            CriteriaConfig(low_density_threshold_kg_m3=120)
+
+        sxx_norm = np.array([[0.5], [0.5], [0.5], [1.5]])
+        top_unbroken_result = CriteriaEvaluator(
+            CriteriaConfig()
         )._calculate_maximal_stresses(system=system)  # pylint: disable=protected-access
 
-        self.assertAlmostEqual(default_result.slab_tensile_criterion, 1 / 3)
-        self.assertAlmostEqual(tuned_result.slab_tensile_criterion, 2 / 3)
+        self.assertAlmostEqual(top_broken_result.slab_tensile_criterion, 1 / 2)
+        self.assertAlmostEqual(top_unbroken_result.slab_tensile_criterion, 1 / 4)
 
     def test_find_minimum_force_convergence(self):
         """Test the convergence of find_minimum_force."""

tests/analysis/test_slab_tensile_comparisons.py

@@ -80,8 +80,8 @@ def _setup_from_cm(*layers: tuple[float, float]) -> SetupDefinition:
     ),
     ComparisonCase(
         name="case_2",
-        setup_a=_setup_from_cm((50, 75), (20, 225)),
-        setup_b=_setup_from_cm((30, 75), (20, 225)),
+        setup_a=_setup_from_cm((30, 75), (20, 225)),
+        setup_b=_setup_from_cm((50, 75), (20, 225)),
     ),
     ComparisonCase(
         name="case_3",
@@ -115,8 +115,13 @@ def _setup_from_cm(*layers: tuple[float, float]) -> SetupDefinition:
     ),
     ComparisonCase(
         name="case_9",
-        setup_a=_setup_from_cm((15, 275), (40, 75)),
-        setup_b=_setup_from_cm((40, 75), (15, 275)),
+        setup_a=_setup_from_cm((40, 75), (15, 275)),
+        setup_b=_setup_from_cm((15, 275), (40, 75)),
+    ),
+    ComparisonCase(
+        name="case_10",
+        setup_a=_setup_from_cm((30, 75), (20, 275)),
+        setup_b=_setup_from_cm((50, 75), (20, 275)),
     ),
 )
 

tests/core/test_field_quantities.py

@@ -272,7 +272,7 @@ def test_weak_layer_shear_stress(self):
         H = self.fq.es.slab.H
         u_surface = self.fq.u(self.Z, h0=H / 2)
 
-        expected = -self.fq.es.weak_layer.kt * (self.Z[3, :] * h / 2 - u_surface)
+        expected = self.fq.es.weak_layer.kt * (self.Z[3, :] * h / 2 - u_surface)
         np.testing.assert_array_almost_equal(
             tau,
             expected,