adtzlr · adtzlr · May 14, 2024 · May 14, 2024 · May 14, 2024 · May 14, 2024
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -11,6 +11,7 @@ All notable changes to this project will be documented in this file. The format
 - 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()`.
 - Add the MORPH material formulation `morph(p)` to be used in `Hyperelastic()`.
+- Add a decorator `@total_lagrange` for Total Lagrange material formulations to be used in `Hyperelastic`.
 
 ### Changed
 - Recfactor the `constitution` module.

diff --git a/docs/felupe/constitution.rst b/docs/felupe/constitution.rst
@@ -67,6 +67,8 @@ There are many different pre-defined constitutive material formulations availabl
 
    Hyperelastic
    MaterialAD
+   total_lagrange
+   updated_lagrange
 
 **Material Models (Strain Energy Functions) for** :class:`~felupe.Hyperelastic`
 
@@ -78,7 +80,6 @@ There are many different pre-defined constitutive material formulations availabl
    finite_strain_viscoelastic
    miehe_goektepe_lulei
    mooney_rivlin
-   morph
    morph_representative_directions
    neo_hooke
    ogden
@@ -88,6 +89,12 @@ There are many different pre-defined constitutive material formulations availabl
    van_der_waals
    yeoh
 
+**Material Models (Strain Energy Functions) for** :class:`~felupe.MaterialAD`
+
+.. autosummary::
+
+   morph
+
 **Small Strain-based User Materials**
 
 .. autosummary::
@@ -244,6 +251,10 @@ There are many different pre-defined constitutive material formulations availabl
 
 .. autofunction:: felupe.third_order_deformation
 
+.. autofunction:: felupe.total_lagrange
+
+.. autofunction:: felupe.updated_lagrange
+
 .. autofunction:: felupe.van_der_waals
 
 .. autofunction:: felupe.yeoh

diff --git a/examples/ex13_morph-rubber-wheel.py b/examples/ex13_morph-rubber-wheel.py
@@ -6,7 +6,7 @@
 [1]_. While the rotation is increased, a constant vertical compression is applied to the
 rubber wheel by a frictionless contact on the bottom. The vertical reaction force is
 then carried out for the rotation angles. The MORPH material model is implemented as a
-first Piola-Kirchhoff stress-based formulation with automatic differentiation. The
+second Piola-Kirchhoff stress-based formulation with automatic differentiation. The
 Tresca invariant of the distortional part of the right Cauchy-Green deformation tensor
 is used as internal state variable, see Eq. :eq:`morph-state`.
 
@@ -102,6 +102,7 @@
 import felupe as fem
 
 
+@fem.total_lagrange
 def morph(F, statevars_old, p):
     "MORPH material model formulation."
 
@@ -111,7 +112,7 @@ def morph(F, statevars_old, p):
     # extract old state variables
     CTSn = tm.array(statevars_old[0], like=C[0, 0])
     Cn = tm.special.from_triu_1d(statevars_old[1:7], like=C)
-    SAn = tm.special.from_triu_1d(statevars_old[7:], like=C)
+    SAn = tm.special.from_triu_1d(statevars_old[7:13], like=C)
 
     # distortional part of right Cauchy-Green deformation tensor
     I3 = tm.linalg.det(C)
@@ -152,12 +153,14 @@ def f(x):
     S = 2 * α * tm.special.dev(CG) @ invC + tm.special.dev(SA @ C) @ invC
 
     try:  # update state variables
-        statevars_new = tm.stack([CTS, *tm.special.triu_1d(C), *tm.special.triu_1d(SA)])
+        statevars_new = np.stack(
+            [CTS.x, *tm.special.triu_1d(C).x, *tm.special.triu_1d(SA).x]
+        )
     except:
         # not possible (and not necessary) during AD-based hessian evaluation
         statevars_new = statevars_old
 
-    return F @ S, statevars_new
+    return S, statevars_new
 
 
 umat = fem.MaterialAD(

diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py
@@ -53,6 +53,8 @@
     ogden_roxburgh,
     saint_venant_kirchhoff,
     third_order_deformation,
+    total_lagrange,
+    updated_lagrange,
     van_der_waals,
     yeoh,
 )
@@ -174,6 +176,8 @@
     "Material",
     "MaterialStrain",
     "Hyperelastic",
+    "total_lagrange",
+    "updated_lagrange",
     "MaterialAD",
     "ViewMaterial",
     "ViewMaterialIncompressible",

diff --git a/src/felupe/constitution/__init__.py b/src/felupe/constitution/__init__.py
@@ -1,7 +1,6 @@
 from ._base import CompositeMaterial, ConstitutiveMaterial, constitutive_material
 from ._kinematics import AreaChange, LineChange, VolumeChange
 from ._material import Material
-from ._material_ad import MaterialAD
 from ._mixed import NearlyIncompressible, ThreeFieldVariation
 from ._view import ViewMaterial, ViewMaterialIncompressible
 from .hyperelasticity import Hyperelastic
@@ -19,7 +18,6 @@
     isochoric_volumetric_split,
     miehe_goektepe_lulei,
     mooney_rivlin,
-    morph,
     morph_representative_directions,
     neo_hooke,
     ogden,
@@ -29,6 +27,8 @@
     van_der_waals,
     yeoh,
 )
+from .lagrange import MaterialAD, total_lagrange, updated_lagrange
+from .lagrange.models import morph
 from .linear_elasticity import (
     LinearElastic,
     LinearElasticLargeStrain,
@@ -45,7 +45,6 @@
 )
 
 __all__ = [
-    "hyperelastic",
     "alexander",
     "arruda_boyce",
     "extended_tube",
@@ -76,6 +75,8 @@
     "MaterialStrain",
     "MaterialAD",
     "Hyperelastic",
+    "total_lagrange",
+    "updated_lagrange",
     "AreaChange",
     "LineChange",
     "VolumeChange",

diff --git a/src/felupe/constitution/hyperelasticity/models/__init__.py b/src/felupe/constitution/hyperelasticity/models/__init__.py
@@ -16,7 +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, morph_representative_directions
+from ._morph import morph_representative_directions
 from ._neo_hooke import neo_hooke
 from ._ogden import ogden
 from ._ogden_roxburgh import ogden_roxburgh
@@ -33,7 +33,6 @@
     "isochoric_volumetric_split",
     "miehe_goektepe_lulei",
     "mooney_rivlin",
-    "morph",
     "morph_representative_directions",
     "neo_hooke",
     "ogden",

diff --git a/src/felupe/constitution/hyperelasticity/models/_morph.py b/src/felupe/constitution/hyperelasticity/models/_morph.py
@@ -15,216 +15,11 @@
 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,
-    if_else,
-    real_to_dual,
-)
-from tensortrax.math.special import try_stack, from_triu_1d, triu_1d, sym, dev, ddot
-from tensortrax.math.linalg import det, inv, eigvalsh, expm
-
-
-def morph(C, statevars, p):
-    r"""Strain energy function of the
-    `MORPH <https://doi.org/10.1016/s0749-6419(02)00091-8>`_ model formulation [1]_.
-
-    Parameters
-    ----------
-    C : tensortrax.Tensor
-        Right Cauchy-Green deformation tensor.
-    statevars : array
-        Vector of stacked state variables (CTS, C, SA).
-    p : list of float
-        A list which contains the 8 material parameters.
-
-    Notes
-    -----
-    The MORPH material model is implemented as a second Piola-Kirchhoff stress-based
-    formulation with automatic differentiation. The Tresca invariant of the distortional
-    part of the right Cauchy-Green deformation tensor is used as internal state
-    variable, see Eq. :eq:`morph-state`.
-
-    ..  warning::
-        While the `MORPH <https://doi.org/10.1016/s0749-6419(02)00091-8>`_-material
-        formulation captures the Mullins effect and quasi-static hysteresis effects of
-        rubber mixtures very nicely, it has been observed to be unstable for medium- to
-        highly-distorted states of deformation.
-
-    ..  math::
-        :label: morph-state
-
-        \boldsymbol{C} &= \boldsymbol{F}^T \boldsymbol{F}
-
-        I_3 &= \det (\boldsymbol{C})
-
-        \hat{\boldsymbol{C}} &= I_3^{-1/3} \boldsymbol{C}
-
-        \hat{\lambda}^2_\alpha &= \text{eigvals}(\hat{\boldsymbol{C}})
-
-        \hat{C}_T &= \max \left( \hat{\lambda}^2_\alpha - \hat{\lambda}^2_\beta \right)
-
-        \hat{C}_T^S &= \max \left( \hat{C}_T, \hat{C}_{T,n}^S \right)
-
-    A sigmoid-function is used inside the deformation-dependent variables
-    :math:`\alpha`, :math:`\beta` and :math:`\gamma`, see Eq. :eq:`morph-sigmoid`.
-
-    ..  math::
-        :label: morph-sigmoid
-
-        f(x) &= \frac{1}{\sqrt{1 + x^2}}
-
-        \alpha &= p_1 + p_2 \ f(p_3\ C_T^S)
+from tensortrax.math import abs as tensor_abs
+from tensortrax.math import array, exp, maximum, sqrt
+from tensortrax.math.special import try_stack
 
