diff --git a/src/hyper/elasticity.py b/src/hyper/elasticity.py
index bfb7393..202096e 100644
--- a/src/hyper/elasticity.py
+++ b/src/hyper/elasticity.py
@@ -1,9 +1,11 @@
 # -*- coding: utf-8 -*-
 
 from abc import ABC, abstractmethod  # For abstract classes
-from . import tensor
+
 import numpy as np
 
+from . import tensor
+
 
 class Elasticity(ABC):
     """
@@ -11,10 +13,11 @@ class Elasticity(ABC):
     1ST_lamdaE_CONSTANT is the frist lamé parameter lambda
     2ST_lamdaE_CONSTANT is the second lamé parameter mu
     """
+
     prop = dict()
 
     def __init__(self, E, nu):
-        G = E / (2 * (1. + nu))
+        G = E / (2 * (1.0 + nu))
         K = E / (3 * (1 - 2 * nu))
         lamda = K - (2 / 3) * G
         mu = G
@@ -88,7 +91,9 @@ def potential(self, C):
         lamda = self.get1LAME()
         mu = self.get2LAME()
         EL = 0.5 * (C - tensor.I(len(C)))  # Lagrangian strain E
-        phi = lamda / 2. * (tensor.trace(EL))**2 + mu * np.tensordot(EL, EL, 2)
+        phi = lamda / 2.0 * (tensor.trace(EL)) ** 2 + mu * np.tensordot(
+            EL, EL, 2
+        )
         return phi
 
     def stress(self, C):
@@ -139,9 +144,9 @@ def potential(self, C):
             C = np.zeros((3, 3))
             C[:2, :2] = K[:, :]
         J = np.sqrt(tensor.det(C))  # J = det(F) and det(C) = J^2
-        part1 = (mu / 2) * (tensor.trace(C) - 3.)
+        part1 = (mu / 2) * (tensor.trace(C) - 3.0)
         part2 = mu * np.log(J)
-        part3 = (lamda / 2) * (np.log(J))**2
+        part3 = (lamda / 2) * (np.log(J)) ** 2
         phi = part1 - part2 + part3
         return phi
 
diff --git a/src/hyper/fem.py b/src/hyper/fem.py
index 185cba8..7bc79e2 100644
--- a/src/hyper/fem.py
+++ b/src/hyper/fem.py
@@ -5,15 +5,15 @@
 
 
 import numpy as np
-from . import tensor
-from . import geometry
 from scipy.linalg import polar
 
+from . import geometry, tensor
 
 ########################################################################
 #          Definition of elements for Finite Elements Method           #
 ########################################################################
 
+
 class FiniteElement:
     """
     Data structure for finite element
@@ -26,13 +26,13 @@ def __init__(self, t, xNod):
         #
         self.calculate_shapes(t, xNod)
 
-        self.F = []        # deformation gradient F = Grad(u) + I
-        self.hencky = []   # hencky strain ln(V) with F = V.R
-        self.E_GL = []     # lagrangian strain E = 1/2 (C - I)
-        self.E_EA = []     # E_EA strain e = 1/2 (I - b^-1)
-        self.PK1 = []      # piola kirchoff I : P = F*S
-        self.sigma = []    # cauchy stress : \sigma
-        self.K = []        # lagrangian tangent operator dP/dF
+        self.F = []  # deformation gradient F = Grad(u) + I
+        self.hencky = []  # hencky strain ln(V) with F = V.R
+        self.E_GL = []  # lagrangian strain E = 1/2 (C - I)
+        self.E_EA = []  # E_EA strain e = 1/2 (I - b^-1)
+        self.PK1 = []  # piola kirchoff I : P = F*S
+        self.sigma = []  # cauchy stress : \sigma
+        self.K = []  # lagrangian tangent operator dP/dF
         d = self.getDim()  # space dimension
         npg = self.getNIntPts()
         for n in range(npg):
@@ -133,14 +133,14 @@ def update(self, uNod, mater):
             self.E_GL[i] = 0.5 * (C - tensor.I(len(C)))
             # compute spatial description
             # from scipy.linalg.polar() method with u=R and p=V,
-            _, V = polar(F, side='left')
+            _, V = polar(F, side="left")
             # replace pure zeros by very low values to prevent "nan" in np.log(V)
             V[V < 1e-10] = 1e-15
             self.hencky[i] = np.log(V)  # ln(V) with F = V.R, "true" strain
             b = tensor.leftCauchyGreen(F)
             self.E_EA[i] = 0.5 * (tensor.I(len(b)) - tensor.inv(b))
             ###
-            if (mater == 0):  # skip next lines: do not compute stress
+            if mater == 0:  # skip next lines: do not compute stress
                 # print("Whahat")
                 continue
             # compute PK2 stress tensor and material tangent operator M=2*dS/dC
@@ -224,7 +224,7 @@ def __init__(self, theMesh, m=0):
                 nodei = theMesh.getNode(label_nodei)
                 xnodei = nodei.getX()
                 # Since we don't need z, we cut it out
-                xnodei = xnodei[:self.dim]
+                xnodei = xnodei[: self.dim]
                 xNod.append(xnodei)
             xNod = np.array(xNod)
             new_FiniteElement = FiniteElement(elem.type, xNod)
@@ -280,7 +280,8 @@ def extract(self, U, e):
         """
         if not isinstance(e, (int)):
             raise Exception(
-                "The element 'e' in 'extract' function is not a integer")
+                "The element 'e' in 'extract' function is not a integer"
+            )
         # print("index = " + str(index))
         connect_e = self.connect[e]
         # print("connect_e = " + str(connect_e))
@@ -409,7 +410,7 @@ def getStressPK1(self, n):
         return avg
 
     def getStrainEulerAlmansi(self, n):
-        """ Return average of EA strain at all integration points of element n"""
+        """Return average of EA strain at all integration points of element n"""
         avg = tensor.tensor(self.dim)
         for tens in self.elems[n].E_EA:
             avg += tens
diff --git a/src/hyper/geometry.py b/src/hyper/geometry.py
index dc64c95..ddab984 100644
--- a/src/hyper/geometry.py
+++ b/src/hyper/geometry.py
@@ -1,6 +1,5 @@
 import numpy as np
 
-
 ########################################################################
 #                           Gauss Quadrature                           #
 ########################################################################
@@ -33,6 +32,7 @@
 #   |           |  (-b, -b), (0, -b), (b, -b) | 25/81, 40/81, 25/81  |
 #    ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾
 
+
 class IPTri:
     """
     Numerical integration rules for triangles
@@ -41,15 +41,20 @@ class IPTri:
         P2 = (1, 0)
         P3 = (0, 1)
     """
+
     @staticmethod
     def getX(npg=3):
         if npg == 1:
-            X = [(1. / 3., 1 / 3.)]
+            X = [(1.0 / 3.0, 1 / 3.0)]
         elif npg == 3:
             X = [(2 / 3, 1 / 6), (1 / 6, 1 / 6), (1 / 6, 2 / 3)]
         elif npg == 4:
-            X = [(1 / 3, 1 / 3), (3 / 5, 1 / 5),
-                 (1 / 5, 1 / 5), (1 / 5, 3 / 5)]
+            X = [
+                (1 / 3, 1 / 3),
+                (3 / 5, 1 / 5),
+                (1 / 5, 1 / 5),
+                (1 / 5, 3 / 5),
+            ]
         else:
             raise ValueError("Not found npg = " + str(npg))
         return X
@@ -57,7 +62,7 @@ def getX(npg=3):
     @staticmethod
     def getW(npg=3):
         if npg == 1:
-            W = [1. / 2.]
+            W = [1.0 / 2.0]
         elif npg == 3:
             W = [1 / 6, 1 / 6, 1 / 6]
         elif npg == 4:
@@ -86,9 +91,17 @@ def getX(npg=3):
             X = [(-a, a), (-a, -a), (a, -a), (a, a)]
         elif npg == 9:
             b = np.sqrt(3 / 5)
-            X = [(-b, b), (0, b), (b, b),
-                 (-b, 0), (0, 0), (b, 0),
-                 (-b, -b), (0, -b), (b, -b)]
+            X = [
+                (-b, b),
+                (0, b),
+                (b, b),
+                (-b, 0),
+                (0, 0),
+                (b, 0),
+                (-b, -b),
+                (0, -b),
+                (b, -b),
+            ]
         else:
             raise ValueError("Not found npg = " + str(npg))
         return X
