Merge pull request #14 from cmccomb/codex/fix-divide-by-zero-warnings

Fix divide-by-zero handling in truss computations

Author: cmccomb

tests/test_member_math.py

@@ -0,0 +1,54 @@
+import warnings
+
+import numpy
+import pytest
+
+import trussme
+
+from tests.helpers import build_triangle_truss
+
+
+def test_member_safety_factors_handle_zero_force_without_runtime_warnings() -> None:
+    truss = build_triangle_truss()
+    member = truss.members[0]
+
+    with warnings.catch_warnings(record=True) as captured:
+        warnings.simplefilter("always")
+        yielding = member.fos_yielding
+        buckling = member.fos_buckling
+
+    runtime_warnings = [
+        warning for warning in captured if issubclass(warning.category, RuntimeWarning)
+    ]
+
+    assert numpy.isinf(yielding)
+    assert numpy.isinf(buckling)
+    assert not runtime_warnings
+
+
+def test_zero_length_member_direction_returns_zero_vector() -> None:
+    truss = trussme.Truss()
+    truss.add_free_joint([0.0, 0.0, 0.0])
+    truss.add_free_joint([0.0, 0.0, 0.0])
+    truss.add_member(0, 1)
+
+    with warnings.catch_warnings(record=True) as captured:
+        warnings.simplefilter("always")
+        direction = truss.members[0].direction
+
+    runtime_warnings = [
+        warning for warning in captured if issubclass(warning.category, RuntimeWarning)
+    ]
+
+    assert numpy.allclose(direction, numpy.zeros(3))
+    assert not runtime_warnings
+
+
+def test_zero_length_member_stiffness_raises_value_error() -> None:
+    truss = trussme.Truss()
+    truss.add_free_joint([0.0, 0.0, 0.0])
+    truss.add_free_joint([0.0, 0.0, 0.0])
+    truss.add_member(0, 1)
+
+    with pytest.raises(ValueError, match="length must be greater than zero"):
+        _ = truss.members[0].stiffness

tests/test_optimization.py

