Merge pull request #108 from adtzlr/add-morph

Add MORPH material model

matadi/__about__.py

@@ -1 +1 @@
-__version__ = "0.1.6"
+__version__ = "0.1.7"

matadi/math.py

@@ -6,6 +6,8 @@
     transpose,
     trace,
     diag,
+    triu,
+    tril,
     adj,
     cofactor,
     cross,
@@ -121,3 +123,61 @@ def dev(T):
 def ddot(A, B):
 
     return trace(transpose(A) @ B)
+
+
+def tresca(C):
+    "Tresca Invariant as maximum difference of two eigenvalues."
+    wC = eigvals(C, 8e-5)
+    return mmax(fabs(wC[[0, 1, 2]] - wC[[1, 2, 0]]))
+
+
+def mexp(C, eps=8e-5):
+    "Exponential Function of a Matrix."
+    w = eigvals(C, eps=eps)
+    I = SX.eye(3)
+
+    M1 = (C - w[1] * I) * (C - w[2] * I) / (w[0] - w[1]) / (w[0] - w[2])
+    M2 = (C - w[2] * I) * (C - w[0] * I) / (w[1] - w[2]) / (w[1] - w[0])
+    M3 = (C - w[0] * I) * (C - w[1] * I) / (w[2] - w[0]) / (w[2] - w[1])
+
+    return exp(w[0]) * M1 + exp(w[1]) * M2 + exp(w[2]) * M3
+
+
+def asvoigt(A, scale=1):
+
+    if A.shape == (3, 3):
+        return vertcat(
+            A[0, 0], 
+            A[1, 1], 
+            A[2, 2], 
+            A[0, 1] * scale, 
+            A[1, 2] * scale, 
+            A[0, 2] * scale,
+        )
+
+    elif A.shape == (2, 2):
+        return vertcat(
+            A[0, 0], 
+            A[1, 1], 
+            A[0, 1] * scale, 
+        )
+
+    else:
+        raise ValueError("Unknown shape of input.")
+
+
+def astensor(A, scale=1):
+
+    if A.shape == (6, 1):
+        A0 = vertcat(A[0] / scale, A[3] / scale, A[5])
+        A1 = vertcat(A[3] / scale, A[1], A[4] / scale)
+        A2 = vertcat(A[5] / scale, A[4] / scale, A[2])
+        return horzcat(A0, A1, A2)
+
+    elif A.shape == (3, 1):
+        A0 = vertcat(A[0], A[2] / scale)
+        A1 = vertcat(A[2] / scale, A[1])
+        return horzcat(A0, A1)
+
+    else:
+        raise ValueError("Unknown shape of input.")

matadi/models/__init__.py

@@ -24,6 +24,7 @@
 )
 
 from ._pseudo_elasticity import ogden_roxburgh
+from ._misc import morph
 
 from . import microsphere
 from .microsphere.nonaffine import miehe_goektepe_lulei

matadi/models/_misc.py