@@ -100,9 +113,17 @@ def getW(npg=4):
         elif npg == 4:
             W = [1, 1, 1, 1]
         elif npg == 9:
-            W = [25 / 81, 40 / 81, 25 / 81,
-                 40 / 81, 64 / 81, 40 / 81,
-                 25 / 81, 40 / 81, 25 / 81]
+            W = [
+                25 / 81,
+                40 / 81,
+                25 / 81,
+                40 / 81,
+                64 / 81,
+                40 / 81,
+                25 / 81,
+                40 / 81,
+                25 / 81,
+            ]
         else:
             raise ValueError("Not found npg = " + str(npg))
         return W
@@ -131,6 +152,7 @@ def getW(npg=4):
 #  (-1, -1)‾‾‾‾‾‾‾‾‾(1, -1)
 #
 
+
 class SFT3:
     """
     Shape functions for 3-node triangle
@@ -139,6 +161,7 @@ class SFT3:
         P2 = (1, 0)
         P3 = (0, 1)
     """
+
     @staticmethod
     def N(x):
         N1 = 1 - x[0] - x[1]
@@ -178,6 +201,7 @@ class SFT6:
         P5 = (0, 1)
         P6 = (0, 0.5)
     """
+
     @staticmethod
     def N(x):
         # Ni(Pj) = delta_ij
@@ -219,6 +243,7 @@ class SFQ4:
         P3 = (1, 1)
         P4 = (-1, 1)
     """
+
     @staticmethod
     def N(x):
         N1 = (x[0] - 1) * (x[1] - 1) / 4
@@ -257,6 +282,7 @@ class SFQ8:
         P7 = (-1, 1)
         P8 = (-1, 0)
     """
+
     @staticmethod
     def N(x):
         N1 = -(x[0] - 1) * (x[1] - 1) * (x[0] + x[1] + 1) / 4
@@ -271,18 +297,26 @@ def N(x):
 
     @staticmethod
     def dN(x):
-        dN1 = [(2 * x[0] + x[1]) * (1 - x[1]) / 4,
-               (1 - x[0]) * (x[0] + 2 * x[1]) / 4]
-        dN2 = [x[0] * (x[1] - 1), (x[0]**2 - 1) / 2]
-        dN3 = [(-2 * x[0] + x[1]) * (x[1] - 1) / 4,
-               (-x[0] + 2 * x[1]) * (x[0] + 1) / 4]
-        dN4 = [(1 - x[1]**2) / 2, -x[1] * (x[0] + 1)]
-        dN5 = [(2 * x[0] + x[1]) * (x[1] + 1) / 4,
-               (x[0] + 1) * (x[0] + 2 * x[1]) / 4]
-        dN6 = [-x[0] * (x[1] + 1), (1 - x[0]**2) / 2]
-        dN7 = [(2 * x[0] - x[1]) * (x[1] + 1) / 4,
-               (x[0] - 1) * (x[0] - 2 * x[1]) / 4]
-        dN8 = [(x[1]**2 - 1) / 2, x[1] * (x[0] - 1)]
+        dN1 = [
+            (2 * x[0] + x[1]) * (1 - x[1]) / 4,
+            (1 - x[0]) * (x[0] + 2 * x[1]) / 4,
+        ]
+        dN2 = [x[0] * (x[1] - 1), (x[0] ** 2 - 1) / 2]
+        dN3 = [
+            (-2 * x[0] + x[1]) * (x[1] - 1) / 4,
+            (-x[0] + 2 * x[1]) * (x[0] + 1) / 4,
+        ]
+        dN4 = [(1 - x[1] ** 2) / 2, -x[1] * (x[0] + 1)]
+        dN5 = [
+            (2 * x[0] + x[1]) * (x[1] + 1) / 4,
+            (x[0] + 1) * (x[0] + 2 * x[1]) / 4,
+        ]
+        dN6 = [-x[0] * (x[1] + 1), (1 - x[0] ** 2) / 2]
+        dN7 = [
+            (2 * x[0] - x[1]) * (x[1] + 1) / 4,
+            (x[0] - 1) * (x[0] - 2 * x[1]) / 4,
+        ]
+        dN8 = [(x[1] ** 2 - 1) / 2, x[1] * (x[0] - 1)]
         return [dN1, dN2, dN3, dN4, dN5, dN6, dN7, dN8]
 
     @staticmethod
diff --git a/src/hyper/gmsh2.py b/src/hyper/gmsh2.py
index 7f9666f..2a9ac85 100644
--- a/src/hyper/gmsh2.py
+++ b/src/hyper/gmsh2.py
@@ -3,6 +3,7 @@
 # Functions to write Gmsh output files
 #
 import numpy as np
+
 from . import mesh as msh
 
 #
@@ -32,7 +33,7 @@ def gmshInput_mesh(inp, d=2):
     # --- version 2
     #
     if fmt == 2:
-        while (l < len(lines)):
+        while l < len(lines):
             if lines[l].startswith("$Nodes"):
                 l += 1
                 nNodes = int(lines[l].strip())
@@ -47,8 +48,12 @@ def gmshInput_mesh(inp, d=2):
                     elif d == 2:
                         mesh.addNode(lbl, float(words[1]), float(words[2]))
                     elif d == 3:
-                        mesh.addNode(lbl, float(words[1]), float(
-                            words[2]), float(words[3]))
+                        mesh.addNode(
+                            lbl,
+                            float(words[1]),
+                            float(words[2]),
+                            float(words[3]),
+                        )
             elif lines[l].startswith("$Elements"):
                 l += 1
                 nElems = int(lines[l].strip())
@@ -59,16 +64,16 @@ def gmshInput_mesh(inp, d=2):
                     typ = int(words[1])
                     ntg = int(words[2])
                     iNodes = []
-                    if (d == 1 and typ == 1):  # 2-node edge
+                    if d == 1 and typ == 1:  # 2-node edge
                         iNodes.append(iTab[int(words[ntg + 3])])
                         iNodes.append(iTab[int(words[ntg + 4])])
                         mesh.addElement(lbl, typ, iNodes)
-                    elif (d == 2 and typ == 2):  # 3-node triangle
+                    elif d == 2 and typ == 2:  # 3-node triangle
                         iNodes.append(iTab[int(words[ntg + 3])])
                         iNodes.append(iTab[int(words[ntg + 4])])
                         iNodes.append(iTab[int(words[ntg + 5])])
                         mesh.addElement(lbl, typ, iNodes)
-                    elif (d == 2 and typ == 3):  # 4-node quadrangle
+                    elif d == 2 and typ == 3:  # 4-node quadrangle
                         iNodes.append(iTab[int(words[ntg + 3])])
                         iNodes.append(iTab[int(words[ntg + 4])])
                         iNodes.append(iTab[int(words[ntg + 5])])
@@ -80,7 +85,7 @@ def gmshInput_mesh(inp, d=2):
     # --- version 4
     #
     elif fmt == 4:
-        while (l < len(lines)):
+        while l < len(lines):
             if lines[l].startswith("$Nodes"):
                 l += 1
                 words = lines[l].split()
@@ -97,7 +102,8 @@ def gmshInput_mesh(inp, d=2):
                     if param != 0:
                         # != instead of <> in Python 3
                         raise ValueError(
-                            'parametric nodes are not handled for now')
+                            "parametric nodes are not handled for now"
+                        )
                     numNodes = int(words[3])
                     for n in range(numNodes):
                         l += 1
@@ -110,8 +116,12 @@ def gmshInput_mesh(inp, d=2):
                         elif d == 2:
                             mesh.addNode(lbl, float(words[1]), float(words[2]))
                         elif d == 3:
-                            mesh.addNode(lbl, float(words[1]), float(
-                                words[2]), float(words[3]))
+                            mesh.addNode(
+                                lbl,
+                                float(words[1]),
+                                float(words[2]),
+                                float(words[3]),
+                            )
             elif lines[l].startswith("$Elements"):
                 l += 1
                 words = lines[l].split()
@@ -129,16 +139,16 @@ def gmshInput_mesh(inp, d=2):
                         words = lines[l].split()
                         lbl = int(words[0])
                         iNodes = []
-                        if (d == 1 and typ == 1):  # 2-node edge
+                        if d == 1 and typ == 1:  # 2-node edge
                             iNodes.append(iTab[int(words[1])])
                             iNodes.append(iTab[int(words[2])])
                             mesh.addElement(lbl, typ, iNodes)
-                        elif (d == 2 and typ == 2):  # 3-node triangle
+                        elif d == 2 and typ == 2:  # 3-node triangle
                             iNodes.append(iTab[int(words[1])])
                             iNodes.append(iTab[int(words[2])])
                             iNodes.append(iTab[int(words[3])])
                             mesh.addElement(lbl, typ, iNodes)
-                        elif (d == 2 and typ == 3):  # 4-node quadrangle
+                        elif d == 2 and typ == 3:  # 4-node quadrangle
                             iNodes.append(iTab[int(words[1])])
                             iNodes.append(iTab[int(words[2])])
                             iNodes.append(iTab[int(words[3])])
@@ -149,9 +159,11 @@ def gmshInput_mesh(inp, d=2):
         assert numNode + 1 == nNodes
     else:
         raise ValueError(
-            'Format of the meshfile is incorrect. Verify the header $MeshFormat')
+            "Format of the meshfile is incorrect. Verify the header $MeshFormat"
+        )
     return mesh
 
+
 #
 # Write mesh data in Gmsh v.2 format
 # INPUT: - out : output stream
@@ -181,14 +193,16 @@ def gmshOutput_mesh(out, mesh):
     out.write("$Elements\n")
     out.write("%d\n" % mesh.nElements())
     for n in range(mesh.nElements()):
-        out.write("%d %d %d %d %d" %
-                  (n + 1, mesh.getElement(n).getType(), 2, 1, 1))
+        out.write(
+            "%d %d %d %d %d" % (n + 1, mesh.getElement(n).getType(), 2, 1, 1)
+        )
         for i in range(mesh.getElement(n).nNodes()):
             out.write(" %d" % (mesh.getElement(n).getNode(i) + 1))
         out.write("\n")
     out.write("$EndElements\n")
     out.flush()
 
+
 #
 # Write nodal field defined on a mesh, in Gmsh V.2 format
 # INPUT: - out : output stream
@@ -205,15 +219,15 @@ def gmshOutput_nodal(out, label, V, it, t):
     """Write nodal field in Gmsh format"""
     out.write("$NodeData\n")
     out.write("%d\n" % 1)
-    out.write("\"%s\"\n" % label)
+    out.write('"%s"\n' % label)
     out.write("%d\n" % 1)
     out.write("%f\n" % t)
     out.write("%d\n" % 3)
     out.write("%d\n" % it)
     sz0 = np.shape(V)[1]
-    if (sz0 == 2):
+    if sz0 == 2:
         sz1 = 3
-    elif(sz0 > 3):
+    elif sz0 > 3:
         sz1 = 9
     else:
         sz1 = 1
@@ -244,15 +258,15 @@ def gmshOutput_element(out, label, V, it, t):
     """Write element field in Gmsh format"""
     out.write("$ElementData\n")
     out.write("%d\n" % 1)
-    out.write("\"%s\"\n" % label)
+    out.write('"%s"\n' % label)
     out.write("%d\n" % 1)
     out.write("%f\n" % t)
     out.write("%d\n" % 3)
     out.write("%d\n" % it)
     sz0 = np.shape(V)[1]
-    if (sz0 == 2 or sz0 == 3):
+    if sz0 == 2 or sz0 == 3:
         sz1 = 3
-    elif(sz0 > 3):
+    elif sz0 > 3:
         sz1 = 9
     else:
         sz1 = 1
@@ -260,15 +274,17 @@ def gmshOutput_element(out, label, V, it, t):
     out.write("%d\n" % np.shape(V)[0])
     for n in range(np.shape(V)[0]):
         out.write("%d" % (n + 1))
-        if (sz1 == 3):
+        if sz1 == 3:
             for i in range(sz0):
                 out.write(" %13.6e" % V[n][i])
             for i in range(sz0, sz1):
                 out.write(" %13.6e" % 0.0)
-        elif (sz1 == 9):
-            if (sz0 == 4):
-                out.write(" %13.6e %13.6e %13.6e %13.6e %13.6e" %
-                          (V[n][0], V[n][1], 0.0, V[n][2], V[n][3]))
+        elif sz1 == 9:
+            if sz0 == 4:
+                out.write(
+                    " %13.6e %13.6e %13.6e %13.6e %13.6e"
+                    % (V[n][0], V[n][1], 0.0, V[n][2], V[n][3])
+                )
                 for i in range(4):
                     out.write(" %13.6e" % 0.0)
             else:
diff --git a/src/hyper/mesh.py b/src/hyper/mesh.py
index f310e40..096bbf7 100644
--- a/src/hyper/mesh.py
+++ b/src/hyper/mesh.py
@@ -45,7 +45,7 @@ class Element:
     def __init__(self, i, t, n):
         """Create element of label i, type t, with list of nodes n"""
         self.id = i
-        self.type = t    # 2:Triangle, 3:Quadrangle
+        self.type = t  # 2:Triangle, 3:Quadrangle
         self.nodes = []  # list of nodes labels (int)
         self.nodes.extend(n)
 
@@ -75,7 +75,9 @@ def __init__(self, d=2):
         """Create mesh of dimension d"""
         self.dim = d  # dimension of the mesh as an integer, 2 by default
         self.nodes = []  # list of all the nodes in the mesh as Node instances
-        self.elems = []  # list of all the elements in the mesh as Element instances
+        self.elems = (
+            []
+        )  # list of all the elements in the mesh as Element instances
 
     def getDimension(self):
         # Returns the dimention of the Mesh
@@ -84,7 +86,7 @@ def getDimension(self):
     def getDim(self):
         return self.getDimension()
 
-    def addNode(self, i, x, y=0., z=0.):
+    def addNode(self, i, x, y=0.0, z=0.0):
         self.nodes.append(Node(i, x, y, z))
 
     def getNode(self, i):
diff --git a/src/hyper/plot.py b/src/hyper/plot.py
index 54a4239..8edc7a3 100644
--- a/src/hyper/plot.py
+++ b/src/hyper/plot.py
@@ -37,4 +37,3 @@ def plot_mesh(mesh, label=None, color=None):
 
     plt.triplot(X_all, Y_all, all_connections, color=color, label=label)
     plt.axis("equal")
-
diff --git a/src/hyper/tensor.py b/src/hyper/tensor.py
index 08e49df..d99b049 100644
--- a/src/hyper/tensor.py
+++ b/src/hyper/tensor.py
@@ -10,11 +10,11 @@
 import numpy as np
 import numpy.linalg as la
 
-
 ###############################################################################
 #                                   Vectors                                   #
 ###############################################################################
 
+
 def vector(d=2):
     """
     Constructor of a vector object (dimension d)
