From 3ab32a5fb9b22ac6c660d5cff0af12255f8784d5 Mon Sep 17 00:00:00 2001
From: Andreas Dutzler <a.dutzler@gmail.com>
Date: Fri, 25 Oct 2024 13:30:44 +0200
Subject: [PATCH] Fix wrong result of `ConstitutiveMaterial.plot()` if any
 stretch <= 0

this raises an error now to avoid the wrong results
---
 CHANGELOG.md                     |  3 +++
 src/felupe/constitution/_view.py | 22 ++++++++++++++++++++++
 tests/test_constitution.py       | 16 ++++++++++++++++
 3 files changed, 41 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 35ffbfb7..ab3e74e7 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/src/felupe/constitution/_view.py b/src/felupe/constitution/_view.py
index 194a1744..a4070df4 100644
--- a/src/felupe/constitution/_view.py
+++ b/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)
diff --git a/tests/test_constitution.py b/tests/test_constitution.py
index 3d209f2b..660dab08 100644
--- a/tests/test_constitution.py
+++ b/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()