Add more hyperelastic models (#903)

* Add more hyperelastic models

SVK-ortho (tensortrax), extended-tube (JAX), neo-hooke (JAX)

* Aktualisieren von CHANGELOG.md

* Aktualisieren von test_constitution_jax.py

CHANGELOG.md

@@ -12,12 +12,13 @@ All notable changes to this project will be documented in this file. The format
 - Add `math.inplane(A, vectors)` to return the in-plane components of a symmetric tensor `A`, where the plane is defined by its standard unit vectors.
 - Add `constitution.jax.Hyperelastic` as a feature-equivalent alternative to `Hyperelastic` with `jax` as backend.
 - Add `constitution.jax.Material(..., jacobian=None)` with JAX as backend. A custom jacobian-callable may be passed to switch between forward- and backward-mode automatic differentiation.
-- Add material models for JAX-based materials. Hyperelastic models available at `constitution.jax.models.hyperelastic`: `miehe_goektepe_lulei()`, `mooney_rivlin()`, `yeoh()`, `third_order_deformation()`, `van_der_waals()`. Lagrange (stress-based) models available at `constitution.jax.models.lagrange`: `morph()`, `morph_representative_directions()`.
+- Add material models for JAX-based materials. Hyperelastic models available at `constitution.jax.models.hyperelastic`: `extended_tube()`, `miehe_goektepe_lulei()`, `mooney_rivlin()`, `neo_hooke()`, `yeoh()`, `third_order_deformation()`, `van_der_waals()`. Lagrange (stress-based) models available at `constitution.jax.models.lagrange`: `morph()`, `morph_representative_directions()`.
 - Add `constitution.jax.total_lagrange()`, `constitution.jax.updated_lagrange()` and `constitution.jax.isochoric_volumetric_split()` function decorators for the JAX hyperelastic material class.
 - Add an optional keyword-argument `Region.astype(copy=True)` to modify the data types of the arrays of the region in-place if `copy=False`.
 - Add `FieldContainer.evaluate.right_cauchy_green_deformation()` and `math.right_cauchy_green_deformation(field)` to evaluate the right Cauchy-Green deformation tensor.
 - Add `math.strain(field, C=None, ..., **kwargs)` to use a given right Cauchy-Green deformation tensor for the evaluation of the strain tensor.
 - Add the hyperelastic foam models `storakers()` and `blatz_ko()` for both AD-backends.
+- Add the `saint_venant_kirchhoff_orthotropic()` hyperelastic model (tensortrax).
 
 ### Changed
 - Change default `np.einsum(..., order="K")` to `np.einsum(..., order="C")` in the methods of `Field`, `FieldAxisymmetric`, `FieldPlaneStrain` and `FieldContainer`.

docs/felupe/constitution/autodiff/jax.rst

@@ -48,8 +48,10 @@ These material model formulations are defined by a strain energy density functio
 .. autosummary::
 
    blatz_ko
+   extended_tube
    miehe_goektepe_lulei
    mooney_rivlin
+   neo_hooke
    storakers
    third_order_deformation
    van_der_waals
@@ -94,10 +96,14 @@ formulations in :class:`~felupe.constitution.jax.Material`.
 
 .. autofunction:: felupe.constitution.jax.models.hyperelastic.blatz_ko
 
+.. autofunction:: felupe.constitution.jax.models.hyperelastic.extended_tube
+
 .. autofunction:: felupe.constitution.jax.models.hyperelastic.miehe_goektepe_lulei
 
 .. autofunction:: felupe.constitution.jax.models.hyperelastic.mooney_rivlin
 
+.. autofunction:: felupe.constitution.jax.models.hyperelastic.neo_hooke
+
 .. autofunction:: felupe.constitution.jax.models.hyperelastic.storakers
 
 .. autofunction:: felupe.constitution.jax.models.hyperelastic.third_order_deformation

docs/felupe/constitution/autodiff/tensortrax.rst

@@ -42,6 +42,7 @@ These material model formulations are defined by a strain energy density functio
    ogden
    ogden_roxburgh
    saint_venant_kirchhoff
+   saint_venant_kirchhoff_orthotropic
    storakers
    third_order_deformation
    van_der_waals
@@ -117,6 +118,8 @@ and :func:`~felupe.constitution.tensortrax.updated_lagrange` material formulatio
 
 .. autofunction:: felupe.constitution.tensortrax.models.hyperelastic.saint_venant_kirchhoff
 
+.. autofunction:: felupe.constitution.tensortrax.models.hyperelastic.saint_venant_kirchhoff_orthotropic
+
 .. autofunction:: felupe.constitution.tensortrax.models.hyperelastic.storakers
 
 .. autofunction:: felupe.constitution.tensortrax.models.hyperelastic.third_order_deformation

src/felupe/__init__.py

@@ -263,6 +263,7 @@
         ogden,
         ogden_roxburgh,
         saint_venant_kirchhoff,
+        saint_venant_kirchhoff_orthotropic,
         storakers,
         third_order_deformation,
         total_lagrange,
@@ -292,6 +293,7 @@
         "ogden",
         "ogden_roxburgh",
         "saint_venant_kirchhoff",
+        "saint_venant_kirchhoff_orthotropic",
         "storakers",
         "third_order_deformation",
         "van_der_waals",

src/felupe/constitution/__init__.py

@@ -68,6 +68,7 @@
         ogden,
         ogden_roxburgh,
         saint_venant_kirchhoff,
+        saint_venant_kirchhoff_orthotropic,
         storakers,
         third_order_deformation,
         van_der_waals,
@@ -92,6 +93,7 @@
         "ogden",
         "ogden_roxburgh",
         "saint_venant_kirchhoff",
+        "saint_venant_kirchhoff_orthotropic",
         "storakers",
         "third_order_deformation",
         "van_der_waals",

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

@@ -1,25 +1,31 @@
 from ._blatz_ko import blatz_ko
+from ._extended_tube import extended_tube
 from ._miehe_goektepe_lulei import miehe_goektepe_lulei
 from ._mooney_rivlin import mooney_rivlin
+from ._neo_hooke import neo_hooke
 from ._storakers import storakers
 from ._third_order_deformation import third_order_deformation
 from ._van_der_waals import van_der_waals
 from ._yeoh import yeoh
 
 __all__ = [
     "blatz_ko",
+    "extended_tube",
     "miehe_goektepe_lulei",
     "mooney_rivlin",
+    "neo_hooke",
     "storakers",
     "third_order_deformation",
     "van_der_waals",
     "yeoh",
 ]
 
 # default (stable) material parameters
-blatz_ko.kwargs = dict(mu=1)
+blatz_ko.kwargs = dict(mu=0)
+extended_tube.kwargs = dict(Gc=0, Ge=0, beta=1, delta=0)
 miehe_goektepe_lulei.kwargs = dict(mu=0, N=100, U=0, p=2, q=2)
 mooney_rivlin.kwargs = dict(C10=0, C01=0)
+neo_hooke.kwargs = dict(mu=0)
 storakers.kwargs = dict(mu=[0], alpha=[2], beta=[1])
 third_order_deformation.kwargs = dict(C10=0, C01=0, C11=0, C20=0, C30=0)
 van_der_waals.kwargs = dict(mu=0, beta=0, a=0, limit=100)

src/felupe/constitution/jax/models/hyperelastic/_blatz_ko.py

@@ -25,7 +25,6 @@
 
 @wraps(blatz_ko_docstring)
 def blatz_ko(C, mu):
-
     I1 = trace(C)
     I2 = (I1**2 - trace(C @ C)) / 2
     I3 = det(C)

src/felupe/constitution/jax/models/hyperelastic/_extended_tube.py

@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+"""
+This file is part of FElupe.
+
+FElupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+FElupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
+"""
+from functools import wraps
+
+from jax.numpy import array, diag, log, sqrt, trace
+from jax.numpy.linalg import det, eigvalsh
+
+from ....tensortrax.models.hyperelastic import extended_tube as extended_tube_docstring
+
+
+@wraps(extended_tube_docstring)
+def extended_tube(C, Gc, delta, Ge, beta):
+    J3 = det(C) ** (-1 / 3)
+    D = J3 * trace(C)
+    λ1, λ2, λ3 = sqrt(J3 * eigvalsh(C + diag(array([0, 1e-4, -1e-4]))))
+    β = beta
+    δ = delta
+    γ = (1 - δ**2) * (D - 3) / (1 - δ**2 * (D - 3))
+    Wc = Gc / 2 * (γ + log(1 - δ**2 * (D - 3)))
+    We = 2 * Ge / β**2 * (λ1**-β + λ2**-β + λ3**-β - 3)
+    return Wc + We

src/felupe/constitution/jax/models/hyperelastic/_neo_hooke.py

@@ -0,0 +1,28 @@
+# -*- coding: utf-8 -*-
+"""
+This file is part of FElupe.
+
+FElupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+FElupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
+"""
+from functools import wraps
+
+from jax.numpy import trace
+from jax.numpy.linalg import det
+
+from ....tensortrax.models.hyperelastic import neo_hooke as neo_hooke_docstring
+
+
+@wraps(neo_hooke_docstring)
+def neo_hooke(C, mu):
+    return mu / 2 * (det(C) ** (-1 / 3) * trace(C) - 3)

src/felupe/constitution/jax/models/hyperelastic/_storakers.py

@@ -17,7 +17,7 @@
 """
 from functools import wraps
 
-from jax.numpy import array, sqrt
+from jax.numpy import array, diag, sqrt
 from jax.numpy import sum as asum
 from jax.numpy.linalg import eigvalsh
 
@@ -26,11 +26,13 @@
 
 @wraps(storakers_docstring)
 def storakers(C, mu, alpha, beta):
-    λ1, λ2, λ3 = sqrt(eigvalsh(C))
+    λ1, λ2, λ3 = sqrt(eigvalsh(C + diag(array([0, -1e-4, 1e-4]))))
     J = λ1 * λ2 * λ3
 
     μ = array(mu)
     α = array(alpha)
     β = array(beta)
 
-    return asum(2 * μ / α**2 * (λ1**α + λ2**α + λ3**α - 3 + (J ** (-α * β) - 1) / β))
+    return asum(
+        2 * μ / α**2 * (λ1**α + λ2**α + λ3**α - 3 + (J ** (-α * β) - 1) / β)
+    )

src/felupe/constitution/jax/models/hyperelastic/_van_der_waals.py

@@ -17,8 +17,7 @@
 """
 from functools import wraps
 
-from jax.lax import select
-from jax.numpy import isclose, log, sqrt, trace
+from jax.numpy import log, sqrt, trace
 from jax.numpy.linalg import det
 
 from ....tensortrax.models.hyperelastic import van_der_waals as van_der_waals_docstring
@@ -30,7 +29,7 @@ def van_der_waals(C, mu, limit, a, beta):
     I1 = J3 * trace(C)
     I2 = (trace(C) ** 2 - J3**2 * trace(C @ C)) / 2
     Im = (1 - beta) * I1 + beta * I2
-    Im = select(isclose(Im, 3), Im + 1e-6, Im)
+    Im += 1e-4
     eta = sqrt((Im - 3) / (limit**2 - 3))
     return mu * (
         -(limit**2 - 3) * (log(1 - eta) + eta) - 2 / 3 * a * ((Im - 3) / 2) ** (3 / 2)

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

@@ -23,6 +23,7 @@
 from ._ogden import ogden
 from ._ogden_roxburgh import ogden_roxburgh
 from ._saint_venant_kirchhoff import saint_venant_kirchhoff
+from ._saint_venant_kirchhoff_orthotropic import saint_venant_kirchhoff_orthotropic
 from ._storakers import storakers
 from ._third_order_deformation import third_order_deformation
 from ._van_der_waals import van_der_waals
@@ -43,6 +44,7 @@
     "ogden",
     "ogden_roxburgh",
     "saint_venant_kirchhoff",
+    "saint_venant_kirchhoff_orthotropic",
     "storakers",
     "third_order_deformation",
     "van_der_waals",
@@ -63,6 +65,13 @@
 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_orthotropic.kwargs = dict(
+    mu=[0.0, 0.0, 0.0],
+    lmbda=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+    r1=[1.0, 0.0, 0.0],
+    r2=[0.0, 1.0, 0.0],
+    r3=[0.0, 0.0, 1.0],
+)
 storakers.kwargs = dict(mu=[0], alpha=[2], beta=[1])
 third_order_deformation.kwargs = dict(C10=0, C01=0, C11=0, C20=0, C30=0)
 van_der_waals.kwargs = dict(mu=0, beta=0, a=0, limit=100)

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

@@ -28,7 +28,7 @@ def extended_tube(C, Gc, delta, Ge, beta):
 
     Parameters
     ----------
-    C : tensortrax.Tensor
+    C : tensortrax.Tensor or jax.Array
         Right Cauchy-Green deformation tensor.
     Gc : float
         Cross-link contribution to the initial shear modulus.
@@ -82,14 +82,27 @@ def extended_tube(C, Gc, delta, Ge, beta):
 
     Examples
     --------
+    First, choose the desired automatic differentiation backend
+
+    ..  pyvista-plot::
+        :context:
+
+        >>> # import felupe.constitution.jax as mat
+        >>> import felupe.constitution.tensortrax as mat
+
+    and create the hyperelastic material.
 
     ..  pyvista-plot::
         :context:
 
         >>> import felupe as fem
         >>>
-        >>> umat = fem.Hyperelastic(
-        ...     fem.extended_tube, Gc=0.1867, Ge=0.2169, beta=0.2, delta=0.09693
+        >>> umat = mat.Hyperelastic(
+        ...     mat.models.hyperelastic.extended_tube,
+        ...     Gc=0.1867,
+        ...     Ge=0.2169,
+        ...     beta=0.2,
+        ...     delta=0.09693,
         ... )
         >>> ax = umat.plot(incompressible=True)
 

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

@@ -26,7 +26,7 @@ def neo_hooke(C, mu):
 
     Parameters
     ----------
-    C : tensortrax.Tensor
+    C : tensortrax.Tensor or jax.Array
         Right Cauchy-Green deformation tensor.
     mu : float
         Shear modulus.
@@ -42,13 +42,22 @@ def neo_hooke(C, mu):
 
     Examples
     --------
+    First, choose the desired automatic differentiation backend
+
+    ..  pyvista-plot::
+        :context:
+
+        >>> # import felupe.constitution.jax as mat
+        >>> import felupe.constitution.tensortrax as mat
+
+    and create the hyperelastic material.
 
     ..  pyvista-plot::
         :context:
 
         >>> import felupe as fem
         >>>
-        >>> umat = fem.Hyperelastic(fem.neo_hooke, mu=1.0)
+        >>> umat = mat.Hyperelastic(mat.models.hyperelastic.neo_hooke, mu=1.0)
         >>> ax = umat.plot(incompressible=True)
 
     ..  pyvista-plot::