@@ -26,6 +26,7 @@ def vector(d=2):
 #                             2nd Order Tensors                               #
 ###############################################################################
 
+
 def tensor(d=2):
     """
     Constructor of 2nd-order tensor (dimension d)
@@ -113,6 +114,7 @@ def PK1toCauchy(F, P):
 def check_symmetric(a, rtol=1e-05, atol=1e-08):
     return np.allclose(a, a.T, rtol=rtol, atol=atol)
 
+
 ###############################################################################
 #                             4th Order Tensors                               #
 ###############################################################################
@@ -192,7 +194,7 @@ def MaterialToLagrangian(F, S, M):
     # print("M = " + str(M.shape))
 
     # value1 = np.einsum('iI,IJKL,kK->iJkL', F, M, F)
-    dPdF += np.einsum('iI,IJKL,kK->iJkL', F, M, F)
+    dPdF += np.einsum("iI,IJKL,kK->iJkL", F, M, F)
     # ident = np.eye(d)
     # w = outerProd4(ident, S)
     dPdF += np.einsum("ik,JL->iJkL", np.eye(d), S)
diff --git a/tests/test_elasticity.py b/tests/test_elasticity.py
index ded08db..9d1bbc8 100644
--- a/tests/test_elasticity.py
+++ b/tests/test_elasticity.py
@@ -1,11 +1,11 @@
-#, -*- coding: utf-8, -*-
+# , -*- coding: utf-8, -*-
 """
 """
 
-import pytest
-from hyper import elasticity
-from hyper import tensor
 import numpy as np
+import pytest
+
+from hyper import elasticity, tensor
 
 
 def test_StVenantKirchhoffValues():
@@ -91,24 +91,21 @@ def test_StVenantKirchhoffStress():
 
     C = np.ones((2, 2))
     PK2test = material.stress(C)
-    PK2good = 10 * np.array([[0, 1],
-                             [1, 0]])
+    PK2good = 10 * np.array([[0, 1], [1, 0]])
     assert PK2test.ndim == 2
     assert PK2test.shape == (2, 2)
     np.testing.assert_almost_equal(PK2test, PK2good)
 
     C = 2 * np.eye(2) + 1
     PK2test = material.stress(C)
-    PK2good = 10 * np.array([[5, 1],
-                             [1, 5]])
+    PK2good = 10 * np.array([[5, 1], [1, 5]])
     assert PK2test.ndim == 2
     assert PK2test.shape == (2, 2)
     np.testing.assert_almost_equal(PK2test, PK2good)
 
     C = np.eye(2) + 1
     PK2test = material.stress(C)
-    PK2good = 5 * np.array([[5, 2],
-                            [2, 5]])
+    PK2good = 5 * np.array([[5, 2], [2, 5]])
     assert PK2test.ndim == 2
     assert PK2test.shape == (2, 2)
     np.testing.assert_almost_equal(PK2test, PK2good)
@@ -119,10 +116,10 @@ def test_StVenantKirchhoffStiffness():
     nu = 0.3
     material = elasticity.StVenantKirchhoffElasticity(E, nu)
 
-    Mgood = [[[[7, 0], [0, 3]],
-              [[0, 2], [2, 0]]],
-             [[[0, 2], [2, 0]],
-              [[3, 0], [0, 7]]]]
+    Mgood = [
+        [[[7, 0], [0, 3]], [[0, 2], [2, 0]]],
+        [[[0, 2], [2, 0]], [[3, 0], [0, 7]]],
+    ]
     Mgood = 5 * np.array(Mgood)
 
     C = np.eye(2)  # No deformation at all
@@ -150,15 +147,23 @@ def test_StVenantKirchhoffStiffness():
     np.testing.assert_almost_equal(Mtest, Mgood)
 
     C = np.eye(3)  # No deformation at all
-    Mgood = [[[[7, 0, 0], [0, 3, 0], [0, 0, 3]],
-              [[0, 2, 0], [2, 0, 0], [0, 0, 0]],
-              [[0, 0, 2], [0, 0, 0], [2, 0, 0]]],
-             [[[0, 2, 0], [2, 0, 0], [0, 0, 0]],
-              [[3, 0, 0], [0, 7, 0], [0, 0, 3]],
-              [[0, 0, 0], [0, 0, 2], [0, 2, 0]]],
-             [[[0, 0, 2], [0, 0, 0], [2, 0, 0]],
-              [[0, 0, 0], [0, 0, 2], [0, 2, 0]],
-              [[3, 0, 0], [0, 3, 0], [0, 0, 7]]]]
+    Mgood = [
+        [
+            [[7, 0, 0], [0, 3, 0], [0, 0, 3]],
+            [[0, 2, 0], [2, 0, 0], [0, 0, 0]],
+            [[0, 0, 2], [0, 0, 0], [2, 0, 0]],
+        ],
+        [
+            [[0, 2, 0], [2, 0, 0], [0, 0, 0]],
+            [[3, 0, 0], [0, 7, 0], [0, 0, 3]],
+            [[0, 0, 0], [0, 0, 2], [0, 2, 0]],
+        ],
+        [
+            [[0, 0, 2], [0, 0, 0], [2, 0, 0]],
+            [[0, 0, 0], [0, 0, 2], [0, 2, 0]],
+            [[3, 0, 0], [0, 3, 0], [0, 0, 7]],
+        ],
+    ]
     Mgood = 5 * np.array(Mgood)
     Mtest = material.stiffness(C)
     assert Mtest.ndim == 4
@@ -247,24 +252,21 @@ def test_NeoHookeanStress():
 
     C = 0.1 * np.eye(2) + 1
     PK2test = material.stress(C)
-    PK2good = [[-103.692, 103.356],
-               [103.356, -103.692]]
+    PK2good = [[-103.692, 103.356], [103.356, -103.692]]
     assert PK2test.ndim == 2
     assert PK2test.shape == (2, 2)
     np.testing.assert_almost_equal(PK2test, PK2good, decimal=3)
 
     C = 2 * np.eye(2) + 1
     PK2test = material.stress(C)
-    PK2good = [[12.098, -0.699],
-               [-0.699, 12.098]]
+    PK2good = [[12.098, -0.699], [-0.699, 12.098]]
     assert PK2test.ndim == 2
     assert PK2test.shape == (2, 2)
     np.testing.assert_almost_equal(PK2test, PK2good, decimal=3)
 
     C = np.eye(2) + 1
     PK2test = material.stress(C)
-    PK2good = [[8.826, 0.586],
-               [0.586, 8.826]]
+    PK2good = [[8.826, 0.586], [0.586, 8.826]]
     assert PK2test.ndim == 2
     assert PK2test.shape == (2, 2)
     np.testing.assert_almost_equal(PK2test, PK2good, decimal=3)
@@ -276,10 +278,10 @@ def test_NeoHookeanStiffness():
     material = elasticity.NeoHookeanElasticity(E, nu)
 
     C = np.eye(2)  # No deformation at all
-    Mgood = [[[[7, 0], [0, 3]],
-              [[0, 2], [2, 0]]],
-             [[[0, 2], [2, 0]],
-              [[3, 0], [0, 7]]]]
+    Mgood = [
+        [[[7, 0], [0, 3]], [[0, 2], [2, 0]]],
+        [[[0, 2], [2, 0]], [[3, 0], [0, 7]]],
+    ]
     Mgood = 5 * np.array(Mgood)
     Mtest = material.stiffness(C)
     assert Mtest.ndim == 4
@@ -287,45 +289,71 @@ def test_NeoHookeanStiffness():
     np.testing.assert_almost_equal(Mtest, Mgood)
 
     C = 0.1 * np.eye(2) + 1
-    Mgood = [[[[1602.624, -1456.93], [-1456.93, 1395.911]],
-              [[-1456.93, 1427.839], [1427.839, -1456.93]]],
-             [[[-1456.93, 1427.839], [1427.839, -1456.93]],
-              [[1395.911, -1456.93], [-1456.93, 1602.624]]]]
+    Mgood = [
+        [
+            [[1602.624, -1456.93], [-1456.93, 1395.911]],
+            [[-1456.93, 1427.839], [1427.839, -1456.93]],
+        ],
+        [
+            [[-1456.93, 1427.839], [1427.839, -1456.93]],
+            [[1395.911, -1456.93], [-1456.93, 1602.624]],
+        ],
+    ]
     Mtest = material.stiffness(C)
     assert Mtest.ndim == 4
     assert Mtest.shape == (2, 2, 2, 2)
     np.testing.assert_almost_equal(Mtest, Mgood, decimal=2)
 
     C = 2 * np.eye(2) + 1
-    Mgood = [[[[0.535, -0.178], [-0.178, 1.934]],
-              [[-0.178, -0.639], [-0.639, -0.178]]],
-             [[[-0.178, -0.639], [-0.639, -0.178]],
-              [[1.934, -0.178], [-0.178, 0.535]]]]
+    Mgood = [
+        [
+            [[0.535, -0.178], [-0.178, 1.934]],
+            [[-0.178, -0.639], [-0.639, -0.178]],
+        ],
+        [
+            [[-0.178, -0.639], [-0.639, -0.178]],
+            [[1.934, -0.178], [-0.178, 0.535]],
+        ],
+    ]
     Mtest = material.stiffness(C)
     assert Mtest.ndim == 4
     assert Mtest.shape == (2, 2, 2, 2)
     np.testing.assert_almost_equal(Mtest, Mgood, decimal=3)
 
     C = np.eye(2) + 1
-    Mgood = [[[[8.231, -4.115], [-4.115, 7.057]],
-              [[-4.115, 2.644], [2.644, -4.115]]],
-             [[[-4.115, 2.644], [2.644, -4.115]],
-              [[7.057, -4.115], [-4.115, 8.231]]]]
+    Mgood = [
+        [
+            [[8.231, -4.115], [-4.115, 7.057]],
+            [[-4.115, 2.644], [2.644, -4.115]],
+        ],
+        [
+            [[-4.115, 2.644], [2.644, -4.115]],
+            [[7.057, -4.115], [-4.115, 8.231]],
+        ],
+    ]
     Mtest = material.stiffness(C)
     assert Mtest.ndim == 4
     assert Mtest.shape == (2, 2, 2, 2)
     np.testing.assert_almost_equal(Mtest, Mgood, decimal=3)
 
     C = np.eye(3)  # No deformation at all
-    Mgood = [[[[7, 0, 0], [0, 3, 0], [0, 0, 3]],
-              [[0, 2, 0], [2, 0, 0], [0, 0, 0]],
-              [[0, 0, 2], [0, 0, 0], [2, 0, 0]]],
-             [[[0, 2, 0], [2, 0, 0], [0, 0, 0]],
-              [[3, 0, 0], [0, 7, 0], [0, 0, 3]],
-              [[0, 0, 0], [0, 0, 2], [0, 2, 0]]],
-             [[[0, 0, 2], [0, 0, 0], [2, 0, 0]],
-              [[0, 0, 0], [0, 0, 2], [0, 2, 0]],
-              [[3, 0, 0], [0, 3, 0], [0, 0, 7]]]]
+    Mgood = [
+        [
+            [[7, 0, 0], [0, 3, 0], [0, 0, 3]],
+            [[0, 2, 0], [2, 0, 0], [0, 0, 0]],
+            [[0, 0, 2], [0, 0, 0], [2, 0, 0]],
+        ],
+        [
+            [[0, 2, 0], [2, 0, 0], [0, 0, 0]],
+            [[3, 0, 0], [0, 7, 0], [0, 0, 3]],
+            [[0, 0, 0], [0, 0, 2], [0, 2, 0]],
+        ],
+        [
+            [[0, 0, 2], [0, 0, 0], [2, 0, 0]],
+            [[0, 0, 0], [0, 0, 2], [0, 2, 0]],
+            [[3, 0, 0], [0, 3, 0], [0, 0, 7]],
+        ],
+    ]
     Mgood = 5 * np.array(Mgood)
     Mtest = material.stiffness(C)
     assert Mtest.ndim == 4
@@ -333,33 +361,23 @@ def test_NeoHookeanStiffness():
     np.testing.assert_almost_equal(Mtest, Mgood)
 
     C = 2 * np.eye(3) + 0.5
-    Mgood = [[[[-0.84, 0.14, 0.14],
-               [0.14, 2.65, -0.55],
-               [0.14, -0.55, 2.65]],
-              [[0.14, -1.77, 0.33],
-               [-1.77, 0.14, 0.33],
-               [0.33, 0.33, -0.56]],
-              [[0.14, 0.33, -1.77],
-               [0.33, -0.56, 0.33],
-               [-1.77, 0.33, 0.14]]],
-             [[[0.14, -1.77, 0.33],
-               [-1.77, 0.14, 0.33],
-               [0.33, 0.33, -0.56]],
-              [[2.66, 0.14, -0.56],
-               [0.14, -0.84, 0.14],
-               [-0.56, 0.14, 2.66]],
-              [[-0.56, 0.33, 0.33],
-               [0.33, 0.14, -1.77],
-               [0.33, -1.77, 0.14]]],
-             [[[0.14, 0.33, -1.77],
-               [0.33, -0.56, 0.33],
-               [-1.77, 0.33, 0.14]],
-              [[-0.56, 0.33, 0.33],
-               [0.33, 0.14, -1.77],
-               [0.33, -1.77, 0.14]],
-              [[2.66, -0.56, 0.14],
-               [-0.56, 2.66, 0.14],
-               [0.14, 0.14, -0.84]]]]
+    Mgood = [
+        [
+            [[-0.84, 0.14, 0.14], [0.14, 2.65, -0.55], [0.14, -0.55, 2.65]],
+            [[0.14, -1.77, 0.33], [-1.77, 0.14, 0.33], [0.33, 0.33, -0.56]],
+            [[0.14, 0.33, -1.77], [0.33, -0.56, 0.33], [-1.77, 0.33, 0.14]],
+        ],
+        [
+            [[0.14, -1.77, 0.33], [-1.77, 0.14, 0.33], [0.33, 0.33, -0.56]],
+            [[2.66, 0.14, -0.56], [0.14, -0.84, 0.14], [-0.56, 0.14, 2.66]],
+            [[-0.56, 0.33, 0.33], [0.33, 0.14, -1.77], [0.33, -1.77, 0.14]],
+        ],
+        [
+            [[0.14, 0.33, -1.77], [0.33, -0.56, 0.33], [-1.77, 0.33, 0.14]],
+            [[-0.56, 0.33, 0.33], [0.33, 0.14, -1.77], [0.33, -1.77, 0.14]],
+            [[2.66, -0.56, 0.14], [-0.56, 2.66, 0.14], [0.14, 0.14, -0.84]],
+        ],
+    ]
     Mtest = material.stiffness(C)
     assert Mtest.ndim == 4
     assert Mtest.shape == (3, 3, 3, 3)
diff --git a/tests/test_fem.py b/tests/test_fem.py
index 6b741a6..9618c80 100644
--- a/tests/test_fem.py
+++ b/tests/test_fem.py
@@ -21,8 +21,9 @@
 """
 
 import numpy as np