@@ -276,6 +276,46 @@ def test_make_truss_generator_function_scaled_updates_other_shapes(
         assert generated.members[0].shape._params[key] == pytest.approx(value)
 
 
+@pytest.mark.parametrize(
+    ("shape_factory", "expected"),
+    [
+        (
+            lambda: trussme.Pipe(r=0.0, t=0.0),
+            {"r": 0.04, "t": 0.0},
+        ),
+        (
+            lambda: trussme.Box(w=0.0, h=0.0, t=0.0),
+            {"w": 0.04, "h": 0.0, "t": 0.0},
+        ),
+        (
+            lambda: trussme.Square(w=0.0, h=0.0),
+            {"w": 0.04, "h": 0.0},
+        ),
+    ],
+)
+def test_make_truss_generator_function_scaled_handles_zero_reference_dimensions(
+    shape_factory,
+    expected,
+) -> None:
+    truss = build_triangle_truss(shape=shape_factory)
+    generator = trussme.make_truss_generator_function(
+        truss,
+        joint_optimization=None,
+        member_optimization="scaled",
+    )
+    x0 = optimize.make_x0(
+        truss,
+        joint_optimization=None,
+        member_optimization="scaled",
+    )
+    x0[0] = 0.04
+
+    generated = generator(x0)
+
+    for key, value in expected.items():
+        assert generated.members[0].shape._params[key] == pytest.approx(value)
+
+
 def test_make_truss_generator_function_full_updates_shapes() -> None:
     truss = build_optimization_truss()
     generator = trussme.make_truss_generator_function(

tests/test_visualize.py

@@ -1,3 +1,5 @@
+import warnings
+
 import matplotlib
 
 matplotlib.use("Agg")
@@ -38,6 +40,22 @@ def test_plot_truss_supports_force_coloring() -> None:
     assert figure.axes
 
 
+def test_plot_truss_supports_force_coloring_with_zero_forces() -> None:
+    truss = build_triangle_truss()
+
+    with warnings.catch_warnings(record=True) as captured:
+        warnings.simplefilter("always")
+        figure = plot_truss(truss, starting_shape="force", deflected_shape="force")
+
+    runtime_warnings = [
+        warning for warning in captured if issubclass(warning.category, RuntimeWarning)
+    ]
+
+    assert isinstance(figure, Figure)
+    assert figure.axes
+    assert not runtime_warnings
+
+
 def test_plot_truss_allows_hidden_shapes() -> None:
     truss = build_triangle_truss()
     truss.analyze()

trussme/components.py

@@ -516,12 +516,18 @@ def direction(self) -> NDArray[numpy.float64]:
         vector_length: NDArray[numpy.float64] = numpy.array(
             self.end_joint.coordinates
         ) - numpy.array(self.begin_joint.coordinates)
-        return vector_length / numpy.linalg.norm(vector_length)
+        magnitude = float(numpy.linalg.norm(vector_length))
+        if magnitude == 0.0:
+            return numpy.zeros_like(vector_length, dtype=numpy.float64)
+        return cast(NDArray[numpy.float64], vector_length / magnitude)
 
     @property
     def stiffness(self) -> float:
         """float: The axial stiffness of the member"""
-        return self.elastic_modulus * self.area / self.length
+        length = self.length
+        if length == 0.0:
+            raise ValueError("Member length must be greater than zero")
+        return self.elastic_modulus * self.area / length
 
     @property
     def stiffness_vector(self) -> NDArray[numpy.float64]:
@@ -554,19 +560,28 @@ def force(self, new_force: float):
     @property
     def fos_yielding(self) -> float:
         """float: The factor of safety against yielding"""
-        return self.yield_strength / abs(self.force / self.area)
+        applied_force = abs(self.force)
+        if applied_force == 0.0:
+            return numpy.inf
+        return self.yield_strength * abs(self.area) / applied_force
 
     @property
     def fos_buckling(self) -> float:
         """float: The factor of safety against buckling"""
+        compressive_force = -self.force
+        if compressive_force <= 0.0:
+            return numpy.inf
+
+        length = self.length
+        if length == 0.0:
+            return numpy.inf
+
         fos = (
-            -(
-                (numpy.pi**2)
-                * self.elastic_modulus
-                * self.moment_of_inertia
-                / (self.length**2)
-            )
-            / self.force
+            (numpy.pi**2)
+            * self.elastic_modulus
+            * self.moment_of_inertia
+            / (length**2)
+            / compressive_force
         )
 
         return fos if fos > 0 else numpy.inf

trussme/optimize.py

@@ -5,6 +5,12 @@
 from trussme import Truss, Goals, read_json, Pipe, Box, Square, Bar
 
 
+def _safe_ratio(numerator: float, denominator: float) -> float:
+    if denominator == 0.0:
+        return 0.0
+    return numerator / denominator
+
+
 def make_x0(
     truss: Truss,
     joint_optimization: Optional[Literal["full"]] = "full",
@@ -223,21 +229,27 @@ def truss_generator(x: list[float]) -> Truss:
                 shape_name: str = configured_truss.members[i].shape.name()
                 p = configured_truss.members[i].shape._params
                 if shape_name == "pipe":
+                    thickness_ratio = _safe_ratio(p["t"], p["r"])
                     configured_truss.members[i].shape = Pipe(
-                        r=x[idx], t=x[idx] * p["t"] / p["r"]
+                        r=x[idx], t=x[idx] * thickness_ratio
                     )
                     idx += 1
                 elif shape_name == "box":
+                    height_ratio = _safe_ratio(p["h"], p["w"])
+                    thickness_ratio = _safe_ratio(p["t"], p["w"])
                     configured_truss.members[i].shape = Box(
-                        w=x[idx], h=x[idx] * p["h"] / p["w"], t=x[idx] * p["t"] / p["w"]
+                        w=x[idx],
+                        h=x[idx] * height_ratio,
+                        t=x[idx] * thickness_ratio,
                     )
                     idx += 1
                 elif shape_name == "bar":
                     configured_truss.members[i].shape = Bar(r=x[idx])
                     idx += 1
                 elif shape_name == "square":
+                    height_ratio = _safe_ratio(p["h"], p["w"])
                     configured_truss.members[i].shape = Square(
-                        w=x[idx], h=x[idx] * p["h"] / p["w"]
+                        w=x[idx], h=x[idx] * height_ratio
                     )
                     idx += 1
 

trussme/visualize.py

@@ -52,7 +52,9 @@ def plot_truss(
     ax.axis("equal")
     ax.set_axis_off()
 
-    scaler: float = numpy.max(numpy.abs([member.force for member in truss.members]))
+    scaler = max((abs(member.force) for member in truss.members), default=0.0)
+    if scaler == 0.0:
+        scaler = 1.0
 
     force_colormap = matplotlib.colors.LinearSegmentedColormap.from_list(
         "force",