Merge pull request #122 from adtzlr/unify-material-methods

Unify material methods

matadi/_lab_compressible.py

@@ -103,7 +103,7 @@ def stability(stretch, stretch_3, stretch_3_eps):
             dl = (np.linalg.inv(B) @ df)[0]
 
             # check volume ratio
-            J = stretch ** 2 * stretch_3
+            J = stretch**2 * stretch_3
 
             # check slope of force
             Q = self.material.gradient([G])[0][0, 0]

matadi/_lab_incompressible.py

@@ -21,7 +21,7 @@ def _ux(self, stretch):
 
     def _bx(self, stretch):
         "Principal stretches of incompressible equi-biaxial load case."
-        return stretch, stretch, 1 / stretch ** 2
+        return stretch, stretch, 1 / stretch**2
 
     def _ps(self, stretch):
         "Principal stretches of incompressible planar shear load case."

matadi/_templates.py

@@ -11,20 +11,22 @@ def __init__(self, material):
         self.x = [self.material.x[0], p]
         self.W = Material(self.x, self._fun)
         self.function = self.W.function
-        self.gradient = self.W.gradient
         self.hessian = self.W.hessian
         self.gradient_vector_product = self.W.gradient_vector_product
         self.hessian_vector_product = self.W.hessian_vector_product
 
     def _fun(self, x):
-        F, p = x
+        F, p = x[:2]
         J = det(F)
         W = self.material.fun(F, **self.material.kwargs)
         U = self.material.fun(J ** (1 / 3) * eye(3), **self.material.kwargs)
         dWdF = grad(W, F)
         dWdJ = trace(dWdF @ F.T) / J / 3
         d2WdJdJ = trace(grad(dWdJ, F) @ F.T) / J / 3
-        return W - U + 1 / d2WdJdJ * (p * dWdJ - p ** 2 / 2)
+        return W - U + 1 / d2WdJdJ * (p * dWdJ - p**2 / 2)
+
+    def gradient(self, x, *args, **kwargs):
+        return [*self.W.gradient(x, *args, **kwargs), None]
 
 
 class TwoFieldVariationPlaneStrain:
@@ -34,13 +36,12 @@ def __init__(self, material):
         self.x = [self.material.x[0], p]
         self.W = Material(self.x, self._fun)
         self.function = self.W.function
-        self.gradient = self.W.gradient
         self.hessian = self.W.hessian
         self.gradient_vector_product = self.W.gradient_vector_product
         self.hessian_vector_product = self.W.hessian_vector_product
 
     def _fun(self, x):
-        F, p = x
+        F, p = x[:2]
         F = horzcat(vertcat(x[0], zeros(1, 2)), zeros(3, 1))
         F[2, 2] = 1  # fixed thickness ratio `h / H = 1`
         J = det(F)
@@ -57,7 +58,10 @@ def _fun(self, x):
         dWdJ = Function("w", [f], [dwdj])(F)
         d2WdJdJ = Function("w", [f], [d2wdjdj])(F)
 
-        return W - U + 1 / d2WdJdJ * (p * dWdJ - p ** 2 / 2)
+        return W - U + 1 / d2WdJdJ * (p * dWdJ - p**2 / 2)
+
+    def gradient(self, x, *args, **kwargs):
+        return [*self.W.gradient(x, *args, **kwargs), None]
 
 
 class ThreeFieldVariation:
@@ -68,17 +72,19 @@ def __init__(self, material):
         self.x = [self.material.x[0], p, J]
         self.W = Material(self.x, self._fun)
         self.function = self.W.function
-        self.gradient = self.W.gradient
         self.hessian = self.W.hessian
         self.gradient_vector_product = self.W.gradient_vector_product
         self.hessian_vector_product = self.W.hessian_vector_product
 
     def _fun(self, x):
