From e6135d0eafb677ee49d3623e375dbf2f77fa80ed Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 13 Nov 2024 16:07:46 +0100 Subject: [PATCH] Add math-function for the right Cauchy-Green deformation tensor (#898) * Add `math.right_cauchy_green_deformation()` * format black * Update CHANGELOG.md --- CHANGELOG.md | 2 ++ src/felupe/constitution/jax/_helpers.py | 1 - src/felupe/constitution/jax/_hyperelastic.py | 1 - src/felupe/constitution/jax/_material.py | 1 - .../jax/models/hyperelastic/_mooney_rivlin.py | 1 - .../hyperelastic/_third_order_deformation.py | 1 - .../jax/models/hyperelastic/_yeoh.py | 1 - .../jax/models/lagrange/_morph.py | 1 - .../jax/models/lagrange/_morph_uniaxial.py | 1 - src/felupe/field/_evaluate.py | 17 +++++++++ src/felupe/math/__init__.py | 2 ++ src/felupe/math/_field.py | 24 ++++++++++--- src/felupe/region/_boundary.py | 35 +++++++++++++++++++ tests/test_constitution_jax.py | 4 --- tests/test_field.py | 1 + 15 files changed, 76 insertions(+), 17 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8a225b7e..342e4d94 100644 --- a/CHANGELOG.md +++ b/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`. diff --git a/src/felupe/constitution/jax/_helpers.py b/src/felupe/constitution/jax/_helpers.py index bd362070..fba57c7c 100644 --- a/src/felupe/constitution/jax/_helpers.py +++ b/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) diff --git a/src/felupe/constitution/jax/_hyperelastic.py b/src/felupe/constitution/jax/_hyperelastic.py index f18d83f9..8460d8cc 100644 --- a/src/felupe/constitution/jax/_hyperelastic.py +++ b/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) diff --git a/src/felupe/constitution/jax/_material.py b/src/felupe/constitution/jax/_material.py index 44317b3e..722979db 100644 --- a/src/felupe/constitution/jax/_material.py +++ b/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 diff --git a/src/felupe/constitution/jax/models/hyperelastic/_mooney_rivlin.py b/src/felupe/constitution/jax/models/hyperelastic/_mooney_rivlin.py index fd04eddc..9f28837c 100644 --- a/src/felupe/constitution/jax/models/hyperelastic/_mooney_rivlin.py +++ b/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 diff --git a/src/felupe/constitution/jax/models/hyperelastic/_third_order_deformation.py b/src/felupe/constitution/jax/models/hyperelastic/_third_order_deformation.py index 429ce33f..b1fe2ce7 100644 --- a/src/felupe/constitution/jax/models/hyperelastic/_third_order_deformation.py +++ b/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 diff --git a/src/felupe/constitution/jax/models/hyperelastic/_yeoh.py b/src/felupe/constitution/jax/models/hyperelastic/_yeoh.py index 6ac3236a..010bec7f 100644 --- a/src/felupe/constitution/jax/models/hyperelastic/_yeoh.py +++ b/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 diff --git a/src/felupe/constitution/jax/models/lagrange/_morph.py b/src/felupe/constitution/jax/models/lagrange/_morph.py index 47c1b213..cf7d8c3c 100644 --- a/src/felupe/constitution/jax/models/lagrange/_morph.py +++ b/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 diff --git a/src/felupe/constitution/jax/models/lagrange/_morph_uniaxial.py b/src/felupe/constitution/jax/models/lagrange/_morph_uniaxial.py index b168539c..dd176f39 100644 --- a/src/felupe/constitution/jax/models/lagrange/_morph_uniaxial.py +++ b/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] diff --git a/src/felupe/field/_evaluate.py b/src/felupe/field/_evaluate.py index 3250e918..7aa5c482 100644 --- a/src/felupe/field/_evaluate.py +++ b/src/felupe/field/_evaluate.py @@ -16,6 +16,7 @@ along with FElupe. If not, see . """ +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. diff --git a/src/felupe/math/__init__.py b/src/felupe/math/__init__.py index bdc3599e..2eefa61a 100644 --- a/src/felupe/math/__init__.py +++ b/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", diff --git a/src/felupe/math/_field.py b/src/felupe/math/_field.py index dd9fb8d3..ab05ca00 100644 --- a/src/felupe/math/_field.py +++ b/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) diff --git a/src/felupe/region/_boundary.py b/src/felupe/region/_boundary.py index 39622e72..15076b4c 100644 --- a/src/felupe/region/_boundary.py +++ b/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__( diff --git a/tests/test_constitution_jax.py b/tests/test_constitution_jax.py index 32c066a7..8cf6dc06 100644 --- a/tests/test_constitution_jax.py +++ b/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, diff --git a/tests/test_field.py b/tests/test_field.py index 3af14b1e..d6a27d86 100644 --- a/tests/test_field.py +++ b/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)