Materials with Automatic Differentiation (JAX)#

This page contains material model formulations with automatic differentiation using jax.

Note

JAX uses single-precision (32bit) data types by default. This requires to relax the tolerance of newtonrhapson() to tol=1e-4. If required, JAX may be enforced to use double-precision at startup with jax.config.update("jax_enable_x64", True).

Note

The number of local XLA devices available must be greater or equal the number of the parallel-mapped axis, i.e. the number of quadrature points per cell when used in Material and Hyperelastic along with parallel=True. To use the multiple cores of a CPU device as multiple local XLA devices, the XLA device count must be defined at startup.

import os

os.environ["XLA_FLAGS"] = "--xla_force_host_platform_device_count=4"

Frameworks

Hyperelastic(fun[, nstatevars, jit, parallel])

A hyperelastic material definition with a given function for the strain energy density function per unit undeformed volume with Automatic Differentiation provided by jax.

Material(fun[, nstatevars, jit, parallel, ...])

A material definition with a given function for the partial derivative of the strain energy function w.r.t.

total_lagrange(material)

Decorate a second Piola-Kirchhoff stress Total-Lagrange material formulation as a first Piola-Kirchoff stress function.

updated_lagrange(material)

Decorate a Cauchy-stress Updated-Lagrange material formulation as a first Piola- Kirchoff stress function.

Material Models for felupe.constitution.jax.Hyperelastic

These material model formulations are defined by a strain energy density function.

blatz_ko(C, mu)

Strain energy function of the Blatz-Ko isotropic hyperelastic foam material formulation [1]_.

extended_tube(C, Gc, delta, Ge, beta)

Strain energy function of the isotropic hyperelastic Extended Tube [1]_ material formulation.

miehe_goektepe_lulei(C, mu, N, U, p, q)

Strain energy function of the isotropic hyperelastic micro-sphere model formulation [1]_.

mooney_rivlin(C, C10, C01)

Strain energy function of the isotropic hyperelastic Mooney-Rivlin material formulation.

neo_hooke(C, mu)

Strain energy function of the isotropic hyperelastic Neo-Hookean material formulation.

storakers(C, mu, alpha, beta)

third_order_deformation(C, C10, C01, C11, ...)

Strain energy function of the isotropic hyperelastic Third-Order-Deformation material formulation.

van_der_waals(C, mu, limit, a, beta)

Strain energy function of the Van der Waals [1]_ material formulation.

yeoh(C, C10, C20, C30)

Strain energy function of the isotropic hyperelastic Yeoh material formulation.

Material Models for felupe.constitution.jax.Material

The material model formulations are defined by the first Piola-Kirchhoff stress tensor. Function-decorators are available to use Total-Lagrange and Updated-Lagrange material formulations in Material.

morph(F, statevars, p)

Second Piola-Kirchhoff stress tensor of the MORPH model formulation [1]_.

morph_representative_directions(F, ...)

Tools

vmap(fun[, in_axes, out_axes, method])

Vectorizing map.

Detailed API Reference

class felupe.constitution.jax.Hyperelastic(fun, nstatevars=0, jit=True, parallel=False, **kwargs)[source]#

A hyperelastic material definition with a given function for the strain energy density function per unit undeformed volume with Automatic Differentiation provided by jax.

Parameters:
  • fun (callable) – A strain energy density function in terms of the right Cauchy-Green deformation tensor \(\boldsymbol{C}\). Function signature must be fun = lambda C, **kwargs: psi for functions without state variables and fun = lambda C, statevars, **kwargs: [psi, statevars_new] for functions with state variables. It is important to use only differentiable math-functions from jax.

  • nstatevars (int, optional) – Number of state variables (default is 0).

  • jit (bool, optional) – A flag to invoke just-in-time compilation (default is True).

  • parallel (bool, optional) – A flag to invoke parallel strain energy density function evaluations (default is False). If True, the quadrature points are executed in parallel. The number of devices must be greater or equal the number of quadrature points per cell.

  • **kwargs (dict, optional) – Optional keyword-arguments for the strain energy density function.

Notes

The strain energy density function \(\psi\) must be given in terms of the right Cauchy-Green deformation tensor \(\boldsymbol{C} = \boldsymbol{F}^T \boldsymbol{F}\).

Warning

It is important to use only differentiable math-functions from jax!

Take this minimal code-block as template

\[\psi = \psi(\boldsymbol{C})\]
import felupe as fem
import felupe.constitution.jax as mat
import jax.numpy as jnp

def neo_hooke(C, mu):
    "Strain energy function of the Neo-Hookean material formulation."
    return mu / 2 * (jnp.linalg.det(C) ** (-1/3) * jnp.trace(C) - 3)

umat = mat.Hyperelastic(neo_hooke, mu=1)

and this code-block for material formulations with state variables.

\[\psi = \psi(\boldsymbol{C}, \boldsymbol{\zeta})\]
import felupe as fem
import felupe.constitution.jax as mat
import jax.numpy as jnp

