feat : track velocity in magnetic example

examples/MagneticRodCases/ConstantMagneticField/magnetic_beam_case.py

@@ -13,7 +13,7 @@
 PLOT_FIGURE = True
 
 
-class MagneticBeamSimulator(BaseSystemCollection, Constraints, Forcing):
+class MagneticBeamSimulator(BaseSystemCollection, Constraints, Forcing, CallBacks):
     pass
 
 
@@ -28,7 +28,7 @@ class MagneticBeamSimulator(BaseSystemCollection, Constraints, Forcing):
 base_radius = 0.15
 base_area = np.pi * base_radius**2
 density = 5000
-nu = 10
+nu = 10.0
 E = 1e6
 I = np.pi / 4 * base_radius**4
 poisson_ratio = 0.5
@@ -66,7 +66,7 @@ class MagneticBeamSimulator(BaseSystemCollection, Constraints, Forcing):
     [np.cos(magnetic_field_angle), np.sin(magnetic_field_angle), 0]
 )
 magnetic_field_object = ConstantMagneticField(
-    magnetic_field_amplitude, ramp_interval=100.0, start_time=0.0, end_time=100000
+    magnetic_field_amplitude, ramp_interval=500.0, start_time=0.0, end_time=100000
 )
 
 # Apply magnetic forces
@@ -75,9 +75,31 @@ class MagneticBeamSimulator(BaseSystemCollection, Constraints, Forcing):
     external_magnetic_field=magnetic_field_object,
 )
 
+# Add callbacks
+class MagneticBeamCallBack(CallBackBaseClass):
+    def __init__(self, step_skip: int, callback_params: dict):
+        CallBackBaseClass.__init__(self)
+        self.every = step_skip
+        self.callback_params = callback_params
+
+    def make_callback(self, system, time, current_step: int):
+        if current_step % self.every == 0:
+            self.callback_params["time"].append(time)
+            self.callback_params["step"].append(current_step)
+            self.callback_params["velocity_norm"].append(
+                np.linalg.norm(system.velocity_collection)
+            )
+
+
+# Add call back for plotting time history of the rod
+post_processing_dict = defaultdict(list)
+magnetic_beam_sim.collect_diagnostics(magnetic_rod).using(
+    MagneticBeamCallBack, step_skip=500, callback_params=post_processing_dict
+)
+
 magnetic_beam_sim.finalize()
 timestepper = PositionVerlet()
-final_time = 250
+final_time = 1000.0
 dl = base_length / n_elem
 dt = 0.01 * dl
 total_steps = int(final_time / dt)
@@ -95,6 +117,20 @@ class MagneticBeamSimulator(BaseSystemCollection, Constraints, Forcing):
         )
         ax.set_xlabel("X")
         ax.set_ylabel("Y")
-        plt.show()
+        fig.show()
+
+        fig = plt.figure(figsize=(10, 8), frameon=True, dpi=150)
+        ax = fig.add_subplot(111)
+        ax.semilogy(
+            np.asarray(post_processing_dict["time"]),
+            np.asarray(post_processing_dict["velocity_norm"]),
+            lw=2,
+            c=to_rgb("xkcd:bluish"),
+        )
+        ax.set_xlabel("t")
+        ax.set_ylabel("|v|")
+        fig.show()
+
+        plt.show() # block
     if SAVE_FIGURE:
         fig.savefig("Magnetic_beam_profile: N=" + str(magnetic_rod.n_elems) + ".png")