Merge pull request #5 from adtzlr/fix-equal-eigvals

Fix equal eigvals

README.md

@@ -13,7 +13,7 @@ pip install matadi
 ```
 
 ## Usage
-First, we have to define a symbolic variable on which our strain energy function will be based on.
+First, a symbolic variable on which our strain energy function will be based on has to be created.
 
 **Note**: *A variable of matADi is an instance of a symbolic variable of casADi (`casadi.SX.sym`). All `matadi.math` functions are simple links to (symbolic) casADi-functions.*
 
@@ -40,7 +40,7 @@ def neohooke(x, mu=1.0, bulk=200.0):
     return mu * (I1_iso - 3) + bulk * (J - 1) ** 2 / 2
 ```
 
-With this simple Python function we create an instance of a **Material**, which allows extra `args` and `kwargs` to be passed to our strain energy function. This instance now enables the evaluation of both **gradient** (stress) and **hessian** (elasticity) via methods based on automatic differentiation - optionally also on input data containing trailing axes. If necessary, the strain energy density function itself will be evaluated on input data with optional trailing axes by the **function** method.
+With this simple Python function at hand, we create an instance of a **Material**, which allows extra `args` and `kwargs` to be passed to our strain energy function. This instance now enables the evaluation of both **gradient** (stress) and **hessian** (elasticity) via methods based on automatic differentiation - optionally also on input data containing trailing axes. If necessary, the strain energy density function itself will be evaluated on input data with optional trailing axes by the **function** method.
 
 ```python
 Mat = Material(
@@ -61,7 +61,7 @@ A = Mat.hessian([defgrad])[0]
 ```
 
 ## Template classes for hyperelasticity
-matADi provides several template classes for hyperelastic materials. Some common isotropic hyperelastic material formulations are located in `matadi.models`. These strain energy functions have to be passed into an instance of `MaterialHyperelastic`. Usage is exactly the same as described above. To convert a hyperelastic material based on the deformation gradient into a mixed three-field formulation suitable for nearly-incompressible behavior (*displacements*, *pressure* and *volume ratio*) an instance of a `MaterialHyperelastic` class has to be passed to the `ThreeFieldVariation` class.
+matADi provides several simple template classes suitable for simple hyperelastic materials. Some common isotropic hyperelastic material formulations are located in `matadi.models` (see list below). These strain energy functions have to be passed as the `fun` argument into an instance of `MaterialHyperelastic`. Usage is exactly the same as described above. To convert a hyperelastic material based on the deformation gradient into a mixed three-field formulation suitable for nearly-incompressible behavior (*displacements*, *pressure* and *volume ratio*) an instance of a `MaterialHyperelastic` class has to be passed to `ThreeFieldVariation`.
 
 ```python
 
@@ -72,7 +72,9 @@ from matadi.models import neo_hooke
 pressure = np.random.rand(5, 100)
 volratio = np.random.rand(5, 100) / 10 + 1
 
-NH = MaterialHyperelastic(fun=neo_hooke, C10=0.5, bulk=20.0)
+kwargs = {"C10": 0.5, "bulk": 20.0}
+
+NH = MaterialHyperelastic(fun=neo_hooke, **kwargs)
 
 W = NH.function([defgrad])[0]
 P = NH.gradient([defgrad])[0]
@@ -83,15 +85,15 @@ P_upJ = ThreeFieldVariation(NH).gradient([defgrad, pressure, volratio])
 A_upJ = ThreeFieldVariation(NH).hessian([defgrad, pressure, volratio])
 ```
 