def viscoelastic(C, Cin, mu, eta, dtime):
    "Finite strain viscoelastic material formulation."

    # unimodular part of the right Cauchy-Green deformation tensor
    Cu = jnp.linalg.det(C) ** (-1 / 3) * C

    # update of state variables by evolution equation
    Ci = Cin.reshape(3, 3) + mu / eta * dtime * Cu
    Ci = jnp.linalg.det(Ci) ** (-1 / 3) * Ci

    # first invariant of elastic part of right Cauchy-Green deformation tensor
    I1 = jnp.trace(Cu @ jnp.linalg.inv(Ci))

    # strain energy function and state variable
    return mu / 2 * (I1 - 3), Ci.ravel()

umat = mat.Hyperelastic(viscoelastic, mu=1, eta=1, dtime=1, nstatevars=9)

Note

See the documentation of JAX for further details. JAX uses single-precision (32bit) data types by default. This requires to relax the tolerance of newtonrhapson() to tol=1e-4. If required, JAX may be enforced to use double-precision at startup with jax.config.update("jax_enable_x64", True).

Examples

View force-stretch curves on elementary incompressible deformations.

>>> import felupe as fem
>>> import felupe.constitution.jax as mat
>>> import jax.numpy as jnp
>>>
>>> def neo_hooke(C, mu):
...     "Strain energy function of the Neo-Hookean material formulation."
...     return mu / 2 * (jnp.linalg.det(C) ** (-1/3) * jnp.trace(C) - 3)
>>>
>>> umat = mat.Hyperelastic(neo_hooke, mu=1)
>>> ax = umat.plot(incompressible=True)
../../../_images/jax-2_00_00.png
copy()#

Return a deep-copy of the constitutive material.

gradient(x)#

Return the evaluated gradient of the strain energy density function.

Parameters:

x (list of ndarray) – The list with input arguments. These contain the extracted fields of a FieldContainer along with the old vector of state variables, [*field.extract(), statevars_old].

Returns:

A list with the evaluated gradient(s) of the strain energy density function and the updated vector of state variables.

Return type:

list of ndarray

hessian(x)#

Return the evaluated upper-triangle components of the hessian(s) of the strain energy density function.

Parameters:

x (list of ndarray) – The list with input arguments. These contain the extracted fields of a FieldContainer along with the old vector of state variables, [*field.extract(), statevars_old].

Returns:

A list with the evaluated upper-triangle components of the hessian(s) of the strain energy density function.

Return type:

list of ndarray

optimize(ux=None, ps=None, bx=None, incompressible=False, relative=False, **kwargs)#

Optimize the material parameters by a least-squares fit on experimental stretch-stress data.

Parameters:
  • ux (array of shape (2, ...) or None, optional) – Experimental uniaxial stretch and force-per-undeformed-area data (default is None).

  • ps (array of shape (2, ...) or None, optional) – Experimental planar-shear stretch and force-per-undeformed-area data (default is None).

  • bx (array of shape (2, ...) or None, optional) – Experimental biaxial stretch and force-per-undeformed-area data (default is None).

  • incompressible (bool, optional) – A flag to enforce incompressible deformations (default is False).

  • relative (bool, optional) – A flag to optimize relative instead of absolute residuals, i.e. (predicted - observed) / observed instead of predicted - observed (default is False).

  • **kwargs (dict, optional) – Optional keyword arguments are passed to scipy.optimize.least_squares().

Returns:

  • ConstitutiveMaterial – A copy of the constitutive material with the optimized material parameters.

  • scipy.optimize.OptimizeResult – Represents the optimization result.

Notes

Warning

At least one load case, i.e. one of the arguments ux, ps or bx must not be None.

The vector of residuals is given in Eq. (3) in case of absolute residuals

(1)#\[ \begin{align}\begin{aligned}\begin{split}\boldsymbol{r} &= \begin{bmatrix} \boldsymbol{r}_\text{ux} \\ \boldsymbol{r}_\text{ps} \\ \boldsymbol{r}_\text{bx} \end{bmatrix}\end{split}\\r_\text{ux}(\lambda_i) &= P_\text{ux}(\lambda_i) - P_\text{ux, observed}(\lambda_i)\\r_\text{ps}(\lambda_i) &= P_\text{ps}(\lambda_i) - P_\text{ps, observed}(\lambda_i)\\r_\text{bx}(\lambda_i) &= P_\text{bx}(\lambda_i) - P_\text{bx, observed}(\lambda_i)\end{aligned}\end{align} \]

and in Eq. (4) in case of relative residuals.

(2)#\[ \begin{align}\begin{aligned}\begin{split}\boldsymbol{r} &= \begin{bmatrix} \boldsymbol{r}_\text{ux} \\ \boldsymbol{r}_\text{ps} \\ \boldsymbol{r}_\text{bx} \end{bmatrix}\end{split}\\r_\text{ux}(\lambda_i) &= \frac{ P_\text{ux}(\lambda_i) - P_\text{ux, observed}(\lambda_i)}{ P_\text{ux, observed}(\lambda_i) }\\r_\text{ps}(\lambda_i) &= \frac{ P_\text{ps}(\lambda_i) - P_\text{ps, observed}(\lambda_i)}{ P_\text{ps, observed}(\lambda_i) }\\r_\text{bx}(\lambda_i) &= \frac{ P_\text{bx}(\lambda_i) - P_\text{bx, observed}(\lambda_i)}{ P_\text{bx, observed}(\lambda_i) }\end{aligned}\end{align} \]

Examples

The Anssari-Benam Bucchi material model formulation is best-fitted on Treloar’s uniaxial and biaxial tension data [1]_.