-        F, p, J = x
+        F, p, J = x[:3]
         detF = det(F)
         Fmod = (J / detF) ** (1 / 3) * F
         return self.material.fun(Fmod, **self.material.kwargs) + p * (detF - J)
 
+    def gradient(self, x, *args, **kwargs):
+        return [*self.W.gradient(x, *args, **kwargs), None]
+
 
 class ThreeFieldVariationPlaneStrain:
     def __init__(self, material):
@@ -88,19 +94,21 @@ def __init__(self, material):
         self.x = [self.material.x[0], p, J]
         self.W = Material(self.x, self._fun)
         self.function = self.W.function
-        self.gradient = self.W.gradient
         self.hessian = self.W.hessian
         self.gradient_vector_product = self.W.gradient_vector_product
         self.hessian_vector_product = self.W.hessian_vector_product
 
     def _fun(self, x):
-        F, p, J = x
+        F, p, J = x[:3]
         F = horzcat(vertcat(x[0], zeros(1, 2)), zeros(3, 1))
         F[2, 2] = 1  # fixed thickness ratio `h / H = 1`
         detF = det(F)
         Fmod = (J / detF) ** (1 / 3) * F
         return self.material.fun(Fmod, **self.material.kwargs) + p * (detF - J)
 
+    def gradient(self, x, *args, **kwargs):
+        return [*self.W.gradient(x, *args, **kwargs), None]
+
 
 class MaterialHyperelastic:
     def __init__(self, fun, **kwargs):
@@ -109,15 +117,21 @@ def __init__(self, fun, **kwargs):
         self.fun = fun
         self.kwargs = kwargs
         self.W = Material(self.x, self._fun_wrapper, kwargs=self.kwargs)
-        self.function = self.W.function
-        self.gradient = self.W.gradient
-        self.hessian = self.W.hessian
         self.gradient_vector_product = self.W.gradient_vector_product
         self.hessian_vector_product = self.W.hessian_vector_product
 
     def _fun_wrapper(self, x, **kwargs):
         return self.fun(x[0], **kwargs)
 
+    def function(self, x, *args, **kwargs):
+        return self.W.function(x[:1], *args, **kwargs)
+
+    def gradient(self, x, *args, **kwargs):
+        return [*self.W.gradient(x[:1], *args, **kwargs), None]
+
+    def hessian(self, x, *args, **kwargs):
+        return self.W.hessian(x[:1], *args, **kwargs)
+
 
 class MaterialHyperelasticPlaneStrain:
     def __init__(self, fun, **kwargs):
@@ -126,9 +140,6 @@ def __init__(self, fun, **kwargs):
         self.fun = fun
         self.kwargs = kwargs
         self.W = Material(self.x, self._fun_wrapper, kwargs=self.kwargs)
-        self.function = self.W.function
-        self.gradient = self.W.gradient
-        self.hessian = self.W.hessian
         self.gradient_vector_product = self.W.gradient_vector_product
         self.hessian_vector_product = self.W.hessian_vector_product
 
@@ -137,6 +148,15 @@ def _fun_wrapper(self, x, **kwargs):
         F[2, 2] = 1  # fixed thickness ratio `h / H = 1`
         return self.fun(F, **kwargs)
 
+    def function(self, x, *args, **kwargs):
+        return self.W.function(x[:1], *args, **kwargs)
+
+    def gradient(self, x, *args, **kwargs):
+        return [*self.W.gradient(x[:1], *args, **kwargs), None]
+
+    def hessian(self, x, *args, **kwargs):
+        return self.W.hessian(x[:1], *args, **kwargs)
+
 
 class MaterialHyperelasticPlaneStressIncompressible(MaterialHyperelasticPlaneStrain):
     def __init__(self, fun, **kwargs):
@@ -166,24 +186,28 @@ def _fun_wrapper(self, x, **kwargs):
 
 class MaterialComposite:
     def __init__(self, materials):
