From 51d2ba92257ae67e94cd91b3ff1002681f62ef0b Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Mon, 24 Aug 2020 08:08:44 -0500 Subject: [PATCH 01/14] enhencement:numba version of rod-rod and self cont This commit adds, numba versions of rod-rod and self contact implementations and also test cases --- elastica/_elastica_numba/_joint.py | 432 ++++++++++++++++-- elastica/joint.py | 4 +- .../rod_rod_contact_inclined_validation.py | 164 +++++++ .../rod_rod_contact_parallel_validation.py | 161 +++++++ .../rod_self_contact_validation.py | 128 ++++++ examples/RodContactCase/post_processing.py | 366 +++++++++++++++ 6 files changed, 1220 insertions(+), 35 deletions(-) create mode 100644 examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py create mode 100644 examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py create mode 100644 examples/RodContactCase/RodSelfContact/rod_self_contact_validation.py create mode 100644 examples/RodContactCase/post_processing.py diff --git a/elastica/_elastica_numba/_joint.py b/elastica/_elastica_numba/_joint.py index 4235f974..ced26912 100644 --- a/elastica/_elastica_numba/_joint.py +++ b/elastica/_elastica_numba/_joint.py @@ -5,6 +5,7 @@ from elastica.joint import FreeJoint from math import sqrt +from elastica._linalg import _batch_product_k_ik_to_ik @numba.njit(cache=True) @@ -97,7 +98,7 @@ def _find_min_dist(x1, e1, x2, e2): @numba.njit(cache=True) -def _calculate_contact_forces( +def _calculate_contact_forces_rod_rigid_body( x_collection_rod, edge_collection_rod, x_cylinder, @@ -184,6 +185,252 @@ def _calculate_contact_forces( external_forces_cylinder[..., 0] += 2.0 * net_contact_force +@numba.njit(cache=True) +def _calculate_contact_forces_rod_rod( + x_collection_rod_one, + radius_rod_one, + length_rod_one, + tangent_rod_one, + velocity_rod_one, + internal_forces_rod_one, + external_forces_rod_one, + x_collection_rod_two, + radius_rod_two, + length_rod_two, + tangent_rod_two, + velocity_rod_two, + internal_forces_rod_two, + external_forces_rod_two, + contact_k, + contact_nu, +): + # We already pass in only the first n_elem x + n_points_rod_one = x_collection_rod_one.shape[1] + n_points_rod_two = x_collection_rod_two.shape[1] + edge_collection_rod_one = _batch_product_k_ik_to_ik(length_rod_one, tangent_rod_one) + edge_collection_rod_two = _batch_product_k_ik_to_ik(length_rod_two, tangent_rod_two) + + for i in range(n_points_rod_one): + for j in range(n_points_rod_two): + radii_sum = radius_rod_one[i] + radius_rod_two[j] + length_sum = length_rod_one[i] + length_rod_two[j] + # Element-wise bounding box + x_selected_rod_one = x_collection_rod_one[..., i] + x_selected_rod_two = x_collection_rod_two[..., j] + + del_x = x_selected_rod_one - x_selected_rod_two + norm_del_x = _norm(del_x) + + # If outside then don't process + if norm_del_x >= (radii_sum + length_sum): + continue + + # find the shortest line segment between the two centerline + # segments : differs from normal cylinder-cylinder intersection + distance_vector = _find_min_dist( + x_selected_rod_one, + edge_collection_rod_one[..., i], + x_selected_rod_two, + edge_collection_rod_two[..., j], + ) + distance_vector_length = _norm(distance_vector) + distance_vector /= distance_vector_length + + gamma = radii_sum - distance_vector_length + + # If distance is large, don't worry about it + if gamma < -1e-5: + continue + + rod_one_elemental_forces = 0.5 * ( + external_forces_rod_one[..., i] + + external_forces_rod_one[..., i + 1] + + internal_forces_rod_one[..., i] + + internal_forces_rod_one[..., i + 1] + ) + + rod_two_elemental_forces = 0.5 * ( + external_forces_rod_two[..., j] + + external_forces_rod_two[..., j + 1] + + internal_forces_rod_two[..., j] + + internal_forces_rod_two[..., j + 1] + ) + + equilibrium_forces = -rod_one_elemental_forces + rod_two_elemental_forces + + normal_force = _dot_product(equilibrium_forces, distance_vector) + # Following line same as np.where(normal_force < 0.0, -normal_force, 0.0) + normal_force = abs(min(normal_force, 0.0)) + + # CHECK FOR GAMMA > 0.0, heaviside but we need to overload it in numba + # As a quick fix, use this instead + mask = (gamma > 0.0) * 1.0 + + contact_force = contact_k * gamma + interpenetration_velocity = 0.5 * ( + (velocity_rod_one[..., i] + velocity_rod_one[..., i + 1]) + - (velocity_rod_two[..., j] + velocity_rod_two[..., j + 1]) + ) + contact_damping_force = contact_nu * _dot_product( + interpenetration_velocity, distance_vector + ) + + # magnitude* direction + net_contact_force = ( + normal_force + 0.5 * mask * (contact_damping_force + contact_force) + ) * distance_vector + + # Add it to the rods at the end of the day + # if i == 0: + # external_forces_rod_one[..., i] -= 0.5 * net_contact_force + # external_forces_rod_one[..., i + 1] -= net_contact_force + # + # if j==0: + # external_forces_rod_two[..., j] += 0.5 * net_contact_force + # external_forces_rod_two[..., j + 1] += net_contact_force + # elif j== n_points_rod_two: + # external_forces_rod_two[..., j] += net_contact_force + # external_forces_rod_two[..., j + 1] += 0.5 * net_contact_force + # else: + # external_forces_rod_two[..., j] += 0.75 * net_contact_force + # external_forces_rod_two[..., j + 1] += 0.75 * net_contact_force + # + # elif i == n_points_rod_one: + # external_forces_rod_one[..., i] -= net_contact_force + # external_forces_rod_one[..., i + 1] -= 0.5 * net_contact_force + # + # if j == 0: + # external_forces_rod_two[..., j] += 0.5 * net_contact_force + # external_forces_rod_two[..., j + 1] += net_contact_force + # elif j == n_points_rod_two: + # external_forces_rod_two[..., j] += net_contact_force + # external_forces_rod_two[..., j + 1] += 0.5 * net_contact_force + # else: + # external_forces_rod_two[..., j] += 0.75 * net_contact_force + # external_forces_rod_two[..., j + 1] += 0.75 * net_contact_force + # else: + # external_forces_rod_one[..., i] -= net_contact_force + # external_forces_rod_one[..., i + 1] -= net_contact_force + # + # if j==0: + # external_forces_rod_two[..., j] += 0.5 * net_contact_force + # external_forces_rod_two[..., j + 1] += 1.5 * net_contact_force + # elif j== n_points_rod_two: + # external_forces_rod_two[..., j] += 1.5*net_contact_force + # external_forces_rod_two[..., j + 1] += 0.5 * net_contact_force + # else: + # external_forces_rod_two[..., j] += net_contact_force + # external_forces_rod_two[..., j + 1] += net_contact_force + + if i == 0: + external_forces_rod_one[..., i] -= net_contact_force * 2 / 3 + external_forces_rod_one[..., i + 1] -= net_contact_force * 4 / 3 + elif i == n_points_rod_one: + external_forces_rod_one[..., i] -= net_contact_force * 4 / 3 + external_forces_rod_one[..., i + 1] -= net_contact_force * 2 / 3 + else: + external_forces_rod_one[..., i] -= net_contact_force + external_forces_rod_one[..., i + 1] -= net_contact_force + + if j == 0: + external_forces_rod_two[..., j] += net_contact_force * 2 / 3 + external_forces_rod_two[..., j + 1] += net_contact_force * 4 / 3 + elif j == n_points_rod_two: + external_forces_rod_two[..., j] += net_contact_force * 4 / 3 + external_forces_rod_two[..., j + 1] += net_contact_force * 2 / 3 + else: + external_forces_rod_two[..., j] += net_contact_force + external_forces_rod_two[..., j + 1] += net_contact_force + + +@numba.njit(cache=True) +def _calculate_contact_forces_self_rod( + x_collection_rod, + radius_rod, + length_rod, + tangent_rod, + velocity_rod, + external_forces_rod, + contact_k, + contact_nu, +): + # We already pass in only the first n_elem x + n_points_rod = x_collection_rod.shape[1] + edge_collection_rod_one = _batch_product_k_ik_to_ik(length_rod, tangent_rod) + + for i in range(n_points_rod): + skip = 1 + np.ceil(0.8 * np.pi * radius_rod[i] / length_rod[i]) + for j in range(i - skip, -1, -1): + radii_sum = radius_rod[i] + radius_rod[j] + length_sum = length_rod[i] + length_rod[j] + # Element-wise bounding box + x_selected_rod_index_i = x_collection_rod[..., i] + x_selected_rod_index_j = x_collection_rod[..., j] + + del_x = x_selected_rod_index_i - x_selected_rod_index_j + norm_del_x = _norm(del_x) + + # If outside then don't process + if norm_del_x >= (radii_sum + length_sum): + continue + + # find the shortest line segment between the two centerline + # segments : differs from normal cylinder-cylinder intersection + distance_vector = _find_min_dist( + x_selected_rod_index_i, + edge_collection_rod_one[..., i], + x_selected_rod_index_j, + edge_collection_rod_one[..., j], + ) + distance_vector_length = _norm(distance_vector) + distance_vector /= distance_vector_length + + gamma = radii_sum - distance_vector_length + + # If distance is large, don't worry about it + if gamma < -1e-5: + continue + + # CHECK FOR GAMMA > 0.0, heaviside but we need to overload it in numba + # As a quick fix, use this instead + mask = (gamma > 0.0) * 1.0 + + contact_force = contact_k * gamma + interpenetration_velocity = 0.5 * ( + (velocity_rod[..., i] + velocity_rod[..., i + 1]) + - (velocity_rod[..., j] + velocity_rod[..., j + 1]) + ) + contact_damping_force = contact_nu * _dot_product( + interpenetration_velocity, distance_vector + ) + + # magnitude* direction + net_contact_force = ( + 0.5 * mask * (contact_damping_force + contact_force) + ) * distance_vector + + # Add it to the rods at the end of the day + # if i == 0: + # external_forces_rod[...,i] -= net_contact_force *2/3 + # external_forces_rod[...,i+1] -= net_contact_force * 4/3 + if i == n_points_rod: + external_forces_rod[..., i] -= net_contact_force * 4 / 3 + external_forces_rod[..., i + 1] -= net_contact_force * 2 / 3 + else: + external_forces_rod[..., i] -= net_contact_force + external_forces_rod[..., i + 1] -= net_contact_force + + if j == 0: + external_forces_rod[..., j] += net_contact_force * 2 / 3 + external_forces_rod[..., j + 1] += net_contact_force * 4 / 3 + # elif j == n_points_rod: + # external_forces_rod[..., j] += net_contact_force * 4/3 + # external_forces_rod[..., j+1] += net_contact_force * 2/3 + else: + external_forces_rod[..., j] += net_contact_force + external_forces_rod[..., j + 1] += net_contact_force + + @numba.njit(cache=True) def _aabbs_not_intersecting(aabb_one, aabb_two): """ Returns true if not intersecting else false""" @@ -198,8 +445,8 @@ def _aabbs_not_intersecting(aabb_one, aabb_two): @numba.njit(cache=True) -def _prune_using_aabbs( - rod_one_position, +def _prune_using_aabbs_rod_rigid_body( + rod_one_position_collection, rod_one_radius_collection, rod_one_length_collection, cylinder_position, @@ -214,8 +461,12 @@ def _prune_using_aabbs( rod_one_length_collection ) for i in range(3): - aabb_rod[i, 0] = np.min(rod_one_position[i]) - max_possible_dimension[i] - aabb_rod[i, 1] = np.max(rod_one_position[i]) + max_possible_dimension[i] + aabb_rod[i, 0] = ( + np.min(rod_one_position_collection[i]) - max_possible_dimension[i] + ) + aabb_rod[i, 1] = ( + np.max(rod_one_position_collection[i]) + max_possible_dimension[i] + ) # Is actually Q^T * d but numba complains about performance so we do # d^T @ Q @@ -235,6 +486,44 @@ def _prune_using_aabbs( return _aabbs_not_intersecting(aabb_cylinder, aabb_rod) +@numba.njit(cache=True) +def _prune_using_aabbs_rod_rod( + rod_one_position_collection, + rod_one_radius_collection, + rod_one_length_collection, + rod_two_position_collection, + rod_two_radius_collection, + rod_two_length_collection, +): + max_possible_dimension = np.zeros((3,)) + aabb_rod_one = np.empty((3, 2)) + aabb_rod_two = np.empty((3, 2)) + max_possible_dimension[...] = np.max(rod_one_radius_collection) + np.max( + rod_one_length_collection + ) + for i in range(3): + aabb_rod_one[i, 0] = ( + np.min(rod_one_position_collection[i]) - max_possible_dimension[i] + ) + aabb_rod_one[i, 1] = ( + np.max(rod_one_position_collection[i]) + max_possible_dimension[i] + ) + + max_possible_dimension[...] = np.max(rod_two_radius_collection) + np.max( + rod_two_length_collection + ) + + for i in range(3): + aabb_rod_two[i, 0] = ( + np.min(rod_two_position_collection[i]) - max_possible_dimension[i] + ) + aabb_rod_two[i, 1] = ( + np.max(rod_two_position_collection[i]) + max_possible_dimension[i] + ) + + return _aabbs_not_intersecting(aabb_rod_two, aabb_rod_one) + + class ExternalContact(FreeJoint): """ Assumes that the second entity is a rigid body for now, can be @@ -253,39 +542,116 @@ class ExternalContact(FreeJoint): def __init__(self, k, nu): super().__init__(k, nu) - def apply_forces(self, rod_one, index_one, cylinder_two, index_two): + def apply_forces(self, rod_one, index_one, rod_two, index_two): + # del index_one, index_two + + # TODO: raise error during the initialization if rod one is rigid body. + + # If rod two has one element, then it is rigid body. + if rod_two.n_elems == 1: + cylinder_two = rod_two + # First, check for a global AABB bounding box, and see whether that + # intersects + if _prune_using_aabbs_rod_rigid_body( + rod_one.position_collection, + rod_one.radius, + rod_one.lengths, + cylinder_two.position_collection, + cylinder_two.director_collection, + cylinder_two.radius, + cylinder_two.length, + ): + return + + x_cyl = ( + cylinder_two.position_collection[..., 0] + - 0.5 * cylinder_two.length * cylinder_two.director_collection[2, :, 0] + ) + + _calculate_contact_forces_rod_rigid_body( + rod_one.position_collection[..., :-1], + rod_one.lengths * rod_one.tangents, + x_cyl, + cylinder_two.length * cylinder_two.director_collection[2, :, 0], + rod_one.radius + cylinder_two.radius, + rod_one.lengths + cylinder_two.length, + rod_one.internal_forces, + rod_one.external_forces, + cylinder_two.external_forces, + rod_one.velocity_collection, + cylinder_two.velocity_collection, + self.k, + self.nu, + ) + + else: + # First, check for a global AABB bounding box, and see whether that + # intersects + + if _prune_using_aabbs_rod_rod( + rod_one.position_collection, + rod_one.radius, + rod_one.lengths, + rod_two.position_collection, + rod_two.radius, + rod_two.lengths, + ): + return + + _calculate_contact_forces_rod_rod( + rod_one.position_collection[ + ..., :-1 + ], # Discount last node, we want element start position + rod_one.radius, + rod_one.lengths, + rod_one.tangents, + rod_one.velocity_collection, + rod_one.internal_forces, + rod_one.external_forces, + rod_two.position_collection[ + ..., :-1 + ], # Discount last node, we want element start position + rod_two.radius, + rod_two.lengths, + rod_two.tangents, + rod_two.velocity_collection, + rod_two.internal_forces, + rod_two.external_forces, + self.k, + self.nu, + ) + + +class SelfContact(FreeJoint): + """ + Assumes that the second entity is a rigid body for now, can be + changed at a later time + + Most of the cylinder-cylinder contact SHOULD be implemented + as given in this paper: + http://larochelle.sdsmt.edu/publications/2005-2009/Collision%20Detection%20of%20Cylindrical%20Rigid%20Bodies%20Using%20Line%20Geometry.pdf + + but, it isn't (the elastica-cpp kernels are implented)! + 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): + super().__init__(k, nu) + + def apply_forces(self, rod_one, index_one, rod_two, index_two): # del index_one, index_two - # First, check for a global AABB bounding box, and see whether that - # intersects - if _prune_using_aabbs( - rod_one.position_collection, + _calculate_contact_forces_self_rod( + rod_one.position_collection[ + ..., :-1 + ], # Discount last node, we want element start position rod_one.radius, rod_one.lengths, - cylinder_two.position_collection, - cylinder_two.director_collection, - cylinder_two.radius, - cylinder_two.length, - ): - return - - x_cyl = ( - cylinder_two.position_collection[..., 0] - - 0.5 * cylinder_two.length * cylinder_two.director_collection[2, :, 0] - ) - - _calculate_contact_forces( - rod_one.position_collection[..., :-1], - rod_one.lengths * rod_one.tangents, - x_cyl, - cylinder_two.length * cylinder_two.director_collection[2, :, 0], - rod_one.radius + cylinder_two.radius, - rod_one.lengths + cylinder_two.length, - rod_one.internal_forces, - rod_one.external_forces, - cylinder_two.external_forces, + rod_one.tangents, rod_one.velocity_collection, - cylinder_two.velocity_collection, + rod_one.external_forces, self.k, self.nu, ) diff --git a/elastica/joint.py b/elastica/joint.py index 32661ffc..a7b4f885 100644 --- a/elastica/joint.py +++ b/elastica/joint.py @@ -1,5 +1,5 @@ __doc__ = """ Module containing joint classes to connect multiple rods together. """ -__all__ = ["FreeJoint", "HingeJoint", "FixedJoint", "ExternalContact"] +__all__ = ["FreeJoint", "HingeJoint", "FixedJoint", "ExternalContact", "SelfContact"] import numpy as np from elastica.utils import Tolerance, MaxDimension from elastica import IMPORT_NUMBA @@ -885,6 +885,6 @@ def apply_torques(self, rod_one, index_one, rod_two, index_two): # ) if IMPORT_NUMBA: - from elastica._elastica_numba._joint import ExternalContact + from elastica._elastica_numba._joint import ExternalContact, SelfContact else: from elastica._elastica_numpy._joint import ExternalContact diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py new file mode 100644 index 00000000..204c1d0a --- /dev/null +++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py @@ -0,0 +1,164 @@ +import sys + +sys.path.append("../../../") +from elastica import * +from examples.RodContactCase.post_processing import ( + plot_video_with_surface, + plot_velocity, +) + + +class InclinedRodRodContact( + BaseSystemCollection, Constraints, Connections, Forcing, CallBacks +): + pass + + +inclined_rod_rod_contact_sim = InclinedRodRodContact() + +# Simulation parameters +dt = 5e-5 +final_time = 20 +total_steps = int(final_time / dt) +time_step = np.float64(final_time / total_steps) +rendering_fps = 20 +step_skip = int(1.0 / (rendering_fps * time_step)) + +# Rod parameters +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 + +# Rod orientations +start = np.zeros(3,) +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)]) + + +# Rod 1 +n_elem_rod_one = 50 +start_rod_one = start + normal * 0.2 + +rod_one = CosseratRod.straight_rod( + n_elem_rod_one, + start_rod_one, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + poisson_ratio, +) + +rod_one.velocity_collection[:] += 0.05 * -normal.reshape(3, 1) + +inclined_rod_rod_contact_sim.append(rod_one) + +# Rod 2 +n_elem_rod_two = 60 + +inclination = np.deg2rad(0) +direction = np.array([0.0, np.cos(inclination), np.sin(inclination)]) +normal = np.array([0.0, -np.sin(inclination), np.cos(inclination)]) +start_rod_two = start + +rod_two = CosseratRod.straight_rod( + n_elem_rod_two, + start_rod_two, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + poisson_ratio, +) + +inclined_rod_rod_contact_sim.append(rod_two) + +# Contact between two rods +inclined_rod_rod_contact_sim.connect(rod_one, rod_two).using( + ExternalContact, k=1e3, nu=0.0 +) + +# Add call backs +class RodCallBack(CallBackBaseClass): + """ + + """ + + def __init__(self, step_skip: int, callback_params: dict): + CallBackBaseClass.__init__(self) + self.every = step_skip + self.callback_params = callback_params + + def make_callback(self, system, time, current_step: int): + if current_step % self.every == 0: + self.callback_params["time"].append(time) + self.callback_params["step"].append(current_step) + self.callback_params["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()) + 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 + + +post_processing_dict_rod1 = defaultdict( + list +) # list which collected data will be append +# set the diagnostics for rod and collect data +inclined_rod_rod_contact_sim.collect_diagnostics(rod_one).using( + RodCallBack, step_skip=step_skip, callback_params=post_processing_dict_rod1, +) + +post_processing_dict_rod2 = defaultdict( + list +) # list which collected data will be append +# set the diagnostics for rod and collect data +inclined_rod_rod_contact_sim.collect_diagnostics(rod_two).using( + RodCallBack, step_skip=step_skip, callback_params=post_processing_dict_rod2, +) + +inclined_rod_rod_contact_sim.finalize() +# Do the simulation + +timestepper = PositionVerlet() +integrate(timestepper, inclined_rod_rod_contact_sim, final_time, total_steps) + +# plotting the videos +filename_video = "inclined_rods_contact.mp4" +plot_video_with_surface( + [post_processing_dict_rod1, post_processing_dict_rod2], + video_name=filename_video, + fps=rendering_fps, + step=1, + vis3D=True, + vis2D=True, +) + +filaname = "inclined_rods_velocity.png" +plot_velocity( + post_processing_dict_rod1, + post_processing_dict_rod2, + filename=filaname, + SAVE_FIGURE=True, +) diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py new file mode 100644 index 00000000..a9637a0c --- /dev/null +++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py @@ -0,0 +1,161 @@ +import sys + +sys.path.append("../../../") +from elastica import * +from examples.RodContactCase.post_processing import ( + plot_video_with_surface, + plot_velocity, +) + + +class ParallelRodRodContact( + BaseSystemCollection, Constraints, Connections, Forcing, CallBacks +): + pass + + +parallel_rod_rod_contact_sim = ParallelRodRodContact() + +# Simulation parameters +dt = 5e-5 +final_time = 10 +total_steps = int(final_time / dt) +time_step = np.float64(final_time / total_steps) +rendering_fps = 20 +step_skip = int(1.0 / (rendering_fps * time_step)) + +# Rod parameters +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 + +# Rod orientations +start = np.zeros(3,) +inclination = np.deg2rad(0) +direction = np.array([0.0, np.cos(inclination), np.sin(inclination)]) +normal = np.array([0.0, -np.sin(inclination), np.cos(inclination)]) + + +# Rod 1 +n_elem_rod_one = 50 +start_rod_one = start + normal * 0.2 + +rod_one = CosseratRod.straight_rod( + n_elem_rod_one, + start_rod_one, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + poisson_ratio, +) + +rod_one.velocity_collection[:] += 0.05 * -normal.reshape(3, 1) + +parallel_rod_rod_contact_sim.append(rod_one) + +# Rod 2 +n_elem_rod_two = 50 + +start_rod_two = start + +rod_two = CosseratRod.straight_rod( + n_elem_rod_two, + start_rod_two, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + poisson_ratio, +) + +parallel_rod_rod_contact_sim.append(rod_two) + +# Contact between two rods +parallel_rod_rod_contact_sim.connect(rod_one, rod_two).using( + ExternalContact, k=1e3, nu=0.001 +) + +# Add call backs +class RodCallBack(CallBackBaseClass): + """ + + """ + + def __init__(self, step_skip: int, callback_params: dict): + CallBackBaseClass.__init__(self) + self.every = step_skip + self.callback_params = callback_params + + def make_callback(self, system, time, current_step: int): + if current_step % self.every == 0: + self.callback_params["time"].append(time) + self.callback_params["step"].append(current_step) + self.callback_params["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()) + 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 + + +post_processing_dict_rod1 = defaultdict( + list +) # list which collected data will be append +# set the diagnostics for rod and collect data +parallel_rod_rod_contact_sim.collect_diagnostics(rod_one).using( + RodCallBack, step_skip=step_skip, callback_params=post_processing_dict_rod1, +) + +post_processing_dict_rod2 = defaultdict( + list +) # list which collected data will be append +# set the diagnostics for rod and collect data +parallel_rod_rod_contact_sim.collect_diagnostics(rod_two).using( + RodCallBack, step_skip=step_skip, callback_params=post_processing_dict_rod2, +) + +parallel_rod_rod_contact_sim.finalize() +# Do the simulation + +timestepper = PositionVerlet() +integrate(timestepper, parallel_rod_rod_contact_sim, final_time, total_steps) + +# plotting the videos +filename_video = "parallel_rods_contact.mp4" +plot_video_with_surface( + [post_processing_dict_rod1, post_processing_dict_rod2], + video_name=filename_video, + fps=rendering_fps, + step=1, + vis3D=True, + vis2D=True, +) + +filaname = "parallel_rods_velocity.png" +plot_velocity( + post_processing_dict_rod1, + post_processing_dict_rod2, + filename=filaname, + SAVE_FIGURE=True, +) diff --git a/examples/RodContactCase/RodSelfContact/rod_self_contact_validation.py b/examples/RodContactCase/RodSelfContact/rod_self_contact_validation.py new file mode 100644 index 00000000..8004ee22 --- /dev/null +++ b/examples/RodContactCase/RodSelfContact/rod_self_contact_validation.py @@ -0,0 +1,128 @@ +import sys + +sys.path.append("../../../") +from elastica import * +from examples.RodContactCase.post_processing import ( + plot_video_with_surface, + plot_velocity, +) + + +class RodSelfContact( + BaseSystemCollection, Constraints, Connections, Forcing, CallBacks +): + pass + + +rod_self_contact = RodSelfContact() + +# Simulation parameters +dt = 5e-5 +final_time = 10 +total_steps = int(final_time / dt) +time_step = np.float64(final_time / total_steps) +rendering_fps = 20 +step_skip = int(1.0 / (rendering_fps * time_step)) + +# Rod parameters +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 + +# Rod orientations +start = np.zeros(3,) +inclination = np.deg2rad(0) +direction = np.array([0.0, np.cos(inclination), np.sin(inclination)]) +normal = np.array([0.0, -np.sin(inclination), np.cos(inclination)]) + + +# Rod 1 +n_elem_rod_one = 50 +start_rod_one = start + normal * 0.2 + +rod_one = CosseratRod.straight_rod( + n_elem_rod_one, + start_rod_one, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + poisson_ratio, +) + +# rod_one.velocity_collection[:] += 0.05 * -normal.reshape(3,1) + +rod_self_contact.append(rod_one) + + +# Contact between two rods +rod_self_contact.connect(rod_one, rod_one).using(SelfContact, k=1e3, nu=0.001) + +# Add call backs +class RodCallBack(CallBackBaseClass): + """ + + """ + + def __init__(self, step_skip: int, callback_params: dict): + CallBackBaseClass.__init__(self) + self.every = step_skip + self.callback_params = callback_params + + def make_callback(self, system, time, current_step: int): + if current_step % self.every == 0: + self.callback_params["time"].append(time) + self.callback_params["step"].append(current_step) + self.callback_params["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()) + 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 + + +post_processing_dict_rod1 = defaultdict( + list +) # list which collected data will be append +# set the diagnostics for rod and collect data +rod_self_contact.collect_diagnostics(rod_one).using( + RodCallBack, step_skip=step_skip, callback_params=post_processing_dict_rod1, +) + + +rod_self_contact.finalize() +# Do the simulation + +timestepper = PositionVerlet() +integrate(timestepper, rod_self_contact, final_time, total_steps) + +# plotting the videos +filename_video = "self_contact.mp4" +plot_video_with_surface( + [post_processing_dict_rod1], + video_name=filename_video, + fps=rendering_fps, + step=1, + vis3D=True, + vis2D=True, +) + +# filaname = "parallel_rods_velocity.png" +# plot_velocity(post_processing_dict_rod1, post_processing_dict_rod2, filename=filaname, SAVE_FIGURE=True) diff --git a/examples/RodContactCase/post_processing.py b/examples/RodContactCase/post_processing.py new file mode 100644 index 00000000..8fde465b --- /dev/null +++ b/examples/RodContactCase/post_processing.py @@ -0,0 +1,366 @@ +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.colors import to_rgb +from mpl_toolkits.mplot3d import proj3d, Axes3D +from tqdm import tqdm + +from typing import Dict, Sequence + + +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}) + + # 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] + + # 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) + + rod_scatters[rod_idx] = ax.scatter( + inst_position[0], + inst_position[1], + inst_position[2], + s=np.pi * (scaling_factor * inst_radius) ** 2, + ) + + ax.set_aspect("equal") + video_name_3D = "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 + ) + 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 + ) + + 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): + 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) + 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, + ) + + ax.set_aspect("equal") + video_name_2D = "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 + ) + 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 + ) + + 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 + + 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) + 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, + ) + + ax.set_aspect("equal") + video_name_2D = "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 + ) + 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 + ) + + 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(-1.0, 1.0) + ax.set_ylim(-1.0, 1.0) + + 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) + 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, + ) + + ax.set_aspect("equal") + video_name_2D = "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 + ) + 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 + ) + + 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()) + + +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].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].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].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) + plt.savefig(filename) From 92e28412826b2db483f212b885b5b3a7d24cabbf Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Tue, 25 Aug 2020 11:16:24 -0500 Subject: [PATCH 02/14] enhencement:Plectone and Solenoid cases added --- .../PlectonemesCase/plectoneme_case.py | 206 ++++++++++++++++ .../SolenoidsCase/solenoid_case.py | 224 ++++++++++++++++++ 2 files changed, 430 insertions(+) create mode 100644 examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py create mode 100644 examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py diff --git a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py new file mode 100644 index 00000000..2c94316a --- /dev/null +++ b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py @@ -0,0 +1,206 @@ +import sys + +sys.path.append("../../../../") +from elastica import * +from examples.RodContactCase.post_processing import ( + plot_video_with_surface, + plot_velocity, +) + + +class PlectonemesCase( + BaseSystemCollection, Constraints, Connections, Forcing, CallBacks +): + pass + + +plectonemes_sim = PlectonemesCase() + +# Simulation parameters +time_twist = 75 +time_relax = 50 +final_time = time_relax + time_twist + +base_length = 1.0 +n_elem = 150 + +dt = 0.0025 * base_length / n_elem # 1E-2 +total_steps = int(final_time / dt) +time_step = np.float64(final_time / total_steps) +rendering_fps = 20 +step_skip = int(1.0 / (rendering_fps * time_step)) + +# Rest of the rod parameters and construct rod +base_radius = 0.0125 +base_area = np.pi * base_radius ** 2 +volume = base_area * base_length +mass = 1.0 +density = mass / volume +nu = 2.0 +E = 1e6 +poisson_ratio = 0.5 + +direction = np.array([0.0, 1.0, 0]) +normal = np.array([0.0, 0.0, 1.0]) +start = np.zeros(3,) + +sherable_rod = CosseratRod.straight_rod( + n_elem, + start, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + poisson_ratio, +) + +sherable_rod.velocity_collection[..., int(n_elem / 2)] -= 1e-4 + + +plectonemes_sim.append(sherable_rod) + +# boundary condition +from elastica._rotations import _get_rotation_matrix + + +class SelonoidsBC(FreeRod): + """ + + """ + + def __init__( + self, + position_start, + position_end, + director_start, + director_end, + twisting_time, + number_of_rotations, + ): + FreeRod.__init__(self) + self.twisting_time = twisting_time + + theta = 2.0 * number_of_rotations * np.pi + + angel_vel_scalar = theta / self.twisting_time + + direction = -(position_end - position_start) / np.linalg.norm( + position_end - position_start + ) + + axis_of_rotation_in_material_frame = director_end @ direction + axis_of_rotation_in_material_frame /= np.linalg.norm( + axis_of_rotation_in_material_frame + ) + + self.final_end_directors = ( + _get_rotation_matrix( + -theta, axis_of_rotation_in_material_frame.reshape(3, 1) + ).reshape(3, 3) + @ director_end + ) # rotation_matrix wants vectors 3,1 + + self.ang_vel = angel_vel_scalar * axis_of_rotation_in_material_frame + + self.position_start = position_start + self.director_start = director_start + + def constrain_values(self, rod, time): + if time > self.twisting_time: + rod.position_collection[..., 0] = self.position_start + rod.position_collection[0, -1] = 0.0 + rod.position_collection[2, -1] = 0.0 + + rod.director_collection[..., 0] = self.director_start + rod.director_collection[..., -1] = self.final_end_directors + + def constrain_rates(self, rod, time): + if time > self.twisting_time: + rod.velocity_collection[..., 0] = 0.0 + rod.omega_collection[..., 0] = 0.0 + + rod.velocity_collection[..., -1] = 0.0 + rod.omega_collection[..., -1] = 0.0 + + else: + rod.velocity_collection[..., 0] = 0.0 + rod.omega_collection[..., 0] = 0.0 + + rod.velocity_collection[0, -1] = 0.0 + rod.velocity_collection[2, -1] = 0.0 + rod.omega_collection[..., -1] = -self.ang_vel + + +plectonemes_sim.constrain(sherable_rod).using( + SelonoidsBC, + constrained_position_idx=(0, -1), + constrained_director_idx=(0, -1), + twisting_time=time_twist, + number_of_rotations=10, +) + +# Add self contact to prevent penetration +plectonemes_sim.connect(sherable_rod, sherable_rod).using(SelfContact, k=1e4, nu=10) + +# Add callback functions for plotting position of the rod later on +class RodCallBack(CallBackBaseClass): + """ + + """ + + def __init__(self, step_skip: int, callback_params: dict): + CallBackBaseClass.__init__(self) + self.every = step_skip + self.callback_params = callback_params + + def make_callback(self, system, time, current_step: int): + if current_step % self.every == 0: + self.callback_params["time"].append(time) + self.callback_params["step"].append(current_step) + self.callback_params["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()) + 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 + + +post_processing_dict = defaultdict(list) # list which collected data will be append +# set the diagnostics for rod and collect data +plectonemes_sim.collect_diagnostics(sherable_rod).using( + RodCallBack, step_skip=step_skip, callback_params=post_processing_dict, +) + +# finalize simulation +plectonemes_sim.finalize() + +# Run the simulation +time_stepper = PositionVerlet() +integrate(time_stepper, plectonemes_sim, final_time, total_steps) + +# plotting the videos +filename_video = "plectonemes.mp4" +plot_video_with_surface( + [post_processing_dict], + video_name=filename_video, + fps=rendering_fps, + step=1, + vis3D=True, + vis2D=True, + x_limits=[-0.5, 0.5], + y_limits=[-0.1, 1.1], + z_limits=[-0.5, 0.5], +) diff --git a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py new file mode 100644 index 00000000..ef6cab6a --- /dev/null +++ b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py @@ -0,0 +1,224 @@ +import sys + +sys.path.append("../../../../") +from elastica import * +from examples.RodContactCase.post_processing import ( + plot_video_with_surface, + plot_velocity, +) + + +class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, CallBacks): + pass + + +solenoid_sim = SolenoidCase() + +# Simulation parameters +number_of_rotations = 13 +time_twist = 5 * number_of_rotations +time_relax = 50 +final_time = 30 # time_relax + time_twist + +base_length = 1.0 +n_elem = 100 + +dt = 0.0025 * base_length / n_elem # 1E-2 +total_steps = int(final_time / dt) +time_step = np.float64(final_time / total_steps) +rendering_fps = 20 +step_skip = int(1.0 / (rendering_fps * time_step)) + +# Rest of the rod parameters and construct rod +base_radius = 0.0125 +base_area = np.pi * base_radius ** 2 +volume = base_area * base_length +mass = 1.0 +density = mass / volume +nu = 2.0 +E = 1e6 +poisson_ratio = 0.5 + +direction = np.array([0.0, 1.0, 0]) +normal = np.array([0.0, 0.0, 1.0]) +start = np.zeros(3,) + +F_pulling_scalar = 300 + +# fractional compression of rod needed to resist force F +c = F_pulling_scalar / (base_area * E) + +end = start + direction * (1 - c) * base_length + +rod_initial_positions = np.zeros((3, n_elem + 1)) +for i in range(3): + rod_initial_positions[i, ...] = np.linspace(start[i], end[i], n_elem + 1) + +sherable_rod = CosseratRod.straight_rod( + n_elem, + start, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + poisson_ratio, + position=rod_initial_positions, +) + +sherable_rod.rest_kappa[:] = sherable_rod.kappa +sherable_rod.rest_sigma[:] = sherable_rod.sigma + +sherable_rod.velocity_collection[..., int(n_elem / 2)] -= 1e-4 + + +solenoid_sim.append(sherable_rod) + +# boundary condition +from elastica._rotations import _get_rotation_matrix + + +class SelonoidsBC(FreeRod): + """ + + """ + + def __init__( + self, + position_start, + position_end, + director_start, + director_end, + twisting_time, + number_of_rotations, + ): + FreeRod.__init__(self) + self.twisting_time = twisting_time + + theta = 2.0 * number_of_rotations * np.pi + + angel_vel_scalar = theta / self.twisting_time + + direction = -(position_end - position_start) / np.linalg.norm( + position_end - position_start + ) + + axis_of_rotation_in_material_frame = director_end @ direction + axis_of_rotation_in_material_frame /= np.linalg.norm( + axis_of_rotation_in_material_frame + ) + + self.final_end_directors = ( + _get_rotation_matrix( + -theta, axis_of_rotation_in_material_frame.reshape(3, 1) + ).reshape(3, 3) + @ director_end + ) # rotation_matrix wants vectors 3,1 + + self.ang_vel = angel_vel_scalar * axis_of_rotation_in_material_frame + + self.position_start = position_start + self.director_start = director_start + + def constrain_values(self, rod, time): + if time > self.twisting_time: + rod.position_collection[..., 0] = self.position_start + rod.position_collection[0, -1] = 0.0 + rod.position_collection[2, -1] = 0.0 + + rod.director_collection[..., 0] = self.director_start + rod.director_collection[..., -1] = self.final_end_directors + + def constrain_rates(self, rod, time): + if time > self.twisting_time: + rod.velocity_collection[..., 0] = 0.0 + rod.omega_collection[..., 0] = 0.0 + + rod.velocity_collection[..., -1] = 0.0 + rod.omega_collection[..., -1] = 0.0 + + else: + rod.velocity_collection[..., 0] = 0.0 + rod.omega_collection[..., 0] = 0.0 + + rod.velocity_collection[0, -1] = 0.0 + rod.velocity_collection[2, -1] = 0.0 + rod.omega_collection[..., -1] = -self.ang_vel + + +solenoid_sim.constrain(sherable_rod).using( + SelonoidsBC, + constrained_position_idx=(0, -1), + constrained_director_idx=(0, -1), + twisting_time=time_twist, + number_of_rotations=number_of_rotations, +) + +solenoid_sim.add_forcing_to(sherable_rod).using( + EndpointForces, np.zeros(3,), F_pulling_scalar * direction, ramp_up_time=time_twist +) + +# Add self contact to prevent penetration +# solenoid_sim.connect(sherable_rod, sherable_rod).using(SelfContact, k=1E4, nu=10) + +# Add callback functions for plotting position of the rod later on +class RodCallBack(CallBackBaseClass): + """ + + """ + + def __init__(self, step_skip: int, callback_params: dict): + CallBackBaseClass.__init__(self) + self.every = step_skip + self.callback_params = callback_params + + def make_callback(self, system, time, current_step: int): + if current_step % self.every == 0: + self.callback_params["time"].append(time) + self.callback_params["step"].append(current_step) + self.callback_params["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()) + 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 + + +post_processing_dict = defaultdict(list) # list which collected data will be append +# set the diagnostics for rod and collect data +solenoid_sim.collect_diagnostics(sherable_rod).using( + RodCallBack, step_skip=step_skip, callback_params=post_processing_dict, +) + +# finalize simulation +solenoid_sim.finalize() + +# Run the simulation +time_stepper = PositionVerlet() +integrate(time_stepper, solenoid_sim, final_time, total_steps) + +# plotting the videos +filename_video = "solenoid.mp4" +plot_video_with_surface( + [post_processing_dict], + video_name=filename_video, + fps=rendering_fps, + step=1, + vis3D=True, + vis2D=True, + x_limits=[-0.5, 0.5], + y_limits=[-0.1, 1.1], + z_limits=[-0.5, 0.5], +) From 5425d07dc1c07fca5e744145346d879e8e82ef6f Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Wed, 26 Aug 2020 11:01:31 -0500 Subject: [PATCH 03/14] work in progress for plectoneme and solenoid --- .../PlectonemesCase/plectoneme_case.py | 8 ++--- .../SolenoidsCase/solenoid_case.py | 35 ++++++++++++------- 2 files changed, 27 insertions(+), 16 deletions(-) diff --git a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py index 2c94316a..a6222cf1 100644 --- a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py +++ b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py @@ -17,9 +17,9 @@ class PlectonemesCase( plectonemes_sim = PlectonemesCase() # Simulation parameters -time_twist = 75 +time_twist = 50 time_relax = 50 -final_time = time_relax + time_twist +final_time = 30#time_relax + time_twist base_length = 1.0 n_elem = 150 @@ -57,7 +57,7 @@ class PlectonemesCase( poisson_ratio, ) -sherable_rod.velocity_collection[..., int(n_elem / 2)] -= 1e-4 +sherable_rod.velocity_collection[2, int(n_elem / 2)] -= 1e-4 plectonemes_sim.append(sherable_rod) @@ -143,7 +143,7 @@ def constrain_rates(self, rod, time): ) # Add self contact to prevent penetration -plectonemes_sim.connect(sherable_rod, sherable_rod).using(SelfContact, k=1e4, nu=10) +# plectonemes_sim.connect(sherable_rod, sherable_rod).using(SelfContact, k=1e4, nu=10) # Add callback functions for plotting position of the rod later on class RodCallBack(CallBackBaseClass): diff --git a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py index ef6cab6a..ac9ecd79 100644 --- a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py +++ b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py @@ -16,22 +16,24 @@ class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, Call # Simulation parameters number_of_rotations = 13 +time_start_twist = 10 time_twist = 5 * number_of_rotations time_relax = 50 -final_time = 30 # time_relax + time_twist +final_time = time_relax + time_twist + time_start_twist base_length = 1.0 n_elem = 100 -dt = 0.0025 * base_length / n_elem # 1E-2 +dt = 0.008 * base_length / n_elem # 1E-2 total_steps = int(final_time / dt) time_step = np.float64(final_time / total_steps) rendering_fps = 20 step_skip = int(1.0 / (rendering_fps * time_step)) # Rest of the rod parameters and construct rod -base_radius = 0.0125 +base_radius = 0.025 base_area = np.pi * base_radius ** 2 +I = np.pi/4* base_radius**4 volume = base_area * base_length mass = 1.0 density = mass / volume @@ -65,13 +67,13 @@ class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, Call nu, E, poisson_ratio, - position=rod_initial_positions, + # position=rod_initial_positions, ) -sherable_rod.rest_kappa[:] = sherable_rod.kappa -sherable_rod.rest_sigma[:] = sherable_rod.sigma +# sherable_rod.rest_kappa[:] = sherable_rod.kappa +# sherable_rod.rest_sigma[:] = sherable_rod.sigma -sherable_rod.velocity_collection[..., int(n_elem / 2)] -= 1e-4 +# sherable_rod.velocity_collection[..., int(n_elem / 2)] -= 1e-3 solenoid_sim.append(sherable_rod) @@ -92,10 +94,12 @@ def __init__( director_start, director_end, twisting_time, + time_twis_start, number_of_rotations, ): FreeRod.__init__(self) self.twisting_time = twisting_time + self.time_twis_start = time_twis_start theta = 2.0 * number_of_rotations * np.pi @@ -123,7 +127,7 @@ def __init__( self.director_start = director_start def constrain_values(self, rod, time): - if time > self.twisting_time: + if time > self.twisting_time+self.time_twis_start: rod.position_collection[..., 0] = self.position_start rod.position_collection[0, -1] = 0.0 rod.position_collection[2, -1] = 0.0 @@ -132,13 +136,17 @@ def constrain_values(self, rod, time): rod.director_collection[..., -1] = self.final_end_directors def constrain_rates(self, rod, time): - if time > self.twisting_time: + if time > self.twisting_time+self.time_twis_start: rod.velocity_collection[..., 0] = 0.0 rod.omega_collection[..., 0] = 0.0 rod.velocity_collection[..., -1] = 0.0 rod.omega_collection[..., -1] = 0.0 + elif time Date: Fri, 28 Aug 2020 08:51:49 -0500 Subject: [PATCH 04/14] workin progress solenoid --- .../RodSelfContact/SolenoidsCase/solenoid_case.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py index ac9ecd79..eba7dca6 100644 --- a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py +++ b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py @@ -155,7 +155,7 @@ def constrain_rates(self, rod, time): rod.velocity_collection[2, -1] = 0.0 rod.omega_collection[..., -1] = -self.ang_vel - rod.velocity_collection[2,:] -= 1E-4 + rod.velocity_collection[2,int(rod.n_elems/2)] -= 1E-4 solenoid_sim.constrain(sherable_rod).using( @@ -168,11 +168,11 @@ def constrain_rates(self, rod, time): ) solenoid_sim.add_forcing_to(sherable_rod).using( - EndpointForces, np.zeros(3,), F_pulling_scalar * direction, ramp_up_time=time_start_twist + EndpointForces, np.zeros(3,), F_pulling_scalar * direction, ramp_up_time=time_start_twist-1 ) # Add self contact to prevent penetration -# solenoid_sim.connect(sherable_rod, sherable_rod).using(SelfContact, k=1E4, nu=10) +solenoid_sim.connect(sherable_rod, sherable_rod).using(SelfContact, k=1E4, nu=10) # Add callback functions for plotting position of the rod later on class RodCallBack(CallBackBaseClass): From b86ed8fc4f3e5287774347f47e2c749ad803b5ef Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 29 Aug 2020 09:12:10 -0500 Subject: [PATCH 05/14] validationc cases: self-contact validated This commit finizalize self-contact validation for numba and it adds validation cases, plectonome and solenoid. --- .../PlectonemesCase/plectoneme_case.py | 25 +++++++---- .../SolenoidsCase/solenoid_case.py | 41 +++++++------------ 2 files changed, 33 insertions(+), 33 deletions(-) diff --git a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py index a6222cf1..dbea5705 100644 --- a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py +++ b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py @@ -17,9 +17,11 @@ class PlectonemesCase( plectonemes_sim = PlectonemesCase() # Simulation parameters -time_twist = 50 +number_of_rotations = 10 +time_start_twist = 0 +time_twist = 5 * number_of_rotations time_relax = 50 -final_time = 30#time_relax + time_twist +final_time = time_relax + time_twist + time_start_twist base_length = 1.0 n_elem = 150 @@ -31,7 +33,7 @@ class PlectonemesCase( step_skip = int(1.0 / (rendering_fps * time_step)) # Rest of the rod parameters and construct rod -base_radius = 0.0125 +base_radius = 0.025 base_area = np.pi * base_radius ** 2 volume = base_area * base_length mass = 1.0 @@ -78,10 +80,12 @@ def __init__( director_start, director_end, twisting_time, + time_twis_start, number_of_rotations, ): FreeRod.__init__(self) self.twisting_time = twisting_time + self.time_twis_start = time_twis_start theta = 2.0 * number_of_rotations * np.pi @@ -109,7 +113,7 @@ def __init__( self.director_start = director_start def constrain_values(self, rod, time): - if time > self.twisting_time: + if time > self.twisting_time + self.time_twis_start: rod.position_collection[..., 0] = self.position_start rod.position_collection[0, -1] = 0.0 rod.position_collection[2, -1] = 0.0 @@ -118,13 +122,17 @@ def constrain_values(self, rod, time): rod.director_collection[..., -1] = self.final_end_directors def constrain_rates(self, rod, time): - if time > self.twisting_time: + if time > self.twisting_time + self.time_twis_start: rod.velocity_collection[..., 0] = 0.0 rod.omega_collection[..., 0] = 0.0 rod.velocity_collection[..., -1] = 0.0 rod.omega_collection[..., -1] = 0.0 + elif time < self.time_twis_start: + rod.velocity_collection[..., 0] = 0.0 + rod.omega_collection[..., 0] = 0.0 + else: rod.velocity_collection[..., 0] = 0.0 rod.omega_collection[..., 0] = 0.0 @@ -133,17 +141,20 @@ def constrain_rates(self, rod, time): rod.velocity_collection[2, -1] = 0.0 rod.omega_collection[..., -1] = -self.ang_vel + rod.velocity_collection[2, int(rod.n_elems / 2)] -= 1e-4 + plectonemes_sim.constrain(sherable_rod).using( SelonoidsBC, constrained_position_idx=(0, -1), constrained_director_idx=(0, -1), + time_twis_start=time_start_twist, twisting_time=time_twist, - number_of_rotations=10, + number_of_rotations=number_of_rotations, ) # Add self contact to prevent penetration -# plectonemes_sim.connect(sherable_rod, sherable_rod).using(SelfContact, k=1e4, nu=10) +plectonemes_sim.connect(sherable_rod, sherable_rod).using(SelfContact, k=1e4, nu=10) # Add callback functions for plotting position of the rod later on class RodCallBack(CallBackBaseClass): diff --git a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py index eba7dca6..c10cc230 100644 --- a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py +++ b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py @@ -19,12 +19,12 @@ class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, Call time_start_twist = 10 time_twist = 5 * number_of_rotations time_relax = 50 -final_time = time_relax + time_twist + time_start_twist +final_time = time_relax + time_twist + time_start_twist base_length = 1.0 n_elem = 100 -dt = 0.008 * base_length / n_elem # 1E-2 +dt = 0.0025 * base_length / n_elem # 1E-2 total_steps = int(final_time / dt) time_step = np.float64(final_time / total_steps) rendering_fps = 20 @@ -33,7 +33,7 @@ class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, Call # Rest of the rod parameters and construct rod base_radius = 0.025 base_area = np.pi * base_radius ** 2 -I = np.pi/4* base_radius**4 +I = np.pi / 4 * base_radius ** 4 volume = base_area * base_length mass = 1.0 density = mass / volume @@ -47,14 +47,6 @@ class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, Call F_pulling_scalar = 300 -# fractional compression of rod needed to resist force F -c = F_pulling_scalar / (base_area * E) - -end = start + direction * (1 - c) * base_length - -rod_initial_positions = np.zeros((3, n_elem + 1)) -for i in range(3): - rod_initial_positions[i, ...] = np.linspace(start[i], end[i], n_elem + 1) sherable_rod = CosseratRod.straight_rod( n_elem, @@ -67,14 +59,8 @@ class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, Call nu, E, poisson_ratio, - # position=rod_initial_positions, ) -# sherable_rod.rest_kappa[:] = sherable_rod.kappa -# sherable_rod.rest_sigma[:] = sherable_rod.sigma - -# sherable_rod.velocity_collection[..., int(n_elem / 2)] -= 1e-3 - solenoid_sim.append(sherable_rod) @@ -94,7 +80,7 @@ def __init__( director_start, director_end, twisting_time, - time_twis_start, + time_twis_start, number_of_rotations, ): FreeRod.__init__(self) @@ -127,7 +113,7 @@ def __init__( self.director_start = director_start def constrain_values(self, rod, time): - if time > self.twisting_time+self.time_twis_start: + if time > self.twisting_time + self.time_twis_start: rod.position_collection[..., 0] = self.position_start rod.position_collection[0, -1] = 0.0 rod.position_collection[2, -1] = 0.0 @@ -136,14 +122,14 @@ def constrain_values(self, rod, time): rod.director_collection[..., -1] = self.final_end_directors def constrain_rates(self, rod, time): - if time > self.twisting_time+self.time_twis_start: + if time > self.twisting_time + self.time_twis_start: rod.velocity_collection[..., 0] = 0.0 rod.omega_collection[..., 0] = 0.0 rod.velocity_collection[..., -1] = 0.0 rod.omega_collection[..., -1] = 0.0 - elif time Date: Mon, 31 Aug 2020 23:15:36 -0500 Subject: [PATCH 06/14] remove unused file example file --- .../rod_self_contact_validation.py | 128 ------------------ 1 file changed, 128 deletions(-) delete mode 100644 examples/RodContactCase/RodSelfContact/rod_self_contact_validation.py diff --git a/examples/RodContactCase/RodSelfContact/rod_self_contact_validation.py b/examples/RodContactCase/RodSelfContact/rod_self_contact_validation.py deleted file mode 100644 index 8004ee22..00000000 --- a/examples/RodContactCase/RodSelfContact/rod_self_contact_validation.py +++ /dev/null @@ -1,128 +0,0 @@ -import sys - -sys.path.append("../../../") -from elastica import * -from examples.RodContactCase.post_processing import ( - plot_video_with_surface, - plot_velocity, -) - - -class RodSelfContact( - BaseSystemCollection, Constraints, Connections, Forcing, CallBacks -): - pass - - -rod_self_contact = RodSelfContact() - -# Simulation parameters -dt = 5e-5 -final_time = 10 -total_steps = int(final_time / dt) -time_step = np.float64(final_time / total_steps) -rendering_fps = 20 -step_skip = int(1.0 / (rendering_fps * time_step)) - -# Rod parameters -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 - -# Rod orientations -start = np.zeros(3,) -inclination = np.deg2rad(0) -direction = np.array([0.0, np.cos(inclination), np.sin(inclination)]) -normal = np.array([0.0, -np.sin(inclination), np.cos(inclination)]) - - -# Rod 1 -n_elem_rod_one = 50 -start_rod_one = start + normal * 0.2 - -rod_one = CosseratRod.straight_rod( - n_elem_rod_one, - start_rod_one, - direction, - normal, - base_length, - base_radius, - density, - nu, - E, - poisson_ratio, -) - -# rod_one.velocity_collection[:] += 0.05 * -normal.reshape(3,1) - -rod_self_contact.append(rod_one) - - -# Contact between two rods -rod_self_contact.connect(rod_one, rod_one).using(SelfContact, k=1e3, nu=0.001) - -# Add call backs -class RodCallBack(CallBackBaseClass): - """ - - """ - - def __init__(self, step_skip: int, callback_params: dict): - CallBackBaseClass.__init__(self) - self.every = step_skip - self.callback_params = callback_params - - def make_callback(self, system, time, current_step: int): - if current_step % self.every == 0: - self.callback_params["time"].append(time) - self.callback_params["step"].append(current_step) - self.callback_params["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()) - 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 - - -post_processing_dict_rod1 = defaultdict( - list -) # list which collected data will be append -# set the diagnostics for rod and collect data -rod_self_contact.collect_diagnostics(rod_one).using( - RodCallBack, step_skip=step_skip, callback_params=post_processing_dict_rod1, -) - - -rod_self_contact.finalize() -# Do the simulation - -timestepper = PositionVerlet() -integrate(timestepper, rod_self_contact, final_time, total_steps) - -# plotting the videos -filename_video = "self_contact.mp4" -plot_video_with_surface( - [post_processing_dict_rod1], - video_name=filename_video, - fps=rendering_fps, - step=1, - vis3D=True, - vis2D=True, -) - -# filaname = "parallel_rods_velocity.png" -# plot_velocity(post_processing_dict_rod1, post_processing_dict_rod2, filename=filaname, SAVE_FIGURE=True) From 558219d81cfbc4d8c6cebfb49bd0ae6b840623be Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Wed, 23 Jun 2021 15:13:21 -0500 Subject: [PATCH 07/14] enhancement: solenoids case collect data for topology analysis This comment adds a new code to save rod position, radius, directors for topology analysis. --- .../SolenoidsCase/solenoid_case.py | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py index c10cc230..8cdf0d08 100644 --- a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py +++ b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py @@ -192,6 +192,7 @@ def make_callback(self, system, time, current_step: int): + system.compute_shear_energy() ) self.callback_params["total_energy"].append(total_energy) + self.callback_params["directors"].append(system.director_collection.copy()) return @@ -222,3 +223,23 @@ def make_callback(self, system, time, current_step: int): y_limits=[-0.1, 1.4], z_limits=[-0.5, 0.5], ) + +# Save data for post-processing and computing topological quantities +import os + +save_folder = os.path.join(os.getcwd(), "data") +os.makedirs(save_folder, exist_ok=True) + +time = np.array(post_processing_dict["time"]) +position_history = np.array(post_processing_dict["position"]) +radius_history = np.array(post_processing_dict["radius"]) +director_history = np.array(post_processing_dict["directors"]) + +np.savez( + os.path.join(save_folder, "solenoid_case_data.npz"), + time=time, + position_history=position_history, + radius_history=radius_history, + director_history=director_history, + base_length = np.array(base_length), +) From 36abde7f41f8df99509d523ec85da327bf211737 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Wed, 23 Jun 2021 18:26:29 -0500 Subject: [PATCH 08/14] bugfix:post-processing script fix due to new version of matplotlib --- examples/RodContactCase/post_processing.py | 41 +++++++++++++--------- 1 file changed, 25 insertions(+), 16 deletions(-) diff --git a/examples/RodContactCase/post_processing.py b/examples/RodContactCase/post_processing.py index 8fde465b..56432907 100644 --- a/examples/RodContactCase/post_processing.py +++ b/examples/RodContactCase/post_processing.py @@ -69,6 +69,8 @@ def plot_video_with_surface( 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], @@ -77,7 +79,7 @@ def plot_video_with_surface( s=np.pi * (scaling_factor * inst_radius) ** 2, ) - ax.set_aspect("equal") + # ax.set_aspect("equal") video_name_3D = "3D_" + video_name with writer.saving(fig, video_name_3D, dpi): @@ -88,9 +90,10 @@ def plot_video_with_surface( inst_position, inst_radius = rod_history_unpacker( rod_idx, time_idx ) - inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] - ) + 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], @@ -124,7 +127,8 @@ def plot_video_with_surface( for rod_idx in range(n_visualized_rods): inst_position, inst_radius = rod_history_unpacker(rod_idx, time_idx) - inst_position = 0.5 * (inst_position[..., 1:] + inst_position[..., :-1]) + 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] @@ -148,9 +152,10 @@ def plot_video_with_surface( inst_position, inst_radius = rod_history_unpacker( rod_idx, time_idx ) - inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] - ) + 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]) @@ -186,7 +191,8 @@ def plot_video_with_surface( for rod_idx in range(n_visualized_rods): inst_position, inst_radius = rod_history_unpacker(rod_idx, time_idx) - inst_position = 0.5 * (inst_position[..., 1:] + inst_position[..., :-1]) + 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] @@ -210,9 +216,10 @@ def plot_video_with_surface( inst_position, inst_radius = rod_history_unpacker( rod_idx, time_idx ) - inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] - ) + 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]) @@ -249,7 +256,8 @@ def plot_video_with_surface( for rod_idx in range(n_visualized_rods): inst_position, inst_radius = rod_history_unpacker(rod_idx, time_idx) - inst_position = 0.5 * (inst_position[..., 1:] + inst_position[..., :-1]) + 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] @@ -273,9 +281,10 @@ def plot_video_with_surface( inst_position, inst_radius = rod_history_unpacker( rod_idx, time_idx ) - inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] - ) + 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]) From 4117bbd64742281bf13502d81d91432a7f187318 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Wed, 23 Jun 2021 22:06:33 -0500 Subject: [PATCH 09/14] enhancement: compute twist density for post processing. --- .../RodSelfContact/SolenoidsCase/solenoid_case.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py index 8cdf0d08..6a7b5749 100644 --- a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py +++ b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py @@ -235,6 +235,15 @@ def make_callback(self, system, time, current_step: int): radius_history = np.array(post_processing_dict["radius"]) director_history = np.array(post_processing_dict["directors"]) +# Compute twist density +theta = 2.0 * number_of_rotations * np.pi +angel_vel_scalar = theta / time_twist + +twist_time_interval_start_idx = np.where(time>time_start_twist)[0][0] +twist_time_interval_end_idx = np.where(time<(time_relax + time_twist))[0][-1] + +twist_density = (time[twist_time_interval_start_idx:twist_time_interval_end_idx]-time_start_twist)*angel_vel_scalar * base_radius + np.savez( os.path.join(save_folder, "solenoid_case_data.npz"), time=time, @@ -242,4 +251,5 @@ def make_callback(self, system, time, current_step: int): radius_history=radius_history, director_history=director_history, base_length = np.array(base_length), + twist_density = twist_density, ) From d6491e01ff6f556d009c7a4c777165aa3dd2b8b9 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 8 Jan 2022 16:42:45 -0600 Subject: [PATCH 10/14] update:RodContactCase example files for poisson This commit updates the rod initialization of RodContactCase examples, instead of inputting poisson ratio we use shear modulus --- .../RodRodContact/rod_rod_contact_inclined_validation.py | 3 +++ .../RodRodContact/rod_rod_contact_parallel_validation.py | 5 +++-- .../RodSelfContact/PlectonemesCase/plectoneme_case.py | 3 ++- .../RodSelfContact/SolenoidsCase/solenoid_case.py | 3 ++- 4 files changed, 10 insertions(+), 4 deletions(-) diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py index 204c1d0a..139bbe8f 100644 --- a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py +++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py @@ -32,6 +32,7 @@ class InclinedRodRodContact( nu = 0.001 E = 3e5 poisson_ratio = 0.5 +shear_modulus = E / (poisson_ratio + 1.0) # Rod orientations start = np.zeros(3,) @@ -55,6 +56,7 @@ class InclinedRodRodContact( nu, E, poisson_ratio, + shear_modulus = shear_modulus, ) rod_one.velocity_collection[:] += 0.05 * -normal.reshape(3, 1) @@ -80,6 +82,7 @@ class InclinedRodRodContact( nu, E, poisson_ratio, + shear_modulus=shear_modulus, ) inclined_rod_rod_contact_sim.append(rod_two) diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py index a9637a0c..7ba4b448 100644 --- a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py +++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py @@ -32,6 +32,7 @@ class ParallelRodRodContact( nu = 0.001 E = 3e5 poisson_ratio = 0.5 +shear_modulus = E / (poisson_ratio + 1.0) # Rod orientations start = np.zeros(3,) @@ -54,7 +55,7 @@ class ParallelRodRodContact( density, nu, E, - poisson_ratio, + shear_modulus = shear_modulus, ) rod_one.velocity_collection[:] += 0.05 * -normal.reshape(3, 1) @@ -76,7 +77,7 @@ class ParallelRodRodContact( density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) parallel_rod_rod_contact_sim.append(rod_two) diff --git a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py index dbea5705..4166e676 100644 --- a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py +++ b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py @@ -41,6 +41,7 @@ class PlectonemesCase( nu = 2.0 E = 1e6 poisson_ratio = 0.5 +shear_modulus = E / (poisson_ratio + 1.0) direction = np.array([0.0, 1.0, 0]) normal = np.array([0.0, 0.0, 1.0]) @@ -56,7 +57,7 @@ class PlectonemesCase( density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) sherable_rod.velocity_collection[2, int(n_elem / 2)] -= 1e-4 diff --git a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py index 6a7b5749..0ddb0a7b 100644 --- a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py +++ b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py @@ -40,6 +40,7 @@ class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, Call nu = 2.0 E = 1e6 poisson_ratio = 0.5 +shear_modulus = E / (poisson_ratio + 1.0) direction = np.array([0.0, 1.0, 0]) normal = np.array([0.0, 0.0, 1.0]) @@ -58,7 +59,7 @@ class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, Call density, nu, E, - poisson_ratio, + shear_modulus=shear_modulus, ) From 8790e1df78b1b0fc3b31985f9217e06fe0153b9d Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 8 Jan 2022 16:48:57 -0600 Subject: [PATCH 11/14] black formatting --- elastica/_elastica_numba/_joint.py | 8 +- elastica/joint.py | 202 ++++++++++++++--------------- 2 files changed, 108 insertions(+), 102 deletions(-) diff --git a/elastica/_elastica_numba/_joint.py b/elastica/_elastica_numba/_joint.py index 782064e4..89a692e8 100644 --- a/elastica/_elastica_numba/_joint.py +++ b/elastica/_elastica_numba/_joint.py @@ -1,6 +1,12 @@ import warnings -from elastica.joint import FreeJoint, HingeJoint, FixedJoint, ExternalContact, SelfContact +from elastica.joint import ( + FreeJoint, + HingeJoint, + FixedJoint, + ExternalContact, + SelfContact, +) warnings.warn( "The numba-implementation is included in the default elastica module. Please import without _elastica_numba.", diff --git a/elastica/joint.py b/elastica/joint.py index 93f8f57a..f7a6f5b9 100644 --- a/elastica/joint.py +++ b/elastica/joint.py @@ -361,19 +361,19 @@ def _find_min_dist(x1, e1, x2, e2): @numba.njit(cache=True) def _calculate_contact_forces_rod_rigid_body( - x_collection_rod, - edge_collection_rod, - x_cylinder, - edge_cylinder, - radii_sum, - length_sum, - internal_forces_rod, - external_forces_rod, - external_forces_cylinder, - velocity_rod, - velocity_cylinder, - contact_k, - contact_nu, + x_collection_rod, + edge_collection_rod, + x_cylinder, + edge_cylinder, + radii_sum, + length_sum, + internal_forces_rod, + external_forces_rod, + external_forces_cylinder, + velocity_rod, + velocity_cylinder, + contact_k, + contact_nu, ): # We already pass in only the first n_elem x n_points = x_collection_rod.shape[1] @@ -403,10 +403,10 @@ def _calculate_contact_forces_rod_rigid_body( continue rod_elemental_forces = 0.5 * ( - external_forces_rod[..., i] - + external_forces_rod[..., i + 1] - + internal_forces_rod[..., i] - + internal_forces_rod[..., i + 1] + external_forces_rod[..., i] + + external_forces_rod[..., i + 1] + + internal_forces_rod[..., i] + + internal_forces_rod[..., i + 1] ) equilibrium_forces = -rod_elemental_forces + external_forces_cylinder[..., 0] @@ -420,8 +420,8 @@ def _calculate_contact_forces_rod_rigid_body( contact_force = contact_k * gamma interpenetration_velocity = ( - 0.5 * (velocity_rod[..., i] + velocity_rod[..., i + 1]) - - velocity_cylinder[..., 0] + 0.5 * (velocity_rod[..., i] + velocity_rod[..., i + 1]) + - velocity_cylinder[..., 0] ) contact_damping_force = contact_nu * _dot_product( interpenetration_velocity, distance_vector @@ -429,8 +429,8 @@ def _calculate_contact_forces_rod_rigid_body( # magnitude* direction net_contact_force = ( - normal_force + 0.5 * mask * (contact_damping_force + contact_force) - ) * distance_vector + normal_force + 0.5 * mask * (contact_damping_force + contact_force) + ) * distance_vector # Add it to the rods at the end of the day if i == 0: @@ -449,22 +449,22 @@ def _calculate_contact_forces_rod_rigid_body( @numba.njit(cache=True) def _calculate_contact_forces_rod_rod( - x_collection_rod_one, - radius_rod_one, - length_rod_one, - tangent_rod_one, - velocity_rod_one, - internal_forces_rod_one, - external_forces_rod_one, - x_collection_rod_two, - radius_rod_two, - length_rod_two, - tangent_rod_two, - velocity_rod_two, - internal_forces_rod_two, - external_forces_rod_two, - contact_k, - contact_nu, + x_collection_rod_one, + radius_rod_one, + length_rod_one, + tangent_rod_one, + velocity_rod_one, + internal_forces_rod_one, + external_forces_rod_one, + x_collection_rod_two, + radius_rod_two, + length_rod_two, + tangent_rod_two, + velocity_rod_two, + internal_forces_rod_two, + external_forces_rod_two, + contact_k, + contact_nu, ): # We already pass in only the first n_elem x n_points_rod_one = x_collection_rod_one.shape[1] @@ -505,17 +505,17 @@ def _calculate_contact_forces_rod_rod( continue rod_one_elemental_forces = 0.5 * ( - external_forces_rod_one[..., i] - + external_forces_rod_one[..., i + 1] - + internal_forces_rod_one[..., i] - + internal_forces_rod_one[..., i + 1] + external_forces_rod_one[..., i] + + external_forces_rod_one[..., i + 1] + + internal_forces_rod_one[..., i] + + internal_forces_rod_one[..., i + 1] ) rod_two_elemental_forces = 0.5 * ( - external_forces_rod_two[..., j] - + external_forces_rod_two[..., j + 1] - + internal_forces_rod_two[..., j] - + internal_forces_rod_two[..., j + 1] + external_forces_rod_two[..., j] + + external_forces_rod_two[..., j + 1] + + internal_forces_rod_two[..., j] + + internal_forces_rod_two[..., j + 1] ) equilibrium_forces = -rod_one_elemental_forces + rod_two_elemental_forces @@ -530,8 +530,8 @@ def _calculate_contact_forces_rod_rod( contact_force = contact_k * gamma interpenetration_velocity = 0.5 * ( - (velocity_rod_one[..., i] + velocity_rod_one[..., i + 1]) - - (velocity_rod_two[..., j] + velocity_rod_two[..., j + 1]) + (velocity_rod_one[..., i] + velocity_rod_one[..., i + 1]) + - (velocity_rod_two[..., j] + velocity_rod_two[..., j + 1]) ) contact_damping_force = contact_nu * _dot_product( interpenetration_velocity, distance_vector @@ -539,8 +539,8 @@ def _calculate_contact_forces_rod_rod( # magnitude* direction net_contact_force = ( - normal_force + 0.5 * mask * (contact_damping_force + contact_force) - ) * distance_vector + normal_force + 0.5 * mask * (contact_damping_force + contact_force) + ) * distance_vector # Add it to the rods at the end of the day if i == 0: @@ -566,14 +566,14 @@ def _calculate_contact_forces_rod_rod( @numba.njit(cache=True) def _calculate_contact_forces_self_rod( - x_collection_rod, - radius_rod, - length_rod, - tangent_rod, - velocity_rod, - external_forces_rod, - contact_k, - contact_nu, + x_collection_rod, + radius_rod, + length_rod, + tangent_rod, + velocity_rod, + external_forces_rod, + contact_k, + contact_nu, ): # We already pass in only the first n_elem x n_points_rod = x_collection_rod.shape[1] @@ -618,8 +618,8 @@ def _calculate_contact_forces_self_rod( contact_force = contact_k * gamma interpenetration_velocity = 0.5 * ( - (velocity_rod[..., i] + velocity_rod[..., i + 1]) - - (velocity_rod[..., j] + velocity_rod[..., j + 1]) + (velocity_rod[..., i] + velocity_rod[..., i + 1]) + - (velocity_rod[..., j] + velocity_rod[..., j + 1]) ) contact_damping_force = contact_nu * _dot_product( interpenetration_velocity, distance_vector @@ -627,8 +627,8 @@ def _calculate_contact_forces_self_rod( # magnitude* direction net_contact_force = ( - 0.5 * mask * (contact_damping_force + contact_force) - ) * distance_vector + 0.5 * mask * (contact_damping_force + contact_force) + ) * distance_vector # Add it to the rods at the end of the day # if i == 0: @@ -654,7 +654,7 @@ def _calculate_contact_forces_self_rod( @numba.njit(cache=True) def _aabbs_not_intersecting(aabb_one, aabb_two): - """ Returns true if not intersecting else false""" + """Returns true if not intersecting else false""" if (aabb_one[0, 1] < aabb_two[0, 0]) | (aabb_one[0, 0] > aabb_two[0, 1]): return 1 if (aabb_one[1, 1] < aabb_two[1, 0]) | (aabb_one[1, 0] > aabb_two[1, 1]): @@ -667,13 +667,13 @@ def _aabbs_not_intersecting(aabb_one, aabb_two): @numba.njit(cache=True) def _prune_using_aabbs_rod_rigid_body( - rod_one_position_collection, - rod_one_radius_collection, - rod_one_length_collection, - cylinder_position, - cylinder_director, - cylinder_radius, - cylinder_length, + rod_one_position_collection, + rod_one_radius_collection, + rod_one_length_collection, + cylinder_position, + cylinder_director, + cylinder_radius, + cylinder_length, ): max_possible_dimension = np.zeros((3,)) aabb_rod = np.empty((3, 2)) @@ -683,10 +683,10 @@ def _prune_using_aabbs_rod_rigid_body( ) for i in range(3): aabb_rod[i, 0] = ( - np.min(rod_one_position_collection[i]) - max_possible_dimension[i] + np.min(rod_one_position_collection[i]) - max_possible_dimension[i] ) aabb_rod[i, 1] = ( - np.max(rod_one_position_collection[i]) + max_possible_dimension[i] + np.max(rod_one_position_collection[i]) + max_possible_dimension[i] ) # Is actually Q^T * d but numba complains about performance so we do @@ -698,7 +698,7 @@ def _prune_using_aabbs_rod_rigid_body( for i in range(3): for j in range(3): cylinder_dimensions_in_world_FOR[i] += ( - cylinder_director[j, i, 0] * cylinder_dimensions_in_local_FOR[j] + cylinder_director[j, i, 0] * cylinder_dimensions_in_local_FOR[j] ) max_possible_dimension = np.abs(cylinder_dimensions_in_world_FOR) @@ -709,12 +709,12 @@ def _prune_using_aabbs_rod_rigid_body( @numba.njit(cache=True) def _prune_using_aabbs_rod_rod( - rod_one_position_collection, - rod_one_radius_collection, - rod_one_length_collection, - rod_two_position_collection, - rod_two_radius_collection, - rod_two_length_collection, + rod_one_position_collection, + rod_one_radius_collection, + rod_one_length_collection, + rod_two_position_collection, + rod_two_radius_collection, + rod_two_length_collection, ): max_possible_dimension = np.zeros((3,)) aabb_rod_one = np.empty((3, 2)) @@ -724,10 +724,10 @@ def _prune_using_aabbs_rod_rod( ) for i in range(3): aabb_rod_one[i, 0] = ( - np.min(rod_one_position_collection[i]) - max_possible_dimension[i] + np.min(rod_one_position_collection[i]) - max_possible_dimension[i] ) aabb_rod_one[i, 1] = ( - np.max(rod_one_position_collection[i]) + max_possible_dimension[i] + np.max(rod_one_position_collection[i]) + max_possible_dimension[i] ) max_possible_dimension[...] = np.max(rod_two_radius_collection) + np.max( @@ -736,10 +736,10 @@ def _prune_using_aabbs_rod_rod( for i in range(3): aabb_rod_two[i, 0] = ( - np.min(rod_two_position_collection[i]) - max_possible_dimension[i] + np.min(rod_two_position_collection[i]) - max_possible_dimension[i] ) aabb_rod_two[i, 1] = ( - np.max(rod_two_position_collection[i]) + max_possible_dimension[i] + np.max(rod_two_position_collection[i]) + max_possible_dimension[i] ) return _aabbs_not_intersecting(aabb_rod_two, aabb_rod_one) @@ -774,19 +774,19 @@ def apply_forces(self, rod_one, index_one, rod_two, index_two): # First, check for a global AABB bounding box, and see whether that # intersects if _prune_using_aabbs_rod_rigid_body( - rod_one.position_collection, - rod_one.radius, - rod_one.lengths, - cylinder_two.position_collection, - cylinder_two.director_collection, - cylinder_two.radius, - cylinder_two.length, + rod_one.position_collection, + rod_one.radius, + rod_one.lengths, + cylinder_two.position_collection, + cylinder_two.director_collection, + cylinder_two.radius, + cylinder_two.length, ): return x_cyl = ( - cylinder_two.position_collection[..., 0] - - 0.5 * cylinder_two.length * cylinder_two.director_collection[2, :, 0] + cylinder_two.position_collection[..., 0] + - 0.5 * cylinder_two.length * cylinder_two.director_collection[2, :, 0] ) _calculate_contact_forces_rod_rigid_body( @@ -803,25 +803,25 @@ def apply_forces(self, rod_one, index_one, rod_two, index_two): cylinder_two.velocity_collection, self.k, self.nu, - ) + ) else: # First, check for a global AABB bounding box, and see whether that # intersects if _prune_using_aabbs_rod_rod( - rod_one.position_collection, - rod_one.radius, - rod_one.lengths, - rod_two.position_collection, - rod_two.radius, - rod_two.lengths, + rod_one.position_collection, + rod_one.radius, + rod_one.lengths, + rod_two.position_collection, + rod_two.radius, + rod_two.lengths, ): return _calculate_contact_forces_rod_rod( rod_one.position_collection[ - ..., :-1 + ..., :-1 ], # Discount last node, we want element start position rod_one.radius, rod_one.lengths, @@ -830,7 +830,7 @@ def apply_forces(self, rod_one, index_one, rod_two, index_two): rod_one.internal_forces, rod_one.external_forces, rod_two.position_collection[ - ..., :-1 + ..., :-1 ], # Discount last node, we want element start position rod_two.radius, rod_two.lengths, @@ -866,7 +866,7 @@ def apply_forces(self, rod_one, index_one, rod_two, index_two): _calculate_contact_forces_self_rod( rod_one.position_collection[ - ..., :-1 + ..., :-1 ], # Discount last node, we want element start position rod_one.radius, rod_one.lengths, From dd9a47568c3ac56756b0c48114ea62ac3901316e Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 8 Jan 2022 17:00:07 -0600 Subject: [PATCH 12/14] update:examples/README for RodContactCase --- examples/README.md | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/examples/README.md b/examples/README.md index 02070d75..91923f38 100644 --- a/examples/README.md +++ b/examples/README.md @@ -45,6 +45,17 @@ Examples can serve as a starting template for customized usages. * [Visualization](./Visualization) * __Purpose__: Include simple examples of raytrace rendering data. * __Features__: POVray +* [RodContactCase](./RodContactCase) + * [RodRodContact](./RodContactCase/RodRodContact) + * __Purpose__: Demonstrates contact between two rods, for different initial conditions. + * __Features__: CosseratRod, ExternalContact + * [RodSelfContact](./RodContactCase/RodSelfContact) + * [PlectonemesCase](./RodContactCase/RodSelfContact/PlectonemesCase) + * __Purpose__: Demonstrates rod self contact with Plectoneme example. + * __Features__: CosseratRod, SelonoidsBC, SelfContact + * [SolenoidsCase](./RodContactCase/RodSelfContact/SolenoidsCase) + * __Purpose__: Demonstrates rod self contact with Solenoid example. + * __Features__: CosseratRod, SelonoidsBC, SelfContact ## Advanced Cases From 0ac0596e0e3fba8e8867ee46bdeeaa3ca013c879 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 8 Jan 2022 17:12:42 -0600 Subject: [PATCH 13/14] update: RELEASE.md --- RELEASE.md | 27 +++++++++------------------ 1 file changed, 9 insertions(+), 18 deletions(-) diff --git a/RELEASE.md b/RELEASE.md index 7c950992..3e1f7269 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -1,31 +1,22 @@ -# Release Note (version 0.2) +# Release Note (version 0.2.1) ## Developer Note -Good luck! If it explode, increase nu. :rofl: If it doesn't explode, thoroughly check for the bug. +Contact model between two different rods and rod with itself is implemented. +Testing the contact model is done through simulations. These simulation scripts can be found under +[RodContactCase](./RodContactCase). +However, in future releases we have to add unit tests for contact model functions to test them and increase code coverage. ## Notable Changes -- #84: Block implementation -- #75: Poisson ratio and definition of modulus [PR #26 in public](https://github.com/GazzolaLab/PyElastica/pull/26) -- #95: MuscularFlagella example case is added -- #100: ExportCallBack is added to export the rod-data into file. -- #109: Numpy-only version is now removed. Numba-implementation is set to default. -- #112: Save and load implementation with the example are added. +- #31: Merge contact model to master [PR #40 in public](https://github.com/GazzolaLab/PyElastica/pull/40) ### Release Note
Click to expand -- Block structure is included as part of optimization strategy. -- Different Poisson ratio is supported. -- Contributing guideline is added. -- Update readme -- Add MuscularFlagella example case -- Minimum requirement for dependencies is specified. -- Shear coefficient is corrected. -- Connection index assertion fixed. -- Remove numpy-only version. -- Save/Load example +- Rod-Rod contact and Rod self contact is added. +- Two example cases for rod-rod contact is added, i.e. two rods colliding to each other in space. +- Two example cases for rod self contact is added, i.e. plectonemes and solenoids.
From 75b635149aac30e6712940f2a39369b23a76ab4c Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sat, 8 Jan 2022 22:11:49 -0600 Subject: [PATCH 14/14] Update RELEASE.md --- RELEASE.md | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/RELEASE.md b/RELEASE.md index 3e1f7269..242dcdce 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -20,3 +20,36 @@ However, in future releases we have to add unit tests for contact model function +--- + +# Release Note (version 0.2) + +## Developer Note + +Good luck! If it explode, increase nu. :rofl: If it doesn't explode, thoroughly check for the bug. + +## Notable Changes +- #84: Block implementation +- #75: Poisson ratio and definition of modulus [PR #26 in public](https://github.com/GazzolaLab/PyElastica/pull/26) +- #95: MuscularFlagella example case is added +- #100: ExportCallBack is added to export the rod-data into file. +- #109: Numpy-only version is now removed. Numba-implementation is set to default. +- #112: Save and load implementation with the example are added. + +### Release Note +
+ Click to expand + +- Block structure is included as part of optimization strategy. +- Different Poisson ratio is supported. +- Contributing guideline is added. +- Update readme +- Add MuscularFlagella example case +- Minimum requirement for dependencies is specified. +- Shear coefficient is corrected. +- Connection index assertion fixed. +- Remove numpy-only version. +- Save/Load example + +
+