adtzlr · adtzlr · May 13, 2024 · May 13, 2024 · May 13, 2024 · May 13, 2024
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -9,6 +9,7 @@ All notable changes to this project will be documented in this file. The format
 - Add the pseudo-elastic `ogden_roxburgh(r, m, beta, material, **kwargs)` material model formulation to be used in `Hyperelastic()`.
 - Add an optional relative-residuals argument to `ConstitutiveMaterial.optimize(relative=False)`.
 - Add a class-decorator `constitutive_material(Msterial, name=None)`.
+- Add the MORPH material formulation implemented by the concept of representative directions `morph_representative_directions(p)` to be used in `Hyperelastic()`.
 
 ### Changed
 - Recfactor the `constitution` module.

diff --git a/docs/felupe/constitution.rst b/docs/felupe/constitution.rst
@@ -78,6 +78,7 @@ There are many different pre-defined constitutive material formulations availabl
    finite_strain_viscoelastic
    miehe_goektepe_lulei
    mooney_rivlin
+   morph_representative_directions
    neo_hooke
    ogden
    ogden_roxburgh
@@ -228,6 +229,8 @@ There are many different pre-defined constitutive material formulations availabl
 
 .. autofunction:: felupe.mooney_rivlin
 
+.. autofunction:: morph_representative_directions
+
 .. autofunction:: felupe.neo_hooke
 
 .. autofunction:: felupe.ogden

diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py
@@ -46,6 +46,7 @@
     linear_elastic_plastic_isotropic_hardening,
     miehe_goektepe_lulei,
     mooney_rivlin,
+    morph_representative_directions,
     neo_hooke,
     ogden,
     ogden_roxburgh,
@@ -187,6 +188,7 @@
     "isochoric_volumetric_split",
     "miehe_goektepe_lulei",
     "mooney_rivlin",
+    "morph_representative_directions",
     "neo_hooke",
     "ogden",
     "saint_venant_kirchhoff",

diff --git a/src/felupe/constitution/__init__.py b/src/felupe/constitution/__init__.py
@@ -19,6 +19,7 @@
     isochoric_volumetric_split,
     miehe_goektepe_lulei,
     mooney_rivlin,
+    morph_representative_directions,
     neo_hooke,
     ogden,
     ogden_roxburgh,
@@ -51,6 +52,7 @@
     "isochoric_volumetric_split",
     "miehe_goektepe_lulei",
     "mooney_rivlin",
+    "morph_representative_directions",
     "neo_hooke",
     "ogden",
     "ogden_roxburgh",

diff --git a/src/felupe/constitution/hyperelasticity/models/__init__.py b/src/felupe/constitution/hyperelasticity/models/__init__.py
@@ -16,6 +16,7 @@
 from ._helpers import isochoric_volumetric_split
 from ._miehe_goektepe_lulei import miehe_goektepe_lulei
 from ._mooney_rivlin import mooney_rivlin
+from ._morph import morph_representative_directions
 from ._neo_hooke import neo_hooke
 from ._ogden import ogden
 from ._ogden_roxburgh import ogden_roxburgh
@@ -32,6 +33,7 @@
     "isochoric_volumetric_split",
     "miehe_goektepe_lulei",
     "mooney_rivlin",
+    "morph_representative_directions",
     "neo_hooke",
     "ogden",
     "ogden_roxburgh",
@@ -47,6 +49,7 @@
 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)
+morph_representative_directions.kwargs = dict(p=[0, 0, 0, 0, 0, 1, 0, 0])
 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)

