Add morph() for JAX-based Material

.github/workflows/upload-codecov.yml

@@ -14,7 +14,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
     steps:
     - uses: actions/setup-python@v4
       with:

CHANGELOG.md

@@ -12,9 +12,11 @@ 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.autodiff.jax.Hyperelastic` as a feature-equivalent alternative to `Hyperelastic` with `jax` as backend.
 - Add `constitution.autodiff.jax.Material` as a feature-equivalent alternative to `MaterialAD` with `jax` as backend.
+- Add the MORPH-material formulation for a JAX-based material `felupe.constitution.autodiff.jax.models.lagrange.morph()`.
 
 ### Changed
 - Change default `np.einsum(..., order="K")` to `np.einsum(..., order="C")` in the methods of `Field`, `FieldAxisymmetric`, `FieldPlaneStrain` and `FieldContainer`.
+- Change supported Python versions to 3.9 - 3.12.
 
 ### Fixed
 - 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.

README.md

@@ -5,7 +5,7 @@
 
 [![FElupe](https://img.shields.io/badge/%F0%9F%94%8D-FElupe-white)](https://felupe.readthedocs.io) [![PyPI version shields.io](https://img.shields.io/pypi/v/felupe.svg)](https://pypi.python.org/pypi/felupe/) [![Documentation Status](https://readthedocs.org/projects/felupe/badge/?version=latest)](https://felupe.readthedocs.io/en/latest/?badge=latest) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![codecov](https://codecov.io/gh/adtzlr/felupe/branch/main/graph/badge.svg?token=J2QP6Y6LVH)](https://codecov.io/gh/adtzlr/felupe) [![DOI](https://zenodo.org/badge/360657894.svg)](https://zenodo.org/badge/latestdoi/360657894) ![Codestyle black](https://img.shields.io/badge/code%20style-black-black) ![PyPI - Downloads](https://img.shields.io/pypi/dm/felupe) [![lite-badge](https://jupyterlite.rtfd.io/en/latest/_static/badge.svg)](https://adtzlr.github.io/felupe-web/lab?path=01_Getting-Started.ipynb) <a target="_blank" href="https://colab.research.google.com/github/adtzlr/felupe-web/blob/main/notebooks/colab/01_Getting-Started.ipynb"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>
 
-FElupe is a Python 3.8+ 🐍 finite element analysis package 📦 focusing on the formulation and numerical solution of nonlinear problems in continuum mechanics of solid bodies 🔧. This package is intended for scientific research 💻, but is also suitable for running nonlinear simulations 🚂 in general 🏎️. In addition to the transformation of general weak forms into sparse vectors and matrices, FElupe provides an efficient high-level abstraction layer for the simulation of the deformation of solid bodies.
+FElupe is a Python 3.9+ 🐍 finite element analysis package 📦 focusing on the formulation and numerical solution of nonlinear problems in continuum mechanics of solid bodies 🔧. This package is intended for scientific research 💻, but is also suitable for running nonlinear simulations 🚂 in general 🏎️. In addition to the transformation of general weak forms into sparse vectors and matrices, FElupe provides an efficient high-level abstraction layer for the simulation of the deformation of solid bodies.
 
 <p align="center">
   <a href="https://felupe.readthedocs.io/en/latest/examples/ex01_beam.html"><img src="https://felupe.readthedocs.io/en/latest/_images/sphx_glr_ex01_beam_thumb.png" height="100px"/></a> <a href="https://felupe.readthedocs.io/en/latest/examples/ex02_plate-with-hole.html"><img src="https://felupe.readthedocs.io/en/latest/_images/sphx_glr_ex02_plate-with-hole_thumb.png" height="100px"/></a> <a

docs/felupe/constitution/jax.rst

@@ -14,6 +14,12 @@ This page contains material model formulations with automatic differentiation us
    constitution.autodiff.jax.Hyperelastic
    constitution.autodiff.jax.Material
 
+**Material Models for** :class:`felupe.constitution.autodiff.jax.Material`
+
+.. autosummary::
+
+   felupe.constitution.autodiff.jax.models.lagrange.morph
+
 **Tools**
 
 .. autosummary::
@@ -32,4 +38,6 @@ This page contains material model formulations with automatic differentiation us
    :undoc-members:
    :inherited-members:
 
+.. autofunction:: felupe.constitution.autodiff.jax.models.lagrange.morph
+
 .. autofunction:: felupe.constitution.autodiff.jax.vmap

docs/index.rst

@@ -1,7 +1,7 @@
 FElupe documentation
 ====================
 
-FElupe is a Python 3.8+ 🐍 finite element analysis package 📦 focusing on the formulation and numerical solution of nonlinear problems in continuum mechanics of solid bodies 🔧. This package is intended for scientific research 💻, but is also suitable for running nonlinear simulations 🚂 in general 🏎️. In addition to the transformation of general weak forms into sparse vectors and matrices, FElupe provides an efficient high-level abstraction layer for the simulation of the deformation of solid bodies.
+FElupe is a Python 3.9+ 🐍 finite element analysis package 📦 focusing on the formulation and numerical solution of nonlinear problems in continuum mechanics of solid bodies 🔧. This package is intended for scientific research 💻, but is also suitable for running nonlinear simulations 🚂 in general 🏎️. In addition to the transformation of general weak forms into sparse vectors and matrices, FElupe provides an efficient high-level abstraction layer for the simulation of the deformation of solid bodies.
 
 .. grid::
 
@@ -72,6 +72,7 @@ where ``[all]`` is a combination of ``[io,parallel,plot,progress,view]`` and ins
 In order to make use of all features of FElupe 💎💰💍👑💎, it is suggested to install all optional dependencies.
 
 * `einsumt <https://github.com/mrkwjc/einsumt>`_ for parallel (threaded) assembly