-        "Composite Material as a sum of a list of materials."
+        "Composite Material as a sum of a list of hyperelastic materials."
         self.materials = materials
         self.fun = self.composite
 
+        # get number of variables defined in the first material
+        self._n = len(self.materials[0].x)
+
     def composite(self):
         "Dummy function for plot title."
         return
 
     def function(self, x, **kwargs):
-        fun = [m.function(x, **kwargs) for m in self.materials]
+        fun = [m.function(x[: self._n], **kwargs) for m in self.materials]
         return [np.sum([f[a] for f in fun], 0) for a in range(len(fun[0]))]
 
     def gradient(self, x, **kwargs):
-        grad = [m.gradient(x, **kwargs) for m in self.materials]
-        return [np.sum([g[a] for g in grad], 0) for a in range(len(grad[0]))]
+        grad = [m.gradient(x[: self._n], **kwargs)[: len(x)] for m in self.materials]
+        res = [np.sum([g[a] for g in grad], 0) for a in range(len(grad[0]))]
+        return [*res, None]
 
     def hessian(self, x, **kwargs):
-        hess = [m.hessian(x, **kwargs) for m in self.materials]
+        hess = [m.hessian(x[: self._n], **kwargs) for m in self.materials]
         return [np.sum([h[a] for h in hess], 0) for a in range(len(hess[0]))]
 
 

matadi/math.py

@@ -76,19 +76,19 @@
 
 
 def zeros_like(T):
-    
+
     return zeros(T.shape)
 
 
 def ones_like(T):
-    
+
     return ones(T.shape)
 
 
 def invariants(T):
 
     I1 = trace(T)
-    I2 = (I1 ** 2 - trace(T @ T)) / 2
+    I2 = (I1**2 - trace(T @ T)) / 2
     I3 = det(T)
 
     return I1, I2, I3

matadi/models/_hyperelasticity_anisotropic.py

@@ -32,12 +32,12 @@ def fiber(F, E, angle, k=1, axis=2, compression=False):
     if k == 0:
         strain = log(stretch)
     else:
-        strain = 1 / k * (stretch ** k - 1)
+        strain = 1 / k * (stretch**k - 1)
 
     if not compression:
         strain = if_else(strain < 0, 0, strain)
 
-    return E * strain ** 2 / 2
+    return E * strain**2 / 2
 
 
 @isochoric_volumetric_split

matadi/models/_hyperelasticity_isotropic.py

@@ -11,7 +11,7 @@ def saint_venant_kirchhoff(F, mu, lmbda):
     C = dot(transpose(F), F)
     I1 = trace(C) / 2 - 3 / 2
     I2 = trace(C @ C) / 4 - trace(C) / 2 + 3 / 4
-    return mu * I2 + lmbda * I1 ** 2 / 2
+    return mu * I2 + lmbda * I1**2 / 2
 
 
 @isochoric_volumetric_split
@@ -59,7 +59,7 @@ def ogden(F, mu, alpha):
     out = 0
     for m, a in zip(mu, alpha):
         wk = wC ** (a / 2)
-        out += 2 * m / a ** 2 * (sum1(wk)[0, 0] - 3)
+        out += 2 * m / a**2 * (sum1(wk)[0, 0] - 3)
 
     return out
 
@@ -70,12 +70,12 @@ def arruda_boyce(F, C1, limit):
     I1 = trace(C)
 
     alpha = [1 / 2, 1 / 20, 11 / 1050, 19 / 7000, 519 / 673750]
-    beta = 1 / limit ** 2
+    beta = 1 / limit**2
 
     out = 0
     for i, a in enumerate(alpha):
         j = i + 1
-        out += a * beta ** (2 * j - 2) * (I1 ** j - 3 ** j)
+        out += a * beta ** (2 * j - 2) * (I1**j - 3**j)
 
     return C1 * out
 
@@ -85,9 +85,9 @@ def extended_tube(F, Gc, delta, Ge, beta):
     C = transpose(F) @ F
     D = trace(C)
     wC = eigvals(C)
