Fix wrong result of ConstitutiveMaterial.plot() if any stretch <= 0

CHANGELOG.md

@@ -13,6 +13,9 @@ All notable changes to this project will be documented in this file. The format
 - Change default `np.einsum(..., order="K")` to `np.einsum(..., order="C")` in the methods of `Field`, `FieldAxisymmetric`, `FieldPlaneStrain` and `FieldContainer`.
 - Fix the number of points for non-disconnected dual meshes. This reduces the assembled (sparse) vector- and matrix-shapes, which are defined on mixed-fields.
 
+### Fixed
+- Fix wrong results of `ConstitutiveMaterial.plot()` if any stretch is non-physical, i.e. lower or equal zero. This raises an error now.
+
 ## [9.0.0] - 2024-09-06
 
 ### Added

src/felupe/constitution/_view.py

@@ -24,6 +24,15 @@
 from ..math import det, linsteps
 
 
+def check_stretches(stretches):
+    "Check if any stretch is lower or equal zero."
+
+    if np.any(stretches <= 0.0):
+        raise ValueError("All stretches must greater than 0.")
+
+    return
+
+
 class PlotMaterial:
     "Plot force-stretch curves of constitutive material formulations."
 
@@ -203,7 +212,10 @@ def uniaxial(self, stretches=None):
         if stretches is None:
             stretches = self.ux
 
+        check_stretches(stretches)
+
         λ1 = stretches
+
         λ2 = λ3 = 1 / np.sqrt(λ1)
         eye = np.eye(3).reshape(3, 3, 1, 1)
 
@@ -265,6 +277,8 @@ def planar(self, stretches=None):
         if stretches is None:
             stretches = self.ps
 
+        check_stretches(stretches)
+
         λ1 = stretches
         λ2 = np.ones_like(λ1)
         λ3 = 1 / λ1
@@ -327,6 +341,8 @@ def biaxial(self, stretches=None):
         if stretches is None:
             stretches = self.bx
 
+        check_stretches(stretches)
+
         λ1 = λ2 = stretches
         λ3 = 1 / λ1**2
         eye = np.eye(3).reshape(3, 3, 1, 1)
@@ -475,6 +491,8 @@ def uniaxial(self, stretches=None):
         if stretches is None:
             stretches = self.ux
 
+        check_stretches(stretches)
+
         λ1 = stretches
         λ2 = λ3 = 1 / np.sqrt(λ1)
         eye = np.eye(3).reshape(3, 3, 1, 1)
@@ -512,6 +530,8 @@ def planar(self, stretches=None):
         if stretches is None:
             stretches = self.ps
 
+        check_stretches(stretches)
+
         λ1 = stretches
         λ2 = np.ones_like(λ1)
         λ3 = 1 / λ1
@@ -554,6 +574,8 @@ def biaxial(self, stretches=None):
         if stretches is None:
             stretches = self.bx
 
+        check_stretches(stretches)
+
         λ1 = λ2 = stretches
         λ3 = 1 / λ1**2
         eye = np.eye(3).reshape(3, 3, 1, 1)

tests/test_constitution.py

@@ -707,6 +707,21 @@ def test_laplace():
     assert A[0].shape == (3, 3, 3, 3, 1, 1)
 
 
+def test_plot_negative_stretches():
+    stretches = np.linspace(-0.5, 1, 16)
+    umat = fem.NeoHooke(mu=1.0, bulk=2.0)
+
+    for incompressible in [False, True]:
+        with pytest.raises(ValueError):
+            umat.plot(ux=stretches, incompressible=incompressible)
+
+        with pytest.raises(ValueError):
+            umat.plot(bx=stretches, incompressible=incompressible)
+
+        with pytest.raises(ValueError):
+            umat.plot(ps=stretches, incompressible=incompressible)
+
+
 if __name__ == "__main__":
     test_nh()
     test_linear()
@@ -728,3 +743,4 @@ def test_laplace():
     test_lagrange()
     test_lagrange_statevars()
     test_laplace()
+    test_plot_negative_stretches()