alder version of meshzoo

conmech/dynamics/contact/contact_law.py

@@ -37,9 +37,7 @@ def normal_bound(var_nu: float, static_displacement_nu: float, dt: float) -> flo
         return 1.0
 
     @staticmethod
-    def tangential_bound(
-        var_nu: float, static_displacement_nu: float, dt: float
-    ) -> float:
+    def tangential_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float:
         """
         Friction bound
 

conmech/dynamics/contact/damped_normal_compliance.py

@@ -20,16 +20,12 @@
 from conmech.dynamics.contact.interior_contact_law import InteriorContactLaw
 
 
-def make_damped_norm_compl(
-    obstacle_level: float, kappa: float, beta: float, interior=False
-):
+def make_damped_norm_compl(obstacle_level: float, kappa: float, beta: float, interior=False):
     superclass = InteriorContactLaw if interior else PotentialOfContactLaw
 
     class DampedNormalCompliance(superclass):
         @staticmethod
-        def normal_bound(
-            var_nu: float, static_displacement_nu: float, dt: float
-        ) -> float:
+        def normal_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float:
             """
             Since multiply by var_nu
             """

conmech/dynamics/factory/_abstract_dynamics_factory.py

@@ -23,14 +23,10 @@ def get_relaxation_tensor(self, W: np.ndarray, coeff: np.ndarray):
     def calculate_acceleration(self, U: np.ndarray, density: float) -> np.ndarray:
         raise NotImplementedError()
 
-    def calculate_thermal_expansion(
-        self, V: np.ndarray, coeff: np.ndarray
-    ) -> np.ndarray:
+    def calculate_thermal_expansion(self, V: np.ndarray, coeff: np.ndarray) -> np.ndarray:
         raise NotImplementedError()
 
-    def calculate_thermal_conductivity(
-        self, W: np.ndarray, coeff: np.ndarray
-    ) -> np.ndarray:
+    def calculate_thermal_conductivity(self, W: np.ndarray, coeff: np.ndarray) -> np.ndarray:
         raise NotImplementedError()
 
     def get_piezoelectric_tensor(self, W: np.ndarray, coeff: np.ndarray) -> np.ndarray:

conmech/dynamics/factory/_dynamics_factory_2d.py

@@ -43,9 +43,7 @@ def get_edges_features_matrix_numba(elements, nodes):
         )
         en0[:] = element_nodes[:, 0]
         en1[:] = element_nodes[:, 1]
-        int_xy = (
-            en0 @ np.array([[2.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 2.0]]) @ en1
-        )
+        int_xy = en0 @ np.array([[2.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 2.0]]) @ en1
         int_matrix = np.array(
             [
                 [
@@ -103,18 +101,11 @@ def get_edges_features_matrix_numba(elements, nodes):
 
                 v = [INT_PH * j_d_phi for j_d_phi in j_d_phi_vec]
 
-                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
-                ]
+                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)
-                )
+                local_stifness_matrices[:, :, element_index, i, j] = element_volume * np.asarray(w)
 
-                edges_features_matrix[
-                    :, element[i], element[j]
-                ] += element_volume * np.array(
+                edges_features_matrix[:, element[i], element[j]] += element_volume * np.array(
                     [
                         volume_at_nodes,
                         u,
@@ -169,9 +160,7 @@ def get_edges_features_matrix(self, elements, nodes):
     def dimension(self) -> int:
         return DIMENSION
 
-    def calculate_constitutive_matrices(
-        self, W: np.ndarray, mu: float, lambda_: float
-    ) -> SM2:
+    def calculate_constitutive_matrices(self, W: np.ndarray, mu: float, lambda_: float) -> SM2:
         A_11 = (2 * mu + lambda_) * W[0, 0] + mu * W[1, 1]
         A_12 = mu * W[1, 0] + lambda_ * W[0, 1]
         A_21 = lambda_ * W[1, 0] + mu * W[0, 1]

conmech/dynamics/factory/_dynamics_factory_3d.py

@@ -48,14 +48,9 @@ def get_edges_features_matrix_numba(elements, nodes):
 
                 v = [INT_PH * j_d_phi for j_d_phi in j_d_phi_vec]
 
-                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
-                ]
+                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]
 
-                edges_features_matrix[
-                    :, element[i], element[j]
-                ] += element_volume * np.array(
+                edges_features_matrix[:, element[i], element[j]] += element_volume * np.array(
                     [
                         volume_at_nodes,
                         u,
@@ -81,9 +76,7 @@ def get_edges_features_matrix_numba(elements, nodes):
 @numba.njit
 def get_integral_parts_numba(element_nodes, element_index):
     x_i = element_nodes[element_index]
-    x_j1, x_j2, x_j3 = list(
-        element_nodes[np.arange(ELEMENT_NODES_COUNT) != element_index]
-    )
+    x_j1, x_j2, x_j3 = list(element_nodes[np.arange(ELEMENT_NODES_COUNT) != element_index])
 
     dm = denominator_numba(x_i, x_j1, x_j2, x_j3)
     element_volume = np.abs(dm) / VOLUME_DIVIDER

conmech/dynamics/factory/dynamics_factory_method.py

@@ -60,23 +60,17 @@ def get_dynamics(elements: np.ndarray, body_prop: BodyProperties, U, V, W):
         poisson_operator = None
 
     if isinstance(body_prop, ElasticProperties):
-        elasticity = factory.calculate_constitutive_matrices(
-            W, body_prop.mu, body_prop.lambda_
-        )
+        elasticity = factory.calculate_constitutive_matrices(W, body_prop.mu, body_prop.lambda_)
     else:
         elasticity = None
 
     if isinstance(body_prop, ViscoelasticProperties):
-        viscosity = factory.calculate_constitutive_matrices(
-            W, body_prop.theta, body_prop.zeta
-        )
+        viscosity = factory.calculate_constitutive_matrices(W, body_prop.theta, body_prop.zeta)
     else:
         viscosity = None
 
     if isinstance(body_prop, TemperatureBodyProperties):
-        thermal_expansion = factory.calculate_thermal_expansion(
-            V, body_prop.thermal_expansion
-        )
+        thermal_expansion = factory.calculate_thermal_expansion(V, body_prop.thermal_expansion)
         thermal_conductivity = factory.calculate_thermal_conductivity(
             W, body_prop.thermal_conductivity
         )
@@ -85,9 +79,7 @@ def get_dynamics(elements: np.ndarray, body_prop: BodyProperties, U, V, W):
         thermal_conductivity = None
 
     if isinstance(body_prop, PiezoelectricBodyProperties):
-        piezoelectricity = factory.get_piezoelectric_tensor(
-            W, body_prop.piezoelectricity
-        )
+        piezoelectricity = factory.get_piezoelectric_tensor(W, body_prop.piezoelectricity)
         permittivity = factory.get_permittivity_tensor(W, body_prop.permittivity)
     else:
         piezoelectricity = None

conmech/dynamics/statement.py

@@ -99,11 +99,7 @@ def update_right_hand_side(self, var):
 
         A = -1 * self.body.dynamics.poisson_operator @ var.displacement[:ind]
 
-        A += (
-            (1 / var.time_step)
-            * self.body.dynamics.acceleration_operator.SM1
-            @ var.velocity[:ind]
-        )
+        A += (1 / var.time_step) * self.body.dynamics.acceleration_operator.SM1 @ var.velocity[:ind]
 
         self.right_hand_side = self.body.dynamics.force.integrate(time=var.time) + A
 
@@ -183,11 +179,7 @@ def update_right_hand_side(self, var):
 
         A = -1 * self.body.dynamics.elasticity @ var.displacement
 
-        A += (
-            (1 / var.time_step)
-            * self.body.dynamics.acceleration_operator
-            @ var.velocity
-        )
+        A += (1 / var.time_step) * self.body.dynamics.acceleration_operator @ var.velocity
 
         self.right_hand_side = self.body.dynamics.force.integrate(time=var.time) + A
 

conmech/helpers/nph.py

@@ -39,6 +39,5 @@ def get_node_index_numba(node, nodes):
 @numba.njit(inline="always")
 def length(edge, nodes):
     return np.sqrt(
-        (nodes[edge[0]][0] - nodes[edge[1]][0]) ** 2
-        + (nodes[edge[0]][1] - nodes[edge[1]][1]) ** 2
+        (nodes[edge[0]][0] - nodes[edge[1]][0]) ** 2 + (nodes[edge[0]][1] - nodes[edge[1]][1]) ** 2
     )

conmech/mesh/boundaries.py

@@ -61,9 +61,7 @@ def get_boundary_surfaces_normals(nodes, boundary_surfaces, boundary_internal_in
 
     internal_nodes = nodes[boundary_internal_indices]
     external_orientation = (-1) * np.sign(
-        nph.elementwise_dot(
-            internal_nodes - tail_nodes, unoriented_normals, keepdims=True
-        )
+        nph.elementwise_dot(internal_nodes - tail_nodes, unoriented_normals, keepdims=True)
     )
     return unoriented_normals * external_orientation
 
@@ -86,19 +84,13 @@ def __init__(self, nodes, boundary_internal_indices: np.ndarray, **kwargs):
     @property
     def boundary_surfaces(self):
         return np.unique(
-            np.vstack(
-                (self.contact_boundary, self.neumann_boundary, self.dirichlet_boundary)
-            ),
+            np.vstack((self.contact_boundary, self.neumann_boundary, self.dirichlet_boundary)),
             axis=1,
         )
 
     @property
     def boundary_nodes_count(self):
-        return (
-            self.contact_nodes_count
-            + self.neumann_nodes_count
-            + self.dirichlet_nodes_count
-        )
+        return self.contact_nodes_count + self.neumann_nodes_count + self.dirichlet_nodes_count
 
     @property
     def boundary_indices(self):
@@ -139,26 +131,20 @@ def get_all_boundary_indices(self, boundary, total_node_count, dimension):
                 condition_start = 0
                 condition_stop = condition_start + i.stop - i.start
                 yield (
-                    slice(
-                        i.start + d * total_node_count, i.stop + d * total_node_count
-                    ),
+                    slice(i.start + d * total_node_count, i.stop + d * total_node_count),
                     slice(condition_start, condition_stop),
                 )
         else:
             # assuming node_indices are sorted
-            discontinuities = np.concatenate(
-                ([0], np.nonzero(np.diff(i) - 1)[0] + 1, [len(i)])
-            )
+            discontinuities = np.concatenate(([0], np.nonzero(np.diff(i) - 1)[0] + 1, [len(i)]))
             starts = i[discontinuities[:-1]]
             stops = i[(discontinuities - 1)[1:]] + 1
             for d in range(dimension):
                 condition_start = 0
                 for start, stop in zip(starts, stops):
                     condition_stop = condition_start + stop - start
                     yield (
-                        slice(
-                            start + d * total_node_count, stop + d * total_node_count
-                        ),
+                        slice(start + d * total_node_count, stop + d * total_node_count),
                         slice(condition_start, condition_stop),
                     )
                     condition_start = condition_stop

conmech/mesh/boundaries_factory.py

@@ -21,9 +21,7 @@ def identify_surfaces_numba(sorted_elements):
     for j in range(element_size):
         # exclude each node from sorted elements and get all combinations to obtain surfaces
         surfaces[i : i + elements_count, :j] = sorted_elements[:, :j]
-        surfaces[i : i + elements_count, j:dim] = sorted_elements[
-            :, j + 1 : element_size
-        ]
+        surfaces[i : i + elements_count, j:dim] = sorted_elements[:, j + 1 : element_size]
         opposing_indices[i : i + elements_count] = sorted_elements[:, j]
         i += elements_count
     return surfaces, opposing_indices
@@ -45,17 +43,13 @@ def extract_unique_indices(surfaces):
 
 
 def extract_unique_elements(elements, opposing_indices):
-    _, indices, count = np.unique(
-        elements, axis=0, return_index=True, return_counts=True
-    )
+    _, indices, count = np.unique(elements, axis=0, return_index=True, return_counts=True)
     unique_indices = indices[count == 1]
     return elements[unique_indices], opposing_indices[unique_indices]
 
 
 def apply_predicate_to_surfaces(surfaces, nodes, predicate: Callable):
-    mask = [
-        predicate(m) for m in np.mean(nodes[surfaces], axis=1)
-    ]  # TODO: #65 Use numba (?)
+    mask = [predicate(m) for m in np.mean(nodes[surfaces], axis=1)]  # TODO: #65 Use numba (?)
     return surfaces[mask]
 
 
@@ -71,9 +65,7 @@ def reorder_boundary_nodes(nodes, elements, is_contact, is_dirichlet):
     Dirichlet override other
     """
     # move boundary nodes to the top
