diff --git a/conmech/dynamics/dynamics.py b/conmech/dynamics/dynamics.py
index a7e0a3b6..3b1953cc 100644
--- a/conmech/dynamics/dynamics.py
+++ b/conmech/dynamics/dynamics.py
@@ -115,7 +115,7 @@ def reinitialize_matrices(self, elements_density: Optional[np.ndarray] = None):
             U,
             V,
             self._w_matrix,
-            self._local_stifness_matrices
+            self._local_stifness_matrices,
         ) = get_basic_matrices(elements=self.mesh.elements, nodes=self.moved_nodes)
 
         if elements_density is not None:
@@ -125,7 +125,10 @@ def reinitialize_matrices(self, elements_density: Optional[np.ndarray] = None):
                     for j, global_j in enumerate(element):
                         for x in (0, 1):
                             for y in (0, 1):
-                                self._w_matrix[x][y][global_i, global_j] += elements_density[element_index] * self._local_stifness_matrices[2*x + y, element_index, i, j]
+                                self._w_matrix[x][y][global_i, global_j] += (
+                                    elements_density[element_index]
+                                    * self._local_stifness_matrices[2 * x + y, element_index, i, j]
+                                )
 
         (
             self.acceleration_operator,
diff --git a/conmech/dynamics/factory/_dynamics_factory_2d.py b/conmech/dynamics/factory/_dynamics_factory_2d.py
index 03ed9560..97f17075 100644
--- a/conmech/dynamics/factory/_dynamics_factory_2d.py
+++ b/conmech/dynamics/factory/_dynamics_factory_2d.py
@@ -87,8 +87,10 @@ def get_edges_features_matrix_numba(elements, nodes):
 
                 w = [[i_d_phi * j_d_phi for j_d_phi in j_d_phi_vec] for i_d_phi in i_d_phi_vec]
 
-                local_stifness_matrices[:, element_index, i, j] = element_volume * np.asarray(w).flatten()
-                
+                local_stifness_matrices[:, element_index, i, j] = (
+                    element_volume * np.asarray(w).flatten()
+                )
+
                 edges_features_matrix[:, element[i], element[j]] += element_volume * np.array(
                     [
                         volume_at_nodes,
diff --git a/conmech/dynamics/factory/dynamics_factory_method.py b/conmech/dynamics/factory/dynamics_factory_method.py
index 4a5684bb..59729bc4 100644
--- a/conmech/dynamics/factory/dynamics_factory_method.py
+++ b/conmech/dynamics/factory/dynamics_factory_method.py
@@ -26,9 +26,11 @@ def get_basic_matrices(elements: np.ndarray, nodes: np.ndarray):
     dimension = len(elements[0]) - 1
     factory = get_factory(dimension)
 
-    edges_features_matrix, element_initial_volume, local_stifness_matrices = factory.get_edges_features_matrix(
-        elements, nodes
-    )
+    (
+        edges_features_matrix,
+        element_initial_volume,
+        local_stifness_matrices,
+    ) = factory.get_edges_features_matrix(elements, nodes)
 
     volume_at_nodes = edges_features_matrix[0]
     U = edges_features_matrix[1]
diff --git a/conmech/simulations/problem_solver.py b/conmech/simulations/problem_solver.py
index 3486a359..81ee099b 100644
--- a/conmech/simulations/problem_solver.py
+++ b/conmech/simulations/problem_solver.py
@@ -381,7 +381,6 @@ def solve(self, *, initial_displacement: Callable, verbose: bool = False, **kwar
 
 
 class NonHomogenousSolver(StaticSolver):
-    
     def update_density(self, density: np.ndarray):
         self.body.reinitialize_matrices(density)
         self.refresh_solvers()
diff --git a/examples/example_density.py b/examples/example_density.py
index 5760ea68..d8f35e25 100644
--- a/examples/example_density.py
+++ b/examples/example_density.py
@@ -34,8 +34,7 @@ def friction_bound(u_nu: float) -> float:
         return 0
 
     boundaries: ... = BoundariesDescription(
-        contact=lambda x: x[1] == 0 and x[0] < 0.2, 
-        dirichlet=lambda x: x[0] == 0
+        contact=lambda x: x[1] == 0 and x[0] < 0.2, dirichlet=lambda x: x[0] == 0
     )