Add math-function for the right Cauchy-Green deformation tensor (#898)

* Add `math.right_cauchy_green_deformation()`

* format black

* Update CHANGELOG.md

CHANGELOG.md

@@ -15,6 +15,8 @@ All notable changes to this project will be documented in this file. The format
 - Add material models for JAX-based materials: `felupe.constitution.jax.models.hyperelastic.miehe_goektepe_lulei()`, `felupe.constitution.jax.models.hyperelastic.mooney_rivlin()`, `felupe.constitution.jax.models.hyperelastic.yeoh()`, `felupe.constitution.jax.models.hyperelastic.third_order_deformation()`, `felupe.constitution.jax.models.lagrange.morph()`, `felupe.constitution.jax.models.lagrange.morph_representative_directions()`.
 - Add `felupe.constitution.jax.total_lagrange()`, `felupe.constitution.jax.updated_lagrange()` and `felupe.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.
 
 ### Changed
 - Change default `np.einsum(..., order="K")` to `np.einsum(..., order="C")` in the methods of `Field`, `FieldAxisymmetric`, `FieldPlaneStrain` and `FieldContainer`.

src/felupe/constitution/jax/_helpers.py

@@ -95,7 +95,6 @@ def vmap2(fun, in_axes=[0, 0], out_axes=[0, 0], methods=[jax.vmap, jax.vmap], **
 
 
 def as_total_lagrange(fun):
-
     @wraps(fun)
     def evaluate(F, *args, **kwargs):
         i, j = jnp.triu_indices(3)

src/felupe/constitution/jax/_hyperelastic.py

@@ -146,7 +146,6 @@ def viscoelastic(C, Cin, mu, eta, dtime):
     """
 
     def __init__(self, fun, nstatevars=0, jit=True, parallel=False, **kwargs):
-
         has_aux = nstatevars > 0
         self.fun = as_total_lagrange(fun)
 

src/felupe/constitution/jax/_material.py

@@ -160,7 +160,6 @@ def viscoelastic(F, Cin, mu, eta, dtime):
     def __init__(
         self, fun, nstatevars=0, jit=True, parallel=False, jacobian=None, **kwargs
     ):
-
         has_aux = nstatevars > 0
         self.fun = fun
 

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

@@ -25,7 +25,6 @@
 
 @wraps(mooney_rivlin_docstring)
 def mooney_rivlin(C, C10, C01):
-
     J3 = det(C) ** (-1 / 3)
     I1 = J3 * trace(C)
     I2 = (I1**2 - J3**2 * trace(C @ C)) / 2

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

@@ -25,7 +25,6 @@
 
 @wraps(tod_docstring)
 def third_order_deformation(C, C10, C01, C11, C20, C30):
-
     J3 = det(C) ** (-1 / 3)
     I1 = J3 * trace(C)
     I2 = (I1**2 - J3**2 * trace(C @ C)) / 2

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

@@ -25,6 +25,5 @@
 
 @wraps(yeoh_docstring)
 def yeoh(C, C10, C20, C30):
-
     I1 = det(C) ** (-1 / 3) * trace(C)
     return C10 * (I1 - 3) + C20 * (I1 - 3) ** 2 + C30 * (I1 - 3) ** 3

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

@@ -28,7 +28,6 @@
 @wraps(morph_docstring)
 @total_lagrange
 def morph(F, statevars, p):
-
     # right Cauchy-Green deformation tensor
     C = F.T @ F
 

src/felupe/constitution/jax/models/lagrange/_morph_uniaxial.py

@@ -25,7 +25,6 @@
 
 @wraps(morph_ux)
 def morph_uniaxial(λ, statevars, p, ε=1e-6):
-
     CTSn = statevars[:21]
     λn = 1 + statevars[21:42]
     SA1n = statevars[42:63]

src/felupe/field/_evaluate.py

@@ -16,6 +16,7 @@
 along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
 """
 
+from ..math import right_cauchy_green_deformation as right_cauchy_green
 from ..math import strain, strain_stretch_1d
 
 
@@ -71,6 +72,22 @@ def deformation_gradient(self):
         """
         return self.field[0].extract()
 
+    def right_cauchy_green_deformation(self):
+        r"""Return the right Cauchy-Green deformation tensor.
+
+        .. math::
+           :label:right-cauchy-green-deformation-tensor
+
+           \boldsymbol{F} &= \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{X}}
+
+           \boldsymbol{C} &= \boldsymbol{F}^T \boldsymbol{F}
+
+           \boldsymbol{C} &= \sum_\alpha \lambda^2_\alpha
+               \ \boldsymbol{N}_\alpha \otimes \boldsymbol{N}_\alpha
+
+        """
+        return right_cauchy_green(self.field)
+
     def strain(self, fun=strain_stretch_1d, tensor=True, asvoigt=False, n=0, **kwargs):
         r"""Return the Lagrangian strain tensor or its principal values.
 

src/felupe/math/__init__.py

@@ -8,6 +8,7 @@
     hess,
     interpolate,
     norm,
+    right_cauchy_green_deformation,
     strain,
     strain_stretch_1d,
     values,
@@ -80,6 +81,7 @@
     "majortranspose",
     "ravel",
     "reshape",
+    "right_cauchy_green_deformation",
     "sym",
     "tovoigt",
     "trace",

src/felupe/math/_field.py

@@ -35,6 +35,12 @@ def deformation_gradient(field, n=0):
     return field[n].extract(grad=True, sym=False, add_identity=True)
 
 
+def right_cauchy_green_deformation(field, n=0):
+    "Return the right Cauchy-Green deformation tensor of the n-th field."
+    F = deformation_gradient(field, n=n)
+    return dot(transpose(F), F)
+
+
 def strain_stretch_1d(stretch, k=0):
     r"""Compute the Seth-Hill strains.
 
@@ -66,13 +72,20 @@ def strain_stretch_1d(stretch, k=0):
     return strain
 
 
-def strain(field, fun=strain_stretch_1d, tensor=True, asvoigt=False, n=0, **kwargs):
+def strain(
+    field, C=None, fun=strain_stretch_1d, tensor=True, asvoigt=False, n=0, **kwargs
+):
     r"""Return Lagrangian strain tensor or its principal values of the n-th field.
 
     Parameters
     ----------
-    field : FieldContainer
-        A field container with a displacement field.
+    field : FieldContainer or None
+        A field container with a displacement field from which the deformation gradient
+        tensor is extracted if ``def_grad`` is None.
+    C : ndarray of shape (N, N, ...) or None, optional
+        Optional array of right Cauchy-Green deformation tensors. If None, the
+        right Cauchy-Green deformation tensor is obtained from the field. Default is
+        None.
     fun : callable, optional
         A callable for the one-dimensional strain-stretch relation. Function signature
         must be ``lambda stretch, **kwargs: strain`` (default is the log. strain,
@@ -132,8 +145,9 @@ def strain(field, fun=strain_stretch_1d, tensor=True, asvoigt=False, n=0, **kwar
     math.strain_stretch_1d : Compute the Seth-Hill strains.
     """
 
-    F = deformation_gradient(field)
-    C = dot(transpose(F), F)
+    if C is None:
+        F = deformation_gradient(field)
+        C = dot(transpose(F), F)
 
     if tensor:
         w, N = eigh(C)

src/felupe/region/_boundary.py

@@ -465,6 +465,41 @@ class RegionBoundary(Region):
         >>> plotter.add_legend()
         >>> mesh.plot(plotter=plotter, style="wireframe").show()
 
+    This examples shows how to evaluate the strain on the faces of the boundary region.
+    A field can't be used to evaluate the strain of a quad-cell in 3D space. However,
+    it is possible to evaluate the in-plane components of the right Cauchy-Green
+    deformation tensor on the boundary. With this tensor at hand, it is possible to
+    evaluate any in-plane strain tensor.
+
+    ..  pyvista-plot::
+        :force_static:
+
+        >>> import felupe as fem
+        >>>
+        >>> mesh = fem.Cube(n=3)
+        >>> region = fem.RegionHexahedron(mesh)
+        >>> field = fem.FieldContainer([fem.Field(region, dim=3)])
+        >>>
+        >>> field[0].values[mesh.x == 1, 0] += 0.5
+        >>> field[0].values[mesh.x == 1, 2] += 0.5
+        >>> field[0].values += mesh.rotate(-30, axis=1).points - mesh.points
+        >>>
+        >>> face = fem.RegionHexahedronBoundary(mesh, mask=mesh.x == mesh.x.max())
+        >>> u = fem.FieldContainer([fem.Field(face, dim=3)])
+        >>> u.link(field)
+        >>>
+        >>> mesh_faces = face.mesh_faces()
+        >>> mesh_faces.points += field[0].values
+        >>>
+        >>> C = fem.math.inplane(u.evaluate.right_cauchy_green_deformation(), face.tangents)
+        >>> e = fem.math.strain(u, C=C, k=0)
+        >>>
+        >>> view = mesh_faces.view(
+        ...     cell_data={"ep": fem.math.eigvalsh(e).max(axis=0).mean(axis=-2).T}
+        ... )
+        >>> view.plot(
+        ...     "ep", label="Strain (Max. Principal)", plotter=field.plot(style="wireframe")
+        ... ).show()
     """
 
     def __init__(

tests/test_constitution_jax.py

@@ -78,7 +78,6 @@ def test_hyperelastic_jax():
 
 
 def test_hyperelastic_jax_statevars():
-
     def W(C, statevars, C10, K):
         I3 = jnp.linalg.det(C)
         J = jnp.sqrt(I3)
@@ -99,7 +98,6 @@ def W(C, statevars, C10, K):
 
 
 def test_material_jax():
-
     def dWdF(F, C10, K):
         J = jnp.linalg.det(F)
         C = F.T @ F
@@ -124,7 +122,6 @@ def dWdF(F, C10, K):
 
 
 def test_material_jax_statevars():
-
     def dWdF(F, statevars, C10, K):
         J = jnp.linalg.det(F)
         C = F.T @ F
@@ -149,7 +146,6 @@ def dWdF(F, statevars, C10, K):
 
 
 def test_material_included_jax_statevars():
-
     for fun, nstatevars in zip(
         [
             mat.models.lagrange.morph,

tests/test_field.py

@@ -163,6 +163,7 @@ def test_3d_mixed():
     f.evaluate.strain()
     f.evaluate.log_strain()
     f.evaluate.green_lagrange_strain()
+    f.evaluate.right_cauchy_green_deformation()
 
     f.extract()
     f.extract(grad=False)