>>> import numpy as np
>>> import felupe as fem
>>>
>>> λ, P = np.array(
...     [
...         [1.000, 0.00],
...         [1.240, 2.30],
...         [1.585, 4.16],
...         [2.180, 6.00],
...         [3.020, 8.80],
...         [4.030, 12.5],
...         [4.760, 16.2],
...         [5.750, 23.6],
...         [6.850, 38.5],
...         [7.250, 49.6],
...         [7.600, 64.4],
...     ]
... ).T * np.array([[1.0], [0.0980665]])
>>>
>>> umat = fem.Hyperelastic(fem.anssari_benam_bucchi)
>>> umat_new, res = umat.optimize(
...     ux=[λ, P], incompressible=True, relative=True
... )
>>>
>>> ux = np.linspace(λ.min(), λ.max(), num=50)
>>> ax = umat_new.plot(incompressible=True, ux=ux, bx=None, ps=None)
>>> ax.plot(λ, P, "C0x")
../../../_images/jax-4_00_00.png

See also

scipy.optimize.least_squares

Solve a nonlinear least-squares problem with bounds on the variables.

References

plot(incompressible=False, **kwargs)#

Return a plot with normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

Parameters:
Return type:

matplotlib.axes.Axes

See also

felupe.ViewMaterial

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

felupe.ViewMaterialIncompressible

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous incompressible deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

screenshot(filename='umat.png', incompressible=False, **kwargs)#

Save a screenshot with normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

Parameters:
  • filename (str, optional) – The filename of the screenshot (default is “umat.png”).

  • incompressible (bool, optional) – A flag to enforce views on incompressible deformations (default is False).

  • **kwargs (dict, optional) – Optional keyword-arguments for ViewMaterial or ViewMaterialIncompressible.

Return type:

matplotlib.axes.Axes

See also

felupe.ViewMaterial

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

felupe.ViewMaterialIncompressible

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous incompressible deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

view(incompressible=False, **kwargs)#

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

Parameters:
Return type:

felupe.ViewMaterial or felupe.ViewMaterialIncompressible

See also

felupe.ViewMaterial

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

felupe.ViewMaterialIncompressible

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous incompressible deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

class felupe.constitution.jax.Material(fun, nstatevars=0, jit=True, parallel=False, jacobian=None, **kwargs)[source]#

A material definition with a given function for the partial derivative of the strain energy function w.r.t. the deformation gradient tensor with Automatic Differentiation provided by jax.

Parameters:
  • fun (callable) – A gradient of the strain energy density function w.r.t. the deformation gradient tensor \(\boldsymbol{F}\). Function signature must be fun = lambda F, **kwargs: P for functions without state variables and fun = lambda F, statevars, **kwargs: [P, statevars_new] for functions with state variables. It is important to use only differentiable math-functions from jax.

  • nstatevars (int, optional) – Number of state variables (default is 0).

  • jit (bool, optional) – A flag to invoke just-in-time compilation (default is True).

  • parallel (bool, optional) – A flag to invoke parallel function evaluations (default is False). If True, the quadrature points are executed in parallel. The number of devices must be greater or equal the number of quadrature points per cell.

  • jacobian (callable or None, optional) – A callable for the Jacobian. Default is None, where jax.jacobian() is used. This may be used to switch to forward-mode differentian jax.jacfwd().

  • **kwargs (dict, optional) – Optional keyword-arguments for the gradient of the strain energy density function.

Notes

The gradient of the strain energy density function \(\frac{\partial \psi}{\partial \boldsymbol{F}}\) must be given in terms of the deformation gradient tensor \(\boldsymbol{F}\).

Warning

It is important to use only differentiable math-functions from jax!

Take this code-block as template

import felupe as fem
import felupe.constitution.jax as mat
import jax.numpy as jnp

def neo_hooke(F, mu):
    "First Piola-Kirchhoff stress of the Neo-Hookean material formulation."

    C = F.T @ F
    Cu = jnp.linalg.det(C) ** (-1/3) * C
    dev = lambda C: C - jnp.trace(C) / 3 * jnp.eye(3)

    return mu * F @ dev(Cu) @ jnp.linalg.inv(C)

umat = mat.Material(neo_hooke, mu=1)

and this code-block for material formulations with state variables:

import felupe as fem
import felupe.constitution.jax as mat
import jax.numpy as jnp

def viscoelastic(F, Cin, mu, eta, dtime):
    "Finite strain viscoelastic material formulation."

    # unimodular part of the right Cauchy-Green deformation tensor
    C = F.T @ F
    Cu = jnp.linalg.det(C) ** (-1 / 3) * C

    # update of state variables by evolution equation
    from_triu = lambda C: C[jnp.array([[0, 1, 2], [1, 3, 4], [2, 4, 5]])]
    Ci = from_triu(Cin) + mu / eta * dtime * Cu
    Ci = jnp.linalg.det(Ci) ** (-1 / 3) * Ci

    # second Piola-Kirchhoff stress tensor
    dev = lambda C: C - jnp.trace(C) / 3 * jnp.eye(3)
    S = mu * dev(Cu @ jnp.linalg.inv(Ci)) @ jnp.linalg.inv(C)

    # first Piola-Kirchhoff stress tensor and state variable
    i, j = jnp.triu_indices(3)
    to_triu = lambda C: C[i, j]
    return F @ S, to_triu(Ci)

