Merge branch 'update-0.3.3' into 359_dev_mesh_rigid_body

.flake8

@@ -1,5 +1,5 @@
 [flake8]
-ignore = E203, E266, E501, W503, F403, F401, W191
+ignore = E203, E266, E501, W503, F403, F401, W191, E704
 max-line-length = 88
 max-complexity = 18
 select = B,C,E,F,W,T4,B9

elastica/_rotations.py

@@ -161,7 +161,7 @@ def _inv_rotate(director_collection):
 def _generate_skew_map(dim: int):
     # TODO Documentation
     # Preallocate
-    mapping_list = [None] * ((dim ** 2 - dim) // 2)
+    mapping_list = [None] * ((dim**2 - dim) // 2)
     # Indexing (i,j), j is the fastest changing
     # r = 2, r here is rank, we deal with only matrices
     for index, (i, j) in enumerate(combinations(range(dim), r=2)):

elastica/collision/AABBCollection.py

@@ -1,5 +1,6 @@
 """ Axis Aligned Bounding Boxes for coarse collision detection
 """
+
 import numpy as np
 from elastica.utils import MaxDimension
 
@@ -129,8 +130,8 @@ def __init__(
         # Check which is the closest and use that as the number of levels
 
         # Else check differences
-        if (4 ** n_levels_bound_above - potential_n_aabbs_in_final_level) > (
-            potential_n_aabbs_in_final_level - 4 ** n_levels_bound_below
+        if (4**n_levels_bound_above - potential_n_aabbs_in_final_level) > (
+            potential_n_aabbs_in_final_level - 4**n_levels_bound_below
         ):
             self.n_levels = n_levels_bound_below + 1
         else:

elastica/contact_utils.py

@@ -46,7 +46,7 @@ def _find_min_dist(x1, e1, x2, e2):
     s = 0.0
     t = 0.0
 
-    parallel = abs(1.0 - e1e2 ** 2 / (e1e1 * e2e2)) < 1e-6
+    parallel = abs(1.0 - e1e2**2 / (e1e1 * e2e2)) < 1e-6
     if parallel:
         # Some are parallel, so do processing
         t = (x2e1 - x1e1) / e1e1  # Comes from taking dot of e1 with a normal

elastica/experimental/connection_contact_joint/parallel_connection.py

@@ -97,7 +97,11 @@ def __init__(
     def apply_forces(self, rod_one, index_one, rod_two, index_two):
         # TODO: documentation
 
-        (self.rod_one_rd2, self.rod_two_rd2, self.spring_force,) = self._apply_forces(
+        (
+            self.rod_one_rd2,
+            self.rod_two_rd2,
+            self.spring_force,
+        ) = self._apply_forces(
             self.k,
             self.nu,
             self.k_repulsive,

elastica/memory_block/memory_block_rod.py

@@ -146,16 +146,16 @@ def __init__(self, systems: Sequence, system_idx_list):
             self.start_idx_in_rod_nodes[n_straight_rods:] = (
                 self.periodic_boundary_nodes_idx[0, 0::3] + 1
             )
-            self.end_idx_in_rod_nodes[
-                n_straight_rods:
-            ] = self.periodic_boundary_nodes_idx[0, 1::3]
+            self.end_idx_in_rod_nodes[n_straight_rods:] = (
+                self.periodic_boundary_nodes_idx[0, 1::3]
+            )
 
             self.start_idx_in_rod_elems[n_straight_rods:] = (
                 self.periodic_boundary_elems_idx[0, 0::2] + 1
             )
-            self.end_idx_in_rod_elems[
-                n_straight_rods:
-            ] = self.periodic_boundary_elems_idx[0, 1::2]
+            self.end_idx_in_rod_elems[n_straight_rods:] = (
+                self.periodic_boundary_elems_idx[0, 1::2]
+            )
 
             self.start_idx_in_rod_voronoi[n_straight_rods:] = (
                 self.periodic_boundary_voronoi_idx[0, :] + 1

elastica/memory_block/memory_block_rod_base.py

@@ -156,7 +156,6 @@ def make_block_memory_periodic_boundary_metadata(n_elems_in_rods):
 
 
 class MemoryBlockRodBase:
-
     """
     This is the base class for memory blocks for rods.
     """

elastica/modules/base_system.py

@@ -48,9 +48,9 @@ def __init__(self):
         self._feature_group_synchronize: Iterable[Callable[[float], None]] = []
         self._feature_group_constrain_values: Iterable[Callable[[float], None]] = []
         self._feature_group_constrain_rates: Iterable[Callable[[float], None]] = []
-        self._feature_group_callback: Iterable[
-            Callable[[float, int, AnyStr], None]
-        ] = []
+        self._feature_group_callback: Iterable[Callable[[float, int, AnyStr], None]] = (
+            []
+        )
         self._feature_group_finalize: Iterable[Callable] = []
         # We need to initialize our mixin classes
         super(BaseSystemCollection, self).__init__()

elastica/rigidbody/sphere.py

@@ -26,7 +26,7 @@ def __init__(self, center, base_radius, density):
         self.density = density
         self.length = 2 * base_radius
         # This is for a rigid body cylinder
-        self.volume = 4.0 / 3.0 * np.pi * base_radius ** 3
+        self.volume = 4.0 / 3.0 * np.pi * base_radius**3
         self.mass = np.array([self.volume * self.density])
         normal = np.array([1.0, 0.0, 0.0]).reshape(3, 1)
         tangents = np.array([0.0, 0.0, 1.0]).reshape(3, 1)
@@ -37,7 +37,7 @@ def __init__(self, center, base_radius, density):
             (MaxDimension.value(), MaxDimension.value()), np.float64
         )
         np.fill_diagonal(
-            mass_second_moment_of_inertia, 2.0 / 5.0 * self.mass * self.radius ** 2
+            mass_second_moment_of_inertia, 2.0 / 5.0 * self.mass * self.radius**2
         )
 
         self.mass_second_moment_of_inertia = mass_second_moment_of_inertia.reshape(

elastica/rod/cosserat_rod.py

@@ -991,7 +991,7 @@ def _compute_internal_torques(
     )
 
     # FIXME: change memory overload instead for the below calls!
-    voronoi_dilatation_inv_cube_cached = 1.0 / voronoi_dilatation ** 3
+    voronoi_dilatation_inv_cube_cached = 1.0 / voronoi_dilatation**3
     # Delta(\tau_L / \Epsilon^3)
     bend_twist_couple_2D = difference_kernel_for_block_structure(
         internal_couple * voronoi_dilatation_inv_cube_cached, ghost_voronoi_idx

elastica/rod/factory_function.py

@@ -212,7 +212,7 @@ def allocate(
     ) / (rest_lengths_temp_for_voronoi[1:] + rest_lengths_temp_for_voronoi[:-1])
 
     # Compute volume of elements
-    volume = np.pi * radius ** 2 * rest_lengths
+    volume = np.pi * radius**2 * rest_lengths
 
     # Compute mass of elements
     mass = np.zeros(n_nodes)

elastica/rod/knot_theory.py

@@ -35,20 +35,16 @@ class KnotTheoryCompatibleProtocol(Protocol):
     """
 
     @property
-    def position_collection(self) -> np.ndarray:
-        ...
+    def position_collection(self) -> np.ndarray: ...
 
     @property
-    def director_collection(self) -> np.ndarray:
-        ...
+    def director_collection(self) -> np.ndarray: ...
 
     @property
-    def radius(self) -> np.ndarray:
-        ...
+    def radius(self) -> np.ndarray: ...
 
     @property
-    def base_length(self) -> np.ndarray:
-        ...
+    def base_length(self) -> np.ndarray: ...
 
 
 class KnotTheory:

elastica/systems/analytical.py

@@ -98,7 +98,7 @@ def __init__(self, omega=2.0 * np.pi, init_val=np.array([1.0, 0.0])):
         self.omega = omega
         self.initial_value = init_val.copy()
         self._state = init_val.copy()
-        self.A_matrix = np.array([[0.0, 1.0], [-self.omega ** 2, 0.0]])
+        self.A_matrix = np.array([[0.0, 1.0], [-self.omega**2, 0.0]])
 
     def analytical_solution(self, time):
         # http://scipp.ucsc.edu/~haber/ph5B/sho09.pdf
@@ -155,7 +155,7 @@ def compute_energy(self, time):
         analytical_state = self.analytical_solution(time)
 
         def energy(st):
-            return self.omega ** 2 * st[0] ** 2 + st[1] ** 2
+            return self.omega**2 * st[0] ** 2 + st[1] ** 2
 
         anal_energy = energy(analytical_state)
         current_energy = energy(self._state)
@@ -174,10 +174,8 @@ class DampedSimpleHarmonicOscillatorSystem(
     def __init__(self, damping=0.5, omega=2.0 * np.pi, init_val=np.array([1.0, 0.0])):
         super(DampedSimpleHarmonicOscillatorSystem, self).__init__(omega, init_val)
         self.zeta = 0.5 * damping
-        self.modified_omega = np.sqrt(
-            self.omega ** 2 - self.zeta ** 2, dtype=np.complex128
-        )
-        self.A_matrix = np.array([[0.0, 1.0], [-self.omega ** 2, -damping]])
+        self.modified_omega = np.sqrt(self.omega**2 - self.zeta**2, dtype=np.complex128)
+        self.A_matrix = np.array([[0.0, 1.0], [-self.omega**2, -damping]])
 
     def analytical_solution(self, time):
         # https://en.wikipedia.org/wiki/Harmonic_oscillator#Transient_solution
@@ -420,15 +418,15 @@ def analytical_solution(self, type, time):
         if type == "Positions":
             analytical_solution = (
                 np.hstack((self.init_pos, self.final_pos))
-                + np.sin(np.pi * time) / np.pi ** 2
+                + np.sin(np.pi * time) / np.pi**2
             )
         elif type == "Velocity":
             analytical_solution = (
                 0.0 * np.hstack((self.init_pos, self.final_pos))
                 + np.cos(np.pi * time) / np.pi
             )
         elif type == "Directors":
-            final_angle = self.omega_value * time + 0.5 * 0.1 * np.pi * time ** 2
+            final_angle = self.omega_value * time + 0.5 * 0.1 * np.pi * time**2
             axis = np.array([0.0, 0.0, 1.0]).reshape(3, 1)  # There is only one director
             # Reshaping done to prevent numba equivalent to complain
             # While we can prevent it here, its' done to make the front end testing scripts "look"

elastica/timestepper/symplectic_steppers.py

@@ -31,7 +31,7 @@ class _SystemInstanceStepper:
     def do_step(
         TimeStepper, _steps_and_prefactors, System, time: np.float64, dt: np.float64
     ):
-        for (kin_prefactor, kin_step, dyn_step) in _steps_and_prefactors[:-1]:
+        for kin_prefactor, kin_step, dyn_step in _steps_and_prefactors[:-1]:
             kin_step(TimeStepper, System, time, dt)
             time += kin_prefactor(TimeStepper, dt)
             System.update_internal_forces_and_torques(time)
@@ -71,7 +71,7 @@ def do_step(
         -------
 
         """
-        for (kin_prefactor, kin_step, dyn_step) in _steps_and_prefactors[:-1]:
+        for kin_prefactor, kin_step, dyn_step in _steps_and_prefactors[:-1]:
 
             for system in SystemCollection._memory_blocks:
                 kin_step(TimeStepper, system, time, dt)

elastica/transformations.py

@@ -79,7 +79,7 @@ def format_matrix_shape(matrix_collection):
     # we need to convert the matrix in that case.
     def assert_proper_square(num1):
         sqrt_num = isqrt(num1)
-        assert sqrt_num ** 2 == num1, "Matrix dimension passed is not a perfect square"
+        assert sqrt_num**2 == num1, "Matrix dimension passed is not a perfect square"
         return sqrt_num
 
     if n_dim == 1:

examples/AxialStretchingCase/axial_stretching.py

@@ -23,6 +23,7 @@
 
     isort:skip_file
 """
+
 import numpy as np
 from matplotlib import pyplot as plt
 
@@ -50,7 +51,7 @@ class StretchingBeamSimulator(
 normal = np.array([0.0, 1.0, 0.0])
 base_length = 1.0
 base_radius = 0.025
-base_area = np.pi * base_radius ** 2
+base_area = np.pi * base_radius**2
 density = 1000
 youngs_modulus = 1e4
 # For shear modulus of 1e4, nu is 99!
@@ -90,6 +91,7 @@ class StretchingBeamSimulator(
     time_step=dt,
 )
 
+
 # Add call backs
 class AxialStretchingCallBack(ea.CallBackBaseClass):
     """

examples/BoundaryConditionsCases/general_constraint_allow_yaw.py

@@ -31,7 +31,7 @@ class GeneralConstraintSimulator(
 normal = np.array([0.0, 1.0, 0.0])
 base_length = 0.2
 base_radius = 0.007
-base_area = np.pi * base_radius ** 2
+base_area = np.pi * base_radius**2
 density = 1750
 E = 3e7
 poisson_ratio = 0.5

examples/ButterflyCase/butterfly.py

@@ -38,7 +38,7 @@ class ButterflySimulator(ea.BaseSystemCollection, ea.CallBacks):
 
 total_length = 3.0
 base_radius = 0.25
-base_area = np.pi * base_radius ** 2
+base_area = np.pi * base_radius**2
 density = 5000
 youngs_modulus = 1e4
 poisson_ratio = 0.5
@@ -75,6 +75,7 @@ class ButterflySimulator(ea.BaseSystemCollection, ea.CallBacks):
 
 butterfly_sim.append(butterfly_rod)
 
+
 # Add call backs
 class VelocityCallBack(ea.CallBackBaseClass):
     """

examples/CatenaryCase/catenary.py

@@ -29,7 +29,7 @@ class CatenarySimulator(BaseSystemCollection, Constraints, Forcing, Damping, Cal
 # catenary parameters
 base_length = 1.0
 base_radius = 0.01
-base_area = np.pi * (base_radius ** 2)
+base_area = np.pi * (base_radius**2)
 volume = base_area * base_length
 mass = 0.2
 density = mass / volume

examples/DynamicCantileverCase/analytical_dynamic_cantilever.py

@@ -74,10 +74,8 @@ def __init__(
 
         betas = np.array([0.596864, 1.49418, 2.50025, 3.49999]) * np.pi / base_length
         self.beta = betas[mode]
-        self.omega = (self.beta ** 2) * np.sqrt(
-            youngs_modulus
-            * moment_of_inertia
-            / (density * base_area * base_length ** 4)
+        self.omega = (self.beta**2) * np.sqrt(
+            youngs_modulus * moment_of_inertia / (density * base_area * base_length**4)
         )
 
         nonparametrized_mode_at_end = self._compute_nonparametrized_mode(

examples/DynamicCantileverCase/dynamic_cantilever.py

@@ -50,10 +50,10 @@ class DynamicCantileverSimulator(
     normal = np.array([0.0, 1.0, 0.0])
     base_length = 1
     base_radius = 0.02
-    base_area = np.pi * base_radius ** 2
+    base_area = np.pi * base_radius**2
     youngs_modulus = 1e5
 
-    moment_of_inertia = np.pi / 4 * base_radius ** 4
+    moment_of_inertia = np.pi / 4 * base_radius**4
 
     dl = base_length / n_elem
     dt = dl * 0.05