diff --git a/src/felupe/constitution/hyperelasticity/models/_morph.py b/src/felupe/constitution/hyperelasticity/models/_morph.py
@@ -0,0 +1,154 @@
+# -*- 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 .microsphere import affine_stretch_statevars
+from tensortrax.math import array, abs as tensor_abs, maximum, sqrt, exp
+from tensortrax.math.special import try_stack
+
+
+def morph_representative_directions(C, statevars, p, ε=1e-8):
+    """Strain energy function of the
+    `MORPH <https://doi.org/10.1016/s0749-6419(02)00091-8>`_ model formulation [1]_,
+    implemented by the concept of
+    `representative directions <https://nbn-resolving.org/urn:nbn:de:bsz:ch1-qucosa-114428>`_
+    [2]_, [3]_.
+
+    Parameters
+    ----------
+    C : tensortrax.Tensor
+        Right Cauchy-Green deformation tensor.
+    statevars : array
+        Vector of stacked state variables (CTS, λ - 1, SA1, SA2).
+    p : list of float
+        A list which contains the 8 material parameters.
+    ε : float, optional
+        A small stabilization parameter (default is 1e-8).
+
+    Examples
+    --------
+    ..  pyvista-plot::
+        :context:
+
+        >>> import felupe as fem
+        >>>
+        >>> umat = fem.Hyperelastic(
+        ...     fem.morph_representative_directions,
+        ...     p=[0.011, 0.408, 0.421, 6.85, 0.0056, 5.54, 5.84, 0.117],
+        ...     nstatevars=84,
+        ... )
+        >>> ax = umat.plot(
+        ...    incompressible=True,
+        ...    ux=fem.math.linsteps(
+        ...        [1, 2, 1, 2.75, 1, 3.5, 1, 4.2, 1, 4.8, 1, 4.8, 1],
+        ...        num=50,
+        ...    ),
+        ...    ps=None,
+        ...    bx=None,
+        ... )
+
+    ..  pyvista-plot::
+        :include-source: False
+        :context:
+        :force_static:
+
+        >>> import pyvista as pv
+        >>>
+        >>> fig = ax.get_figure()
+        >>> chart = pv.ChartMPL(fig)
+        >>> chart.show()
+
+    References
+    ----------
+    .. [1] D. Besdo and J. Ihlemann, "A phenomenological constitutive model for
+       rubberlike materials and its numerical applications", International Journal
+       of Plasticity, vol. 19, no. 7. Elsevier BV, pp. 1019–1036, Jul. 2003. doi:
+       `10.1016/s0749-6419(02)00091-8 <https://doi.org/10.1016/s0749-6419(02)00091-8>`_.
+
+    .. [2] M. Freund, "Verallgemeinerung eindimensionaler Materialmodelle für die
+       Finite-Elemente-Methode", Dissertation, Technische Universität Chemnitz,
+       Chemnitz, 2013.
+
+    .. [3] C. Miehe, S. Göktepe and F. Lulei, "A micro-macro approach to rubber-like
+       materials - Part I: the non-affine micro-sphere model of rubber elasticity",
+       Journal of the Mechanics and Physics of Solids, vol. 52, no. 11. Elsevier BV, pp.
+       2617–2660, Nov. 2004. doi:
+       `10.1016/j.jmps.2004.03.011 <https://doi.org/10.1016/j.jmps.2004.03.011>`_.
+    """
+
+    def morph_uniaxial(λ, statevars, p, ε=1e-8):
+        """Return the force (per undeformed area) for a given longitudinal stretch in
+        uniaxial incompressible tension or compression for the MORPH material
+        formulation [1]_, [2]_.
+
+        Parameters
+        ----------
+        λ : tensortrax.Tensor
+            Longitudinal stretch of uniaxial incompressible deformation.
+        statevars : array
+            Vector of stacked state variables (CTS, λ - 1, SA1, SA2).
+        p : list of float
+            A list which contains the 8 material parameters.
+        ε : float, optional
+            A small stabilization parameter (default is 1e-8).
+
+        References
+        ----------
+        .. [1] D. Besdo and J. Ihlemann, "A phenomenological constitutive model for
+           rubberlike materials and its numerical applications", International Journal
+           of Plasticity, vol. 19, no. 7. Elsevier BV, pp. 1019–1036, Jul. 2003. doi:
+           `10.1016/s0749-6419(02)00091-8 <https://doi.org/10.1016/s0749-6419(02)00091-8>`_.
+
+        .. [2] M. Freund, "Verallgemeinerung eindimensionaler Materialmodelle für die
+           Finite-Elemente-Methode", Dissertation, Technische Universität Chemnitz,
+           Chemnitz, 2013.
+
+        """
+        CTSn = array(statevars[:21], like=C, shape=(21,))
+        λn = array(statevars[21:42], like=C, shape=(21,)) + 1
+        SA1n = array(statevars[42:63], like=C, shape=(21,))
+        SA2n = array(statevars[63:], like=C, shape=(21,))
+
+        CT = tensor_abs(λ**2 - 1 / λ)
+        CTS = maximum(CT, CTSn)
+
+        L1 = 2 * (λ**3 / λn - λn**2) / 3
+        L2 = (λn**2 / λ**3 - 1 / λn) / 3
+        LT = tensor_abs(L1 - L2)
+
+        sigmoid = lambda x: 1 / sqrt(1 + x**2)
+        α = p[0] + p[1] * sigmoid(p[2] * CTS)
+        β = p[3] * sigmoid(p[2] * CTS)
+        γ = p[4] * CTS * (1 - sigmoid(CTS / p[5]))
+
+        L1_LT = L1 / (ε + LT)
+        L2_LT = L2 / (ε + LT)
+        CT_CTS = CT / (ε + CTS)
+
+        SL1 = (γ * exp(p[6] * L1_LT * CT_CTS) + p[7] * L1_LT) / λ**2
+        SL2 = (γ * exp(p[6] * L2_LT * CT_CTS) + p[7] * L2_LT) * λ
+
+        SA1 = (SA1n + β * LT * SL1) / (1 + β * LT)
+        SA2 = (SA2n + β * LT * SL2) / (1 + β * LT)
+
+        dψdλ = (2 * α + SA1) * λ - (2 * α + SA2) / λ**2
+        statevars_new = try_stack([CTS, (λ - 1), SA1, SA2], fallback=statevars)
+
+        return 5 * dψdλ.real_to_dual(λ), statevars_new
+
+    return affine_stretch_statevars(
+        C, statevars, f=morph_uniaxial, kwargs={"p": p, "ε": ε}
+    )
diff --git a/src/felupe/constitution/hyperelasticity/models/_ogden_roxburgh.py b/src/felupe/constitution/hyperelasticity/models/_ogden_roxburgh.py
@@ -16,7 +16,7 @@
 along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
 """
 from tensortrax import Tensor, Δ, Δδ, f, δ
-from tensortrax.math import array, maximum, tanh
+from tensortrax.math import array, maximum, real_to_dual, tanh
 
 
 def ogden_roxburgh(C, Wmax_n, material, r, m, beta, **kwargs):
@@ -106,19 +106,11 @@ def ogden_roxburgh(C, Wmax_n, material, r, m, beta, **kwargs):
     W = material(C, **kwargs)
     Wmax = maximum(W, array(Wmax_n, like=W))
 
-    def pseudo_elastic_strain_energy(W, Wmax):
-        "Strain energy density function with custom variations."
+    # evolution equation
+    η = lambda W: 1 - 1 / r * tanh((Wmax - W) / (m + beta * Wmax))
 
-        # evolution equation
-        η = 1 - 1 / r * tanh((Wmax - W) / (m + beta * Wmax))
+    # custom first- and second-partial derivatives
+    # set the variation to δη * W (and the linearization to Δη * δW + η * ΔδW)
+    dWdF = real_to_dual(η(W), W)
 
-        # custom first- and second-partial derivatives
-        return Tensor(
-            x=f(η) * f(W),
-            δx=f(η) * δ(W),
-            Δx=f(η) * Δ(W),
-            Δδx=δ(η) * Δ(W) + f(η) * Δδ(W),
-            ntrax=W.ntrax,
-        )
-
-    return pseudo_elastic_strain_energy(W, Wmax), Wmax
+    return dWdF, Wmax
diff --git a/src/felupe/constitution/hyperelasticity/models/microsphere/__init__.py b/src/felupe/constitution/hyperelasticity/models/microsphere/__init__.py
@@ -1,9 +1,10 @@
 from ._chain import langevin, linear
-from ._framework_affine import affine_stretch, affine_tube
+from ._framework_affine import affine_stretch, affine_stretch_statevars, affine_tube
 from ._framework_nonaffine import nonaffine_stretch, nonaffine_tube
 
 __all__ = [
     "affine_stretch",
+    "affine_stretch_statevars",
     "affine_tube",
     "linear",
     "langevin",

diff --git a/src/felupe/constitution/hyperelasticity/models/microsphere/_framework_affine.py b/src/felupe/constitution/hyperelasticity/models/microsphere/_framework_affine.py
@@ -16,6 +16,19 @@ def affine_stretch(C, f, kwargs, quadrature=BazantOh(n=21)):
     return einsum("a...,a->...", f(λ, **kwargs), w)
 
 
+def affine_stretch_statevars(C, statevars, f, kwargs, quadrature=BazantOh(n=21)):
+    "Micro-sphere model: Affine stretch part."
+
+    r = quadrature.points
+    w = quadrature.weights
+
+    # affine stretches (distortional part)
+    λ = det(C) ** (-1 / 6) * sqrt(einsum("ai,ij...,aj->a...", r, C, r))
+    ψ, statevars_new = f(λ, statevars, **kwargs)
+
+    return einsum("a...,a->...", ψ, w), statevars_new
+
+
 def affine_tube(C, f, kwargs, quadrature=BazantOh(n=21)):
     "Micro-sphere model: Affine area-stretch part."
 

diff --git a/src/felupe/mesh/_geometry.py b/src/felupe/mesh/_geometry.py
@@ -295,11 +295,11 @@ def __init__(self, *xi, indexing="ij", **kwargs):
 class RectangleArbitraryOrderQuad(Rectangle):
     """A rectangular 2d-mesh with an arbitrarr-order Lagrange quad between ``a`` and
     ``b``.