umat = mat.Material(viscoelastic, mu=1, eta=1, dtime=1, nstatevars=6)

Note

See the documentation of JAX for further details. JAX uses single-precision (32bit) data types by default. This requires to relax the tolerance of newtonrhapson() to tol=1e-4. If required, JAX may be enforced to use double-precision at startup with jax.config.update("jax_enable_x64", True).

Examples

View force-stretch curves on elementary incompressible deformations.

>>> import felupe as fem
>>> import felupe.constitution.jax as mat
>>> import jax.numpy as jnp
>>>
>>> def neo_hooke(F, mu):
...     "First Piola-Kirchhoff stress of the Neo-Hookean material formulation."
...
...     C = F.T @ F
...     Cu = jnp.linalg.det(C) ** (-1/3) * C
...     dev = lambda C: C - jnp.trace(C) / 3 * jnp.eye(3)
...
...     return mu * F @ dev(Cu) @ jnp.linalg.inv(C)
>>>
>>> umat = mat.Material(neo_hooke, mu=1)
>>> ax = umat.plot(incompressible=True)
../../../_images/jax-6_00_00.png
copy()#

Return a deep-copy of the constitutive material.

gradient(x)[source]#

Return the evaluated gradient of the strain energy density function.

Parameters:

x (list of ndarray) – The list with input arguments. These contain the extracted fields of a FieldContainer along with the old vector of state variables, [*field.extract(), statevars_old].

Returns:

A list with the evaluated gradient(s) of the strain energy density function and the updated vector of state variables.

Return type:

list of ndarray

hessian(x)[source]#

Return the evaluated upper-triangle components of the hessian(s) of the strain energy density function.

Parameters:

x (list of ndarray) – The list with input arguments. These contain the extracted fields of a FieldContainer along with the old vector of state variables, [*field.extract(), statevars_old].

Returns:

A list with the evaluated upper-triangle components of the hessian(s) of the strain energy density function.

Return type:

list of ndarray

optimize(ux=None, ps=None, bx=None, incompressible=False, relative=False, **kwargs)#

Optimize the material parameters by a least-squares fit on experimental stretch-stress data.

Parameters:
  • ux (array of shape (2, ...) or None, optional) – Experimental uniaxial stretch and force-per-undeformed-area data (default is None).

  • ps (array of shape (2, ...) or None, optional) – Experimental planar-shear stretch and force-per-undeformed-area data (default is None).

  • bx (array of shape (2, ...) or None, optional) – Experimental biaxial stretch and force-per-undeformed-area data (default is None).

  • incompressible (bool, optional) – A flag to enforce incompressible deformations (default is False).

  • relative (bool, optional) – A flag to optimize relative instead of absolute residuals, i.e. (predicted - observed) / observed instead of predicted - observed (default is False).

  • **kwargs (dict, optional) – Optional keyword arguments are passed to scipy.optimize.least_squares().

Returns:

  • ConstitutiveMaterial – A copy of the constitutive material with the optimized material parameters.

  • scipy.optimize.OptimizeResult – Represents the optimization result.

Notes

Warning

At least one load case, i.e. one of the arguments ux, ps or bx must not be None.

The vector of residuals is given in Eq. (3) in case of absolute residuals

(3)#\[ \begin{align}\begin{aligned}\begin{split}\boldsymbol{r} &= \begin{bmatrix} \boldsymbol{r}_\text{ux} \\ \boldsymbol{r}_\text{ps} \\ \boldsymbol{r}_\text{bx} \end{bmatrix}\end{split}\\r_\text{ux}(\lambda_i) &= P_\text{ux}(\lambda_i) - P_\text{ux, observed}(\lambda_i)\\r_\text{ps}(\lambda_i) &= P_\text{ps}(\lambda_i) - P_\text{ps, observed}(\lambda_i)\\r_\text{bx}(\lambda_i) &= P_\text{bx}(\lambda_i) - P_\text{bx, observed}(\lambda_i)\end{aligned}\end{align} \]

and in Eq. (4) in case of relative residuals.

(4)#\[ \begin{align}\begin{aligned}\begin{split}\boldsymbol{r} &= \begin{bmatrix} \boldsymbol{r}_\text{ux} \\ \boldsymbol{r}_\text{ps} \\ \boldsymbol{r}_\text{bx} \end{bmatrix}\end{split}\\r_\text{ux}(\lambda_i) &= \frac{ P_\text{ux}(\lambda_i) - P_\text{ux, observed}(\lambda_i)}{ P_\text{ux, observed}(\lambda_i) }\\r_\text{ps}(\lambda_i) &= \frac{ P_\text{ps}(\lambda_i) - P_\text{ps, observed}(\lambda_i)}{ P_\text{ps, observed}(\lambda_i) }\\r_\text{bx}(\lambda_i) &= \frac{ P_\text{bx}(\lambda_i) - P_\text{bx, observed}(\lambda_i)}{ P_\text{bx, observed}(\lambda_i) }\end{aligned}\end{align} \]

Examples

The Anssari-Benam Bucchi material model formulation is best-fitted on Treloar’s uniaxial and biaxial tension data [1]_.