-        \beta &= p_4\ f(p_3\ C_T^S)
-
-        \gamma &= p_5\ C_T^S\ \left( 1 - f\left(\frac{C_T^S}{p_6}\right) \right)
-
-    The rate of deformation is described by the Lagrangian tensor and its Tresca-
-    invariant, see Eq. :eq:`morph-rate-of-deformation`.
-
-    ..  note::
-        It is important to evaluate the incremental right Cauchy-Green tensor by the
-        difference of the final and the previous state of deformation, not by its
-        variation with respect to the deformation gradient tensor.
-
-    ..  math::
-        :label: morph-rate-of-deformation
-
-        \hat{\boldsymbol{L}} &= \text{sym}\left(
-                \text{dev}(\boldsymbol{C}^{-1} \Delta\boldsymbol{C})
-            \right) \hat{\boldsymbol{C}}
-
-        \lambda_{\hat{\boldsymbol{L}}, \alpha} &= \text{eigvals}(\hat{\boldsymbol{L}})
-
-        \hat{L}_T &= \max \left(
-            \lambda_{\hat{\boldsymbol{L}}, \alpha}-\lambda_{\hat{\boldsymbol{L}}, \beta}
-        \right)
-
-        \Delta\boldsymbol{C} &= \boldsymbol{C} - \boldsymbol{C}_n
-
-    The additional stresses evolve between the limiting stresses, see Eq.
-    :eq:`morph-stresses`. The additional deviatoric-enforcement terms [1]_ are neglected
-    in this implementation.
-
-    ..  math::
-        :label: morph-stresses
-
-        \boldsymbol{S}_L &= \left(
-            \gamma \exp \left(p_7 \frac{\hat{\boldsymbol{L}}}{\hat{L}_T}
-                \frac{\hat{C}_T}{\hat{C}_T^S} \right) +
-                p8 \frac{\hat{\boldsymbol{L}}}{\hat{L}_T}
-        \right) \boldsymbol{C}^{-1}
-
-        \boldsymbol{S}_A &= \frac{
-            \boldsymbol{S}_{A,n} + \beta\ \hat{L}_T\ \boldsymbol{S}_L
-        }{1 + \beta\ \hat{L}_T}
-
-        \boldsymbol{S} &= 2 \alpha\ \text{dev}( \hat{\boldsymbol{C}} )
-            \boldsymbol{C}^{-1}+\text{dev}\left(\boldsymbol{S}_A\ \boldsymbol{C}\right)
-            \boldsymbol{C}^{-1}
-
-    ..  note::
-        Only the upper-triangle entries of the symmetric stress-tensor state
-        variables are stored in the solid body. Hence, it is necessary to extract such
-        variables with :func:`tm.special.from_triu_1d` and export them as
-        :func:`tm.special.triu_1d`.
-
-    Examples
-    --------
-    ..  pyvista-plot::
-        :context:
-
-        >>> import felupe as fem
-        >>>
-        >>> umat = fem.Hyperelastic(
-        ...     fem.morph,
-        ...     p=[0.039, 0.371, 0.174, 2.41, 0.0094, 6.84, 5.65, 0.244],
-        ...     nstatevars=13,
-        ... )
-        >>> 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],
-        ...        [1, 2.75, 1, 2.75],
-        ...        num=20,
-        ...    ),
-        ...    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>`_.
-
-    See Also
-    --------
-    felupe.morph_representative_directions : Strain energy function of the MORPH model
-        formulation, implemented by the concept of representative directions.
-    """
-
-    # extract old state variables
-    CTSn = array(statevars[0], like=C[0, 0])
-    Cn = from_triu_1d(statevars[1:7], like=C)
-    SAn = from_triu_1d(statevars[7:13], like=C)
-
-    # distortional part of right Cauchy-Green deformation tensor
-    I3 = det(C)
-    CG = C * I3 ** (-1 / 3)
-
-    # inverse of and incremental right Cauchy-Green deformation tensor
-    invC = inv(C)
-    dC = C - Cn
-
-    # eigenvalues of right Cauchy-Green deformation tensor (sorted in ascending order)
-    λCG = eigvalsh(CG)
-
-    # Tresca invariant of distortional part of right Cauchy-Green deformation tensor
-    CTG = λCG[-1] - λCG[0]
-
-    # maximum Tresca invariant in load history
-    CTS = maximum(CTG, CTSn)
-
-    def sigmoid(x):
-        "Algebraic sigmoid function."
-        return 1 / sqrt(1 + x**2)
-
-    # material parameters
-    α = p[0] + p[1] * sigmoid(p[2] * CTS)
-    β = p[3] * sigmoid(p[2] * CTS)
-    γ = p[4] * CTS * (1 - sigmoid(CTS / p[5]))
-
-    LG = sym(dev(invC @ dC)) @ CG
-    λLG = eigvalsh(LG)
-    LTG = λLG[-1] - λLG[0]
-
-    # limiting stresses "L" and additional stresses "A"
-    SL = (γ * expm(p[6] * LG / LTG * CTG / CTS) + p[7] * LG / LTG) @ invC
-    SA = (SAn + β * LTG * SL) / (1 + β * LTG)
-
-    # second Piola-Kirchhoff stress tensor
-    dψdC = α * dev(CG) @ invC + dev(SA @ C) @ invC / 2
-    statevars_new = try_stack([[CTS], triu_1d(C), triu_1d(SA)], fallback=statevars)
-
-    return real_to_dual(dψdC, C, mul=ddot), statevars_new
+from .microsphere import affine_stretch_statevars
 
 
 def morph_representative_directions(C, statevars, p, ε=1e-8):