-    nodes, elements, boundary_nodes_count = reorder(
-        nodes, elements, lambda _: True, to_top=True
-    )
+    nodes, elements, boundary_nodes_count = reorder(nodes, elements, lambda _: True, to_top=True)
     # then move contact nodes to the top
     nodes, elements, contact_nodes_count = reorder(
         nodes,
@@ -82,9 +74,7 @@ def reorder_boundary_nodes(nodes, elements, is_contact, is_dirichlet):
         to_top=True,
     )
     # finally move dirichlet nodes to the bottom
-    nodes, elements, dirichlet_nodes_count = reorder(
-        nodes, elements, is_dirichlet, to_top=False
-    )
+    nodes, elements, dirichlet_nodes_count = reorder(nodes, elements, is_dirichlet, to_top=False)
     return (
         nodes,
         elements,
@@ -173,17 +163,11 @@ def identify_boundaries_and_reorder_nodes(
             is_dirichlet=is_dirichlet,
         )
 
-        neumann_nodes_count = (
-            boundary_nodes_count - contact_nodes_count - dirichlet_nodes_count
-        )
+        neumann_nodes_count = boundary_nodes_count - contact_nodes_count - dirichlet_nodes_count
 
-        (boundary_surfaces, boundary_internal_indices, *_) = get_boundary_surfaces(
-            elements
-        )
+        (boundary_surfaces, boundary_internal_indices, *_) = get_boundary_surfaces(elements)
 
-        contact_boundary = apply_predicate_to_surfaces(
-            boundary_surfaces, initial_nodes, is_contact
-        )
+        contact_boundary = apply_predicate_to_surfaces(boundary_surfaces, initial_nodes, is_contact)
         dirichlet_boundary = apply_predicate_to_surfaces(
             boundary_surfaces, initial_nodes, is_dirichlet
         )
@@ -200,14 +184,10 @@ def identify_boundaries_and_reorder_nodes(
         )
         neumann_boundary = Boundary(
             surfaces=neumann_boundary,
-            node_indices=slice(
-                contact_nodes_count, contact_nodes_count + neumann_nodes_count
-            ),
+            node_indices=slice(contact_nodes_count, contact_nodes_count + neumann_nodes_count),
             node_count=neumann_nodes_count,
         )
-        dirichlet_indices = slice(
-            len(initial_nodes) - dirichlet_nodes_count, len(initial_nodes)
-        )
+        dirichlet_indices = slice(len(initial_nodes) - dirichlet_nodes_count, len(initial_nodes))
         dirichlet_boundary = Boundary(
             surfaces=dirichlet_boundary,
             node_indices=dirichlet_indices,
@@ -220,9 +200,7 @@ def identify_boundaries_and_reorder_nodes(
         other_boundaries = {}
         for name, indicator in boundaries_description.indicators.items():
             if name not in ("contact", "dirichlet"):
-                surfaces = apply_predicate_to_surfaces(
-                    boundary_surfaces, initial_nodes, indicator
-                )
+                surfaces = apply_predicate_to_surfaces(boundary_surfaces, initial_nodes, indicator)
                 node_indices = np.unique(surfaces)
                 if len(node_indices):
                     other_boundaries[name] = Boundary(

conmech/mesh/mesh.py

@@ -105,9 +105,7 @@ def __init__(
         ) = BoundariesFactory.identify_boundaries_and_reorder_nodes(
             unordered_nodes, unordered_elements, boundaries_description
         )
-        edges_matrix = get_edges_matrix(
-            nodes_count=len(self.nodes), elements=self.elements
-        )
+        edges_matrix = get_edges_matrix(nodes_count=len(self.nodes), elements=self.elements)
         self.edges = get_edges_list_numba(edges_matrix)
 
     @property
@@ -173,11 +171,7 @@ def independent_nodes_count(self):
     @property
     def free_nodes_count(self):
         # TODO: #65 CHECK
-        return (
-            self.independent_nodes_count
-            - self.contact_nodes_count
-            - self.dirichlet_nodes_count
-        )
+        return self.independent_nodes_count - self.contact_nodes_count - self.dirichlet_nodes_count
 
     @property
     def boundary_indices(self):

conmech/mesh/zoo/barboteu_2008.py

@@ -44,9 +44,7 @@ def __init__(self, mesh_descr: Barboteu2008MeshDescription):
 
     @staticmethod
     def _set_mesh_size(geom, mesh_descr: Barboteu2008MeshDescription):
-        geom.set_mesh_size_callback(
-            lambda dim, tag, x, y, z, *_: mesh_descr.max_element_perimeter
-        )
+        geom.set_mesh_size_callback(lambda dim, tag, x, y, z, *_: mesh_descr.max_element_perimeter)
 
     @staticmethod
     def _get_nodes_and_elements(geom, dim):