>>> import numpy as np
>>> import felupe as fem
>>>
>>> λ, P = np.array(
...     [
...         [1.000, 0.00],
...         [1.240, 2.30],
...         [1.585, 4.16],
...         [2.180, 6.00],
...         [3.020, 8.80],
...         [4.030, 12.5],
...         [4.760, 16.2],
...         [5.750, 23.6],
...         [6.850, 38.5],
...         [7.250, 49.6],
...         [7.600, 64.4],
...     ]
... ).T * np.array([[1.0], [0.0980665]])
>>>
>>> umat = fem.Hyperelastic(fem.anssari_benam_bucchi)
>>> umat_new, res = umat.optimize(
...     ux=[λ, P], incompressible=True, relative=True
... )
>>>
>>> ux = np.linspace(λ.min(), λ.max(), num=50)
>>> ax = umat_new.plot(incompressible=True, ux=ux, bx=None, ps=None)
>>> ax.plot(λ, P, "C0x")
../../../_images/jax-8_00_00.png

See also

scipy.optimize.least_squares

Solve a nonlinear least-squares problem with bounds on the variables.

References

plot(incompressible=False, **kwargs)#

Return a plot with normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

Parameters:
Return type:

matplotlib.axes.Axes

See also

felupe.ViewMaterial

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

felupe.ViewMaterialIncompressible

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous incompressible deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

screenshot(filename='umat.png', incompressible=False, **kwargs)#

Save a screenshot with normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

Parameters:
  • filename (str, optional) – The filename of the screenshot (default is “umat.png”).

  • incompressible (bool, optional) – A flag to enforce views on incompressible deformations (default is False).

  • **kwargs (dict, optional) – Optional keyword-arguments for ViewMaterial or ViewMaterialIncompressible.

Return type:

matplotlib.axes.Axes

See also

felupe.ViewMaterial

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

felupe.ViewMaterialIncompressible

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous incompressible deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

view(incompressible=False, **kwargs)#

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

Parameters:
Return type:

felupe.ViewMaterial or felupe.ViewMaterialIncompressible

See also

felupe.ViewMaterial

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

felupe.ViewMaterialIncompressible

Create views on normal force per undeformed area vs. stretch curves for the elementary homogeneous incompressible deformations uniaxial tension/compression, planar shear and biaxial tension of a given isotropic material formulation.

felupe.constitution.jax.total_lagrange(material)[source]#

Decorate a second Piola-Kirchhoff stress Total-Lagrange material formulation as a first Piola-Kirchoff stress function.

Notes

\[\delta \psi = \boldsymbol{F} \boldsymbol{S} : \delta \boldsymbol{F}\]

Examples

>>> import felupe as fem
>>> import felupe.constitution.jax as mat
>>> import jax.numpy as jnp
>>>
>>> @mat.total_lagrange
>>> def neo_hooke_total_lagrange(F, mu=1):
>>>     C = F.T @ F
>>>     dev = lambda C: C - jnp.trace(C) / 3 * jnp.eye(3)
>>>     S = mu * dev(jnp.linalg.det(C)**(-1/3) * C) @ jnp.linalg.inv(C)
>>>     return S
>>>
>>> umat = mat.Material(neo_hooke_total_lagrange, mu=1)

See also

felupe.constitution.jax.Hyperelastic

A hyperelastic material definition with a given function for the strain energy density function per unit undeformed volume with Automatic Differentiation provided by jax.

felupe.constitution.jax.Material

A material definition with a given function for the partial derivative of the strain energy function w.r.t. the deformation gradient tensor with Automatic Differentiation provided by jax.

felupe.constitution.jax.updated_lagrange(material)[source]#

Decorate a Cauchy-stress Updated-Lagrange material formulation as a first Piola- Kirchoff stress function.

Notes

\[\delta \psi = J \boldsymbol{\sigma} \boldsymbol{F}^{-T} : \delta \boldsymbol{F}\]

Examples

>>> import felupe as fem
>>> import felupe.constitution.jax as mat
>>> import jax.numpy as jnp
>>>
>>> @fem.updated_lagrange
>>> def neo_hooke_updated_lagrange(F, mu=1):
>>>     J = jnp.linalg.det(F)
>>>     b = F @ F.T
>>>     dev = lambda b: b - jnp.trace(b) / 3 * jnp.eye(3)
>>>     τ = mu * dev(J**(-2/3) * b)
>>>     return τ / J
>>>
>>> umat = mat.Material(neo_hooke_updated_lagrange, mu=1)

See also

felupe.constitution.jax.Hyperelastic

A hyperelastic material definition with a given function for the strain energy density function per unit undeformed volume with Automatic Differentiation provided by jax.

felupe.constitution.jax.Material

A material definition with a given function for the partial derivative of the strain energy function w.r.t. the deformation gradient tensor with Automatic Differentiation provided by jax.

felupe.constitution.jax.models.hyperelastic.blatz_ko(C, mu)[source]#

Strain energy function of the Blatz-Ko isotropic hyperelastic foam material formulation [1]_.

Parameters:

Notes

The Poisson ratio of the Blatz-Ko model formulation is \(\nu = 0.25\). The strain energy function is given in Eq. (17)

(5)#\[\psi = \frac{\mu}{2} \left(\frac{I_2}{I_3} + 2 \sqrt{I_3} - 5 \right)\]