-    g = (1 - delta ** 2) * (D - 3) / (1 - delta ** 2 * (D - 3))
-    Wc = Gc / 2 * (g + log(1 - delta ** 2 * (D - 3)))
-    We = 2 * Ge / beta ** 2 * sum1(wC ** (-beta / 2) - 1)
+    g = (1 - delta**2) * (D - 3) / (1 - delta**2 * (D - 3))
+    Wc = Gc / 2 * (g + log(1 - delta**2 * (D - 3)))
+    We = 2 * Ge / beta**2 * sum1(wC ** (-beta / 2) - 1)
     return Wc + We
 
 
@@ -97,7 +97,7 @@ def van_der_waals(F, mu, limit, a, beta):
     I1 = trace(C)
     I2 = (trace(C) ** 2 - trace(C @ C)) / 2
     I = (1 - beta) * I1 + beta * I2
-    eta = sqrt((I - 3) / (limit ** 2 - 3))
+    eta = sqrt((I - 3) / (limit**2 - 3))
     return mu * (
-        -(limit ** 2 - 3) * (log(1 - eta) + eta) - 2 / 3 * a * ((I - 3) / 2) ** (3 / 2)
+        -(limit**2 - 3) * (log(1 - eta) + eta) - 2 / 3 * a * ((I - 3) / 2) ** (3 / 2)
     )

matadi/models/_misc.py

@@ -51,7 +51,7 @@ def morph(x, p1, p2, p3, p4, p5, p6, p7, p8):
     CT_CTS = if_else(CTS > 0, CT / CTS, CT)
 
     # MORPH deformation-dependent material parameters
-    f = lambda x: 1 / sqrt(1 + x ** 2)
+    f = lambda x: 1 / sqrt(1 + x**2)
     a = p1 + p2 * f(p3 * CTS)
     b = p4 * f(p3 * CTS)
     c = p5 * CTS * (1 - f(CTS / p6))

matadi/models/_templates.py

@@ -28,7 +28,9 @@ def fun(x, C10, r, m, beta, bulk):
             # for a pseudo-elastic material formulation
             return eta * gradient(W, F) + gradient(U, F), Wmax
 
-        super().__init__(fun=fun, statevars_shape=(1, 1), C10=C10, r=r, m=m, beta=beta, bulk=bulk)
+        super().__init__(
+            fun=fun, statevars_shape=(1, 1), C10=C10, r=r, m=m, beta=beta, bulk=bulk
+        )
 
 
 class Morph(MaterialTensorGeneral):

matadi/models/microsphere/_chain.py

@@ -7,7 +7,7 @@ def langevin(stretch, mu, N):
     Langevin function is defined by a Padé approximation."""
 
     x = stretch / sqrt(N)
-    L = x * (3 - x ** 2) / (1 - x ** 2)
+    L = x * (3 - x**2) / (1 - x**2)
 
     return mu * N * (x * L + log(L / sinh(L)))
 
@@ -22,14 +22,14 @@ def langevin2(stretch, mu, N):
     w.r.t. the stretch gives a real-valued result for the region of interest
     `N > stretch ** 2`."""
 
-    return mu * (stretch ** 2 / 2 - N * log(stretch ** 2 - N))
+    return mu * (stretch**2 / 2 - N * log(stretch**2 - N))
 
 
 def gauss(stretch, mu):
     """Gaussian model given by the free energy
     of a single chain as a function of the stretch."""
 
-    return 3 * mu / 2 * (stretch ** 2 - 1)
+    return 3 * mu / 2 * (stretch**2 - 1)
 
 
 def linear(stretch, mu):

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "matadi"
-version = "0.1.10"
+version = "0.1.11"
 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.10
+version = 0.1.11
 author = Andreas Dutzler
 author_email = [email protected]
 description = Material Definition with Automatic Differentiation