.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "_gallery/KnotCase/run_knot_simulation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr__gallery_KnotCase_run_knot_simulation.py: Knot Simulation =============== This script simulates the formation of an overhand knot in a soft rod. It demonstrates how to create a controller to manipulate a node on the rod, which can be used for tasks like trajectory tracing or proportional control. .. GENERATED FROM PYTHON SOURCE LINES 9-26 .. code-block:: Python from typing import Any, TypeAlias from numpy.typing import NDArray from elastica.typing import RodType import numpy as np import matplotlib.pyplot as plt from collections import defaultdict import elastica as ea from knot_forcing import TargetPoseProportionalControl from knot_visualization import plot_video3D Position: TypeAlias = NDArray[np.float64] # vector (3) Orientation: TypeAlias = NDArray[np.float64] # SO3 matrix (3, 3) Pose: TypeAlias = tuple[Position, Orientation] .. GENERATED FROM PYTHON SOURCE LINES 27-30 Simulation Setup ---------------- We define a simulator class that inherits from the necessary mixins. .. GENERATED FROM PYTHON SOURCE LINES 30-48 .. code-block:: Python class SoftRodSimulator( ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping, ea.CallBacks, ea.Contact, ): pass simulator = SoftRodSimulator() final_time = 5 dt = 0.0002 .. GENERATED FROM PYTHON SOURCE LINES 49-53 Callback Setup -------------- We also define a callback class to record the position of the rod during the simulation. .. GENERATED FROM PYTHON SOURCE LINES 53-80 .. code-block:: Python class Callback(ea.CallBackBaseClass): """ Records the position of the rod """ def __init__(self, callback_params: dict) -> None: ea.CallBackBaseClass.__init__(self) self.every = 200 self.callback_params = callback_params def make_callback(self, system: RodType, time: float, current_step: int) -> None: if current_step % self.every == 0: self.callback_params["time"].append(time) self.callback_params["step"].append(current_step) self.callback_params["radius"].append(system.radius.copy()) self.callback_params["position"].append(system.position_collection.copy()) self.callback_params["orientation"].append( system.director_collection.copy() ) return recorded_history: dict[str, list[Any]] = defaultdict(list) .. GENERATED FROM PYTHON SOURCE LINES 81-84 Rod Setup --------- Next, we set up the parameters for the rod. .. GENERATED FROM PYTHON SOURCE LINES 84-115 .. code-block:: Python # setting up test params n_elem = 50 start = np.zeros((3,)) direction = np.array([1.0, 0.0, 0.0]) normal = np.array([0.0, 1.0, 0.0]) base_length = 1.2 base_radius = 0.025 density = 2000 youngs_modulus = 1e6 poisson_ratio = 0.5 shear_modulus = youngs_modulus / (2 * (poisson_ratio + 1.0)) # We create the `CosseratRod` object and add it to the simulator. stretchable_rod = ea.CosseratRod.straight_rod( n_elem, start, direction, normal, base_length, base_radius, density, youngs_modulus=youngs_modulus, shear_modulus=shear_modulus, ) simulator.append(stretchable_rod) simulator.collect_diagnostics(stretchable_rod).using( Callback, callback_params=recorded_history ) .. GENERATED FROM PYTHON SOURCE LINES 116-121 Controller Setup ---------------- We define a function that returns the target pose (position and orientation) for the controller at a given time. This function creates the trajectory for the end of the rod to follow to tie the knot. .. GENERATED FROM PYTHON SOURCE LINES 121-161 .. code-block:: Python activation_time = 4 def base_target(t: float, rod: RodType) -> Pose: target_position = direction * base_length - 5 * base_radius * normal if t <= activation_time / 2: ratio = min(2 * t / activation_time, 1.0) angular_ratio = ratio * np.pi * 2 position = target_position * ratio orientation_twist = np.array( [ [0, np.cos(angular_ratio), np.sin(angular_ratio)], [0, -np.sin(angular_ratio), np.cos(angular_ratio)], [1, 0, 0], ], dtype=float, ) else: ratio = min(2 * (t - activation_time / 2) / activation_time, 1.0) R = 8 position = np.array( [ target_position[0] * (1 - ratio), -R * base_radius * np.cos(2 * ratio * 12) * (1 - ratio), -R * base_radius * np.sin(2 * ratio * 12) * (1 - ratio), ] ) angular_ratio = (1 - ratio) * np.pi * 2 orientation_twist = np.array( [ [0, np.cos(angular_ratio), -np.sin(angular_ratio)], [0, np.sin(angular_ratio), np.cos(angular_ratio)], [1, 0, 0], ], dtype=float, ) return position, orientation_twist .. GENERATED FROM PYTHON SOURCE LINES 162-165 We add a `TargetPoseProportionalControl` forcing to the rod. This controller applies forces and torques to drive a specific node of the rod to the target pose. The class is defined in `knot_forcing.py`. .. GENERATED FROM PYTHON SOURCE LINES 165-179 .. code-block:: Python # Control point p = 3e3 pt = 5e0 simulator.add_forcing_to(stretchable_rod).using( TargetPoseProportionalControl, elem_index=0, p_linear_value=p, p_angular_value=pt, target=base_target, ramp_up_time=1e-6, target_history=recorded_history["base_pose"], ) .. GENERATED FROM PYTHON SOURCE LINES 180-183 Boundary Conditions ------------------- We apply boundary conditions to fix the other end of the rod. .. GENERATED FROM PYTHON SOURCE LINES 183-189 .. code-block:: Python # Boundary conditions simulator.constrain(stretchable_rod).using( ea.FixedConstraint, constrained_position_idx=(-1, -20) ) .. GENERATED FROM PYTHON SOURCE LINES 190-194 Contact Setup ------------- We enable self-contact detection for the rod to prevent it from passing through itself. .. GENERATED FROM PYTHON SOURCE LINES 194-200 .. code-block:: Python # Self contact simulator.detect_contact_between(stretchable_rod, stretchable_rod).using( ea.RodSelfContact, k=1e4, nu=3 ) .. GENERATED FROM PYTHON SOURCE LINES 201-204 Environmental Forcing and Damping --------------------------------- We add gravity and damping to the system. .. GENERATED FROM PYTHON SOURCE LINES 204-221 .. code-block:: Python # Gravity simulator.add_forcing_to(stretchable_rod).using( ea.GravityForces, acc_gravity=np.array([0.0, 0.0, -9.80665]) ) # Damping damping_constant = 5.0 simulator.dampen(stretchable_rod).using( ea.AnalyticalLinearDamper, translational_damping_constant=damping_constant, rotational_damping_constant=damping_constant * 0.01, time_step=dt, ) simulator.dampen(stretchable_rod).using(ea.LaplaceDissipationFilter, filter_order=5) .. GENERATED FROM PYTHON SOURCE LINES 222-225 Finalize and Run ---------------- We finalize the simulator and create the time-stepper. .. GENERATED FROM PYTHON SOURCE LINES 225-237 .. code-block:: Python # Finalize and run the simulation simulator.finalize() timestepper = ea.PositionVerlet() total_steps = int(final_time / dt) print("Total steps", total_steps) dt = final_time / total_steps time = 0.0 for i in range(total_steps): time = timestepper.step(simulator, time, dt) .. rst-class:: sphx-glr-script-out .. code-block:: none Contact features should be instantiated lastly. Total steps 25000 /home/docs/checkouts/readthedocs.org/user_builds/pyelastica-mirror/checkouts/mod-examples-timestepper/examples/KnotCase/knot_forcing.py:107: NumbaPerformanceWarning: '@' is faster on contiguous arrays, called on (Array(float64, 2, 'A', False, aligned=True), Array(float64, 2, 'C', False, aligned=True)) rotation = orientation.T @ target_orientation /home/docs/checkouts/readthedocs.org/user_builds/pyelastica-mirror/checkouts/mod-examples-timestepper/examples/KnotCase/knot_forcing.py:118: NumbaPerformanceWarning: '@' is faster on contiguous arrays, called on (Array(float64, 2, 'A', False, aligned=True), Array(float64, 1, 'C', False, aligned=True)) external_torques[..., index] -= orientation @ torque .. GENERATED FROM PYTHON SOURCE LINES 238-242 Post-Processing --------------- After the simulation, we can generate a 3D video of the knot tying process. .. GENERATED FROM PYTHON SOURCE LINES 242-246 .. code-block:: Python filename_video = "knot3D.mp4" plot_video3D(recorded_history, video_name=filename_video, margin=0.2, fps=10) .. rst-class:: sphx-glr-script-out .. code-block:: none Video saved as knot3D.mp4 .. GENERATED FROM PYTHON SOURCE LINES 247-252 .. video:: ../../../examples/KnotCase/knot3D.mp4 :width: 720 :autoplay: :muted: :loop: .. GENERATED FROM PYTHON SOURCE LINES 255-257 We can also plot the topological quantities of the knot, such as twist, writhe, and link, as a function of time. .. GENERATED FROM PYTHON SOURCE LINES 257-281 .. code-block:: Python # Plot knot topological quantities timestep = np.asarray(recorded_history["time"]) positions = np.asarray(recorded_history["position"]) orientations = np.asarray(recorded_history["orientation"]) radii = np.asarray(recorded_history["radius"]) total_twist, _ = ea.compute_twist(positions, orientations[:, 0, ...]) total_writhe = ea.compute_writhe(positions, np.float64(base_length), "next_tangent") total_link = ea.compute_link( positions, orientations[:, 0, ...], radii, np.float64(base_length), "next_tangent", ) plt.figure() plt.plot(timestep, total_twist, label="twist") plt.plot(timestep, total_writhe, label="writhe") plt.plot(timestep, total_link, label="link") plt.legend() plt.xlabel("time") plt.ylabel("link-writhe-twist quantity") plt.show() .. image-sg:: /_gallery/KnotCase/images/sphx_glr_run_knot_simulation_001.png :alt: run knot simulation :srcset: /_gallery/KnotCase/images/sphx_glr_run_knot_simulation_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/pyelastica-mirror/checkouts/mod-examples-timestepper/elastica/rod/knot_theory.py:244: NumbaPerformanceWarning: np.dot() is faster on contiguous arrays, called on (Array(float64, 1, 'A', False, aligned=True), Array(float64, 1, 'A', False, aligned=True)) dot_product = np.dot( /home/docs/checkouts/readthedocs.org/user_builds/pyelastica-mirror/checkouts/mod-examples-timestepper/elastica/rod/knot_theory.py:254: NumbaPerformanceWarning: np.dot() is faster on contiguous arrays, called on (Array(float64, 1, 'A', False, aligned=True), Array(float64, 1, 'C', False, aligned=True)) angle *= np.sign(np.dot(tangent[:, i], cross)) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 26.716 seconds) .. _sphx_glr_download__gallery_KnotCase_run_knot_simulation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: run_knot_simulation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: run_knot_simulation.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: run_knot_simulation.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_