-    
+
     Examples
     --------
     ..  pyvista-plot::
-        
+
         >>> import felupe as fem
         >>>
         >>> mesh = fem.mesh.RectangleArbitraryOrderQuad(order=5).add_runouts()
@@ -317,11 +317,11 @@ def __init__(self, a=(0.0, 0.0), b=(1.0, 1.0), order=2):
 class CubeArbitraryOrderHexahedron(Cube):
     """A rectangular 2d-mesh with an arbitrarr-order Lagrange hexahedron between ``a``
     and ``b``.
-    
+
     Examples
     --------
     ..  pyvista-plot::
-        
+
         >>> import felupe as fem
         >>>
         >>> mesh = fem.mesh.CubeArbitraryOrderHexahedron(order=5).add_runouts()

diff --git a/tests/test_constitution.py b/tests/test_constitution.py
@@ -299,35 +299,38 @@ def elasticity(x, mu, lmbda):
     dsde = linear_elastic.hessian([F, None])
 
 
-def test_umat_ogden_roxburgh():
+def test_umat_hyperelastic_statevars():
     r, x = pre(sym=False, add_identity=True)
     F = x[0]
 
-    import matplotlib.pyplot as plt
-    import tensortrax.math as tm
-
     for model, kwargs, incompressible in [
         (
             fem.constitution.ogden_roxburgh,
             dict(r=3, m=1, beta=0, material=fem.neo_hooke, mu=1, nstatevars=1),
             True,
-        )
+        ),
+        (
+            fem.constitution.morph_representative_directions,
+            dict(
+                p=[0.011, 0.408, 0.421, 6.85, 0.0056, 5.54, 5.84, 0.117], nstatevars=84
+            ),
+            True,
+        ),
     ]:
         umat = fem.Hyperelastic(model, **kwargs)
 
-        statevars = np.zeros((8, 1))
+        statevars = np.zeros((kwargs["nstatevars"], 8, 1))
         s, statevars_new = umat.gradient([F, statevars])
         dsde = umat.hessian([F, statevars])
 
-        ux = fem.math.linsteps([1, 1.5, 1, 2, 1, 2.5, 1], num=[10, 10, 20, 20, 30, 30])
+        ux = fem.math.linsteps([1, 2, 1], num=10)
         ax = umat.plot(ux=ux, bx=None, ps=None, incompressible=True)
 
 
 def test_umat_hyperelastic():
     r, x = pre(sym=False, add_identity=True)
     F = x[0]
 
-    import matplotlib.pyplot as plt
     import tensortrax.math as tm
 
     def neo_hooke(C, mu=1):
@@ -609,7 +612,7 @@ def test_optimize():
     test_umat()
     test_umat_hyperelastic()
     test_umat_hyperelastic2()
-    test_umat_ogden_roxburgh()
+    test_umat_hyperelastic_statevars()
     test_umat_viscoelastic()
     test_umat_viscoelastic2()
     test_umat_strain()