The shear modulus \(\mu\) is related to young’s modulus as denoted in Eq. (18).

(6)#\[\mu = \frac{2 E}{5}\]

Examples

First, choose the desired automatic differentiation backend

>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the hyperelastic material.

>>> import felupe as fem
>>>
>>> umat = mat.Hyperelastic(mat.models.hyperelastic.blatz_ko, mu=1.0)
>>> ax = umat.plot()
../../../_images/jax-11_00_00.png

References

felupe.constitution.jax.models.hyperelastic.extended_tube(C, Gc, delta, Ge, beta)[source]#

Strain energy function of the isotropic hyperelastic Extended Tube [1]_ material formulation.

Parameters:
  • C (tensortrax.Tensor or jax.Array) – Right Cauchy-Green deformation tensor.

  • Gc (float) – Cross-link contribution to the initial shear modulus.

  • delta (float) – Finite extension parameter of the polymer strands.

  • Ge (float) – Constraint contribution to the initial shear modulus.

  • beta (float) – Global rearrangements of cross-links upon deformation (release of topological constraints).

Notes

The strain energy function is given in Eq. (19)

(7)#\[\psi = \frac{G_c}{2} \left[ \frac{\left( 1 - \delta^2 \right) \left( \hat{I}_1 - 3 \right)}{1 - \delta^2 \left( \hat{I}_1 - 3 \right)} + \ln \left( 1 - \delta^2 \left( \hat{I}_1 - 3 \right) \right) \right] + \frac{2 G_e}{\beta^2} \left( \hat{\lambda}_1^{-\beta} + \hat{\lambda}_2^{-\beta} + \hat{\lambda}_3^{-\beta} - 3 \right)\]

with the first main invariant of the distortional part of the right Cauchy-Green deformation tensor as given in Eq. (20)

(8)#\[\hat{I}_1 = J^{-2/3} \text{tr}\left( \boldsymbol{C} \right)\]

and the principal stretches, obtained from the distortional part of the right Cauchy-Green deformation tensor, see Eq. (21).

(9)#\[ \begin{align}\begin{aligned}\lambda^2_\alpha &= \text{eigvals}\left( \boldsymbol{C} \right)\\\hat{\lambda}_\alpha &= J^{-1/3} \lambda_\alpha\end{aligned}\end{align} \]

The initial shear modulus results from the sum of the cross-link and the constraint contributions to the total initial shear modulus as denoted in Eq. (22).

(10)#\[\mu = G_e + G_c\]

Examples

First, choose the desired automatic differentiation backend

>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the hyperelastic material.

>>> import felupe as fem
>>>
>>> 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)
../../../_images/jax-14_00_00.png

References

felupe.constitution.jax.models.hyperelastic.miehe_goektepe_lulei(C, mu, N, U, p, q)[source]#

Strain energy function of the isotropic hyperelastic micro-sphere model formulation [1]_.

Parameters:
  • C (tensortrax.Tensor or jax.Array) – Right Cauchy-Green deformation tensor.

  • mu (float) – Shear modulus (ground state stifness).

  • N (float) – Number of chain segments (chain locking response).

  • U (float) – Tube geometry parameter (3D locking characteristics).

  • p (float) – Non-affine stretch parameter (additional constraint stifness).

  • q (float) – Non-affine tube parameter (shape of constraint stress).

Examples

First, choose the desired automatic differentiation backend

>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the hyperelastic material.

>>> import felupe as fem
>>>
>>> umat = mat.Hyperelastic(
...     mat.models.hyperelastic.miehe_goektepe_lulei,
...     mu=0.1475,
...     N=3.273,
...     p=9.31,
...     U=9.94,
...     q=0.567,
... )
>>> ux = ps = fem.math.linsteps([1, 2], num=50)
>>> bx = fem.math.linsteps([1, 1.5], num=50)
>>> ax = umat.plot(ux=ux, ps=ps, bx=bx, incompressible=True)
../../../_images/jax-17_00_00.png

References

felupe.constitution.jax.models.hyperelastic.mooney_rivlin(C, C10, C01)[source]#

Strain energy function of the isotropic hyperelastic Mooney-Rivlin material formulation.

Parameters:
  • C (tensortrax.Tensor or jax.Array) – Right Cauchy-Green deformation tensor.

  • C10 (float) – First material parameter associated to the first invariant.

  • C01 (float) – Second material parameter associated to the second invariant.

Notes

The strain energy function is given in Eq. (29)

(11)#\[\psi = C_{10} \left(\hat{I}_1 - 3 \right) + C_{01} \left(\hat{I}_2 - 3 \right)\]

with the first and second main invariant of the distortional part of the right Cauchy-Green deformation tensor, see Eq. (30).

(12)#\[ \begin{align}\begin{aligned}\hat{I}_1 &= J^{-2/3} \text{tr}\left( \boldsymbol{C} \right)\\\hat{I}_2 &= J^{-4/3} \frac{1}{2} \left( \text{tr}\left(\boldsymbol{C}\right)^2 - \text{tr}\left(\boldsymbol{C}^2\right) \right)\end{aligned}\end{align} \]

The doubled sum of both material parameters is equal to the shear modulus \(\mu\) as denoted in Eq. (31).

(13)#\[\mu = 2 \left( C_{10} + C_{01} \right)\]

Examples

First, choose the desired automatic differentiation backend

