Skip to content

Commit

Permalink
Merge pull request #667 from adtzlr/use-tensortrax-official-api
Browse files Browse the repository at this point in the history 
Use only the official API of `tensortrax.math`
  • Loading branch information
@adtzlr
adtzlr authored Feb 28, 2024
2 parents 42d69b1 + b5976c2 commit 33bacdd
Show file tree
Hide file tree
Showing 4 changed files with 73 additions and 44 deletions.
2 changes: 2 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,13 @@ All notable changes to this project will be documented in this file. The format
- Add `mesh.stack(meshes)` as method to `MeshContainer.stack()`. Note that this only supports mesh containers with meshes of same cell-types.
- Add `NeoHooke.gradient(out=None)` and `NeoHooke.hessian(out=None)` for a location to store the results. Also for `NeoHookeCompressible`.
- Add `out`-keyword to `gradient()` and `hessian` of `NearlyIncompressible` and `ThreeFieldVariation`.
- Add optional initial state variables in `ViewMaterial(statevars=None)` and `ViewMaterialIncompressible(statevars=None)`.

### Changed
- Rename `Mesh.save()` to `Mesh.write()` and add `Mesh.save()` as an alias to `Mesh.write()`.
- Enhance the performance of `NeoHooke`, `NeoHookeCompressible`, `SolidBody` and `SolidBodyNearlyIncompressible`.
- Enhance the performance of `math.inv(out=None)` and `math.det(out=None)`.
- Use only the offical API of `tensortrax`. A workaround is used to ensure compatibility with `tensortrax` <= v0.17.1.

# Fixed
- Fix missing support for third-order- and second-order tensor combinations to `math.dot(A, B, mode=(2,3))` and `math.ddot(A, B, mode=(2,3))`.
Expand Down
11 changes: 8 additions & 3 deletions src/felupe/constitution/_models_hyperelasticity_ad.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,11 +19,16 @@
from functools import wraps

import numpy as np
from tensortrax.math import log, sqrt
from tensortrax.math import linalg, log, special, sqrt
from tensortrax.math import sum as sum1
from tensortrax.math import trace
from tensortrax.math._linalg import det, eigvalsh, inv
from tensortrax.math._special import from_triu_1d, triu_1d

det = linalg.det
inv = linalg.inv
eigvalsh = linalg.eigvalsh

from_triu_1d = special.from_triu_1d
triu_1d = special.triu_1d