+* `jax <https://github.com/jax-ml/jax>`_ for JAX-based material formulations
 * `h5py <https://github.com/h5py/h5py>`_ for writing XDMF result files
 * `matplotlib <https://github.com/matplotlib/matplotlib>`_ for plotting graphs
 * `meshio <https://github.com/nschloe/meshio>`_ for mesh-related I/O

pyproject.toml

@@ -14,7 +14,7 @@ authors = [
   {email = "a.dutzler@gmail.com"},
 
 ]
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 license = {file = "LICENSE"}
 keywords = [
   "computational-mechanics",
@@ -38,7 +38,6 @@ classifiers = [
   "Operating System :: OS Independent",
   "Programming Language :: Python",
   "Programming Language :: Python :: 3",
-  "Programming Language :: Python :: 3.8",
   "Programming Language :: Python :: 3.9",
   "Programming Language :: Python :: 3.10",
   "Programming Language :: Python :: 3.11",

src/felupe/constitution/autodiff/jax/__init__.py

@@ -1,5 +1,6 @@
+from . import models
 from ._hyperelastic import Hyperelastic
 from ._material import Material
 from ._tools import vmap
 
-__all__ = ["Hyperelastic", "Material", "vmap"]
+__all__ = ["Hyperelastic", "Material", "models", "vmap"]

src/felupe/constitution/autodiff/jax/_material.py

@@ -100,7 +100,8 @@ def viscoelastic(F, Cin, mu, eta, dtime):
             S = mu * dev(Cu @ jnp.linalg.inv(Ci)) @ jnp.linalg.inv(C)
 
             # first Piola-Kirchhoff stress tensor and state variable
-            to_triu = lambda C: C[*jnp.triu_indices(3)]
+            i, j = triu_indices(3)
+            to_triu = lambda C: C[i, j]
             return F @ S, to_triu(Ci)
 
         umat = fem.constitution.autodiff.jax.Material(

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

@@ -0,0 +1,3 @@
+from . import hyperelastic, lagrange
+
+__all__ = ["hyperelastic", "lagrange"]

src/felupe/constitution/autodiff/jax/models/lagrange/__init__.py

@@ -0,0 +1,8 @@
+from ._morph import morph
+
+__all__ = [
+    "morph",
+]
+
+# default (stable) material parameters
+morph.kwargs = dict(p=[0, 0, 0, 0, 0, 1, 0, 0])

src/felupe/constitution/autodiff/jax/models/lagrange/_morph.py

@@ -0,0 +1,228 @@
+# -*- 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/>.
+"""
+
+
+def morph(F, statevars, p):
+    r"""Second Piola-Kirchhoff stress tensor 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)
+
+        \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}
+
+    Examples
+    --------
+    ..  pyvista-plot::
+        :context:
+
+        >>> import felupe as fem
+        >>>
+        >>> umat = fem.constitution.autodiff.jax.Material(
+        ...     fem.constitution.autodiff.jax.models.lagrange.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>`_.
+    """
+
+    from jax.numpy import (
+        array,
+        concatenate,
+        diag,
+        eye,
+        maximum,
+        sqrt,
+        trace,
+        triu_indices,
+    )
+    from jax.numpy.linalg import det, eigvalsh, inv
+    from jax.scipy.linalg import expm
+
+    # right Cauchy-Green deformation tensor
+    C = F.T @ F
+
+    # extract old state variables
+    CTSn = statevars[0]
+    from_triu = lambda C: C[array([[0, 1, 2], [1, 3, 4], [2, 4, 5]])]
+    Cn = from_triu(statevars[1:7])
+    SAn = from_triu(statevars[7:13])
+
+    # 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)
+    eigvalsh2 = lambda C: eigvalsh(C + diag(array([1e-4, -1e-4, 0])))
+    λCG = eigvalsh2(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]))
+
+    dev = lambda C: C - trace(C) / 3 * eye(3)
+    sym = lambda C: (C + C.T) / 2
+
+    LG = sym(dev(invC @ dC)) @ CG
+    λLG = eigvalsh2(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
+    S = 2 * α * dev(CG) @ invC + dev(SA @ C) @ invC
+
+    i, j = triu_indices(3)
+    to_triu = lambda C: C[i, j]
+    statevars_new = concatenate([array([CTS]), to_triu(C), to_triu(SA)])
+
+    return F @ S, statevars_new

tests/test_constitution_jax.py

@@ -193,9 +193,34 @@ def dWdF(F, statevars, C10, K):
         pass
 
 
+def test_material_included_jax_statevars():
+    try:
+        umat = fem.constitution.autodiff.jax.Material(
+            fem.constitution.autodiff.jax.models.lagrange.morph,
+            p=[0.039, 0.371, 0.174, 2.41, 0.0094, 6.84, 5.65, 0.244],
+            nstatevars=13,
+        )
+        mesh = fem.Cube(n=2)
+        region = fem.RegionHexahedron(mesh)
+        field = fem.FieldContainer([fem.Field(region, dim=3)])
+
+        boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)
+        solid = fem.SolidBody(umat=umat, field=field)
+
+        move = fem.math.linsteps([0, 1], num=3)
+        ramp = {boundaries["move"]: move}
+        step = fem.Step(items=[solid], ramp=ramp, boundaries=boundaries)
+        job = fem.Job(steps=[step])
+        job.evaluate(tol=1e-4)
+
+    except ModuleNotFoundError:
+        pass
+
+
 if __name__ == "__main__":
     test_vmap()
     test_hyperelastic_jax()
     test_hyperelastic_jax_statevars()
     test_material_jax()
     test_material_jax_statevars()
+    test_material_included_jax_statevars()