>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the hyperelastic material.

>>> umat = mat.Hyperelastic(
...     mat.models.hyperelastic.mooney_rivlin, C10=0.3, C01=0.8
... )
>>> ax = umat.plot(incompressible=True)
../../../_images/jax-20_00_00.png
felupe.constitution.jax.models.hyperelastic.neo_hooke(C, mu)[source]#

Strain energy function of the isotropic hyperelastic Neo-Hookean material formulation.

Parameters:

Notes

The strain energy function is given in Eq. (32).

(14)#\[\psi = \frac{\mu}{2} \left(\text{tr}\left(\hat{\boldsymbol{C}}\right) - 3\right)\]

Examples

First, choose the desired automatic differentiation backend

>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the hyperelastic material.

>>> import felupe as fem
>>>
>>> umat = mat.Hyperelastic(mat.models.hyperelastic.neo_hooke, mu=1.0)
>>> ax = umat.plot(incompressible=True)
../../../_images/jax-23_00_00.png
felupe.constitution.jax.models.hyperelastic.storakers(C, mu, alpha, beta)[source]#

Strain energy function of the Storåkers isotropic hyperelastic foam material formulation [1]_.

Parameters:

Notes

The strain energy function is given in Eq. (39)

(15)#\[\psi = \sum_i \frac{2 \mu_i}{\alpha^2_i} \left[ \lambda_1^{\alpha_i} + \lambda_2^{\alpha_i} + \lambda_3^{\alpha_i} - 3 + \frac{1}{\beta_i} \left( J^{-\alpha_i \beta_i} - 1 \right) \right]\]

The sum of the moduli \(\mu_i\) is equal to the initial shear modulus \(\mu\), see Eq. (40),

(16)#\[\mu = \sum_i \mu_i\]

and the initial bulk modulus is given in Eq. (41).

(17)#\[K = \sum_i 2 \mu_i \left( \frac{1}{3} + \beta_i \right)\]

Examples

First, choose the desired automatic differentiation backend

>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the hyperelastic material.

>>> import felupe as fem
>>>
>>> umat = mat.Hyperelastic(
...     mat.models.hyperelastic.storakers,
...     mu=[4.5 * (1.85 / 2), -4.5 * (-9.2 / 2)],
...     alpha=[1.85, -9.2],
...     beta=[0.92, 0.92],
... )
>>> ax = umat.plot(
...     ux=fem.math.linsteps([1, 2], 15),
...     ps=fem.math.linsteps([1, 1], 15),
...     bx=fem.math.linsteps([1, 1], 9),
... )
../../../_images/jax-26_00_00.png

References

felupe.constitution.jax.models.hyperelastic.third_order_deformation(C, C10, C01, C11, C20, C30)[source]#

Strain energy function of the isotropic hyperelastic Third-Order-Deformation material formulation.

Parameters:
  • C (tensortrax.Tensor or jax.Array) – Right Cauchy-Green deformation tensor.

  • C10 (float) – Material parameter associated to the linear term of the first invariant.

  • C01 (float) – Material parameter associated to the linear term of the second invariant.

  • C11 (float) – Material parameter associated to the mixed term of the first and second invariant.

  • C20 (float) – Material parameter associated to the quadratic term of the first invariant.

  • C30 (float) – Material parameter associated to the cubic term of the first invariant.

Notes

The strain energy function is given in Eq. (42)

(18)#\[ \begin{align}\begin{aligned}\psi &= C_{10} \left(\hat{I}_1 - 3 \right) + C_{01} \left(\hat{I}_2 - 3 \right) + C_{11} \left(\hat{I}_1 - 3 \right) \left(\hat{I}_2 - 3 \right)\\ &+ C_{20} \left(\hat{I}_1 - 3 \right)^2 + C_{30} \left(\hat{I}_1 - 3 \right)^3\end{aligned}\end{align} \]

with the first and second main invariant of the distortional part of the right Cauchy-Green deformation tensor, see Eq. (43).

(19)#\[ \begin{align}\begin{aligned}\hat{I}_1 &= J^{-2/3} \text{tr}\left( \boldsymbol{C} \right)\\\hat{I}_2 &= J^{-4/3} \frac{1}{2} \left( \text{tr}\left(\boldsymbol{C}\right)^2 - \text{tr}\left(\boldsymbol{C}^2\right) \right)\end{aligned}\end{align} \]

The doubled sum of the material parameters \(C_{10}\) and \(C_{01}\) is equal to the initial shear modulus \(\mu\) as denoted in Eq. (44).

(20)#\[\mu = 2 \left( C_{10} + C_{01} \right)\]

Examples

First, choose the desired automatic differentiation backend

>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the hyperelastic material.

>>> umat = mat.Hyperelastic(
...     mat.models.hyperelastic.third_order_deformation,
...     C10=0.5,
...     C01=0.1,
...     C11=0.01,
...     C20=-0.1,
...     C30=0.02,
... )
>>> ax = umat.plot(incompressible=True)
../../../_images/jax-29_00_00.png
felupe.constitution.jax.models.hyperelastic.van_der_waals(C, mu, limit, a, beta)[source]#

Strain energy function of the Van der Waals [1]_ material formulation.