-from hyper import gmsh2 as gmsh
+
 from hyper import fem
+from hyper import gmsh2 as gmsh
 
 
 def vector_fct(x):
@@ -63,7 +64,7 @@ def test_Tri14():
         E.append(FE_model.getStrainGreenLagrange(n).flatten())
 
     outputfilename = inputfilename.replace(".msh", "-val.msh")
-    with open(outputfilename, 'w') as outputfile:
+    with open(outputfilename, "w") as outputfile:
         gmsh.gmshOutput_mesh(outputfile, mesh)
 
         fieldname = "U: displacement field"
@@ -101,7 +102,7 @@ def test_Tri56():
         E.append(FE_model.getStrainGreenLagrange(n).flatten())
 
     outputfilename = inputfilename.replace(".msh", "-val.msh")
-    with open(outputfilename, 'w') as outputfile:
+    with open(outputfilename, "w") as outputfile:
         gmsh.gmshOutput_mesh(outputfile, mesh)
 
         fieldname = "U: displacement field"
@@ -139,7 +140,7 @@ def test_Quad11():
         E.append(FE_model.getStrainGreenLagrange(n).flatten())
 
     outputfilename = inputfilename.replace(".msh", "-val.msh")
-    with open(outputfilename, 'w') as outputfile:
+    with open(outputfilename, "w") as outputfile:
         gmsh.gmshOutput_mesh(outputfile, mesh)
 
         fieldname = "U: displacement field"
