Skip to content

Commit

Permalink
Merge pull request #123 from adtzlr/enhance-tensor-materials
Browse files Browse the repository at this point in the history
Tensor-Materials: Rename `function` to `gradient` and `gradient` to `hessian`
  • Loading branch information
adtzlr authored Nov 10, 2022
2 parents 74738a6 + 0455b09 commit 5390b66
Show file tree
Hide file tree
Showing 8 changed files with 29 additions and 37 deletions.
6 changes: 3 additions & 3 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -195,7 +195,7 @@ c) the sign of the resulting stretch from a small superposed force in one direct
For tensor-valued material definitions use `MaterialTensor` (e.g. any stress-strain relation). Also, please have a look at [casADi's documentation](https://web.casadi.org/). It is very powerful but unfortunately does not support all the Python stuff you would expect. For example Python's default if-else-statements can't be used in combination with symbolic conditions (use `math.if_else(cond, if_true, if_false)` instead). Contrary to [casADi](https://web.casadi.org/), the gradient of the eigenvalue function is stabilized by a perturbation of the diagonal components.

### A **Material** with state variables
A generalized material model with optional state variables for the (u/p)-formulation is created by an instance of `MaterialTensor`. If the argument `triu` is set to `True` the gradient method returns only the upper triangle entries of the gradient components. If some of the input variables are internal state variables the number of these variables have to be passed to the optional argument `statevars`. While the hyperelastic material classes are defined by a strain energy function, this one is defined by the first Piola-Kirchhoff stress tensor. Internally, state variables are equal to default variables but they are excluded from gradient calculations. State variables may also be used as placeholders for additional quantities, e.g. the initial deformation gradient at the beginning of an increment or the time increment. Hence, it is a very flexible class not restricted to hyperelasticity.
A generalized material model with optional state variables for the (u/p)-formulation is created by an instance of `MaterialTensor`. If the argument `triu` is set to `True` the gradient method returns only the upper triangle entries of the gradient components. If some of the input variables are internal state variables the number of these variables have to be passed to the optional argument `statevars`. While the hyperelastic material classes are defined by a strain energy function, this one is defined by the first Piola-Kirchhoff stress tensor. Internally, state variables are equal to default variables but they are excluded from gradient calculations. State variables may also be used as placeholders for additional quantities, e.g. the initial deformation gradient at the beginning of an increment or the time increment. Hence, it is a very flexible class not restricted to hyperelasticity. For consistency, the methods `gradient` and `hessian` of a tensor-based material refer to the gradient and hessian of the strain energy function.

Included [pseudo-elastic material models](https://github.com/adtzlr/matadi/blob/main/matadi/models/_pseudo_elasticity.py):
- [Ogden-Roxburgh](https://doi.org/10.1098%2Frspa.1999.0431) ([code](https://github.com/adtzlr/matadi/blob/main/matadi/models/_pseudo_elasticity.py#L4-L16))
Expand Down Expand Up @@ -248,8 +248,8 @@ statevars = np.random.rand(1, 1, 5, 100)
for a in range(3):
defgrad[a, a] += 1.0

dWdF, dWdp, statevars_new = NH.function([defgrad, pressure, statevars])
d2WdFdF, d2WdFdp, d2Wdpdp = NH.gradient([defgrad, pressure, statevars])
dWdF, dWdp, statevars_new = NH.gradient([defgrad, pressure, statevars])
d2WdFdF, d2WdFdp, d2Wdpdp = NH.hessian([defgrad, pressure, statevars])
```

**Hint**: *The above Neo-Hooke as well as the MORPH material model formulation, both within the u/p-framework, are available as ready-to-go materials in `matadi.models` as [`NeoHookeOgdenRoxburgh()`](https://github.com/adtzlr/matadi/blob/main/matadi/models/_templates.py) and [`Morph()`](https://github.com/adtzlr/matadi/blob/main/matadi/models/_templates.py).*
Expand Down
6 changes: 2 additions & 4 deletions matadi/_material.py
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,7 @@ def __init__(

# init Function
super().__init__(x=x, fun=fun, args=args, kwargs=kwargs)
self.gradient = self.function

# no. of active variables
n = len(self.x) - statevars
Expand All @@ -224,9 +225,6 @@ def __init__(
for f in self._f[:n]
]

# alias
self.jacobian = self.gradient

# store only upper-triangle entries of gradients
if triu:
i, j = np.triu_indices(len(self.x[:n]))
Expand Down Expand Up @@ -266,7 +264,7 @@ def __init__(
else:
self._idx_gradient.append((*a, *b))

def gradient(self, x, threads=cpu_count()):
def hessian(self, x, threads=cpu_count()):
"Return list of gradients."
return apply(
x,
Expand Down
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[project]
name = "matadi"
version = "0.1.11"
version = "0.1.12"
description = "Material Definition with Automatic Differentiation"
readme = "README.md"
requires-python = ">=3.6"
Expand Down
2 changes: 1 addition & 1 deletion setup.cfg
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[metadata]
name = matadi
version = 0.1.11
version = 0.1.12
author = Andreas Dutzler
author_email = a.dutzler@gmail.com
description = Material Definition with Automatic Differentiation
Expand Down
8 changes: 4 additions & 4 deletions tests/test_displacement-pressure.py
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,8 @@ def test_up_state():
for a in range(3):
FF[a, a] += 1

P = M.function([FF, pp, zz]) # stress, constraint, statevars_new
A = M.gradient([FF, pp, zz])
P = M.gradient([FF, pp, zz]) # stress, constraint, statevars_new
A = M.hessian([FF, pp, zz])

assert len(P) == 3
assert len(A) == 3
Expand Down Expand Up @@ -122,8 +122,8 @@ def fun(x):
for a in range(3):
FF[a, a] += 1

P = M.function([FF, pp]) # stress, constraint
A = M.gradient([FF, pp])
P = M.gradient([FF, pp]) # stress, constraint
A = M.hessian([FF, pp])

assert len(P) == 2
assert len(A) == 3
Expand Down
6 changes: 2 additions & 4 deletions tests/test_eigvals_grad.py
Original file line number Diff line number Diff line change
Expand Up @@ -109,15 +109,13 @@ def g(x):
Eye = np.eye(3)
Eye4 = np.einsum("ij,kl->ikjl", Eye, Eye)

WW = W.function([FF])
dW = W.gradient([FF])
DW = W.jacobian([FF])
WW = W.gradient([FF])
dW = W.hessian([FF])

dW_gvp = W.gradient_vector_product([FF, dF])

assert np.allclose(np.einsum("ij...,ij...->...", dW[0], dF), dW_gvp)

assert np.allclose(dW, DW)
assert np.allclose(dW[0][:, :, :, :, 0, 0], Eye4)


Expand Down
28 changes: 12 additions & 16 deletions tests/test_simple.py
Original file line number Diff line number Diff line change
Expand Up @@ -92,19 +92,15 @@ def test_tensor():

W = MaterialTensor(x=[F], fun=fun)

W0 = W.function([FF])
dW = W.gradient([FF])
DW = W.jacobian([FF])

W0 = W.gradient([FF])
dW = W.hessian([FF])
# dW and DW are always lists...
assert W0[0].shape == (3, 3, 8, 1000)
assert dW[0].shape == (3, 3, 3, 3, 8, 1000)
assert DW[0].shape == (3, 3, 3, 3, 8, 1000)

# check output of parallel evaluation
W0 = W.function([FF], threads=2)
dW = W.gradient([FF], threads=2)
DW = W.jacobian([FF], threads=2)
W0 = W.gradient([FF], threads=2)
dW = W.hessian([FF], threads=2)

assert W0[0].shape == (3, 3, 8, 1000)
assert dW[0].shape == (3, 3, 3, 3, 8, 1000)
Expand All @@ -113,37 +109,37 @@ def test_tensor():
pp = np.random.rand(8, 1000)

W = MaterialTensor(x=[p], fun=lambda x: x[0], compress=True)
W0 = W.function([pp], threads=2)
W0 = W.gradient([pp], threads=2)

assert W0[0].shape == (8, 1000)

pp = np.random.rand(1, 1, 8, 1000)

W = MaterialTensor(x=[p], fun=lambda x: x[0])
W0 = W.function([pp])
W0 = W.gradient([pp])

assert W0[0].shape == (1, 1, 8, 1000)

# init mixed Material with upper triangle gradient
W = MaterialTensor(x=[F, p], fun=lambda x: x, triu=True)
P = W.function([FF, pp])
A = W.gradient([FF, pp])
P = W.gradient([FF, pp])
A = W.hessian([FF, pp])

assert len(P) == 2
assert len(A) == 3

# init mixed Material with full gradient
W = MaterialTensor(x=[F, p], fun=lambda x: x, triu=False)
P = W.function([FF, pp])
A = W.gradient([FF, pp])
P = W.gradient([FF, pp])
A = W.hessian([FF, pp])

assert len(P) == 2
assert len(A) == 4

# init mixed Material with upper triangle gradient and state variables
W = MaterialTensor(x=[F, p, z], fun=lambda x: [*x, z], triu=True, statevars=1)
P = W.function([FF, pp, zz])
A = W.gradient([FF, pp, zz])
P = W.gradient([FF, pp, zz])
A = W.hessian([FF, pp, zz])

assert len(P) == 4
assert len(A) == 3
Expand Down
8 changes: 4 additions & 4 deletions tests/test_template-materials.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,8 @@ def test_u_templates():
for a in range(3):
FF[a, a] += 1

P = M.function([FF, zz]) # stress, constraint, statevars_new
A = M.gradient([FF, zz])
P = M.gradient([FF, zz]) # stress, constraint, statevars_new
A = M.hessian([FF, zz])

assert len(P) == 2
assert len(A) == 1
Expand All @@ -54,8 +54,8 @@ def test_up_templates():
for a in range(3):
FF[a, a] += 1

P = M.function([FF, pp, zz]) # stress, constraint, statevars_new
A = M.gradient([FF, pp, zz])
P = M.gradient([FF, pp, zz]) # stress, constraint, statevars_new
A = M.hessian([FF, pp, zz])

assert len(P) == 3
assert len(A) == 3
Expand Down

0 comments on commit 5390b66

Please sign in to comment.