Parameters:
  • C (tensortrax.Tensor or jax.Array) – Right Cauchy-Green deformation tensor.

  • mu (float) – Initial shear modulus.

  • limit (float) – Limiting stretch \(\lambda_m\) at which the polymer chain network becomes locked.

  • a (float) – Attractive interactions between the quasi-particles.

  • beta (float) – Mixed-Invariant factor: 0 for pure I1- and 1 for pure I2-contribution.

Examples

First, choose the desired automatic differentiation backend

>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the hyperelastic material.

>>> import felupe as fem
>>>
>>> umat = mat.Hyperelastic(
...     mat.models.hyperelastic.van_der_waals,
...     mu=1.0,
...     beta=0.1,
...     a=0.5,
...     limit=5.0
... )
>>> ax = umat.plot(incompressible=True)
../../../_images/jax-32_00_00.png

References

felupe.constitution.jax.models.hyperelastic.yeoh(C, C10, C20, C30)[source]#

Strain energy function of the isotropic hyperelastic Yeoh material formulation.

Parameters:
  • C (tensortrax.Tensor or jax.Array) – Right Cauchy-Green deformation tensor.

  • C10 (float) – Material parameter associated to the linear term of the first invariant.

  • C20 (float) – Material parameter associated to the quadratic term of the first invariant.

  • C30 (float) – Material parameter associated to the cubic term of the first invariant.

Notes

The strain energy function is given in Eq. (45)

(21)#\[\psi = C_{10} \left(\hat{I}_1 - 3 \right) + C_{20} \left(\hat{I}_1 - 3 \right)^2 + C_{30} \left(\hat{I}_1 - 3 \right)^3\]

with the first main invariant of the distortional part of the right Cauchy-Green deformation tensor, see Eq. (46).

(22)#\[\hat{I}_1 = J^{-2/3} \text{tr}\left( \boldsymbol{C} \right)\]

The \(C_{10}\) material parameter is equal to half the initial shear modulus \(\mu\) as denoted in Eq. (47).

(23)#\[\mu = 2 C_{10}\]

Examples

First, choose the desired automatic differentiation backend

>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the hyperelastic material.

>>> umat = mat.Hyperelastic(
...     mat.models.hyperelastic.yeoh, C10=0.5, C20=-0.1, C30=0.02
... )
>>> ax = umat.plot(incompressible=True)
../../../_images/jax-35_00_00.png
felupe.constitution.jax.models.lagrange.morph(F, statevars, p)[source]#

Second Piola-Kirchhoff stress tensor of the MORPH model formulation [1]_.

Parameters:
  • F (tensortrax.Tensor or jax.Array) – Deformation gradient 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. (48).

Warning

While the MORPH-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.

(24)#\[ \begin{align}\begin{aligned}\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)\end{aligned}\end{align} \]

A sigmoid-function is used inside the deformation-dependent variables \(\alpha\), \(\beta\) and \(\gamma\), see Eq. (49).

(25)#\[ \begin{align}\begin{aligned}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)\end{aligned}\end{align} \]

The rate of deformation is described by the Lagrangian tensor and its Tresca- invariant, see Eq. (50).

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.

(26)#\[ \begin{align}\begin{aligned}\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\end{aligned}\end{align} \]

The additional stresses evolve between the limiting stresses, see Eq. (51). The additional deviatoric-enforcement terms [1]_ are neglected in this implementation.

(27)#\[ \begin{align}\begin{aligned}\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}\end{aligned}\end{align} \]

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 tm.special.from_triu_1d() and export them as tm.special.triu_1d().

Examples

First, choose the desired automatic differentiation backend

>>> import felupe as fem
>>>
>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the material.

>>> umat = mat.Material(
...     mat.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,
... )
../../../_images/jax-38_00_00.png

References

See also

felupe.constitution.tensortrax.models.lagrange.morph_representative_directions

Strain energy function of the MORPH model formulation, implemented by the concept of representative directions.

felupe.constitution.jax.models.lagrange.morph_representative_directions

Strain energy function of the MORPH model formulation, implemented by the concept of representative directions.

felupe.constitution.jax.models.lagrange.morph_representative_directions(F, statevars, p, ε=1e-06)[source]#

First Piola-Kirchhoff stress tensor of the MORPH model formulation [1]_, implemented by the concept of representative directions [2], [3].

Parameters:
  • F (tensortrax.Tensor or jax.Array) – Deformation gradient 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-6).

Examples

First, choose the desired automatic differentiation backend

>>> import felupe as fem
>>>
>>> # import felupe.constitution.jax as mat
>>> import felupe.constitution.tensortrax as mat

and create the material.

>>> umat = mat.Material(
...     mat.models.lagrange.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],
...        [1, 2.75, 1, 2.75],
...        num=20,
...    ),
...    ps=None,
...    bx=None,
... )
../../../_images/jax-41_00_00.png

References

See also

felupe.constitution.tensortrax.models.lagrange.morph

Strain energy function of the MORPH model formulation.

felupe.constitution.jax.models.lagrange.morph

Strain energy function of the MORPH model formulation.

felupe.constitution.jax.vmap(fun, in_axes=0, out_axes=0, method=<function vmap>, **kwargs)[source]#

Vectorizing map. Creates a function which maps fun over argument axes. This decorator treats all non-specified arguments and keyword-arguments as static.

See also

jax.vmap

Vectorizing map. Creates a function which maps fun over argument axes.