@@ -177,7 +178,7 @@ def test_Quad44():
         E.append(FE_model.getStrainGreenLagrange(n).flatten())
 
     outputfilename = inputfilename.replace(".msh", "-val.msh")
-    with open(outputfilename, 'w') as outputfile:
+    with open(outputfilename, "w") as outputfile:
         gmsh.gmshOutput_mesh(outputfile, mesh)
 
         fieldname = "U: displacement field"
diff --git a/tests/test_geometry.py b/tests/test_geometry.py
index f624229..9ae9576 100644
--- a/tests/test_geometry.py
+++ b/tests/test_geometry.py
@@ -1,7 +1,8 @@
 # -*- coding: utf-8 -*-
 
-import pytest
 import numpy as np
+import pytest
+
 from hyper import geometry
 
 
@@ -135,8 +136,14 @@ def test_derivativeSFT6():
     assert dNtest == dNgood
 
     P = (0.5, 0.5)
-    dNgood = [[1.0, 1.0], [-2.0, -2.0], [1.0, 0],
-              [2.0, 2.0], [0, 1.0], [-2.0, -2.0]]
+    dNgood = [
+        [1.0, 1.0],
+        [-2.0, -2.0],
+        [1.0, 0],
+        [2.0, 2.0],
+        [0, 1.0],
+        [-2.0, -2.0],
+    ]
     dNtest = geometry.SFT6.dN(P)
     assert dNtest == dNgood
 
@@ -146,9 +153,14 @@ def test_derivativeSFT6():
     assert dNtest == dNgood
 
     P = (1 / 3, 1 / 3)
-    dNgood = [[-1 / 3, -1 / 3], [0, -4 / 3],
-              [1 / 3, 0], [4 / 3, 4 / 3],
-              [0, 1 / 3], [-4 / 3, 0]]
+    dNgood = [
+        [-1 / 3, -1 / 3],
+        [0, -4 / 3],
+        [1 / 3, 0],
+        [4 / 3, 4 / 3],
+        [0, 1 / 3],
+        [-4 / 3, 0],
+    ]
     dNtest = geometry.SFT6.dN(P)
     np.testing.assert_almost_equal(dNtest, dNgood)
 