-Available isotropic hyperelastic models:
-- [x] Neo-Hooke
-- [x] Mooney-Rivlin
-- [x] Yeoh
-- [x] Third-Order-Deformation
-- [x] Ogden
-- [x] Arruda-Boyce
+Available isotropic hyperelastic material models:
+- [Neo-Hooke](https://en.wikipedia.org/wiki/Neo-Hookean_solid)
+- [Mooney-Rivlin](https://en.wikipedia.org/wiki/Mooney%E2%80%93Rivlin_solid)
+- [Yeoh](https://en.wikipedia.org/wiki/Yeoh_(hyperelastic_model))
+- [Third-Order-Deformation (James-Green-Simpson)](https://onlinelibrary.wiley.com/doi/abs/10.1002/app.1975.070190723)
+- [Ogden](https://en.wikipedia.org/wiki/Ogden_(hyperelastic_model))
+- [Arruda-Boyce](https://en.wikipedia.org/wiki/Arruda%E2%80%93Boyce_model)
 
-Any user-defined isotropic hyperelastic strain energy density function may be passed as the `fun` argument of an instance of `MaterialHyperelastic` by using the following template:
+Any user-defined isotropic hyperelastic strain energy density function may be passed as the `fun` argument of `MaterialHyperelastic` by using the following template:
 
 ```python
 def fun(F, **kwargs):
@@ -100,4 +102,4 @@ def fun(F, **kwargs):
 ```
 
 ## Hints
-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-statements can't be used in combination with a symbolic boolean operation. If you use `eigvals` to symbolically calculate eigenvalues and their corresponding gradients please call gradient and hessian methods with `Mat.gradient([defgrad], modify=True, eps=1e-5)` to avoid the gradient to be filled with NaN's. This is because the gradient of the implemented eigenvalue calculation is not defined for the case of repeated equal eigenvalues. The `modify` argument adds a small number `eps=1e-5` to the diagonal entries of the input data.
+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).

matadi/__about__.py

@@ -1 +1 @@
-__version__ = "0.0.4"
+__version__ = "0.0.5"

matadi/_apply.py

@@ -36,13 +36,3 @@ def rshape(z):
             fun_shape = [(1,)]
 
     return [np.array(o).reshape(*f, *ax, order="F") for o, f in zip(out, fun_shape)]
-
-
-def modify(x, mdify=True, eps=1e-5):
-    if mdify:
-        y = x.copy()
-        y[0, 0] += eps
-        y[1, 1] -= eps
-    else:
-        y = x
-    return y

matadi/_material.py

@@ -1,7 +1,7 @@
 import numpy as np
 import casadi as ca
 
-from ._apply import apply, modify as mdify
+from ._apply import apply
 
 
 class Material:
@@ -60,40 +60,23 @@ def __init__(self, x, fun, args=(), kwargs={}, compress=False):
                 if j >= i:
                     self._idx_hessian.append((*a, *b))
 
-    def function(self, x, modify=None, eps=1e-5):
+    def function(self, x):
         "Return the function."
-        if modify is not None:
-            y = [mdify(y, m, eps) for y, m in zip(x, modify)]
-        else:
-            y = x
-        return apply(
-            y,
-            fun=self._function,
-            x_shape=self._idx_gradient,
-            fun_shape=[()],
-        )
+        return apply(x, fun=self._function, x_shape=self._idx_gradient, fun_shape=[()],)
 
-    def gradient(self, x, modify=None, eps=1e-5):
+    def gradient(self, x):
         "Return list of gradients."
-        if modify is not None:
-            y = [mdify(y, m, eps) for y, m in zip(x, modify)]
-        else:
-            y = x
         return apply(
-            y,
+            x,
             fun=self._gradient,
             x_shape=self._idx_gradient,
             fun_shape=self._idx_gradient,
         )
 
-    def hessian(self, x, modify=None, eps=1e-5):
+    def hessian(self, x):
         "Return upper-triangle entries of hessian."
-        if modify is not None:
-            y = [mdify(y, m, eps) for y, m in zip(x, modify)]
-        else:
-            y = x
         return apply(
-            y,
+            x,
             fun=self._hessian,
             x_shape=self._idx_gradient,
             fun_shape=self._idx_hessian,

matadi/math.py

@@ -10,7 +10,6 @@
     cross,
     times,
     eig_symbolic,
-    eig_symbolic as eigvals,
     ldl,
     # trig
     sin,
@@ -30,4 +29,31 @@
     sum1,
     sum2,
     sumsqr,
+    pi,
+    fabs,
+    linspace,
+    #
+    if_else,
+    SX,
+    DM,
+    MX,
 )
+
+
+def invariants(T):
+
+    I1 = trace(T)
+    I2 = (I1 ** 2 - trace(T @ T)) / 2
+    I3 = det(T)
+
+    return I1, I2, I3
+
+
+def eigvals(T, eps=1e-5):
+
+    T[0, 0] += eps
+    T[1, 1] -= eps
+
+    wT = eig_symbolic(T)
+
+    return wT

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = matadi
-version = 0.0.4
+version = 0.0.5
 author = Andreas Dutzler
 author_email = a.dutzler@gmail.com
 description = Material Definition with Automatic Differentiation

tests/test_eigvals_grad.py

@@ -4,29 +4,32 @@
 from matadi.math import transpose, eigvals, sum1, trace
 
 
+def fun(x):
+    F = x[0]
+    return (sum1(eigvals(transpose(F) @ F))[0, 0] - 3) / 2
+
+
 def test_eigvals():
 
     # variables
     F = Variable("F", 3, 3)
 
-    def fun(x):
-        F = x[0]
-        return (sum1(eigvals(transpose(F) @ F)) - 3) / 2
-
     # data
     np.random.seed(2345537)
     FF = np.random.rand(3, 3, 5, 100)
     for a in range(3):
         FF[a, a] += 1
 
-    FF[:, :, 0, 1] = np.diag([1.2, 0.7, 1.2])
+    # input with repeated equal eigenvalues
+    FF[:, :, 0, 1] = np.diag(2.4 * np.ones(3))
+    FF[:, :, 0, 2] = np.diag([2.4, 1.2, 2.4])
 
     # init Material
     W = Material(x=[F], fun=fun)
 
-    WW = W.function([FF], modify=[True], eps=1e-6)
-    dW = W.gradient([FF], modify=[True], eps=1e-6)
-    DW = W.hessian([FF], modify=[True], eps=1e-6)
+    WW = W.function([FF])
+    dW = W.gradient([FF])
+    DW = W.hessian([FF])
 
     Eye = np.eye(3)
     Eye4 = np.einsum("ij,kl->ikjl", Eye, Eye)
@@ -39,32 +42,29 @@ def fun(x):
     # check gradient
     assert np.allclose(dW[0][:, :, 0, 0], FF[:, :, 0, 0])
     assert np.allclose(dW[0][:, :, 0, 1], FF[:, :, 0, 1])
+    assert np.allclose(dW[0][:, :, 0, 2], FF[:, :, 0, 2])
 
     # check hessian
     assert np.allclose(DW[0][:, :, :, :, 0, 0], Eye4)
     assert np.allclose(DW[0][:, :, :, :, 0, 1], Eye4)
+    assert np.allclose(DW[0][:, :, :, :, 0, 2], Eye4)
 
 
 def test_eigvals_single():
 
     # variables
     F = Variable("F", 3, 3)
 
-    def fun(x):
-        F = x[0]
-
-        return (sum1(eigvals(transpose(F) @ F))[0, 0] - 3) / 2
-
     # data
-    FF = np.diag([1.2, 0.7, 1.2])
+    FF = np.diag(2.4 * np.ones(3))
 
     # init Material
     W = Material(x=[F], fun=fun)
 
     WW = W.function([FF])
-    WW = W.function([FF], modify=[True], eps=1e-6)
-    dW = W.gradient([FF], modify=[True], eps=1e-6)
-    DW = W.hessian([FF], modify=[True], eps=1e-6)
+    WW = W.function([FF])
+    dW = W.gradient([FF])
+    DW = W.hessian([FF])
 
     Eye = np.eye(3)
     Eye4 = np.einsum("ij,kl->ikjl", Eye, Eye)

tests/test_modify.py → tests/test_simple.py

@@ -1,51 +1,44 @@
 import numpy as np
 
 from matadi import Variable, Material
-from matadi.math import det, transpose, trace
+from matadi.math import det, transpose, trace, invariants, sqrt
 
 
-def test_modify():
-
-    # variables
-    F = Variable("F", 3, 3)
-
-    def neohooke(x, mu=1.0, bulk=200.0):
+def neohooke(x, mu=1.0, bulk=200.0):
         """Strain energy density function of nearly-incompressible
         Neo-Hookean isotropic hyperelastic material formulation."""
 
         F = x[0]
-
-        J = det(F)
         C = transpose(F) @ F
-        I1_iso = J ** (-2 / 3) * trace(C)
+
+        I1, I2, I3 = invariants(C)
+        J = sqrt(I3)
+
+        I1_iso = I3 ** (-1 / 3) * trace(C)
 
         return mu * (I1_iso - 3) + bulk * (J - 1) ** 2 / 2
 
+
+def test_simple():
+
+    # variables
+    F = Variable("F", 3, 3)
+
     # data
     FF = np.random.rand(3, 3, 5, 100)
     for a in range(3):
         FF[a, a] += 1
 
     # init Material
-    W = Material(
-        x=[F],
-        fun=neohooke,
-        kwargs={"mu": 1.0, "bulk": 10.0},
-    )
+    W = Material(x=[F], fun=neohooke, kwargs={"mu": 1.0, "bulk": 10.0},)
 
     dW = W.gradient([FF])
     DW = W.hessian([FF])
 
-    dW2 = W.gradient([FF], modify=[True])
-    DW2 = W.hessian([FF], modify=[True])
-
     # dW and DW are always lists...
     assert dW[0].shape == (3, 3, 5, 100)
     assert DW[0].shape == (3, 3, 3, 3, 5, 100)
 
-    assert dW2[0].shape == (3, 3, 5, 100)
-    assert DW2[0].shape == (3, 3, 3, 3, 5, 100)
-
 
 if __name__ == "__main__":
-    test_modify()
+    test_simple()