@@ -0,0 +1,66 @@
+from ._helpers import displacement_pressure_split, volumetric
+from ._hyperelasticity_isotropic import neo_hooke
+from ..math import (
+    det, inv, dev, sym, sqrt, if_else, vertsplit, vertcat, 
+    asvoigt, astensor, tresca, mexp, gradient
+)
+
+
+@displacement_pressure_split
+def morph(x, param, bulk):
+    "MORPH consitututive material formulation."
+
+    # split input into the deformation gradient and the vector of state variables
+    F, statevars = x[0], x[-1]
+
+    # split state variables
+    CTSn, Cn, SZn = vertsplit(statevars, [0, 1, 7, 13])
+    Cn, SZn = astensor(Cn), astensor(SZn)
+
+    # right Cauchy-Green deformation tensor
+    C = F.T @ F
+    J = det(F)
+    dC = C - Cn
+
+    # (Incremental) Inverse  of right Cauchy-Green deformation tensor
+    invC = inv(C)
+
+    # isochoric part of C and L
+    CG = J ** (-2 / 3) * C
+    LG = dev(sym(dC @ invC)) @ CG
+
+    # # von mises invariants of C and L
+    CT = tresca(CG)
+    LT = tresca(LG)
+
+    # # maximum historical von mises invariant of CG
+    CTS = if_else(CT > CTSn, CT, CTSn)
+
+    # # stable LG / LT and CT / CTS
+    LGLT = if_else(LT > 0, LG / LT, LG)
+    CTCTS = if_else(CTS > 0, CT / CTS, CT)
+
+    # MORPH deformation-dependent material parameters
+    f = lambda x: 1 / sqrt(1 + x ** 2)
+    a = param[0] + param[1] * f(param[2] * CTS)
+    b = param[3] * f(param[2] * CTS)
+    c = param[4] * CTS * (1 - f(CTS / param[5]))
+
+    # # Hull stress
+    SH = (c * mexp(param[6] * LGLT * CTCTS) + param[7] * LGLT) @ invC
+
+    # # helper "kappa" for implict euler update of overstress
+    SZ = (SZn + b * LT * SH) / (1 + b * LT)
+
+    # hyperelastic distortional and volumetric parts of strain energy function
+    W = neo_hooke(F, a)
+    U = volumetric(J, bulk)
+
+    # overstress as second Piola-Kirchhoff stress
+    S = dev(SZ @ C) @ invC
+
+    # update statevars
+    statevars_new = vertcat(CTS, asvoigt(C), asvoigt(SZ))
+
+    # first Piola-Kirchhoff stress
+    return [gradient(W, F) + gradient(U, F) + F @ S, statevars_new]

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "matadi"
-version = "0.1.6"
+version = "0.1.7"
 description = "Material Definition with Automatic Differentiation"
 readme = "README.md"
 requires-python = ">=3.6"

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = matadi
-version = 0.1.6
+version = 0.1.7
 author = Andreas Dutzler
 author_email = [email protected]
 description = Material Definition with Automatic Differentiation

tests/test_displacement-pressure.py

@@ -8,7 +8,11 @@
     volumetric,
 )
 
-from matadi.math import det, gradient, dev, inv
+from matadi.math import (
+    det, gradient, dev, inv, asvoigt, astensor, vertcat, triu, zeros
+)
+
+import pytest
 
 
 def test_up_state():
@@ -25,6 +29,26 @@ def fun(x, C10=0.5, bulk=5000, r=3, m=1, beta=0):
 
         # split `x` into the deformation gradient and the state variable
         F, Wmaxn = x[0], x[-1]
+
+        # (begin) test `asvoigt()` and `astensor()`
+        C_3D = F.T @ F
+        C_2D = vertcat(C_3D[0, 0], C_3D[1, 0], C_3D[0, 1], C_3D[1, 1]).reshape((2, 2))
+
+        C_6 = asvoigt(C_3D)
+        C_4 = asvoigt(C_2D)
+
+        C_from_C_6 = astensor(C_6)
+        C_2D_from_C_4 = astensor(C_4)
+
+        assert C_3D[0, 1] == C_from_C_6[0, 1]
+        assert C_2D[0, 1] == C_2D_from_C_4[0, 1]
+
+        with pytest.raises(ValueError):
+            asvoigt(C_6)
+
+        with pytest.raises(ValueError):
+            astensor(C_3D)
+        # (end) test `asvoigt()` and `astensor()`
 
         # isochoric and volumetric parts of the hyperelastic strain energy function
         W = neo_hooke(F, C10)

tests/test_morph.py

@@ -0,0 +1,40 @@
+import numpy as np
+
+from matadi import Variable, MaterialTensor
+from matadi.models import (
+    morph
+)
+
+
+def test_up_morph():
+
+    # deformation gradient
+    F = Variable("F", 3, 3)
+
+    # state variables
+    z = Variable("z", 13, 1)
+
+    kwargs = {
+        "param": [0.035, 0.37, 0.17, 2.4, 0.01, 6.4, 5.5, 0.24],
+        "bulk": 5000
+    }
+
+    p = morph.p
+    M = MaterialTensor(x=[F, p, z], fun=morph, triu=True, statevars=1, kwargs=kwargs)
+
+    FF = np.random.rand(3, 3, 8, 100)
+    pp = np.random.rand(1, 8, 100)
+    zz = np.random.rand(13, 1, 8, 100)
+
+    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])
+
+    assert len(P) == 3
+    assert len(A) == 3
+
+
+if __name__ == "__main__":
+    test_up_morph()