@@ -203,72 +215,126 @@ def test_derivativeSFQ4():
 def test_derivativeSFQ8():
     P = (-1, -1)
     dNtest = geometry.SFQ8.dN(P)
-    dNgood = [[-1.5, -1.5], [2, 0],
-              [-0.5, 0], [0, 0],
-              [0, 0], [0, 0],
-              [0, -0.5], [0, 2]]
+    dNgood = [
+        [-1.5, -1.5],
+        [2, 0],
+        [-0.5, 0],
+        [0, 0],
+        [0, 0],
+        [0, 0],
+        [0, -0.5],
+        [0, 2],
+    ]
     assert dNtest == dNgood
 
     P = (0, -1)
     dNtest = geometry.SFQ8.dN(P)
-    dNgood = [[-0.5, -0.5], [0, -0.5],
-              [0.5, -0.5], [0, 1],
-              [0, -0.5], [0, 0.5],
-              [0, -0.5], [0, 1]]
+    dNgood = [
+        [-0.5, -0.5],
+        [0, -0.5],
+        [0.5, -0.5],
+        [0, 1],
+        [0, -0.5],
+        [0, 0.5],
+        [0, -0.5],
+        [0, 1],
+    ]
     assert dNtest == dNgood
 
     P = (1, -1)
     dNtest = geometry.SFQ8.dN(P)
-    dNgood = [[0.5, 0], [-2, 0],
-              [1.5, -1.5], [0, 2],
-              [0, -0.5], [0, 0],
-              [0, 0], [0, 0]]
+    dNgood = [
+        [0.5, 0],
+        [-2, 0],
+        [1.5, -1.5],
+        [0, 2],
+        [0, -0.5],
+        [0, 0],
+        [0, 0],
+        [0, 0],
+    ]
     assert dNtest == dNgood
 
     P = (-1, 0)
     dNtest = geometry.SFQ8.dN(P)
-    dNgood = [[-0.5, -0.5], [1, 0],
-              [-0.5, 0], [0.5, 0],
-              [-0.5, 0], [1, 0],
-              [-0.5, 0.5], [-0.5, 0]]
+    dNgood = [
+        [-0.5, -0.5],
+        [1, 0],
+        [-0.5, 0],
+        [0.5, 0],
+        [-0.5, 0],
+        [1, 0],
+        [-0.5, 0.5],
+        [-0.5, 0],
+    ]
     assert dNtest == dNgood
 
     P = (0, 0)
     dNtest = geometry.SFQ8.dN(P)
-    dNgood = [[0, 0], [0, -0.5],
-              [0, 0], [0.5, 0],
-              [0, 0], [0, 0.5],
-              [0, 0], [-0.5, 0]]
+    dNgood = [
+        [0, 0],
+        [0, -0.5],
+        [0, 0],
+        [0.5, 0],
+        [0, 0],
+        [0, 0.5],
+        [0, 0],
+        [-0.5, 0],
+    ]
     assert dNtest == dNgood
 
     P = (1, 0)
     dNtest = geometry.SFQ8.dN(P)
-    dNgood = [[0.5, 0], [-1, 0],
-              [0.5, -0.5], [0.5, 0],
-              [0.5, 0.5], [-1, 0],
-              [0.5, 0], [-0.5, 0]]
+    dNgood = [
+        [0.5, 0],
+        [-1, 0],
+        [0.5, -0.5],
+        [0.5, 0],
+        [0.5, 0.5],
+        [-1, 0],
+        [0.5, 0],
+        [-0.5, 0],
+    ]
     assert dNtest == dNgood
 
     P = (-1, 1)
     dNtest = geometry.SFQ8.dN(P)
-    dNgood = [[0, 0.5], [0, 0],
-              [0, 0], [0, 0],
-              [-0.5, 0], [2, 0],
-              [-1.5, 1.5], [0, -2]]
+    dNgood = [
+        [0, 0.5],
+        [0, 0],
+        [0, 0],
+        [0, 0],
+        [-0.5, 0],
+        [2, 0],
+        [-1.5, 1.5],
+        [0, -2],
+    ]
     assert dNtest == dNgood
 
     P = (0, 1)
     dNtest = geometry.SFQ8.dN(P)
-    dNgood = [[0, 0.5], [0, -0.5],
-              [0, 0.5], [0, -1],
-              [0.5, 0.5], [0, 0.5],
-              [-0.5, 0.5], [0, -1]]
+    dNgood = [
+        [0, 0.5],
+        [0, -0.5],
+        [0, 0.5],
+        [0, -1],
+        [0.5, 0.5],
+        [0, 0.5],
+        [-0.5, 0.5],
+        [0, -1],
+    ]
     assert dNtest == dNgood
 
     P = (1, 1)
     dNtest = geometry.SFQ8.dN(P)
