fix: update example cases to use consistent syntax for custom forcing/couples

Minor changes:
Removed Seaborn version specification
Checked all example cases are running and producing correct result
Bug in dissipation module

Author: skim0119

elastica/experimental/connection_contact_joint/generic_system_type_connection.py

@@ -107,6 +107,7 @@ def apply_forces(
         index_one: int,
         system_two: SystemType,
         index_two: int,
+        time: np.float64,
     ):
         """
         Apply joint force to the connected systems.
@@ -196,6 +197,7 @@ def apply_torques(
         index_one: int,
         system_two: SystemType,
         index_two: int,
+        time: np.float64,
     ):
         """
         Apply restoring joint torques to the connected systems.
@@ -214,7 +216,6 @@ def apply_torques(
             Index of second system for joint.
 
         """
-        pass
 
 
 class GenericSystemTypeFixedJoint(GenericSystemTypeFreeJoint):

elastica/experimental/connection_contact_joint/parallel_connection.py

@@ -82,7 +82,7 @@ def __init__(
         offset_btw_rods,
         **kwargs,
     ):
-        super().__init__(np.array(k), np.array(nu))
+        super().__init__(k=np.array(k), nu=np.array(nu))
 
         self.k_repulsive = np.array(k_repulsive)
 
@@ -96,7 +96,7 @@ def __init__(
         )
 
     # Apply force is same as free joint
-    def apply_forces(self, system_one, index_one, system_two, index_two):
+    def apply_forces(self, system_one, index_one, system_two, index_two, time):
         # TODO: documentation
 
         (
@@ -191,7 +191,7 @@ def _apply_forces(
 
         # Compute spring force between two rods
         distance_vector = surface_position_rod_two - surface_position_rod_one
-        np.round_(distance_vector, 12, distance_vector)
+        np.round(distance_vector, 12, distance_vector)
         spring_force = k * (distance_vector)
 
         # Damping force
@@ -220,7 +220,7 @@ def _apply_forces(
             + offset_rod_two
         )
 
-        np.round_(penetration_strain, 12, penetration_strain)
+        np.round(penetration_strain, 12, penetration_strain)
         # Contact present only if rods penetrate to each other
         idx_penetrate = np.where(penetration_strain < 0)[0]
         k_contact = np.zeros(index_one.shape[0])
@@ -256,7 +256,7 @@ def _apply_forces(
             spring_force,
         )
 
-    def apply_torques(self, system_one, index_one, system_two, index_two):
+    def apply_torques(self, system_one, index_one, system_two, index_two, time):
         # pass
 
         self._apply_torques(

elastica/modules/base_system.py

@@ -18,6 +18,7 @@
 )
 
 import numpy as np
+from itertools import chain
 
 from collections.abc import MutableSequence
 
@@ -267,7 +268,9 @@ def constrain_rates(self, time: np.float64) -> None:
         Call constrain rates functions for all features.
         Features are registered in _feature_group_constrain_rates.
         """
-        for func in self._feature_group_constrain_rates:
+        for func in chain(
+            self._feature_group_constrain_rates, self._feature_group_damping
+        ):
             func(time=time)
 
     @final

examples/Binder/2_Slithering_Snake.ipynb

@@ -630,12 +630,11 @@
     "    ax.set_xlim([0 - margin, 3 + margin])\n",
     "    ax.set_ylim([-1.5 - margin, 1.5 + margin])\n",
     "    with writer.saving(fig, video_name, dpi=100):\n",
-    "        with plt.style.context(\"seaborn-v0_8-whitegrid\"):\n",
-    "            for time in range(1, len(t), step):\n",
-    "                rod_lines_2d.set_xdata([positions_over_time[time][2]])\n",
-    "                rod_lines_2d.set_ydata([positions_over_time[time][0]])\n",
+    "        for time in range(1, len(t), step):\n",
+    "            rod_lines_2d.set_xdata([positions_over_time[time][2]])\n",
+    "            rod_lines_2d.set_ydata([positions_over_time[time][0]])\n",
     "\n",
-    "                writer.grab_frame()\n",
+    "            writer.grab_frame()\n",
     "    plt.close(fig)\n",
     "\n",
     "\n",
@@ -696,13 +695,12 @@
     "        linewidth=3,\n",
     "    )[0]\n",
     "    with writer.saving(fig, video_name, dpi=100):\n",
-    "        with plt.style.context(\"seaborn-v0_8-whitegrid\"):\n",
-    "            for time in range(1, len(t), step):\n",
-    "                rod_lines_3d.set_xdata([positions_over_time[time][2]])\n",
-    "                rod_lines_3d.set_ydata([positions_over_time[time][0]])\n",
-    "                rod_lines_3d.set_3d_properties([positions_over_time[time][1]])\n",
+    "        for time in range(1, len(t), step):\n",
+    "            rod_lines_3d.set_xdata([positions_over_time[time][2]])\n",
+    "            rod_lines_3d.set_ydata([positions_over_time[time][0]])\n",
+    "            rod_lines_3d.set_3d_properties(positions_over_time[time][1])\n",
     "\n",
-    "                writer.grab_frame()\n",
+    "            writer.grab_frame()\n",
     "    plt.close(fig)\n",
     "\n",
     "\n",

examples/CantileverDistributedLoad/cantilever_distrubuted_load_postprecessing.py

@@ -241,41 +241,38 @@ def plot_video_with_surface(
         video_name_3D = folder_name + "3D_" + video_name
 
         with writer.saving(fig, video_name_3D, dpi):
-            with plt.style.context("seaborn-v0_8-whitegrid"):
-                for time_idx in tqdm(range(0, sim_time.shape[0], int(step))):
+            for time_idx in tqdm(range(0, sim_time.shape[0], int(step))):
 
-                    for rod_idx in range(n_visualized_rods):
-                        inst_position, inst_radius = rod_history_unpacker(
-                            rod_idx, time_idx
-                        )
-                        if not inst_position.shape[1] == inst_radius.shape[0]:
-                            inst_position = 0.5 * (
-                                inst_position[..., 1:] + inst_position[..., :-1]
-                            )
-
-                        rod_scatters[rod_idx]._offsets3d = (
-                            inst_position[0],
-                            inst_position[1],
-                            inst_position[2],
+                for rod_idx in range(n_visualized_rods):
+                    inst_position, inst_radius = rod_history_unpacker(rod_idx, time_idx)
+                    if not inst_position.shape[1] == inst_radius.shape[0]:
+                        inst_position = 0.5 * (
+                            inst_position[..., 1:] + inst_position[..., :-1]
                         )
 
-                        # rod_scatters[rod_idx].set_offsets(inst_position[:2].T)
-                        rod_scatters[rod_idx].set_sizes(
-                            np.pi * (scaling_factor * inst_radius) ** 2
+                    rod_scatters[rod_idx]._offsets3d = (
+                        inst_position[0],
+                        inst_position[1],
+                        inst_position[2],
+                    )
+
+                    # rod_scatters[rod_idx].set_offsets(inst_position[:2].T)
+                    rod_scatters[rod_idx].set_sizes(
+                        np.pi * (scaling_factor * inst_radius) ** 2
+                    )
+
+                if sphere_flag:
+                    for sphere_idx in range(n_visualized_spheres):
+                        sphere_position, _ = sphere_history_unpacker(
+                            sphere_idx, time_idx
+                        )
+                        sphere_artists[sphere_idx]._offsets3d = (
+                            sphere_position[0],
+                            sphere_position[1],
+                            sphere_position[2],
                         )
 
-                    if sphere_flag:
-                        for sphere_idx in range(n_visualized_spheres):
-                            sphere_position, _ = sphere_history_unpacker(
-                                sphere_idx, time_idx
-                            )
-                            sphere_artists[sphere_idx]._offsets3d = (
-                                sphere_position[0],
-                                sphere_position[1],
-                                sphere_position[2],
-                            )
-
-                    writer.grab_frame()
+                writer.grab_frame()
 
         # Be a good boy and close figures
         # https://stackoverflow.com/a/37451036

examples/CantileverDistributedLoad/cantilever_nonconservative_distributed_load.py

@@ -11,6 +11,12 @@
 )
 
 
+class SquareRodSimulator(
+    ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping, ea.CallBacks
+):
+    pass
+
+
 def cantilever_subjected_to_a_nonconservative_load(
     n_elem,
     base_length,
@@ -21,11 +27,6 @@ def cantilever_subjected_to_a_nonconservative_load(
     animation=False,
     plot_figure_equilibrium=False,
 ):
-    class SquareRodSimulator(
-        ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping, ea.CallBacks
-    ):
-        pass
-
     square_rod_sim = SquareRodSimulator()
 
     square_rod = ea.CosseratRod.straight_rod(
@@ -118,7 +119,6 @@ def make_callback(self, system, time, current_step: int):
     timestepper = ea.PositionVerlet()
 
     total_steps = int(final_time / dt)
-    print(square_rod_sim)
     print("Total steps", total_steps)
     ea.integrate(timestepper, square_rod_sim, final_time, total_steps)
 
@@ -229,15 +229,15 @@ def make_callback(self, system, time, current_step: int):
         load_on_rod,
         x_tip_experiment,
         color="blue",
-        marker="s",
+        marker="o",
         linestyle="None",
         label="x_tip/L",
     )
     plt.scatter(
         load_on_rod,
         y_tip_experiment,
         color="red",
-        marker="s",
+        marker="o",
         linestyle="None",
         label="y_tip/L",
     )

examples/CantileverTransversalLoadCase/cantilever_transversal_load_postprocessing.py

@@ -121,40 +121,37 @@ def plot_video_with_surface(
         video_name_3D = folder_name + "3D_" + video_name
 
         with writer.saving(fig, video_name_3D, dpi):
-            with plt.style.context("seaborn-v0_8-whitegrid"):
-                for time_idx in tqdm(range(0, sim_time.shape[0], int(step))):
+            for time_idx in tqdm(range(0, sim_time.shape[0], int(step))):
 
-                    for rod_idx in range(n_visualized_rods):
-                        inst_position, inst_radius = rod_history_unpacker(
-                            rod_idx, time_idx
-                        )
-                        if not inst_position.shape[1] == inst_radius.shape[0]:
-                            inst_position = 0.5 * (
-                                inst_position[..., 1:] + inst_position[..., :-1]
-                            )
-
-                        rod_scatters[rod_idx]._offsets3d = (
-                            inst_position[0],
-                            inst_position[1],
-                            inst_position[2],
+                for rod_idx in range(n_visualized_rods):
+                    inst_position, inst_radius = rod_history_unpacker(rod_idx, time_idx)
+                    if not inst_position.shape[1] == inst_radius.shape[0]:
+                        inst_position = 0.5 * (
+                            inst_position[..., 1:] + inst_position[..., :-1]
                         )
 
-                        rod_scatters[rod_idx].set_sizes(
-                            np.pi * (scaling_factor * inst_radius) ** 2
+                    rod_scatters[rod_idx]._offsets3d = (
+                        inst_position[0],
+                        inst_position[1],
+                        inst_position[2],
+                    )
+
+                    rod_scatters[rod_idx].set_sizes(
+                        np.pi * (scaling_factor * inst_radius) ** 2
+                    )
+
+                if sphere_flag:
+                    for sphere_idx in range(n_visualized_spheres):
+                        sphere_position, _ = sphere_history_unpacker(
+                            sphere_idx, time_idx
+                        )
+                        sphere_artists[sphere_idx]._offsets3d = (
+                            sphere_position[0],
+                            sphere_position[1],
+                            sphere_position[2],
                         )
 
-                    if sphere_flag:
-                        for sphere_idx in range(n_visualized_spheres):
-                            sphere_position, _ = sphere_history_unpacker(
-                                sphere_idx, time_idx
-                            )
-                            sphere_artists[sphere_idx]._offsets3d = (
-                                sphere_position[0],
-                                sphere_position[1],
-                                sphere_position[2],
-                            )
-
-                    writer.grab_frame()
+                writer.grab_frame()
 
         # Be a good boy and close figures
         # https://stackoverflow.com/a/37451036

examples/CatenaryCase/post_processing.py

@@ -1,8 +1,8 @@
 import numpy as np
 import matplotlib
 
-matplotlib.use("Agg")  # Must be before importing matplotlib.pyplot or pylab!
 from matplotlib import pyplot as plt
+import matplotlib.animation as manimation
 from tqdm import tqdm
 import scipy as sci
 
@@ -14,8 +14,6 @@ def plot_video(
     xlim: tuple[float, float] = (0, 4),
     ylim: tuple[float, float] = (-1, 1),
 ) -> None:
-    import matplotlib.animation as manimation
-
     positions_over_time = np.array(plot_params["position"])
 
     print("plot video")
@@ -52,7 +50,6 @@ def plot_catenary(
     """
     Catenary analytical solution from Routh, Edward John (1891). "Chapter X: On Strings". A Treatise on Analytical Statics. University Press.
     """
-    matplotlib.use("TkAgg")
     position = np.array(plot_params["position"])
     lowest_point = np.min(position[-1][2])
     x_catenary = np.linspace(0, base_length, 100)

examples/ContinuumSnakeWithLiftingWaveCase/snake_contact.py

@@ -293,7 +293,9 @@ def _allowed_system_two(self) -> list[Type]:
         # Modify this list to include the allowed system types for contact
         return [SurfaceBase]
 
-    def apply_contact(self, system_one: RodType, system_two: SystemType) -> None:
+    def apply_contact(
+        self, system_one: RodType, system_two: SystemType, time: float
+    ) -> None:
         """
         In the case of contact with the plane, this function computes the plane reaction force on the element.
 

examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py

@@ -120,22 +120,21 @@ def apply_forces(self, system, time: np.float64 = 0.0):
     rod_one, rod_two, (0, n_elem), (0, n_elem)
 )
 
-for i in range(n_elem):
-    parallel_connection_sim.connect(
-        first_rod=rod_one, second_rod=rod_two, first_connect_idx=i, second_connect_idx=i
-    ).using(
-        SurfaceJointSideBySide,
-        k=1e2,
-        nu=1e-5,
-        k_repulsive=1e3,
-        rod_one_direction_vec_in_material_frame=rod_one_direction_vec_in_material_frame[
-            :, i
-        ],
-        rod_two_direction_vec_in_material_frame=rod_two_direction_vec_in_material_frame[
-            :, i
-        ],
-        offset_btw_rods=offset_btw_rods[i],
-    )  # k=kg/s2 nu=kg/s 1e-2
+indices = np.arange(n_elem)
+parallel_connection_sim.connect(
+    first_rod=rod_one,
+    second_rod=rod_two,
+    first_connect_idx=indices,
+    second_connect_idx=indices,
+).using(
+    SurfaceJointSideBySide,
+    k=1e2,
+    nu=1e-5,
+    k_repulsive=1e3,
+    rod_one_direction_vec_in_material_frame=rod_one_direction_vec_in_material_frame,
+    rod_two_direction_vec_in_material_frame=rod_two_direction_vec_in_material_frame,
+    offset_btw_rods=offset_btw_rods,
+)  # k=kg/s2 nu=kg/s 1e-2
 
 
 # add damping