diff --git a/matadi/_lab_compressible.py b/matadi/_lab_compressible.py
index 347d424..cbfe36a 100644
--- a/matadi/_lab_compressible.py
+++ b/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]
diff --git a/matadi/_lab_incompressible.py b/matadi/_lab_incompressible.py
index 483b47e..cfcd123 100644
--- a/matadi/_lab_incompressible.py
+++ b/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."
diff --git a/matadi/_templates.py b/matadi/_templates.py
index fb24c0b..e54644b 100644
--- a/matadi/_templates.py
+++ b/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]))]
 
 
diff --git a/matadi/math.py b/matadi/math.py
index caeab8c..6a17c59 100644
--- a/matadi/math.py
+++ b/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
diff --git a/matadi/models/_hyperelasticity_anisotropic.py b/matadi/models/_hyperelasticity_anisotropic.py
index ef9df5c..afc0a8c 100644
--- a/matadi/models/_hyperelasticity_anisotropic.py
+++ b/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
diff --git a/matadi/models/_hyperelasticity_isotropic.py b/matadi/models/_hyperelasticity_isotropic.py
index 1eff63a..6387b7e 100644
--- a/matadi/models/_hyperelasticity_isotropic.py
+++ b/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)
     )
diff --git a/matadi/models/_misc.py b/matadi/models/_misc.py
index 898cd0f..0ab54b1 100644
--- a/matadi/models/_misc.py
+++ b/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))
diff --git a/matadi/models/_templates.py b/matadi/models/_templates.py
index 9ef0c2d..560d9c3 100644
--- a/matadi/models/_templates.py
+++ b/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):
diff --git a/matadi/models/microsphere/_chain.py b/matadi/models/microsphere/_chain.py
index 865c971..29f638a 100644
--- a/matadi/models/microsphere/_chain.py
+++ b/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):
diff --git a/pyproject.toml b/pyproject.toml
index 4152e2b..cfb404d 100644
--- a/pyproject.toml
+++ b/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"
diff --git a/setup.cfg b/setup.cfg
index b423b6c..68afa33 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = matadi
-version = 0.1.10
+version = 0.1.11
 author = Andreas Dutzler
 author_email = a.dutzler@gmail.com
 description = Material Definition with Automatic Differentiation
diff --git a/tests/test_eigvals_grad.py b/tests/test_eigvals_grad.py
index 4742701..af6b0fb 100644
--- a/tests/test_eigvals_grad.py
+++ b/tests/test_eigvals_grad.py
@@ -15,43 +15,43 @@ def fun_mexp(x):
 
 
 def test_eigvals():
-    
+
     # test several repeated principal stretches
     for w in [1, 1.1, 0.1, 2.4, 12]:
 
         # variables
         F = Variable("F", 3, 3)
-    
+
         # data
         np.random.seed(2345537)
         FF = np.random.rand(3, 3, 5, 100)
         for a in range(3):
             FF[a, a] += 1
-    
+
         # input with repeated equal eigenvalues
         FF[:, :, 0, 1] = np.diag(w * np.ones(3))
         FF[:, :, 0, 2] = np.diag([w, w * 1.01, w])
-    
+
         # init Material
         W = Material(x=[F], fun=fun)
-    
+
         WW = W.function([FF])
         dW = W.gradient([FF])
         DW = W.hessian([FF])
-    
+
         Eye = np.eye(3)
         Eye4 = np.einsum("ij,kl->ikjl", Eye, Eye)
-    
+
         # check function
         f = FF[:, :, 0, 0]
         c = f.T @ f
         assert np.isclose(WW[0][0, 0, 0, 0], (np.linalg.eigvals(c).sum() - 3) / 2)
-    
+
         # 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)
@@ -122,24 +122,24 @@ def g(x):
 
 
 def test_mexp():
-    
+
     # test several repeated principal stretches
     for w in [1, 1.1, 0.1, 2.4, 12]:
 
         # variables
         F = Variable("F", 3, 3)
-    
+
         # data
         FF = np.diag([w * 1.01, w, w])
         FF = FF.reshape(3, 3, 1, 1)
-    
+
         # init Material
         W = Material(x=[F], fun=fun_mexp)
-    
+
         WW = W.function([FF])
         dW = W.gradient([FF])
         DW = W.hessian([FF])
-    
+
         assert not np.any(np.isnan(WW))
         assert not np.any(np.isnan(dW))
         assert not np.any(np.isnan(DW))
diff --git a/tests/test_models.py b/tests/test_models.py
index 56ab9aa..aad25b3 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -104,7 +104,7 @@ def test_models():
         DW = HM.hessian_vector_product([FF], [dF], [DF])
 
         assert len(W) == 1
-        assert len(P) == 1
+        assert len(P) == 1 + 1
         assert len(A) == 1
 
         dW_check = np.einsum("ij...,ij...->...", P[0], dF)
@@ -118,7 +118,7 @@ def test_models():
         A = HM_mixed.hessian([FF, pp, JJ])
 
         assert len(W) == 1
-        assert len(P) == 3
+        assert len(P) == 3 + 1
         assert len(A) == 6
 
         dW = HM_mixed.gradient_vector_product([FF, pp, JJ], [dF, dp, dJ])
@@ -156,7 +156,7 @@ def test_models():
         A = HM_mixed2.hessian([FF, pp])
 
         assert len(W) == 1
-        assert len(P) == 2
+        assert len(P) == 2 + 1
         assert len(A) == 3
 
         dW = HM_mixed2.gradient_vector_product([FF, pp], [dF, dp])
@@ -194,7 +194,7 @@ def test_models():
     A = comp.hessian([FF])
 
     assert len(W) == 1
-    assert len(P) == 1
+    assert len(P) == 1 + 1
     assert len(A) == 1
 
     nh_mixed = matadi.ThreeFieldVariation(nh)
@@ -207,7 +207,7 @@ def test_models():
     A = comp.hessian([FF, pp, JJ])
 
     assert len(W) == 1
-    assert len(P) == 3
+    assert len(P) == 3 + 1
     assert len(A) == 6
 
     nh_mixed2 = matadi.TwoFieldVariation(nh)
@@ -220,7 +220,7 @@ def test_models():
     A = comp.hessian([FF, pp])
 
     assert len(W) == 1
-    assert len(P) == 2
+    assert len(P) == 2 + 1
     assert len(A) == 3
 
 
diff --git a/tests/test_template-materials.py b/tests/test_template-materials.py
index f9f5539..367d2ff 100644
--- a/tests/test_template-materials.py
+++ b/tests/test_template-materials.py
@@ -12,9 +12,9 @@ def fun(x, C10=0.5, bulk=5000):
 
     W = neo_hooke(F, C10)
     U = volumetric(J, bulk)
-    
+
     statevars_old = x[-1]
-    statevars_new = ones_like(statevars_old) # only for testing
+    statevars_new = ones_like(statevars_old)  # only for testing
     statevars_new = zeros_like(statevars_old)
 
     return gradient(W, F) + gradient(U, F), statevars_new