-    dNgood = [[0, 0], [0, 0],
-              [0, 0.5], [0, -2],
-              [1.5, 1.5], [-2, 0],
-              [0.5, 0], [0, 0]]
+    dNgood = [
+        [0, 0],
+        [0, 0],
+        [0, 0.5],
+        [0, -2],
+        [1.5, 1.5],
+        [-2, 0],
+        [0.5, 0],
+        [0, 0],
+    ]
     assert dNtest == dNgood
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index d70bbcc..2ee5e06 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -18,11 +18,12 @@
 """
 
 import numpy as np
+
 import hyper.gmsh2 as gmsh
 
 
 def scalar_fct(x):
-    f = x[0]**2 + x[1]**2
+    f = x[0] ** 2 + x[1] ** 2
     return f
 
 
@@ -107,6 +108,7 @@ def test_ReadQuad44():
 #    TESTING THE POSITION OF NODES
 #################################################
 
+
 def test_PositionNodesTri8():
     # Reading the mesh
     inputfilename = "tests/msh/meshfile-tri8.msh"
@@ -120,8 +122,10 @@ def test_PositionNodesTri8():
     for i in range(nNodes):
         Node_i = mesh.getNode(i)
         X[i] = Node_i.getX()
-    Xgood = [[0, 1, 1, 0, 0.5, 0.5, 0.25, 0.75],
-             [0, 0, 0.3, 0.3, 0, 0.3, 0.15, 0.15]]
+    Xgood = [
+        [0, 1, 1, 0, 0.5, 0.5, 0.25, 0.75],
+        [0, 0, 0.3, 0.3, 0, 0.3, 0.15, 0.15],
+    ]
     Xgood = np.array(Xgood).T
     np.testing.assert_almost_equal(X, Xgood)
 
@@ -138,8 +142,10 @@ def test_PositionNodesQuad8():
     for i in range(nNodes):
         Node_i = mesh.getNode(i)
         X[i] = Node_i.getX()
-    Xgood = [[0, 20, 20, 0, 5, 10, 15, 20, 15, 10, 5, 0, 10, 5, 15],
-             [0, 0, 6, 6, 0, 0, 0, 3, 6, 6, 6, 3, 3, 3, 3]]
+    Xgood = [
+        [0, 20, 20, 0, 5, 10, 15, 20, 15, 10, 5, 0, 10, 5, 15],
+        [0, 0, 6, 6, 0, 0, 0, 3, 6, 6, 6, 3, 3, 3, 3],
+    ]
     Xgood = np.array(Xgood).T / 20
     np.testing.assert_almost_equal(X, Xgood)
 
@@ -148,6 +154,7 @@ def test_PositionNodesQuad8():
 # #  TESTING FIELDS
 # ################################################
 
+
 def test_FieldsTri8():
     inputfilename = "tests/msh/meshfile-tri8.msh"
     with open(inputfilename, "r") as meshfile:
@@ -171,13 +178,14 @@ def test_FieldsTri8():
 
     # Writting the results in a mesh
     outputfilename = inputfilename.replace(".msh", "-val.msh")
-    with open(outputfilename, 'w') as outputfile:
+    with open(outputfilename, "w") as outputfile:
         # Write the base mesh
         gmsh.gmshOutput_mesh(outputfile, mesh)
 
         # Write scalar field
-        gmsh.gmshOutput_nodal(outputfile, "scalar field", V,
-                              it=0, t=0.0)  # write nodal field 'V'
+        gmsh.gmshOutput_nodal(
+            outputfile, "scalar field", V, it=0, t=0.0
+        )  # write nodal field 'V'
 
         # Write vectorial field
         gmsh.gmshOutput_nodal(outputfile, "VectorField", U, it=0, t=0)
@@ -206,13 +214,14 @@ def test_FieldsQuad8():
 
     # Writting the results in a mesh
     outputfilename = inputfilename.replace(".msh", "-val.msh")
-    with open(outputfilename, 'w') as outputfile:
+    with open(outputfilename, "w") as outputfile:
         # Write the base mesh
         gmsh.gmshOutput_mesh(outputfile, mesh)
 
         # Write scalar field
-        gmsh.gmshOutput_nodal(outputfile, "scalar field", V,
-                              it=0, t=0.0)  # write nodal field 'V'
+        gmsh.gmshOutput_nodal(
+            outputfile, "scalar field", V, it=0, t=0.0
+        )  # write nodal field 'V'
 
         # Write vectorial field
         gmsh.gmshOutput_nodal(outputfile, "VectorField", U, it=0, t=0)
diff --git a/tests/test_residual.py b/tests/test_residual.py
index 8260327..cb26e16 100644
--- a/tests/test_residual.py
+++ b/tests/test_residual.py
@@ -4,9 +4,9 @@
 
 
 import numpy as np
+
+from hyper import elasticity, fem
 from hyper import gmsh2 as gmsh
-from hyper import fem
-from hyper import elasticity
 
 
 def vector_fct(x):
@@ -41,24 +41,24 @@ def computeTestValues(FEmodel, U):
 
 
 def writeFieldsOnGmshFile(fields, mesh, filename):
-    with open(filename, 'w') as file:
+    with open(filename, "w") as file:
         gmsh.gmshOutput_mesh(file, mesh)
 
-        for (fieldname, field) in fields["nodal"]:
+        for fieldname, field in fields["nodal"]:
             gmsh.gmshOutput_nodal(file, fieldname, field, 0, 0.0)
 
-        for (fieldname, field) in fields["element"]:
+        for fieldname, field in fields["element"]:
             gmsh.gmshOutput_element(file, fieldname, field, 0, 0.0)
 
 
 def test_Tri56StVenantKirchhoff():
     inputfilename = "tests/msh/triangle-tri56.msh"
     outputfilename = "tests/msh/triangle-tri56-StVenant-val.msh"
-    E = 210.e9
+    E = 210.0e9
     nu = 0.3
     material = elasticity.StVenantKirchhoffElasticity(E, nu)
 
-    with open(inputfilename, 'r') as meshfile:
+    with open(inputfilename, "r") as meshfile:
         mesh = gmsh.gmshInput_mesh(meshfile)
 
     FEmodel = fem.FEModel(mesh, material)
@@ -66,21 +66,24 @@ def test_Tri56StVenantKirchhoff():
 
     R, F, P = computeTestValues(FEmodel, U)
 
-    fields = {"nodal": [("U: Displacement", U),
-                        ("R: Residual", R)],
-              "element": [("F: Deformation Gradient", F),
-                          ("P: Engineering stress", P)]}
+    fields = {
+        "nodal": [("U: Displacement", U), ("R: Residual", R)],
+        "element": [
+            ("F: Deformation Gradient", F),
+            ("P: Engineering stress", P),
+        ],
+    }
     writeFieldsOnGmshFile(fields, mesh, outputfilename)
 
 
 def test_Quad44StVenantKirchhoff():
     inputfilename = "tests/msh/triangle-quad44.msh"
     outputfilename = "tests/msh/triangle-quad44-StVenant-val.msh"
-    E = 210.e9
+    E = 210.0e9
     nu = 0.3
     material = elasticity.StVenantKirchhoffElasticity(E, nu)
 
-    with open(inputfilename, 'r') as meshfile:
+    with open(inputfilename, "r") as meshfile:
         mesh = gmsh.gmshInput_mesh(meshfile)
 
     FEmodel = fem.FEModel(mesh, material)
@@ -88,21 +91,24 @@ def test_Quad44StVenantKirchhoff():
 
     R, F, P = computeTestValues(FEmodel, U)
 
-    fields = {"nodal": [("U: Displacement", U),
-                        ("R: Residual", R)],
-              "element": [("F: Deformation Gradient", F),
-                          ("P: Engineering stress", P)]}
+    fields = {
+        "nodal": [("U: Displacement", U), ("R: Residual", R)],
+        "element": [
+            ("F: Deformation Gradient", F),
+            ("P: Engineering stress", P),
+        ],
+    }
     writeFieldsOnGmshFile(fields, mesh, outputfilename)
 
 
 def test_Tri56NeoHookean():
     inputfilename = "tests/msh/triangle-tri56.msh"
     outputfilename = "tests/msh/triangle-tri56-NeoHookean-val.msh"
-    E = 10.e6
+    E = 10.0e6
     nu = 0.45
     material = elasticity.StVenantKirchhoffElasticity(E, nu)
 
-    with open(inputfilename, 'r') as meshfile:
+    with open(inputfilename, "r") as meshfile:
         mesh = gmsh.gmshInput_mesh(meshfile)
 
     FEmodel = fem.FEModel(mesh, material)
@@ -110,21 +116,24 @@ def test_Tri56NeoHookean():
 
     R, F, P = computeTestValues(FEmodel, U)
 
-    fields = {"nodal": [("U: Displacement", U),
-                        ("R: Residual", R)],
-              "element": [("F: Deformation Gradient", F),
-                          ("P: Engineering stress", P)]}
+    fields = {
+        "nodal": [("U: Displacement", U), ("R: Residual", R)],
+        "element": [
+            ("F: Deformation Gradient", F),
+            ("P: Engineering stress", P),
+        ],
+    }
     writeFieldsOnGmshFile(fields, mesh, outputfilename)
 
 
 def test_Quad44NeoHookean():
     inputfilename = "tests/msh/triangle-tri56.msh"
     outputfilename = "tests/msh/triangle-tri56-NeoHookean-val.msh"
-    E = 10.e6
+    E = 10.0e6
     nu = 0.45
     material = elasticity.StVenantKirchhoffElasticity(E, nu)
 
-    with open(inputfilename, 'r') as meshfile:
+    with open(inputfilename, "r") as meshfile:
         mesh = gmsh.gmshInput_mesh(meshfile)
 
     FEmodel = fem.FEModel(mesh, material)
@@ -132,10 +141,13 @@ def test_Quad44NeoHookean():
 
     R, F, P = computeTestValues(FEmodel, U)
 
-    fields = {"nodal": [("U: Displacement", U),
-                        ("R: Residual", R)],
-              "element": [("F: Deformation Gradient", F),
-                          ("P: Engineering stress", P)]}
+    fields = {
+        "nodal": [("U: Displacement", U), ("R: Residual", R)],
+        "element": [
+            ("F: Deformation Gradient", F),
+            ("P: Engineering stress", P),
+        ],
+    }
     writeFieldsOnGmshFile(fields, mesh, outputfilename)
 
 
diff --git a/tests/test_tangent.py b/tests/test_tangent.py
index e290cf5..0b69005 100644
--- a/tests/test_tangent.py
+++ b/tests/test_tangent.py
@@ -5,9 +5,9 @@
 """
 
 import numpy as np
+
+from hyper import elasticity, fem
 from hyper import gmsh2 as gmsh
-from hyper import elasticity
-from hyper import fem
 
 
 def vector_fct(x):
@@ -34,7 +34,7 @@ def computeKnum(FEmodel, U, DU=1e-5):
 
     for i in range(nNodes):
         for j in range(dim):
-            Ud[i, j] += DU / 2.
+            Ud[i, j] += DU / 2.0
             # R(U_nj + DU/2), residual "plus"
             Rp = FEmodel.computeResidual(Ud)
 
@@ -50,11 +50,11 @@ def computeKnum(FEmodel, U, DU=1e-5):
 
 def test_Quad11StVenantKirchhoff():
     inputfilename = "tests/msh/triangle-quad11.msh"
-    E = 210.e9
+    E = 210.0e9
     nu = 0.3
     material = elasticity.StVenantKirchhoffElasticity(E, nu)
 
-    with open(inputfilename, 'r') as meshfile:
+    with open(inputfilename, "r") as meshfile:
         mesh = gmsh.gmshInput_mesh(meshfile)
 
     FEmodel = fem.FEModel(mesh, material)
@@ -67,11 +67,11 @@ def test_Quad11StVenantKirchhoff():
 
 def test_Quad44StVenantKirchhoff():
     inputfilename = "tests/msh/triangle-quad44.msh"
-    E = 210.e9
+    E = 210.0e9
     nu = 0.3
     material = elasticity.StVenantKirchhoffElasticity(E, nu)
 
