Add an optional strain exponent for the SVK-material (#905)

* Add an optional strain exponent for the SVK-material

* Make the third orthotropic-direction optional

in `saint_venant_kirchhoff_orthotropic()`

src/felupe/constitution/tensortrax/models/hyperelastic/__init__.py

@@ -64,7 +64,7 @@
 neo_hooke.kwargs = dict(mu=0)
 ogden.kwargs = dict(mu=[0, 0], alpha=[2, -2])
 ogden_roxburgh.kwargs = dict(r=100, m=1, beta=0, material=neo_hooke, mu=0)
-saint_venant_kirchhoff.kwargs = dict(mu=0.0, lmbda=0.0)
+saint_venant_kirchhoff.kwargs = dict(mu=0.0, lmbda=0.0, k=2)
 saint_venant_kirchhoff_orthotropic.kwargs = dict(
     mu=[0.0, 0.0, 0.0],
     lmbda=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0],

src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff.py

@@ -16,10 +16,13 @@
 along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
 """
 
+from tensortrax.math import base, log
+from tensortrax.math import sum as tsum
 from tensortrax.math import trace
+from tensortrax.math.linalg import eigvalsh
 
 
-def saint_venant_kirchhoff(C, mu, lmbda):
+def saint_venant_kirchhoff(C, mu, lmbda, k=2):
     r"""Strain energy function of the isotropic hyperelastic
     `Saint-Venant Kirchhoff <https://en.wikipedia.org/wiki/Hyperelastic_material#Saint_Venant-Kirchhoff_model>`_
     material formulation.
@@ -32,6 +35,9 @@ def saint_venant_kirchhoff(C, mu, lmbda):
         Second Lamé constant (shear modulus).
     lmbda : float
         First Lamé constant (shear modulus).
+    k : float, optional
+        Strain exponent (default is 2). If 2, the Green-Lagrange strain measure is used.
+        For any other value, the family of Seth-Hill strains is used.
 
     Notes
     -----
@@ -79,6 +85,22 @@ def saint_venant_kirchhoff(C, mu, lmbda):
 
 
     """
-    I1 = trace(C) / 2 - 3 / 2
-    I2 = trace(C @ C) / 4 - trace(C) / 2 + 3 / 4
+    eye = base.eye
+
+    if k == 2:
+        E = (C - eye(C)) / 2
+
+        I1 = trace(E)
+        I2 = trace(E @ E)
+
+    else:
+        λ2 = eigvalsh(C)
+        if k == 0:
+            Ek = log(λ2) / 2
+        else:
+            Ek = (λ2 ** (k / 2) - 1) / k
+
+        I1 = tsum(Ek)
+        I2 = tsum(Ek**2)
+
     return mu * I2 + lmbda * I1**2 / 2

src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff_orthotropic.py

@@ -16,11 +16,15 @@
 along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
 """
 from numpy import array
-from tensortrax.math import base, einsum
+from tensortrax.math import base, einsum, log
+from tensortrax.math import sum as tsum
+from tensortrax.math.linalg import eigh
 from tensortrax.math.special import from_triu_1d
 
+from .....math import cross
 
-def saint_venant_kirchhoff_orthotropic(C, mu, lmbda, r1, r2, r3):
+
+def saint_venant_kirchhoff_orthotropic(C, mu, lmbda, r1, r2, r3=None, k=2):
     r"""Strain energy function of the orthotropic hyperelastic
     `Saint-Venant Kirchhoff <https://en.wikipedia.org/wiki/Hyperelastic_material#Saint_Venant-Kirchhoff_model>`_
     material formulation.
@@ -30,16 +34,20 @@ def saint_venant_kirchhoff_orthotropic(C, mu, lmbda, r1, r2, r3):
     C : tensortrax.Tensor
         Right Cauchy-Green deformation tensor.
     mu : list of float
-        List of the three second Lamé parameters :math:`\mu_1,\mu_2, \mu3`.
+        List of the three second Lamé parameters :math:`\mu_1,\mu_2, \mu_3`.
     lmbda : list of float
-        List of six (upper triangle) first Lamé parameters
-        :math:`\lambda_11, \lambda_12, \lambda_13, \lambda_22, \lambda_23, \lambda_33`.
+        List of six (upper triangle) first Lamé parameters :math:`\lambda_{11},
+        \lambda_{12}, \lambda_{13}, \lambda_{22}, \lambda_{23}, \lambda_{33}`.
     r1 : list of float
         First normal vector of planes of symmetry.
     r2 : list of float
         Second normal vector of planes of symmetry.
-    r3 : list of float
-        Third normal vector of planes of symmetry.
+    r3 : list of float or None, optional
+        Third normal vector of planes of symmetry. If None, the third normal vector
+        is evaluated as :math:`r_1 \times r_2`. Default is None.
+    k : float, optional
+        Strain exponent (default is 2). If 2, the Green-Lagrange strain measure is used.
+        For any other value, the family of Seth-Hill strains is used.
 
     Notes
     -----
@@ -96,12 +104,24 @@ def saint_venant_kirchhoff_orthotropic(C, mu, lmbda, r1, r2, r3):
 
     """
     eye = base.eye
+
+    if k == 2:
+        E = (C - eye(C)) / 2
+
+    else:
+        λ2, M = eigh(C)
+        if k == 0:
+            E = tsum(log(λ2) / 2 * M, axis=0)
+        else:
+            E = tsum((λ2 ** (k / 2) - 1) / k * M, axis=0)
+
     μ = array(mu)
     λ = from_triu_1d(array(lmbda))
 
-    r = array([r1, r2, r3])
+    if r3 is None:
+        r3 = cross(r1, r2)
 
-    E = (C - eye(C)) / 2
+    r = array([r1, r2, r3])
     Err = einsum("ai...,ij...,aj...->a...", r, E, r)
 
     λI1 = einsum("ab,a...,b...->...", λ / 2, Err, Err)

tests/test_constitution.py

@@ -342,6 +342,8 @@ def neo_hooke(C, mu=1):
     for model, kwargs, incompressible in [
         (neo_hooke, {"mu": 1}, True),
         (fem.saint_venant_kirchhoff, {"mu": 1, "lmbda": 20.0}, False),
+        (fem.saint_venant_kirchhoff, {"mu": 1, "lmbda": 20.0, "k": 0}, False),
+        (fem.saint_venant_kirchhoff, {"mu": 1, "lmbda": 20.0, "k": 1}, False),
         (
             fem.saint_venant_kirchhoff_orthotropic,
             {
@@ -353,6 +355,30 @@ def neo_hooke(C, mu=1):
             },
             False,
         ),
+        (
+            fem.saint_venant_kirchhoff_orthotropic,
+            {
+                "mu": [1, 1, 1],
+                "lmbda": [20, 20, 20, 20, 20, 20],
+                "r1": np.eye(3)[:, 0],
+                "r2": np.eye(3)[:, 1],
+                "r3": None,
+                "k": 0,
+            },
+            False,
+        ),
+        (
+            fem.saint_venant_kirchhoff_orthotropic,
+            {
+                "mu": [1, 1, 1],
+                "lmbda": [20, 20, 20, 20, 20, 20],
+                "r1": np.eye(3)[:, 0],
+                "r2": np.eye(3)[:, 1],
+                "r3": np.eye(3)[:, 2],
+                "k": 1,
+            },
+            False,
+        ),
         (fem.neo_hooke, {"mu": 1}, True),
         (fem.mooney_rivlin, {"C10": 0.3, "C01": 0.8}, True),
         (fem.yeoh, {"C10": 0.5, "C20": -0.1, "C30": 0.02}, True),