From a511e49af1808de81eef6a49e986373b326785b2 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 8 Jun 2022 04:28:04 -0500 Subject: [PATCH 01/37] resolve #106: avoid 0.0 default for rampup --- elastica/external_forces.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/elastica/external_forces.py b/elastica/external_forces.py index 6ef9645c..b0de600d 100644 --- a/elastica/external_forces.py +++ b/elastica/external_forces.py @@ -139,7 +139,7 @@ class EndpointForces(NoForces): """ - def __init__(self, start_force, end_force, ramp_up_time=0.0): + def __init__(self, start_force, end_force, ramp_up_time): """ Parameters @@ -157,15 +157,10 @@ def __init__(self, start_force, end_force, ramp_up_time=0.0): super(EndpointForces, self).__init__() self.start_force = start_force self.end_force = end_force - assert ramp_up_time >= 0.0 + assert ramp_up_time > 0.0 self.ramp_up_time = ramp_up_time def apply_forces(self, system, time=0.0): - # factor = min(1.0, time / self.ramp_up_time) - # - # system.external_forces[..., 0] += self.start_force * factor - # system.external_forces[..., -1] += self.end_force * factor - self.compute_end_point_forces( system.external_forces, self.start_force, @@ -307,7 +302,7 @@ def __init__( phase_shift, direction, rest_lengths, - ramp_up_time=0.0, + ramp_up_time, with_spline=False, ): """ @@ -340,7 +335,7 @@ def __init__( self.wave_number = wave_number self.phase_shift = phase_shift - assert ramp_up_time >= 0.0 + assert ramp_up_time > 0.0 self.ramp_up_time = ramp_up_time # s is the position of nodes on the rod, we go from node=1 to node=nelem-1, because there is no From 51a5fcf873af4e1ff8e1a791412b692436956650 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 11 Jun 2022 14:02:08 -0500 Subject: [PATCH 02/37] bugfix:numpy import for rod contact cases --- .../RodRodContact/rod_rod_contact_inclined_validation.py | 1 + .../RodRodContact/rod_rod_contact_parallel_validation.py | 1 + .../RodSelfContact/PlectonemesCase/plectoneme_case.py | 1 + .../RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py | 1 + 4 files changed, 4 insertions(+) diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py index 8f71addb..bb4e13e0 100644 --- a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py +++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py @@ -1,3 +1,4 @@ +import numpy as np import sys sys.path.append("../../../") diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py index 09b18ca3..e5f81969 100644 --- a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py +++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py @@ -1,3 +1,4 @@ +import numpy as np import sys sys.path.append("../../../") diff --git a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py index 1f7e817b..7eb8e64d 100644 --- a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py +++ b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py @@ -1,3 +1,4 @@ +import numpy as np import sys sys.path.append("../../../../") diff --git a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py index 0ce2fdce..ab9e6371 100644 --- a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py +++ b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py @@ -1,3 +1,4 @@ +import numpy as np import sys sys.path.append("../../../../") From 8b7cc395a08a8b0fae9ab8d50233b4870a1918c2 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 11 Jun 2022 14:03:01 -0500 Subject: [PATCH 03/37] enhancement:adding example cases for contact between rigid cylinder and rod. --- .../rod_cylinder_contact_validation.py | 290 ++++++++++++++++++ .../rod_cylinder_contact_with_y_normal.py | 266 ++++++++++++++++ 2 files changed, 556 insertions(+) create mode 100644 examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py create mode 100644 examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py new file mode 100644 index 00000000..cd9f5442 --- /dev/null +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py @@ -0,0 +1,290 @@ +import numpy as np + +# FIXME without appending sys.path make it more generic +import sys + +sys.path.append("../../../") +from elastica import * +from matplotlib import pyplot as plt + + +class SingleRodSingleCylinderInteractionSimulator( + BaseSystemCollection, Constraints, Connections, Forcing, CallBacks +): + pass + + +single_rod_sim = SingleRodSingleCylinderInteractionSimulator() +# setting up test params +n_elem = 50 + +inclination = np.deg2rad(30) +direction = np.array([0.0, np.cos(inclination), np.sin(inclination)]) +normal = np.array([0.0, -np.sin(inclination), np.cos(inclination)]) + +# can be y or z too, meant for testing purposes of rod-body contact in different planes +action_plane_key = "x" + +# can be set to True, checks collision at tips of rod +tip_collision = True + +_roll_key = 0 if action_plane_key == "x" else (1 if action_plane_key == "y" else 2) +if action_plane_key == "x": + global_rot_mat = np.eye(3) +elif action_plane_key == "y": + # Rotate +ve 90 about z + global_rot_mat = np.zeros((3, 3)) + global_rot_mat[0, 1] = -1.0 + global_rot_mat[1, 0] = 1.0 + global_rot_mat[2, 2] = 1.0 +else: + # Rotate -ve 90 abuot y + global_rot_mat = np.zeros((3, 3)) + global_rot_mat[1, 1] = 1.0 + global_rot_mat[0, 2] = 1.0 + global_rot_mat[2, 0] = 1.0 + + +direction = global_rot_mat @ direction +normal = global_rot_mat @ normal + +base_length = 0.5 +base_radius = 0.01 +base_area = np.pi * base_radius ** 2 +density = 1750 +nu = 0.001 +E = 3e5 +poisson_ratio = 0.5 +shear_modulus = E / (1 + poisson_ratio) + +cylinder_start = global_rot_mat @ np.array([0.3, 0.0, 0.0]) +cylinder_direction = global_rot_mat @ np.array([0.0, 0.0, 1.0]) +cylinder_normal = global_rot_mat @ np.array([0.0, 1.0, 0.0]) + +cylinder_height = 0.4 +cylinder_radius = 10.0 * base_radius + +# Cylinder surface starts at 0.2 +tip_offset = 0.0 +if tip_collision: + # The random choice decides which tip of the rod intersects with cylinder + tip_choice = np.random.choice([1, -1]) + tip_offset = 0.5 * tip_choice * base_length * np.cos(inclination) + +start_rod_1 = np.array( + [ + 0.15, + -0.5 * base_length * np.cos(inclination) + tip_offset, + 0.5 * cylinder_height - 0.5 * base_length * np.sin(inclination), + ] +) +start_rod_1 = global_rot_mat @ start_rod_1 + +rod1 = CosseratRod.straight_rod( + n_elem, + start_rod_1, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + shear_modulus=shear_modulus, +) +# Give it an initial push +rod1.velocity_collection[_roll_key, ...] = 0.05 +single_rod_sim.append(rod1) + + +cylinder = Cylinder( + cylinder_start, + cylinder_direction, + cylinder_normal, + cylinder_height, + cylinder_radius, + density, +) +single_rod_sim.append(cylinder) + +single_rod_sim.connect(rod1, cylinder).using(ExternalContact, 1e2, 0.1) + + +# Add call backs +class PositionCollector(CallBackBaseClass): + """ + Call back function for continuum snake + """ + + 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) + # Collect only x + self.callback_params["position"].append(system.position_collection.copy()) + self.callback_params["com"].append(system.compute_position_center_of_mass()) + return + + +recorded_rod_history = defaultdict(list) +single_rod_sim.collect_diagnostics(rod1).using( + PositionCollector, step_skip=200, callback_params=recorded_rod_history +) +recorded_cyl_history = defaultdict(list) +single_rod_sim.collect_diagnostics(cylinder).using( + PositionCollector, step_skip=200, callback_params=recorded_cyl_history +) + +single_rod_sim.finalize() +timestepper = PositionVerlet() +final_time = 2.0 +dl = base_length / n_elem +dt = 1e-4 +total_steps = int(final_time / dt) +print("Total steps", total_steps) + +integrate(timestepper, single_rod_sim, final_time, total_steps) + +if True: + + def make_data_for_cylinder_along_z(cstart, cradius, cheight): + center_x, center_y = cstart[0], cstart[1] + z = np.linspace(0, cheight, 5) + theta = np.linspace(0, 2 * np.pi, 20) + theta_grid, z_grid = np.meshgrid(theta, z) + x_grid = cradius * np.cos(theta_grid) + center_x + y_grid = cradius * np.sin(theta_grid) + center_y + z_grid += cstart[2] + return np.roll([x_grid, y_grid, z_grid], _roll_key) + + XC, YC, ZC = make_data_for_cylinder_along_z( + cylinder_start, cylinder_radius, cylinder_height + ) + + def plot_video( + plot_params: dict, + cylinder_history: dict, + video_name="video.mp4", + margin=0.2, + fps=60, + step=1, + *args, + **kwargs + ): # (time step, x/y/z, node) + import matplotlib.animation as manimation + + plt.rcParams.update({"font.size": 22}) + + # Should give a (n_time, 3, n_elem) array + positions = np.array(plot_params["position"]) + # (n_time, 3) array + com = np.array(plot_params["com"]) + + cylinder_com = np.array(cylinder_history["com"]) + cylinder_origin = cylinder_com - 0.5 * cylinder_height * cylinder_direction + + print("plot video") + FFMpegWriter = manimation.writers["ffmpeg"] + metadata = dict( + title="Movie Test", artist="Matplotlib", comment="Movie support!" + ) + writer = FFMpegWriter(fps=fps, metadata=metadata) + dpi = 50 + + # min_limits = np.roll(np.array([0.0, -0.5 * cylinder_height, 0.0]), _roll_key) + if True: + from mpl_toolkits.mplot3d import Axes3D + + fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) + ax = plt.axes(projection="3d") # fig.add_subplot(111) + ax.grid(b=True, which="minor", color="k", linestyle="--") + ax.grid(b=True, which="major", color="k", linestyle="-") + # plt.axis("square") + i = 0 + (rod_line,) = ax.plot( + positions[i, 0], positions[i, 1], positions[i, 2], lw=3.0 + ) + XC, YC, ZC = make_data_for_cylinder_along_z( + cylinder_origin[i, ...], cylinder_radius, cylinder_height + ) + surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + + min_limits = global_rot_mat @ np.array([0.0, -0.5 * cylinder_height, 0.0]) + min_limits = -np.abs(min_limits) + max_limits = min_limits + cylinder_height + + ax.set_xlim([min_limits[0], max_limits[0]]) + ax.set_ylim([min_limits[1], max_limits[1]]) + ax.set_zlim([min_limits[2], max_limits[2]]) + with writer.saving(fig, video_name, dpi): + with plt.style.context("seaborn-white"): + for i in range(0, positions.shape[0], int(step)): + rod_line.set_xdata(positions[i, 0]) + rod_line.set_ydata(positions[i, 1]) + rod_line.set_3d_properties(positions[i, 2]) + XC, YC, ZC = make_data_for_cylinder_along_z( + cylinder_origin[i, ...], cylinder_radius, cylinder_height + ) + surf.remove() + surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) + writer.grab_frame() + if True: + from matplotlib.patches import Circle + + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + i = 0 + positions = np.roll(positions, -_roll_key, axis=1) + com = np.roll(com, -_roll_key, axis=1) + cstart = np.roll(cylinder_origin, -_roll_key, axis=1) + (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], lw=3.0) + (tip_line,) = ax.plot(com[i, 0], com[i, 1], "k--") + + min_limits = np.array([0.0, -0.5 * cylinder_height, 0.0]) + max_limits = min_limits + cylinder_height + + ax.set_xlim([min_limits[0], max_limits[0]]) + ax.set_ylim([min_limits[1], max_limits[1]]) + + circle_artist = Circle( + (cstart[i, 0], cstart[i, 1]), cylinder_radius, color="g" + ) + ax.add_artist(circle_artist) + ax.set_aspect("equal") + video_name = "2D_" + video_name + with writer.saving(fig, video_name, dpi): + with plt.style.context("fivethirtyeight"): + for i in range(0, positions.shape[0], int(step)): + rod_line.set_xdata(positions[i, 0]) + rod_line.set_ydata(positions[i, 1]) + tip_line.set_xdata(com[:i, 0]) + tip_line.set_ydata(com[:i, 1]) + circle_artist.center = cstart[i, 0], cstart[i, 1] + writer.grab_frame() + + plot_video(recorded_rod_history, recorded_cyl_history, "cylinder_rod_collision.mp4") + + positions = np.array(recorded_rod_history["position"]) + sim_time = np.array(recorded_rod_history["time"]) + colliding_element_idx = n_elem // 2 + if tip_collision: + colliding_element_idx = 0 if tip_choice == 1 else -1 + colliding_element_history = positions[:, :, colliding_element_idx] + fig = plt.figure(3, figsize=(8, 5)) + ax = fig.add_subplot(111) + ax.plot(sim_time, colliding_element_history[:, _roll_key]) + ax.hlines( + cylinder_start[_roll_key] - cylinder_radius - base_radius, + sim_time[0], + sim_time[-1], + "k", + linestyle="dashed", + ) + plt.show() diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py new file mode 100644 index 00000000..b06172b1 --- /dev/null +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py @@ -0,0 +1,266 @@ +import numpy as np + +# FIXME without appending sys.path make it more generic +import sys + +sys.path.append("../../../") +from elastica import * +from matplotlib import pyplot as plt + + +class SingleRodSingleCylinderInteractionSimulator( + BaseSystemCollection, Constraints, Connections, Forcing, CallBacks +): + pass + + +single_rod_sim = SingleRodSingleCylinderInteractionSimulator() +# setting up test params +n_elem = 50 + +inclination = np.deg2rad(30) +direction = np.array([0.0, np.sin(inclination), -np.cos(inclination)]) +normal = np.array([0.0, np.cos(inclination), np.sin(inclination)]) + +base_length = 0.5 +base_radius = 0.01 +base_area = np.pi * base_radius ** 2 +density = 1750 +nu = 0.001 +E = 3e5 +poisson_ratio = 0.5 +shear_modulus = E / (1 + poisson_ratio) + +cylinder_start = np.array([0.3, 0.0, 0.0]) +cylinder_direction = np.array([0.0, 1.0, 0.0]) +cylinder_normal = np.array([1.0, 0.0, 0.0]) + +cylinder_height = 0.4 +cylinder_radius = 10.0 * base_radius + +# can be set to True, checks collision at tips of rod +tip_collision = False + +# Cylinder surface starts at 0.2 +tip_offset = 0.0 + +if tip_collision: + # The random choice decides which tip of the rod intersects with cylinder + tip_choice = np.random.choice([1, -1]) + tip_offset = 0.5 * tip_choice * base_length * np.cos(inclination) + +start_rod_1 = np.array( + [ + 0.15, + 0.5 * cylinder_height - 0.5 * base_length * np.sin(inclination), + 0.5 * base_length * np.cos(inclination) + tip_offset, + ] +) +# start_rod_1[2] = cylinder_radius + base_length + +rod1 = CosseratRod.straight_rod( + n_elem, + start_rod_1, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + shear_modulus=shear_modulus, +) +# Give it an initial push +rod1.velocity_collection[0, ...] = 0.05 +single_rod_sim.append(rod1) + + +cylinder = Cylinder( + cylinder_start, + cylinder_direction, + cylinder_normal, + cylinder_height, + cylinder_radius, + density, +) +single_rod_sim.append(cylinder) + +single_rod_sim.connect(rod1, cylinder).using(ExternalContact, 1e2, 0.1) + + +# Add call backs +class PositionCollector(CallBackBaseClass): + """ + Call back function for continuum snake + """ + + 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) + # Collect only x + self.callback_params["position"].append(system.position_collection.copy()) + self.callback_params["com"].append(system.compute_position_center_of_mass()) + return + + +recorded_rod_history = defaultdict(list) +single_rod_sim.collect_diagnostics(rod1).using( + PositionCollector, step_skip=200, callback_params=recorded_rod_history +) +recorded_cyl_history = defaultdict(list) +single_rod_sim.collect_diagnostics(cylinder).using( + PositionCollector, step_skip=200, callback_params=recorded_cyl_history +) + +single_rod_sim.finalize() +timestepper = PositionVerlet() +final_time = 2.0 +dl = base_length / n_elem +dt = 1e-4 +total_steps = int(final_time / dt) +print("Total steps", total_steps) + +integrate(timestepper, single_rod_sim, final_time, total_steps) + +if True: + + def make_data_for_cylinder_along_y(cstart, cradius, cheight): + center_x, center_z = cstart[0], cstart[1] + y = np.linspace(0, cheight, 5) + theta = np.linspace(0, 2 * np.pi, 20) + theta_grid, y_grid = np.meshgrid(theta, y) + x_grid = cradius * np.cos(theta_grid) + center_x + z_grid = cradius * np.sin(theta_grid) + center_z + y_grid += cstart[2] + return [x_grid, y_grid, z_grid] + + XC, YC, ZC = make_data_for_cylinder_along_y( + cylinder_start, cylinder_radius, cylinder_height + ) + + def plot_video( + plot_params: dict, + cylinder_history: dict, + video_name="video.mp4", + margin=0.2, + fps=60, + step=1, + *args, + **kwargs + ): # (time step, x/y/z, node) + import matplotlib.animation as manimation + + plt.rcParams.update({"font.size": 22}) + + # Should give a (n_time, 3, n_elem) array + positions = np.array(plot_params["position"]) + # (n_time, 3) array + com = np.array(plot_params["com"]) + + cylinder_com = np.array(cylinder_history["com"]) + cylinder_origin = cylinder_com - 0.5 * cylinder_height * cylinder_direction + + print("plot video") + FFMpegWriter = manimation.writers["ffmpeg"] + metadata = dict( + title="Movie Test", artist="Matplotlib", comment="Movie support!" + ) + writer = FFMpegWriter(fps=fps, metadata=metadata) + dpi = 50 + + # min_limits = np.roll(np.array([0.0, -0.5 * cylinder_height, 0.0]), _roll_key) + if True: + from mpl_toolkits.mplot3d import Axes3D + + fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) + ax = plt.axes(projection="3d") # fig.add_subplot(111) + ax.grid(b=True, which="minor", color="k", linestyle="--") + ax.grid(b=True, which="major", color="k", linestyle="-") + # plt.axis("square") + i = 0 + (rod_line,) = ax.plot( + positions[i, 0], positions[i, 1], positions[i, 2], lw=3.0 + ) + XC, YC, ZC = make_data_for_cylinder_along_y( + cylinder_origin[i, ...], cylinder_radius, cylinder_height + ) + surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + + min_limits = np.array([0.0, 0.0, -0.5 * cylinder_height]) + min_limits = -np.abs(min_limits) + max_limits = min_limits + cylinder_height + + ax.set_xlim([min_limits[0], max_limits[0]]) + ax.set_ylim([min_limits[1], max_limits[1]]) + ax.set_zlim([min_limits[2], max_limits[2]]) + with writer.saving(fig, video_name, dpi): + with plt.style.context("seaborn-white"): + for i in range(0, positions.shape[0], int(step)): + rod_line.set_xdata(positions[i, 0]) + rod_line.set_ydata(positions[i, 1]) + rod_line.set_3d_properties(positions[i, 2]) + XC, YC, ZC = make_data_for_cylinder_along_y( + cylinder_origin[i, ...], cylinder_radius, cylinder_height + ) + surf.remove() + surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) + writer.grab_frame() + if True: + from matplotlib.patches import Circle + + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + i = 0 + cstart = cylinder_origin + (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], lw=3.0) + (tip_line,) = ax.plot(com[i, 0], com[i, 1], "k--") + + min_limits = np.array([0.0, 0.0, -0.5 * cylinder_height]) + max_limits = min_limits + cylinder_height + + ax.set_xlim([min_limits[0], max_limits[0]]) + ax.set_ylim([min_limits[2], max_limits[2]]) + + circle_artist = Circle( + (cstart[i, 0], cstart[i, 2]), cylinder_radius, color="g" + ) + ax.add_artist(circle_artist) + ax.set_aspect("equal") + video_name = "2D_" + video_name + with writer.saving(fig, video_name, dpi): + with plt.style.context("fivethirtyeight"): + for i in range(0, positions.shape[0], int(step)): + rod_line.set_xdata(positions[i, 0]) + rod_line.set_ydata(positions[i, 2]) + tip_line.set_xdata(com[:i, 0]) + tip_line.set_ydata(com[:i, 2]) + circle_artist.center = cstart[i, 0], cstart[i, 2] + writer.grab_frame() + + plot_video(recorded_rod_history, recorded_cyl_history, "cylinder_rod_collision.mp4") + + positions = np.array(recorded_rod_history["position"]) + sim_time = np.array(recorded_rod_history["time"]) + colliding_element_idx = n_elem // 2 + if tip_collision: + colliding_element_idx = -1 if tip_choice == 1 else 0 + colliding_element_history = positions[:, :, colliding_element_idx] + fig = plt.figure(3, figsize=(8, 5)) + ax = fig.add_subplot(111) + ax.plot(sim_time, colliding_element_history[:, 0]) + ax.hlines( + cylinder_start[0] - cylinder_radius - base_radius, + sim_time[0], + sim_time[-1], + "k", + linestyle="dashed", + ) + plt.show() From fcb45e351485a2c1f07b5ca3dd2f7a5d142dfd01 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 11 Jun 2022 14:32:32 -0500 Subject: [PATCH 04/37] fix: flag names --- .../rod_cylinder_contact_validation.py | 143 +++++++++--------- .../rod_cylinder_contact_with_y_normal.py | 140 +++++++++-------- 2 files changed, 139 insertions(+), 144 deletions(-) diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py index cd9f5442..2fcc2025 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py @@ -14,6 +14,9 @@ class SingleRodSingleCylinderInteractionSimulator( pass +# Options +PLOT_FIGURE = True + single_rod_sim = SingleRodSingleCylinderInteractionSimulator() # setting up test params n_elem = 50 @@ -149,7 +152,7 @@ def make_callback(self, system, time, current_step: int): integrate(timestepper, single_rod_sim, final_time, total_steps) -if True: +if PLOT_FIGURE: def make_data_for_cylinder_along_z(cstart, cradius, cheight): center_x, center_y = cstart[0], cstart[1] @@ -196,78 +199,72 @@ def plot_video( dpi = 50 # min_limits = np.roll(np.array([0.0, -0.5 * cylinder_height, 0.0]), _roll_key) - if True: - from mpl_toolkits.mplot3d import Axes3D - - fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) - ax = plt.axes(projection="3d") # fig.add_subplot(111) - ax.grid(b=True, which="minor", color="k", linestyle="--") - ax.grid(b=True, which="major", color="k", linestyle="-") - # plt.axis("square") - i = 0 - (rod_line,) = ax.plot( - positions[i, 0], positions[i, 1], positions[i, 2], lw=3.0 - ) - XC, YC, ZC = make_data_for_cylinder_along_z( - cylinder_origin[i, ...], cylinder_radius, cylinder_height - ) - surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_zlabel("z") - - min_limits = global_rot_mat @ np.array([0.0, -0.5 * cylinder_height, 0.0]) - min_limits = -np.abs(min_limits) - max_limits = min_limits + cylinder_height - - ax.set_xlim([min_limits[0], max_limits[0]]) - ax.set_ylim([min_limits[1], max_limits[1]]) - ax.set_zlim([min_limits[2], max_limits[2]]) - with writer.saving(fig, video_name, dpi): - with plt.style.context("seaborn-white"): - for i in range(0, positions.shape[0], int(step)): - rod_line.set_xdata(positions[i, 0]) - rod_line.set_ydata(positions[i, 1]) - rod_line.set_3d_properties(positions[i, 2]) - XC, YC, ZC = make_data_for_cylinder_along_z( - cylinder_origin[i, ...], cylinder_radius, cylinder_height - ) - surf.remove() - surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) - writer.grab_frame() - if True: - from matplotlib.patches import Circle - - fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) - ax = fig.add_subplot(111) - i = 0 - positions = np.roll(positions, -_roll_key, axis=1) - com = np.roll(com, -_roll_key, axis=1) - cstart = np.roll(cylinder_origin, -_roll_key, axis=1) - (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], lw=3.0) - (tip_line,) = ax.plot(com[i, 0], com[i, 1], "k--") - - min_limits = np.array([0.0, -0.5 * cylinder_height, 0.0]) - max_limits = min_limits + cylinder_height - - ax.set_xlim([min_limits[0], max_limits[0]]) - ax.set_ylim([min_limits[1], max_limits[1]]) - - circle_artist = Circle( - (cstart[i, 0], cstart[i, 1]), cylinder_radius, color="g" - ) - ax.add_artist(circle_artist) - ax.set_aspect("equal") - video_name = "2D_" + video_name - with writer.saving(fig, video_name, dpi): - with plt.style.context("fivethirtyeight"): - for i in range(0, positions.shape[0], int(step)): - rod_line.set_xdata(positions[i, 0]) - rod_line.set_ydata(positions[i, 1]) - tip_line.set_xdata(com[:i, 0]) - tip_line.set_ydata(com[:i, 1]) - circle_artist.center = cstart[i, 0], cstart[i, 1] - writer.grab_frame() + from mpl_toolkits.mplot3d import Axes3D + + fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) + ax = plt.axes(projection="3d") # fig.add_subplot(111) + ax.grid(b=True, which="minor", color="k", linestyle="--") + ax.grid(b=True, which="major", color="k", linestyle="-") + # plt.axis("square") + i = 0 + (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], positions[i, 2], lw=3.0) + XC, YC, ZC = make_data_for_cylinder_along_z( + cylinder_origin[i, ...], cylinder_radius, cylinder_height + ) + surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + + min_limits = global_rot_mat @ np.array([0.0, -0.5 * cylinder_height, 0.0]) + min_limits = -np.abs(min_limits) + max_limits = min_limits + cylinder_height + + ax.set_xlim([min_limits[0], max_limits[0]]) + ax.set_ylim([min_limits[1], max_limits[1]]) + ax.set_zlim([min_limits[2], max_limits[2]]) + with writer.saving(fig, video_name, dpi): + with plt.style.context("seaborn-white"): + for i in range(0, positions.shape[0], int(step)): + rod_line.set_xdata(positions[i, 0]) + rod_line.set_ydata(positions[i, 1]) + rod_line.set_3d_properties(positions[i, 2]) + XC, YC, ZC = make_data_for_cylinder_along_z( + cylinder_origin[i, ...], cylinder_radius, cylinder_height + ) + surf.remove() + surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) + writer.grab_frame() + from matplotlib.patches import Circle + + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + i = 0 + positions = np.roll(positions, -_roll_key, axis=1) + com = np.roll(com, -_roll_key, axis=1) + cstart = np.roll(cylinder_origin, -_roll_key, axis=1) + (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], lw=3.0) + (tip_line,) = ax.plot(com[i, 0], com[i, 1], "k--") + + min_limits = np.array([0.0, -0.5 * cylinder_height, 0.0]) + max_limits = min_limits + cylinder_height + + ax.set_xlim([min_limits[0], max_limits[0]]) + ax.set_ylim([min_limits[1], max_limits[1]]) + + circle_artist = Circle((cstart[i, 0], cstart[i, 1]), cylinder_radius, color="g") + ax.add_artist(circle_artist) + ax.set_aspect("equal") + video_name = "2D_" + video_name + with writer.saving(fig, video_name, dpi): + with plt.style.context("fivethirtyeight"): + for i in range(0, positions.shape[0], int(step)): + rod_line.set_xdata(positions[i, 0]) + rod_line.set_ydata(positions[i, 1]) + tip_line.set_xdata(com[:i, 0]) + tip_line.set_ydata(com[:i, 1]) + circle_artist.center = cstart[i, 0], cstart[i, 1] + writer.grab_frame() plot_video(recorded_rod_history, recorded_cyl_history, "cylinder_rod_collision.mp4") diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py index b06172b1..359fab9b 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py @@ -14,6 +14,9 @@ class SingleRodSingleCylinderInteractionSimulator( pass +# Options +PLOT_FIGURE = True + single_rod_sim = SingleRodSingleCylinderInteractionSimulator() # setting up test params n_elem = 50 @@ -127,7 +130,7 @@ def make_callback(self, system, time, current_step: int): integrate(timestepper, single_rod_sim, final_time, total_steps) -if True: +if PLOT_FIGURE: def make_data_for_cylinder_along_y(cstart, cradius, cheight): center_x, center_z = cstart[0], cstart[1] @@ -174,76 +177,71 @@ def plot_video( dpi = 50 # min_limits = np.roll(np.array([0.0, -0.5 * cylinder_height, 0.0]), _roll_key) - if True: - from mpl_toolkits.mplot3d import Axes3D - - fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) - ax = plt.axes(projection="3d") # fig.add_subplot(111) - ax.grid(b=True, which="minor", color="k", linestyle="--") - ax.grid(b=True, which="major", color="k", linestyle="-") - # plt.axis("square") - i = 0 - (rod_line,) = ax.plot( - positions[i, 0], positions[i, 1], positions[i, 2], lw=3.0 - ) - XC, YC, ZC = make_data_for_cylinder_along_y( - cylinder_origin[i, ...], cylinder_radius, cylinder_height - ) - surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_zlabel("z") - - min_limits = np.array([0.0, 0.0, -0.5 * cylinder_height]) - min_limits = -np.abs(min_limits) - max_limits = min_limits + cylinder_height - - ax.set_xlim([min_limits[0], max_limits[0]]) - ax.set_ylim([min_limits[1], max_limits[1]]) - ax.set_zlim([min_limits[2], max_limits[2]]) - with writer.saving(fig, video_name, dpi): - with plt.style.context("seaborn-white"): - for i in range(0, positions.shape[0], int(step)): - rod_line.set_xdata(positions[i, 0]) - rod_line.set_ydata(positions[i, 1]) - rod_line.set_3d_properties(positions[i, 2]) - XC, YC, ZC = make_data_for_cylinder_along_y( - cylinder_origin[i, ...], cylinder_radius, cylinder_height - ) - surf.remove() - surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) - writer.grab_frame() - if True: - from matplotlib.patches import Circle - - fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) - ax = fig.add_subplot(111) - i = 0 - cstart = cylinder_origin - (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], lw=3.0) - (tip_line,) = ax.plot(com[i, 0], com[i, 1], "k--") - - min_limits = np.array([0.0, 0.0, -0.5 * cylinder_height]) - max_limits = min_limits + cylinder_height - - ax.set_xlim([min_limits[0], max_limits[0]]) - ax.set_ylim([min_limits[2], max_limits[2]]) - - circle_artist = Circle( - (cstart[i, 0], cstart[i, 2]), cylinder_radius, color="g" - ) - ax.add_artist(circle_artist) - ax.set_aspect("equal") - video_name = "2D_" + video_name - with writer.saving(fig, video_name, dpi): - with plt.style.context("fivethirtyeight"): - for i in range(0, positions.shape[0], int(step)): - rod_line.set_xdata(positions[i, 0]) - rod_line.set_ydata(positions[i, 2]) - tip_line.set_xdata(com[:i, 0]) - tip_line.set_ydata(com[:i, 2]) - circle_artist.center = cstart[i, 0], cstart[i, 2] - writer.grab_frame() + from mpl_toolkits.mplot3d import Axes3D + + fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) + ax = plt.axes(projection="3d") # fig.add_subplot(111) + ax.grid(b=True, which="minor", color="k", linestyle="--") + ax.grid(b=True, which="major", color="k", linestyle="-") + # plt.axis("square") + i = 0 + (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], positions[i, 2], lw=3.0) + XC, YC, ZC = make_data_for_cylinder_along_y( + cylinder_origin[i, ...], cylinder_radius, cylinder_height + ) + surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + + min_limits = np.array([0.0, 0.0, -0.5 * cylinder_height]) + min_limits = -np.abs(min_limits) + max_limits = min_limits + cylinder_height + + ax.set_xlim([min_limits[0], max_limits[0]]) + ax.set_ylim([min_limits[1], max_limits[1]]) + ax.set_zlim([min_limits[2], max_limits[2]]) + with writer.saving(fig, video_name, dpi): + with plt.style.context("seaborn-white"): + for i in range(0, positions.shape[0], int(step)): + rod_line.set_xdata(positions[i, 0]) + rod_line.set_ydata(positions[i, 1]) + rod_line.set_3d_properties(positions[i, 2]) + XC, YC, ZC = make_data_for_cylinder_along_y( + cylinder_origin[i, ...], cylinder_radius, cylinder_height + ) + surf.remove() + surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) + writer.grab_frame() + + from matplotlib.patches import Circle + + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + i = 0 + cstart = cylinder_origin + (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], lw=3.0) + (tip_line,) = ax.plot(com[i, 0], com[i, 1], "k--") + + min_limits = np.array([0.0, 0.0, -0.5 * cylinder_height]) + max_limits = min_limits + cylinder_height + + ax.set_xlim([min_limits[0], max_limits[0]]) + ax.set_ylim([min_limits[2], max_limits[2]]) + + circle_artist = Circle((cstart[i, 0], cstart[i, 2]), cylinder_radius, color="g") + ax.add_artist(circle_artist) + ax.set_aspect("equal") + video_name = "2D_" + video_name + with writer.saving(fig, video_name, dpi): + with plt.style.context("fivethirtyeight"): + for i in range(0, positions.shape[0], int(step)): + rod_line.set_xdata(positions[i, 0]) + rod_line.set_ydata(positions[i, 2]) + tip_line.set_xdata(com[:i, 0]) + tip_line.set_ydata(com[:i, 2]) + circle_artist.center = cstart[i, 0], cstart[i, 2] + writer.grab_frame() plot_video(recorded_rod_history, recorded_cyl_history, "cylinder_rod_collision.mp4") From 2e635d8ada0cac7ad672a6352a22c86e952e6a37 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 11 Jun 2022 14:34:35 -0500 Subject: [PATCH 05/37] fix:capitalize flag names --- .../RodRigidBodyContact/rod_cylinder_contact_validation.py | 6 +++--- .../rod_cylinder_contact_with_y_normal.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py index 2fcc2025..ff297a2c 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py @@ -29,7 +29,7 @@ class SingleRodSingleCylinderInteractionSimulator( action_plane_key = "x" # can be set to True, checks collision at tips of rod -tip_collision = True +TIP_COLLISION = True _roll_key = 0 if action_plane_key == "x" else (1 if action_plane_key == "y" else 2) if action_plane_key == "x": @@ -69,7 +69,7 @@ class SingleRodSingleCylinderInteractionSimulator( # Cylinder surface starts at 0.2 tip_offset = 0.0 -if tip_collision: +if TIP_COLLISION: # The random choice decides which tip of the rod intersects with cylinder tip_choice = np.random.choice([1, -1]) tip_offset = 0.5 * tip_choice * base_length * np.cos(inclination) @@ -271,7 +271,7 @@ def plot_video( positions = np.array(recorded_rod_history["position"]) sim_time = np.array(recorded_rod_history["time"]) colliding_element_idx = n_elem // 2 - if tip_collision: + if TIP_COLLISION: colliding_element_idx = 0 if tip_choice == 1 else -1 colliding_element_history = positions[:, :, colliding_element_idx] fig = plt.figure(3, figsize=(8, 5)) diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py index 359fab9b..b8e22c32 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py @@ -42,12 +42,12 @@ class SingleRodSingleCylinderInteractionSimulator( cylinder_radius = 10.0 * base_radius # can be set to True, checks collision at tips of rod -tip_collision = False +TIP_COLLISION = False # Cylinder surface starts at 0.2 tip_offset = 0.0 -if tip_collision: +if TIP_COLLISION: # The random choice decides which tip of the rod intersects with cylinder tip_choice = np.random.choice([1, -1]) tip_offset = 0.5 * tip_choice * base_length * np.cos(inclination) @@ -248,7 +248,7 @@ def plot_video( positions = np.array(recorded_rod_history["position"]) sim_time = np.array(recorded_rod_history["time"]) colliding_element_idx = n_elem // 2 - if tip_collision: + if TIP_COLLISION: colliding_element_idx = -1 if tip_choice == 1 else 0 colliding_element_history = positions[:, :, colliding_element_idx] fig = plt.figure(3, figsize=(8, 5)) From 2997d100d09364aaff4c003d1e90c9eec39bd7b4 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 11 Jun 2022 15:12:47 -0500 Subject: [PATCH 06/37] move:post-processing scripts to post_processing.py --- .../RodRigidBodyContact/post_processing.py | 158 ++++++++++++++++++ .../rod_cylinder_contact_validation.py | 151 +++-------------- .../rod_cylinder_contact_with_y_normal.py | 149 ++--------------- 3 files changed, 195 insertions(+), 263 deletions(-) create mode 100644 examples/RigidBodyCases/RodRigidBodyContact/post_processing.py diff --git a/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py b/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py new file mode 100644 index 00000000..0c519698 --- /dev/null +++ b/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py @@ -0,0 +1,158 @@ +import numpy as np +from matplotlib import pyplot as plt + + +def make_data_for_cylinder_along_y(cstart, cradius, cheight): + center_x, center_z = cstart[0], cstart[1] + y = np.linspace(0, cheight, 5) + theta = np.linspace(0, 2 * np.pi, 20) + theta_grid, y_grid = np.meshgrid(theta, y) + x_grid = cradius * np.cos(theta_grid) + center_x + z_grid = cradius * np.sin(theta_grid) + center_z + y_grid += cstart[2] + return [x_grid, y_grid, z_grid] + + +def plot_video( + rod_history: dict, + cylinder_history: dict, + video_name="video.mp4", + margin=0.2, + fps=60, + step=1, + *args, + **kwargs +): # (time step, x/y/z, node) + + cylinder_start = np.array(cylinder_history["position"])[0, ...] + cylinder_radius = kwargs.get("cylinder_radius") + cylinder_height = kwargs.get("cylinder_height") + cylinder_direction = kwargs.get("cylinder_direction") + + XC, YC, ZC = make_data_for_cylinder_along_y( + cylinder_start, cylinder_radius, cylinder_height + ) + + import matplotlib.animation as manimation + + plt.rcParams.update({"font.size": 22}) + + # Should give a (n_time, 3, n_elem) array + positions = np.array(rod_history["position"]) + # (n_time, 3) array + com = np.array(rod_history["com"]) + + cylinder_com = np.array(cylinder_history["com"]) + cylinder_origin = cylinder_com - 0.5 * cylinder_height * cylinder_direction + + print("plot video") + FFMpegWriter = manimation.writers["ffmpeg"] + metadata = dict(title="Movie Test", artist="Matplotlib", comment="Movie support!") + writer = FFMpegWriter(fps=fps, metadata=metadata) + dpi = 50 + + # min_limits = np.roll(np.array([0.0, -0.5 * cylinder_height, 0.0]), _roll_key) + from mpl_toolkits.mplot3d import Axes3D + + fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) + ax = plt.axes(projection="3d") # fig.add_subplot(111) + ax.grid(b=True, which="minor", color="k", linestyle="--") + ax.grid(b=True, which="major", color="k", linestyle="-") + # plt.axis("square") + i = 0 + (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], positions[i, 2], lw=3.0) + XC, YC, ZC = make_data_for_cylinder_along_y( + cylinder_origin[i, ...], cylinder_radius, cylinder_height + ) + surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + + min_limits = np.array([0.0, 0.0, -0.5 * cylinder_height]) + min_limits = -np.abs(min_limits) + max_limits = min_limits + cylinder_height + + ax.set_xlim([min_limits[0], max_limits[0]]) + ax.set_ylim([min_limits[1], max_limits[1]]) + ax.set_zlim([min_limits[2], max_limits[2]]) + with writer.saving(fig, video_name, dpi): + with plt.style.context("seaborn-white"): + for i in range(0, positions.shape[0], int(step)): + rod_line.set_xdata(positions[i, 0]) + rod_line.set_ydata(positions[i, 1]) + rod_line.set_3d_properties(positions[i, 2]) + XC, YC, ZC = make_data_for_cylinder_along_y( + cylinder_origin[i, ...], cylinder_radius, cylinder_height + ) + surf.remove() + surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) + writer.grab_frame() + + from matplotlib.patches import Circle + + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + i = 0 + cstart = cylinder_origin + (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], lw=3.0) + (tip_line,) = ax.plot(com[i, 0], com[i, 1], "k--") + + min_limits = np.array([0.0, 0.0, -0.5 * cylinder_height]) + max_limits = min_limits + cylinder_height + + ax.set_xlim([min_limits[0], max_limits[0]]) + ax.set_ylim([min_limits[2], max_limits[2]]) + + circle_artist = Circle((cstart[i, 0], cstart[i, 2]), cylinder_radius, color="g") + ax.add_artist(circle_artist) + ax.set_aspect("equal") + video_name = "2D_" + video_name + with writer.saving(fig, video_name, dpi): + with plt.style.context("fivethirtyeight"): + for i in range(0, positions.shape[0], int(step)): + rod_line.set_xdata(positions[i, 0]) + rod_line.set_ydata(positions[i, 2]) + tip_line.set_xdata(com[:i, 0]) + tip_line.set_ydata(com[:i, 2]) + circle_artist.center = cstart[i, 0], cstart[i, 2] + writer.grab_frame() + + +def plot_cylinder_rod_position( + rod_history, + cylinder_history, + cylinder_radius, + rod_base_radius, + TIP_COLLISION, + TIP_CHOICE, + _roll_key=0, +): + cylinder_start = np.array(cylinder_history["position"])[0, ...] + positions = np.array(rod_history["position"]) + sim_time = np.array(rod_history["time"]) + + n_elem = positions.shape[-1] + + fig = plt.figure(figsize=(10, 8), frameon=True, dpi=150) + plt.rcParams.update({"font.size": 18}) + ax = fig.add_subplot(111) + colliding_element_idx = n_elem // 2 + if TIP_COLLISION: + colliding_element_idx = 0 if TIP_CHOICE == 1 else -1 + colliding_element_history = positions[:, :, colliding_element_idx] + # fig = plt.figure(3, figsize=(8, 5)) + # ax = fig.add_subplot(111) + ax.plot(sim_time, colliding_element_history[:, _roll_key], label="rod") + ax.hlines( + cylinder_start[_roll_key] - cylinder_radius - rod_base_radius, + sim_time[0], + sim_time[-1], + "k", + linestyle="dashed", + label="cylinder", + ) + plt.xlabel("Time [s]", fontsize=20) + plt.ylabel("Position", fontsize=20) + fig.legend(prop={"size": 20}) + plt.show() diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py index ff297a2c..ad538395 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_validation.py @@ -5,7 +5,7 @@ sys.path.append("../../../") from elastica import * -from matplotlib import pyplot as plt +from post_processing import plot_video, plot_cylinder_rod_position class SingleRodSingleCylinderInteractionSimulator( @@ -30,6 +30,7 @@ class SingleRodSingleCylinderInteractionSimulator( # can be set to True, checks collision at tips of rod TIP_COLLISION = True +TIP_CHOICE = 1 _roll_key = 0 if action_plane_key == "x" else (1 if action_plane_key == "y" else 2) if action_plane_key == "x": @@ -71,8 +72,8 @@ class SingleRodSingleCylinderInteractionSimulator( tip_offset = 0.0 if TIP_COLLISION: # The random choice decides which tip of the rod intersects with cylinder - tip_choice = np.random.choice([1, -1]) - tip_offset = 0.5 * tip_choice * base_length * np.cos(inclination) + TIP_CHOICE = np.random.choice([1, -1]) + tip_offset = 0.5 * TIP_CHOICE * base_length * np.cos(inclination) start_rod_1 = np.array( [ @@ -153,135 +154,21 @@ def make_callback(self, system, time, current_step: int): integrate(timestepper, single_rod_sim, final_time, total_steps) if PLOT_FIGURE: - - def make_data_for_cylinder_along_z(cstart, cradius, cheight): - center_x, center_y = cstart[0], cstart[1] - z = np.linspace(0, cheight, 5) - theta = np.linspace(0, 2 * np.pi, 20) - theta_grid, z_grid = np.meshgrid(theta, z) - x_grid = cradius * np.cos(theta_grid) + center_x - y_grid = cradius * np.sin(theta_grid) + center_y - z_grid += cstart[2] - return np.roll([x_grid, y_grid, z_grid], _roll_key) - - XC, YC, ZC = make_data_for_cylinder_along_z( - cylinder_start, cylinder_radius, cylinder_height + plot_video( + recorded_rod_history, + recorded_cyl_history, + "cylinder_rod_collision.mp4", + cylinder_direction=cylinder_direction, + cylinder_height=cylinder_height, + cylinder_radius=cylinder_radius, ) - def plot_video( - plot_params: dict, - cylinder_history: dict, - video_name="video.mp4", - margin=0.2, - fps=60, - step=1, - *args, - **kwargs - ): # (time step, x/y/z, node) - import matplotlib.animation as manimation - - plt.rcParams.update({"font.size": 22}) - - # Should give a (n_time, 3, n_elem) array - positions = np.array(plot_params["position"]) - # (n_time, 3) array - com = np.array(plot_params["com"]) - - cylinder_com = np.array(cylinder_history["com"]) - cylinder_origin = cylinder_com - 0.5 * cylinder_height * cylinder_direction - - print("plot video") - FFMpegWriter = manimation.writers["ffmpeg"] - metadata = dict( - title="Movie Test", artist="Matplotlib", comment="Movie support!" - ) - writer = FFMpegWriter(fps=fps, metadata=metadata) - dpi = 50 - - # min_limits = np.roll(np.array([0.0, -0.5 * cylinder_height, 0.0]), _roll_key) - from mpl_toolkits.mplot3d import Axes3D - - fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) - ax = plt.axes(projection="3d") # fig.add_subplot(111) - ax.grid(b=True, which="minor", color="k", linestyle="--") - ax.grid(b=True, which="major", color="k", linestyle="-") - # plt.axis("square") - i = 0 - (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], positions[i, 2], lw=3.0) - XC, YC, ZC = make_data_for_cylinder_along_z( - cylinder_origin[i, ...], cylinder_radius, cylinder_height - ) - surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_zlabel("z") - - min_limits = global_rot_mat @ np.array([0.0, -0.5 * cylinder_height, 0.0]) - min_limits = -np.abs(min_limits) - max_limits = min_limits + cylinder_height - - ax.set_xlim([min_limits[0], max_limits[0]]) - ax.set_ylim([min_limits[1], max_limits[1]]) - ax.set_zlim([min_limits[2], max_limits[2]]) - with writer.saving(fig, video_name, dpi): - with plt.style.context("seaborn-white"): - for i in range(0, positions.shape[0], int(step)): - rod_line.set_xdata(positions[i, 0]) - rod_line.set_ydata(positions[i, 1]) - rod_line.set_3d_properties(positions[i, 2]) - XC, YC, ZC = make_data_for_cylinder_along_z( - cylinder_origin[i, ...], cylinder_radius, cylinder_height - ) - surf.remove() - surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) - writer.grab_frame() - from matplotlib.patches import Circle - - fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) - ax = fig.add_subplot(111) - i = 0 - positions = np.roll(positions, -_roll_key, axis=1) - com = np.roll(com, -_roll_key, axis=1) - cstart = np.roll(cylinder_origin, -_roll_key, axis=1) - (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], lw=3.0) - (tip_line,) = ax.plot(com[i, 0], com[i, 1], "k--") - - min_limits = np.array([0.0, -0.5 * cylinder_height, 0.0]) - max_limits = min_limits + cylinder_height - - ax.set_xlim([min_limits[0], max_limits[0]]) - ax.set_ylim([min_limits[1], max_limits[1]]) - - circle_artist = Circle((cstart[i, 0], cstart[i, 1]), cylinder_radius, color="g") - ax.add_artist(circle_artist) - ax.set_aspect("equal") - video_name = "2D_" + video_name - with writer.saving(fig, video_name, dpi): - with plt.style.context("fivethirtyeight"): - for i in range(0, positions.shape[0], int(step)): - rod_line.set_xdata(positions[i, 0]) - rod_line.set_ydata(positions[i, 1]) - tip_line.set_xdata(com[:i, 0]) - tip_line.set_ydata(com[:i, 1]) - circle_artist.center = cstart[i, 0], cstart[i, 1] - writer.grab_frame() - - plot_video(recorded_rod_history, recorded_cyl_history, "cylinder_rod_collision.mp4") - - positions = np.array(recorded_rod_history["position"]) - sim_time = np.array(recorded_rod_history["time"]) - colliding_element_idx = n_elem // 2 - if TIP_COLLISION: - colliding_element_idx = 0 if tip_choice == 1 else -1 - colliding_element_history = positions[:, :, colliding_element_idx] - fig = plt.figure(3, figsize=(8, 5)) - ax = fig.add_subplot(111) - ax.plot(sim_time, colliding_element_history[:, _roll_key]) - ax.hlines( - cylinder_start[_roll_key] - cylinder_radius - base_radius, - sim_time[0], - sim_time[-1], - "k", - linestyle="dashed", + plot_cylinder_rod_position( + recorded_rod_history, + recorded_cyl_history, + cylinder_radius=cylinder_radius, + rod_base_radius=base_radius, + TIP_COLLISION=TIP_COLLISION, + TIP_CHOICE=TIP_CHOICE, + _roll_key=_roll_key, ) - plt.show() diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py index b8e22c32..e001c062 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_with_y_normal.py @@ -5,7 +5,7 @@ sys.path.append("../../../") from elastica import * -from matplotlib import pyplot as plt +from post_processing import plot_video, plot_cylinder_rod_position class SingleRodSingleCylinderInteractionSimulator( @@ -47,10 +47,11 @@ class SingleRodSingleCylinderInteractionSimulator( # Cylinder surface starts at 0.2 tip_offset = 0.0 +TIP_CHOICE = 1 if TIP_COLLISION: # The random choice decides which tip of the rod intersects with cylinder - tip_choice = np.random.choice([1, -1]) - tip_offset = 0.5 * tip_choice * base_length * np.cos(inclination) + TIP_CHOICE = np.random.choice([1, -1]) + tip_offset = 0.5 * TIP_CHOICE * base_length * np.cos(inclination) start_rod_1 = np.array( [ @@ -131,134 +132,20 @@ def make_callback(self, system, time, current_step: int): integrate(timestepper, single_rod_sim, final_time, total_steps) if PLOT_FIGURE: - - def make_data_for_cylinder_along_y(cstart, cradius, cheight): - center_x, center_z = cstart[0], cstart[1] - y = np.linspace(0, cheight, 5) - theta = np.linspace(0, 2 * np.pi, 20) - theta_grid, y_grid = np.meshgrid(theta, y) - x_grid = cradius * np.cos(theta_grid) + center_x - z_grid = cradius * np.sin(theta_grid) + center_z - y_grid += cstart[2] - return [x_grid, y_grid, z_grid] - - XC, YC, ZC = make_data_for_cylinder_along_y( - cylinder_start, cylinder_radius, cylinder_height + plot_video( + recorded_rod_history, + recorded_cyl_history, + "cylinder_rod_collision.mp4", + cylinder_direction=cylinder_direction, + cylinder_height=cylinder_height, + cylinder_radius=cylinder_radius, ) - def plot_video( - plot_params: dict, - cylinder_history: dict, - video_name="video.mp4", - margin=0.2, - fps=60, - step=1, - *args, - **kwargs - ): # (time step, x/y/z, node) - import matplotlib.animation as manimation - - plt.rcParams.update({"font.size": 22}) - - # Should give a (n_time, 3, n_elem) array - positions = np.array(plot_params["position"]) - # (n_time, 3) array - com = np.array(plot_params["com"]) - - cylinder_com = np.array(cylinder_history["com"]) - cylinder_origin = cylinder_com - 0.5 * cylinder_height * cylinder_direction - - print("plot video") - FFMpegWriter = manimation.writers["ffmpeg"] - metadata = dict( - title="Movie Test", artist="Matplotlib", comment="Movie support!" - ) - writer = FFMpegWriter(fps=fps, metadata=metadata) - dpi = 50 - - # min_limits = np.roll(np.array([0.0, -0.5 * cylinder_height, 0.0]), _roll_key) - from mpl_toolkits.mplot3d import Axes3D - - fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) - ax = plt.axes(projection="3d") # fig.add_subplot(111) - ax.grid(b=True, which="minor", color="k", linestyle="--") - ax.grid(b=True, which="major", color="k", linestyle="-") - # plt.axis("square") - i = 0 - (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], positions[i, 2], lw=3.0) - XC, YC, ZC = make_data_for_cylinder_along_y( - cylinder_origin[i, ...], cylinder_radius, cylinder_height - ) - surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_zlabel("z") - - min_limits = np.array([0.0, 0.0, -0.5 * cylinder_height]) - min_limits = -np.abs(min_limits) - max_limits = min_limits + cylinder_height - - ax.set_xlim([min_limits[0], max_limits[0]]) - ax.set_ylim([min_limits[1], max_limits[1]]) - ax.set_zlim([min_limits[2], max_limits[2]]) - with writer.saving(fig, video_name, dpi): - with plt.style.context("seaborn-white"): - for i in range(0, positions.shape[0], int(step)): - rod_line.set_xdata(positions[i, 0]) - rod_line.set_ydata(positions[i, 1]) - rod_line.set_3d_properties(positions[i, 2]) - XC, YC, ZC = make_data_for_cylinder_along_y( - cylinder_origin[i, ...], cylinder_radius, cylinder_height - ) - surf.remove() - surf = ax.plot_surface(XC, YC, ZC, color="g", alpha=0.5) - writer.grab_frame() - - from matplotlib.patches import Circle - - fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) - ax = fig.add_subplot(111) - i = 0 - cstart = cylinder_origin - (rod_line,) = ax.plot(positions[i, 0], positions[i, 1], lw=3.0) - (tip_line,) = ax.plot(com[i, 0], com[i, 1], "k--") - - min_limits = np.array([0.0, 0.0, -0.5 * cylinder_height]) - max_limits = min_limits + cylinder_height - - ax.set_xlim([min_limits[0], max_limits[0]]) - ax.set_ylim([min_limits[2], max_limits[2]]) - - circle_artist = Circle((cstart[i, 0], cstart[i, 2]), cylinder_radius, color="g") - ax.add_artist(circle_artist) - ax.set_aspect("equal") - video_name = "2D_" + video_name - with writer.saving(fig, video_name, dpi): - with plt.style.context("fivethirtyeight"): - for i in range(0, positions.shape[0], int(step)): - rod_line.set_xdata(positions[i, 0]) - rod_line.set_ydata(positions[i, 2]) - tip_line.set_xdata(com[:i, 0]) - tip_line.set_ydata(com[:i, 2]) - circle_artist.center = cstart[i, 0], cstart[i, 2] - writer.grab_frame() - - plot_video(recorded_rod_history, recorded_cyl_history, "cylinder_rod_collision.mp4") - - positions = np.array(recorded_rod_history["position"]) - sim_time = np.array(recorded_rod_history["time"]) - colliding_element_idx = n_elem // 2 - if TIP_COLLISION: - colliding_element_idx = -1 if tip_choice == 1 else 0 - colliding_element_history = positions[:, :, colliding_element_idx] - fig = plt.figure(3, figsize=(8, 5)) - ax = fig.add_subplot(111) - ax.plot(sim_time, colliding_element_history[:, 0]) - ax.hlines( - cylinder_start[0] - cylinder_radius - base_radius, - sim_time[0], - sim_time[-1], - "k", - linestyle="dashed", + plot_cylinder_rod_position( + recorded_rod_history, + recorded_cyl_history, + cylinder_radius=cylinder_radius, + rod_base_radius=base_radius, + TIP_COLLISION=TIP_COLLISION, + TIP_CHOICE=TIP_CHOICE, ) - plt.show() From 44a441777cb25c76c5d3fcc312107d8d84d84733 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 11 Jun 2022 15:34:15 -0500 Subject: [PATCH 07/37] refactor EndPointForcesSinusoidal example and test cases --- elastica/external_forces.py | 71 ++++++++++++++++++ .../external_force_class_for_joint_test.py | 72 ------------------- examples/JointCases/fixed_joint.py | 3 - examples/JointCases/hinge_joint.py | 3 - examples/JointCases/spherical_joint.py | 3 - tests/test_external_forces.py | 4 +- 6 files changed, 72 insertions(+), 84 deletions(-) delete mode 100644 examples/JointCases/external_force_class_for_joint_test.py diff --git a/elastica/external_forces.py b/elastica/external_forces.py index b0de600d..a3a44919 100644 --- a/elastica/external_forces.py +++ b/elastica/external_forces.py @@ -6,6 +6,7 @@ "UniformTorques", "UniformForces", "MuscleTorques", + "EndpointForcesSinusoidal", ] import numpy as np @@ -470,3 +471,73 @@ def inplace_substraction(external_force_or_torque, force_or_torque): for i in range(3): for k in range(blocksize): external_force_or_torque[i, k] -= force_or_torque[i, k] + + +# Force class. This force class is used only for joint test cases +class EndpointForcesSinusoidal(NoForces): + """Applies sinusoidal forces on endpoints""" + + def __init__( + self, + start_force_mag, + end_force_mag, + ramp_up_time=0.0, + tangent_direction=np.array([0, 0, 1]), + normal_direction=np.array([0, 1, 0]), + ): + super(EndpointForcesSinusoidal, self).__init__() + # Start force + self.start_force_mag = start_force_mag + self.end_force_mag = end_force_mag + + # Applied force directions + self.normal_direction = normal_direction + # self.roll_direction = np.cross(tangent_direction, normal_direction) + self.roll_direction = np.cross(normal_direction, tangent_direction) + + assert ramp_up_time >= 0.0 + self.ramp_up_time = ramp_up_time + + def apply_forces(self, system, time=0.0): + + if time < self.ramp_up_time: + # When time smaller than ramp up time apply the force in normal direction + # First pull the rod upward or downward direction some time. + start_force = -2.0 * self.start_force_mag * self.normal_direction + end_force = -2.0 * self.end_force_mag * self.normal_direction + + system.external_forces[..., 0] += start_force + system.external_forces[..., -1] += end_force + + else: + # When time is greater than ramp up time, forces are applied in normal + # and roll direction or forces are in a plane perpendicular to the + # direction. + + # First force applied to start of the rod + roll_forces_start = ( + self.start_force_mag + * np.cos(0.5 * np.pi * (time - self.ramp_up_time)) + * self.roll_direction + ) + normal_forces_start = ( + self.start_force_mag + * np.sin(0.5 * np.pi * (time - self.ramp_up_time)) + * self.normal_direction + ) + start_force = roll_forces_start + normal_forces_start + # Now force applied to end of the rod + roll_forces_end = ( + self.end_force_mag + * np.cos(0.5 * np.pi * (time - self.ramp_up_time)) + * self.roll_direction + ) + normal_forces_end = ( + self.end_force_mag + * np.sin(0.5 * np.pi * (time - self.ramp_up_time)) + * self.normal_direction + ) + end_force = roll_forces_end + normal_forces_end + # Update external forces + system.external_forces[..., 0] += start_force + system.external_forces[..., -1] += end_force diff --git a/examples/JointCases/external_force_class_for_joint_test.py b/examples/JointCases/external_force_class_for_joint_test.py deleted file mode 100644 index 3fa31877..00000000 --- a/examples/JointCases/external_force_class_for_joint_test.py +++ /dev/null @@ -1,72 +0,0 @@ -import numpy as np -from elastica.external_forces import NoForces - - -# Force class. This force class is used only for joint test cases -class EndpointForcesSinusoidal(NoForces): - """Applies sinusoidal forces on endpoints""" - - def __init__( - self, - start_force_mag, - end_force_mag, - ramp_up_time=0.0, - tangent_direction=np.array([0, 0, 1]), - normal_direction=np.array([0, 1, 0]), - ): - super(EndpointForcesSinusoidal, self).__init__() - # Start force - self.start_force_mag = start_force_mag - self.end_force_mag = end_force_mag - - # Applied force directions - self.normal_direction = normal_direction - # self.roll_direction = np.cross(tangent_direction, normal_direction) - self.roll_direction = np.cross(normal_direction, tangent_direction) - - assert ramp_up_time >= 0.0 - self.ramp_up_time = ramp_up_time - - def apply_forces(self, system, time=0.0): - - if time < self.ramp_up_time: - # When time smaller than ramp up time apply the force in normal direction - # First pull the rod upward or downward direction some time. - start_force = -2.0 * self.start_force_mag * self.normal_direction - end_force = -2.0 * self.end_force_mag * self.normal_direction - - system.external_forces[..., 0] += start_force - system.external_forces[..., -1] += end_force - - else: - # When time is greater than ramp up time, forces are applied in normal - # and roll direction or forces are in a plane perpendicular to the - # direction. - - # First force applied to start of the rod - roll_forces_start = ( - self.start_force_mag - * np.cos(0.5 * np.pi * (time - self.ramp_up_time)) - * self.roll_direction - ) - normal_forces_start = ( - self.start_force_mag - * np.sin(0.5 * np.pi * (time - self.ramp_up_time)) - * self.normal_direction - ) - start_force = roll_forces_start + normal_forces_start - # Now force applied to end of the rod - roll_forces_end = ( - self.end_force_mag - * np.cos(0.5 * np.pi * (time - self.ramp_up_time)) - * self.roll_direction - ) - normal_forces_end = ( - self.end_force_mag - * np.sin(0.5 * np.pi * (time - self.ramp_up_time)) - * self.normal_direction - ) - end_force = roll_forces_end + normal_forces_end - # Update external forces - system.external_forces[..., 0] += start_force - system.external_forces[..., -1] += end_force diff --git a/examples/JointCases/fixed_joint.py b/examples/JointCases/fixed_joint.py index e41d434f..eaad4134 100644 --- a/examples/JointCases/fixed_joint.py +++ b/examples/JointCases/fixed_joint.py @@ -6,9 +6,6 @@ # FIXME without appending sys.path make it more generic sys.path.append("../../") from elastica import * -from examples.JointCases.external_force_class_for_joint_test import ( - EndpointForcesSinusoidal, -) from examples.JointCases.joint_cases_postprocessing import ( plot_position, plot_video, diff --git a/examples/JointCases/hinge_joint.py b/examples/JointCases/hinge_joint.py index 19f5da99..82f1b399 100644 --- a/examples/JointCases/hinge_joint.py +++ b/examples/JointCases/hinge_joint.py @@ -6,9 +6,6 @@ # FIXME without appending sys.path make it more generic sys.path.append("../../") from elastica import * -from examples.JointCases.external_force_class_for_joint_test import ( - EndpointForcesSinusoidal, -) from examples.JointCases.joint_cases_postprocessing import ( plot_position, plot_video, diff --git a/examples/JointCases/spherical_joint.py b/examples/JointCases/spherical_joint.py index fce947f5..bf7b88dc 100644 --- a/examples/JointCases/spherical_joint.py +++ b/examples/JointCases/spherical_joint.py @@ -7,9 +7,6 @@ # FIXME without appending sys.path make it more generic sys.path.append("../../") from elastica import * -from examples.JointCases.external_force_class_for_joint_test import ( - EndpointForcesSinusoidal, -) from examples.JointCases.joint_cases_postprocessing import ( plot_position, plot_video, diff --git a/tests/test_external_forces.py b/tests/test_external_forces.py index 67a80f98..417bda8a 100644 --- a/tests/test_external_forces.py +++ b/tests/test_external_forces.py @@ -15,11 +15,9 @@ MuscleTorques, inplace_addition, inplace_substraction, -) -from elastica.utils import Tolerance -from examples.JointCases.external_force_class_for_joint_test import ( EndpointForcesSinusoidal, ) +from elastica.utils import Tolerance def mock_rod_init(self): From 82d1dac7d5a2b89102b2a978b876dbebbf63e337 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 11 Jun 2022 22:04:39 -0500 Subject: [PATCH 08/37] enhancement:rod-rod contact compute total velocity --- examples/RodContactCase/post_processing.py | 23 ++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/examples/RodContactCase/post_processing.py b/examples/RodContactCase/post_processing.py index 93c5f349..45e27b04 100644 --- a/examples/RodContactCase/post_processing.py +++ b/examples/RodContactCase/post_processing.py @@ -330,6 +330,13 @@ def plot_velocity( axs[0].plot(time[:], avg_velocity_rod_one[:, 0], linewidth=3, label="rod_one") axs[0].plot(time[:], avg_velocity_rod_two[:, 0], linewidth=3, label="rod_two") + axs[0].plot( + time[:], + avg_velocity_rod_one[:, 0] + avg_velocity_rod_two[:, 0], + "--", + linewidth=3, + label="total", + ) axs[0].set_ylabel("x velocity", fontsize=20) axs[1].plot( @@ -342,6 +349,12 @@ def plot_velocity( avg_velocity_rod_two[:, 1], linewidth=3, ) + axs[1].plot( + time[:], + avg_velocity_rod_one[:, 1] + avg_velocity_rod_two[:, 1], + "--", + linewidth=3, + ) axs[1].set_ylabel("y velocity", fontsize=20) axs[2].plot( @@ -354,6 +367,12 @@ def plot_velocity( avg_velocity_rod_two[:, 2], linewidth=3, ) + axs[2].plot( + time[:], + avg_velocity_rod_one[:, 2] + avg_velocity_rod_two[:, 2], + "--", + linewidth=3, + ) axs[2].set_ylabel("z velocity", fontsize=20) axs[3].semilogy( @@ -383,8 +402,8 @@ def plot_velocity( plt.close(plt.gcf()) if SAVE_FIGURE: - # fig.savefig(filename) - plt.savefig(filename) + fig.savefig(filename) + # plt.savefig(filename) def plot_link_writhe_twist(twist_density, total_twist, total_writhe, total_link): From d3f91f2429f9640fbe707dee1a5bd0f584167bb6 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 11 Jun 2022 22:05:07 -0500 Subject: [PATCH 09/37] hotfix:rod-rod contact examples remove poisson ratio from init This commit removes poisson ratio from the initialization of rods in rod-rod contact examples. Since initializing poisson ratio is depreciated. --- .../RodRodContact/rod_rod_contact_inclined_validation.py | 3 +-- .../RodRodContact/rod_rod_contact_parallel_validation.py | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py index bb4e13e0..4c364c23 100644 --- a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py +++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py @@ -30,7 +30,7 @@ class InclinedRodRodContact( base_radius = 0.01 base_area = np.pi * base_radius ** 2 density = 1750 -nu = 0.001 +nu = 0.0 E = 3e5 poisson_ratio = 0.5 shear_modulus = E / (poisson_ratio + 1.0) @@ -58,7 +58,6 @@ class InclinedRodRodContact( density, nu, E, - poisson_ratio, shear_modulus=shear_modulus, ) diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py index e5f81969..6f983ab2 100644 --- a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py +++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py @@ -30,7 +30,7 @@ class ParallelRodRodContact( base_radius = 0.01 base_area = np.pi * base_radius ** 2 density = 1750 -nu = 0.001 +nu = 0.0 E = 3e5 poisson_ratio = 0.5 shear_modulus = E / (poisson_ratio + 1.0) From e52ec2b12eb0feb94f49cf1a1da382ba6aee0b30 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 11 Jun 2022 22:45:27 -0500 Subject: [PATCH 10/37] bugfix:remove poisson ratio from initialization --- examples/HelicalBucklingCase/convergence_helicalbuckling.py | 2 +- examples/MuscularFlagella/muscular_flagella.py | 5 +++-- .../ContinuumSnakeVisualization/continuum_snake.py | 3 ++- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/examples/HelicalBucklingCase/convergence_helicalbuckling.py b/examples/HelicalBucklingCase/convergence_helicalbuckling.py index ca728381..44b94944 100644 --- a/examples/HelicalBucklingCase/convergence_helicalbuckling.py +++ b/examples/HelicalBucklingCase/convergence_helicalbuckling.py @@ -45,6 +45,7 @@ def simulate_helicalbucklin_beam_with( number_of_rotations = 27 # For shear modulus of 1e4, nu is 99! poisson_ratio = 99 + shear_modulus = E / (poisson_ratio + 1.0) shear_matrix = np.repeat(1e5 * np.identity((3))[:, :, np.newaxis], n_elem, axis=2) temp_bend_matrix = np.zeros((3, 3)) np.fill_diagonal(temp_bend_matrix, [1.345, 1.345, 0.789]) @@ -60,7 +61,6 @@ def simulate_helicalbucklin_beam_with( density, nu, E, - poisson_ratio, ) # TODO: CosseratRod has to be able to take shear matrix as input, we should change it as done below diff --git a/examples/MuscularFlagella/muscular_flagella.py b/examples/MuscularFlagella/muscular_flagella.py index 427feb08..a9e49220 100644 --- a/examples/MuscularFlagella/muscular_flagella.py +++ b/examples/MuscularFlagella/muscular_flagella.py @@ -40,6 +40,7 @@ class MuscularFlagellaSimulator( base_length_body = 1.927 # mm E = 3.86e6 # MPa poisson_ratio = 0.5 +shear_modulus = E / (poisson_ratio + 1.0) base_radius_head = 0.02 # mm base_radius_tail = 0.007 # mm @@ -67,7 +68,7 @@ class MuscularFlagellaSimulator( density_body, nu_body, E, - poisson_ratio, + shear_modulus=shear_modulus, ) # In order to match bending stiffness of the tail as given in below reference, recompute and @@ -152,7 +153,7 @@ class MuscularFlagellaSimulator( density_muscle, nu_muscle, E_muscle, - poisson_ratio, + shear_modulus=shear_modulus, ) muscular_flagella_sim.append(flagella_muscle) diff --git a/examples/Visualization/ContinuumSnakeVisualization/continuum_snake.py b/examples/Visualization/ContinuumSnakeVisualization/continuum_snake.py index d76da61c..46695225 100644 --- a/examples/Visualization/ContinuumSnakeVisualization/continuum_snake.py +++ b/examples/Visualization/ContinuumSnakeVisualization/continuum_snake.py @@ -36,6 +36,7 @@ def run_snake(b_coeff, SAVE_RESULTS=False): nu = 5.0 E = 1e7 poisson_ratio = 0.5 + shear_modulus = E / (poisson_ratio + 1.0) shearable_rod = CosseratRod.straight_rod( n_elem, @@ -47,7 +48,7 @@ def run_snake(b_coeff, SAVE_RESULTS=False): density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) snake_sim.append(shearable_rod) From ce5e12acf8c7df6aadc8ce7867e525b38e1ab1a0 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 12 Jun 2022 13:55:42 -0500 Subject: [PATCH 11/37] hotfix:rod-rod contact post-processing vis2D flag fix --- examples/RodContactCase/post_processing.py | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/RodContactCase/post_processing.py b/examples/RodContactCase/post_processing.py index 45e27b04..2ebf30ab 100644 --- a/examples/RodContactCase/post_processing.py +++ b/examples/RodContactCase/post_processing.py @@ -12,7 +12,6 @@ def plot_video_with_surface( video_name="video.mp4", fps=60, step=1, - vis2D=True, **kwargs, ): plt.rcParams.update({"font.size": 22}) From 312a7066e1024946489bbd1f70e76360f6ce0898 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 12 Jun 2022 19:49:06 -0500 Subject: [PATCH 12/37] docs: documentation is added for Endpointforcessinusoidal --- elastica/external_forces.py | 44 ++++++++++++++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 3 deletions(-) diff --git a/elastica/external_forces.py b/elastica/external_forces.py index a3a44919..a2a2f801 100644 --- a/elastica/external_forces.py +++ b/elastica/external_forces.py @@ -473,9 +473,31 @@ def inplace_substraction(external_force_or_torque, force_or_torque): external_force_or_torque[i, k] -= force_or_torque[i, k] -# Force class. This force class is used only for joint test cases class EndpointForcesSinusoidal(NoForces): - """Applies sinusoidal forces on endpoints""" + """ + This class applies sinusoidally varying forces to the ends of a rod. + Forces are applied in a plane, which is defined by the tangent_direction and normal_direction. + + Attributes + ---------- + start_force_mag: float + Magnitude of the force that is applied to the start of the rod (node 0). + end_force_mag: float + Magnitude of the force that is applied to the end of the rod (node -1). + ramp_up_time: float + Applied forces are applied in the normal direction until time reaches ramp_up_time. + normal_direction: np.ndarray + An array (3,) contains type float. + This is the normal direction of the rod. + roll_direction: np.ndarray + An array (3,) contains type float. + This is the direction perpendicular to rod tangent, and rod normal. + + Notes + ----- + In order to see example how to use this class, see joint examples. + + """ def __init__( self, @@ -485,6 +507,23 @@ def __init__( tangent_direction=np.array([0, 0, 1]), normal_direction=np.array([0, 1, 0]), ): + """ + + Parameters + ---------- + start_force_mag: float + Magnitude of the force that is applied to the start of the rod (node 0). + end_force_mag: float + Magnitude of the force that is applied to the end of the rod (node -1). + ramp_up_time: float + Applied muscle torques are ramped up until ramp up time. + tangent_direction: np.ndarray + An array (3,) contains type float. + This is the tangent direction of the rod, or normal of the plane that forces applied. + normal_direction: np.ndarray + An array (3,) contains type float. + This is the normal direction of the rod. + """ super(EndpointForcesSinusoidal, self).__init__() # Start force self.start_force_mag = start_force_mag @@ -492,7 +531,6 @@ def __init__( # Applied force directions self.normal_direction = normal_direction - # self.roll_direction = np.cross(tangent_direction, normal_direction) self.roll_direction = np.cross(normal_direction, tangent_direction) assert ramp_up_time >= 0.0 From b6fb6f723b0c2a3c0b2c707a7da484195530f045 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Mon, 13 Jun 2022 14:47:23 -0500 Subject: [PATCH 13/37] update:documentation external_forces.rst --- docs/api/external_forces.rst | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/docs/api/external_forces.rst b/docs/api/external_forces.rst index fdb55b39..8f84f66f 100644 --- a/docs/api/external_forces.rst +++ b/docs/api/external_forces.rst @@ -22,6 +22,7 @@ External force and environmental interaction are represented as force/torque bou UniformForces UniformTorques MuscleTorques + EndpointForcesSinusoidal .. rubric:: Available Interaction @@ -40,11 +41,12 @@ Compatibility Forcing Rod Rigid Body ========================== ======= ============ NoForces ✅ ✅ -EndpointForces ✅ ✅ +EndpointForces ✅ ❌ GravityForces ✅ ✅ UniformForces ✅ ✅ UniformTorques ✅ ✅ MuscleTorques ✅ ❌ +EndpointForcesSinusoidal ✅ ❌ ========================== ======= ============ ========================== ======= ============ @@ -78,6 +80,9 @@ Built-in External Forces .. autoclass:: MuscleTorques :special-members: __init__ +.. autoclass:: EndpointForcesSinusoidal + :special-members: __init__ + Built-in Environment Interactions --------------------------------- .. automodule:: elastica.interaction From 6becfdf4e2be05de54d75bdec5ae34b0d8bed99c Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Thu, 23 Jun 2022 23:36:12 -0500 Subject: [PATCH 14/37] enhancement:rod-cylinder contact add contact torques on cylinder This commit fixes the rod-cylinder contact by adding torques due to the contact force on cylidner. Previously, we were only dealing with forces acting on cylinders however any contact force acting away from the cg will create torque. --- elastica/joint.py | 68 +++++++++++++++++++++++++++++++++++------------ 1 file changed, 51 insertions(+), 17 deletions(-) diff --git a/elastica/joint.py b/elastica/joint.py index 0f9eb841..98b37d0d 100644 --- a/elastica/joint.py +++ b/elastica/joint.py @@ -355,20 +355,24 @@ def _find_min_dist(x1, e1, x2, e2): s = potential_s t = 1.0 - return x2 + s * e2 - x1 - t * e1 + # Return distance, contact point of system 2, contact point of system 1 + return x2 + s * e2 - x1 - t * e1, x2 + s * e2, x1 - t * e1 @numba.njit(cache=True) def _calculate_contact_forces_rod_rigid_body( x_collection_rod, edge_collection_rod, - x_cylinder, + x_cylinder_center, + x_cylinder_tip, edge_cylinder, radii_sum, length_sum, internal_forces_rod, external_forces_rod, external_forces_cylinder, + external_torques_cylinder, + cylinder_director_collection, velocity_rod, velocity_cylinder, contact_k, @@ -376,11 +380,13 @@ def _calculate_contact_forces_rod_rigid_body( ): # We already pass in only the first n_elem x n_points = x_collection_rod.shape[1] + cylinder_total_contact_forces = np.zeros((3)) + cylinder_total_contact_torques = np.zeros((3)) for i in range(n_points): # Element-wise bounding box x_selected = x_collection_rod[..., i] # x_cylinder is already a (,) array from outised - del_x = x_selected - x_cylinder + del_x = x_selected - x_cylinder_tip norm_del_x = _norm(del_x) # If outside then don't process @@ -389,8 +395,8 @@ def _calculate_contact_forces_rod_rigid_body( # find the shortest line segment between the two centerline # segments : differs from normal cylinder-cylinder intersection - distance_vector = _find_min_dist( - x_selected, edge_collection_rod[..., i], x_cylinder, edge_cylinder + distance_vector, x_cylinder_contact_point, _ = _find_min_dist( + x_selected, edge_collection_rod[..., i], x_cylinder_tip, edge_cylinder ) distance_vector_length = _norm(distance_vector) distance_vector /= distance_vector_length @@ -428,22 +434,43 @@ def _calculate_contact_forces_rod_rigid_body( # magnitude* direction net_contact_force = ( - normal_force + 0.5 * mask * (contact_damping_force + contact_force) - ) * distance_vector + 0.5 + * mask + * (normal_force + (contact_damping_force + contact_force)) + * distance_vector + ) + + # Torques acting on the cylinder + moment_arm = x_cylinder_contact_point - x_cylinder_center # Add it to the rods at the end of the day if i == 0: - external_forces_rod[..., i] -= 0.5 * net_contact_force - external_forces_rod[..., i + 1] -= net_contact_force - external_forces_cylinder[..., 0] += 1.5 * net_contact_force + external_forces_rod[..., i] -= 2 / 3 * net_contact_force + external_forces_rod[..., i + 1] -= 4 / 3 * net_contact_force + cylinder_total_contact_forces += 2.0 * net_contact_force + cylinder_total_contact_torques += np.cross( + moment_arm, 2.0 * net_contact_force + ) elif i == n_points: - external_forces_rod[..., i] -= net_contact_force - external_forces_rod[..., i + 1] -= 0.5 * net_contact_force - external_forces_cylinder[..., 0] += 1.5 * net_contact_force + external_forces_rod[..., i] -= 4 / 3 * net_contact_force + external_forces_rod[..., i + 1] -= 2 / 3 * net_contact_force + cylinder_total_contact_forces += 2.0 * net_contact_force + cylinder_total_contact_torques += np.cross( + moment_arm, 2.0 * net_contact_force + ) else: external_forces_rod[..., i] -= net_contact_force external_forces_rod[..., i + 1] -= net_contact_force - external_forces_cylinder[..., 0] += 2.0 * net_contact_force + cylinder_total_contact_forces += 2.0 * net_contact_force + cylinder_total_contact_torques += np.cross( + moment_arm, 2.0 * net_contact_force + ) + + # Update the cylinder external forces and torques + external_forces_cylinder[..., 0] += cylinder_total_contact_forces + external_torques_cylinder[..., 0] += ( + cylinder_director_collection @ cylinder_total_contact_torques + ) @numba.njit(cache=True) @@ -488,7 +515,7 @@ def _calculate_contact_forces_rod_rod( # find the shortest line segment between the two centerline # segments : differs from normal cylinder-cylinder intersection - distance_vector = _find_min_dist( + distance_vector, _, _ = _find_min_dist( x_selected_rod_one, edge_collection_rod_one[..., i], x_selected_rod_two, @@ -596,7 +623,7 @@ def _calculate_contact_forces_self_rod( # find the shortest line segment between the two centerline # segments : differs from normal cylinder-cylinder intersection - distance_vector = _find_min_dist( + distance_vector, _, _ = _find_min_dist( x_selected_rod_index_i, edge_collection_rod_one[..., i], x_selected_rod_index_j, @@ -787,9 +814,14 @@ def apply_forces(self, rod_one, index_one, rod_two, index_two): - 0.5 * cylinder_two.length * cylinder_two.director_collection[2, :, 0] ) + rod_element_position = 0.5 * ( + rod_one.position_collection[..., 1:] + + rod_one.position_collection[..., :-1] + ) _calculate_contact_forces_rod_rigid_body( - rod_one.position_collection[..., :-1], + rod_element_position, rod_one.lengths * rod_one.tangents, + cylinder_two.position_collection[..., 0], x_cyl, cylinder_two.length * cylinder_two.director_collection[2, :, 0], rod_one.radius + cylinder_two.radius, @@ -797,6 +829,8 @@ def apply_forces(self, rod_one, index_one, rod_two, index_two): rod_one.internal_forces, rod_one.external_forces, cylinder_two.external_forces, + cylinder_two.external_torques, + cylinder_two.director_collection[:, :, 0], rod_one.velocity_collection, cylinder_two.velocity_collection, self.k, From 8e7614a42e9443020b09cd98db5f9c8e6439319b Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Thu, 23 Jun 2022 23:42:15 -0500 Subject: [PATCH 15/37] enhancement:add parallel rod-cylinder contact simulation --- examples/README.md | 3 + .../RodRigidBodyContact/post_processing.py | 528 +++++++++++++++++- .../rod_cylinder_parallel_contact.py | 222 ++++++++ 3 files changed, 752 insertions(+), 1 deletion(-) create mode 100644 examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_parallel_contact.py diff --git a/examples/README.md b/examples/README.md index 52800518..713d9f83 100644 --- a/examples/README.md +++ b/examples/README.md @@ -36,6 +36,9 @@ Examples can serve as a starting template for customized usages. * [RigidBodyCases](./RigidBodyCases) * __Purpose__: Demonstrate usage of rigid body on simulation. * __Features__: Cylinder, Sphere + * [RodRigidBodyContact](./RigidBodyCases/RodRigidBodyContact) + * __Purpose__: Demonstrate contact between cylinder and rod, for different intial conditions. + * __Features__: Cylinder, CosseratRods, ExternalContact * [HelicalBucklingCase](./HelicalBucklingCase) * __Purpose__: Demonstrate helical buckling with extreme twisting boundary condition. * __Features__: HelicalBucklingBC diff --git a/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py b/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py index 0c519698..5a0684e1 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py @@ -1,5 +1,8 @@ import numpy as np from matplotlib import pyplot as plt +from matplotlib import cm +from typing import Dict, Sequence +from tqdm import tqdm def make_data_for_cylinder_along_y(cstart, cradius, cheight): @@ -21,7 +24,7 @@ def plot_video( fps=60, step=1, *args, - **kwargs + **kwargs, ): # (time step, x/y/z, node) cylinder_start = np.array(cylinder_history["position"])[0, ...] @@ -156,3 +159,526 @@ def plot_cylinder_rod_position( plt.ylabel("Position", fontsize=20) fig.legend(prop={"size": 20}) plt.show() + + +def plot_velocity( + plot_params_rod_one: dict, + plot_params_rod_two: dict, + filename="velocity.png", + SAVE_FIGURE=False, +): + time = np.array(plot_params_rod_one["time"]) + avg_velocity_rod_one = np.array(plot_params_rod_one["com_velocity"]) + avg_velocity_rod_two = np.array(plot_params_rod_two["com_velocity"]) + total_energy_rod_one = np.array(plot_params_rod_one["total_energy"]) + total_energy_rod_two = np.array(plot_params_rod_two["total_energy"]) + + fig = plt.figure(figsize=(12, 10), frameon=True, dpi=150) + axs = [] + axs.append(plt.subplot2grid((4, 1), (0, 0))) + axs.append(plt.subplot2grid((4, 1), (1, 0))) + axs.append(plt.subplot2grid((4, 1), (2, 0))) + axs.append(plt.subplot2grid((4, 1), (3, 0))) + + axs[0].plot(time[:], avg_velocity_rod_one[:, 0], linewidth=3, label="rod_one") + axs[0].plot(time[:], avg_velocity_rod_two[:, 0], linewidth=3, label="rod_two") + axs[0].plot( + time[:], + avg_velocity_rod_one[:, 0] + avg_velocity_rod_two[:, 0], + "--", + linewidth=3, + label="total", + ) + axs[0].set_ylabel("x velocity", fontsize=20) + + axs[1].plot( + time[:], + avg_velocity_rod_one[:, 1], + linewidth=3, + ) + axs[1].plot( + time[:], + avg_velocity_rod_two[:, 1], + linewidth=3, + ) + axs[1].plot( + time[:], + avg_velocity_rod_one[:, 1] + avg_velocity_rod_two[:, 1], + "--", + linewidth=3, + ) + axs[1].set_ylabel("y velocity", fontsize=20) + + axs[2].plot( + time[:], + avg_velocity_rod_one[:, 2], + linewidth=3, + ) + axs[2].plot( + time[:], + avg_velocity_rod_two[:, 2], + linewidth=3, + ) + axs[2].plot( + time[:], + avg_velocity_rod_one[:, 2] + avg_velocity_rod_two[:, 2], + "--", + linewidth=3, + ) + axs[2].set_ylabel("z velocity", fontsize=20) + + axs[3].semilogy( + time[:], + total_energy_rod_one[:], + linewidth=3, + ) + axs[3].semilogy( + time[:], + total_energy_rod_two[:], + linewidth=3, + ) + axs[3].semilogy( + time[:], + np.abs(total_energy_rod_one[:] - total_energy_rod_two[:]), + "--", + linewidth=3, + ) + axs[3].set_ylabel("total_energy", fontsize=20) + axs[3].set_xlabel("time [s]", fontsize=20) + + plt.tight_layout() + # fig.align_ylabels() + fig.legend(prop={"size": 20}) + # fig.savefig(filename) + # plt.show() + plt.close(plt.gcf()) + + if SAVE_FIGURE: + fig.savefig(filename) + + +def plot_video_with_surface( + rods_history: Sequence[Dict], + video_name="video.mp4", + fps=60, + step=1, + vis2D=True, + **kwargs, +): + plt.rcParams.update({"font.size": 22}) + + folder_name = kwargs.get("folder_name", "") + + # 2d case + import matplotlib.animation as animation + + # simulation time + sim_time = np.array(rods_history[0]["time"]) + + # Rod + n_visualized_rods = len(rods_history) # should be one for now + # Rod info + rod_history_unpacker = lambda rod_idx, t_idx: ( + rods_history[rod_idx]["position"][t_idx], + rods_history[rod_idx]["radius"][t_idx], + ) + # Rod center of mass + com_history_unpacker = lambda rod_idx, t_idx: rods_history[rod_idx]["com"][time_idx] + + # Generate target sphere data + sphere_flag = False + if kwargs.__contains__("sphere_history"): + sphere_flag = True + sphere_history = kwargs.get("sphere_history") + n_visualized_spheres = len(sphere_history) # should be one for now + sphere_history_unpacker = lambda sph_idx, t_idx: ( + sphere_history[sph_idx]["position"][t_idx], + sphere_history[sph_idx]["radius"][t_idx], + ) + # color mapping + sphere_cmap = cm.get_cmap("Spectral", n_visualized_spheres) + + # video pre-processing + print("plot scene visualization video") + FFMpegWriter = animation.writers["ffmpeg"] + metadata = dict(title="Movie Test", artist="Matplotlib", comment="Movie support!") + writer = FFMpegWriter(fps=fps, metadata=metadata) + dpi = kwargs.get("dpi", 100) + + xlim = kwargs.get("x_limits", (-1.0, 1.0)) + ylim = kwargs.get("y_limits", (-1.0, 1.0)) + zlim = kwargs.get("z_limits", (-0.05, 1.0)) + + difference = lambda x: x[1] - x[0] + max_axis_length = max(difference(xlim), difference(ylim)) + # The scaling factor from physical space to matplotlib space + scaling_factor = (2 * 0.1) / max_axis_length # Octopus head dimension + scaling_factor *= 2.6e3 # Along one-axis + + if kwargs.get("vis3D", True): + fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) + ax = plt.axes(projection="3d") + + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + ax.set_zlim(*zlim) + + time_idx = 0 + rod_lines = [None for _ in range(n_visualized_rods)] + rod_com_lines = [None for _ in range(n_visualized_rods)] + rod_scatters = [None for _ in range(n_visualized_rods)] + + 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] = ax.scatter( + inst_position[0], + inst_position[1], + inst_position[2], + s=np.pi * (scaling_factor * inst_radius) ** 2, + ) + + if sphere_flag: + sphere_artists = [None for _ in range(n_visualized_spheres)] + for sphere_idx in range(n_visualized_spheres): + sphere_position, sphere_radius = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx] = ax.scatter( + sphere_position[0], + sphere_position[1], + sphere_position[2], + s=np.pi * (scaling_factor * sphere_radius) ** 2, + ) + # sphere_radius, + # color=sphere_cmap(sphere_idx),) + ax.add_artist(sphere_artists[sphere_idx]) + + # ax.set_aspect("equal") + video_name_3D = folder_name + "3D_" + video_name + + with writer.saving(fig, video_name_3D, dpi): + with plt.style.context("seaborn-whitegrid"): + 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], + ) + + # 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], + ) + + writer.grab_frame() + + # Be a good boy and close figures + # https://stackoverflow.com/a/37451036 + # plt.close(fig) alone does not suffice + # See https://github.com/matplotlib/matplotlib/issues/8560/ + plt.close(plt.gcf()) + + if kwargs.get("vis2D", True): + max_axis_length = max(difference(xlim), difference(ylim)) + # The scaling factor from physical space to matplotlib space + scaling_factor = (2 * 0.1) / max_axis_length # Octopus head dimension + scaling_factor *= 2.6e3 # Along one-axis + + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + + time_idx = 0 + rod_lines = [None for _ in range(n_visualized_rods)] + rod_com_lines = [None for _ in range(n_visualized_rods)] + rod_scatters = [None for _ in range(n_visualized_rods)] + + 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_lines[rod_idx] = ax.plot( + inst_position[0], inst_position[1], "r", lw=0.5 + )[0] + inst_com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx] = ax.plot(inst_com[0], inst_com[1], "k--", lw=2.0)[0] + + rod_scatters[rod_idx] = ax.scatter( + inst_position[0], + inst_position[1], + s=np.pi * (scaling_factor * inst_radius) ** 2, + ) + + if sphere_flag: + sphere_artists = [None for _ in range(n_visualized_spheres)] + for sphere_idx in range(n_visualized_spheres): + sphere_position, sphere_radius = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx] = Circle( + (sphere_position[0], sphere_position[1]), + sphere_radius, + color=sphere_cmap(sphere_idx), + ) + ax.add_artist(sphere_artists[sphere_idx]) + + ax.set_aspect("equal") + video_name_2D = folder_name + "2D_xy_" + video_name + + with writer.saving(fig, video_name_2D, dpi): + with plt.style.context("seaborn-whitegrid"): + 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_lines[rod_idx].set_xdata(inst_position[0]) + rod_lines[rod_idx].set_ydata(inst_position[1]) + + com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx].set_xdata(com[0]) + rod_com_lines[rod_idx].set_ydata(com[1]) + + 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].center = ( + sphere_position[0], + sphere_position[1], + ) + + writer.grab_frame() + + # Be a good boy and close figures + # https://stackoverflow.com/a/37451036 + # plt.close(fig) alone does not suffice + # See https://github.com/matplotlib/matplotlib/issues/8560/ + plt.close(plt.gcf()) + + # Plot zy + max_axis_length = max(difference(zlim), difference(ylim)) + # The scaling factor from physical space to matplotlib space + scaling_factor = (2 * 0.1) / max_axis_length # Octopus head dimension + scaling_factor *= 2.6e3 # Along one-axis + + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + ax.set_xlim(*zlim) + ax.set_ylim(*ylim) + + time_idx = 0 + rod_lines = [None for _ in range(n_visualized_rods)] + rod_com_lines = [None for _ in range(n_visualized_rods)] + rod_scatters = [None for _ in range(n_visualized_rods)] + + 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_lines[rod_idx] = ax.plot( + inst_position[2], inst_position[1], "r", lw=0.5 + )[0] + inst_com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx] = ax.plot(inst_com[2], inst_com[1], "k--", lw=2.0)[0] + + rod_scatters[rod_idx] = ax.scatter( + inst_position[2], + inst_position[1], + s=np.pi * (scaling_factor * inst_radius) ** 2, + ) + + if sphere_flag: + sphere_artists = [None for _ in range(n_visualized_spheres)] + for sphere_idx in range(n_visualized_spheres): + sphere_position, sphere_radius = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx] = Circle( + (sphere_position[2], sphere_position[1]), + sphere_radius, + color=sphere_cmap(sphere_idx), + ) + ax.add_artist(sphere_artists[sphere_idx]) + + ax.set_aspect("equal") + video_name_2D = folder_name + "2D_zy_" + video_name + + with writer.saving(fig, video_name_2D, dpi): + with plt.style.context("seaborn-whitegrid"): + 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_lines[rod_idx].set_xdata(inst_position[2]) + rod_lines[rod_idx].set_ydata(inst_position[1]) + + com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx].set_xdata(com[2]) + rod_com_lines[rod_idx].set_ydata(com[1]) + + rod_scatters[rod_idx].set_offsets( + np.vstack((inst_position[2], inst_position[1])).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].center = ( + sphere_position[2], + sphere_position[1], + ) + + writer.grab_frame() + + # Be a good boy and close figures + # https://stackoverflow.com/a/37451036 + # plt.close(fig) alone does not suffice + # See https://github.com/matplotlib/matplotlib/issues/8560/ + plt.close(plt.gcf()) + + # Plot xz + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + ax.set_xlim(*xlim) + ax.set_ylim(*zlim) + + # The scaling factor from physical space to matplotlib space + max_axis_length = max(difference(zlim), difference(xlim)) + scaling_factor = (2 * 0.1) / (max_axis_length) # Octopus head dimension + scaling_factor *= 2.6e3 # Along one-axis + + time_idx = 0 + rod_lines = [None for _ in range(n_visualized_rods)] + rod_com_lines = [None for _ in range(n_visualized_rods)] + rod_scatters = [None for _ in range(n_visualized_rods)] + + 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_lines[rod_idx] = ax.plot( + inst_position[0], inst_position[2], "r", lw=0.5 + )[0] + inst_com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx] = ax.plot(inst_com[0], inst_com[2], "k--", lw=2.0)[0] + + rod_scatters[rod_idx] = ax.scatter( + inst_position[0], + inst_position[2], + s=np.pi * (scaling_factor * inst_radius) ** 2, + ) + + if sphere_flag: + sphere_artists = [None for _ in range(n_visualized_spheres)] + for sphere_idx in range(n_visualized_spheres): + sphere_position, sphere_radius = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx] = Circle( + (sphere_position[0], sphere_position[2]), + sphere_radius, + color=sphere_cmap(sphere_idx), + ) + ax.add_artist(sphere_artists[sphere_idx]) + + ax.set_aspect("equal") + video_name_2D = folder_name + "2D_xz_" + video_name + + with writer.saving(fig, video_name_2D, dpi): + with plt.style.context("seaborn-whitegrid"): + 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_lines[rod_idx].set_xdata(inst_position[0]) + rod_lines[rod_idx].set_ydata(inst_position[2]) + + com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx].set_xdata(com[0]) + rod_com_lines[rod_idx].set_ydata(com[2]) + + rod_scatters[rod_idx].set_offsets( + np.vstack((inst_position[0], 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].center = ( + sphere_position[0], + sphere_position[2], + ) + + writer.grab_frame() + + # Be a good boy and close figures + # https://stackoverflow.com/a/37451036 + # plt.close(fig) alone does not suffice + # See https://github.com/matplotlib/matplotlib/issues/8560/ + plt.close(plt.gcf()) diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_parallel_contact.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_parallel_contact.py new file mode 100644 index 00000000..7cd63776 --- /dev/null +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_parallel_contact.py @@ -0,0 +1,222 @@ +import numpy as np +import sys + +sys.path.append("../../../") +from elastica import * +from post_processing import plot_velocity, plot_video_with_surface + + +class RodCylinderParallelContact( + BaseSystemCollection, Constraints, Connections, CallBacks, Forcing +): + pass + + +rod_cylinder_parallel_contact_simulator = RodCylinderParallelContact() + +# time step etc +final_time = 5.0 +time_step = 1e-4 +total_steps = int(final_time / time_step) + 1 +rendering_fps = 30 # 20 * 1e1 +step_skip = int(1.0 / (rendering_fps * time_step)) + +base_length = 0.5 +base_radius = 0.1 +density = 1750 +E = 3e5 +poisson_ratio = 0.5 +shear_modulus = E / (2 * (1 + poisson_ratio)) +n_elem = 50 +nu = 0.0 +start = np.zeros((3,)) +direction = np.array([0.0, 0.0, 1.0]) +normal = np.array([0.0, 1.0, 0.0]) +binormal = np.cross(direction, normal) + +rod = CosseratRod.straight_rod( + n_elem, + start, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + shear_modulus=shear_modulus, +) + +rod_cylinder_parallel_contact_simulator.append(rod) +rod.velocity_collection[0, :] -= 0.2 + + +cylinder_height = base_length +cylinder_radius = base_radius + +cylinder_start = start + binormal +cylinder_direction = direction +cylinder_normal = normal + +rigid_body = Cylinder( + start=cylinder_start, + direction=cylinder_direction, + normal=cylinder_normal, + base_length=cylinder_height, + base_radius=cylinder_radius, + density=density, +) +rod_cylinder_parallel_contact_simulator.append(rigid_body) + +# Add contact between rigid body and rod +rod_cylinder_parallel_contact_simulator.connect(rod, rigid_body).using( + ExternalContact, k=1e3, nu=0.1 +) + +# Add callbacks +post_processing_dict_list = [] +# For rod +class StraightRodCallBack(CallBackBaseClass): + """ + Call back function for two arm octopus + """ + + 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["position"].append(system.position_collection.copy()) + self.callback_params["radius"].append(system.radius.copy()) + self.callback_params["com"].append(system.compute_position_center_of_mass()) + if current_step == 0: + self.callback_params["lengths"].append(system.rest_lengths.copy()) + else: + self.callback_params["lengths"].append(system.lengths.copy()) + + self.callback_params["com_velocity"].append( + system.compute_velocity_center_of_mass() + ) + + total_energy = ( + system.compute_translational_energy() + + system.compute_rotational_energy() + + system.compute_bending_energy() + + system.compute_shear_energy() + ) + self.callback_params["total_energy"].append(total_energy) + + return + + +class RigidCylinderCallBack(CallBackBaseClass): + """ + Call back function for two arm octopus + """ + + def __init__( + self, step_skip: int, callback_params: dict, resize_cylinder_elems: int + ): + CallBackBaseClass.__init__(self) + self.every = step_skip + self.callback_params = callback_params + self.n_elem_cylinder = resize_cylinder_elems + self.n_node_cylinder = self.n_elem_cylinder + 1 + + 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) + + cylinder_center_position = system.position_collection + cylinder_length = system.length + cylinder_direction = system.director_collection[2, :, :].reshape(3, 1) + cylinder_radius = system.radius + + # Expand cylinder data. Create multiple points on cylinder later to use for rendering. + + start_position = ( + cylinder_center_position - cylinder_length / 2 * cylinder_direction + ) + + cylinder_position_collection = ( + start_position + + np.linspace(0, cylinder_length[0], self.n_node_cylinder) + * cylinder_direction + ) + cylinder_radius_collection = ( + np.ones((self.n_elem_cylinder)) * cylinder_radius + ) + cylinder_length_collection = ( + np.ones((self.n_elem_cylinder)) * cylinder_length + ) + cylinder_velocity_collection = ( + np.ones((self.n_node_cylinder)) * system.velocity_collection + ) + + self.callback_params["position"].append(cylinder_position_collection.copy()) + self.callback_params["velocity"].append(cylinder_velocity_collection.copy()) + self.callback_params["radius"].append(cylinder_radius_collection.copy()) + self.callback_params["com"].append(system.compute_position_center_of_mass()) + + self.callback_params["lengths"].append(cylinder_length_collection.copy()) + self.callback_params["com_velocity"].append( + system.velocity_collection[..., 0].copy() + ) + + total_energy = ( + system.compute_translational_energy() + + system.compute_rotational_energy() + ) + self.callback_params["total_energy"].append(total_energy[..., 0].copy()) + + return + + +post_processing_dict_list.append(defaultdict(list)) +rod_cylinder_parallel_contact_simulator.collect_diagnostics(rod).using( + StraightRodCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_list[0], +) +# For rigid body +post_processing_dict_list.append(defaultdict(list)) +rod_cylinder_parallel_contact_simulator.collect_diagnostics(rigid_body).using( + RigidCylinderCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_list[1], + resize_cylinder_elems=n_elem, +) + +rod_cylinder_parallel_contact_simulator.finalize() +timestepper = PositionVerlet() + +integrate(timestepper, rod_cylinder_parallel_contact_simulator, final_time, total_steps) + + +# Plot the rods +plot_video_with_surface( + post_processing_dict_list, + video_name="rod_cylinder_contact.mp4", + fps=rendering_fps, + step=1, + # The following parameters are optional + x_limits=(-base_length * 5, base_length * 5), # Set bounds on x-axis + y_limits=(-base_length * 5, base_length * 5), # Set bounds on y-axis + z_limits=(-base_length * 5, base_length * 5), # Set bounds on z-axis + dpi=100, # Set the quality of the image + vis3D=True, # Turn on 3D visualization + vis2D=True, # Turn on projected (2D) visualization +) + +filaname = "rod_rigid_velocity.png" +plot_velocity( + post_processing_dict_list[0], + post_processing_dict_list[1], + filename=filaname, + SAVE_FIGURE=True, +) From 96f74fd50cd08d076ea5c363932ffa314948cc32 Mon Sep 17 00:00:00 2001 From: Bhosale Date: Thu, 30 Jun 2022 10:57:10 -0500 Subject: [PATCH 16/37] fix: plotting rod radius --- .../RodRigidBodyContact/post_processing.py | 61 ++++++++++++++++--- ...lel_contact.py => rod_cylinder_contact.py} | 0 2 files changed, 52 insertions(+), 9 deletions(-) rename examples/RigidBodyCases/RodRigidBodyContact/{rod_cylinder_parallel_contact.py => rod_cylinder_contact.py} (100%) diff --git a/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py b/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py index 5a0684e1..14b81311 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py @@ -337,11 +337,18 @@ def plot_video_with_surface( if not inst_position.shape[1] == inst_radius.shape[0]: inst_position = 0.5 * (inst_position[..., 1:] + inst_position[..., :-1]) + # for reference see + # https://stackoverflow.com/questions/48172928/scale-matplotlib-pyplot + # -axes-scatter-markersize-by-x-scale/48174228#48174228 + scaling_factor = ( + ax.get_window_extent().width / (max_axis_length) * 72.0 / fig.dpi + ) rod_scatters[rod_idx] = ax.scatter( inst_position[0], inst_position[1], inst_position[2], - s=np.pi * (scaling_factor * inst_radius) ** 2, + # for circle s = 4/pi*area = 4 * r^2 + s=4 * (scaling_factor * inst_radius) ** 2, ) if sphere_flag: @@ -350,11 +357,14 @@ def plot_video_with_surface( sphere_position, sphere_radius = sphere_history_unpacker( sphere_idx, time_idx ) + scaling_factor = ( + ax.get_window_extent().width / (max_axis_length) * 72.0 / fig.dpi + ) sphere_artists[sphere_idx] = ax.scatter( sphere_position[0], sphere_position[1], sphere_position[2], - s=np.pi * (scaling_factor * sphere_radius) ** 2, + s=4 * (scaling_factor * sphere_radius) ** 2, ) # sphere_radius, # color=sphere_cmap(sphere_idx),) @@ -382,9 +392,15 @@ def plot_video_with_surface( inst_position[2], ) + scaling_factor = ( + ax.get_window_extent().width + / (max_axis_length) + * 72.0 + / fig.dpi + ) # rod_scatters[rod_idx].set_offsets(inst_position[:2].T) rod_scatters[rod_idx].set_sizes( - np.pi * (scaling_factor * inst_radius) ** 2 + 4 * (scaling_factor * inst_radius) ** 2 ) if sphere_flag: @@ -432,10 +448,13 @@ def plot_video_with_surface( inst_com = com_history_unpacker(rod_idx, time_idx) rod_com_lines[rod_idx] = ax.plot(inst_com[0], inst_com[1], "k--", lw=2.0)[0] + scaling_factor = ( + ax.get_window_extent().width / (max_axis_length) * 72.0 / fig.dpi + ) rod_scatters[rod_idx] = ax.scatter( inst_position[0], inst_position[1], - s=np.pi * (scaling_factor * inst_radius) ** 2, + s=4 * (scaling_factor * inst_radius) ** 2, ) if sphere_flag: @@ -475,8 +494,14 @@ def plot_video_with_surface( rod_com_lines[rod_idx].set_ydata(com[1]) rod_scatters[rod_idx].set_offsets(inst_position[:2].T) + scaling_factor = ( + ax.get_window_extent().width + / (max_axis_length) + * 72.0 + / fig.dpi + ) rod_scatters[rod_idx].set_sizes( - np.pi * (scaling_factor * inst_radius) ** 2 + 4 * (scaling_factor * inst_radius) ** 2 ) if sphere_flag: @@ -523,10 +548,13 @@ def plot_video_with_surface( inst_com = com_history_unpacker(rod_idx, time_idx) rod_com_lines[rod_idx] = ax.plot(inst_com[2], inst_com[1], "k--", lw=2.0)[0] + scaling_factor = ( + ax.get_window_extent().width / (max_axis_length) * 72.0 / fig.dpi + ) rod_scatters[rod_idx] = ax.scatter( inst_position[2], inst_position[1], - s=np.pi * (scaling_factor * inst_radius) ** 2, + s=4 * (scaling_factor * inst_radius) ** 2, ) if sphere_flag: @@ -568,8 +596,14 @@ def plot_video_with_surface( rod_scatters[rod_idx].set_offsets( np.vstack((inst_position[2], inst_position[1])).T ) + scaling_factor = ( + ax.get_window_extent().width + / (max_axis_length) + * 72.0 + / fig.dpi + ) rod_scatters[rod_idx].set_sizes( - np.pi * (scaling_factor * inst_radius) ** 2 + 4 * (scaling_factor * inst_radius) ** 2 ) if sphere_flag: @@ -616,10 +650,13 @@ def plot_video_with_surface( inst_com = com_history_unpacker(rod_idx, time_idx) rod_com_lines[rod_idx] = ax.plot(inst_com[0], inst_com[2], "k--", lw=2.0)[0] + scaling_factor = ( + ax.get_window_extent().width / (max_axis_length) * 72.0 / fig.dpi + ) rod_scatters[rod_idx] = ax.scatter( inst_position[0], inst_position[2], - s=np.pi * (scaling_factor * inst_radius) ** 2, + s=4 * (scaling_factor * inst_radius) ** 2, ) if sphere_flag: @@ -661,8 +698,14 @@ def plot_video_with_surface( rod_scatters[rod_idx].set_offsets( np.vstack((inst_position[0], inst_position[2])).T ) + scaling_factor = ( + ax.get_window_extent().width + / (max_axis_length) + * 72.0 + / fig.dpi + ) rod_scatters[rod_idx].set_sizes( - np.pi * (scaling_factor * inst_radius) ** 2 + 4 * (scaling_factor * inst_radius) ** 2 ) if sphere_flag: diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_parallel_contact.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py similarity index 100% rename from examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_parallel_contact.py rename to examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py From a241f650785f9dd97896e9c4aeda72353fc597d4 Mon Sep 17 00:00:00 2001 From: Bhosale Date: Thu, 30 Jun 2022 10:57:38 -0500 Subject: [PATCH 17/37] feat: generic inclination angle collision functionality --- .../rod_cylinder_contact.py | 435 +++++++++--------- .../rod_cylinder_inclined_contact.py | 7 + .../rod_cylinder_parallel_contact.py | 6 + 3 files changed, 235 insertions(+), 213 deletions(-) create mode 100644 examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_inclined_contact.py create mode 100644 examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_parallel_contact.py diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py index 7cd63776..c362f566 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py @@ -6,217 +6,226 @@ from post_processing import plot_velocity, plot_video_with_surface -class RodCylinderParallelContact( - BaseSystemCollection, Constraints, Connections, CallBacks, Forcing -): - pass - - -rod_cylinder_parallel_contact_simulator = RodCylinderParallelContact() - -# time step etc -final_time = 5.0 -time_step = 1e-4 -total_steps = int(final_time / time_step) + 1 -rendering_fps = 30 # 20 * 1e1 -step_skip = int(1.0 / (rendering_fps * time_step)) - -base_length = 0.5 -base_radius = 0.1 -density = 1750 -E = 3e5 -poisson_ratio = 0.5 -shear_modulus = E / (2 * (1 + poisson_ratio)) -n_elem = 50 -nu = 0.0 -start = np.zeros((3,)) -direction = np.array([0.0, 0.0, 1.0]) -normal = np.array([0.0, 1.0, 0.0]) -binormal = np.cross(direction, normal) - -rod = CosseratRod.straight_rod( - n_elem, - start, - direction, - normal, - base_length, - base_radius, - density, - nu, - E, - shear_modulus=shear_modulus, -) - -rod_cylinder_parallel_contact_simulator.append(rod) -rod.velocity_collection[0, :] -= 0.2 - - -cylinder_height = base_length -cylinder_radius = base_radius - -cylinder_start = start + binormal -cylinder_direction = direction -cylinder_normal = normal - -rigid_body = Cylinder( - start=cylinder_start, - direction=cylinder_direction, - normal=cylinder_normal, - base_length=cylinder_height, - base_radius=cylinder_radius, - density=density, -) -rod_cylinder_parallel_contact_simulator.append(rigid_body) - -# Add contact between rigid body and rod -rod_cylinder_parallel_contact_simulator.connect(rod, rigid_body).using( - ExternalContact, k=1e3, nu=0.1 -) - -# Add callbacks -post_processing_dict_list = [] -# For rod -class StraightRodCallBack(CallBackBaseClass): - """ - Call back function for two arm octopus - """ - - 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["position"].append(system.position_collection.copy()) - self.callback_params["radius"].append(system.radius.copy()) - self.callback_params["com"].append(system.compute_position_center_of_mass()) - if current_step == 0: - self.callback_params["lengths"].append(system.rest_lengths.copy()) - else: - self.callback_params["lengths"].append(system.lengths.copy()) - - self.callback_params["com_velocity"].append( - system.compute_velocity_center_of_mass() - ) - - total_energy = ( - system.compute_translational_energy() - + system.compute_rotational_energy() - + system.compute_bending_energy() - + system.compute_shear_energy() - ) - self.callback_params["total_energy"].append(total_energy) - - return - - -class RigidCylinderCallBack(CallBackBaseClass): - """ - Call back function for two arm octopus - """ - - def __init__( - self, step_skip: int, callback_params: dict, resize_cylinder_elems: int +def rod_cylinder_contact_case(inclination_angle=0.0): + class RodCylinderParallelContact( + BaseSystemCollection, Constraints, Connections, CallBacks, Forcing ): - CallBackBaseClass.__init__(self) - self.every = step_skip - self.callback_params = callback_params - self.n_elem_cylinder = resize_cylinder_elems - self.n_node_cylinder = self.n_elem_cylinder + 1 - - 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) - - cylinder_center_position = system.position_collection - cylinder_length = system.length - cylinder_direction = system.director_collection[2, :, :].reshape(3, 1) - cylinder_radius = system.radius - - # Expand cylinder data. Create multiple points on cylinder later to use for rendering. - - start_position = ( - cylinder_center_position - cylinder_length / 2 * cylinder_direction - ) - - cylinder_position_collection = ( - start_position - + np.linspace(0, cylinder_length[0], self.n_node_cylinder) - * cylinder_direction - ) - cylinder_radius_collection = ( - np.ones((self.n_elem_cylinder)) * cylinder_radius - ) - cylinder_length_collection = ( - np.ones((self.n_elem_cylinder)) * cylinder_length - ) - cylinder_velocity_collection = ( - np.ones((self.n_node_cylinder)) * system.velocity_collection - ) - - self.callback_params["position"].append(cylinder_position_collection.copy()) - self.callback_params["velocity"].append(cylinder_velocity_collection.copy()) - self.callback_params["radius"].append(cylinder_radius_collection.copy()) - self.callback_params["com"].append(system.compute_position_center_of_mass()) - - self.callback_params["lengths"].append(cylinder_length_collection.copy()) - self.callback_params["com_velocity"].append( - system.velocity_collection[..., 0].copy() - ) - - total_energy = ( - system.compute_translational_energy() - + system.compute_rotational_energy() - ) - self.callback_params["total_energy"].append(total_energy[..., 0].copy()) - - return - - -post_processing_dict_list.append(defaultdict(list)) -rod_cylinder_parallel_contact_simulator.collect_diagnostics(rod).using( - StraightRodCallBack, - step_skip=step_skip, - callback_params=post_processing_dict_list[0], -) -# For rigid body -post_processing_dict_list.append(defaultdict(list)) -rod_cylinder_parallel_contact_simulator.collect_diagnostics(rigid_body).using( - RigidCylinderCallBack, - step_skip=step_skip, - callback_params=post_processing_dict_list[1], - resize_cylinder_elems=n_elem, -) - -rod_cylinder_parallel_contact_simulator.finalize() -timestepper = PositionVerlet() - -integrate(timestepper, rod_cylinder_parallel_contact_simulator, final_time, total_steps) - - -# Plot the rods -plot_video_with_surface( - post_processing_dict_list, - video_name="rod_cylinder_contact.mp4", - fps=rendering_fps, - step=1, - # The following parameters are optional - x_limits=(-base_length * 5, base_length * 5), # Set bounds on x-axis - y_limits=(-base_length * 5, base_length * 5), # Set bounds on y-axis - z_limits=(-base_length * 5, base_length * 5), # Set bounds on z-axis - dpi=100, # Set the quality of the image - vis3D=True, # Turn on 3D visualization - vis2D=True, # Turn on projected (2D) visualization -) - -filaname = "rod_rigid_velocity.png" -plot_velocity( - post_processing_dict_list[0], - post_processing_dict_list[1], - filename=filaname, - SAVE_FIGURE=True, -) + pass + + rod_cylinder_parallel_contact_simulator = RodCylinderParallelContact() + + # time step etc + final_time = 5.0 + time_step = 1e-4 + total_steps = int(final_time / time_step) + 1 + rendering_fps = 30 # 20 * 1e1 + step_skip = int(1.0 / (rendering_fps * time_step)) + + base_length = 0.5 + base_radius = 0.1 + density = 1750 + E = 3e5 + poisson_ratio = 0.5 + shear_modulus = E / (2 * (1 + poisson_ratio)) + n_elem = 50 + nu = 0.0 + start = np.zeros((3,)) + direction = np.array([np.sin(inclination_angle), 0.0, np.cos(inclination_angle)]) + normal = np.array([0.0, 1.0, 0.0]) + + rod = CosseratRod.straight_rod( + n_elem, + start, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + shear_modulus=shear_modulus, + ) + + rod_cylinder_parallel_contact_simulator.append(rod) + rod.velocity_collection[0, :] -= 0.2 + + cylinder_height = base_length + cylinder_radius = base_radius + + cylinder_start = start + np.array([-1.0, 0.0, 0.0]) + cylinder_direction = np.array([0.0, 0.0, 1.0]) + cylinder_normal = np.array([0.0, 1.0, 0.0]) + + rigid_body = Cylinder( + start=cylinder_start, + direction=cylinder_direction, + normal=cylinder_normal, + base_length=cylinder_height, + base_radius=cylinder_radius, + density=density, + ) + rod_cylinder_parallel_contact_simulator.append(rigid_body) + + # Add contact between rigid body and rod + rod_cylinder_parallel_contact_simulator.connect(rod, rigid_body).using( + ExternalContact, k=1e3, nu=0.1 + ) + + # Add callbacks + post_processing_dict_list = [] + # For rod + class StraightRodCallBack(CallBackBaseClass): + """ + Call back function for two arm octopus + """ + + 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["position"].append( + system.position_collection.copy() + ) + self.callback_params["radius"].append(system.radius.copy()) + self.callback_params["com"].append( + system.compute_position_center_of_mass() + ) + if current_step == 0: + self.callback_params["lengths"].append(system.rest_lengths.copy()) + else: + self.callback_params["lengths"].append(system.lengths.copy()) + + self.callback_params["com_velocity"].append( + system.compute_velocity_center_of_mass() + ) + + total_energy = ( + system.compute_translational_energy() + + system.compute_rotational_energy() + + system.compute_bending_energy() + + system.compute_shear_energy() + ) + self.callback_params["total_energy"].append(total_energy) + + return + + class RigidCylinderCallBack(CallBackBaseClass): + """ + Call back function for two arm octopus + """ + + def __init__( + self, step_skip: int, callback_params: dict, resize_cylinder_elems: int + ): + CallBackBaseClass.__init__(self) + self.every = step_skip + self.callback_params = callback_params + self.n_elem_cylinder = resize_cylinder_elems + self.n_node_cylinder = self.n_elem_cylinder + 1 + + 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) + + cylinder_center_position = system.position_collection + cylinder_length = system.length + cylinder_direction = system.director_collection[2, :, :].reshape(3, 1) + cylinder_radius = system.radius + + # Expand cylinder data. Create multiple points on cylinder later to use for rendering. + + start_position = ( + cylinder_center_position - cylinder_length / 2 * cylinder_direction + ) + + cylinder_position_collection = ( + start_position + + np.linspace(0, cylinder_length[0], self.n_node_cylinder) + * cylinder_direction + ) + cylinder_radius_collection = ( + np.ones((self.n_elem_cylinder)) * cylinder_radius + ) + cylinder_length_collection = ( + np.ones((self.n_elem_cylinder)) * cylinder_length + ) + cylinder_velocity_collection = ( + np.ones((self.n_node_cylinder)) * system.velocity_collection + ) + + self.callback_params["position"].append( + cylinder_position_collection.copy() + ) + self.callback_params["velocity"].append( + cylinder_velocity_collection.copy() + ) + self.callback_params["radius"].append(cylinder_radius_collection.copy()) + self.callback_params["com"].append( + system.compute_position_center_of_mass() + ) + + self.callback_params["lengths"].append( + cylinder_length_collection.copy() + ) + self.callback_params["com_velocity"].append( + system.velocity_collection[..., 0].copy() + ) + + total_energy = ( + system.compute_translational_energy() + + system.compute_rotational_energy() + ) + self.callback_params["total_energy"].append(total_energy[..., 0].copy()) + + return + + post_processing_dict_list.append(defaultdict(list)) + rod_cylinder_parallel_contact_simulator.collect_diagnostics(rod).using( + StraightRodCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_list[0], + ) + # For rigid body + post_processing_dict_list.append(defaultdict(list)) + rod_cylinder_parallel_contact_simulator.collect_diagnostics(rigid_body).using( + RigidCylinderCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_list[1], + resize_cylinder_elems=n_elem, + ) + + rod_cylinder_parallel_contact_simulator.finalize() + timestepper = PositionVerlet() + + integrate( + timestepper, rod_cylinder_parallel_contact_simulator, final_time, total_steps + ) + + # Plot the rods + plot_video_with_surface( + post_processing_dict_list, + video_name="rod_cylinder_contact.mp4", + fps=rendering_fps, + step=1, + # The following parameters are optional + x_limits=(-base_length * 5, base_length * 5), # Set bounds on x-axis + y_limits=(-base_length * 5, base_length * 5), # Set bounds on y-axis + z_limits=(-base_length * 5, base_length * 5), # Set bounds on z-axis + dpi=100, # Set the quality of the image + vis3D=True, # Turn on 3D visualization + vis2D=True, # Turn on projected (2D) visualization + ) + + filaname = "rod_rigid_velocity.png" + plot_velocity( + post_processing_dict_list[0], + post_processing_dict_list[1], + filename=filaname, + SAVE_FIGURE=True, + ) diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_inclined_contact.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_inclined_contact.py new file mode 100644 index 00000000..3c50e6a2 --- /dev/null +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_inclined_contact.py @@ -0,0 +1,7 @@ +from rod_cylinder_contact import rod_cylinder_contact_case + +import numpy as np + +if __name__ == "__main__": + inclination_angle = np.deg2rad(-30.0) + rod_cylinder_contact_case(inclination_angle=inclination_angle) diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_parallel_contact.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_parallel_contact.py new file mode 100644 index 00000000..127cff45 --- /dev/null +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_parallel_contact.py @@ -0,0 +1,6 @@ +from rod_cylinder_contact import rod_cylinder_contact_case + + +if __name__ == "__main__": + inclination_angle = 0.0 + rod_cylinder_contact_case(inclination_angle=inclination_angle) From a1c7d5360c6ed4ebd4faadcfa992f90beedeabef Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Fri, 1 Jul 2022 16:18:48 -0500 Subject: [PATCH 18/37] fix:remove normal force from rigid body rod contact --- elastica/joint.py | 16 +++++++--------- .../RodRigidBodyContact/rod_cylinder_contact.py | 8 +++++--- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/elastica/joint.py b/elastica/joint.py index 98b37d0d..7b589a12 100644 --- a/elastica/joint.py +++ b/elastica/joint.py @@ -423,22 +423,20 @@ def _calculate_contact_forces_rod_rigid_body( # As a quick fix, use this instead mask = (gamma > 0.0) * 1.0 - contact_force = contact_k * gamma + # Compute contact spring force + contact_force = contact_k * gamma * distance_vector interpenetration_velocity = ( 0.5 * (velocity_rod[..., i] + velocity_rod[..., i + 1]) - velocity_cylinder[..., 0] ) - contact_damping_force = contact_nu * _dot_product( - interpenetration_velocity, distance_vector + # Compute contact damping + normal_interpenetration_velocity = ( + _dot_product(interpenetration_velocity, distance_vector) * distance_vector ) + contact_damping_force = contact_nu * normal_interpenetration_velocity # magnitude* direction - net_contact_force = ( - 0.5 - * mask - * (normal_force + (contact_damping_force + contact_force)) - * distance_vector - ) + net_contact_force = 0.5 * mask * (contact_damping_force + contact_force) # Torques acting on the cylinder moment_arm = x_cylinder_contact_point - x_cylinder_center diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py index c362f566..97c94a8d 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact.py @@ -15,7 +15,7 @@ class RodCylinderParallelContact( rod_cylinder_parallel_contact_simulator = RodCylinderParallelContact() # time step etc - final_time = 5.0 + final_time = 10.0 time_step = 1e-4 total_steps = int(final_time / time_step) + 1 rendering_fps = 30 # 20 * 1e1 @@ -28,7 +28,7 @@ class RodCylinderParallelContact( poisson_ratio = 0.5 shear_modulus = E / (2 * (1 + poisson_ratio)) n_elem = 50 - nu = 0.0 + nu = 0.5 start = np.zeros((3,)) direction = np.array([np.sin(inclination_angle), 0.0, np.cos(inclination_angle)]) normal = np.array([0.0, 1.0, 0.0]) @@ -68,7 +68,9 @@ class RodCylinderParallelContact( # Add contact between rigid body and rod rod_cylinder_parallel_contact_simulator.connect(rod, rigid_body).using( - ExternalContact, k=1e3, nu=0.1 + ExternalContact, + k=5e4, + nu=0.1, ) # Add callbacks From a69def4056d8cdd29c3098d2fe9c11241878c4e1 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 3 Jul 2022 23:08:12 -0500 Subject: [PATCH 19/37] enhancement:implement cylinder-rod contact friction model --- elastica/joint.py | 100 +++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 90 insertions(+), 10 deletions(-) diff --git a/elastica/joint.py b/elastica/joint.py index 7b589a12..7390d6b1 100644 --- a/elastica/joint.py +++ b/elastica/joint.py @@ -377,6 +377,8 @@ def _calculate_contact_forces_rod_rigid_body( velocity_cylinder, contact_k, contact_nu, + velocity_damping_coefficient, + kinetic_friction_coefficient, ): # We already pass in only the first n_elem x n_points = x_collection_rod.shape[1] @@ -425,19 +427,49 @@ def _calculate_contact_forces_rod_rigid_body( # Compute contact spring force contact_force = contact_k * gamma * distance_vector - interpenetration_velocity = ( - 0.5 * (velocity_rod[..., i] + velocity_rod[..., i + 1]) - - velocity_cylinder[..., 0] + interpenetration_velocity = velocity_cylinder[..., 0] - 0.5 * ( + velocity_rod[..., i] + velocity_rod[..., i + 1] ) # Compute contact damping normal_interpenetration_velocity = ( _dot_product(interpenetration_velocity, distance_vector) * distance_vector ) - contact_damping_force = contact_nu * normal_interpenetration_velocity + contact_damping_force = -contact_nu * normal_interpenetration_velocity # magnitude* direction net_contact_force = 0.5 * mask * (contact_damping_force + contact_force) + # Compute friction + slip_interpenetration_velocity = ( + interpenetration_velocity - normal_interpenetration_velocity + ) + slip_interpenetration_velocity_mag = np.linalg.norm( + slip_interpenetration_velocity + ) + if slip_interpenetration_velocity_mag >= 1e-12: + slip_interpenetration_velocity_unitized = ( + slip_interpenetration_velocity / slip_interpenetration_velocity_mag + ) + else: + slip_interpenetration_velocity_unitized = np.zeros( + 3, + ) + # Compute friction force in the slip direction. + damping_force_in_slip_direction = ( + velocity_damping_coefficient * slip_interpenetration_velocity_mag + ) + # Compute kinetic friction + kinetic_friction_force = kinetic_friction_coefficient * np.linalg.norm( + net_contact_force + ) + # Compare damping force in slip direction and kinetic friction and minimum is the friction force. + friction_force = ( + -min(damping_force_in_slip_direction, kinetic_friction_force) + * slip_interpenetration_velocity_unitized + ) + # Update contact force + net_contact_force += friction_force + # Torques acting on the cylinder moment_arm = x_cylinder_contact_point - x_cylinder_center @@ -771,20 +803,66 @@ def _prune_using_aabbs_rod_rod( class ExternalContact(FreeJoint): """ - Assumes that the second entity is a rigid body for now, can be - changed at a later time + This class is for applying contact forces between rod-cylinder and rod-rod. + If you are want to apply contact forces between rod and cylinder, first system is always rod and second system + is always cylinder. + In addition to the contact forces, user can define apply friction forces between rod and cylinder that + are in contact. For details on friction model refer to this `paper `_. - Most of the cylinder-cylinder contact SHOULD be implemented - as given in this `paper `_. + TODO: Currently friction force is between rod-cylinder, in future implement friction forces between rod-rod. + + Examples + -------- + How to define contact between rod and cylinder. + + >>> simulator.connect(rod, cylidner).using( + ... ExternalContact, + ... k=1e4, + ... nu=10, + ... velocity_damping_coefficient=10, + ... kinetic_friction_coefficient=10, + ... ) - but, it isn't (the elastica-cpp kernels are implented)! + How to define contact between rod and rod. + + >>> simulator.connect(rod, rod).using( + ... ExternalContact, + ... k=1e4, + ... nu=10, + ... ) + + + Developer Note + -------------- + Most of the cylinder-cylinder contact SHOULD be implemented + as given in this `paper `, + but the elastica-cpp kernels are implemented. This is maybe to speed-up the kernel, but it's potentially dangerous as it does not deal with "end" conditions correctly. + """ - def __init__(self, k, nu): + def __init__( + self, k, nu, velocity_damping_coefficient=0, kinetic_friction_coefficient=0 + ): + """ + + Parameters + ---------- + k : float + Contact spring constant. + nu : float + Contact damping constant or coulomb coefficient of friction. + velocity_damping_coefficient : float + Velocity damping coefficient between rigid-body and rod contact is used to apply friction force in the + slip direction. + kinetic_friction_coefficient : float + Kinetic friction coefficient for rigid-body and rod contact. + """ super().__init__(k, nu) + self.velocity_damping_coefficient = velocity_damping_coefficient + self.kinetic_friction_coefficient = kinetic_friction_coefficient def apply_forces(self, rod_one, index_one, rod_two, index_two): # del index_one, index_two @@ -833,6 +911,8 @@ def apply_forces(self, rod_one, index_one, rod_two, index_two): cylinder_two.velocity_collection, self.k, self.nu, + self.velocity_damping_coefficient, + self.kinetic_friction_coefficient, ) else: From f53761c448231ff896540b91f01d922ebca38658 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 3 Jul 2022 23:09:45 -0500 Subject: [PATCH 20/37] enhancement:add rod-cylinder contact friction example case. --- .../RodRigidBodyContact/post_processing.py | 30 ++ .../rod_cylinder_contact_friction.py | 279 ++++++++++++++++++ 2 files changed, 309 insertions(+) create mode 100644 examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction.py diff --git a/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py b/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py index 14b81311..ebe9d83f 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py @@ -725,3 +725,33 @@ def plot_video_with_surface( # plt.close(fig) alone does not suffice # See https://github.com/matplotlib/matplotlib/issues/8560/ plt.close(plt.gcf()) + + +def plot_force_vs_energy( + force, + total_final_energy, + filename="energy_vs_force.png", + SAVE_FIGURE=False, +): + + fig = plt.figure(figsize=(12, 10), frameon=True, dpi=150) + axs = [] + axs.append(plt.subplot2grid((1, 1), (0, 0))) + + axs[0].plot( + force, + total_final_energy, + linewidth=3, + ) + axs[0].set_ylabel("total energy", fontsize=20) + axs[0].set_xlabel("force ", fontsize=20) + + plt.tight_layout() + # fig.align_ylabels() + fig.legend(prop={"size": 20}) + # fig.savefig(filename) + # plt.show() + plt.close(plt.gcf()) + + if SAVE_FIGURE: + fig.savefig(filename) diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction.py new file mode 100644 index 00000000..b314eccd --- /dev/null +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction.py @@ -0,0 +1,279 @@ +import numpy as np +import sys + +sys.path.append("../../../") +from elastica import * +from post_processing import plot_velocity, plot_video_with_surface + + +def rod_cylinder_contact_friction_case(pulling_force=0.0, POST_PROCESSING=False): + class RodCylinderParallelContact( + BaseSystemCollection, Constraints, Connections, CallBacks, Forcing + ): + pass + + rod_cylinder_parallel_contact_simulator = RodCylinderParallelContact() + + # time step etc + final_time = 10.0 + time_step = 1e-4 + total_steps = int(final_time / time_step) + 1 + rendering_fps = 30 # 20 * 1e1 + step_skip = int(1.0 / (rendering_fps * time_step)) + + base_length = 0.5 + base_radius = 0.1 + density = 1750 + E = 3e5 + poisson_ratio = 0.5 + shear_modulus = E / (2 * (1 + poisson_ratio)) + n_elem = 50 + nu = 0.5 + start = np.zeros((3,)) + direction = np.array([0, 0.0, 1.0]) + normal = np.array([0.0, 1.0, 0.0]) + + rod = CosseratRod.straight_rod( + n_elem, + start, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + shear_modulus=shear_modulus, + ) + + rod_cylinder_parallel_contact_simulator.append(rod) + # rod.velocity_collection[0, :] -= 0.2 + + # Apply uniform forces on the rod + rod_cylinder_parallel_contact_simulator.add_forcing_to(rod).using( + UniformForces, force=1.0 * pulling_force, direction=direction + ) + # Push the rod towards the cylinder to make sure contact is there + rod_cylinder_parallel_contact_simulator.add_forcing_to(rod).using( + UniformForces, force=0.1, direction=np.array([-1.0, 0.0, 0.0]) + ) + + cylinder_height = 3 * base_length + cylinder_radius = base_radius + + cylinder_start = start + np.array([-1.0, 0.0, 0.0]) * 2 * base_radius + cylinder_direction = np.array([0.0, 0.0, 1.0]) + cylinder_normal = np.array([0.0, 1.0, 0.0]) + + rigid_body = Cylinder( + start=cylinder_start, + direction=cylinder_direction, + normal=cylinder_normal, + base_length=cylinder_height, + base_radius=cylinder_radius, + density=density, + ) + rod_cylinder_parallel_contact_simulator.append(rigid_body) + + # Constrain the rigid body position and directors + rod_cylinder_parallel_contact_simulator.constrain(rigid_body).using( + OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,) + ) + + # Add contact between rigid body and rod + rod_cylinder_parallel_contact_simulator.connect(rod, rigid_body).using( + ExternalContact, + k=1e4, + nu=10, + velocity_damping_coefficient=10, + kinetic_friction_coefficient=10, + ) + + # Add callbacks + post_processing_dict_list = [] + # For rod + class StraightRodCallBack(CallBackBaseClass): + """ + Call back function for two arm octopus + """ + + 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["position"].append( + system.position_collection.copy() + ) + self.callback_params["radius"].append(system.radius.copy()) + self.callback_params["com"].append( + system.compute_position_center_of_mass() + ) + if current_step == 0: + self.callback_params["lengths"].append(system.rest_lengths.copy()) + else: + self.callback_params["lengths"].append(system.lengths.copy()) + + self.callback_params["com_velocity"].append( + system.compute_velocity_center_of_mass() + ) + + total_energy = ( + system.compute_translational_energy() + + system.compute_rotational_energy() + + system.compute_bending_energy() + + system.compute_shear_energy() + ) + self.callback_params["total_energy"].append(total_energy) + + return + + class RigidCylinderCallBack(CallBackBaseClass): + """ + Call back function for two arm octopus + """ + + def __init__( + self, step_skip: int, callback_params: dict, resize_cylinder_elems: int + ): + CallBackBaseClass.__init__(self) + self.every = step_skip + self.callback_params = callback_params + self.n_elem_cylinder = resize_cylinder_elems + self.n_node_cylinder = self.n_elem_cylinder + 1 + + 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) + + cylinder_center_position = system.position_collection + cylinder_length = system.length + cylinder_direction = system.director_collection[2, :, :].reshape(3, 1) + cylinder_radius = system.radius + + # Expand cylinder data. Create multiple points on cylinder later to use for rendering. + + start_position = ( + cylinder_center_position - cylinder_length / 2 * cylinder_direction + ) + + cylinder_position_collection = ( + start_position + + np.linspace(0, cylinder_length[0], self.n_node_cylinder) + * cylinder_direction + ) + cylinder_radius_collection = ( + np.ones((self.n_elem_cylinder)) * cylinder_radius + ) + cylinder_length_collection = ( + np.ones((self.n_elem_cylinder)) * cylinder_length + ) + cylinder_velocity_collection = ( + np.ones((self.n_node_cylinder)) * system.velocity_collection + ) + + self.callback_params["position"].append( + cylinder_position_collection.copy() + ) + self.callback_params["velocity"].append( + cylinder_velocity_collection.copy() + ) + self.callback_params["radius"].append(cylinder_radius_collection.copy()) + self.callback_params["com"].append( + system.compute_position_center_of_mass() + ) + + self.callback_params["lengths"].append( + cylinder_length_collection.copy() + ) + self.callback_params["com_velocity"].append( + system.velocity_collection[..., 0].copy() + ) + + total_energy = ( + system.compute_translational_energy() + + system.compute_rotational_energy() + ) + self.callback_params["total_energy"].append(total_energy[..., 0].copy()) + + return + + if POST_PROCESSING: + post_processing_dict_list.append(defaultdict(list)) + rod_cylinder_parallel_contact_simulator.collect_diagnostics(rod).using( + StraightRodCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_list[0], + ) + # For rigid body + post_processing_dict_list.append(defaultdict(list)) + rod_cylinder_parallel_contact_simulator.collect_diagnostics(rigid_body).using( + RigidCylinderCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_list[1], + resize_cylinder_elems=n_elem, + ) + + rod_cylinder_parallel_contact_simulator.finalize() + timestepper = PositionVerlet() + + integrate( + timestepper, rod_cylinder_parallel_contact_simulator, final_time, total_steps + ) + + if POST_PROCESSING: + # Plot the rods + plot_video_with_surface( + post_processing_dict_list, + video_name="rod_cylinder_contact.mp4", + fps=rendering_fps, + step=1, + # The following parameters are optional + x_limits=(-base_length * 5, base_length * 5), # Set bounds on x-axis + y_limits=(-base_length * 5, base_length * 5), # Set bounds on y-axis + z_limits=(-base_length * 5, base_length * 5), # Set bounds on z-axis + dpi=100, # Set the quality of the image + vis3D=True, # Turn on 3D visualization + vis2D=True, # Turn on projected (2D) visualization + ) + + filaname = "rod_rigid_velocity.png" + plot_velocity( + post_processing_dict_list[0], + post_processing_dict_list[1], + filename=filaname, + SAVE_FIGURE=True, + ) + + # Compute final total energy + total_final_energy = ( + rod.compute_translational_energy() + + rod.compute_rotational_energy() + + rod.compute_bending_energy() + + rod.compute_shear_energy() + ) + + return total_final_energy + + +if __name__ == "__main__": + import multiprocessing as mp + from examples.RigidBodyCases.RodRigidBodyContact.post_processing import ( + plot_force_vs_energy, + ) + + # total_energy = rod_cylinder_contact_friction_case(pulling_force=0.5, POST_PROCESSING=True) + + force = list([(float(x)) / 100.0 for x in range(0, 90, 5)]) + + with mp.Pool(mp.cpu_count()) as pool: + result_total_energy = pool.map(rod_cylinder_contact_friction_case, force) + + plot_force_vs_energy( + force, result_total_energy, filename="rod_energy_vs_force.png", SAVE_FIGURE=True + ) From f318dd81eb69bd70d4ad78b2594eca998768c453 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Mon, 4 Jul 2022 09:08:59 -0500 Subject: [PATCH 21/37] fix:connection docs, compatibility table Rigid body is not supported for fixed and free joint methods. Correct the compatibility table. --- docs/api/connections.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/api/connections.rst b/docs/api/connections.rst index 4d7503b3..42b8e3af 100644 --- a/docs/api/connections.rst +++ b/docs/api/connections.rst @@ -25,9 +25,9 @@ Compatibility =============================== ==== =========== Connection / Contact / Joints Rod Rigid Body =============================== ==== =========== -FreeJoint ✅ ✅ +FreeJoint ✅ ❌ ExternalContact ✅ ❌ -FixedJoint ✅ ✅ +FixedJoint ✅ ❌ HingeJoint ✅ ❌ SelfContact ✅ ❌ =============================== ==== =========== From 5b536176b89b762324ef03d4acec34fa30434a05 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20St=C3=B6lzle?= Date: Mon, 4 Jul 2022 17:47:40 +0200 Subject: [PATCH 22/37] Fix #128: save_every condition in ExportCallBack (#125) * Fix bug with `save_every` condition in ExportCallBack * minimizing any potential delays for saving into storage Co-authored-by: Seung Hyun Kim --- elastica/callback_functions.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/elastica/callback_functions.py b/elastica/callback_functions.py index 73fcada4..c200805f 100644 --- a/elastica/callback_functions.py +++ b/elastica/callback_functions.py @@ -200,11 +200,12 @@ def make_callback(self, system, time, current_step: int): + sys.getsizeof(velocity) + sys.getsizeof(director) ) - if ( - self.buffer_size > ExportCallBack.FILE_SIZE_CUTOFF - or (current_step + 1) % self.save_every == 0 - ): - self._dump() + + if ( + self.buffer_size > ExportCallBack.FILE_SIZE_CUTOFF + or (current_step + 1) % self.save_every == 0 + ): + self._dump() def _dump(self, **kwargs): file_path = f"{self.save_path}_{self.file_count}.dat" From 536a49e8a497ddcc37b60e3cf8f0df4e29c96432 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 4 Jul 2022 10:20:10 -0500 Subject: [PATCH 23/37] doc: clarify and refactor step_every, add destructor --- elastica/callback_functions.py | 29 ++++++++++++++++++++++------- tests/test_callback_functions.py | 2 +- 2 files changed, 23 insertions(+), 8 deletions(-) diff --git a/elastica/callback_functions.py b/elastica/callback_functions.py index c200805f..35033cfa 100644 --- a/elastica/callback_functions.py +++ b/elastica/callback_functions.py @@ -114,14 +114,16 @@ def __init__( path: str, method: str, initial_file_count: int = 0, - save_every: int = 1e8, + file_save_interval: int = 1e8, ): """ Parameters ---------- step_skip : int - Collect data at each step_skip steps. + Interval to collect simulation data into buffer. + The data will be collected at every `dt * step_skip` + interval. path : str Path to save the file. If directories are prepended, they must exist. The filename depends on the method. @@ -131,8 +133,9 @@ def __init__( allowed. initial_file_count : int Initial file count index that will be appended - save_every : int - Save the file every save_every steps. (default=1e8) + file_save_interval : int + Interval, in steps, to export/save collected buffer + as file. (default=1e8) """ # Assertions MIN_STEP_SKIP = 100 @@ -141,14 +144,15 @@ def __init__( assert ( method in ExportCallBack.AVAILABLE_METHOD ), f"The exporting method ({method}) is not supported. Please use one of {ExportCallBack.AVAILABLE_METHOD}." - assert os.path.exists(path), "The export path does not exist." + # TODO: Assertion is temporarily disabled. Should be fixed with #127 + # assert os.path.exists(path), "The export path does not exist." # Argument Parameters self.step_skip = step_skip self.save_path = path self.method = method self.file_count = initial_file_count - self.save_every = save_every + self.file_save_interval = file_save_interval # Data collector from collections import defaultdict @@ -203,11 +207,15 @@ def make_callback(self, system, time, current_step: int): if ( self.buffer_size > ExportCallBack.FILE_SIZE_CUTOFF - or (current_step + 1) % self.save_every == 0 + or (current_step + 1) % self.file_save_interval == 0 ): self._dump() def _dump(self, **kwargs): + """ + Dump dictionary buffer (self.buffer) to a file and clear + the buffer. + """ file_path = f"{self.save_path}_{self.file_count}.dat" data = {k: np.array(v) for k, v in self.buffer.items()} if self.method == ExportCallBack.AVAILABLE_METHOD[0]: @@ -225,3 +233,10 @@ def _dump(self, **kwargs): self.file_count += 1 self.buffer_size = 0 self.buffer.clear() + + def __del__(self): + """ + Save residual buffer on exit + """ + if self.buffer_size: + self._dump() diff --git a/tests/test_callback_functions.py b/tests/test_callback_functions.py index 1d578a74..036efa00 100644 --- a/tests/test_callback_functions.py +++ b/tests/test_callback_functions.py @@ -150,7 +150,7 @@ def test_export_call_back_base_class(self, n_elem): "directors": [], } - callback = ExportCallBack(step_skip, ".", "tempfile", save_every=10) + callback = ExportCallBack(step_skip, ".", "tempfile", file_save_interval=10) for i in range(10): callback.make_callback(mock_rod, time[i], current_step[i]) From aa5427a57463de21f31e74b9b92b64b5f1b24735 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Mon, 4 Jul 2022 17:18:11 -0500 Subject: [PATCH 24/37] refactor external contact with friction --- elastica/joint.py | 39 ++++++++++++++++++--------------------- 1 file changed, 18 insertions(+), 21 deletions(-) diff --git a/elastica/joint.py b/elastica/joint.py index 7390d6b1..6cd15a9e 100644 --- a/elastica/joint.py +++ b/elastica/joint.py @@ -378,7 +378,7 @@ def _calculate_contact_forces_rod_rigid_body( contact_k, contact_nu, velocity_damping_coefficient, - kinetic_friction_coefficient, + friction_coefficient, ): # We already pass in only the first n_elem x n_points = x_collection_rod.shape[1] @@ -446,25 +446,20 @@ def _calculate_contact_forces_rod_rigid_body( slip_interpenetration_velocity_mag = np.linalg.norm( slip_interpenetration_velocity ) - if slip_interpenetration_velocity_mag >= 1e-12: - slip_interpenetration_velocity_unitized = ( - slip_interpenetration_velocity / slip_interpenetration_velocity_mag - ) - else: - slip_interpenetration_velocity_unitized = np.zeros( - 3, - ) + slip_interpenetration_velocity_unitized = slip_interpenetration_velocity / ( + slip_interpenetration_velocity_mag + 1e-14 + ) # Compute friction force in the slip direction. damping_force_in_slip_direction = ( velocity_damping_coefficient * slip_interpenetration_velocity_mag ) - # Compute kinetic friction - kinetic_friction_force = kinetic_friction_coefficient * np.linalg.norm( + # Compute Coulombic friction + coulombic_friction_force = friction_coefficient * np.linalg.norm( net_contact_force ) # Compare damping force in slip direction and kinetic friction and minimum is the friction force. friction_force = ( - -min(damping_force_in_slip_direction, kinetic_friction_force) + -min(damping_force_in_slip_direction, coulombic_friction_force) * slip_interpenetration_velocity_unitized ) # Update contact force @@ -808,9 +803,13 @@ class ExternalContact(FreeJoint): is always cylinder. In addition to the contact forces, user can define apply friction forces between rod and cylinder that are in contact. For details on friction model refer to this `paper `_. - TODO: Currently friction force is between rod-cylinder, in future implement friction forces between rod-rod. + Notes + ----- + The `velocity_damping_coefficient` is set to a high value (e.g. 1e4) to minimize slip and simulate stiction + (static friction), while friction_coefficient corresponds to the Coulombic friction coefficient. + Examples -------- How to define contact between rod and cylinder. @@ -843,9 +842,7 @@ class ExternalContact(FreeJoint): """ - def __init__( - self, k, nu, velocity_damping_coefficient=0, kinetic_friction_coefficient=0 - ): + def __init__(self, k, nu, velocity_damping_coefficient=0, friction_coefficient=0): """ Parameters @@ -853,16 +850,16 @@ def __init__( k : float Contact spring constant. nu : float - Contact damping constant or coulomb coefficient of friction. + Contact damping constant. velocity_damping_coefficient : float Velocity damping coefficient between rigid-body and rod contact is used to apply friction force in the slip direction. - kinetic_friction_coefficient : float - Kinetic friction coefficient for rigid-body and rod contact. + friction_coefficient : float + For Coulombic friction coefficient for rigid-body and rod contact. """ super().__init__(k, nu) self.velocity_damping_coefficient = velocity_damping_coefficient - self.kinetic_friction_coefficient = kinetic_friction_coefficient + self.friction_coefficient = friction_coefficient def apply_forces(self, rod_one, index_one, rod_two, index_two): # del index_one, index_two @@ -912,7 +909,7 @@ def apply_forces(self, rod_one, index_one, rod_two, index_two): self.k, self.nu, self.velocity_damping_coefficient, - self.kinetic_friction_coefficient, + self.friction_coefficient, ) else: From fadc8596b37cf1cfc8c8e36dec04f1657914b341 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Mon, 4 Jul 2022 17:19:42 -0500 Subject: [PATCH 25/37] refactor:rigid-body rod contact friction examples --- .../RodRigidBodyContact/post_processing.py | 8 ++-- .../rod_cylinder_contact_friction.py | 44 ++++++------------- .../rod_cylinder_contact_friction_case.py | 8 ++++ ...d_cylinder_contact_friction_phase_space.py | 25 +++++++++++ 4 files changed, 52 insertions(+), 33 deletions(-) create mode 100644 examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction_case.py create mode 100644 examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction_phase_space.py diff --git a/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py b/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py index ebe9d83f..cc089926 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/post_processing.py @@ -728,8 +728,9 @@ def plot_video_with_surface( def plot_force_vs_energy( - force, + normalized_force, total_final_energy, + friction_coefficient, filename="energy_vs_force.png", SAVE_FIGURE=False, ): @@ -739,12 +740,13 @@ def plot_force_vs_energy( axs.append(plt.subplot2grid((1, 1), (0, 0))) axs[0].plot( - force, + normalized_force, total_final_energy, linewidth=3, ) + plt.axvline(x=friction_coefficient, linewidth=3, color="r", label="threshold") axs[0].set_ylabel("total energy", fontsize=20) - axs[0].set_xlabel("force ", fontsize=20) + axs[0].set_xlabel("normalized force", fontsize=20) plt.tight_layout() # fig.align_ylabels() diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction.py index b314eccd..7ad3e569 100644 --- a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction.py +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction.py @@ -6,7 +6,9 @@ from post_processing import plot_velocity, plot_video_with_surface -def rod_cylinder_contact_friction_case(pulling_force=0.0, POST_PROCESSING=False): +def rod_cylinder_contact_friction_case( + force_coefficient=0.1, normal_force_mag=10, POST_PROCESSING=False +): class RodCylinderParallelContact( BaseSystemCollection, Constraints, Connections, CallBacks, Forcing ): @@ -15,7 +17,7 @@ class RodCylinderParallelContact( rod_cylinder_parallel_contact_simulator = RodCylinderParallelContact() # time step etc - final_time = 10.0 + final_time = 20.0 time_step = 1e-4 total_steps = int(final_time / time_step) + 1 rendering_fps = 30 # 20 * 1e1 @@ -47,18 +49,18 @@ class RodCylinderParallelContact( ) rod_cylinder_parallel_contact_simulator.append(rod) - # rod.velocity_collection[0, :] -= 0.2 - # Apply uniform forces on the rod + # Push the rod towards the cylinder to make sure contact is there + normal_force_direction = np.array([-1.0, 0.0, 0.0]) rod_cylinder_parallel_contact_simulator.add_forcing_to(rod).using( - UniformForces, force=1.0 * pulling_force, direction=direction + UniformForces, force=normal_force_mag, direction=normal_force_direction ) - # Push the rod towards the cylinder to make sure contact is there + # Apply uniform forces on the rod rod_cylinder_parallel_contact_simulator.add_forcing_to(rod).using( - UniformForces, force=0.1, direction=np.array([-1.0, 0.0, 0.0]) + UniformForces, force=normal_force_mag * force_coefficient, direction=direction ) - cylinder_height = 3 * base_length + cylinder_height = 8 * base_length cylinder_radius = base_radius cylinder_start = start + np.array([-1.0, 0.0, 0.0]) * 2 * base_radius @@ -83,10 +85,10 @@ class RodCylinderParallelContact( # Add contact between rigid body and rod rod_cylinder_parallel_contact_simulator.connect(rod, rigid_body).using( ExternalContact, - k=1e4, - nu=10, - velocity_damping_coefficient=10, - kinetic_friction_coefficient=10, + k=1e5, + nu=100, + velocity_damping_coefficient=1e5, + friction_coefficient=0.5, ) # Add callbacks @@ -259,21 +261,3 @@ def make_callback(self, system, time, current_step: int): ) return total_final_energy - - -if __name__ == "__main__": - import multiprocessing as mp - from examples.RigidBodyCases.RodRigidBodyContact.post_processing import ( - plot_force_vs_energy, - ) - - # total_energy = rod_cylinder_contact_friction_case(pulling_force=0.5, POST_PROCESSING=True) - - force = list([(float(x)) / 100.0 for x in range(0, 90, 5)]) - - with mp.Pool(mp.cpu_count()) as pool: - result_total_energy = pool.map(rod_cylinder_contact_friction_case, force) - - plot_force_vs_energy( - force, result_total_energy, filename="rod_energy_vs_force.png", SAVE_FIGURE=True - ) diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction_case.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction_case.py new file mode 100644 index 00000000..c8dfe251 --- /dev/null +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction_case.py @@ -0,0 +1,8 @@ +if __name__ == "__main__": + from examples.RigidBodyCases.RodRigidBodyContact.rod_cylinder_contact_friction import ( + rod_cylinder_contact_friction_case, + ) + + total_energy = rod_cylinder_contact_friction_case( + force_coefficient=0.1, normal_force_mag=10, POST_PROCESSING=True + ) diff --git a/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction_phase_space.py b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction_phase_space.py new file mode 100644 index 00000000..7b7ae1d1 --- /dev/null +++ b/examples/RigidBodyCases/RodRigidBodyContact/rod_cylinder_contact_friction_phase_space.py @@ -0,0 +1,25 @@ +if __name__ == "__main__": + import multiprocessing as mp + from examples.RigidBodyCases.RodRigidBodyContact.rod_cylinder_contact_friction import ( + rod_cylinder_contact_friction_case, + ) + from examples.RigidBodyCases.RodRigidBodyContact.post_processing import ( + plot_force_vs_energy, + ) + + # total_energy = rod_cylinder_contact_friction_case(friction_coefficient=1.0, normal_force_mag=10, velocity_damping_coefficient=1E4, POST_PROCESSING=True) + + friction_coefficient = list([(float(x)) / 100.0 for x in range(0, 100, 5)]) + + with mp.Pool(mp.cpu_count()) as pool: + result_total_energy = pool.map( + rod_cylinder_contact_friction_case, friction_coefficient + ) + + plot_force_vs_energy( + friction_coefficient, + result_total_energy, + friction_coefficient=0.5, + filename="rod_energy_vs_force.png", + SAVE_FIGURE=True, + ) From f11be7f2bdbecbd02f527c4a7999d6f7b6ec973b Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sun, 10 Jul 2022 14:56:22 -0500 Subject: [PATCH 26/37] update: export callback parameter fix --- elastica/callback_functions.py | 35 +++++++++++++++++++++------------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/elastica/callback_functions.py b/elastica/callback_functions.py index 35033cfa..206c0275 100644 --- a/elastica/callback_functions.py +++ b/elastica/callback_functions.py @@ -5,6 +5,9 @@ import os import sys import numpy as np +import logging + +from collections import defaultdict class CallBackBaseClass: @@ -111,23 +114,25 @@ class ExportCallBack(CallBackBaseClass): def __init__( self, step_skip: int, - path: str, + filename: str, + directory: str, method: str, initial_file_count: int = 0, file_save_interval: int = 1e8, ): """ - Parameters ---------- step_skip : int Interval to collect simulation data into buffer. The data will be collected at every `dt * step_skip` interval. - path : str - Path to save the file. If directories are prepended, - they must exist. The filename depends on the method. - The path is not expected to include extension. + filename : str + Name of the file without extension. The extension will be + determined depend on the method. + directory : str + Path to save the file. If directory doesn't exist, it will + be created. Any existing files will be overwritten. method : str Method name. Only the name in AVAILABLE_METHOD is allowed. @@ -135,7 +140,7 @@ def __init__( Initial file count index that will be appended file_save_interval : int Interval, in steps, to export/save collected buffer - as file. (default=1e8) + as file. (default = 1e8) """ # Assertions MIN_STEP_SKIP = 100 @@ -144,19 +149,20 @@ def __init__( assert ( method in ExportCallBack.AVAILABLE_METHOD ), f"The exporting method ({method}) is not supported. Please use one of {ExportCallBack.AVAILABLE_METHOD}." - # TODO: Assertion is temporarily disabled. Should be fixed with #127 - # assert os.path.exists(path), "The export path does not exist." + + # Create directory + if os.path.exists(directory): + logging.warning(f"Save file already exists in {directory}. Files will be overwritten") + os.makedirs(directory, exist_ok=True) # Argument Parameters self.step_skip = step_skip - self.save_path = path + self.save_path = os.path.join(directory, filename) + "_{:02d}.{}" self.method = method self.file_count = initial_file_count self.file_save_interval = file_save_interval # Data collector - from collections import defaultdict - self.buffer = defaultdict(list) self.buffer_size = 0 @@ -165,16 +171,19 @@ def __init__( import pickle self._pickle = pickle + self._ext = "pkl" elif method == ExportCallBack.AVAILABLE_METHOD[1]: from numpy import savez self._savez = savez + self._ext = "npz" elif method == ExportCallBack.AVAILABLE_METHOD[2]: import tempfile import pickle self._tempfile = tempfile.NamedTemporaryFile(delete=False) self._pickle = pickle + self._ext = "pkl" def make_callback(self, system, time, current_step: int): """ @@ -216,7 +225,7 @@ def _dump(self, **kwargs): Dump dictionary buffer (self.buffer) to a file and clear the buffer. """ - file_path = f"{self.save_path}_{self.file_count}.dat" + file_path = self.save_path.format(self.file_count, self._ext) data = {k: np.array(v) for k, v in self.buffer.items()} if self.method == ExportCallBack.AVAILABLE_METHOD[0]: # pickle From 350c54b83cf8cb7939be4f7261f845d7bf17bc89 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 11 Jul 2022 00:42:34 -0500 Subject: [PATCH 27/37] test: add more test for export file --- elastica/callback_functions.py | 28 ++- tests/test_callback_functions.py | 284 ++++++++++++++++++++++++------- 2 files changed, 247 insertions(+), 65 deletions(-) diff --git a/elastica/callback_functions.py b/elastica/callback_functions.py index 206c0275..63cbc794 100644 --- a/elastica/callback_functions.py +++ b/elastica/callback_functions.py @@ -129,7 +129,7 @@ def __init__( interval. filename : str Name of the file without extension. The extension will be - determined depend on the method. + determined depend on the method. directory : str Path to save the file. If directory doesn't exist, it will be created. Any existing files will be overwritten. @@ -152,7 +152,9 @@ def __init__( # Create directory if os.path.exists(directory): - logging.warning(f"Save file already exists in {directory}. Files will be overwritten") + logging.warning( + f"Save file already exists in {directory}. Files will be overwritten" + ) os.makedirs(directory, exist_ok=True) # Argument Parameters @@ -215,7 +217,7 @@ def make_callback(self, system, time, current_step: int): ) if ( - self.buffer_size > ExportCallBack.FILE_SIZE_CUTOFF + self.buffer_size > self.FILE_SIZE_CUTOFF or (current_step + 1) % self.file_save_interval == 0 ): self._dump() @@ -243,9 +245,25 @@ def _dump(self, **kwargs): self.buffer_size = 0 self.buffer.clear() - def __del__(self): + def get_last_saved_path(self) -> str: + """ + Return last saved file path. If no file has been saved, + return None + """ + if self.file_count == 0: + return None + else: + return self.save_path.format(self.file_count - 1, self._ext) + + def close(self): """ - Save residual buffer on exit + Save residual buffer """ if self.buffer_size: self._dump() + + def clear(self): + """ + Alias to `close` + """ + self.close() diff --git a/tests/test_callback_functions.py b/tests/test_callback_functions.py index 036efa00..38c25faa 100644 --- a/tests/test_callback_functions.py +++ b/tests/test_callback_functions.py @@ -1,24 +1,33 @@ __doc__ = """ Call back functions for rod test module """ -import sys +import os, sys # System imports import numpy as np from numpy.testing import assert_allclose, assert_array_equal from elastica.callback_functions import CallBackBaseClass, MyCallBack, ExportCallBack from elastica.utils import Tolerance +import tempfile import pytest -def mock_rod_init(self): - self.n_elems = 0.0 - self.position_collection = 0.0 - self.velocity_collection = 0.0 - self.director_collection = 0.0 - self.external_forces = 0.0 - self.external_torques = 0.0 +class MockRod: + def __init__(self): + self.n_elems = 0.0 + self.position_collection = 0.0 + self.velocity_collection = 0.0 + self.director_collection = 0.0 + self.external_forces = 0.0 + self.external_torques = 0.0 -MockRod = type("MockRod", (object,), {"__init__": mock_rod_init}) +class MockRodWithElements: + def __init__(self, n_elems): + self.n_elems = n_elems + self.position_collection = np.random.rand(3, n_elems) + self.velocity_collection = np.random.rand(3, n_elems) + self.director_collection = np.random.rand(3, 3, n_elems) + self.external_forces = np.random.rand(3, n_elems) + self.external_torques = np.random.rand(3, n_elems) class TestCallBackBaseClass: @@ -26,9 +35,6 @@ def test_call_back_base_class(self): """ This test case tests, base class for call functions. make_callback does not do anything, but for completion this test case is here. - Returns - ------- - """ mock_rod = MockRod() @@ -46,32 +52,15 @@ def test_call_back_base_class(self): class TestMyCallBackClass: - @pytest.mark.parametrize("n_elem", [2, 4, 16]) - def test_my_call_back_base_class(self, n_elem): + @pytest.mark.parametrize("n_elems", [2, 4, 16]) + def test_my_call_back_base_class(self, n_elems): """ This test case is for testing MyCallBack function. - Parameters - ---------- - n_elem - - Returns - ------- - """ - - mock_rod = MockRod() + mock_rod = MockRodWithElements(n_elems) time = np.random.rand(10) current_step = list(range(10)) - position_collection = np.random.rand(3, 10) - velocity_collection = np.random.rand(3, 10) - director_collection = np.random.rand(3, 3, 10) - - # set arrays in mock rod - mock_rod.n_elems = n_elem - mock_rod.position_collection = position_collection - mock_rod.velocity_collection = velocity_collection - mock_rod.director_collection = director_collection step_skip = 1 list_test = { @@ -95,9 +84,9 @@ def test_my_call_back_base_class(self, n_elem): list_correct["time"].append(time[i]) list_correct["step"].append(current_step[i]) - list_correct["position"].append(position_collection) - list_correct["velocity"].append(velocity_collection) - list_correct["directors"].append(director_collection) + list_correct["position"].append(mock_rod.position_collection.copy()) + list_correct["velocity"].append(mock_rod.velocity_collection.copy()) + list_correct["directors"].append(mock_rod.director_collection.copy()) assert_allclose(list_test["time"], list_correct["time"], atol=Tolerance.atol()) assert_allclose(list_test["step"], list_correct["step"], atol=Tolerance.atol()) @@ -113,33 +102,99 @@ def test_my_call_back_base_class(self, n_elem): class TestExportCallBackClass: - @pytest.mark.parametrize("n_elem", [2, 4, 16]) - def test_export_call_back_base_class(self, n_elem): + @pytest.mark.parametrize("method", ["0", 1, "numba", "test", "some string", None]) + def test_export_call_back_unavailable_save_methods(self, method): + with pytest.raises(AssertionError) as excinfo: + callback = ExportCallBack(1, "rod", "tempdir", method) + + @pytest.mark.parametrize("method", ExportCallBack.AVAILABLE_METHOD) + def test_export_call_back_available_save_methods(self, method): + try: + callback = ExportCallBack(1, "rod", "tempdir", method) + except Error: + pytest.fail( + f"Could not create callback module with available method {method}" + ) + + def test_export_call_back_interval_by_filesize(self): + mock_rod = MockRodWithElements(5) + with tempfile.TemporaryDirectory() as temp_dir_path: + callback = ExportCallBack(1, "rod", temp_dir_path, "npz") + callback.FILE_SIZE_CUTOFF = 1 + callback.make_callback(mock_rod, 1, 1) + + saved_path_name = callback.save_path.format(0, "npz") + assert os.path.exists(saved_path_name), "File is not saved." + + @pytest.mark.parametrize("file_save_interval", [1, 5, 10, 15]) + def test_export_call_back_file_save_interval_param(self, file_save_interval): + mock_rod = MockRodWithElements(5) + with tempfile.TemporaryDirectory() as temp_dir_path: + callback = ExportCallBack( + 1, "rod", temp_dir_path, "npz", file_save_interval=file_save_interval + ) + for step in range(file_save_interval): + callback.make_callback(mock_rod, 1, step) + + saved_path_name = callback.get_last_saved_path() + assert os.path.exists(saved_path_name), "File is not saved." + + @pytest.mark.parametrize("file_save_interval", [5, 10]) + def test_export_call_back_file_save_interval_param_ext(self, file_save_interval): + mock_rod = MockRodWithElements(5) + n_repeat = 3 + with tempfile.TemporaryDirectory() as temp_dir_path: + callback = ExportCallBack( + 1, "rod", temp_dir_path, "npz", file_save_interval=file_save_interval + ) + for step in range(file_save_interval * n_repeat): + callback.make_callback(mock_rod, 1, step) + + saved_path_name = callback.get_last_saved_path() + assert os.path.exists(saved_path_name), "File is not saved." + + def test_export_call_back_file_not_saved(self): + with tempfile.TemporaryDirectory() as temp_dir_path: + callback = ExportCallBack( + 1, "rod", temp_dir_path, "npz", file_save_interval=10 + ) + saved_path_name = callback.get_last_saved_path() + assert saved_path_name is None, f"{saved_path_name} should be None" + + def test_export_call_back_close_test(self): + mock_rod = MockRodWithElements(5) + with tempfile.TemporaryDirectory() as temp_dir_path: + callback = ExportCallBack( + 1, "rod", temp_dir_path, "npz", file_save_interval=50 + ) + for step in range(10): + callback.make_callback(mock_rod, 1, step) + callback.close() + saved_path_name = callback.get_last_saved_path() + assert os.path.exists(saved_path_name), "File is not saved." + + def test_export_call_back_clear_test(self): + mock_rod = MockRodWithElements(5) + with tempfile.TemporaryDirectory() as temp_dir_path: + callback = ExportCallBack( + 1, "rod", temp_dir_path, "npz", file_save_interval=50 + ) + for step in range(10): + callback.make_callback(mock_rod, 1, step) + callback.clear() + saved_path_name = callback.get_last_saved_path() + assert os.path.exists(saved_path_name), "File is not saved." + + @pytest.mark.parametrize("n_elems", [2, 4, 16]) + def test_export_call_back_class_tempfile_option(self, n_elems): """ - This test case is for testing ExportCallBack function. - Parameters - ---------- - n_elem - - Returns - ------- - + This test case is for testing ExportCallBack function, saving into temporary files. """ import pickle - mock_rod = MockRod() - + mock_rod = MockRodWithElements(n_elems) time = np.random.rand(10) current_step = list(range(10)) - position_collection = np.random.rand(3, 10) - velocity_collection = np.random.rand(3, 10) - director_collection = np.random.rand(3, 3, 10) - - # set arrays in mock rod - mock_rod.n_elems = n_elem - mock_rod.position_collection = position_collection - mock_rod.velocity_collection = velocity_collection - mock_rod.director_collection = director_collection step_skip = 1 list_correct = { @@ -150,15 +205,17 @@ def test_export_call_back_base_class(self, n_elem): "directors": [], } - callback = ExportCallBack(step_skip, ".", "tempfile", file_save_interval=10) + callback = ExportCallBack( + step_skip, "rod", "tempdir", "tempfile", file_save_interval=10 + ) for i in range(10): callback.make_callback(mock_rod, time[i], current_step[i]) list_correct["time"].append(time[i]) list_correct["step"].append(current_step[i]) - list_correct["position"].append(position_collection) - list_correct["velocity"].append(velocity_collection) - list_correct["directors"].append(director_collection) + list_correct["position"].append(mock_rod.position_collection) + list_correct["velocity"].append(mock_rod.velocity_collection) + list_correct["directors"].append(mock_rod.director_collection) file = open(callback._tempfile.name, "rb") list_test = pickle.load(file) @@ -175,3 +232,110 @@ def test_export_call_back_base_class(self, n_elem): ) callback._tempfile.close() + + @pytest.mark.parametrize("n_elems", [2, 4, 16]) + def test_export_call_back_class_npz_option(self, n_elems): + """ + This test case is for testing ExportCallBack function, saving into numpy files. + """ + filename = "test_rod" + mock_rod = MockRodWithElements(n_elems) + time = np.random.rand(10) + current_step = list(range(10)) + + step_skip = 1 + list_correct = { + "time": [], + "step": [], + "position": [], + "velocity": [], + "directors": [], + } + + with tempfile.TemporaryDirectory() as temp_dir_path: + callback = ExportCallBack( + step_skip, filename, temp_dir_path, "npz", file_save_interval=10 + ) + for i in range(10): + callback.make_callback(mock_rod, time[i], current_step[i]) + + list_correct["time"].append(time[i]) + list_correct["step"].append(current_step[i]) + list_correct["position"].append(mock_rod.position_collection) + list_correct["velocity"].append(mock_rod.velocity_collection) + list_correct["directors"].append(mock_rod.director_collection) + + saved_path_name = callback.get_last_saved_path() + assert os.path.exists(saved_path_name), "File does not exist" + list_test = np.load(saved_path_name) + + assert_allclose( + list_test["time"], list_correct["time"], atol=Tolerance.atol() + ) + assert_allclose( + list_test["step"], list_correct["step"], atol=Tolerance.atol() + ) + assert_allclose( + list_test["position"], list_correct["position"], atol=Tolerance.atol() + ) + assert_allclose( + list_test["velocity"], list_correct["velocity"], atol=Tolerance.atol() + ) + assert_allclose( + list_test["directors"], list_correct["directors"], atol=Tolerance.atol() + ) + + @pytest.mark.parametrize("n_elems", [2, 4, 16]) + def test_export_call_back_class_pickle_option(self, n_elems): + """ + This test case is for testing ExportCallBack function, saving into pickle files. + """ + import pickle + + filename = "test_rod" + mock_rod = MockRodWithElements(n_elems) + time = np.random.rand(10) + current_step = list(range(10)) + + step_skip = 1 + list_correct = { + "time": [], + "step": [], + "position": [], + "velocity": [], + "directors": [], + } + + with tempfile.TemporaryDirectory() as temp_dir_path: + callback = ExportCallBack( + step_skip, filename, temp_dir_path, "pickle", file_save_interval=10 + ) + for i in range(10): + callback.make_callback(mock_rod, time[i], current_step[i]) + + list_correct["time"].append(time[i]) + list_correct["step"].append(current_step[i]) + list_correct["position"].append(mock_rod.position_collection) + list_correct["velocity"].append(mock_rod.velocity_collection) + list_correct["directors"].append(mock_rod.director_collection) + + saved_path_name = callback.get_last_saved_path() + assert os.path.exists(saved_path_name), "File does not exist" + file = open(saved_path_name, "rb") + list_test = pickle.load(file) + + assert_allclose( + list_test["time"], list_correct["time"], atol=Tolerance.atol() + ) + assert_allclose( + list_test["step"], list_correct["step"], atol=Tolerance.atol() + ) + assert_allclose( + list_test["position"], list_correct["position"], atol=Tolerance.atol() + ) + assert_allclose( + list_test["velocity"], list_correct["velocity"], atol=Tolerance.atol() + ) + assert_allclose( + list_test["directors"], list_correct["directors"], atol=Tolerance.atol() + ) From 6ed52399862219991dd38e1b313e9aa7d504cb7a Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 11 Jul 2022 00:43:30 -0500 Subject: [PATCH 28/37] add global path for mock import --- tests/__init__.py | 0 tests/test_boundary_conditions.py | 2 +- tests/test_interaction.py | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) create mode 100644 tests/__init__.py diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/test_boundary_conditions.py b/tests/test_boundary_conditions.py index 7e3144c5..af289e77 100644 --- a/tests/test_boundary_conditions.py +++ b/tests/test_boundary_conditions.py @@ -3,7 +3,7 @@ # System imports import numpy as np -from test_rod.test_rods import MockTestRod +from tests.test_rod.test_rods import MockTestRod from elastica.boundary_conditions import ( ConstraintBase, FreeBC, diff --git a/tests/test_interaction.py b/tests/test_interaction.py index 82ac01a8..8f60494e 100644 --- a/tests/test_interaction.py +++ b/tests/test_interaction.py @@ -12,7 +12,7 @@ SlenderBodyTheory, ) -from test_rod.test_rods import MockTestRod +from tests.test_rod.test_rods import MockTestRod class BaseRodClass(MockTestRod): From 0faa029886c48868e1d634567dffec38411a465e Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 11 Jul 2022 01:07:30 -0500 Subject: [PATCH 29/37] test: add warning msg test --- elastica/callback_functions.py | 3 +-- tests/test_callback_functions.py | 22 ++++++++++++++++++++++ 2 files changed, 23 insertions(+), 2 deletions(-) diff --git a/elastica/callback_functions.py b/elastica/callback_functions.py index 63cbc794..098f5052 100644 --- a/elastica/callback_functions.py +++ b/elastica/callback_functions.py @@ -1,7 +1,6 @@ __doc__ = """ Module contains callback classes to save simulation data for rod-like objects """ __all__ = ["CallBackBaseClass", "MyCallBack", "ExportCallBack"] -import warnings import os import sys import numpy as np @@ -145,7 +144,7 @@ def __init__( # Assertions MIN_STEP_SKIP = 100 if step_skip <= MIN_STEP_SKIP: - warnings.warn(f"We recommend step_skip at least {MIN_STEP_SKIP}") + logging.warning(f"We recommend ({step_skip=}) at least {MIN_STEP_SKIP}") assert ( method in ExportCallBack.AVAILABLE_METHOD ), f"The exporting method ({method}) is not supported. Please use one of {ExportCallBack.AVAILABLE_METHOD}." diff --git a/tests/test_callback_functions.py b/tests/test_callback_functions.py index 38c25faa..fff6cf64 100644 --- a/tests/test_callback_functions.py +++ b/tests/test_callback_functions.py @@ -2,6 +2,7 @@ import os, sys # System imports +import logging import numpy as np from numpy.testing import assert_allclose, assert_array_equal from elastica.callback_functions import CallBackBaseClass, MyCallBack, ExportCallBack @@ -116,6 +117,27 @@ def test_export_call_back_available_save_methods(self, method): f"Could not create callback module with available method {method}" ) + @pytest.mark.parametrize("step_skip", [1, 5, 10, 99]) + def test_export_call_back_small_stepsize_warning(self, caplog, step_skip): + mock_rod = MockRodWithElements(5) + with tempfile.TemporaryDirectory() as temp_dir_path: + ExportCallBack(step_skip, "rod", temp_dir_path, "npz") + record_tuple = caplog.record_tuples[0] + assert record_tuple[0] == "root" + assert record_tuple[1] == logging.WARNING + assert str(100) in record_tuple[2] + assert f"recommend ({step_skip=}) at least" in record_tuple[2] + + def test_export_call_back_file_recreate_warning(self, caplog): + mock_rod = MockRodWithElements(5) + with tempfile.TemporaryDirectory() as temp_dir_path: + ExportCallBack(1000, "rod", temp_dir_path, "npz") + ExportCallBack(1000, "rod", temp_dir_path, "npz") + record_tuple = caplog.record_tuples[0] + assert record_tuple[0] == "root" + assert record_tuple[1] == logging.WARNING + assert "already exists" in record_tuple[2] + def test_export_call_back_interval_by_filesize(self): mock_rod = MockRodWithElements(5) with tempfile.TemporaryDirectory() as temp_dir_path: From 38f66930ad4e5e7ab92a4a214b76ab9bc35b12e8 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 11 Jul 2022 01:14:10 -0500 Subject: [PATCH 30/37] test: add step_skip test --- tests/test_callback_functions.py | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/tests/test_callback_functions.py b/tests/test_callback_functions.py index fff6cf64..6dbd7b47 100644 --- a/tests/test_callback_functions.py +++ b/tests/test_callback_functions.py @@ -117,7 +117,7 @@ def test_export_call_back_available_save_methods(self, method): f"Could not create callback module with available method {method}" ) - @pytest.mark.parametrize("step_skip", [1, 5, 10, 99]) + @pytest.mark.parametrize("step_skip", [2, 5, 20, 50, 99]) def test_export_call_back_small_stepsize_warning(self, caplog, step_skip): mock_rod = MockRodWithElements(5) with tempfile.TemporaryDirectory() as temp_dir_path: @@ -161,6 +161,34 @@ def test_export_call_back_file_save_interval_param(self, file_save_interval): saved_path_name = callback.get_last_saved_path() assert os.path.exists(saved_path_name), "File is not saved." + @pytest.mark.parametrize("step_skip", [2, 5, 10, 15]) + def test_export_call_back_step_skip_param(self, step_skip): + mock_rod = MockRodWithElements(5) + with tempfile.TemporaryDirectory() as temp_dir_path: + callback = ExportCallBack(step_skip, "rod", temp_dir_path, "npz") + callback.make_callback(mock_rod, 1, step_skip - 1) + # Check empty + callback.clear() + saved_path_name = callback.get_last_saved_path() + assert saved_path_name is None, "No file should be saved." + + # Check saved + callback.make_callback(mock_rod, 1, step_skip) + callback.clear() + saved_path_name = callback.get_last_saved_path() + assert saved_path_name is not None, "File should be saved." + assert os.path.exists(saved_path_name), "File should be saved" + + # Check saved file number + callback.make_callback(mock_rod, 1, step_skip * 2) + callback.clear() + callback.make_callback(mock_rod, 1, step_skip * 5) + callback.clear() + saved_path_name = callback.get_last_saved_path() + assert ( + str(2) in saved_path_name + ), f"Total 3 file should be saved: {saved_path_name}" + @pytest.mark.parametrize("file_save_interval", [5, 10]) def test_export_call_back_file_save_interval_param_ext(self, file_save_interval): mock_rod = MockRodWithElements(5) From 9009fa69dc0a062d2ecf6819e8ca973a044e401c Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 11 Jul 2022 01:17:47 -0500 Subject: [PATCH 31/37] test: remove unsupported f-string format --- elastica/callback_functions.py | 4 +++- tests/test_callback_functions.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/elastica/callback_functions.py b/elastica/callback_functions.py index 098f5052..976458e7 100644 --- a/elastica/callback_functions.py +++ b/elastica/callback_functions.py @@ -144,7 +144,9 @@ def __init__( # Assertions MIN_STEP_SKIP = 100 if step_skip <= MIN_STEP_SKIP: - logging.warning(f"We recommend ({step_skip=}) at least {MIN_STEP_SKIP}") + logging.warning( + f"We recommend (step_skip={step_skip}) at least {MIN_STEP_SKIP}" + ) assert ( method in ExportCallBack.AVAILABLE_METHOD ), f"The exporting method ({method}) is not supported. Please use one of {ExportCallBack.AVAILABLE_METHOD}." diff --git a/tests/test_callback_functions.py b/tests/test_callback_functions.py index 6dbd7b47..80300b14 100644 --- a/tests/test_callback_functions.py +++ b/tests/test_callback_functions.py @@ -126,7 +126,7 @@ def test_export_call_back_small_stepsize_warning(self, caplog, step_skip): assert record_tuple[0] == "root" assert record_tuple[1] == logging.WARNING assert str(100) in record_tuple[2] - assert f"recommend ({step_skip=}) at least" in record_tuple[2] + assert f"recommend (step_skip={step_skip}) at least" in record_tuple[2] def test_export_call_back_file_recreate_warning(self, caplog): mock_rod = MockRodWithElements(5) From dc8b5bacf4b17e1d219f12728876705836b8f52b Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Tue, 12 Jul 2022 17:31:01 -0500 Subject: [PATCH 32/37] Update RELEASE.md Add release note for v0.2.4 --- RELEASE.md | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/RELEASE.md b/RELEASE.md index c2248082..f6493fc3 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -1,5 +1,21 @@ -# Release Note (version 0.2.3) +# Release Note (version 0.2.4) + +## What's Changed + +* Refactor EndPointForcesSinusoidal example and test cases by @armantekinalp in https://github.com/GazzolaLab/PyElastica/pull/110 +* Fix save_every condition in ExportCallBack by @mstoelzle in https://github.com/GazzolaLab/PyElastica/pull/125 +* Fix and update contact examples by @armantekinalp in https://github.com/GazzolaLab/PyElastica/pull/109 +* Update rigid body rod contact by @armantekinalp in https://github.com/GazzolaLab/PyElastica/pull/117 +* Update rigid body rod contact friction by @armantekinalp in https://github.com/GazzolaLab/PyElastica/pull/124 +* Update ExportCallback by @skim0119 in https://github.com/GazzolaLab/PyElastica/pull/130 +## New Contributors + +* @mstoelzle made their first contribution in https://github.com/GazzolaLab/PyElastica/pull/125 + +--- + +# Release Note (version 0.2.3) ## Developer Note The major updates are knot theory module added to the Cosserat rod as *mixin*, and muscular snake example is added. @@ -8,6 +24,7 @@ The major updates are knot theory module added to the Cosserat rod as *mixin*, a - #70: Knot theory module to compute topological quantities. - #71: Reorganize rod constructor warning messages and collect messages in log. - #72: Muscular snake example is added. + --- # Release Note (version 0.2.2) From 427a1bf89e5b5827b310b23a84726757cd6a1f42 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Tue, 12 Jul 2022 21:36:55 -0500 Subject: [PATCH 33/37] bump version number to 0.2.4 --- elastica/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elastica/version.py b/elastica/version.py index ae50f1f6..4b51ef6f 100644 --- a/elastica/version.py +++ b/elastica/version.py @@ -1 +1 @@ -VERSION = "0.2.3" +VERSION = "0.2.4" From 71f6b73d7fddf7940d0fda28483d17463bfdd6fb Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Tue, 12 Jul 2022 21:46:58 -0500 Subject: [PATCH 34/37] fix:bump version number in setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 142d458e..386b3461 100644 --- a/setup.py +++ b/setup.py @@ -14,7 +14,7 @@ EMAIL = "armant2@illinois.edu" AUTHOR = "GazzolaLab" REQUIRES_PYTHON = ">=3.6.0" -VERSION = "0.2.3" +VERSION = "0.2.4" # What packages are required for this module to be executed? REQUIRED = ["numpy>=1.19.2", "matplotlib>=3.3.2", "scipy>=1.5.2", "tqdm>=4.61.1", "numba>=0.51.0"] From 91e5371c91ca7f2ffc4c72f2df4fac78579cfe40 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Tue, 12 Jul 2022 22:09:07 -0500 Subject: [PATCH 35/37] update:gitworkflow to run tests on macos and windows --- .github/workflows/main.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 538296b7..1ba10769 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -19,7 +19,7 @@ jobs: strategy: matrix: - os: [ubuntu-latest] # [ubuntu-latest, macos-latest] # Run macos tests if really required, since they charge 10 times more for macos + os: [ubuntu-latest, macos-latest, windows-latest] # Run macos tests if really required, since they charge 10 times more for macos python-version: [3.6, 3.7, 3.8] # Steps represent a sequence of tasks that will be executed as part of the job From 42ece6b292bc8930567ee02f4b42f0c210abb8ad Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Wed, 13 Jul 2022 13:58:19 -0500 Subject: [PATCH 36/37] fix:github workflow remove windows from the os list --- .github/workflows/main.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 1ba10769..c85b14c8 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -19,7 +19,7 @@ jobs: strategy: matrix: - os: [ubuntu-latest, macos-latest, windows-latest] # Run macos tests if really required, since they charge 10 times more for macos + os: [ubuntu-latest, macos-latest] #, windows-latest] # Run macos tests if really required, since they charge 10 times more for macos python-version: [3.6, 3.7, 3.8] # Steps represent a sequence of tasks that will be executed as part of the job From 9e16e3d484d9c30673cc00897e5d06b091db08c4 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Thu, 14 Jul 2022 23:03:32 -0500 Subject: [PATCH 37/37] Update callback_functions.py --- elastica/callback_functions.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/elastica/callback_functions.py b/elastica/callback_functions.py index 976458e7..09f5b8aa 100644 --- a/elastica/callback_functions.py +++ b/elastica/callback_functions.py @@ -128,10 +128,12 @@ def __init__( interval. filename : str Name of the file without extension. The extension will be - determined depend on the method. + determined depend on the method. File will be saved with the + name _.. directory : str - Path to save the file. If directory doesn't exist, it will - be created. Any existing files will be overwritten. + Directory to save the file. If directory doesn't exist, it will + be created. During the save, any existing files in this directory + could be overwritten. method : str Method name. Only the name in AVAILABLE_METHOD is allowed. @@ -154,7 +156,7 @@ def __init__( # Create directory if os.path.exists(directory): logging.warning( - f"Save file already exists in {directory}. Files will be overwritten" + f"The directory ({directory}) already exists. Previously saved data could be overwritten." ) os.makedirs(directory, exist_ok=True)