-    with open(inputfilename, 'r') as meshfile:
+    with open(inputfilename, "r") as meshfile:
         mesh = gmsh.gmshInput_mesh(meshfile)
 
     FEmodel = fem.FEModel(mesh, material)
@@ -84,11 +84,11 @@ def test_Quad44StVenantKirchhoff():
 
 def test_Quad11NeoHookean():
     inputfilename = "tests/msh/triangle-quad11.msh"
-    E = 210.e9
+    E = 210.0e9
     nu = 0.3
     material = elasticity.StVenantKirchhoffElasticity(E, nu)
 
-    with open(inputfilename, 'r') as meshfile:
+    with open(inputfilename, "r") as meshfile:
         mesh = gmsh.gmshInput_mesh(meshfile)
 
     FEmodel = fem.FEModel(mesh, material)
@@ -101,11 +101,11 @@ def test_Quad11NeoHookean():
 
 def test_Quad44NeoHookean():
     inputfilename = "tests/msh/triangle-quad44.msh"
-    E = 210.e9
+    E = 210.0e9
     nu = 0.3
     material = elasticity.StVenantKirchhoffElasticity(E, nu)
 
-    with open(inputfilename, 'r') as meshfile:
+    with open(inputfilename, "r") as meshfile:
         mesh = gmsh.gmshInput_mesh(meshfile)
 
     FEmodel = fem.FEModel(mesh, material)
diff --git a/tests/test_tensor.py b/tests/test_tensor.py
index 5096874..a0f8d9e 100644
--- a/tests/test_tensor.py
+++ b/tests/test_tensor.py
@@ -1,5 +1,6 @@
 # -*- coding: utf-8 -*-
 import numpy as np
+
 from hyper import tensor
 
 
@@ -13,7 +14,7 @@ def test_vector():
     vectest[1] = 2.0
     assert type(vectest) is np.ndarray
     assert vectest.ndim == 1
-    assert vectest.shape == (2, )
+    assert vectest.shape == (2,)
     np.testing.assert_array_equal(vectest, [1, 2])
 
     # Testing attribution of array values
@@ -48,16 +49,13 @@ def test_tensor():
 
 def test_I():
     Itest = tensor.I()
-    Igood = [[1, 0],
-             [0, 1]]
+    Igood = [[1, 0], [0, 1]]
     assert type(Itest) is np.ndarray
     assert Itest.shape == (2, 2)
     np.testing.assert_array_equal(Itest, Igood)
 
     Itest = tensor.I(3)
-    Igood = [[1, 0, 0],
-             [0, 1, 0],
-             [0, 0, 1]]
+    Igood = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
     assert type(Itest) is np.ndarray
     assert Itest.shape == (3, 3)
     np.testing.assert_array_equal(Itest, Igood)
@@ -65,8 +63,7 @@ def test_I():
 
 def test_outerProd():
     v = [1, 2]
-    Mgood = [[1, 2],
-             [2, 4]]
+    Mgood = [[1, 2], [2, 4]]
     Mtest = tensor.outerProd(v, v)
     assert Mtest.ndim == 2
     assert Mtest.shape == (2, 2)
@@ -74,47 +71,39 @@ def test_outerProd():
 
 
 def test_det():
-    M = [[1, 2],
-         [2, 4]]
+    M = [[1, 2], [2, 4]]
     detgood = 0
     dettest = tensor.det(M)
     assert dettest == detgood
 
-    M = [[1, 2],
-         [2, 5]]
+    M = [[1, 2], [2, 5]]
     detgood = 1
     dettest = tensor.det(M)
     assert dettest == detgood
 
 
 def test_trace():
-    M = [[1, 2],
-         [2, 4]]
+    M = [[1, 2], [2, 4]]
     trgood = 5
     trtest = tensor.trace(M)
     assert trtest == trgood
 
-    M = [[1, 2],
-         [2, 5]]
+    M = [[1, 2], [2, 5]]
     trgood = 6
     trtest = tensor.trace(M)
     assert trtest == trgood
 
 
 def test_inv():
-    M = [[1, 0],
-         [0, 1]]
-    Vgood = [[1, 0],
-             [0, 1]]
+    M = [[1, 0], [0, 1]]
+    Vgood = [[1, 0], [0, 1]]
     Vtest = tensor.inv(M)
     assert Vtest.ndim == 2
     assert Vtest.shape == (2, 2)
     np.testing.assert_array_equal(Vtest, Vgood)
 
-    M = [[2, 3],
-         [2, 5]]
-    Vgood = [[1.25, -0.75],
-             [-0.5, 0.5]]
+    M = [[2, 3], [2, 5]]
+    Vgood = [[1.25, -0.75], [-0.5, 0.5]]
     Vtest = tensor.inv(M)
     assert Vtest.ndim == 2
     assert Vtest.shape == (2, 2)
@@ -122,10 +111,8 @@ def test_inv():
 
 
 def test_rightCauchyGreen():
-    F = [[2, 3],
-         [2, 5]]
-    Cgood = [[8, 16],
-             [16, 34]]
+    F = [[2, 3], [2, 5]]
+    Cgood = [[8, 16], [16, 34]]
     Ctest = tensor.rightCauchyGreen(F)
     assert Ctest.ndim == 2
     assert Ctest.shape == (2, 2)
@@ -133,10 +120,8 @@ def test_rightCauchyGreen():
 
 
 def test_leftCauchyGreen():
-    F = [[2, 3],
-         [2, 5]]
-    Bgood = [[13, 19],
-             [19, 29]]
+    F = [[2, 3], [2, 5]]
+    Bgood = [[13, 19], [19, 29]]
     Btest = tensor.leftCauchyGreen(F)
     assert Btest.ndim == 2
     assert Btest.shape == (2, 2)
@@ -156,10 +141,10 @@ def test_tensor4():
 
 
 def test_II():
-    IIgood = [[[[1, 0], [0, 0]],
-               [[0, 1], [0, 0]]],
-              [[[0, 0], [1, 0]],
-               [[0, 0], [0, 1]]]]
+    IIgood = [
+        [[[1, 0], [0, 0]], [[0, 1], [0, 0]]],
+        [[[0, 0], [1, 0]], [[0, 0], [0, 1]]],
+    ]
     IItest = tensor.II()
     assert IItest.ndim == 4
     assert IItest.shape == (2, 2, 2, 2)
@@ -168,10 +153,10 @@ def test_II():
 
 def test_IISym():
     IItest = tensor.IISym()
-    IIgood = [[[[1, 0], [0, 0]],
-               [[0, 0.5], [0.5, 0]]],
-              [[[0, 0.5], [0.5, 0]],
-               [[0, 0], [0, 1]]]]
+    IIgood = [
+        [[[1, 0], [0, 0]], [[0, 0.5], [0.5, 0]]],
+        [[[0, 0.5], [0.5, 0]], [[0, 0], [0, 1]]],
+    ]
     assert IItest.ndim == 4
     assert IItest.shape == (2, 2, 2, 2)
     np.testing.assert_array_equal(IItest, IIgood)
@@ -179,22 +164,21 @@ def test_IISym():
 
 def test_KK():
     KKtest = tensor.KK()
-    KKgood = [[[[1, 0], [0, 1]],
-               [[0, 0], [0, 0]]],
-              [[[0, 0], [0, 0]],
-               [[1, 0], [0, 1]]]]
+    KKgood = [
+        [[[1, 0], [0, 1]], [[0, 0], [0, 0]]],
+        [[[0, 0], [0, 0]], [[1, 0], [0, 1]]],
+    ]
     assert KKtest.ndim == 4
     assert KKtest.shape == (2, 2, 2, 2)
     np.testing.assert_array_equal(KKtest, KKgood)
 
 
 def test_outerProd4():
-    C = [[8, 16],
-         [16, 34]]
-    Mgood = [[[[64, 128], [128, 272]],
-              [[128, 256], [256, 544]]],
-             [[[128, 256], [256, 544]],
-              [[272, 544], [544, 1156]]]]
+    C = [[8, 16], [16, 34]]
+    Mgood = [
+        [[[64, 128], [128, 272]], [[128, 256], [256, 544]]],
+        [[[128, 256], [256, 544]], [[272, 544], [544, 1156]]],
+    ]
     Mtest = tensor.outerProd4(C, C)
     assert Mtest.ndim == 4
     assert Mtest.shape == (2, 2, 2, 2)