def isochoric_volumetric_split(fun):
Expand Down
83 changes: 51 additions & 32 deletions src/felupe/constitution/_view.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,9 @@ class ViewMaterial(PlotMaterial):
bx : ndarray, optional
Array with stretches for equi-biaxial tension. Default is
``linsteps([1.0, 1.75], num=15)```.
statevars : ndarray or None, optional
Array with state variables (default is None). If None, the state variables are
assumed to be initially zero.
Examples
--------
Expand All @@ -146,12 +149,14 @@ def __init__(
ux=linsteps([0.7, 2.5], num=36),
ps=linsteps([1, 2.5], num=30),
bx=linsteps([1, 1.75], num=15),
statevars=None,
):
self.umat = umat
self.ux = ux
self.ps = ps
self.bx = bx
self.statevars_included = self.umat.x[-1].size > 0
self.statevars_included = (self.umat.x[-1].size > 0,)
self.statevars = statevars

def uniaxial(self, stretches=None):
"""Normal force per undeformed area vs. stretch curve for a uniaxial
Expand Down Expand Up @@ -180,11 +185,12 @@ def fun(λ3):
λ2 = λ3
F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)
if self.statevars_included:
statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
if self.statevars is None:
self.statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
P = np.zeros_like(F)
for increment, defgrad in enumerate(F.T):
P[..., [increment]], statevars = self.umat.gradient(
[F[..., [increment]], statevars]
[F[..., [increment]], self.statevars]
)
else:
P, statevars = self.umat.gradient([F, None])
Expand All @@ -196,14 +202,15 @@ def fun(λ3):
F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)

if self.statevars_included:
statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
if self.statevars is None:
self.statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
P = np.zeros_like(F)
for increment, defgrad in enumerate(F.T):
P[..., [increment]], statevars = self.umat.gradient(
[F[..., [increment]], statevars]
P[..., [increment]], self.statevars = self.umat.gradient(
[F[..., [increment]], self.statevars]
)
else:
P, statevars = self.umat.gradient([F, None])
P, self.statevars = self.umat.gradient([F, None])

valid = det(F) > np.sqrt(np.finfo(float).eps)
if not np.all(valid):
Expand Down Expand Up @@ -239,11 +246,12 @@ def planar(self, stretches=None):
def fun(λ3):
F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)
if self.statevars_included:
statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
if self.statevars is None:
self.statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
P = np.zeros_like(F)
for increment, defgrad in enumerate(F.T):
P[..., [increment]], statevars = self.umat.gradient(
[F[..., [increment]], statevars]
[F[..., [increment]], self.statevars]
)
else:
P, statevars = self.umat.gradient([F, None])
Expand All @@ -255,14 +263,15 @@ def fun(λ3):
F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)

if self.statevars_included:
statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
if self.statevars is None:
self.statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
P = np.zeros_like(F)
for increment, defgrad in enumerate(F.T):
P[..., [increment]], statevars = self.umat.gradient(
[F[..., [increment]], statevars]
P[..., [increment]], self.statevars = self.umat.gradient(
[F[..., [increment]], self.statevars]
)
else:
P, statevars = self.umat.gradient([F, None])
P, self.statevars = self.umat.gradient([F, None])

valid = det(F) > np.sqrt(np.finfo(float).eps)
if not np.all(valid):
Expand Down Expand Up @@ -297,14 +306,15 @@ def biaxial(self, stretches=None):
def fun(λ3):
F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)
if self.statevars_included:
statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
if self.statevars is None:
self.statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
P = np.zeros_like(F)
for increment, defgrad in enumerate(F.T):
P[..., [increment]], statevars = self.umat.gradient(
[F[..., [increment]], statevars]
[F[..., [increment]], self.statevars]
)
else:
P, statevars = self.umat.gradient([F, None])
P, self.statevars = self.umat.gradient([F, None])
return P[2, 2].ravel()

from scipy.optimize import root
Expand All @@ -313,14 +323,15 @@ def fun(λ3):
F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)

if self.statevars_included:
statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
if self.statevars is None:
self.statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
P = np.zeros_like(F)
for increment, defgrad in enumerate(F.T):
P[..., [increment]], statevars = self.umat.gradient(
[F[..., [increment]], statevars]
P[..., [increment]], self.statevars = self.umat.gradient(
[F[..., [increment]], self.statevars]
)
else:
P, statevars = self.umat.gradient([F, None])
P, self.statevars = self.umat.gradient([F, None])

valid = det(F) > np.sqrt(np.finfo(float).eps)
if not np.all(valid):
Expand Down Expand Up @@ -350,6 +361,9 @@ class ViewMaterialIncompressible(PlotMaterial):
bx : ndarray, optional
Array with stretches for incompressible equi-biaxial tension. Default is
``linsteps([1, 1.75], num=15)``.
statevars : ndarray or None, optional
Array with state variables (default is None). If None, the state variables are
assumed to be initially zero.
Examples
--------
Expand Down Expand Up @@ -382,12 +396,14 @@ def __init__(
ux=linsteps([0.7, 2.5], num=36),
ps=linsteps([1, 2.5], num=30),
bx=linsteps([1, 1.75], num=15),
statevars=None,
):
self.umat = umat
self.ux = ux
self.ps = ps
self.bx = bx
self.statevars_included = self.umat.x[-1].size > 0
self.statevars = statevars

def uniaxial(self, stretches=None):
"""Normal force per undeformed area vs. stretch curve for a uniaxial
Expand All @@ -414,14 +430,15 @@ def uniaxial(self, stretches=None):
F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)

if self.statevars_included:
statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
if self.statevars is None:
self.statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
P = np.zeros_like(F)
for increment, defgrad in enumerate(F.T):
P[..., [increment]], statevars = self.umat.gradient(
[F[..., [increment]], statevars]
P[..., [increment]], self.statevars = self.umat.gradient(
[F[..., [increment]], self.statevars]
)
else:
P, statevars = self.umat.gradient([F, None])
P, self.statevars = self.umat.gradient([F, None])

return λ1, (P[0, 0] - λ3 / λ1 * P[2, 2]).ravel(), "Uniaxial (Incompressible)"

Expand Down Expand Up @@ -451,14 +468,15 @@ def planar(self, stretches=None):
F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)

if self.statevars_included:
statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
if self.statevars is None:
self.statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
P = np.zeros_like(F)
for increment, defgrad in enumerate(F.T):
P[..., [increment]], statevars = self.umat.gradient(
[F[..., [increment]], statevars]
P[..., [increment]], self.statevars = self.umat.gradient(
[F[..., [increment]], self.statevars]
)
else:
P, statevars = self.umat.gradient([F, None])
P, self.statevars = self.umat.gradient([F, None])

return (
λ1,
Expand Down Expand Up @@ -491,13 +509,14 @@ def biaxial(self, stretches=None):
F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)

if self.statevars_included:
statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
if self.statevars is None:
self.statevars = np.zeros((*self.umat.x[-1].shape, 1, 1))
P = np.zeros_like(F)
for increment, defgrad in enumerate(F.T):
P[..., [increment]], statevars = self.umat.gradient(
[F[..., [increment]], statevars]
P[..., [increment]], self.statevars = self.umat.gradient(
[F[..., [increment]], self.statevars]
)
else:
P, statevars = self.umat.gradient([F, None])
P, self.statevars = self.umat.gradient([F, None])

return λ1, (P[0, 0] - λ3 / λ1 * P[2, 2]).ravel(), "Biaxial (Incompressible)"
21 changes: 12 additions & 9 deletions tests/test_constitution.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -307,25 +307,28 @@ def test_umat_hyperelastic(savefig=False):
def neo_hooke(C, mu=1):
return mu / 2 * (tm.linalg.det(C) ** (-1 / 3) * tm.trace(C) - 3)

for model, kwargs in [
(neo_hooke, {"mu": 1}),
(fem.constitution.saint_venant_kirchhoff, {"mu": 1, "lmbda": 20.0}),
(fem.constitution.neo_hooke, {"mu": 1}),
(fem.constitution.mooney_rivlin, {"C10": 0.3, "C01": 0.8}),
(fem.constitution.yeoh, {"C10": 0.5, "C20": -0.1, "C30": 0.02}),
for model, kwargs, incompressible in [
(neo_hooke, {"mu": 1}, True),
(fem.constitution.saint_venant_kirchhoff, {"mu": 1, "lmbda": 20.0}, False),
(fem.constitution.neo_hooke, {"mu": 1}, True),
(fem.constitution.mooney_rivlin, {"C10": 0.3, "C01": 0.8}, True),
(fem.constitution.yeoh, {"C10": 0.5, "C20": -0.1, "C30": 0.02}, True),
(
fem.constitution.third_order_deformation,
{"C10": 0.5, "C01": 0.1, "C11": 0.01, "C20": -0.1, "C30": 0.02},
True,
),
(fem.constitution.ogden, {"mu": [1, 0.2], "alpha": [1.7, -1.5]}),
(fem.constitution.arruda_boyce, {"C1": 1.0, "limit": 3.2}),
(fem.constitution.ogden, {"mu": [1, 0.2], "alpha": [1.7, -1.5]}, True),
(fem.constitution.arruda_boyce, {"C1": 1.0, "limit": 3.2}, True),
(
fem.constitution.extended_tube,
{"Gc": 0.1867, "Ge": 0.2169, "beta": 0.2, "delta": 0.09693},
True,
),
(
fem.constitution.van_der_waals,
{"mu": 1.0, "beta": 0.1, "a": 0.5, "limit": 5.0},
True,
),
]:
umat = fem.Hyperelastic(model, **kwargs)
Expand All @@ -336,7 +339,7 @@ def neo_hooke(C, mu=1):
if savefig:
ax = umat.screenshot(
filename=f"../docs/felupe/images/umat_{umat.fun.__name__}.png",
incompressible=True,
incompressible=incompressible,
)

for incompressible in [False, True]:
Expand Down

0 comments on commit 33bacdd

Please sign in to comment.