From 97a3bbb13822de9b824c18340b4971010228d947 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Fri, 11 Mar 2022 15:32:01 -0600 Subject: [PATCH 01/47] renew documentation design to sphinx-book --- docs/conf.py | 7 ++++++- docs/requirements.txt | 3 ++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 9cd8fae78..1fdeb7ba7 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -78,7 +78,12 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. -html_theme = 'sphinx_rtd_theme' +html_theme = 'sphinx_book_theme' +html_theme_option = { + "repository_url": "https://github.com/SophT-Team/SophT/docs", + "use_repository_button": True, +} +html_title = "PyElastica" html_logo = "https://github.com/GazzolaLab/PyElastica/blob/assets/docs/assets/Logo.png?raw=true" #pygments_style = "sphinx" diff --git a/docs/requirements.txt b/docs/requirements.txt index f7684f579..522ad88f3 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,7 +1,8 @@ # File: docs/requirements.txt sphinx==4.4.0 -sphinx_rtd_theme==1.0.0 +#sphinx_rtd_theme==1.0.0 +sphinx-book-theme readthedocs-sphinx-search==0.1.1 sphinx-autodoc-typehints myst-parser From cbb8a9ec2a4ff303cc05aa9acfcfdeb0ac79ffae Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sat, 12 Mar 2022 11:07:29 -0600 Subject: [PATCH 02/47] edit repo url --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index 1fdeb7ba7..f42db9129 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -80,7 +80,7 @@ html_theme = 'sphinx_book_theme' html_theme_option = { - "repository_url": "https://github.com/SophT-Team/SophT/docs", + "repository_url": "https://github.com/GazzolaLab/PyElastica", "use_repository_button": True, } html_title = "PyElastica" From 0d324cf53a859cec2d590b81ea61c84563605069 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sat, 12 Mar 2022 15:18:54 -0600 Subject: [PATCH 03/47] add link-writhe-twist computation to cosserat_rod module Co-authored-by: Arman Tekinalp --- elastica/rod/cosserat_rod.py | 3 +- elastica/rod/knot_theory.py | 711 +++++++++++++++++++++++++++++++++++ 2 files changed, 713 insertions(+), 1 deletion(-) create mode 100644 elastica/rod/knot_theory.py diff --git a/elastica/rod/cosserat_rod.py b/elastica/rod/cosserat_rod.py index 4b6359501..fa0eb12bf 100644 --- a/elastica/rod/cosserat_rod.py +++ b/elastica/rod/cosserat_rod.py @@ -14,6 +14,7 @@ ) from elastica._rotations import _inv_rotate from elastica.rod.factory_function import allocate +from elastica.rod.knot_theory import KnotTheory from elastica._calculus import ( quadrature_kernel_for_block_structure, difference_kernel_for_block_structure, @@ -31,7 +32,7 @@ def _get_z_vector(): return np.array([0.0, 0.0, 1.0]).reshape(3, -1) -class CosseratRod(RodBase): +class CosseratRod(RodBase, KnotTheory): """ Cosserat Rod class. This is the preferred class for rods because it is derived from some of the essential base classes. diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py new file mode 100644 index 000000000..911d8e692 --- /dev/null +++ b/elastica/rod/knot_theory.py @@ -0,0 +1,711 @@ +__docs__ = """ +This script is for computing the link-writhe-twist of a rod using the method 1a from Klenin & Langowski 2000 paper. + +Following codes are adapted from section S2 of Charles et. al. PRL 2019 paper. +""" + +import sys + +if sys.version_info.minor >= 8: + # typing Protocol is introduced in python 3.8 + from typing import Protocol +elif sys.version_info.minor < 8: + # Protocol is implemented in typing_extensions for previous pythons + from typing_extensions import Protocol + +from typing import Union + +from numba import njit +import numpy as np + +from elastica.rod.rod_base import RodBase +from elastica._linalg import _batch_norm, _batch_dot, _batch_cross + + +class KnotTheoryCompatibleProtocol(Protocol): + @property + def position_collection(self) -> np.ndarray: + ... + + @property + def director_collection(self) -> np.ndarray: + ... + + @property + def radius_collection(self) -> np.ndarray: + ... + + @property + def base_length(self) -> np.ndarray: + ... + + +class KnotTheory: + # This mixin should be used in RodBase-derived class that satisfies KnotCompatibleProtocol + MIXIN_PROTOCOL = Union[RodBase, KnotTheoryCompatibleProtocol] + + def compute_twist(self: MIXIN_PROTOCOL): + total_twist, local_twist = compute_twist( + self.position_collection[None, ...], + self.director_collection[0][None, ...], + ) + return total_twist[0] + + def compute_writhe( + self: MIXIN_PROTOCOL, type_of_additional_segment: str = "next_tangent" + ): + return compute_writhe( + self.position_collection[None, ...], + self.rest_lengths.sum(), + type_of_additional_segment, + )[0] + + def compute_link( + self: MIXIN_PROTOCOL, type_of_additional_segment: str = "next_tangent" + ): + return compute_link( + self.position_collection[None, ...], + self.director_collection[0][None, ...], + self.radius[None, ...], + self.rest_lengths.sum(), + type_of_additional_segment, + )[0] + + +def compute_twist(center_line, normal_collection): + """ + Compute the twist of a rod, using center_line and normal collection. Methods used in this function is + adapted from method 2a Klenin & Langowski 2000 paper. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + normal_collection : numpy.ndarray + 3D (time, 3, n_elems) array containing data with 'float' type. + Time history of rod elements normal direction. + + Warnings + ------- + If center line is straight, although the normals of each element is pointing different direction computed twist + will be zero. + + Returns + ------- + + """ + + total_twist, local_twist = _compute_twist(center_line, normal_collection) + + return total_twist, local_twist + + +@njit(cache=True) +def _compute_twist(center_line, normal_collection): + """ + Compute the twist of a rod, using center_line and normal collection. Methods used in this function is + adapted from method 2a Klenin & Langowski 2000 paper. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + normal_collection : numpy.ndarray + 3D (time, 3, n_elems) array containing data with 'float' type. + Time history of rod elements normal direction. + + Returns + ------- + + """ + + timesize, _, blocksize = center_line.shape + + total_twist = np.zeros((timesize)) + local_twist = np.zeros((timesize, blocksize - 2)) # Only consider interior nodes + + # compute s vector + for k in range(timesize): + # s is a secondary vector field. + s = center_line[k, :, 1:] - center_line[k, :, :-1] + # Compute tangents + tangent = s / _batch_norm(s) + + # Compute the projection of normal collection (d1) on normal-binormal plane. + projection_of_normal_collection = ( + normal_collection[k, :, :] + - _batch_dot(tangent, normal_collection[k, :, :]) * tangent + ) + projection_of_normal_collection /= _batch_norm(projection_of_normal_collection) + + # Eq 27 in Klenin & Langowski 2000 + # p is defined on interior nodes + p = _batch_cross(s[:, :-1], s[:, 1:]) + p /= _batch_norm(p) + + # Compute the angle we need to turn d1 around s to get p + # sign part tells whether d1 must be rotated ccw(+) or cw(-) around s + alpha = np.sign( + _batch_dot( + _batch_cross(projection_of_normal_collection[:, :-1], p), s[:, :-1] + ) + ) * np.arccos(_batch_dot(projection_of_normal_collection[:, :-1], p)) + + gamma = np.sign( + _batch_dot( + _batch_cross(p, projection_of_normal_collection[:, 1:]), s[:, 1:] + ) + ) * np.arccos(_batch_dot(projection_of_normal_collection[:, 1:], p)) + + # An angle 1 is a full rotation, 0.5 is rotation by pi, 0.25 is pi/2 etc. + alpha /= 2 * np.pi + gamma /= 2 * np.pi + twist_temp = alpha + gamma + # Make sure twist is between (-1/2 to 1/2) as defined in pg 313 Klenin & Langowski 2000 + idx = np.where(twist_temp > 0.5)[0] + twist_temp[idx] -= 1 + idx = np.where(twist_temp < -0.5)[0] + twist_temp[idx] += 1 + + # Check if there is any Nan. Nan's appear when rod tangents are parallel to each other. + idx = np.where(np.isnan(twist_temp))[0] + twist_temp[idx] = 0.0 + + local_twist[k, :] = twist_temp + total_twist[k] = np.sum(twist_temp) + + return total_twist, local_twist + + +def compute_writhe(center_line, segment_length, type_of_additional_segment): + """ + This function computes the total writhe history of a rod. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + segment_length : float + Length of added segments. + type_of_additional_segment : str + Determines the method to compute new segments (elements) added to the rod. + Valid inputs are "next_tangent", "end_to_end", "net_tangent", otherwise program uses the center line. + + Returns + ------- + total_writhe + + """ + + center_line_with_added_segments, _, _ = compute_additional_segment( + center_line, segment_length, type_of_additional_segment + ) + + """ + Total writhe of a rod is computed using the method 1a from Klenin & Langowski 2000 + """ + total_writhe = _compute_writhe(center_line_with_added_segments) + + return total_writhe + + +@njit(cache=True) +def _compute_writhe(center_line): + """ + This function computes the total writhe history of a rod. + Equations used are from method 1a from Klenin & Langowski 2000 paper. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + + Returns + ------- + + """ + + time, _, blocksize = center_line.shape + + omega_star = np.zeros((blocksize - 2, blocksize - 1)) + segment_writhe = np.zeros((blocksize - 2, blocksize - 1)) + total_writhe = np.zeros((time)) + + # Compute the writhe between each pair first. + for k in range(time): + for i in range(blocksize - 2): + for j in range(i + 1, blocksize - 1): + + point_one = center_line[k, :, i] + point_two = center_line[k, :, i + 1] + point_three = center_line[k, :, j] + point_four = center_line[k, :, j + 1] + + # Eq 15 in Klenin & Langowski 2000 + r12 = point_two - point_one + r34 = point_four - point_three + r14 = point_four - point_one + r13 = point_three - point_one + r23 = point_three - point_two + r24 = point_four - point_two + + n1 = np.cross(r13, r14) + n1 /= np.linalg.norm(n1) + n2 = np.cross(r14, r24) + n2 /= np.linalg.norm(n2) + n3 = np.cross(r24, r23) + n3 /= np.linalg.norm(n3) + n4 = np.cross(r23, r13) + n4 /= np.linalg.norm(n4) + + # Eq 16a in Klenin & Langowski 2000 + omega_star[i, j] = ( + np.arcsin(np.dot(n1, n2)) + + np.arcsin(np.dot(n2, n3)) + + np.arcsin(np.dot(n3, n4)) + + np.arcsin(np.dot(n4, n1)) + ) + + if np.isnan(omega_star[i, j]): + omega_star[i, j] = 0 + + # Eq 16b in Klenin & Langowski 2000 + segment_writhe[i, j] = ( + omega_star[i, j] + * np.sign(np.dot(np.cross(r34, r12), r13)) + / (4 * np.pi) + ) + + # Compute the total writhe + # Eq 13 in Klenin & Langowski 2000 + total_writhe[k] = 2 * np.sum(segment_writhe) + + return total_writhe + + +def compute_link( + center_line: np.ndarray, + normal_collection: np.ndarray, + radius: np.ndarray, + segment_length: float, + type_of_additional_segment: str, +): + """ + Compute the link of a rod. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + normal_collection : numpy.ndarray + 3D (time, 3, n_elems) array containing data with 'float' type. + Time history of rod elements normal direction. + radius : numpy.ndarray + 2D (time, n_elems) array containing data with 'float' type. + Time history of rod element radius. + segment_length : float + Length of added segments. + type_of_additional_segment : str + Determines the method to compute new segments (elements) added to the rod. + Valid inputs are "next_tangent", "end_to_end", "net_tangent", otherwise program uses the center line. + """ + + # Compute auxiliary line + auxiliary_line = _compute_auxiliary_line(center_line, normal_collection, radius) + + # Add segments at the beginning and end of the rod center line and auxiliary line. + ( + center_line_with_added_segments, + beginning_direction, + end_direction, + ) = compute_additional_segment( + center_line, segment_length, type_of_additional_segment + ) + auxiliary_line_with_added_segments = _compute_auxiliary_line_added_segments( + beginning_direction, end_direction, auxiliary_line, segment_length + ) + + """ + Total link of a rod is computed using the method 1a from Klenin & Langowski 2000 + """ + total_link = _compute_link( + center_line_with_added_segments, auxiliary_line_with_added_segments + ) + + return total_link + + +@njit(cache=True) +def _compute_auxiliary_line(center_line, normal_collection, radius): + """ + This function computes the auxiliary line using rod center line and normal collection. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + normal_collection : numpy.ndarray + 3D (time, 3, n_elems) array containing data with 'float' type. + Time history of rod elements normal direction. + radius : numpy.ndarray + 2D (time, n_elems) array containing data with 'float' type. + Time history of rod element radius. + + Returns + ------- + + """ + time, _, blocksize = center_line.shape + auxiliary_line = np.zeros(center_line.shape) + projection_of_normal_collection = np.zeros((3, blocksize)) + radius_on_nodes = np.zeros((blocksize)) + + for i in range(time): + tangent = center_line[i, :, 1:] - center_line[i, :, :-1] + tangent /= _batch_norm(tangent) + # Compute the projection of normal collection (d1) on normal-binormal plane. + projection_of_normal_collection_temp = ( + normal_collection[i, :, :] + - _batch_dot(tangent, normal_collection[i, :, :]) * tangent + ) + projection_of_normal_collection_temp /= _batch_norm( + projection_of_normal_collection_temp + ) + + # First node have the same direction with second node. They share the same element. + # TODO: Instead of this maybe we should use the trapezoidal rule or averaging operator for normal and radius. + projection_of_normal_collection[:, 0] = projection_of_normal_collection_temp[ + :, 0 + ] + projection_of_normal_collection[:, 1:] = projection_of_normal_collection_temp[:] + radius_on_nodes[0] = radius[i, 0] + radius_on_nodes[1:] = radius[i, :] + + auxiliary_line[i, :, :] = ( + radius_on_nodes * projection_of_normal_collection + center_line[i, :, :] + ) + + return auxiliary_line + + +@njit(cache=True) +def _compute_link(center_line, auxiliary_line): + """ + This function computes the total link history of a rod. + Equations used are from method 1a from Klenin & Langowski 2000 paper. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + auxiliary_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of auxiliary line. + + Returns + ------- + + """ + timesize, _, blocksize_center_line = center_line.shape + blocksize_auxiliary_line = auxiliary_line.shape[-1] + + omega_star = np.zeros((blocksize_center_line - 1, blocksize_auxiliary_line - 1)) + segment_link = np.zeros((blocksize_center_line - 1, blocksize_auxiliary_line - 1)) + total_link = np.zeros((timesize)) + + # Compute the writhe between each pair first. + for k in range(timesize): + for i in range(blocksize_center_line - 1): + for j in range(blocksize_auxiliary_line - 1): + + point_one = center_line[k, :, i] + point_two = center_line[k, :, i + 1] + point_three = auxiliary_line[k, :, j] + point_four = auxiliary_line[k, :, j + 1] + + # Eq 15 in Klenin & Langowski 2000 + r12 = point_two - point_one + r34 = point_four - point_three + r14 = point_four - point_one + r13 = point_three - point_one + r23 = point_three - point_two + r24 = point_four - point_two + + n1 = np.cross(r13, r14) + n1 /= np.linalg.norm(n1) + n2 = np.cross(r14, r24) + n2 /= np.linalg.norm(n2) + n3 = np.cross(r24, r23) + n3 /= np.linalg.norm(n3) + n4 = np.cross(r23, r13) + n4 /= np.linalg.norm(n4) + + # Eq 16a in Klenin & Langowski 2000 + omega_star[i, j] = ( + np.arcsin(np.dot(n1, n2)) + + np.arcsin(np.dot(n2, n3)) + + np.arcsin(np.dot(n3, n4)) + + np.arcsin(np.dot(n4, n1)) + ) + + if np.isnan(omega_star[i, j]): + omega_star[i, j] = 0 + + # Eq 16b in Klenin & Langowski 2000 + segment_link[i, j] = ( + omega_star[i, j] + * np.sign(np.dot(np.cross(r34, r12), r13)) + / (4 * np.pi) + ) + + # Compute the total writhe + # Eq 6 in Klenin & Langowski 2000 + # Unlike the writhe, link computed using two curves so we do not multiply by 2 + total_link[k] = np.sum(segment_link) + + return total_link + + +@njit(cache=True) +def _compute_auxiliary_line_added_segments( + beginning_direction, end_direction, auxiliary_line, segment_length +): + """ + This code is for computing position of added segments to the auxiliary line. + + Parameters + ---------- + beginning_direction : numpy.ndarray + 2D (time, 3) array containing data with 'float' type. + Time history of center line tangent at the beginning. + end_direction : numpy.ndarray + 2D (time, 3) array containing data with 'float' type. + Time history of center line tangent at the end. + auxiliary_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of auxiliary line. + segment_length : float + Length of added segments. + + Returns + ------- + + """ + timesize, _, blocksize = auxiliary_line.shape + + new_auxiliary_line = np.zeros((timesize, 3, blocksize + 2)) + + new_auxiliary_line[:, :, 1:-1] = auxiliary_line + + new_auxiliary_line[:, :, 0] = ( + auxiliary_line[:, :, 0] + beginning_direction * segment_length + ) + + new_auxiliary_line[:, :, -1] = ( + auxiliary_line[:, :, -1] + end_direction * segment_length + ) + + return new_auxiliary_line + + +def compute_additional_segment(center_line, segment_length, type_of_additional_segment): + """ + This function adds two points at the end of center line. Distance from the center line is given by segment_length. + Direction from center line to the new point locations can be computed using 3 methods, which can be selected by + type. For more details section S2 of Charles et. al. PRL 2019 paper. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + segment_length : float + Length of added segments. + type_of_additional_segment : str + Determines the method to compute new segments (elements) added to the rod. + Valid inputs are "next_tangent", "end_to_end", "net_tangent", otherwise it returns the center line. + + Returns + ------- + + """ + + if type_of_additional_segment == "next_tangent": + ( + center_line, + beginning_direction, + end_direction, + ) = _compute_additional_segment_using_next_tangent(center_line, segment_length) + elif type_of_additional_segment == "end_to_end": + ( + center_line, + beginning_direction, + end_direction, + ) = _compute_additional_segment_using_end_to_end(center_line, segment_length) + elif type_of_additional_segment == "net_tangent": + ( + center_line, + beginning_direction, + end_direction, + ) = _compute_additional_segment_using_net_tangent(center_line, segment_length) + else: + print("Additional segments are not used") + beginning_direction = np.zeros((center_line.shape[0], 3)) + end_direction = np.zeros((center_line.shape[0], 3)) + + return center_line, beginning_direction, end_direction + + +@njit(cache=True) +def _compute_additional_segment_using_next_tangent(center_line, segment_length): + """ + This function adds a two new point at the begining and end of the center line. Distance of these points are + given in segment_length. Direction of these points are computed using the rod tangents at the begining and end. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + segment_length : float + Length of added segments. + + Returns + ------- + + """ + + timesize, _, blocksize = center_line.shape + new_center_line = np.zeros( + (timesize, 3, blocksize + 2) + ) # +2 is for added two new points + beginning_direction = np.zeros((timesize, 3)) + end_direction = np.zeros((timesize, 3)) + + for i in range(timesize): + # Direction of the additional point at the beginning of the rod + direction_of_rod_begin = center_line[i, :, 0] - center_line[i, :, 1] + direction_of_rod_begin /= np.linalg.norm(direction_of_rod_begin) + + # Direction of the additional point at the end of the rod + direction_of_rod_end = center_line[i, :, -1] - center_line[i, :, -2] + direction_of_rod_end /= np.linalg.norm(direction_of_rod_end) + + first_point = center_line[i, :, 0] + segment_length * direction_of_rod_begin + last_point = center_line[i, :, -1] + segment_length * direction_of_rod_end + + new_center_line[i, :, 1:-1] = center_line[i, :, :] + new_center_line[i, :, 0] = first_point + new_center_line[i, :, -1] = last_point + beginning_direction[i, :] = direction_of_rod_begin + end_direction[i, :] = direction_of_rod_end + + return new_center_line, beginning_direction, end_direction + + +@njit(cache=True) +def _compute_additional_segment_using_end_to_end(center_line, segment_length): + """ + This function adds a two new point at the begining and end of the center line. Distance of these points are + given in segment_length. Direction of these points are computed using the rod node end positions. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + segment_length : float + Length of added segments. + + Returns + ------- + + """ + timesize, _, blocksize = center_line.shape + new_center_line = np.zeros( + (timesize, 3, blocksize + 2) + ) # +2 is for added two new points + beginning_direction = np.zeros((timesize, 3)) + end_direction = np.zeros((timesize, 3)) + + for i in range(timesize): + # Direction of the additional point at the beginning of the rod + direction_of_rod_begin = center_line[i, :, 0] - center_line[i, :, -1] + direction_of_rod_begin /= np.linalg.norm(direction_of_rod_begin) + + # Direction of the additional point at the end of the rod + direction_of_rod_end = -direction_of_rod_begin + + first_point = center_line[i, :, 0] + segment_length * direction_of_rod_begin + last_point = center_line[i, :, -1] + segment_length * direction_of_rod_end + + new_center_line[i, :, 1:-1] = center_line[i, :, :] + new_center_line[i, :, 0] = first_point + new_center_line[i, :, -1] = last_point + + beginning_direction[i, :] = direction_of_rod_begin + end_direction[i, :] = direction_of_rod_end + + return new_center_line, beginning_direction, end_direction + + +@njit(cache=True) +def _compute_additional_segment_using_net_tangent(center_line, segment_length): + """ + This function adds a two new point at the begining and end of the center line. Distance of these points are + given in segment_length. Direction of these points are point wise avarege of nodes at the first and second half + of the rod. + + Parameters + ---------- + center_line : numpy.ndarray + 3D (time, 3, n_nodes) array containing data with 'float' type. + Time history of rod node positions. + segment_length : float + Length of added segments. + + Returns + ------- + + """ + timesize, _, blocksize = center_line.shape + new_center_line = np.zeros( + (timesize, 3, blocksize + 2) + ) # +2 is for added two new points + beginning_direction = np.zeros((timesize, 3)) + end_direction = np.zeros((timesize, 3)) + + for i in range(timesize): + # Direction of the additional point at the beginning of the rod + n_nodes_begin = int(np.floor(blocksize / 2)) + average_begin = ( + np.sum(center_line[i, :, :n_nodes_begin], axis=1) / n_nodes_begin + ) + n_nodes_end = int(np.ceil(blocksize / 2)) + average_end = np.sum(center_line[i, :, n_nodes_end:], axis=1) / ( + blocksize - n_nodes_end + 1 + ) + + direction_of_rod_begin = average_begin - average_end + direction_of_rod_begin /= np.linalg.norm(direction_of_rod_begin) + + direction_of_rod_end = -direction_of_rod_begin + + first_point = center_line[i, :, 0] + segment_length * direction_of_rod_begin + last_point = center_line[i, :, -1] + segment_length * direction_of_rod_end + + new_center_line[i, :, 1:-1] = center_line[i, :, :] + new_center_line[i, :, 0] = first_point + new_center_line[i, :, -1] = last_point + + beginning_direction[i, :] = direction_of_rod_begin + end_direction[i, :] = direction_of_rod_end + + return new_center_line, beginning_direction, end_direction From 0d9f4968a65816aacab771e7092d0c5a5754268e Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sat, 12 Mar 2022 19:28:08 -0600 Subject: [PATCH 04/47] adjust alias for import --- elastica/rod/__init__.py | 1 + elastica/rod/knot_theory.py | 1 + 2 files changed, 2 insertions(+) diff --git a/elastica/rod/__init__.py b/elastica/rod/__init__.py index 8d846bf60..08263853b 100644 --- a/elastica/rod/__init__.py +++ b/elastica/rod/__init__.py @@ -1,5 +1,6 @@ __doc__ = """Rod classes and its data structures """ +from elastica.rod.knot_theory import * from elastica.rod.data_structures import * from elastica.rod.rod_base import RodBase diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index 911d8e692..6fffe5808 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -3,6 +3,7 @@ Following codes are adapted from section S2 of Charles et. al. PRL 2019 paper. """ +__all__ = ["KnotTheoryCompatibleProtocol", "KnotTheory", "compute_twist", "compute_link", "compute_writhe"] import sys From 25ffb8a723e13bc07001882ba9e1ef8e4e34dbf4 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Tue, 15 Mar 2022 23:34:00 -0500 Subject: [PATCH 05/47] wip: test for knot theory --- Makefile | 2 ++ tests/test_boundary_conditions.py | 2 +- tests/test_interaction.py | 2 +- tests/{ => test_rod}/test_cosserat_rod.py | 0 tests/test_rod/test_knot_theory.py | 30 ++++++++++++++++++++ tests/{test_rod.py => test_rod/test_rods.py} | 0 6 files changed, 34 insertions(+), 2 deletions(-) rename tests/{ => test_rod}/test_cosserat_rod.py (100%) create mode 100644 tests/test_rod/test_knot_theory.py rename tests/{test_rod.py => test_rod/test_rods.py} (100%) diff --git a/Makefile b/Makefile index b45ffe6c2..a16207652 100644 --- a/Makefile +++ b/Makefile @@ -12,6 +12,8 @@ flake8: @flake8 --version @flake8 elastica tests +test: + @python -m pytest all:black flake8 ci:black_check flake8 diff --git a/tests/test_boundary_conditions.py b/tests/test_boundary_conditions.py index d8143d786..7e3144c57 100644 --- a/tests/test_boundary_conditions.py +++ b/tests/test_boundary_conditions.py @@ -3,7 +3,7 @@ # System imports import numpy as np -from test_rod import MockTestRod +from test_rod.test_rods import MockTestRod from elastica.boundary_conditions import ( ConstraintBase, FreeBC, diff --git a/tests/test_interaction.py b/tests/test_interaction.py index eb010344a..82ac01a80 100644 --- a/tests/test_interaction.py +++ b/tests/test_interaction.py @@ -12,7 +12,7 @@ SlenderBodyTheory, ) -from test_rod import MockTestRod +from test_rod.test_rods import MockTestRod class BaseRodClass(MockTestRod): diff --git a/tests/test_cosserat_rod.py b/tests/test_rod/test_cosserat_rod.py similarity index 100% rename from tests/test_cosserat_rod.py rename to tests/test_rod/test_cosserat_rod.py diff --git a/tests/test_rod/test_knot_theory.py b/tests/test_rod/test_knot_theory.py new file mode 100644 index 000000000..a80c4bbec --- /dev/null +++ b/tests/test_rod/test_knot_theory.py @@ -0,0 +1,30 @@ +__doc__ = """ Knot Theory class testing """ + +import pytest +import numpy as np +from numpy.testing import assert_allclose + +from elastica.rod.data_structures import _bootstrap_from_data +from elastica.rod.data_structures import ( + _KinematicState, + _DynamicState, +) +from elastica.utils import MaxDimension + +from test_rod.test_rods import MockTestRod + +from elastica.rod.rod_base import RodBase +from elastica.rod.knot_theory import KnotTheory + +class TestRodWithKnotTheory(MockTestRod, KnotTheory): + def __init__(self): + super().__init__() + self.radius = np.random.randn(MaxDimension.value(), self.n_elem) + +def test_knot_theory_mixin_methods(): + rod = TestRodWithKnotTheory() + assert hasattr(rod, "MIXIN_PROTOCOL"), "Expected to mix-in variables: MIXIN_PROTOCOL" + assert hasattr(rod, "compute_writhe"), "Expected to mix-in functionals into the rod class: compute_writhe" + assert hasattr(rod, "compute_twist"), "Expected to mix-in functionals into the rod class: compute_twist" + assert hasattr(rod, "compute_link"), "Expected to mix-in functionals into the rod class: compute_link" + diff --git a/tests/test_rod.py b/tests/test_rod/test_rods.py similarity index 100% rename from tests/test_rod.py rename to tests/test_rod/test_rods.py From aa3e9b85f71c1f1077f3d9ce19f88409d6733d8b Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Tue, 15 Mar 2022 23:34:17 -0500 Subject: [PATCH 06/47] wip: documentation on knot theory --- elastica/rod/knot_theory.py | 34 +++++++++++++++++++++++++++--- tests/test_rod/test_knot_theory.py | 19 ++++++++++++----- 2 files changed, 45 insertions(+), 8 deletions(-) diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index 6fffe5808..e0de0fd65 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -1,9 +1,15 @@ __docs__ = """ -This script is for computing the link-writhe-twist of a rod using the method 1a from Klenin & Langowski 2000 paper. +This script is for computing the link-writhe-twist of a rod using the method from Klenin & Langowski 2000 paper. Following codes are adapted from section S2 of Charles et. al. PRL 2019 paper. """ -__all__ = ["KnotTheoryCompatibleProtocol", "KnotTheory", "compute_twist", "compute_link", "compute_writhe"] +__all__ = [ + "KnotTheoryCompatibleProtocol", + "KnotTheory", + "compute_twist", + "compute_link", + "compute_writhe", +] import sys @@ -24,6 +30,11 @@ class KnotTheoryCompatibleProtocol(Protocol): + """KnotTheoryCompatibleProtocol + + Required properties to use KnotTheory mixin + """ + @property def position_collection(self) -> np.ndarray: ... @@ -46,6 +57,7 @@ class KnotTheory: MIXIN_PROTOCOL = Union[RodBase, KnotTheoryCompatibleProtocol] def compute_twist(self: MIXIN_PROTOCOL): + """compute_twist""" total_twist, local_twist = compute_twist( self.position_collection[None, ...], self.director_collection[0][None, ...], @@ -55,6 +67,14 @@ def compute_twist(self: MIXIN_PROTOCOL): def compute_writhe( self: MIXIN_PROTOCOL, type_of_additional_segment: str = "next_tangent" ): + """compute_writhe + + Parameters + ---------- + type_of_additional_segment : str + Determines the method to compute new segments (elements) added to the rod. + Valid inputs are "next_tangent", "end_to_end", "net_tangent", otherwise program uses the center line. + """ return compute_writhe( self.position_collection[None, ...], self.rest_lengths.sum(), @@ -64,6 +84,14 @@ def compute_writhe( def compute_link( self: MIXIN_PROTOCOL, type_of_additional_segment: str = "next_tangent" ): + """compute_link + + Parameters + ---------- + type_of_additional_segment : str + Determines the method to compute new segments (elements) added to the rod. + Valid inputs are "next_tangent", "end_to_end", "net_tangent", otherwise program uses the center line. + """ return compute_link( self.position_collection[None, ...], self.director_collection[0][None, ...], @@ -88,7 +116,7 @@ def compute_twist(center_line, normal_collection): Time history of rod elements normal direction. Warnings - ------- + -------- If center line is straight, although the normals of each element is pointing different direction computed twist will be zero. diff --git a/tests/test_rod/test_knot_theory.py b/tests/test_rod/test_knot_theory.py index a80c4bbec..abdca6b7b 100644 --- a/tests/test_rod/test_knot_theory.py +++ b/tests/test_rod/test_knot_theory.py @@ -16,15 +16,24 @@ from elastica.rod.rod_base import RodBase from elastica.rod.knot_theory import KnotTheory + class TestRodWithKnotTheory(MockTestRod, KnotTheory): def __init__(self): super().__init__() self.radius = np.random.randn(MaxDimension.value(), self.n_elem) + def test_knot_theory_mixin_methods(): rod = TestRodWithKnotTheory() - assert hasattr(rod, "MIXIN_PROTOCOL"), "Expected to mix-in variables: MIXIN_PROTOCOL" - assert hasattr(rod, "compute_writhe"), "Expected to mix-in functionals into the rod class: compute_writhe" - assert hasattr(rod, "compute_twist"), "Expected to mix-in functionals into the rod class: compute_twist" - assert hasattr(rod, "compute_link"), "Expected to mix-in functionals into the rod class: compute_link" - + assert hasattr( + rod, "MIXIN_PROTOCOL" + ), "Expected to mix-in variables: MIXIN_PROTOCOL" + assert hasattr( + rod, "compute_writhe" + ), "Expected to mix-in functionals into the rod class: compute_writhe" + assert hasattr( + rod, "compute_twist" + ), "Expected to mix-in functionals into the rod class: compute_twist" + assert hasattr( + rod, "compute_link" + ), "Expected to mix-in functionals into the rod class: compute_link" From ad033c7a5837eb22152b4d1590d96d9f28d580fd Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 16 Mar 2022 16:25:14 -0500 Subject: [PATCH 07/47] Create parallel_connection.py --- .../experimental/connection_contact_joint/parallel_connection.py | 1 + 1 file changed, 1 insertion(+) create mode 100644 elastica/experimental/connection_contact_joint/parallel_connection.py diff --git a/elastica/experimental/connection_contact_joint/parallel_connection.py b/elastica/experimental/connection_contact_joint/parallel_connection.py new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/elastica/experimental/connection_contact_joint/parallel_connection.py @@ -0,0 +1 @@ + From d12acde3b5a63a8bd0cf8e0ff937ad8d5fcdfef0 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 16 Mar 2022 16:25:41 -0500 Subject: [PATCH 08/47] Update version.py --- elastica/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elastica/version.py b/elastica/version.py index 0bcc5d990..ae50f1f69 100644 --- a/elastica/version.py +++ b/elastica/version.py @@ -1 +1 @@ -VERSION = "0.2.2" +VERSION = "0.2.3" From 02480d6fc75e0c724fc37526a1617eee034bb4e9 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Wed, 16 Mar 2022 23:05:02 -0500 Subject: [PATCH 09/47] enhancement:Add parallel_connection class This class is for connecting two parallel rods at their surface. Current implementation is at experimental stage. --- .../parallel_connection.py | 296 ++++++++++++++++++ 1 file changed, 296 insertions(+) diff --git a/elastica/experimental/connection_contact_joint/parallel_connection.py b/elastica/experimental/connection_contact_joint/parallel_connection.py index 8b1378917..d9a74ab67 100644 --- a/elastica/experimental/connection_contact_joint/parallel_connection.py +++ b/elastica/experimental/connection_contact_joint/parallel_connection.py @@ -1 +1,297 @@ +__doc__ = """Contains SurfaceJointSideBySide class which connects two parallel rods .""" +import numpy as np +import numba +from numba import njit +from elastica.joint import FreeJoint +# Join the two rods +from elastica._linalg import ( + _batch_norm, + _batch_cross, + _batch_matvec, + _batch_dot, + _batch_matmul, + _batch_matrix_transpose, +) + + +def get_connection_vector_straight_straight_rod( + rod_one, + rod_two, +): + + # Compute rod element positions + rod_one_element_position = 0.5 * ( + rod_one.position_collection[..., 1:] + rod_one.position_collection[..., :-1] + ) + rod_two_element_position = 0.5 * ( + rod_two.position_collection[..., 1:] + rod_two.position_collection[..., :-1] + ) + + # Lets get the distance between rod elements + distance_vector_rod_one_to_rod_two = ( + rod_two_element_position - rod_one_element_position + ) + distance_vector_rod_one_to_rod_two_norm = _batch_norm( + distance_vector_rod_one_to_rod_two + ) + distance_vector_rod_one_to_rod_two /= distance_vector_rod_one_to_rod_two_norm + + distance_vector_rod_two_to_rod_one = -distance_vector_rod_one_to_rod_two + + rod_one_direction_vec_in_material_frame = _batch_matvec( + rod_one.director_collection, distance_vector_rod_one_to_rod_two + ) + rod_two_direction_vec_in_material_frame = _batch_matvec( + rod_two.director_collection, distance_vector_rod_two_to_rod_one + ) + + offset_btw_rods = distance_vector_rod_one_to_rod_two_norm - ( + rod_one.radius + rod_two.radius + ) + + return ( + rod_one_direction_vec_in_material_frame, + rod_two_direction_vec_in_material_frame, + offset_btw_rods, + ) + + +class SurfaceJointSideBySide(FreeJoint): + """ + TODO: documentation + """ + + def __init__( + self, + k, + nu, + k_repulsive, + rod_one_direction_vec_in_material_frame, + rod_two_direction_vec_in_material_frame, + offset_btw_rods, + **kwargs, + ): + super().__init__(np.array(k), np.array(nu)) + + self.k_repulsive = np.array(k_repulsive) + + self.offset_btw_rods = np.array(offset_btw_rods) + + self.rod_one_direction_vec_in_material_frame = np.array( + rod_one_direction_vec_in_material_frame + ).T + self.rod_two_direction_vec_in_material_frame = np.array( + rod_two_direction_vec_in_material_frame + ).T + + # Apply force is same as free joint + def apply_forces(self, rod_one, index_one, rod_two, index_two): + # TODO: documentation + + (self.rod_one_rd2, self.rod_two_rd2, self.spring_force,) = self._apply_forces( + self.k, + self.nu, + self.k_repulsive, + index_one, + index_two, + self.rod_one_direction_vec_in_material_frame, + self.rod_two_direction_vec_in_material_frame, + self.offset_btw_rods, + rod_one.director_collection, + rod_two.director_collection, + rod_one.position_collection, + rod_two.position_collection, + rod_one.radius, + rod_two.radius, + rod_one.dilatation, + rod_two.dilatation, + rod_one.velocity_collection, + rod_two.velocity_collection, + rod_one.external_forces, + rod_two.external_forces, + ) + + @staticmethod + @njit(cache=True) + def _apply_forces( + k, + nu, + k_repulsive, + index_one, + index_two, + rod_one_direction_vec_in_material_frame, + rod_two_direction_vec_in_material_frame, + rest_offset_btw_rods, + rod_one_director_collection, + rod_two_director_collection, + rod_one_position_collection, + rod_two_position_collection, + rod_one_radius, + rod_two_radius, + rod_one_dilatation, + rod_two_dilatation, + rod_one_velocity_collection, + rod_two_velocity_collection, + rod_one_external_forces, + rod_two_external_forces, + ): + + rod_one_to_rod_two_connection_vec = ( + rod_one_director_collection[:, :, index_one].T + @ rod_one_direction_vec_in_material_frame + ) + rod_two_to_rod_one_connection_vec = ( + rod_two_director_collection[:, :, index_two].T + @ rod_two_direction_vec_in_material_frame + ) + + # Compute element positions + rod_one_element_position = 0.5 * ( + rod_one_position_collection[:, index_one] + + rod_one_position_collection[:, index_one + 1] + ) + rod_two_element_position = 0.5 * ( + rod_two_position_collection[:, index_two] + + rod_two_position_collection[:, index_two + 1] + ) + + # If there is an offset between rod one and rod two surface, then it should change as a function of dilatation. + offset_rod_one = ( + 0.5 * rest_offset_btw_rods / np.sqrt(rod_one_dilatation[index_one]) + ) + offset_rod_two = ( + 0.5 * rest_offset_btw_rods / np.sqrt(rod_two_dilatation[index_two]) + ) + + # Compute vector r*d2 (radius * connection vector) for each rod and element + rod_one_rd2 = rod_one_to_rod_two_connection_vec * ( + rod_one_radius[index_one] + offset_rod_one + ) + rod_two_rd2 = rod_two_to_rod_one_connection_vec * ( + rod_two_radius[index_two] + offset_rod_two + ) + + # Compute connection points on the rod surfaces + surface_position_rod_one = rod_one_element_position + rod_one_rd2 + surface_position_rod_two = rod_two_element_position + rod_two_rd2 + + # Compute spring force between two rods + distance_vector = surface_position_rod_two - surface_position_rod_one + np.round_(distance_vector, 12, distance_vector) + spring_force = k * (distance_vector) + + # Damping force + rod_one_element_velocity = 0.5 * ( + rod_one_velocity_collection[:, index_one] + + rod_one_velocity_collection[:, index_one + 1] + ) + rod_two_element_velocity = 0.5 * ( + rod_two_velocity_collection[:, index_two] + + rod_two_velocity_collection[:, index_two + 1] + ) + + relative_velocity = rod_two_element_velocity - rod_one_element_velocity + + distance = np.linalg.norm(distance_vector) + + # normalized_distance_vector = np.zeros((relative_velocity.shape)) + # + # idx_nonzero_distance = np.where(distance >= 1e-12)[0] + + if distance >= 1e-12: + normalized_distance_vector = distance_vector / distance + else: + normalized_distance_vector = np.zeros( + 3, + ) + + normal_relative_velocity_vector = ( + np.dot(relative_velocity, normalized_distance_vector) + * normalized_distance_vector + ) + + damping_force = -nu * normal_relative_velocity_vector + + # Compute the total force + total_force = spring_force + damping_force + + # Compute contact forces. Contact forces are applied in the case one rod penetrates to the other, in that case + # we apply a repulsive force. + center_distance = rod_two_element_position - rod_one_element_position + center_distance_unit_vec = center_distance / np.linalg.norm(center_distance) + penetration = np.linalg.norm(center_distance) - ( + rod_one_radius[index_one] + + offset_rod_one + + rod_two_radius[index_two] + + offset_rod_two + ) + + round(penetration, 12) + # Contact present only if rods penetrate to each other + if penetration < 0: + # Hertzian contact + contact_force = ( + -k_repulsive * np.abs(penetration) ** (1.5) * center_distance_unit_vec + ) + else: + contact_force = np.zeros( + 3, + ) + + # Add contact forces + total_force += contact_force + + # Re-distribute forces from elements to nodes. + rod_one_external_forces[..., index_one] += 0.5 * total_force + rod_one_external_forces[..., index_one + 1] += 0.5 * total_force + rod_two_external_forces[..., index_two] -= 0.5 * total_force + rod_two_external_forces[..., index_two + 1] -= 0.5 * total_force + + return ( + rod_one_rd2, + rod_two_rd2, + spring_force, + ) + + def apply_torques(self, rod_one, index_one, rod_two, index_two): + # pass + + self._apply_torques( + self.spring_force, + self.rod_one_rd2, + self.rod_two_rd2, + index_one, + index_two, + rod_one.director_collection, + rod_two.director_collection, + rod_one.external_torques, + rod_two.external_torques, + ) + + @staticmethod + @njit(cache=True) + def _apply_torques( + spring_force, + rod_one_rd2, + rod_two_rd2, + index_one, + index_two, + rod_one_director_collection, + rod_two_director_collection, + rod_one_external_torques, + rod_two_external_torques, + ): + # Compute torques due to the connection forces + torque_on_rod_one = np.cross(rod_one_rd2, spring_force) + torque_on_rod_two = np.cross(rod_two_rd2, -spring_force) + + torque_on_rod_one_material_frame = ( + rod_one_director_collection[:, :, index_one] @ torque_on_rod_one + ) + torque_on_rod_two_material_frame = ( + rod_two_director_collection[:, :, index_two] @ torque_on_rod_two + ) + + rod_one_external_torques[..., index_one] += torque_on_rod_one_material_frame + rod_two_external_torques[..., index_two] += torque_on_rod_two_material_frame From d5dad47f940c2b15ab9633dfe3a535386e895861 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Thu, 17 Mar 2022 00:40:31 -0500 Subject: [PATCH 10/47] enhancement:add parallel connection example. --- .../parallel_connection_example.py | 171 ++++++++++++++++++ examples/README.md | 5 + 2 files changed, 176 insertions(+) create mode 100644 examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py diff --git a/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py b/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py new file mode 100644 index 000000000..a0de0c46d --- /dev/null +++ b/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py @@ -0,0 +1,171 @@ +__doc__ = """Parallel connection example""" + +import numpy as np +import sys + +# FIXME without appending sys.path make it more generic +sys.path.append("../../../../") +from elastica import * +from elastica.experimental.connection_contact_joint.parallel_connection import get_connection_vector_straight_straight_rod, SurfaceJointSideBySide +from elastica._calculus import difference_kernel +from examples.JointCases.joint_cases_postprocessing import ( + plot_position, + plot_video, + plot_video_xy, + plot_video_xz, +) + + +class ParallelConnection( + BaseSystemCollection, Constraints, Connections, Forcing, CallBacks +): + pass + + +parallel_connection_sim = ParallelConnection() + +# setting up test params +n_elem = 10 +direction = np.array([0.0, 0.0, 1.0]) +normal = np.array([0.0, 1.0, 0.0]) +binormal = np.cross(direction, normal) +base_length = 0.2 +base_radius = 0.007 +base_area = np.pi * base_radius ** 2 +density = 1750 +nu = 1e-2 +E = 3e4 +poisson_ratio = 0.5 +shear_modulus = E / (poisson_ratio + 1.0) + +start_rod_1 = np.zeros((3,)) + 0.1 * direction +start_rod_2 = start_rod_1 + binormal * 2 * base_radius + +# Create rod 1 +rod_one = CosseratRod.straight_rod( + n_elem, + start_rod_1, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + shear_modulus=shear_modulus, +) +parallel_connection_sim.append(rod_one) +# Create rod 2 +rod_two = CosseratRod.straight_rod( + n_elem, + start_rod_2, + direction, + normal, + base_length, + base_radius, + density, + nu, + E, + shear_modulus=shear_modulus, +) +parallel_connection_sim.append(rod_two) + +# Apply boundary conditions to rod1. +parallel_connection_sim.constrain(rod_one).using( + OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,) +) + +# Apply boundary conditions to rod2. +parallel_connection_sim.constrain(rod_two).using( + OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,) +) + +# Apply a contraction force on rod one. +class ContractionForce(NoForces): + def __init__(self, ramp, force_mag, ): + self.ramp = ramp + self.force_mag = force_mag + + def apply_forces(self, system, time: np.float64 = 0.0): + # Ramp the force + factor = min(1.0, time / self.ramp) + + system.external_forces[:] -= factor * difference_kernel(self.force_mag * system.tangents) + +parallel_connection_sim.add_forcing_to(rod_one).using(ContractionForce, ramp=0.5, force_mag=1. ) + +# Connect rod 1 and rod 2 +( + rod_one_direction_vec_in_material_frame, + rod_two_direction_vec_in_material_frame, + offset_btw_rods, +) = get_connection_vector_straight_straight_rod(rod_one, rod_two) + +for i in range(n_elem): + parallel_connection_sim.connect( + first_rod=rod_one, second_rod=rod_two, first_connect_idx=i, second_connect_idx=i + ).using( + SurfaceJointSideBySide, k=1E2, nu=1E-5, k_repulsive = 1E3, rod_one_direction_vec_in_material_frame=rod_one_direction_vec_in_material_frame[:,i], + rod_two_direction_vec_in_material_frame = rod_two_direction_vec_in_material_frame[:,i], + offset_btw_rods = offset_btw_rods[i] + ) # k=kg/s2 nu=kg/s 1e-2 + +class ParallelConnecitonCallback(CallBackBaseClass): + """ + Call back function for parallel connection + """ + + 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["velocity"].append(system.velocity_collection.copy()) + return + + +pp_list_rod1 = defaultdict(list) +pp_list_rod2 = defaultdict(list) + + +parallel_connection_sim.collect_diagnostics(rod_one).using( + ParallelConnecitonCallback, step_skip=1000, callback_params=pp_list_rod1 +) +parallel_connection_sim.collect_diagnostics(rod_two).using( + ParallelConnecitonCallback, step_skip=1000, callback_params=pp_list_rod2 +) + + +parallel_connection_sim.finalize() +timestepper = PositionVerlet() + +final_time = 5.0 +dl = base_length / n_elem +dt = 1e-5 +total_steps = int(final_time / dt) +print("Total steps", total_steps) +integrate(timestepper, parallel_connection_sim, final_time, total_steps) + +PLOT_FIGURE = True +SAVE_FIGURE = False +PLOT_VIDEO = True + +# plotting results +if PLOT_FIGURE: + filename = "parallel_connection_test_last_node_pos_xy.png" + plot_position(pp_list_rod1, pp_list_rod2, filename, SAVE_FIGURE) + +if PLOT_VIDEO: + filename = "parallel_connection_test.mp4" + plot_video(pp_list_rod1, pp_list_rod2, video_name=filename, margin=0.2, fps=100) + plot_video_xy( + pp_list_rod1, pp_list_rod2, video_name=filename + "_xy.mp4", margin=0.2, fps=100 + ) + plot_video_xz( + pp_list_rod1, pp_list_rod2, video_name=filename + "_xz.mp4", margin=0.2, fps=100 + ) diff --git a/examples/README.md b/examples/README.md index 0c3b67e4c..20ab95948 100644 --- a/examples/README.md +++ b/examples/README.md @@ -67,3 +67,8 @@ Examples can serve as a starting template for customized usages. * [Elastica RL control](https://github.com/GazzolaLab/Elastica-RL-control) - Case presented in [Elastica: A compliant mechanics environment for soft robotic control](https://doi.org/10.1109/LRA.2021.3063698) * [Gym Softrobot](https://github.com/skim0119/gym-softrobot) - Soft-robot control environment developed in OpenAI-gym format to study slender body control with reinforcement learning. + +## Experimental Cases +* [ParallelConnectionExample](./ExperimentalCases/ParallelConnectionExample) + * __Purpose__: Demonstrate the usage of parallel connection. + * __Features__: connect two parallel rods \ No newline at end of file From aa366a64986b4f599a93d71740bef82b07837f3c Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sun, 20 Mar 2022 03:03:46 -0500 Subject: [PATCH 11/47] Update conf.py Fix typo --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index 831462830..905e3d0c7 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -81,7 +81,7 @@ # a list of builtin themes. html_theme = 'sphinx_book_theme' -html_theme_option = { +html_theme_options = { "repository_url": "https://github.com/GazzolaLab/PyElastica", "use_repository_button": True, } From 8c7314fc28dc9d0795d71bd614b0393ae7d475a8 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sun, 20 Mar 2022 23:57:40 -0500 Subject: [PATCH 12/47] wip: knot theory unittest --- tests/test_rod/test_knot_theory.py | 29 +++++++++++++++++++++++------ tests/test_rod/test_rods.py | 12 ++++++++---- 2 files changed, 31 insertions(+), 10 deletions(-) diff --git a/tests/test_rod/test_knot_theory.py b/tests/test_rod/test_knot_theory.py index abdca6b7b..be7361e8d 100644 --- a/tests/test_rod/test_knot_theory.py +++ b/tests/test_rod/test_knot_theory.py @@ -14,16 +14,21 @@ from test_rod.test_rods import MockTestRod from elastica.rod.rod_base import RodBase -from elastica.rod.knot_theory import KnotTheory -class TestRodWithKnotTheory(MockTestRod, KnotTheory): - def __init__(self): - super().__init__() - self.radius = np.random.randn(MaxDimension.value(), self.n_elem) +@pytest.fixture +def knot_theory(): + from elastica.rod import knot_theory + return knot_theory + + +def test_knot_theory_mixin_methods(knot_theory): + class TestRodWithKnotTheory(MockTestRod, knot_theory.KnotTheory): + def __init__(self): + super().__init__() + self.radius = np.random.randn(MaxDimension.value(), self.n_elem) -def test_knot_theory_mixin_methods(): rod = TestRodWithKnotTheory() assert hasattr( rod, "MIXIN_PROTOCOL" @@ -37,3 +42,15 @@ def test_knot_theory_mixin_methods(): assert hasattr( rod, "compute_link" ), "Expected to mix-in functionals into the rod class: compute_link" + + +def test_knot_theory_mixin_methods_with_no_radius(knot_theory): + class TestRodWithKnotTheoryWithoutRadius(MockTestRod, knot_theory.KnotTheory): + def __init__(self): + super().__init__() + + rod = TestRodWithKnotTheoryWithoutRadius() + with pytest.raises(AttributeError) as e_info: + rod.compute_writhe() + with pytest.raises(AttributeError) as e_info: + rod.compute_link() diff --git a/tests/test_rod/test_rods.py b/tests/test_rod/test_rods.py index 9577a3b83..bb1a82544 100644 --- a/tests/test_rod/test_rods.py +++ b/tests/test_rod/test_rods.py @@ -17,14 +17,18 @@ class MockTestRod: def __init__(self): self.n_elem = 32 - self.position_collection = np.random.randn(MaxDimension.value(), self.n_elem) + self.position_collection = np.random.randn( + MaxDimension.value(), self.n_elem + 1 + ) self.director_collection = np.random.randn( MaxDimension.value(), MaxDimension.value(), self.n_elem ) - self.velocity_collection = np.random.randn(MaxDimension.value(), self.n_elem) + self.velocity_collection = np.random.randn( + MaxDimension.value(), self.n_elem + 1 + ) self.omega_collection = np.random.randn(MaxDimension.value(), self.n_elem) - self.mass = np.abs(np.random.randn(self.n_elem)) - self.external_forces = np.zeros(self.n_elem) + self.mass = np.abs(np.random.randn(self.n_elem + 1)) + self.external_forces = np.zeros(self.n_elem + 1) # Choosing 15 and 31 as nelems to reflect common expected From 6b6c2c13825606f20750d811c8702cf952aaed54 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 21 Mar 2022 00:09:45 -0500 Subject: [PATCH 13/47] update: flake8 exclude doc configuration --- .flake8 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/.flake8 b/.flake8 index 34b2930b3..7f925ebc7 100644 --- a/.flake8 +++ b/.flake8 @@ -3,4 +3,8 @@ ignore = E203, E266, E501, W503, F403, F401, W191 max-line-length = 88 max-complexity = 18 select = B,C,E,F,W,T4,B9 -exclude = .git,__pycache__,elastica/rod.py \ No newline at end of file +exclude = + .git, + __pycache__, + docs/conf.py + tests From fff6804d67b9d6006e105c29cc319a0f18235371 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Tue, 22 Mar 2022 12:28:46 -0500 Subject: [PATCH 14/47] update: knot_theory inline documentation --- docs/api/rods.rst | 17 +++++++- docs/conf.py | 3 ++ elastica/rod/knot_theory.py | 78 +++++++++++++++++++++++++++---------- 3 files changed, 75 insertions(+), 23 deletions(-) diff --git a/docs/api/rods.rst b/docs/api/rods.rst index c3248118c..628818898 100644 --- a/docs/api/rods.rst +++ b/docs/api/rods.rst @@ -1,12 +1,12 @@ Rods -===== +==== .. automodule:: elastica.rod.rod_base :members: :exclude-members: __weakref__ Cosserat Rod -~~~~~~~~~~~~ +------------ +------------+-------------------+----------------------------------------+-----------------------------+ | | On Nodes (+1) | On Elements (n_elements) | On Voronoi (-1) | @@ -39,9 +39,22 @@ Cosserat Rod .. automodule:: elastica.rod.cosserat_rod :exclude-members: __weakref__, __init__, update_accelerations, zeroed_out_external_forces_and_torques, compute_internal_forces_and_torques :members: + :inherited-members: .. Constitutive Models .. ~~~~~~~~~~~~~~~~~~~ .. .. automodule:: elastica.rod.constitutive_model .. :members: .. :exclude-members: __weakref__ + + +Knot Theory (Mixin) +~~~~~~~~~~~~~~~~~~~ + +.. autoclass:: elastica.rod.knot_theory.KnotTheory + +.. autoclass:: elastica.rod.knot_theory.KnotTheoryCompatibleProtocol + +.. automodule:: elastica.rod.knot_theory + :exclude-members: __weakref__, __init__, KnotTheory, KnotTheoryCompatibleProtocol + :members: diff --git a/docs/conf.py b/docs/conf.py index e6662bcf2..964ad0a03 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -90,5 +90,8 @@ html_static_path = ['_static'] html_css_files = ['css/*', 'css/logo.css'] +# -- Options for autodoc --------------------------------------------------- +autodoc_member_order = 'bysource' + # -- Options for numpydoc --------------------------------------------------- numpydoc_show_class_members = False diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index e0de0fd65..d55935c7c 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -53,11 +53,28 @@ def base_length(self) -> np.ndarray: class KnotTheory: - # This mixin should be used in RodBase-derived class that satisfies KnotCompatibleProtocol + """ + This mixin should be used in RodBase-derived class that satisfies KnotCompatibleProtocol. + The theory behind this module is based on the method from Klenin & Langowski 2000 paper. + + KnotTheory can be mixed with any rod-class based on RodBase:: + + class MyRod(RodBase, KnotTheory): + def __init__(self): + super().__init__() + rod = MyRod(...) + + total_twist = rod.compute_twist() + total_link = rod.compute_link() + + """ + MIXIN_PROTOCOL = Union[RodBase, KnotTheoryCompatibleProtocol] def compute_twist(self: MIXIN_PROTOCOL): - """compute_twist""" + """ + See :ref:`api/rods:Knot Theory (Mixin)` for the detail. + """ total_twist, local_twist = compute_twist( self.position_collection[None, ...], self.director_collection[0][None, ...], @@ -67,7 +84,8 @@ def compute_twist(self: MIXIN_PROTOCOL): def compute_writhe( self: MIXIN_PROTOCOL, type_of_additional_segment: str = "next_tangent" ): - """compute_writhe + """ + See :ref:`api/rods:Knot Theory (Mixin)` for the detail. Parameters ---------- @@ -84,7 +102,8 @@ def compute_writhe( def compute_link( self: MIXIN_PROTOCOL, type_of_additional_segment: str = "next_tangent" ): - """compute_link + """ + See :ref:`api/rods:Knot Theory (Mixin)` for the detail. Parameters ---------- @@ -106,6 +125,11 @@ def compute_twist(center_line, normal_collection): Compute the twist of a rod, using center_line and normal collection. Methods used in this function is adapted from method 2a Klenin & Langowski 2000 paper. + Warnings + -------- + If center line is straight, although the normals of each element is pointing different direction computed twist + will be zero. + Parameters ---------- center_line : numpy.ndarray @@ -115,13 +139,10 @@ def compute_twist(center_line, normal_collection): 3D (time, 3, n_elems) array containing data with 'float' type. Time history of rod elements normal direction. - Warnings - -------- - If center line is straight, although the normals of each element is pointing different direction computed twist - will be zero. - Returns ------- + total_twist : numpy.ndarray + local_twist : numpy.ndarray """ @@ -133,9 +154,6 @@ def compute_twist(center_line, normal_collection): @njit(cache=True) def _compute_twist(center_line, normal_collection): """ - Compute the twist of a rod, using center_line and normal collection. Methods used in this function is - adapted from method 2a Klenin & Langowski 2000 paper. - Parameters ---------- center_line : numpy.ndarray @@ -147,7 +165,8 @@ def _compute_twist(center_line, normal_collection): Returns ------- - + total_twist : numpy.ndarray + local_twist : numpy.ndarray """ timesize, _, blocksize = center_line.shape @@ -211,6 +230,7 @@ def _compute_twist(center_line, normal_collection): def compute_writhe(center_line, segment_length, type_of_additional_segment): """ This function computes the total writhe history of a rod. + Equations used are from method 1a from Klenin & Langowski 2000 paper. Parameters ---------- @@ -225,7 +245,7 @@ def compute_writhe(center_line, segment_length, type_of_additional_segment): Returns ------- - total_writhe + total_writhe : numpy.ndarray """ @@ -244,9 +264,6 @@ def compute_writhe(center_line, segment_length, type_of_additional_segment): @njit(cache=True) def _compute_writhe(center_line): """ - This function computes the total writhe history of a rod. - Equations used are from method 1a from Klenin & Langowski 2000 paper. - Parameters ---------- center_line : numpy.ndarray @@ -255,7 +272,7 @@ def _compute_writhe(center_line): Returns ------- - + total_writhe : numpy.ndarray """ time, _, blocksize = center_line.shape @@ -324,7 +341,8 @@ def compute_link( type_of_additional_segment: str, ): """ - Compute the link of a rod. + This function computes the total link history of a rod. + Equations used are from method 1a from Klenin & Langowski 2000 paper. Parameters ---------- @@ -342,6 +360,11 @@ def compute_link( type_of_additional_segment : str Determines the method to compute new segments (elements) added to the rod. Valid inputs are "next_tangent", "end_to_end", "net_tangent", otherwise program uses the center line. + + Returns + ------- + total_link : numpy.ndarray + """ # Compute auxiliary line @@ -388,6 +411,7 @@ def _compute_auxiliary_line(center_line, normal_collection, radius): Returns ------- + auxillary_line : numpy.ndarray """ time, _, blocksize = center_line.shape @@ -426,8 +450,6 @@ def _compute_auxiliary_line(center_line, normal_collection, radius): @njit(cache=True) def _compute_link(center_line, auxiliary_line): """ - This function computes the total link history of a rod. - Equations used are from method 1a from Klenin & Langowski 2000 paper. Parameters ---------- @@ -440,6 +462,7 @@ def _compute_link(center_line, auxiliary_line): Returns ------- + total_link : numpy.ndarray """ timesize, _, blocksize_center_line = center_line.shape @@ -525,6 +548,7 @@ def _compute_auxiliary_line_added_segments( Returns ------- + new_auxiliary_line : numpy.ndarray """ timesize, _, blocksize = auxiliary_line.shape @@ -563,6 +587,9 @@ def compute_additional_segment(center_line, segment_length, type_of_additional_s Returns ------- + center_line : numpy.ndarray + beginning_direction : numpy.ndarray + end_direction : numpy.ndarray """ @@ -608,6 +635,9 @@ def _compute_additional_segment_using_next_tangent(center_line, segment_length): Returns ------- + center_line : numpy.ndarray + beginning_direction : numpy.ndarray + end_direction : numpy.ndarray """ @@ -655,6 +685,9 @@ def _compute_additional_segment_using_end_to_end(center_line, segment_length): Returns ------- + center_line : numpy.ndarray + beginning_direction : numpy.ndarray + end_direction : numpy.ndarray """ timesize, _, blocksize = center_line.shape @@ -702,6 +735,9 @@ def _compute_additional_segment_using_net_tangent(center_line, segment_length): Returns ------- + center_line : numpy.ndarray + beginning_direction : numpy.ndarray + end_direction : numpy.ndarray """ timesize, _, blocksize = center_line.shape From 4634b14f4ff7a6dcb4f3d2db932161cf8ca0572e Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Tue, 22 Mar 2022 12:32:32 -0500 Subject: [PATCH 15/47] update: allow free-functions to be accesible --- elastica/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/elastica/__init__.py b/elastica/__init__.py index b9a976d14..98bb897f4 100644 --- a/elastica/__init__.py +++ b/elastica/__init__.py @@ -4,6 +4,7 @@ from collections import defaultdict from elastica.wrappers import * from elastica.rod.cosserat_rod import * +from elastica.rod.knot_theory import * from elastica.rigidbody import * from elastica.boundary_conditions import * from elastica.external_forces import * From dffce47bb8aebb93cc4fb23396b30fb42e80eb21 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Tue, 22 Mar 2022 23:33:54 -0500 Subject: [PATCH 16/47] update: reduce codebase, update documentation --- elastica/rod/knot_theory.py | 236 +++++++++++------------------------- 1 file changed, 74 insertions(+), 162 deletions(-) diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index d55935c7c..4f0508b35 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -67,6 +67,27 @@ def __init__(self): total_twist = rod.compute_twist() total_link = rod.compute_link() + There are few alternative way of handling edge-condition in computing Link and Writhe. + Here, we provide three methods: "next_tangent", "end_to_end", and "net_tangent". + The default *type_of_additional_segment* is set to "next_tangent." + + ========================== ===================================== + type_of_additional_segment Description + ========================== ===================================== + next_tangent | Adds a two new point at the begining and end of the center line. + | Distance of these points are given in segment_length. + | Direction of these points are computed using the rod tangents at + | the begining and end. + end_to_end | Adds a two new point at the begining and end of the center line. + | Distance of these points are given in segment_length. + | Direction of these points are computed using the rod node end + | positions. + net_tangent | Adds a two new point at the begining and end of the center line. + | Distance of these points are given in segment_length. Direction of + | these points are point wise avarege of nodes at the first and + | second half of the rod. + ========================== ===================================== + """ MIXIN_PROTOCOL = Union[RodBase, KnotTheoryCompatibleProtocol] @@ -249,7 +270,7 @@ def compute_writhe(center_line, segment_length, type_of_additional_segment): """ - center_line_with_added_segments, _, _ = compute_additional_segment( + center_line_with_added_segments, _, _ = _compute_additional_segment( center_line, segment_length, type_of_additional_segment ) @@ -375,7 +396,7 @@ def compute_link( center_line_with_added_segments, beginning_direction, end_direction, - ) = compute_additional_segment( + ) = _compute_additional_segment( center_line, segment_length, type_of_additional_segment ) auxiliary_line_with_added_segments = _compute_auxiliary_line_added_segments( @@ -568,12 +589,26 @@ def _compute_auxiliary_line_added_segments( return new_auxiliary_line -def compute_additional_segment(center_line, segment_length, type_of_additional_segment): +@njit(cache=True) +def _compute_additional_segment( + center_line, segment_length, type_of_additional_segment +): """ This function adds two points at the end of center line. Distance from the center line is given by segment_length. Direction from center line to the new point locations can be computed using 3 methods, which can be selected by type. For more details section S2 of Charles et. al. PRL 2019 paper. + next_tangent: + This function adds a two new point at the begining and end of the center line. Distance of these points are + given in segment_length. Direction of these points are computed using the rod tangents at the begining and end. + end_to_end: + This function adds a two new point at the begining and end of the center line. Distance of these points are + given in segment_length. Direction of these points are computed using the rod node end positions. + net_tangent: + This function adds a two new point at the begining and end of the center line. Distance of these points are + given in segment_length. Direction of these points are point wise avarege of nodes at the first and second half + of the rod. + Parameters ---------- center_line : numpy.ndarray @@ -583,7 +618,7 @@ def compute_additional_segment(center_line, segment_length, type_of_additional_s Length of added segments. type_of_additional_segment : str Determines the method to compute new segments (elements) added to the rod. - Valid inputs are "next_tangent", "end_to_end", "net_tangent", otherwise it returns the center line. + Valid inputs are "next_tangent", "end_to_end", "net_tangent". If None, returns the center line. Returns ------- @@ -592,104 +627,11 @@ def compute_additional_segment(center_line, segment_length, type_of_additional_s end_direction : numpy.ndarray """ - - if type_of_additional_segment == "next_tangent": - ( - center_line, - beginning_direction, - end_direction, - ) = _compute_additional_segment_using_next_tangent(center_line, segment_length) - elif type_of_additional_segment == "end_to_end": - ( - center_line, - beginning_direction, - end_direction, - ) = _compute_additional_segment_using_end_to_end(center_line, segment_length) - elif type_of_additional_segment == "net_tangent": - ( - center_line, - beginning_direction, - end_direction, - ) = _compute_additional_segment_using_net_tangent(center_line, segment_length) - else: - print("Additional segments are not used") + if type_of_additional_segment is None: beginning_direction = np.zeros((center_line.shape[0], 3)) end_direction = np.zeros((center_line.shape[0], 3)) + return center_line, beginning_direction, end_direction - return center_line, beginning_direction, end_direction - - -@njit(cache=True) -def _compute_additional_segment_using_next_tangent(center_line, segment_length): - """ - This function adds a two new point at the begining and end of the center line. Distance of these points are - given in segment_length. Direction of these points are computed using the rod tangents at the begining and end. - - Parameters - ---------- - center_line : numpy.ndarray - 3D (time, 3, n_nodes) array containing data with 'float' type. - Time history of rod node positions. - segment_length : float - Length of added segments. - - Returns - ------- - center_line : numpy.ndarray - beginning_direction : numpy.ndarray - end_direction : numpy.ndarray - - """ - - timesize, _, blocksize = center_line.shape - new_center_line = np.zeros( - (timesize, 3, blocksize + 2) - ) # +2 is for added two new points - beginning_direction = np.zeros((timesize, 3)) - end_direction = np.zeros((timesize, 3)) - - for i in range(timesize): - # Direction of the additional point at the beginning of the rod - direction_of_rod_begin = center_line[i, :, 0] - center_line[i, :, 1] - direction_of_rod_begin /= np.linalg.norm(direction_of_rod_begin) - - # Direction of the additional point at the end of the rod - direction_of_rod_end = center_line[i, :, -1] - center_line[i, :, -2] - direction_of_rod_end /= np.linalg.norm(direction_of_rod_end) - - first_point = center_line[i, :, 0] + segment_length * direction_of_rod_begin - last_point = center_line[i, :, -1] + segment_length * direction_of_rod_end - - new_center_line[i, :, 1:-1] = center_line[i, :, :] - new_center_line[i, :, 0] = first_point - new_center_line[i, :, -1] = last_point - beginning_direction[i, :] = direction_of_rod_begin - end_direction[i, :] = direction_of_rod_end - - return new_center_line, beginning_direction, end_direction - - -@njit(cache=True) -def _compute_additional_segment_using_end_to_end(center_line, segment_length): - """ - This function adds a two new point at the begining and end of the center line. Distance of these points are - given in segment_length. Direction of these points are computed using the rod node end positions. - - Parameters - ---------- - center_line : numpy.ndarray - 3D (time, 3, n_nodes) array containing data with 'float' type. - Time history of rod node positions. - segment_length : float - Length of added segments. - - Returns - ------- - center_line : numpy.ndarray - beginning_direction : numpy.ndarray - end_direction : numpy.ndarray - - """ timesize, _, blocksize = center_line.shape new_center_line = np.zeros( (timesize, 3, blocksize + 2) @@ -697,72 +639,42 @@ def _compute_additional_segment_using_end_to_end(center_line, segment_length): beginning_direction = np.zeros((timesize, 3)) end_direction = np.zeros((timesize, 3)) - for i in range(timesize): - # Direction of the additional point at the beginning of the rod - direction_of_rod_begin = center_line[i, :, 0] - center_line[i, :, -1] - direction_of_rod_begin /= np.linalg.norm(direction_of_rod_begin) - - # Direction of the additional point at the end of the rod - direction_of_rod_end = -direction_of_rod_begin - - first_point = center_line[i, :, 0] + segment_length * direction_of_rod_begin - last_point = center_line[i, :, -1] + segment_length * direction_of_rod_end - - new_center_line[i, :, 1:-1] = center_line[i, :, :] - new_center_line[i, :, 0] = first_point - new_center_line[i, :, -1] = last_point - - beginning_direction[i, :] = direction_of_rod_begin - end_direction[i, :] = direction_of_rod_end - - return new_center_line, beginning_direction, end_direction - - -@njit(cache=True) -def _compute_additional_segment_using_net_tangent(center_line, segment_length): - """ - This function adds a two new point at the begining and end of the center line. Distance of these points are - given in segment_length. Direction of these points are point wise avarege of nodes at the first and second half - of the rod. - - Parameters - ---------- - center_line : numpy.ndarray - 3D (time, 3, n_nodes) array containing data with 'float' type. - Time history of rod node positions. - segment_length : float - Length of added segments. - - Returns - ------- - center_line : numpy.ndarray - beginning_direction : numpy.ndarray - end_direction : numpy.ndarray + if type_of_additional_segment == "next_tangent": + for i in range(timesize): + # Direction of the additional point at the beginning of the rod + direction_of_rod_begin = center_line[i, :, 0] - center_line[i, :, 1] + direction_of_rod_begin /= np.linalg.norm(direction_of_rod_begin) + + # Direction of the additional point at the end of the rod + direction_of_rod_end = center_line[i, :, -1] - center_line[i, :, -2] + direction_of_rod_end /= np.linalg.norm(direction_of_rod_end) + elif type_of_additional_segment == "end_to_end": + for i in range(timesize): + # Direction of the additional point at the beginning of the rod + direction_of_rod_begin = center_line[i, :, 0] - center_line[i, :, -1] + direction_of_rod_begin /= np.linalg.norm(direction_of_rod_begin) - """ - timesize, _, blocksize = center_line.shape - new_center_line = np.zeros( - (timesize, 3, blocksize + 2) - ) # +2 is for added two new points - beginning_direction = np.zeros((timesize, 3)) - end_direction = np.zeros((timesize, 3)) + # Direction of the additional point at the end of the rod + direction_of_rod_end = -direction_of_rod_begin + elif type_of_additional_segment == "net_tangent": + for i in range(timesize): + # Direction of the additional point at the beginning of the rod + n_nodes_begin = int(np.floor(blocksize / 2)) + average_begin = ( + np.sum(center_line[i, :, :n_nodes_begin], axis=1) / n_nodes_begin + ) + n_nodes_end = int(np.ceil(blocksize / 2)) + average_end = np.sum(center_line[i, :, n_nodes_end:], axis=1) / ( + blocksize - n_nodes_end + 1 + ) + direction_of_rod_begin = average_begin - average_end + direction_of_rod_begin /= np.linalg.norm(direction_of_rod_begin) + direction_of_rod_end = -direction_of_rod_begin + else: + raise NotImplementedError("unavailable type_of_additional_segment is given") + # Compute new centerline and beginning/end direction for i in range(timesize): - # Direction of the additional point at the beginning of the rod - n_nodes_begin = int(np.floor(blocksize / 2)) - average_begin = ( - np.sum(center_line[i, :, :n_nodes_begin], axis=1) / n_nodes_begin - ) - n_nodes_end = int(np.ceil(blocksize / 2)) - average_end = np.sum(center_line[i, :, n_nodes_end:], axis=1) / ( - blocksize - n_nodes_end + 1 - ) - - direction_of_rod_begin = average_begin - average_end - direction_of_rod_begin /= np.linalg.norm(direction_of_rod_begin) - - direction_of_rod_end = -direction_of_rod_begin - first_point = center_line[i, :, 0] + segment_length * direction_of_rod_begin last_point = center_line[i, :, -1] + segment_length * direction_of_rod_end From 7c3bc5719a5dd9c88c438d486d26358b1ff64dd8 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Tue, 22 Mar 2022 23:52:21 -0500 Subject: [PATCH 17/47] wip: test case for compute_additional_segment --- tests/test_rod/test_knot_theory.py | 49 +++++++++++++++++++++++++++++- 1 file changed, 48 insertions(+), 1 deletion(-) diff --git a/tests/test_rod/test_knot_theory.py b/tests/test_rod/test_knot_theory.py index be7361e8d..eeeaf3dd5 100644 --- a/tests/test_rod/test_knot_theory.py +++ b/tests/test_rod/test_knot_theory.py @@ -11,9 +11,10 @@ ) from elastica.utils import MaxDimension -from test_rod.test_rods import MockTestRod +from test_rods import MockTestRod from elastica.rod.rod_base import RodBase +from elastica.rod.knot_theory import _compute_additional_segment @pytest.fixture @@ -54,3 +55,49 @@ def __init__(self): rod.compute_writhe() with pytest.raises(AttributeError) as e_info: rod.compute_link() + + +@pytest.mark.parametrize("timesteps", [1, 5, 10]) +@pytest.mark.parametrize("n_elem", [1, 3, 8]) +@pytest.mark.parametrize("segment_length", [1.0, 10.0, 100.0]) +def test_knot_theory_compute_additional_segment_next_tangent_case( + timesteps, n_elem, segment_length +): + ... + + +@pytest.mark.parametrize("timesteps", [1, 5, 10]) +@pytest.mark.parametrize("n_elem", [1, 3, 8]) +@pytest.mark.parametrize("segment_length", [1.0, 10.0, 100.0]) +def test_knot_theory_compute_additional_segment_end_to_end_case( + timesteps, n_elem, segment_length +): + ... + + +@pytest.mark.parametrize("timesteps", [1, 5, 10]) +@pytest.mark.parametrize("n_elem", [1, 3, 8]) +@pytest.mark.parametrize("segment_length", [1.0, 10.0, 100.0]) +def test_knot_theory_compute_additional_segment_net_tangent_case( + timesteps, n_elem, segment_length +): + ... + + +@pytest.mark.parametrize("timesteps", [1, 5, 10]) +@pytest.mark.parametrize("n_elem", [1, 3, 8]) +@pytest.mark.parametrize("segment_length", [1.0, 10.0, 100.0]) +def test_knot_theory_compute_additional_segment_none_case( + timesteps, n_elem, segment_length +): + center_line = np.random.random([timesteps, 3, n_elem]) + new_center_line, beginning_direction, end_direction = _compute_additional_segment( + center_line, segment_length, None + ) + + assert_allclose(new_center_line, center_line) + assert_allclose(beginning_direction, 0.0) + assert_allclose(end_direction, 0.0) + assert_allclose(new_center_line.shape, [timesteps, 3, n_elem]) + assert_allclose(beginning_direction.shape[0], timesteps) + assert_allclose(end_direction.shape[0], timesteps) From b41f6ab1a48983e82784ba60528eb89cbc2258a7 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 23 Mar 2022 00:28:42 -0500 Subject: [PATCH 18/47] wip: complete additional_segment test cases --- tests/test_rod/test_knot_theory.py | 101 +++++++++++++++++++++++------ 1 file changed, 81 insertions(+), 20 deletions(-) diff --git a/tests/test_rod/test_knot_theory.py b/tests/test_rod/test_knot_theory.py index eeeaf3dd5..9683e2757 100644 --- a/tests/test_rod/test_knot_theory.py +++ b/tests/test_rod/test_knot_theory.py @@ -14,7 +14,10 @@ from test_rods import MockTestRod from elastica.rod.rod_base import RodBase -from elastica.rod.knot_theory import _compute_additional_segment +from elastica.rod.knot_theory import ( + KnotTheoryCompatibleProtocol, + _compute_additional_segment, +) @pytest.fixture @@ -24,6 +27,13 @@ def knot_theory(): return knot_theory +def test_knot_theory_protocol(): + # To clear the protocol test coverage + with pytest.raises(TypeError) as e_info: + protocol = KnotTheoryCompatibleProtocol() + assert "cannot be instantiated" in e_info + + def test_knot_theory_mixin_methods(knot_theory): class TestRodWithKnotTheory(MockTestRod, knot_theory.KnotTheory): def __init__(self): @@ -57,31 +67,82 @@ def __init__(self): rod.compute_link() -@pytest.mark.parametrize("timesteps", [1, 5, 10]) -@pytest.mark.parametrize("n_elem", [1, 3, 8]) -@pytest.mark.parametrize("segment_length", [1.0, 10.0, 100.0]) -def test_knot_theory_compute_additional_segment_next_tangent_case( - timesteps, n_elem, segment_length -): - ... +@pytest.mark.parametrize("type_str", ["randomstr1", "nextnext_tangent", " "]) +def test_knot_theory_compute_additional_segment_integrity(type_str): + center_line = np.zeros([1, 3, 10]) + center_line[:, 2, :] = np.arange(10) + with pytest.raises(NotImplementedError) as e_info: + _compute_additional_segment(center_line, 10.0, type_str) -@pytest.mark.parametrize("timesteps", [1, 5, 10]) -@pytest.mark.parametrize("n_elem", [1, 3, 8]) +@pytest.mark.parametrize("n_elem", [2, 3, 8]) @pytest.mark.parametrize("segment_length", [1.0, 10.0, 100.0]) -def test_knot_theory_compute_additional_segment_end_to_end_case( - timesteps, n_elem, segment_length +@pytest.mark.parametrize("type_str", ["next_tangent", "end_to_end", "net_tangent"]) +def test_knot_theory_compute_additional_segment_straight_case( + n_elem, segment_length, type_str ): - ... + # If straight rod give, result should be same regardless of type + center_line = np.zeros([1, 3, n_elem]) + center_line[0, 2, :] = np.linspace(0, 5, n_elem) + ncl, bd, ed = _compute_additional_segment(center_line, segment_length, type_str) + assert_allclose(ncl[0, :, 0], np.array([0, 0, -segment_length])) + assert_allclose( + ncl[0, :, -1], np.array([0, 0, center_line[0, 2, -1] + segment_length]) + ) + assert_allclose(bd[0], np.array([0, 0, -1])) + assert_allclose(ed[0], np.array([0, 0, 1])) -@pytest.mark.parametrize("timesteps", [1, 5, 10]) -@pytest.mark.parametrize("n_elem", [1, 3, 8]) -@pytest.mark.parametrize("segment_length", [1.0, 10.0, 100.0]) -def test_knot_theory_compute_additional_segment_net_tangent_case( - timesteps, n_elem, segment_length -): - ... +def test_knot_theory_compute_additional_segment_next_tangent_case(): + center_line = np.array( + [[0, 0, 0], [0, 0, 1], [0, 1, 1], [0, 1, 0]], dtype=np.float64 + ).T[None, ...] + segment_length = 10 + ncl, bd, ed = _compute_additional_segment( + center_line, segment_length, "next_tangent" + ) + # fmt: off + assert_allclose(ncl[0], + np.array([[ 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 1., 1., 1.], + [-10., 0., 1., 1., 0., -10.]])) + assert_allclose(bd[0], np.array([ 0., 0., -1.])) + assert_allclose(ed[0], np.array([ 0., 0., -1.])) + # fmt: on + + +def test_knot_theory_compute_additional_segment_end_to_end_case(): + center_line = np.array( + [[0, 0, 0], [0, 0, 1], [0, 1, 1], [0, 1, 0]], dtype=np.float64 + ).T[None, ...] + segment_length = 10 + ncl, bd, ed = _compute_additional_segment(center_line, segment_length, "end_to_end") + # fmt: off + assert_allclose(ncl[0], + np.array([[ 0., 0., 0., 0., 0., 0.], + [-10., 0., 0., 1., 1., 11.], + [ 0., 0., 1., 1., 0., 0.]])) + assert_allclose(bd[0], np.array([ 0., -1., 0.])) + assert_allclose(ed[0], np.array([-0., 1., -0.])) + # fmt: on + + +def test_knot_theory_compute_additional_segment_net_tangent_case(): + center_line = np.array( + [[0, 0, 0], [0, 0, 1], [0, 1, 1], [0, 1, 0]], dtype=np.float64 + ).T[None, ...] + segment_length = 10 + ncl, bd, ed = _compute_additional_segment( + center_line, segment_length, "net_tangent" + ) + # fmt: off + assert_allclose(ncl[0], + np.array([[ 0. , 0. , 0. , 0. , 0. , 0. ], + [-9.701425 , 0. , 0. , 1. , 1. , 10.701425 ], + [ 2.42535625, 0. , 1. , 1. , 0. , -2.42535625]])) + assert_allclose(bd[0], np.array([ 0. , -0.9701425 , 0.24253563])) + assert_allclose(ed[0], np.array([-0. , 0.9701425 , -0.24253563])) + # fmt: on @pytest.mark.parametrize("timesteps", [1, 5, 10]) From fb657e6680d0fa7eac45023b397f65fadbf4586a Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 23 Mar 2022 00:50:56 -0500 Subject: [PATCH 19/47] update: add assertion to clearify the error message --- elastica/rod/knot_theory.py | 34 +++++++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index 4f0508b35..c160728a0 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -166,6 +166,15 @@ def compute_twist(center_line, normal_collection): local_twist : numpy.ndarray """ + # fmt: off + # Format is turned off because I want the assertion message to display the line. + assert center_line.shape[2] == normal_collection.shape[2] + 1, \ + "Please check the shape (axis-2) of center_line(n_node=n_elems+1) or normal_collection(n_elems)." + assert center_line.shape[0] == normal_collection.shape[0], \ + "The number of timesteps (axis-0) must be equal" + assert center_line.shape[1] == normal_collection.shape[1] == 3, \ + "The dimension (axis-1) must be 3" + # fmt: on total_twist, local_twist = _compute_twist(center_line, normal_collection) @@ -269,14 +278,22 @@ def compute_writhe(center_line, segment_length, type_of_additional_segment): total_writhe : numpy.ndarray """ + # fmt: off + # Format is turned off because I want the assertion message to display the line. + assert type(segment_length) is float, \ + "segment_length is not a float. (not numpy.float)" + assert center_line.shape[2] == normal_collection.shape[2] + 1, \ + "Please check the shape (axis-2) of center_line(n_node=n_elems+1) or normal_collection(n_elems)." + assert center_line.shape[0] == normal_collection.shape[0], \ + "The number of timesteps (axis-0) must be equal" + assert center_line.shape[1] == normal_collection.shape[1] == 3, \ + "The dimension (axis-1) must be 3" + # fmt: on center_line_with_added_segments, _, _ = _compute_additional_segment( center_line, segment_length, type_of_additional_segment ) - """ - Total writhe of a rod is computed using the method 1a from Klenin & Langowski 2000 - """ total_writhe = _compute_writhe(center_line_with_added_segments) return total_writhe @@ -387,6 +404,17 @@ def compute_link( total_link : numpy.ndarray """ + # fmt: off + # Format is turned off because I want the assertion message to display the line. + assert type(segment_length) is float, \ + "segment_length is not a float. (not numpy.float)" + assert center_line.shape[2] == normal_collection.shape[2] + 1 == radius.shape[2] + 1, \ + "Please check the shape (axis-2) of center_line(n_node=n_elems+1) or normal_collection(n_elems)." + assert center_line.shape[0] == normal_collection.shape[0] == radius.shape[0], \ + "The number of timesteps (axis-0) must be equal" + assert center_line.shape[1] == normal_collection.shape[1] == 3, \ + "The dimension (axis-1) for center_line and normal_collection must be 3" + # fmt: on # Compute auxiliary line auxiliary_line = _compute_auxiliary_line(center_line, normal_collection, radius) From 6ba7b3a7a71d3e05fa0e0acccc0b45162defa634 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 23 Mar 2022 01:16:48 -0500 Subject: [PATCH 20/47] wip: arithmetic test for LWT --- elastica/rod/knot_theory.py | 12 ++---- tests/test_rod/test_knot_theory.py | 59 ++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 8 deletions(-) diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index c160728a0..a2b5606a6 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -282,11 +282,7 @@ def compute_writhe(center_line, segment_length, type_of_additional_segment): # Format is turned off because I want the assertion message to display the line. assert type(segment_length) is float, \ "segment_length is not a float. (not numpy.float)" - assert center_line.shape[2] == normal_collection.shape[2] + 1, \ - "Please check the shape (axis-2) of center_line(n_node=n_elems+1) or normal_collection(n_elems)." - assert center_line.shape[0] == normal_collection.shape[0], \ - "The number of timesteps (axis-0) must be equal" - assert center_line.shape[1] == normal_collection.shape[1] == 3, \ + assert center_line.shape[1] == 3, \ "The dimension (axis-1) must be 3" # fmt: on @@ -406,9 +402,9 @@ def compute_link( """ # fmt: off # Format is turned off because I want the assertion message to display the line. - assert type(segment_length) is float, \ - "segment_length is not a float. (not numpy.float)" - assert center_line.shape[2] == normal_collection.shape[2] + 1 == radius.shape[2] + 1, \ + #assert type(segment_length) is float, \ + # "segment_length is not a float. (not numpy.float)" + assert center_line.shape[2] == normal_collection.shape[2] + 1 == radius.shape[1] + 1, \ "Please check the shape (axis-2) of center_line(n_node=n_elems+1) or normal_collection(n_elems)." assert center_line.shape[0] == normal_collection.shape[0] == radius.shape[0], \ "The number of timesteps (axis-0) must be equal" diff --git a/tests/test_rod/test_knot_theory.py b/tests/test_rod/test_knot_theory.py index 9683e2757..1463c2ccd 100644 --- a/tests/test_rod/test_knot_theory.py +++ b/tests/test_rod/test_knot_theory.py @@ -16,6 +16,9 @@ from elastica.rod.rod_base import RodBase from elastica.rod.knot_theory import ( KnotTheoryCompatibleProtocol, + compute_twist, + compute_writhe, + compute_link, _compute_additional_segment, ) @@ -67,6 +70,62 @@ def __init__(self): rod.compute_link() +def test_compute_twist_arithmetic(): + # fmt: off + center_line = np.array( + [[0, 0, 0], [0, 0, 1], [0, 1, 1], [0, 1, 0], [1, 1, 0]], + dtype=np.float64).T[None, ...] + normal_collection = np.array( + [[1, 0, 0], [0, 1, 1], [1, 1, 0], [0,1,0]], + dtype=np.float64).T[None, ...] + # fmt: on + a, b = compute_twist(center_line, normal_collection) + assert np.isclose(a[0], 0.75) + assert_allclose(b[0], np.array([0.25, 0.125, 0.375])) + + +@pytest.mark.parametrize( + "type_str, sol", + [ + ("next_tangent", -0.477268070084), + ("end_to_end", -0.37304522216388), + ("net_tangent", -0.26423311709925), + ], +) +def test_compute_writhe_arithmetic(type_str, sol): + # fmt: off + center_line = np.array( + [[0, 0, 0], [0, 0, 1], [0, 1, 1], [0, 1, 0], [1, 1, 0]], + dtype=np.float64).T[None, ...] + segment_length = 10.0 + # fmt: on + a = compute_writhe(center_line, segment_length, type_str) + assert np.isclose(a[0], sol) + + +@pytest.mark.parametrize( + "type_str, sol", + [ + ("next_tangent", -0.703465518706), + ("end_to_end", -0.4950786438825), + ("net_tangent", -0.321184858244), + ], +) +def test_compute_link_arithmetic(type_str, sol): + # fmt: off + center_line = np.array( + [[0, 0, 0], [0, 0, 1], [0, 1, 1], [0, 1, 0], [1, 1, 0]], + dtype=np.float64).T[None, ...] + normal_collection = np.array( + [[1, 0, 0], [0, 1, 1], [1, 1, 0], [0,1,0]], dtype=np.float64 + ).T[None, ...] + radius = np.array([1, 2, 4, 2], dtype=np.float64)[None,...] + segment_length = 10.0 + # fmt: on + a = compute_link(center_line, normal_collection, radius, segment_length, type_str) + assert np.isclose(a[0], sol) + + @pytest.mark.parametrize("type_str", ["randomstr1", "nextnext_tangent", " "]) def test_knot_theory_compute_additional_segment_integrity(type_str): center_line = np.zeros([1, 3, 10]) From b2faea66086629688f417d0a3f778f0d96da389c Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 23 Mar 2022 01:38:32 -0500 Subject: [PATCH 21/47] update: test for KnotTheory used in rod-mixin --- elastica/rod/knot_theory.py | 5 ++- tests/test_rod/test_knot_theory.py | 49 ++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 3 deletions(-) diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index a2b5606a6..bdf6916c1 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -44,7 +44,7 @@ def director_collection(self) -> np.ndarray: ... @property - def radius_collection(self) -> np.ndarray: + def radius(self) -> np.ndarray: ... @property @@ -132,6 +132,7 @@ def compute_link( Determines the method to compute new segments (elements) added to the rod. Valid inputs are "next_tangent", "end_to_end", "net_tangent", otherwise program uses the center line. """ + print(self.rest_lengths.sum()) return compute_link( self.position_collection[None, ...], self.director_collection[0][None, ...], @@ -280,8 +281,6 @@ def compute_writhe(center_line, segment_length, type_of_additional_segment): """ # fmt: off # Format is turned off because I want the assertion message to display the line. - assert type(segment_length) is float, \ - "segment_length is not a float. (not numpy.float)" assert center_line.shape[1] == 3, \ "The dimension (axis-1) must be 3" # fmt: on diff --git a/tests/test_rod/test_knot_theory.py b/tests/test_rod/test_knot_theory.py index 1463c2ccd..fca73b940 100644 --- a/tests/test_rod/test_knot_theory.py +++ b/tests/test_rod/test_knot_theory.py @@ -70,6 +70,55 @@ def __init__(self): rod.compute_link() +@pytest.mark.parametrize( + "position_collection, director_collection, radius, segment_length, sol_total_twist, sol_total_writhe, sol_total_link", + # fmt: off + [ + ( + np.array([[0, 0, 0], [0, 0, 1], [0, 1, 1], [0, 1, 0], [1, 1, 0]], # position_collection + dtype=np.float64).T, + np.array([[1, 0, 0], [0, 1, 1], [1, 1, 0], [0,1,0]], # director_collection + dtype=np.float64).T[None,...], + np.array([1, 2, 4, 2], dtype=np.float64), # radius + np.array([10.0]), # segment_length + 0.75, # solution total twist + -0.477268070084, # solution total writhe + -0.703465518706 + ), + ], + # solution total link + # fmt: on +) +def test_knot_theory_mixin_methods_arithmetic( + knot_theory, + position_collection, + director_collection, + radius, + segment_length, + sol_total_twist, + sol_total_writhe, + sol_total_link, +): + class TestRod(RodBase, knot_theory.KnotTheory): + def __init__( + self, position_collection, director_collection, radius, segment_length + ): + self.position_collection = position_collection + self.director_collection = director_collection + self.radius = radius + self.rest_lengths = segment_length + + test_rod = TestRod(position_collection, director_collection, radius, segment_length) + + twist = test_rod.compute_twist() + writhe = test_rod.compute_writhe() + link = test_rod.compute_link() + + assert np.isclose(twist, sol_total_twist) + assert np.isclose(writhe, sol_total_writhe) + assert np.isclose(link, sol_total_link) + + def test_compute_twist_arithmetic(): # fmt: off center_line = np.array( From e66c291f812f2fe9b98dc63b8c00a4a7b2d85b6a Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 23 Mar 2022 01:41:32 -0500 Subject: [PATCH 22/47] fmt: flake8 --- elastica/rod/knot_theory.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index bdf6916c1..187aa676e 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -170,11 +170,11 @@ def compute_twist(center_line, normal_collection): # fmt: off # Format is turned off because I want the assertion message to display the line. assert center_line.shape[2] == normal_collection.shape[2] + 1, \ - "Please check the shape (axis-2) of center_line(n_node=n_elems+1) or normal_collection(n_elems)." + "Please check the shape (axis-2) of center_line(n_node=n_elems+1) or normal_collection(n_elems)." assert center_line.shape[0] == normal_collection.shape[0], \ - "The number of timesteps (axis-0) must be equal" + "The number of timesteps (axis-0) must be equal" assert center_line.shape[1] == normal_collection.shape[1] == 3, \ - "The dimension (axis-1) must be 3" + "The dimension (axis-1) must be 3" # fmt: on total_twist, local_twist = _compute_twist(center_line, normal_collection) @@ -282,7 +282,7 @@ def compute_writhe(center_line, segment_length, type_of_additional_segment): # fmt: off # Format is turned off because I want the assertion message to display the line. assert center_line.shape[1] == 3, \ - "The dimension (axis-1) must be 3" + "The dimension (axis-1) must be 3" # fmt: on center_line_with_added_segments, _, _ = _compute_additional_segment( @@ -401,14 +401,12 @@ def compute_link( """ # fmt: off # Format is turned off because I want the assertion message to display the line. - #assert type(segment_length) is float, \ - # "segment_length is not a float. (not numpy.float)" assert center_line.shape[2] == normal_collection.shape[2] + 1 == radius.shape[1] + 1, \ - "Please check the shape (axis-2) of center_line(n_node=n_elems+1) or normal_collection(n_elems)." + "Please check the shape (axis-2) of center_line(n_node=n_elems+1) or normal_collection(n_elems)." assert center_line.shape[0] == normal_collection.shape[0] == radius.shape[0], \ - "The number of timesteps (axis-0) must be equal" + "The number of timesteps (axis-0) must be equal" assert center_line.shape[1] == normal_collection.shape[1] == 3, \ - "The dimension (axis-1) for center_line and normal_collection must be 3" + "The dimension (axis-1) for center_line and normal_collection must be 3" # fmt: on # Compute auxiliary line From 9f26e77fb76893b69c4978c095b45a75c0c8a1ed Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 23 Mar 2022 01:53:57 -0500 Subject: [PATCH 23/47] update coveragerc --- .coveragerc | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/.coveragerc b/.coveragerc index ed02c8e79..52568494e 100644 --- a/.coveragerc +++ b/.coveragerc @@ -4,6 +4,12 @@ exclude_lines = # Don't complain if non-runnable code isn't run: if 0: if __name__ == .__main__.: + ... + pass + def __repr__ + from + import + [run] branch = True From 81a5bcfcd0215b668f17753732a179d3c24d2b75 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 27 Mar 2022 23:41:43 -0500 Subject: [PATCH 24/47] enhancement:Adding how to use LWT func after sim This commit adds an example how to call and use LWT functions after the simulation completed. --- examples/README.md | 4 +- .../SolenoidsCase/solenoid_case.py | 67 +++++++++++++++---- examples/RodContactCase/post_processing.py | 64 +++++++++++++++--- 3 files changed, 109 insertions(+), 26 deletions(-) diff --git a/examples/README.md b/examples/README.md index 0c3b67e4c..18b041f16 100644 --- a/examples/README.md +++ b/examples/README.md @@ -51,8 +51,8 @@ Examples can serve as a starting template for customized usages. * __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 + * __Purpose__: Demonstrates rod self contact with Solenoid example and how to use link-writhe-twist after simulation completed. + * __Features__: CosseratRod, SelonoidsBC, SelfContact, Link-Writhe-Twist ## Functional Examples diff --git a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py index b65b8cd17..2499985f4 100644 --- a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py +++ b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py @@ -5,6 +5,7 @@ from examples.RodContactCase.post_processing import ( plot_video_with_surface, plot_velocity, + plot_link_writhe_twist, ) @@ -44,7 +45,9 @@ class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, Call direction = np.array([0.0, 1.0, 0]) normal = np.array([0.0, 0.0, 1.0]) -start = np.zeros(3,) +start = np.zeros( + 3, +) F_pulling_scalar = 300 @@ -70,9 +73,7 @@ class SolenoidCase(BaseSystemCollection, Constraints, Connections, Forcing, Call class SelonoidsBC(ConstraintBase): - """ - - """ + """ """ def __init__( self, @@ -157,7 +158,9 @@ def constrain_rates(self, rod, time): solenoid_sim.add_forcing_to(sherable_rod).using( EndpointForces, - np.zeros(3,), + np.zeros( + 3, + ), F_pulling_scalar * direction, ramp_up_time=time_start_twist - 1, ) @@ -167,9 +170,7 @@ def constrain_rates(self, rod, time): # 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) @@ -202,7 +203,9 @@ def make_callback(self, system, time, current_step: int): 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, + RodCallBack, + step_skip=step_skip, + callback_params=post_processing_dict, ) # finalize simulation @@ -241,17 +244,53 @@ def make_callback(self, system, time, current_step: int): 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_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 +) + +# Compute link-writhe-twist +normal_history = director_history[:, 0, :, :] +segment_length = 10 * base_length + +type_of_additional_segment = "next_tangent" -twist_density = (time[twist_time_interval_start_idx:twist_time_interval_end_idx]-time_start_twist)*angel_vel_scalar * base_radius +total_twist, local_twist = compute_twist(position_history, normal_history) + +total_link = compute_link( + position_history, + normal_history, + radius_history, + segment_length, + type_of_additional_segment, +) + +total_writhe = compute_writhe( + position_history, segment_length, type_of_additional_segment +) + +# Plot link-writhe-twist +plot_link_writhe_twist( + twist_density, + total_twist[twist_time_interval_start_idx:twist_time_interval_end_idx], + total_writhe[twist_time_interval_start_idx:twist_time_interval_end_idx], + total_link[twist_time_interval_start_idx:twist_time_interval_end_idx], +) +# Save simulation data 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), - twist_density = twist_density, + base_length=np.array(base_length), + twist_density=twist_density, + total_twist=total_twist, + total_writhe=total_writhe, + total_link=total_link, ) diff --git a/examples/RodContactCase/post_processing.py b/examples/RodContactCase/post_processing.py index 564329070..93c5f349e 100644 --- a/examples/RodContactCase/post_processing.py +++ b/examples/RodContactCase/post_processing.py @@ -92,7 +92,7 @@ def plot_video_with_surface( ) if not inst_position.shape[1] == inst_radius.shape[0]: inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] + inst_position[..., 1:] + inst_position[..., :-1] ) rod_scatters[rod_idx]._offsets3d = ( @@ -154,7 +154,7 @@ def plot_video_with_surface( ) if not inst_position.shape[1] == inst_radius.shape[0]: inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] + inst_position[..., 1:] + inst_position[..., :-1] ) rod_lines[rod_idx].set_xdata(inst_position[0]) @@ -218,7 +218,7 @@ def plot_video_with_surface( ) if not inst_position.shape[1] == inst_radius.shape[0]: inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] + inst_position[..., 1:] + inst_position[..., :-1] ) rod_lines[rod_idx].set_xdata(inst_position[2]) @@ -283,7 +283,7 @@ def plot_video_with_surface( ) if not inst_position.shape[1] == inst_radius.shape[0]: inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] + inst_position[..., 1:] + inst_position[..., :-1] ) rod_lines[rod_idx].set_xdata(inst_position[0]) @@ -333,26 +333,38 @@ def plot_velocity( axs[0].set_ylabel("x velocity", fontsize=20) axs[1].plot( - time[:], avg_velocity_rod_one[:, 1], linewidth=3, + time[:], + avg_velocity_rod_one[:, 1], + linewidth=3, ) axs[1].plot( - time[:], avg_velocity_rod_two[:, 1], linewidth=3, + 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, + time[:], + avg_velocity_rod_one[:, 2], + linewidth=3, ) axs[2].plot( - time[:], avg_velocity_rod_two[:, 2], linewidth=3, + 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, + time[:], + total_energy_rod_one[:], + linewidth=3, ) axs[3].semilogy( - time[:], total_energy_rod_two[:], linewidth=3, + time[:], + total_energy_rod_two[:], + linewidth=3, ) axs[3].semilogy( time[:], @@ -373,3 +385,35 @@ def plot_velocity( if SAVE_FIGURE: # fig.savefig(filename) plt.savefig(filename) + + +def plot_link_writhe_twist(twist_density, total_twist, total_writhe, total_link): + + plt.rcParams.update({"font.size": 22}) + fig = plt.figure(figsize=(10, 10), frameon=True, dpi=150) + + axs = [] + axs.append(plt.subplot2grid((1, 1), (0, 0))) + axs[0].plot( + twist_density, + total_twist, + label="twist", + ) + axs[0].plot( + twist_density, + total_writhe, + label="writhe", + ) + axs[0].plot( + twist_density, + total_link, + label="link", + ) + axs[0].set_xlabel("twist density", fontsize=20) + axs[0].set_ylabel("link-twist-writhe", fontsize=20) + axs[0].set_xlim(0, 2.0) + plt.tight_layout() + fig.align_ylabels() + fig.legend(prop={"size": 20}) + fig.savefig("link_twist_writhe.png") + plt.close(plt.gcf()) From ed28f1be7c7c5270d56cda06acabbc297fe34be9 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 27 Mar 2022 23:52:01 -0500 Subject: [PATCH 25/47] update:solenoid case script update. --- .../RodSelfContact/SolenoidsCase/solenoid_case.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py index 2499985f4..0ce2fdcef 100644 --- a/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py +++ b/examples/RodContactCase/RodSelfContact/SolenoidsCase/solenoid_case.py @@ -229,12 +229,7 @@ def make_callback(self, system, time, current_step: int): 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) - +# Compute topological quantities time = np.array(post_processing_dict["time"]) position_history = np.array(post_processing_dict["position"]) radius_history = np.array(post_processing_dict["radius"]) @@ -282,6 +277,10 @@ def make_callback(self, system, time, current_step: int): ) # Save simulation data +import os + +save_folder = os.path.join(os.getcwd(), "data") +os.makedirs(save_folder, exist_ok=True) np.savez( os.path.join(save_folder, "solenoid_case_data.npz"), time=time, From 8e7b3f43f357db74aa5251a5ce20481ae6791b43 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 27 Mar 2022 23:52:28 -0500 Subject: [PATCH 26/47] enhancement:add LWT function usage for plectoneme This commit adds example to show how to use LWT, using the Plectoneme example. --- examples/README.md | 6 +- .../PlectonemesCase/plectoneme_case.py | 83 +++++++++++++++++-- 2 files changed, 78 insertions(+), 11 deletions(-) diff --git a/examples/README.md b/examples/README.md index 18b041f16..e79241e48 100644 --- a/examples/README.md +++ b/examples/README.md @@ -48,10 +48,10 @@ Examples can serve as a starting template for customized usages. * __Features__: CosseratRod, ExternalContact * [RodSelfContact](./RodContactCase/RodSelfContact) * [PlectonemesCase](./RodContactCase/RodSelfContact/PlectonemesCase) - * __Purpose__: Demonstrates rod self contact with Plectoneme example. - * __Features__: CosseratRod, SelonoidsBC, SelfContact + * __Purpose__: Demonstrates rod self contact with Plectoneme example, and how to use link-writhe-twist after simulation completed. + * __Features__: CosseratRod, SelonoidsBC, SelfContact, Link-Writhe-Twist * [SolenoidsCase](./RodContactCase/RodSelfContact/SolenoidsCase) - * __Purpose__: Demonstrates rod self contact with Solenoid example and how to use link-writhe-twist after simulation completed. + * __Purpose__: Demonstrates rod self contact with Solenoid example, and how to use link-writhe-twist after simulation completed. * __Features__: CosseratRod, SelonoidsBC, SelfContact, Link-Writhe-Twist ## Functional Examples diff --git a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py index 54a13175c..1f7e817be 100644 --- a/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py +++ b/examples/RodContactCase/RodSelfContact/PlectonemesCase/plectoneme_case.py @@ -5,6 +5,7 @@ from examples.RodContactCase.post_processing import ( plot_video_with_surface, plot_velocity, + plot_link_writhe_twist, ) @@ -45,7 +46,9 @@ class PlectonemesCase( direction = np.array([0.0, 1.0, 0]) normal = np.array([0.0, 0.0, 1.0]) -start = np.zeros(3,) +start = np.zeros( + 3, +) sherable_rod = CosseratRod.straight_rod( n_elem, @@ -70,9 +73,7 @@ class PlectonemesCase( class SelonoidsBC(ConstraintBase): - """ - - """ + """ """ def __init__( self, @@ -160,9 +161,7 @@ def constrain_rates(self, rod, time): # 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) @@ -187,6 +186,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 @@ -194,7 +194,9 @@ def make_callback(self, system, time, current_step: int): 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, + RodCallBack, + step_skip=step_skip, + callback_params=post_processing_dict, ) # finalize simulation @@ -217,3 +219,68 @@ def make_callback(self, system, time, current_step: int): y_limits=[-0.1, 1.1], z_limits=[-0.5, 0.5], ) + +# Compute topological quantities +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"]) + +# 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 +) + +# Compute link-writhe-twist +normal_history = director_history[:, 0, :, :] +segment_length = 10 * base_length + +type_of_additional_segment = "next_tangent" + +total_twist, local_twist = compute_twist(position_history, normal_history) + +total_link = compute_link( + position_history, + normal_history, + radius_history, + segment_length, + type_of_additional_segment, +) + +total_writhe = compute_writhe( + position_history, segment_length, type_of_additional_segment +) + +# Plot link-writhe-twist +plot_link_writhe_twist( + twist_density, + total_twist[twist_time_interval_start_idx:twist_time_interval_end_idx], + total_writhe[twist_time_interval_start_idx:twist_time_interval_end_idx], + total_link[twist_time_interval_start_idx:twist_time_interval_end_idx], +) + +# Save simulation data +import os + +save_folder = os.path.join(os.getcwd(), "data") +os.makedirs(save_folder, exist_ok=True) +np.savez( + os.path.join(save_folder, "plectoneme_case_data.npz"), + time=time, + position_history=position_history, + radius_history=radius_history, + director_history=director_history, + base_length=np.array(base_length), + twist_density=twist_density, + total_twist=total_twist, + total_writhe=total_writhe, + total_link=total_link, +) From 8cd3770b1bb9e166383e1d216ef7275114a39f6d Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 28 Mar 2022 03:05:33 -0500 Subject: [PATCH 27/47] update: documentation for LWT, suggest batch process --- docs/api/rods.rst | 6 ++-- elastica/rod/knot_theory.py | 63 +++++++++++++++++++++++++++++++------ 2 files changed, 57 insertions(+), 12 deletions(-) diff --git a/docs/api/rods.rst b/docs/api/rods.rst index 628818898..1b620efdc 100644 --- a/docs/api/rods.rst +++ b/docs/api/rods.rst @@ -51,10 +51,10 @@ Cosserat Rod Knot Theory (Mixin) ~~~~~~~~~~~~~~~~~~~ -.. autoclass:: elastica.rod.knot_theory.KnotTheory +.. .. autoclass:: elastica.rod.knot_theory.KnotTheory -.. autoclass:: elastica.rod.knot_theory.KnotTheoryCompatibleProtocol +.. .. autoclass:: elastica.rod.knot_theory.KnotTheoryCompatibleProtocol .. automodule:: elastica.rod.knot_theory - :exclude-members: __weakref__, __init__, KnotTheory, KnotTheoryCompatibleProtocol + :exclude-members: __init__ :members: diff --git a/elastica/rod/knot_theory.py b/elastica/rod/knot_theory.py index 187aa676e..53f1f44ed 100644 --- a/elastica/rod/knot_theory.py +++ b/elastica/rod/knot_theory.py @@ -1,7 +1,13 @@ -__docs__ = """ -This script is for computing the link-writhe-twist of a rod using the method from Klenin & Langowski 2000 paper. +__doc__ = """ +This script is for computing the link-writhe-twist (LWT) of a rod using the method from Klenin & Langowski 2000 paper. +Algorithms are adapted from section S2 of Charles et. al. PRL 2019 paper. -Following codes are adapted from section S2 of Charles et. al. PRL 2019 paper. +Following example cases includes computing LWT quantities to study the bifurcation: + +- `Example case (PlectonemesCase) `_ +- `Example case (SolenoidCase) `_ + +The details discussion is included in `N Charles et. al. PRL (2019) `_. """ __all__ = [ "KnotTheoryCompatibleProtocol", @@ -144,13 +150,23 @@ def compute_link( def compute_twist(center_line, normal_collection): """ - Compute the twist of a rod, using center_line and normal collection. Methods used in this function is - adapted from method 2a Klenin & Langowski 2000 paper. + Compute the twist of a rod, using center_line and normal collection. + + Methods used in this function is adapted from method 2a Klenin & Langowski 2000 paper. - Warnings - -------- - If center line is straight, although the normals of each element is pointing different direction computed twist - will be zero. + .. warning:: If center line is straight, although the normals of each element is pointing different direction computed twist will be zero. + + Typical runtime of this function is longer than simulation steps. While we provide a function to compute + topological quantities at every timesteps, **we highly recommend** to compute LWT during the post-processing + stage.:: + + import elastica + ... + normal_collection = director_collection[:,0,...] # shape of director (time, 3, 3, n_elems) + elastica.compute_twist( + center_line, # shape (time, 3, n_nodes) + normal_collection # shape (time, 3, n_elems) + ) Parameters ---------- @@ -261,8 +277,21 @@ def _compute_twist(center_line, normal_collection): def compute_writhe(center_line, segment_length, type_of_additional_segment): """ This function computes the total writhe history of a rod. + Equations used are from method 1a from Klenin & Langowski 2000 paper. + Typical runtime of this function is longer than simulation steps. While we provide a function to compute + topological quantities at every timesteps, **we highly recommend** to compute LWT during the post-processing + stage.:: + + import elastica + ... + elastica.compute_writhe( + center_line, # shape (time, 3, n_nodes) + segment_length, + type_of_additional_segment="next_tangent" + ) + Parameters ---------- center_line : numpy.ndarray @@ -375,8 +404,24 @@ def compute_link( ): """ This function computes the total link history of a rod. + Equations used are from method 1a from Klenin & Langowski 2000 paper. + Typical runtime of this function is longer than simulation steps. While we provide a function to compute + topological quantities at every timesteps, **we highly recommend** to compute LWT during the post-processing + stage.:: + + import elastica + ... + normal_collection = director_collection[:,0,...] # shape of director (time, 3, 3, n_elems) + elastica.compute_link( + center_line, # shape (time, 3, n_nodes) + normal_collection, # shape (time 3, n_elems) + radius, # shape (time, n_elems) + segment_length, + type_of_additional_segment="next_tangent" + ) + Parameters ---------- center_line : numpy.ndarray From cef4ae40a9661bf0bdb64a09b4652b8c42585f1e Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sun, 10 Apr 2022 00:09:33 -0500 Subject: [PATCH 28/47] update doc config --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index 030965f76..592f12e13 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -86,7 +86,7 @@ "use_repository_button": True, } html_title = "PyElastica" -html_logo = "https://github.com/GazzolaLab/PyElastica/blob/assets/docs/assets/Logo.png?raw=true" +html_logo = "_static/assets/Logo.png" #pygments_style = "sphinx" # Add any paths that contain custom static files (such as style sheets) here, From a06ac6393261b45dcfa61bced1d13484e3edf737 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sun, 10 Apr 2022 00:41:41 -0500 Subject: [PATCH 29/47] update: moduli assignment for Cosserat Rod The logging level of some messages are lowered (warning->info) in order to avoid unnecessary logging. --- elastica/rod/factory_function.py | 53 +++++++++++++++++--------------- 1 file changed, 28 insertions(+), 25 deletions(-) diff --git a/elastica/rod/factory_function.py b/elastica/rod/factory_function.py index 4e92cc89a..c2421df18 100644 --- a/elastica/rod/factory_function.py +++ b/elastica/rod/factory_function.py @@ -1,6 +1,9 @@ __doc__ = """ Factory function to allocate variables for Cosserat Rod""" __all__ = ["allocate"] +import typing +from typing import Optional import warnings +import logging import numpy as np from numpy.testing import assert_allclose from elastica.utils import MaxDimension, Tolerance @@ -16,8 +19,9 @@ def allocate( base_radius, density, nu, - youngs_modulus, - # poisson_ratio, + youngs_modulus: float, + poisson_ratio: Optional[float] = None, + shear_modulus: Optional[float] = None, # alpha_c=0.964, *args, **kwargs @@ -215,7 +219,7 @@ def allocate( ) # Shear/Stretch matrix - shear_modulus = get_shear_modulus(youngs_modulus, kwargs) + shear_modulus = get_shear_modulus(youngs_modulus, poisson_ratio, shear_modulus) # Value taken based on best correlation for Poisson ratio = 0.5, from # "On Timoshenko's correction for shear in vibrating beams" by Kaneko, 1975 @@ -398,31 +402,38 @@ def allocate( ) -def get_shear_modulus(youngs_modulus, kwargs): +def get_shear_modulus( + youngs_modulus: float, + poisson_ratio: Optional[float], + shear_modulus: Optional[float], +): """ - From the kwargs get shear modulus, or compute it. This function contains warnining messages. + From the kwargs get shear modulus, or compute it. This function contains logging messages. Parameters ---------- youngs_modulus : float - - kwargs + poisson_ratio: Optional[float] + shear_modulus: Optional[float] Returns ------- + shear_modulus: float """ + + log = logging.getLogger("rod constructor: get_shear_modulus") + # Shear/Stretch matrix - if kwargs.__contains__("shear_modulus"): + if shear_modulus and not poisson_ratio: # User set shear modulus use that. - shear_modulus = kwargs.get("shear_modulus") + pass - if kwargs.__contains__("shear_modulus") and kwargs.__contains__("poisson_ratio"): + if shear_modulus and poisson_ratio: # User set shear modulus and also a poisson ratio. Do not use poisson ratio and raise warning. - shear_modulus = kwargs.get("shear_modulus") - poisson_ratio = kwargs.get("poisson_ratio") message = ( "Both a Poisson ratio and a shear modulus are provided. " + "In such case, we prioritize the shear modulus." "The Poisson ratio is only used to compute a shear modulus " "so the provided Poisson ratio of ( " + str(poisson_ratio) @@ -432,12 +443,8 @@ def get_shear_modulus(youngs_modulus, kwargs): ) warnings.warn(message, category=UserWarning) - if kwargs.__contains__("poisson_ratio") and ( - not kwargs.__contains__("shear_modulus") - ): - # User set poisson ratio but not the shear modulus. Then use poisson ratio to compute shear modulus, and - # raise a warning. - poisson_ratio = kwargs.get("poisson_ratio") + if poisson_ratio and not shear_modulus: + # Use poisson ratio to compute shear modulus shear_modulus = youngs_modulus / (poisson_ratio + 1.0) message = ( @@ -450,14 +457,11 @@ def get_shear_modulus(youngs_modulus, kwargs): "Use of a Poisson ratio will be depreciated in a future release. " "It is encouraged that you discontinue using a Poisson ratio and instead directly provide the shear_modulus. \n" ) - warnings.warn(message, category=UserWarning) + log.info(message) - if not ( - kwargs.__contains__("poisson_ratio") or kwargs.__contains__("shear_modulus") - ): + if not poisson_ratio and not shear_modulus: # If user does not set poisson ratio or shear modulus, then take poisson ratio as 0.5 and raise warning. poisson_ratio = 0.5 - shear_modulus = youngs_modulus / (poisson_ratio + 1.0) message = ( @@ -469,7 +473,6 @@ def get_shear_modulus(youngs_modulus, kwargs): + ", " "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0)." ) - - warnings.warn(message, category=UserWarning) + log.info(message) return shear_modulus From 6b9dc8c76196f66bcc978cec8b77afd95f9e03d6 Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sun, 10 Apr 2022 13:33:01 -0500 Subject: [PATCH 30/47] remove unused test cases --- tests/test_rod_initialisation.py | 115 ------------------------------- 1 file changed, 115 deletions(-) diff --git a/tests/test_rod_initialisation.py b/tests/test_rod_initialisation.py index 943b13497..43e51ec7b 100644 --- a/tests/test_rod_initialisation.py +++ b/tests/test_rod_initialisation.py @@ -816,121 +816,6 @@ def test_shear_matrix_for_varying_shear_modulus_warning_message_check_if_poisson ) -@pytest.mark.parametrize("n_elems", [5, 10, 50]) -@pytest.mark.parametrize("poisson_ratio", [0.5, 0.3, 1.0]) -def test_shear_matrix_for_varying_poisson_ratio_warning_message_check_if_no_shear_modulus_defined( - n_elems, poisson_ratio -): - """ - This test, is checking if for user defined poisson ratio and validity of shear matrix. - We expect if shear modulus is not defined then a UserWarninig is raised and shear matrix will be computed. - - Returns - ------- - - """ - start = np.array([0.0, 0.0, 0.0]) - direction = np.array([1.0, 0.0, 0.0]) - normal = np.array([0.0, 0.0, 1.0]) - base_length = 1.0 - base_radius = 0.1 - density = 1000 - nu = 0.1 - youngs_modulus = 1e6 - base_area = np.pi * base_radius ** 2 - - with pytest.warns(UserWarning): - mockrod = MockRodForTest.straight_rod( - n_elems, - start, - direction, - normal, - base_length, - base_radius, - density, - nu, - youngs_modulus, - poisson_ratio=poisson_ratio, - ) - - test_shear_matrix = mockrod.shear_matrix - - shear_modulus = youngs_modulus / (1 + poisson_ratio) - correct_shear_matrix = np.zeros((3, 3)) - np.fill_diagonal( - correct_shear_matrix[:], - [ - 0.964 * shear_modulus * base_area, - 0.964 * shear_modulus * base_area, - youngs_modulus * base_area, - ], - ) - - for k in range(n_elems): - assert_allclose( - correct_shear_matrix, - test_shear_matrix[..., k], - atol=Tolerance.atol(), - ) - - -@pytest.mark.parametrize("n_elems", [5, 10, 50]) -def test_shear_matrix_for_no_shear_modulus_or_poisson_ratio_defined_warning_message_check( - n_elems, -): - """ - This test, checks validity of shear matrix if there is no user defined shear modulus and poisson ratio. - Then Elastica uses poisson ratio of 0.5 and computes shear matrix. We expect a UserWarning to be raised. - - Returns - ------- - - """ - start = np.array([0.0, 0.0, 0.0]) - direction = np.array([1.0, 0.0, 0.0]) - normal = np.array([0.0, 0.0, 1.0]) - base_length = 1.0 - base_radius = 0.1 - density = 1000 - nu = 0.1 - youngs_modulus = 1e6 - poisson_ratio = 0.5 - base_area = np.pi * base_radius ** 2 - - with pytest.warns(UserWarning): - mockrod = MockRodForTest.straight_rod( - n_elems, - start, - direction, - normal, - base_length, - base_radius, - density, - nu, - youngs_modulus, - ) - - test_shear_matrix = mockrod.shear_matrix - - shear_modulus = youngs_modulus / (1 + poisson_ratio) - correct_shear_matrix = np.zeros((3, 3)) - np.fill_diagonal( - correct_shear_matrix[:], - [ - 0.964 * shear_modulus * base_area, - 0.964 * shear_modulus * base_area, - youngs_modulus * base_area, - ], - ) - - for k in range(n_elems): - assert_allclose( - correct_shear_matrix, - test_shear_matrix[..., k], - atol=Tolerance.atol(), - ) - - def test_inertia_shear_bend_matrices_for_varying_radius(): """ This test, is checking if for user defined varying radius, validity From 46b42448069f786c859fd8c249c5c9797ed3925d Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 11 Apr 2022 00:13:52 -0500 Subject: [PATCH 31/47] update: remove poisson ratio --- elastica/rod/factory_function.py | 84 +++----------------------------- tests/test_rod_initialisation.py | 24 +-------- 2 files changed, 8 insertions(+), 100 deletions(-) diff --git a/elastica/rod/factory_function.py b/elastica/rod/factory_function.py index c2421df18..6bed11333 100644 --- a/elastica/rod/factory_function.py +++ b/elastica/rod/factory_function.py @@ -20,12 +20,15 @@ def allocate( density, nu, youngs_modulus: float, - poisson_ratio: Optional[float] = None, shear_modulus: Optional[float] = None, # alpha_c=0.964, *args, **kwargs ): + if "poisson_ratio" in kwargs: + raise NameError( + "Poisson's ratio is deprecated for Cosserat Rod for clarity. Please provide shear_modulus instead." + ) # sanity checks here assert n_elements > 1 @@ -219,7 +222,8 @@ def allocate( ) # Shear/Stretch matrix - shear_modulus = get_shear_modulus(youngs_modulus, poisson_ratio, shear_modulus) + if not shear_modulus: + shear_modulus = youngs_modulus / (2.0 * (1.0 + 0.5)) # Value taken based on best correlation for Poisson ratio = 0.5, from # "On Timoshenko's correction for shear in vibrating beams" by Kaneko, 1975 @@ -400,79 +404,3 @@ def allocate( damping_forces, damping_torques, ) - - -def get_shear_modulus( - youngs_modulus: float, - poisson_ratio: Optional[float], - shear_modulus: Optional[float], -): - """ - From the kwargs get shear modulus, or compute it. This function contains logging messages. - - Parameters - ---------- - youngs_modulus : float - poisson_ratio: Optional[float] - shear_modulus: Optional[float] - - Returns - ------- - shear_modulus: float - - """ - - log = logging.getLogger("rod constructor: get_shear_modulus") - - # Shear/Stretch matrix - if shear_modulus and not poisson_ratio: - # User set shear modulus use that. - pass - - if shear_modulus and poisson_ratio: - # User set shear modulus and also a poisson ratio. Do not use poisson ratio and raise warning. - message = ( - "Both a Poisson ratio and a shear modulus are provided. " - "In such case, we prioritize the shear modulus." - "The Poisson ratio is only used to compute a shear modulus " - "so the provided Poisson ratio of ( " - + str(poisson_ratio) - + " ) is being ignored in favor of the provided shear modulus ( " - + str(shear_modulus) - + " ). \n" - ) - warnings.warn(message, category=UserWarning) - - if poisson_ratio and not shear_modulus: - # Use poisson ratio to compute shear modulus - shear_modulus = youngs_modulus / (poisson_ratio + 1.0) - - message = ( - "The given Poisson ratio of " - + str(poisson_ratio) - + " is used only to compute a shear modulus of " - + str(shear_modulus) - + ", " - "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0). " - "Use of a Poisson ratio will be depreciated in a future release. " - "It is encouraged that you discontinue using a Poisson ratio and instead directly provide the shear_modulus. \n" - ) - log.info(message) - - if not poisson_ratio and not shear_modulus: - # If user does not set poisson ratio or shear modulus, then take poisson ratio as 0.5 and raise warning. - poisson_ratio = 0.5 - shear_modulus = youngs_modulus / (poisson_ratio + 1.0) - - message = ( - "Shear modulus cannot be found in kwargs. " - "Poisson ratio " - + str(poisson_ratio) - + " is used to compute shear modulus " - + str(shear_modulus) - + ", " - "using the equation: shear_modulus = youngs_modulus / (poisson_ratio + 1.0)." - ) - log.info(message) - - return shear_modulus diff --git a/tests/test_rod_initialisation.py b/tests/test_rod_initialisation.py index 43e51ec7b..097b2c0fd 100644 --- a/tests/test_rod_initialisation.py +++ b/tests/test_rod_initialisation.py @@ -152,7 +152,6 @@ def straight_rod( density, nu, youngs_modulus, - # poisson_ratio, alpha_c=0.964, *args, **kwargs @@ -758,7 +757,7 @@ def test_shear_matrix_for_varying_shear_modulus(n_elems, shear_modulus): @pytest.mark.parametrize("n_elems", [5, 10, 50]) @pytest.mark.parametrize("shear_modulus", [5e3, 10e3, 50e3]) -def test_shear_matrix_for_varying_shear_modulus_warning_message_check_if_poisson_ratio_defined( +def test_shear_matrix_for_varying_shear_modulus_error_message_check_if_poisson_ratio_defined( n_elems, shear_modulus ): """ @@ -781,7 +780,7 @@ def test_shear_matrix_for_varying_shear_modulus_warning_message_check_if_poisson poisson_ratio = 0.3 base_area = np.pi * base_radius ** 2 - with pytest.warns(UserWarning): + with pytest.raises(NameError): mockrod = MockRodForTest.straight_rod( n_elems, start, @@ -796,25 +795,6 @@ def test_shear_matrix_for_varying_shear_modulus_warning_message_check_if_poisson poisson_ratio=poisson_ratio, ) - test_shear_matrix = mockrod.shear_matrix - - correct_shear_matrix = np.zeros((3, 3)) - np.fill_diagonal( - correct_shear_matrix[:], - [ - 0.964 * shear_modulus * base_area, - 0.964 * shear_modulus * base_area, - youngs_modulus * base_area, - ], - ) - - for k in range(n_elems): - assert_allclose( - correct_shear_matrix, - test_shear_matrix[..., k], - atol=Tolerance.atol(), - ) - def test_inertia_shear_bend_matrices_for_varying_radius(): """ From 8587ac7107afadf802fa8581377cdfa3c87498da Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Mon, 11 Apr 2022 00:16:50 -0500 Subject: [PATCH 32/47] update: add info-level logging message for default poisson ratio --- elastica/rod/factory_function.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/elastica/rod/factory_function.py b/elastica/rod/factory_function.py index 6bed11333..460d22b0b 100644 --- a/elastica/rod/factory_function.py +++ b/elastica/rod/factory_function.py @@ -223,6 +223,11 @@ def allocate( # Shear/Stretch matrix if not shear_modulus: + log = logging.getLogger() + log.info( + """Shear modulus is not explicitly given.\n + In such case, we compute shear_modulus assuming poisson's ratio of 0.5""" + ) shear_modulus = youngs_modulus / (2.0 * (1.0 + 0.5)) # Value taken based on best correlation for Poisson ratio = 0.5, from From 01ada3097f4f4f1472eac54e5cfbc4a1b3a6eb8d Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Fri, 15 Apr 2022 18:04:35 -0500 Subject: [PATCH 33/47] enhancement:initial commit of muscular snake --- .../MuscularSnake/get_connection_vector.py | 54 ++ examples/MuscularSnake/muscle_forces.py | 133 +++++ examples/MuscularSnake/muscular_snake.py | 373 ++++++++++++ examples/MuscularSnake/post_processing.py | 529 ++++++++++++++++++ 4 files changed, 1089 insertions(+) create mode 100644 examples/MuscularSnake/get_connection_vector.py create mode 100644 examples/MuscularSnake/muscle_forces.py create mode 100644 examples/MuscularSnake/muscular_snake.py create mode 100644 examples/MuscularSnake/post_processing.py diff --git a/examples/MuscularSnake/get_connection_vector.py b/examples/MuscularSnake/get_connection_vector.py new file mode 100644 index 000000000..30b852469 --- /dev/null +++ b/examples/MuscularSnake/get_connection_vector.py @@ -0,0 +1,54 @@ +import numpy as np +# Join the two rods +from elastica._linalg import ( + _batch_norm, + _batch_matvec, +) + + + +def get_connection_vector_straight_straight_rod( + rod_one, + rod_two, + rod_one_start_idx, + rod_one_end_idx, +): + + # Compute rod element positions + rod_one_element_position = 0.5 * ( + rod_one.position_collection[..., 1:] + rod_one.position_collection[..., :-1] + ) + rod_one_element_position = rod_one_element_position[:, rod_one_start_idx:rod_one_end_idx] + rod_two_element_position = 0.5 * ( + rod_two.position_collection[..., 1:] + rod_two.position_collection[..., :-1] + ) + + # Lets get the distance between rod elements + distance_vector_rod_one_to_rod_two = ( + rod_two_element_position - rod_one_element_position + ) + distance_vector_rod_one_to_rod_two_norm = _batch_norm( + distance_vector_rod_one_to_rod_two + ) + distance_vector_rod_one_to_rod_two /= distance_vector_rod_one_to_rod_two_norm + + distance_vector_rod_two_to_rod_one = -distance_vector_rod_one_to_rod_two + + rod_one_direction_vec_in_material_frame = _batch_matvec( + rod_one.director_collection[:,:,rod_one_start_idx:rod_one_end_idx], distance_vector_rod_one_to_rod_two + ) + rod_two_direction_vec_in_material_frame = _batch_matvec( + rod_two.director_collection, distance_vector_rod_two_to_rod_one + ) + + offset_btw_rods = distance_vector_rod_one_to_rod_two_norm - ( + rod_one.radius[rod_one_start_idx:rod_one_end_idx] + rod_two.radius + ) + + return ( + rod_one_direction_vec_in_material_frame, + rod_two_direction_vec_in_material_frame, + offset_btw_rods, + ) + + diff --git a/examples/MuscularSnake/muscle_forces.py b/examples/MuscularSnake/muscle_forces.py new file mode 100644 index 000000000..b30c4b6cd --- /dev/null +++ b/examples/MuscularSnake/muscle_forces.py @@ -0,0 +1,133 @@ +__doc__ = """ Muscular snake muscle forces NumPy implementation """ +__all__ = ["MuscleForces"] +import numpy as np +import numba +from numba import njit +from elastica import NoForces +from elastica._calculus import difference_kernel + + +class MuscleForces(NoForces): + """ + This class is for computing muscle forces. Detailed formulation is given in Eq. 2 + Zhang et. al. Nature Comm 2019 paper + + Attributes + ---------- + amplitude : float + Amplitude of muscle forces. + wave_number : float + Wave number for muscle actuation. + side_of_body : int + Depending on which side of body, left or right this variable becomes +1 or -1. + This variable determines the sin wave direction. + time_delay : float + Delay time for muscles. + muscle_start_end_index : numpy.ndarray + 1D (2) array containing data with 'int' type. + Element start and end index of muscle forces. + """ + + def __init__( + self, + amplitude, + wave_number, + arm_length, + time_delay, + side_of_body, + muscle_start_end_index, + step, + post_processing, + ): + """ + + Parameters + ---------- + amplitude : float + Amplitude of muscle forces. + wave_number : float + Wave number for muscle actuation. + arm_length : float + Used to map the torques optimized by CMA into the contraction forces of our muscles. + time_delay : float + Delay time for muscles. + side_of_body : int + Depending on which side of body, left or right this variable becomes +1 or -1. + muscle_start_end_index : numpy.ndarray + 1D (2) array containing data with 'int' type. + Element start and end index of muscle forces. + """ + + self.amplitude = amplitude + self.wave_number = wave_number + self.side_of_body = side_of_body + self.time_delay = time_delay + self.muscle_start_end_index = muscle_start_end_index + + """ + For legacy purposes, the input from the optimizer is given in terms of + torque amplitudes (see Gazzola et al. Royal Society Open Source, 2018). + We then use the same set up and map those values into muscle + contractile forces by dividing them by the arm length. This is captured + through the parameter "factor". This also enables a meaningful + comparison with the continuum snake case of the above reference. + """ + self.amplitude /= arm_length + + self.post_processing = post_processing + self.step = step + self.counter=0 + + def apply_forces(self, system, time: np.float = 0.0): + forces = self._apply_forces( + self.amplitude, + self.wave_number, + self.side_of_body, + time, + self.time_delay, + self.muscle_start_end_index, + system.tangents, + system.external_forces, + ) + + if self.counter % self.step ==0: + self.post_processing["time"].append(time) + self.post_processing["step"].append(self.counter) + self.post_processing["external_forces"].append(forces.copy()) + self.post_processing["element_position"].append(np.cumsum(system.lengths).copy()) + + self.counter += 1 + + + @staticmethod + @njit(cache=True) + def _apply_forces( + amplitude, + wave_number, + side_of_body, + time, + time_delay, + muscle_start_end_index, + tangents, + external_forces, + ): + + real_time = time - time_delay + ramp = real_time if real_time <= 1.0 else 1.0 + factor = 0.0 if real_time <= 0.0 else ramp # max(0.0, ramp) + + forces = np.zeros(external_forces.shape) + + muscle_forces = ( + factor + * amplitude + * (side_of_body * 0.5 * np.sin(wave_number * real_time) + 0.5) + * tangents[..., muscle_start_end_index[0] : muscle_start_end_index[1]] + ) + + forces[ + ..., muscle_start_end_index[0] : muscle_start_end_index[1] + 1 + ] += difference_kernel(muscle_forces) + + external_forces += forces + return forces \ No newline at end of file diff --git a/examples/MuscularSnake/muscular_snake.py b/examples/MuscularSnake/muscular_snake.py new file mode 100644 index 000000000..619535b8d --- /dev/null +++ b/examples/MuscularSnake/muscular_snake.py @@ -0,0 +1,373 @@ +__doc__ = """Muscular snake example from Zhang et. al. Nature Comm 2019 paper.""" +from elastica import * +from examples.MuscularSnake.post_processing import plot_video_with_surface, plot_snake_velocity +from examples.MuscularSnake.muscle_forces import MuscleForces +from elastica.experimental.connection_contact_joint.parallel_connection import SurfaceJointSideBySide +from examples.MuscularSnake.get_connection_vector import get_connection_vector_straight_straight_rod + +# Set base simulator class +class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, Forcing, CallBacks): + pass + +muscular_snake_simulator = MuscularSnakeSimulator() + +# Simulation parameters +final_time = 0.01#16.0 +time_step = 5E-6 +total_steps = int(final_time/time_step) +rendering_fps = 30 +step_skip = int(1.0/(rendering_fps*time_step)) + + +rod_list = [] +# Snake body +n_elem_body = 100 +density_body = 1000 +base_length_body = 1.0 +base_radius_body = 0.025 +E = 1e7 +nu = 4E-3 +shear_modulus = E/ 2*(0.5 + 1.0) +poisson_ratio = 0.5 + +direction = np.array([1.0, 0.0, 0.0]) +normal = np.array([0.0, 0.0, 1.0]) +start = np.array([0.0, 0.0, base_radius_body]) + +snake_body = CosseratRod.straight_rod( + n_elem_body, + start, + direction, + normal, + base_length_body, + base_radius_body, + density_body, + nu, + youngs_modulus=E, + shear_modulus=shear_modulus, +) + +body_elem_length = snake_body.rest_lengths[0] + +# Define muscle fibers +n_muscle_fibers = 8 + +# Muscle force amplitudes +muscle_force_amplitudes = np.array( + [22.96, 22.96, 20.95, 20.95, 9.51, 9.51, 13.7, 13.7] +)[::-1]/2 + +# Set connection index of first node of each muscle with body +muscle_start_connection_index = [4, 4, 33, 33, 23, 23, 61, 61] +muscle_end_connection_index = [] +muscle_glue_connection_index = ( + [] +) # These are the middle node idx of muscles that are glued to body +muscle_rod_list = [] +""" +The muscle density is higher than the physiological one, since +we lump many muscles (SSP-SP, LD and IC) into one actuator. These rods +also represent the two tendons on the sides of the muscle which biologically +have a higher density than the muscle itself. For these reasons,we set the +muscle density to approximately twice the biological value. +""" + +density_muscle = 2000 +E_muscle = 1e4 +nu_muscle = nu +shear_modulus_muscle = E_muscle/ 2*(0.5 + 1.0) + +# Muscle group 1 and 3, define two antagonistic muscle pairs +n_elem_muscle_group_one_to_three = 13*3 +base_length_muscle = 0.39 +""" +In our simulation, we lump many biological tendons into one computational +tendon. As a result, our computational tendon is bigger in size, set as elements other than 4-8 +below. +""" +muscle_radius = np.zeros((n_elem_muscle_group_one_to_three)) +muscle_radius[:] = 0.003 # First set tendon radius for whole rod. +muscle_radius[4*3:9*3] = 0.006 # Change the radius of muscle elements + +for i in range(int(n_muscle_fibers/2)): + + index = muscle_start_connection_index[i] + # Chose which side of body we are attaching the muscles. Note that these muscles are antagonistic pairs. + # So they are at the opposite sides of the body and side_sign determines that. + side_sign = -1 if i % 2 == 0 else 1 + start_muscle = np.array( + [ + index * body_elem_length, + side_sign * (base_radius_body+0.003), + base_radius_body, + ] + ) + + muscle_rod = CosseratRod.straight_rod( + n_elem_muscle_group_one_to_three, + start_muscle, + direction, + normal, + base_length_muscle, + muscle_radius, + density_muscle, + nu_muscle, + youngs_modulus=E_muscle, + shear_modulus=shear_modulus_muscle, + ) + + + """ + The biological tendons have a high Young's modulus E.,but are very slender. + As a result, they resist extension (stretch) but can bend easily. + Due to our decision to lump tendons and in order to mimic the above behavior + of the biological tendons, we use a lower Young's + Modulus and harden the stiffness of the shear and stretch modes only. + Numerically, this is done by putting a pre-factor of 50000 before the + shear/stretch matrix below. The actual value of the prefactor does not matter, + what is important is that it is a high value to high stretch/shear stiffness. + """ + + muscle_rod.shear_matrix[..., :4*3] *= 50000 + muscle_rod.shear_matrix[..., 9*3:] *= 50000 + + muscle_rod_list.append(muscle_rod) + muscle_end_connection_index.append( + index + n_elem_muscle_group_one_to_three + ) + muscle_glue_connection_index.append( + np.hstack((np.arange(0, 4*3, 1, dtype=np.int64), np.arange(9*3, n_elem_muscle_group_one_to_three, 1, dtype=np.int64))) + ) + + +# Muscle group 2 and 4, define two antagonistic muscle pairs +n_elem_muscle_group_two_to_four = 33 +base_length_muscle = 0.33 +""" +In our simulation, we lump many biological tendons into one computational +tendon. As a result, our computational tendon is bigger in size, set as rm_t +below. +""" +muscle_radius = np.zeros((n_elem_muscle_group_two_to_four)) +muscle_radius[:] = 0.003 # First set tendon radius for whole rod. +muscle_radius[4*3:9*3] = 0.006 # Change the radius of muscle elements + +for i in range(int(n_muscle_fibers / 2), n_muscle_fibers): + + index = muscle_start_connection_index[i] + # Chose which side of body we are attaching the muscles. Note that these muscles are antagonistic pairs. + # So they are at the opposite sides of the body and side_sign determines that. + side_sign = -1 if i % 2 == 0 else 1 + start_muscle = np.array( + [ + index * body_elem_length, + side_sign * (base_radius_body+0.003), + base_radius_body, + ] + ) + + muscle_rod = CosseratRod.straight_rod( + n_elem_muscle_group_two_to_four, + start_muscle, + direction, + normal, + base_length_muscle, + muscle_radius, + density_muscle, + nu_muscle, + youngs_modulus=E_muscle, + shear_modulus=shear_modulus_muscle, + ) + + """ + The biological tendons have a high Young's modulus E.,but are very slender. + As a result, they resist extension (stretch) but can bend easily. + Due to our decision to lump tendons and in order to mimic the above behavior + of the biological tendons, we use a lower Young's + Modulus and harden the stiffness of the shear and stretch modes only. + Numerically, this is done by putting a pre-factor of 50000 before the + shear/stretch matrix below. The actual value of the prefactor does not matter, + what is important is that it is a high value to high stretch/shear stiffness. + """ + + muscle_rod.shear_matrix[..., :4*3] *= 50000 + muscle_rod.shear_matrix[..., 9*3:] *= 50000 + + muscle_rod_list.append(muscle_rod) + muscle_end_connection_index.append( + index + n_elem_muscle_group_two_to_four + ) + muscle_glue_connection_index.append( + # np.array([0,1, 2, 3, 9, 10 ], dtype=np.int) + np.hstack((np.arange(0, 4*3, 1, dtype=np.int64), np.arange(9*3, n_elem_muscle_group_two_to_four, 1, dtype=np.int64))) + ) + + +# After initializing the rods append them on to the simulation +rod_list.append(snake_body) +rod_list = rod_list + muscle_rod_list +for _, my_rod in enumerate(rod_list): + muscular_snake_simulator.append(my_rod) + +# Muscle actuation +post_processing_forces_dict_list = [] + +for i in range(n_muscle_fibers): + post_processing_forces_dict_list.append(defaultdict(list)) + muscle_rod = muscle_rod_list[i] + side_of_body = 1 if i % 2 == 0 else -1 + + time_delay = muscle_start_connection_index[::-1][i] * 1.0 / 101.76 + + muscular_snake_simulator.add_forcing_to(muscle_rod).using( + MuscleForces, + amplitude=muscle_force_amplitudes[i], + wave_number=2.0 * np.pi / 1.0, + arm_length=(base_radius_body+0.003), + time_delay=time_delay, + side_of_body=side_of_body, + muscle_start_end_index=np.array([4*3, 9*3], np.int64), + step=step_skip, + post_processing=post_processing_forces_dict_list[i], + ) + + +straight_straight_rod_connection_list = [] +straight_straight_rod_connection_post_processing_dict = defaultdict(list) +for idx, rod_two in enumerate(muscle_rod_list): + rod_one = snake_body + ( + rod_one_direction_vec_in_material_frame, + rod_two_direction_vec_in_material_frame, + offset_btw_rods, + ) = get_connection_vector_straight_straight_rod(rod_one, rod_two,muscle_start_connection_index[idx], muscle_end_connection_index[idx]) + straight_straight_rod_connection_list.append( + [ + rod_one, + rod_two, + rod_one_direction_vec_in_material_frame.copy(), + rod_two_direction_vec_in_material_frame.copy(), + offset_btw_rods.copy(), + ] + ) + for k in range(rod_two.n_elems): + rod_one_index = k + muscle_start_connection_index[idx] + rod_two_index = k + k_conn = rod_one.radius[rod_one_index] * rod_two.radius[rod_two_index] / (rod_one.radius[rod_one_index] + rod_two.radius[rod_two_index]) * body_elem_length * E / (rod_one.radius[rod_one_index] + rod_two.radius[rod_two_index]) + + if k < 12 or k>= 27: + scale = 1*2 + scale_contact = 20 + else: + scale = 0.01*5 + scale_contact = 20 + + muscular_snake_simulator.connect( + first_rod=rod_one, + second_rod=rod_two, + first_connect_idx=rod_one_index, + second_connect_idx=rod_two_index, + ).using( + SurfaceJointSideBySide, + k=k_conn * scale , + nu=1E-4, + k_repulsive=k_conn * scale_contact, + rod_one_direction_vec_in_material_frame=rod_one_direction_vec_in_material_frame[ + ..., k + ], + rod_two_direction_vec_in_material_frame=rod_two_direction_vec_in_material_frame[ + ..., k + ], + offset_btw_rods=offset_btw_rods[k], + post_processing_dict=straight_straight_rod_connection_post_processing_dict, + step_skip=step_skip, + ) + + # Friction forces + # Only apply to the snake body. + gravitational_acc = -9.81 + muscular_snake_simulator.add_forcing_to(snake_body).using( + GravityForces, acc_gravity=np.array([0.0, 0.0, gravitational_acc]) + ) + + origin_plane = np.array([0.0, 0.0, 0.0]) + normal_plane = normal + slip_velocity_tol = 1e-8 + froude = 0.1 + period = 1.0 + mu = base_length_body / (period * period * np.abs(gravitational_acc) * froude) + kinetic_mu_array = np.array( + [1.0*mu, 1.5 * mu, 2.0 * mu] + ) # [forward, backward, sideways] + static_mu_array = 2 * kinetic_mu_array + muscular_snake_simulator.add_forcing_to(snake_body).using( + AnisotropicFrictionalPlane, + k=1E1, + nu=20, + plane_origin=origin_plane, + plane_normal=normal_plane, + slip_velocity_tol=slip_velocity_tol, + static_mu_array=static_mu_array, + kinetic_mu_array=kinetic_mu_array, + ) + +class MuscularSnakeCallBack(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["com"].append( + system.compute_position_center_of_mass() + ) + self.callback_params["radius"].append(system.radius.copy()) + self.callback_params["velocity"].append( + system.velocity_collection.copy() + ) + self.callback_params["avg_velocity"].append( + system.compute_velocity_center_of_mass() + ) + + self.callback_params["center_of_mass"].append( + system.compute_position_center_of_mass() + ) + +post_processing_dict_list = [] + +for idx, rod in enumerate(rod_list): + post_processing_dict_list.append(defaultdict(list)) + muscular_snake_simulator.collect_diagnostics(rod).using( + MuscularSnakeCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_list[idx], + ) + +muscular_snake_simulator.finalize() +timestepper = PositionVerlet() +integrate(timestepper, muscular_snake_simulator, final_time, total_steps) + + +plot_video_with_surface( + post_processing_dict_list, + video_name="muscular_snake.mp4", + fps=rendering_fps, + step=1, + # The following parameters are optional + x_limits=(-0.1, 1.0), # Set bounds on x-axis + y_limits=(-0.3, 0.3), # Set bounds on y-axis + z_limits=(-0.3, 0.3), # Set bounds on z-axis + dpi=100, # Set the quality of the image + vis3D=True, # Turn on 3D visualization + vis2D=True, # Turn on projected (2D) visualization +) + +plot_snake_velocity(post_processing_dict_list[0], period=period, filename="muscular_snake_velocity.png") diff --git a/examples/MuscularSnake/post_processing.py b/examples/MuscularSnake/post_processing.py new file mode 100644 index 000000000..8de47eb14 --- /dev/null +++ b/examples/MuscularSnake/post_processing.py @@ -0,0 +1,529 @@ +import numpy as np +import matplotlib + +matplotlib.use("Agg") # Must be before importing matplotlib.pyplot or pylab! +from matplotlib import pyplot as plt +from matplotlib.colors import to_rgb +from matplotlib import cm +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}) + + folder_name = kwargs.get("folder_name","") + + # 2d case + import matplotlib.animation as animation + + # simulation time + sim_time = np.array(rods_history[0]["time"]) + + # Rod + n_visualized_rods = len(rods_history) # should be one for now + # Rod info + rod_history_unpacker = lambda rod_idx, t_idx: ( + rods_history[rod_idx]["position"][t_idx], + rods_history[rod_idx]["radius"][t_idx], + ) + # Rod center of mass + com_history_unpacker = lambda rod_idx, t_idx: rods_history[rod_idx]["com"][time_idx] + + # Generate target sphere data + sphere_flag = False + if kwargs.__contains__("sphere_history"): + sphere_flag = True + sphere_history = kwargs.get("sphere_history") + n_visualized_spheres = len(sphere_history) # should be one for now + sphere_history_unpacker = lambda sph_idx, t_idx: ( + sphere_history[sph_idx]["position"][t_idx], + sphere_history[sph_idx]["radius"][t_idx], + ) + # color mapping + sphere_cmap = cm.get_cmap("Spectral", n_visualized_spheres) + + # video pre-processing + print("plot scene visualization video") + FFMpegWriter = animation.writers["ffmpeg"] + metadata = dict(title="Movie Test", artist="Matplotlib", comment="Movie support!") + writer = FFMpegWriter(fps=fps, metadata=metadata) + dpi = kwargs.get("dpi", 100) + + xlim = kwargs.get("x_limits", (-1.0, 1.0)) + ylim = kwargs.get("y_limits", (-1.0, 1.0)) + zlim = kwargs.get("z_limits", (-0.05, 1.0)) + + difference = lambda x: x[1] - x[0] + max_axis_length = max(difference(xlim), difference(ylim)) + # The scaling factor from physical space to matplotlib space + scaling_factor = (2 * 0.1) / max_axis_length # Octopus head dimension + scaling_factor *= 2.6e3 # Along one-axis + + if kwargs.get("vis3D", True): + fig = plt.figure(1, figsize=(10, 8), frameon=True, dpi=dpi) + ax = plt.axes(projection="3d") + + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + ax.set_zlim(*zlim) + + time_idx = 0 + rod_lines = [None for _ in range(n_visualized_rods)] + rod_com_lines = [None for _ in range(n_visualized_rods)] + rod_scatters = [None for _ in range(n_visualized_rods)] + + for rod_idx in range(n_visualized_rods): + inst_position, inst_radius = rod_history_unpacker(rod_idx, time_idx) + if not inst_position.shape[1] == inst_radius.shape[0]: + inst_position = 0.5 * (inst_position[..., 1:] + inst_position[..., :-1]) + + rod_scatters[rod_idx] = ax.scatter( + inst_position[0], + inst_position[1], + inst_position[2], + s=np.pi * (scaling_factor * inst_radius) ** 2, + ) + + if sphere_flag: + sphere_artists = [None for _ in range(n_visualized_spheres)] + for sphere_idx in range(n_visualized_spheres): + sphere_position, sphere_radius = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx] = ax.scatter( + sphere_position[0], + sphere_position[1], + sphere_position[2], + s=np.pi * (scaling_factor * sphere_radius) ** 2, + ) + # sphere_radius, + # color=sphere_cmap(sphere_idx),) + ax.add_artist(sphere_artists[sphere_idx]) + + # ax.set_aspect("equal") + video_name_3D = folder_name+"3D_" + video_name + + with writer.saving(fig, video_name_3D, dpi): + with plt.style.context("seaborn-whitegrid"): + for time_idx in tqdm(range(0, sim_time.shape[0], int(step))): + + for rod_idx in range(n_visualized_rods): + inst_position, inst_radius = rod_history_unpacker( + rod_idx, time_idx + ) + if not inst_position.shape[1] == inst_radius.shape[0]: + inst_position = 0.5 * ( + inst_position[..., 1:] + inst_position[..., :-1] + ) + + rod_scatters[rod_idx]._offsets3d = ( + inst_position[0], + inst_position[1], + inst_position[2], + ) + + # rod_scatters[rod_idx].set_offsets(inst_position[:2].T) + rod_scatters[rod_idx].set_sizes( + np.pi * (scaling_factor * inst_radius) ** 2 + ) + + if sphere_flag: + for sphere_idx in range(n_visualized_spheres): + sphere_position, _ = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx]._offsets3d = ( + sphere_position[0], + sphere_position[1], + sphere_position[2], + ) + + writer.grab_frame() + + # Be a good boy and close figures + # https://stackoverflow.com/a/37451036 + # plt.close(fig) alone does not suffice + # See https://github.com/matplotlib/matplotlib/issues/8560/ + plt.close(plt.gcf()) + + if kwargs.get("vis2D", True): + max_axis_length = max(difference(xlim), difference(ylim)) + # The scaling factor from physical space to matplotlib space + scaling_factor = (2 * 0.1) / max_axis_length # Octopus head dimension + scaling_factor *= 2.6e3 # Along one-axis + + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + + time_idx = 0 + rod_lines = [None for _ in range(n_visualized_rods)] + rod_com_lines = [None for _ in range(n_visualized_rods)] + rod_scatters = [None for _ in range(n_visualized_rods)] + + for rod_idx in range(n_visualized_rods): + inst_position, inst_radius = rod_history_unpacker(rod_idx, time_idx) + if not inst_position.shape[1] == inst_radius.shape[0]: + inst_position = 0.5 * (inst_position[..., 1:] + inst_position[..., :-1]) + rod_lines[rod_idx] = ax.plot( + inst_position[0], inst_position[1], "r", lw=0.5 + )[0] + inst_com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx] = ax.plot(inst_com[0], inst_com[1], "k--", lw=2.0)[0] + + rod_scatters[rod_idx] = ax.scatter( + inst_position[0], + inst_position[1], + s=np.pi * (scaling_factor * inst_radius) ** 2, + ) + + if sphere_flag: + sphere_artists = [None for _ in range(n_visualized_spheres)] + for sphere_idx in range(n_visualized_spheres): + sphere_position, sphere_radius = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx] = Circle( + (sphere_position[0], sphere_position[1]), + sphere_radius, + color=sphere_cmap(sphere_idx), + ) + ax.add_artist(sphere_artists[sphere_idx]) + + ax.set_aspect("equal") + video_name_2D = folder_name+"2D_xy_" + video_name + + with writer.saving(fig, video_name_2D, dpi): + with plt.style.context("seaborn-whitegrid"): + for time_idx in tqdm(range(0, sim_time.shape[0], int(step))): + + for rod_idx in range(n_visualized_rods): + inst_position, inst_radius = rod_history_unpacker( + rod_idx, time_idx + ) + if not inst_position.shape[1] == inst_radius.shape[0]: + inst_position = 0.5 * ( + inst_position[..., 1:] + inst_position[..., :-1] + ) + + rod_lines[rod_idx].set_xdata(inst_position[0]) + rod_lines[rod_idx].set_ydata(inst_position[1]) + + com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx].set_xdata(com[0]) + rod_com_lines[rod_idx].set_ydata(com[1]) + + rod_scatters[rod_idx].set_offsets(inst_position[:2].T) + rod_scatters[rod_idx].set_sizes( + np.pi * (scaling_factor * inst_radius) ** 2 + ) + + if sphere_flag: + for sphere_idx in range(n_visualized_spheres): + sphere_position, _ = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx].center = ( + sphere_position[0], + sphere_position[1], + ) + + writer.grab_frame() + + # Be a good boy and close figures + # https://stackoverflow.com/a/37451036 + # plt.close(fig) alone does not suffice + # See https://github.com/matplotlib/matplotlib/issues/8560/ + plt.close(plt.gcf()) + + # Plot zy + max_axis_length = max(difference(zlim), difference(ylim)) + # The scaling factor from physical space to matplotlib space + scaling_factor = (2 * 0.1) / max_axis_length # Octopus head dimension + scaling_factor *= 2.6e3 # Along one-axis + + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + ax.set_xlim(*zlim) + ax.set_ylim(*ylim) + + time_idx = 0 + rod_lines = [None for _ in range(n_visualized_rods)] + rod_com_lines = [None for _ in range(n_visualized_rods)] + rod_scatters = [None for _ in range(n_visualized_rods)] + + for rod_idx in range(n_visualized_rods): + inst_position, inst_radius = rod_history_unpacker(rod_idx, time_idx) + if not inst_position.shape[1] == inst_radius.shape[0]: + inst_position = 0.5 * (inst_position[..., 1:] + inst_position[..., :-1]) + rod_lines[rod_idx] = ax.plot( + inst_position[2], inst_position[1], "r", lw=0.5 + )[0] + inst_com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx] = ax.plot(inst_com[2], inst_com[1], "k--", lw=2.0)[0] + + rod_scatters[rod_idx] = ax.scatter( + inst_position[2], + inst_position[1], + s=np.pi * (scaling_factor * inst_radius) ** 2, + ) + + if sphere_flag: + sphere_artists = [None for _ in range(n_visualized_spheres)] + for sphere_idx in range(n_visualized_spheres): + sphere_position, sphere_radius = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx] = Circle( + (sphere_position[2], sphere_position[1]), + sphere_radius, + color=sphere_cmap(sphere_idx), + ) + ax.add_artist(sphere_artists[sphere_idx]) + + ax.set_aspect("equal") + video_name_2D = folder_name+"2D_zy_" + video_name + + with writer.saving(fig, video_name_2D, dpi): + with plt.style.context("seaborn-whitegrid"): + for time_idx in tqdm(range(0, sim_time.shape[0], int(step))): + + for rod_idx in range(n_visualized_rods): + inst_position, inst_radius = rod_history_unpacker( + rod_idx, time_idx + ) + if not inst_position.shape[1] == inst_radius.shape[0]: + inst_position = 0.5 * ( + inst_position[..., 1:] + inst_position[..., :-1] + ) + + rod_lines[rod_idx].set_xdata(inst_position[2]) + rod_lines[rod_idx].set_ydata(inst_position[1]) + + com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx].set_xdata(com[2]) + rod_com_lines[rod_idx].set_ydata(com[1]) + + rod_scatters[rod_idx].set_offsets( + np.vstack((inst_position[2], inst_position[1])).T + ) + rod_scatters[rod_idx].set_sizes( + np.pi * (scaling_factor * inst_radius) ** 2 + ) + + if sphere_flag: + for sphere_idx in range(n_visualized_spheres): + sphere_position, _ = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx].center = ( + sphere_position[2], + sphere_position[1], + ) + + writer.grab_frame() + + # Be a good boy and close figures + # https://stackoverflow.com/a/37451036 + # plt.close(fig) alone does not suffice + # See https://github.com/matplotlib/matplotlib/issues/8560/ + plt.close(plt.gcf()) + + # Plot xz + fig = plt.figure(2, figsize=(10, 8), frameon=True, dpi=dpi) + ax = fig.add_subplot(111) + ax.set_xlim(*xlim) + ax.set_ylim(*zlim) + + # The scaling factor from physical space to matplotlib space + max_axis_length = max(difference(zlim), difference(xlim)) + scaling_factor = (2 * 0.1) / (max_axis_length) # Octopus head dimension + scaling_factor *= 2.6e3 # Along one-axis + + time_idx = 0 + rod_lines = [None for _ in range(n_visualized_rods)] + rod_com_lines = [None for _ in range(n_visualized_rods)] + rod_scatters = [None for _ in range(n_visualized_rods)] + + for rod_idx in range(n_visualized_rods): + inst_position, inst_radius = rod_history_unpacker(rod_idx, time_idx) + if not inst_position.shape[1] == inst_radius.shape[0]: + inst_position = 0.5 * (inst_position[..., 1:] + inst_position[..., :-1]) + rod_lines[rod_idx] = ax.plot( + inst_position[0], inst_position[2], "r", lw=0.5 + )[0] + inst_com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx] = ax.plot(inst_com[0], inst_com[2], "k--", lw=2.0)[0] + + rod_scatters[rod_idx] = ax.scatter( + inst_position[0], + inst_position[2], + s=np.pi * (scaling_factor * inst_radius) ** 2, + ) + + if sphere_flag: + sphere_artists = [None for _ in range(n_visualized_spheres)] + for sphere_idx in range(n_visualized_spheres): + sphere_position, sphere_radius = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx] = Circle( + (sphere_position[0], sphere_position[2]), + sphere_radius, + color=sphere_cmap(sphere_idx), + ) + ax.add_artist(sphere_artists[sphere_idx]) + + ax.set_aspect("equal") + video_name_2D = folder_name+"2D_xz_" + video_name + + with writer.saving(fig, video_name_2D, dpi): + with plt.style.context("seaborn-whitegrid"): + for time_idx in tqdm(range(0, sim_time.shape[0], int(step))): + + for rod_idx in range(n_visualized_rods): + inst_position, inst_radius = rod_history_unpacker( + rod_idx, time_idx + ) + if not inst_position.shape[1] == inst_radius.shape[0]: + inst_position = 0.5 * ( + inst_position[..., 1:] + inst_position[..., :-1] + ) + + rod_lines[rod_idx].set_xdata(inst_position[0]) + rod_lines[rod_idx].set_ydata(inst_position[2]) + + com = com_history_unpacker(rod_idx, time_idx) + rod_com_lines[rod_idx].set_xdata(com[0]) + rod_com_lines[rod_idx].set_ydata(com[2]) + + rod_scatters[rod_idx].set_offsets( + np.vstack((inst_position[0], inst_position[2])).T + ) + rod_scatters[rod_idx].set_sizes( + np.pi * (scaling_factor * inst_radius) ** 2 + ) + + if sphere_flag: + for sphere_idx in range(n_visualized_spheres): + sphere_position, _ = sphere_history_unpacker( + sphere_idx, time_idx + ) + sphere_artists[sphere_idx].center = ( + sphere_position[0], + sphere_position[2], + ) + + writer.grab_frame() + + # Be a good boy and close figures + # https://stackoverflow.com/a/37451036 + # plt.close(fig) alone does not suffice + # See https://github.com/matplotlib/matplotlib/issues/8560/ + plt.close(plt.gcf()) + +def plot_snake_velocity( + plot_params: dict, + period, + filename="slithering_snake_velocity.png", +): + time_per_period = np.array(plot_params["time"]) / period + avg_velocity = np.array(plot_params["avg_velocity"]) + + [ + velocity_in_direction_of_rod, + velocity_in_rod_roll_dir, + _, + _, + ] = compute_projected_velocity(plot_params, period) + + fig = plt.figure(figsize=(10, 8), frameon=True, dpi=150) + ax = fig.add_subplot(111) + ax.grid(b=True, which="minor", color="k", linestyle="--") + ax.grid(b=True, which="major", color="k", linestyle="-") + ax.plot( + time_per_period[:], velocity_in_direction_of_rod[:, 0], "r-", label="forward" + ) + ax.plot( + time_per_period[:], + velocity_in_rod_roll_dir[:, 1], + c=to_rgb("xkcd:bluish"), + label="lateral", + ) + ax.plot(time_per_period[:], avg_velocity[:, 2], "k-", label="normal") + fig.legend(prop={"size": 20}) + fig.savefig(filename) + +def compute_projected_velocity(plot_params: dict, period): + + time_per_period = np.array(plot_params["time"]) / period + avg_velocity = np.array(plot_params["avg_velocity"]) + center_of_mass = np.array(plot_params["center_of_mass"]) + + # Compute rod velocity in rod direction. We need to compute that because, + # after snake starts to move it chooses an arbitrary direction, which does not + # have to be initial tangent direction of the rod. Thus we need to project the + # snake velocity with respect to its new tangent and roll direction, after that + # we will get the correct forward and lateral speed. After this projection + # lateral velocity of the snake has to be oscillating between + and - values with + # zero mean. + + # Number of steps in one period. + period_step = int(period / (time_per_period[-1] - time_per_period[-2])) + 1 + number_of_period = int(time_per_period.shape[0] / period_step) + # Center of mass position averaged in one period + center_of_mass_averaged_over_one_period = np.zeros((number_of_period - 2, 3)) + for i in range(1, number_of_period - 1): + # position of center of mass averaged over one period + center_of_mass_averaged_over_one_period[i - 1] = np.mean( + center_of_mass[(i + 1) * period_step : (i + 2) * period_step] + - center_of_mass[(i + 0) * period_step : (i + 1) * period_step], + axis=0, + ) + + # Average the rod directions over multiple periods and get the direction of the rod. + direction_of_rod = np.mean(center_of_mass_averaged_over_one_period, axis=0) + direction_of_rod /= np.linalg.norm(direction_of_rod, ord=2) + print("direction of rod " + str(direction_of_rod)) + + # Compute the projected rod velocity in the direction of the rod + velocity_mag_in_direction_of_rod = np.einsum( + "ji,i->j", avg_velocity, direction_of_rod + ) + velocity_in_direction_of_rod = np.einsum( + "j,i->ji", velocity_mag_in_direction_of_rod, direction_of_rod + ) + + # Get the lateral or roll velocity of the rod after subtracting its projected + # velocity in the direction of rod + velocity_in_rod_roll_dir = avg_velocity - velocity_in_direction_of_rod + + # Compute the average velocity over the simulation, this can be used for optimizing snake + # for fastest forward velocity. We start after first period, because of the ramping up happens + # in first period. + average_velocity_over_simulation = np.mean( + velocity_in_direction_of_rod[period_step * 2 :], axis=0 + ) + + return ( + velocity_in_direction_of_rod, + velocity_in_rod_roll_dir, + average_velocity_over_simulation[0], + average_velocity_over_simulation[1], + ) + From a079f22a2fd20bbce6acfc9adba3b773e56766e3 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 16 Apr 2022 09:50:48 -0500 Subject: [PATCH 34/47] fix:path issues and black fmt muscular snake --- examples/MuscularSnake/muscular_snake.py | 141 ++++++++++++++--------- 1 file changed, 84 insertions(+), 57 deletions(-) diff --git a/examples/MuscularSnake/muscular_snake.py b/examples/MuscularSnake/muscular_snake.py index 619535b8d..c49dd0832 100644 --- a/examples/MuscularSnake/muscular_snake.py +++ b/examples/MuscularSnake/muscular_snake.py @@ -1,22 +1,35 @@ __doc__ = """Muscular snake example from Zhang et. al. Nature Comm 2019 paper.""" +import sys + +sys.path.append("../../") from elastica import * -from examples.MuscularSnake.post_processing import plot_video_with_surface, plot_snake_velocity +from examples.MuscularSnake.post_processing import ( + plot_video_with_surface, + plot_snake_velocity, +) from examples.MuscularSnake.muscle_forces import MuscleForces -from elastica.experimental.connection_contact_joint.parallel_connection import SurfaceJointSideBySide -from examples.MuscularSnake.get_connection_vector import get_connection_vector_straight_straight_rod +from elastica.experimental.connection_contact_joint.parallel_connection import ( + SurfaceJointSideBySide, +) +from examples.MuscularSnake.get_connection_vector import ( + get_connection_vector_straight_straight_rod, +) # Set base simulator class -class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, Forcing, CallBacks): +class MuscularSnakeSimulator( + BaseSystemCollection, Constraints, Connections, Forcing, CallBacks +): pass + muscular_snake_simulator = MuscularSnakeSimulator() # Simulation parameters -final_time = 0.01#16.0 -time_step = 5E-6 -total_steps = int(final_time/time_step) +final_time = 16.0 +time_step = 5e-6 +total_steps = int(final_time / time_step) rendering_fps = 30 -step_skip = int(1.0/(rendering_fps*time_step)) +step_skip = int(1.0 / (rendering_fps * time_step)) rod_list = [] @@ -26,8 +39,8 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For base_length_body = 1.0 base_radius_body = 0.025 E = 1e7 -nu = 4E-3 -shear_modulus = E/ 2*(0.5 + 1.0) +nu = 4e-3 +shear_modulus = E / 2 * (0.5 + 1.0) poisson_ratio = 0.5 direction = np.array([1.0, 0.0, 0.0]) @@ -53,9 +66,9 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For n_muscle_fibers = 8 # Muscle force amplitudes -muscle_force_amplitudes = np.array( - [22.96, 22.96, 20.95, 20.95, 9.51, 9.51, 13.7, 13.7] -)[::-1]/2 +muscle_force_amplitudes = ( + np.array([22.96, 22.96, 20.95, 20.95, 9.51, 9.51, 13.7, 13.7])[::-1] / 2 +) # Set connection index of first node of each muscle with body muscle_start_connection_index = [4, 4, 33, 33, 23, 23, 61, 61] @@ -75,10 +88,10 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For density_muscle = 2000 E_muscle = 1e4 nu_muscle = nu -shear_modulus_muscle = E_muscle/ 2*(0.5 + 1.0) +shear_modulus_muscle = E_muscle / 2 * (0.5 + 1.0) # Muscle group 1 and 3, define two antagonistic muscle pairs -n_elem_muscle_group_one_to_three = 13*3 +n_elem_muscle_group_one_to_three = 13 * 3 base_length_muscle = 0.39 """ In our simulation, we lump many biological tendons into one computational @@ -87,9 +100,9 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For """ muscle_radius = np.zeros((n_elem_muscle_group_one_to_three)) muscle_radius[:] = 0.003 # First set tendon radius for whole rod. -muscle_radius[4*3:9*3] = 0.006 # Change the radius of muscle elements +muscle_radius[4 * 3 : 9 * 3] = 0.006 # Change the radius of muscle elements -for i in range(int(n_muscle_fibers/2)): +for i in range(int(n_muscle_fibers / 2)): index = muscle_start_connection_index[i] # Chose which side of body we are attaching the muscles. Note that these muscles are antagonistic pairs. @@ -98,9 +111,9 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For start_muscle = np.array( [ index * body_elem_length, - side_sign * (base_radius_body+0.003), + side_sign * (base_radius_body + 0.003), base_radius_body, - ] + ] ) muscle_rod = CosseratRod.straight_rod( @@ -116,7 +129,6 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For shear_modulus=shear_modulus_muscle, ) - """ The biological tendons have a high Young's modulus E.,but are very slender. As a result, they resist extension (stretch) but can bend easily. @@ -128,15 +140,18 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For what is important is that it is a high value to high stretch/shear stiffness. """ - muscle_rod.shear_matrix[..., :4*3] *= 50000 - muscle_rod.shear_matrix[..., 9*3:] *= 50000 + muscle_rod.shear_matrix[..., : 4 * 3] *= 50000 + muscle_rod.shear_matrix[..., 9 * 3 :] *= 50000 muscle_rod_list.append(muscle_rod) - muscle_end_connection_index.append( - index + n_elem_muscle_group_one_to_three - ) + muscle_end_connection_index.append(index + n_elem_muscle_group_one_to_three) muscle_glue_connection_index.append( - np.hstack((np.arange(0, 4*3, 1, dtype=np.int64), np.arange(9*3, n_elem_muscle_group_one_to_three, 1, dtype=np.int64))) + np.hstack( + ( + np.arange(0, 4 * 3, 1, dtype=np.int64), + np.arange(9 * 3, n_elem_muscle_group_one_to_three, 1, dtype=np.int64), + ) + ) ) @@ -150,7 +165,7 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For """ muscle_radius = np.zeros((n_elem_muscle_group_two_to_four)) muscle_radius[:] = 0.003 # First set tendon radius for whole rod. -muscle_radius[4*3:9*3] = 0.006 # Change the radius of muscle elements +muscle_radius[4 * 3 : 9 * 3] = 0.006 # Change the radius of muscle elements for i in range(int(n_muscle_fibers / 2), n_muscle_fibers): @@ -161,9 +176,9 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For start_muscle = np.array( [ index * body_elem_length, - side_sign * (base_radius_body+0.003), + side_sign * (base_radius_body + 0.003), base_radius_body, - ] + ] ) muscle_rod = CosseratRod.straight_rod( @@ -190,16 +205,19 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For what is important is that it is a high value to high stretch/shear stiffness. """ - muscle_rod.shear_matrix[..., :4*3] *= 50000 - muscle_rod.shear_matrix[..., 9*3:] *= 50000 + muscle_rod.shear_matrix[..., : 4 * 3] *= 50000 + muscle_rod.shear_matrix[..., 9 * 3 :] *= 50000 muscle_rod_list.append(muscle_rod) - muscle_end_connection_index.append( - index + n_elem_muscle_group_two_to_four - ) + muscle_end_connection_index.append(index + n_elem_muscle_group_two_to_four) muscle_glue_connection_index.append( # np.array([0,1, 2, 3, 9, 10 ], dtype=np.int) - np.hstack((np.arange(0, 4*3, 1, dtype=np.int64), np.arange(9*3, n_elem_muscle_group_two_to_four, 1, dtype=np.int64))) + np.hstack( + ( + np.arange(0, 4 * 3, 1, dtype=np.int64), + np.arange(9 * 3, n_elem_muscle_group_two_to_four, 1, dtype=np.int64), + ) + ) ) @@ -223,10 +241,10 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For MuscleForces, amplitude=muscle_force_amplitudes[i], wave_number=2.0 * np.pi / 1.0, - arm_length=(base_radius_body+0.003), + arm_length=(base_radius_body + 0.003), time_delay=time_delay, side_of_body=side_of_body, - muscle_start_end_index=np.array([4*3, 9*3], np.int64), + muscle_start_end_index=np.array([4 * 3, 9 * 3], np.int64), step=step_skip, post_processing=post_processing_forces_dict_list[i], ) @@ -240,7 +258,12 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For rod_one_direction_vec_in_material_frame, rod_two_direction_vec_in_material_frame, offset_btw_rods, - ) = get_connection_vector_straight_straight_rod(rod_one, rod_two,muscle_start_connection_index[idx], muscle_end_connection_index[idx]) + ) = get_connection_vector_straight_straight_rod( + rod_one, + rod_two, + muscle_start_connection_index[idx], + muscle_end_connection_index[idx], + ) straight_straight_rod_connection_list.append( [ rod_one, @@ -253,13 +276,20 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For for k in range(rod_two.n_elems): rod_one_index = k + muscle_start_connection_index[idx] rod_two_index = k - k_conn = rod_one.radius[rod_one_index] * rod_two.radius[rod_two_index] / (rod_one.radius[rod_one_index] + rod_two.radius[rod_two_index]) * body_elem_length * E / (rod_one.radius[rod_one_index] + rod_two.radius[rod_two_index]) + k_conn = ( + rod_one.radius[rod_one_index] + * rod_two.radius[rod_two_index] + / (rod_one.radius[rod_one_index] + rod_two.radius[rod_two_index]) + * body_elem_length + * E + / (rod_one.radius[rod_one_index] + rod_two.radius[rod_two_index]) + ) - if k < 12 or k>= 27: - scale = 1*2 + if k < 12 or k >= 27: + scale = 1 * 2 scale_contact = 20 else: - scale = 0.01*5 + scale = 0.01 * 5 scale_contact = 20 muscular_snake_simulator.connect( @@ -269,8 +299,8 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For second_connect_idx=rod_two_index, ).using( SurfaceJointSideBySide, - k=k_conn * scale , - nu=1E-4, + k=k_conn * scale, + nu=1e-4, k_repulsive=k_conn * scale_contact, rod_one_direction_vec_in_material_frame=rod_one_direction_vec_in_material_frame[ ..., k @@ -297,12 +327,12 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For period = 1.0 mu = base_length_body / (period * period * np.abs(gravitational_acc) * froude) kinetic_mu_array = np.array( - [1.0*mu, 1.5 * mu, 2.0 * mu] + [1.0 * mu, 1.5 * mu, 2.0 * mu] ) # [forward, backward, sideways] static_mu_array = 2 * kinetic_mu_array muscular_snake_simulator.add_forcing_to(snake_body).using( AnisotropicFrictionalPlane, - k=1E1, + k=1e1, nu=20, plane_origin=origin_plane, plane_normal=normal_plane, @@ -311,6 +341,7 @@ class MuscularSnakeSimulator(BaseSystemCollection, Constraints, Connections, For kinetic_mu_array=kinetic_mu_array, ) + class MuscularSnakeCallBack(CallBackBaseClass): def __init__(self, step_skip: int, callback_params: dict): CallBackBaseClass.__init__(self) @@ -319,20 +350,13 @@ def __init__(self, step_skip: int, callback_params: dict): 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["com"].append( - system.compute_position_center_of_mass() - ) + self.callback_params["position"].append(system.position_collection.copy()) + self.callback_params["com"].append(system.compute_position_center_of_mass()) self.callback_params["radius"].append(system.radius.copy()) - self.callback_params["velocity"].append( - system.velocity_collection.copy() - ) + self.callback_params["velocity"].append(system.velocity_collection.copy()) self.callback_params["avg_velocity"].append( system.compute_velocity_center_of_mass() ) @@ -341,6 +365,7 @@ def make_callback(self, system, time, current_step: int): system.compute_position_center_of_mass() ) + post_processing_dict_list = [] for idx, rod in enumerate(rod_list): @@ -370,4 +395,6 @@ def make_callback(self, system, time, current_step: int): vis2D=True, # Turn on projected (2D) visualization ) -plot_snake_velocity(post_processing_dict_list[0], period=period, filename="muscular_snake_velocity.png") +plot_snake_velocity( + post_processing_dict_list[0], period=period, filename="muscular_snake_velocity.png" +) From 0d8f5c29543e9bfc3a0e0a80862e3e14fd920322 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sat, 16 Apr 2022 22:28:24 -0500 Subject: [PATCH 35/47] bugfix:friction forces defined on wrong rods --- examples/MuscularSnake/muscular_snake.py | 52 ++++++++++++------------ 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/examples/MuscularSnake/muscular_snake.py b/examples/MuscularSnake/muscular_snake.py index c49dd0832..da5dc4379 100644 --- a/examples/MuscularSnake/muscular_snake.py +++ b/examples/MuscularSnake/muscular_snake.py @@ -313,33 +313,33 @@ class MuscularSnakeSimulator( step_skip=step_skip, ) - # Friction forces - # Only apply to the snake body. - gravitational_acc = -9.81 - muscular_snake_simulator.add_forcing_to(snake_body).using( - GravityForces, acc_gravity=np.array([0.0, 0.0, gravitational_acc]) - ) +# Friction forces +# Only apply to the snake body. +gravitational_acc = -9.81 +muscular_snake_simulator.add_forcing_to(snake_body).using( + GravityForces, acc_gravity=np.array([0.0, 0.0, gravitational_acc]) +) - origin_plane = np.array([0.0, 0.0, 0.0]) - normal_plane = normal - slip_velocity_tol = 1e-8 - froude = 0.1 - period = 1.0 - mu = base_length_body / (period * period * np.abs(gravitational_acc) * froude) - kinetic_mu_array = np.array( - [1.0 * mu, 1.5 * mu, 2.0 * mu] - ) # [forward, backward, sideways] - static_mu_array = 2 * kinetic_mu_array - muscular_snake_simulator.add_forcing_to(snake_body).using( - AnisotropicFrictionalPlane, - k=1e1, - nu=20, - plane_origin=origin_plane, - plane_normal=normal_plane, - slip_velocity_tol=slip_velocity_tol, - static_mu_array=static_mu_array, - kinetic_mu_array=kinetic_mu_array, - ) +origin_plane = np.array([0.0, 0.0, 0.0]) +normal_plane = normal +slip_velocity_tol = 1e-8 +froude = 0.1 +period = 1.0 +mu = base_length_body / (period * period * np.abs(gravitational_acc) * froude) +kinetic_mu_array = np.array( + [1.0 * mu, 1.5 * mu, 2.0 * mu] +) # [forward, backward, sideways] +static_mu_array = 2 * kinetic_mu_array +muscular_snake_simulator.add_forcing_to(snake_body).using( + AnisotropicFrictionalPlane, + k=1e1, + nu=20, + plane_origin=origin_plane, + plane_normal=normal_plane, + slip_velocity_tol=slip_velocity_tol, + static_mu_array=static_mu_array, + kinetic_mu_array=kinetic_mu_array, +) class MuscularSnakeCallBack(CallBackBaseClass): From 9bc2310259da7c3c68d1e7a01e0ad7426820255b Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 17 Apr 2022 10:29:40 -0500 Subject: [PATCH 36/47] Update: Example readme.md updated and muscular snake description added --- examples/README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/README.md b/examples/README.md index c4af06dc5..52800518d 100644 --- a/examples/README.md +++ b/examples/README.md @@ -21,6 +21,9 @@ Examples can serve as a starting template for customized usages. * [ContinuumSnakeCase](./ContinuumSnakeCase) * __Purpose__: Demonstrate simple case of modeling biological creature using PyElastica. The example use friction to create slithering snake, and optimize the speed using [CMA](https://github.com/CMA-ES/pycma). * __Features__: CosseratRod, MuscleTorques, AnisotropicFrictionalPlane, Gravity, CMA Optimization + * [MuscularSnake](./MuscularSnake) + * __Purpose__: Example of [Parallel connection module](../elastica/experimental/connection_contact_joint/parallel_connection.py) and customized [Force module](./MuscularSnake/muscle_forces.py) to implement muscular snake. + * __Features__: MuscleForces(custom implemented) * [ButterflyCase](./ButterflyCase) * __Purpose__: Demonstrate simple restoration with initial strain. * __Features__: CosseratRod From e5884245aacbaf458fde7f8aac244f160146ee9c Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 17 Apr 2022 18:07:11 -0500 Subject: [PATCH 37/47] update:get_connection_vector function is updated Now it takes the indexes that user want to connect. --- .../parallel_connection.py | 19 +++++-- .../parallel_connection_example.py | 2 +- .../MuscularSnake/get_connection_vector.py | 54 ------------------- examples/MuscularSnake/muscular_snake.py | 9 ++-- 4 files changed, 21 insertions(+), 63 deletions(-) delete mode 100644 examples/MuscularSnake/get_connection_vector.py diff --git a/elastica/experimental/connection_contact_joint/parallel_connection.py b/elastica/experimental/connection_contact_joint/parallel_connection.py index d9a74ab67..5113e07b9 100644 --- a/elastica/experimental/connection_contact_joint/parallel_connection.py +++ b/elastica/experimental/connection_contact_joint/parallel_connection.py @@ -18,15 +18,25 @@ def get_connection_vector_straight_straight_rod( rod_one, rod_two, + rod_one_idx, + rod_two_idx, ): + rod_one_start_idx, rod_one_end_idx = rod_one_idx + rod_two_start_idx, rod_two_end_idx = rod_two_idx # Compute rod element positions rod_one_element_position = 0.5 * ( rod_one.position_collection[..., 1:] + rod_one.position_collection[..., :-1] ) + rod_one_element_position = rod_one_element_position[ + :, rod_one_start_idx:rod_one_end_idx + ] rod_two_element_position = 0.5 * ( rod_two.position_collection[..., 1:] + rod_two.position_collection[..., :-1] ) + rod_two_element_position = rod_two_element_position[ + :, rod_two_start_idx:rod_two_end_idx + ] # Lets get the distance between rod elements distance_vector_rod_one_to_rod_two = ( @@ -40,14 +50,17 @@ def get_connection_vector_straight_straight_rod( distance_vector_rod_two_to_rod_one = -distance_vector_rod_one_to_rod_two rod_one_direction_vec_in_material_frame = _batch_matvec( - rod_one.director_collection, distance_vector_rod_one_to_rod_two + rod_one.director_collection[:, :, rod_one_start_idx:rod_one_end_idx], + distance_vector_rod_one_to_rod_two, ) rod_two_direction_vec_in_material_frame = _batch_matvec( - rod_two.director_collection, distance_vector_rod_two_to_rod_one + rod_two.director_collection[:, :, rod_two_start_idx:rod_two_end_idx], + distance_vector_rod_two_to_rod_one, ) offset_btw_rods = distance_vector_rod_one_to_rod_two_norm - ( - rod_one.radius + rod_two.radius + rod_one.radius[rod_one_start_idx:rod_one_end_idx] + + rod_two.radius[rod_two_start_idx:rod_two_end_idx] ) return ( diff --git a/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py b/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py index a0de0c46d..08db51cd3 100644 --- a/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py +++ b/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py @@ -99,7 +99,7 @@ def apply_forces(self, system, time: np.float64 = 0.0): rod_one_direction_vec_in_material_frame, rod_two_direction_vec_in_material_frame, offset_btw_rods, -) = get_connection_vector_straight_straight_rod(rod_one, rod_two) +) = get_connection_vector_straight_straight_rod(rod_one, rod_two, (0, n_elem), (0, n_elem)) for i in range(n_elem): parallel_connection_sim.connect( diff --git a/examples/MuscularSnake/get_connection_vector.py b/examples/MuscularSnake/get_connection_vector.py deleted file mode 100644 index 30b852469..000000000 --- a/examples/MuscularSnake/get_connection_vector.py +++ /dev/null @@ -1,54 +0,0 @@ -import numpy as np -# Join the two rods -from elastica._linalg import ( - _batch_norm, - _batch_matvec, -) - - - -def get_connection_vector_straight_straight_rod( - rod_one, - rod_two, - rod_one_start_idx, - rod_one_end_idx, -): - - # Compute rod element positions - rod_one_element_position = 0.5 * ( - rod_one.position_collection[..., 1:] + rod_one.position_collection[..., :-1] - ) - rod_one_element_position = rod_one_element_position[:, rod_one_start_idx:rod_one_end_idx] - rod_two_element_position = 0.5 * ( - rod_two.position_collection[..., 1:] + rod_two.position_collection[..., :-1] - ) - - # Lets get the distance between rod elements - distance_vector_rod_one_to_rod_two = ( - rod_two_element_position - rod_one_element_position - ) - distance_vector_rod_one_to_rod_two_norm = _batch_norm( - distance_vector_rod_one_to_rod_two - ) - distance_vector_rod_one_to_rod_two /= distance_vector_rod_one_to_rod_two_norm - - distance_vector_rod_two_to_rod_one = -distance_vector_rod_one_to_rod_two - - rod_one_direction_vec_in_material_frame = _batch_matvec( - rod_one.director_collection[:,:,rod_one_start_idx:rod_one_end_idx], distance_vector_rod_one_to_rod_two - ) - rod_two_direction_vec_in_material_frame = _batch_matvec( - rod_two.director_collection, distance_vector_rod_two_to_rod_one - ) - - offset_btw_rods = distance_vector_rod_one_to_rod_two_norm - ( - rod_one.radius[rod_one_start_idx:rod_one_end_idx] + rod_two.radius - ) - - return ( - rod_one_direction_vec_in_material_frame, - rod_two_direction_vec_in_material_frame, - offset_btw_rods, - ) - - diff --git a/examples/MuscularSnake/muscular_snake.py b/examples/MuscularSnake/muscular_snake.py index da5dc4379..9a0bfcbdb 100644 --- a/examples/MuscularSnake/muscular_snake.py +++ b/examples/MuscularSnake/muscular_snake.py @@ -10,11 +10,10 @@ from examples.MuscularSnake.muscle_forces import MuscleForces from elastica.experimental.connection_contact_joint.parallel_connection import ( SurfaceJointSideBySide, -) -from examples.MuscularSnake.get_connection_vector import ( get_connection_vector_straight_straight_rod, ) + # Set base simulator class class MuscularSnakeSimulator( BaseSystemCollection, Constraints, Connections, Forcing, CallBacks @@ -261,8 +260,8 @@ class MuscularSnakeSimulator( ) = get_connection_vector_straight_straight_rod( rod_one, rod_two, - muscle_start_connection_index[idx], - muscle_end_connection_index[idx], + (muscle_start_connection_index[idx],muscle_end_connection_index[idx]), + (0, rod_two.n_elems) ) straight_straight_rod_connection_list.append( [ @@ -333,7 +332,7 @@ class MuscularSnakeSimulator( muscular_snake_simulator.add_forcing_to(snake_body).using( AnisotropicFrictionalPlane, k=1e1, - nu=20, + nu=40, plane_origin=origin_plane, plane_normal=normal_plane, slip_velocity_tol=slip_velocity_tol, From edba6e7cd576870f7a627d8c077e0f5c186b826b Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Wed, 27 Apr 2022 18:41:20 -0500 Subject: [PATCH 38/47] refactor: try to remove duplicated commentaries --- elastica/rod/factory_function.py | 292 ++++++++++++++----------------- 1 file changed, 128 insertions(+), 164 deletions(-) diff --git a/elastica/rod/factory_function.py b/elastica/rod/factory_function.py index 460d22b0b..267f7a074 100644 --- a/elastica/rod/factory_function.py +++ b/elastica/rod/factory_function.py @@ -1,7 +1,7 @@ __doc__ = """ Factory function to allocate variables for Cosserat Rod""" __all__ = ["allocate"] import typing -from typing import Optional +from typing import Optional, Tuple import warnings import logging import numpy as np @@ -20,12 +20,19 @@ def allocate( density, nu, youngs_modulus: float, + nu_for_torques: Optional[float] = None, shear_modulus: Optional[float] = None, - # alpha_c=0.964, + position: Optional[np.ndarray] = None, + directors: Optional[np.ndarray] = None, + rest_sigma: Optional[np.ndarray] = None, + rest_kappa: Optional[np.ndarray] = None, *args, - **kwargs + **kwargs, ): + log = logging.getLogger() + if "poisson_ratio" in kwargs: + # Deprecation warning for poission_ratio raise NameError( "Poisson's ratio is deprecated for Cosserat Rod for clarity. Please provide shear_modulus instead." ) @@ -36,35 +43,14 @@ def allocate( assert np.sqrt(np.dot(normal, normal)) > Tolerance.atol() assert np.sqrt(np.dot(direction, direction)) > Tolerance.atol() - # Set the position array - position = np.zeros((MaxDimension.value(), n_elements + 1)) - # check if position is in kwargs, if it is use user defined position otherwise generate position - if kwargs.__contains__("position"): - position_temp = np.array(kwargs["position"]) - - # Check the shape of the input position - assert position_temp.shape == (MaxDimension.value(), n_elements + 1), ( - "Given position shape is not correct, it should be " - + str(position.shape) - + " but instead " - + str(position_temp.shape) - ) - # Check if the start position of the rod and first entry of position array are the same - assert_allclose( - position_temp[..., 0], - start, - atol=Tolerance.atol(), - err_msg=str( - "First entry of position" + " (" + str(position_temp[..., 0]) + " ) " - " is different than start " + " (" + str(start) + " ) " - ), - ) - position = position_temp.copy() - - else: + # check if position is given. + if position is None: # Generate straight and uniform rod + # Set the position array + position = np.zeros((MaxDimension.value(), n_elements + 1)) end = start + direction * base_length for i in range(0, 3): position[i, ...] = np.linspace(start[i], end[i], n_elements + 1) + _position_validity_checker(position, start, n_elements) # Compute rest lengths and tangents position_diff = position[..., 1:] - position[..., :-1] @@ -72,66 +58,9 @@ def allocate( tangents = position_diff / rest_lengths normal /= np.linalg.norm(normal) - # Set the directors matrix - directors = np.zeros((MaxDimension.value(), MaxDimension.value(), n_elements)) - # check if directors is in kwargs, if it use user defined directors otherwise generate directors - if kwargs.__contains__("directors"): - directors_temp = np.array(kwargs["directors"]) - - # Check the shape of input directors - assert directors_temp.shape == ( - MaxDimension.value(), - MaxDimension.value(), - n_elements, - ), ( - " Given directors shape is not correct, it should be " - + str(directors.shape) - + " but instead " - + str(directors_temp.shape) - ) - - # Check if d1, d2, d3 are unit vectors - d1 = directors_temp[0, ...] - d2 = directors_temp[1, ...] - d3 = directors_temp[2, ...] - assert_allclose( - _batch_norm(d1), - np.ones((n_elements)), - atol=Tolerance.atol(), - err_msg=(" d1 vector of input director matrix is not unit vector "), - ) - assert_allclose( - _batch_norm(d2), - np.ones((n_elements)), - atol=Tolerance.atol(), - err_msg=(" d2 vector of input director matrix is not unit vector "), - ) - assert_allclose( - _batch_norm(d3), - np.ones((n_elements)), - atol=Tolerance.atol(), - err_msg=(" d3 vector of input director matrix is not unit vector "), - ) - - # Check if d3xd1 = d2 - assert_allclose( - _batch_cross(d3, d1), - d2, - atol=Tolerance.atol(), - err_msg=(" d3 x d1 != d2 of input director matrix"), - ) - - # Check if computed tangents from position is the same with d3 - assert_allclose( - tangents, - d3, - atol=Tolerance.atol(), - err_msg=" Tangent vector computed using node positions is different than d3 vector of input directors", - ) - - directors[:] = directors_temp[:] - - else: + if directors is None: # Generate straight uniform rod + # Set the directors matrix + directors = np.zeros((MaxDimension.value(), MaxDimension.value(), n_elements)) # Construct directors using tangents and normal normal_collection = np.repeat(normal[:, np.newaxis], n_elements, axis=1) # Check if rod normal and rod tangent are perpendicular to each other otherwise @@ -145,42 +74,27 @@ def allocate( directors[0, ...] = normal_collection directors[1, ...] = _batch_cross(tangents, normal_collection) directors[2, ...] = tangents + _directors_validity_checker(directors, tangents, n_elements) # Set radius array radius = np.zeros((n_elements)) # Check if the user input radius is valid radius_temp = np.array(base_radius) - assert radius_temp.ndim < 2, ( - "Input radius shape is not correct " - + str(radius_temp.shape) - + " It should be " - + str(radius.shape) - + " or single floating number " - ) + _assert_dim(radius_temp, 2, "radius") radius[:] = radius_temp # Check if the elements of radius are greater than tolerance - for k in range(n_elements): - assert radius[k] > Tolerance.atol(), ( - " Radius has to be greater than 0" + " Check you radius input!" - ) + assert np.all(radius > Tolerance.atol()), " Radius has to be greater than 0." # Set density array density_array = np.zeros((n_elements)) # Check if the user input density is valid density_temp = np.array(density) - assert density_temp.ndim < 2, ( - "Input density shape is not correct " - + str(density_temp.shape) - + " It should be " - + str(density_array.shape) - + " or single floating number " - ) + _assert_dim(density_temp, 2, "density") density_array[:] = density_temp # Check if the elements of density are greater than tolerance - for k in range(n_elements): - assert density_array[k] > Tolerance.atol(), ( - " Density has to be greater than 0" + " Check you density input!" - ) + assert np.all( + density_array > Tolerance.atol() + ), " Density has to be greater than 0." # Second moment of inertia A0 = np.pi * radius * radius @@ -205,7 +119,7 @@ def allocate( # sanity check of mass second moment of inertia if (mass_second_moment_of_inertia < Tolerance.atol()).all(): message = "Mass moment of inertia matrix smaller than tolerance, please check provided radius, density and length." - warnings.warn(message, category=UserWarning) + log.warning(message) # Inverse of second moment of inertia inv_mass_second_moment_of_inertia = np.zeros( @@ -223,7 +137,6 @@ def allocate( # Shear/Stretch matrix if not shear_modulus: - log = logging.getLogger() log.info( """Shear modulus is not explicitly given.\n In such case, we compute shear_modulus assuming poisson's ratio of 0.5""" @@ -232,8 +145,7 @@ def allocate( # Value taken based on best correlation for Poisson ratio = 0.5, from # "On Timoshenko's correction for shear in vibrating beams" by Kaneko, 1975 - alpha_c = kwargs.get("alpha_c", 0.964) - + alpha_c = 0.964 shear_matrix = np.zeros( (MaxDimension.value(), MaxDimension.value(), n_elements), np.float64 ) @@ -260,9 +172,11 @@ def allocate( shear_modulus * I0_3[i], ], ) - for k in range(n_elements): - for i in range(0, MaxDimension.value()): - assert bend_matrix[i, i, k] > Tolerance.atol() + for i in range(0, MaxDimension.value()): + assert np.all( + bend_matrix[i, i, :] > Tolerance.atol() + ), " Bend matrix has to be greater than 0." + # Compute bend matrix in Voronoi Domain bend_matrix = ( bend_matrix[..., 1:] * rest_lengths[1:] @@ -281,60 +195,32 @@ def allocate( dissipation_constant_for_forces = np.zeros((n_elements)) # Check if the user input nu is valid nu_temp = np.array(nu) - assert nu_temp.ndim < 2, ( - "Input dissipation constant(nu) for forces shape is not correct " - + str(nu_temp.shape) - + " It should be " - + str(dissipation_constant_for_forces.shape) - + " or single floating number " - ) + _assert_dim(nu_temp, 2, "dissipation constant (nu) for forces)") dissipation_constant_for_forces[:] = nu # Check if the elements of dissipation constant greater than tolerance - for k in range(n_elements): - assert dissipation_constant_for_forces[k] >= 0.0, ( - " Dissipation constant has to be equal or greater than 0 " - + " Check your dissipation constant(nu) input!" - ) - - dissipation_constant_for_torques = np.zeros((n_elements)) - if kwargs.__contains__("nu_for_torques"): - temp_nu_for_torques = np.array(kwargs["nu_for_torques"]) - assert temp_nu_for_torques.ndim < 2, ( - "Input dissipation constant(nu) for torques shape is not correct " - + str(temp_nu_for_torques.shape) - + " It should be " - + str(dissipation_constant_for_torques.shape) - + " or single floating number " - ) - dissipation_constant_for_torques[:] = temp_nu_for_torques + assert np.all( + dissipation_constant_for_forces >= 0.0 + ), " Dissipation constant(nu) has to be equal or greater than 0." + # Custom nu for torques + if nu_for_torques is None: + dissipation_constant_for_torques = dissipation_constant_for_forces.copy() else: - dissipation_constant_for_torques[:] = dissipation_constant_for_forces + dissipation_constant_for_torques = np.asarray(nu_for_torques) + _assert_dim( + dissipation_constant_for_torques, 2, "dissipation constant (nu) for torque)" + ) # Generate rest sigma and rest kappa, use user input if defined # set rest strains and curvature to be zero at start # if found in kwargs modify (say for curved rod) - rest_sigma = np.zeros((MaxDimension.value(), n_elements)) - if kwargs.__contains__("rest_sigma"): - temp_rest_sigma = np.array(kwargs["rest_sigma"]) - assert temp_rest_sigma.shape == rest_sigma.shape, ( - "Input rest sigma shape is not correct " - + str(temp_rest_sigma.shape) - + " It should be " - + str(rest_sigma.shape) - ) - rest_sigma[:] = temp_rest_sigma - - rest_kappa = np.zeros((MaxDimension.value(), n_elements - 1)) - if kwargs.__contains__("rest_kappa"): - temp_rest_kappa = np.array(kwargs["rest_kappa"]) - assert temp_rest_kappa.shape == rest_kappa.shape, ( - "Input rest kappa shape is not correct " - + str(temp_rest_kappa.shape) - + " It should be " - + str(rest_kappa.shape) - ) - rest_kappa[:] = temp_rest_kappa + if rest_sigma is None: + rest_sigma = np.zeros((MaxDimension.value(), n_elements)) + _assert_shape(rest_sigma, (MaxDimension.value(), n_elements), "rest_sigma") + + if rest_kappa is None: + rest_kappa = np.zeros((MaxDimension.value(), n_elements - 1)) + _assert_shape(rest_kappa, (MaxDimension.value(), n_elements - 1), "rest_kappa") # Compute rest voronoi length rest_voronoi_lengths = 0.5 * (rest_lengths[1:] + rest_lengths[:-1]) @@ -409,3 +295,81 @@ def allocate( damping_forces, damping_torques, ) + + +def _assert_dim(vector, max_dim: int, name: str): + assert vector.ndim < max_dim, ( + f"Input {name} dimension is not correct {vector.shape}" + + f" It should be maximum {max_dim}D vector or single floating number." + ) + + +def _assert_shape(array: np.ndarray, expected_shape: Tuple[int], name: str): + assert array.shape == expected_shape, ( + f"Given {name} shape is not correct, it should be " + + str(expected_shape) + + " but instead " + + str(array.shape) + ) + + +def _position_validity_checker(position, start, n_elements): + """Checker on user-defined position validity""" + _assert_shape(position, (MaxDimension.value(), n_elements + 1), "position") + + # Check if the start position of the rod and first entry of position array are the same + assert_allclose( + position[..., 0], + start, + atol=Tolerance.atol(), + err_msg=str( + "First entry of position" + " (" + str(position[..., 0]) + " ) " + " is different than start " + " (" + str(start) + " ) " + ), + ) + + +def _directors_validity_checker(directors, tangents, n_elements): + """Checker on user-defined directors validity""" + _assert_shape( + directors, (MaxDimension.value(), MaxDimension.value(), n_elements), "directors" + ) + + # Check if d1, d2, d3 are unit vectors + d1 = directors[0, ...] + d2 = directors[1, ...] + d3 = directors[2, ...] + assert_allclose( + _batch_norm(d1), + np.ones((n_elements)), + atol=Tolerance.atol(), + err_msg=(" d1 vector of input director matrix is not unit vector "), + ) + assert_allclose( + _batch_norm(d2), + np.ones((n_elements)), + atol=Tolerance.atol(), + err_msg=(" d2 vector of input director matrix is not unit vector "), + ) + assert_allclose( + _batch_norm(d3), + np.ones((n_elements)), + atol=Tolerance.atol(), + err_msg=(" d3 vector of input director matrix is not unit vector "), + ) + + # Check if d3xd1 = d2 + assert_allclose( + _batch_cross(d3, d1), + d2, + atol=Tolerance.atol(), + err_msg=(" d3 x d1 != d2 of input director matrix"), + ) + + # Check if computed tangents from position is the same with d3 + assert_allclose( + tangents, + d3, + atol=Tolerance.atol(), + err_msg=" Tangent vector computed using node positions is different than d3 vector of input directors", + ) From 59d85f399c123c73eb97f6fa41c5f99e534f4b0e Mon Sep 17 00:00:00 2001 From: Bhosale Date: Sun, 15 May 2022 21:12:24 -0500 Subject: [PATCH 39/47] fmt: add blacken option for examples --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index a16207652..4c032cbb6 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ black: @black --version - @black --required-version 21.12b0 elastica tests + @black --required-version 21.12b0 elastica tests examples black_check: From be9bdc43e67b3d4259b080bf275745ca31b5fffb Mon Sep 17 00:00:00 2001 From: Bhosale Date: Sun, 15 May 2022 21:12:39 -0500 Subject: [PATCH 40/47] fmt: blacken examples --- .../parallel_connection_example.py | 39 ++++++++++++++---- examples/MuscularSnake/muscle_forces.py | 13 +++--- examples/MuscularSnake/muscular_snake.py | 4 +- examples/MuscularSnake/post_processing.py | 41 ++++++++++--------- examples/RestartExample/restart_example.py | 13 ++++-- .../rod_rod_contact_inclined_validation.py | 18 ++++---- .../rod_rod_contact_parallel_validation.py | 18 ++++---- 7 files changed, 93 insertions(+), 53 deletions(-) diff --git a/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py b/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py index 08db51cd3..5ca8d3aaf 100644 --- a/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py +++ b/examples/ExperimentalCases/ParallelConnectionExample/parallel_connection_example.py @@ -6,7 +6,10 @@ # FIXME without appending sys.path make it more generic sys.path.append("../../../../") from elastica import * -from elastica.experimental.connection_contact_joint.parallel_connection import get_connection_vector_straight_straight_rod, SurfaceJointSideBySide +from elastica.experimental.connection_contact_joint.parallel_connection import ( + get_connection_vector_straight_straight_rod, + SurfaceJointSideBySide, +) from elastica._calculus import difference_kernel from examples.JointCases.joint_cases_postprocessing import ( plot_position, @@ -82,7 +85,11 @@ class ParallelConnection( # Apply a contraction force on rod one. class ContractionForce(NoForces): - def __init__(self, ramp, force_mag, ): + def __init__( + self, + ramp, + force_mag, + ): self.ramp = ramp self.force_mag = force_mag @@ -90,26 +97,42 @@ def apply_forces(self, system, time: np.float64 = 0.0): # Ramp the force factor = min(1.0, time / self.ramp) - system.external_forces[:] -= factor * difference_kernel(self.force_mag * system.tangents) + system.external_forces[:] -= factor * difference_kernel( + self.force_mag * system.tangents + ) + -parallel_connection_sim.add_forcing_to(rod_one).using(ContractionForce, ramp=0.5, force_mag=1. ) +parallel_connection_sim.add_forcing_to(rod_one).using( + ContractionForce, ramp=0.5, force_mag=1.0 +) # Connect rod 1 and rod 2 ( rod_one_direction_vec_in_material_frame, rod_two_direction_vec_in_material_frame, offset_btw_rods, -) = get_connection_vector_straight_straight_rod(rod_one, rod_two, (0, n_elem), (0, n_elem)) +) = get_connection_vector_straight_straight_rod( + rod_one, rod_two, (0, n_elem), (0, n_elem) +) for i in range(n_elem): parallel_connection_sim.connect( first_rod=rod_one, second_rod=rod_two, first_connect_idx=i, second_connect_idx=i ).using( - SurfaceJointSideBySide, k=1E2, nu=1E-5, k_repulsive = 1E3, rod_one_direction_vec_in_material_frame=rod_one_direction_vec_in_material_frame[:,i], - rod_two_direction_vec_in_material_frame = rod_two_direction_vec_in_material_frame[:,i], - offset_btw_rods = offset_btw_rods[i] + SurfaceJointSideBySide, + k=1e2, + nu=1e-5, + k_repulsive=1e3, + rod_one_direction_vec_in_material_frame=rod_one_direction_vec_in_material_frame[ + :, i + ], + rod_two_direction_vec_in_material_frame=rod_two_direction_vec_in_material_frame[ + :, i + ], + offset_btw_rods=offset_btw_rods[i], ) # k=kg/s2 nu=kg/s 1e-2 + class ParallelConnecitonCallback(CallBackBaseClass): """ Call back function for parallel connection diff --git a/examples/MuscularSnake/muscle_forces.py b/examples/MuscularSnake/muscle_forces.py index b30c4b6cd..3a59c6dde 100644 --- a/examples/MuscularSnake/muscle_forces.py +++ b/examples/MuscularSnake/muscle_forces.py @@ -37,7 +37,7 @@ def __init__( side_of_body, muscle_start_end_index, step, - post_processing, + post_processing, ): """ @@ -76,7 +76,7 @@ def __init__( self.post_processing = post_processing self.step = step - self.counter=0 + self.counter = 0 def apply_forces(self, system, time: np.float = 0.0): forces = self._apply_forces( @@ -90,15 +90,16 @@ def apply_forces(self, system, time: np.float = 0.0): system.external_forces, ) - if self.counter % self.step ==0: + if self.counter % self.step == 0: self.post_processing["time"].append(time) self.post_processing["step"].append(self.counter) self.post_processing["external_forces"].append(forces.copy()) - self.post_processing["element_position"].append(np.cumsum(system.lengths).copy()) + self.post_processing["element_position"].append( + np.cumsum(system.lengths).copy() + ) self.counter += 1 - @staticmethod @njit(cache=True) def _apply_forces( @@ -130,4 +131,4 @@ def _apply_forces( ] += difference_kernel(muscle_forces) external_forces += forces - return forces \ No newline at end of file + return forces diff --git a/examples/MuscularSnake/muscular_snake.py b/examples/MuscularSnake/muscular_snake.py index 9a0bfcbdb..d7495b78c 100644 --- a/examples/MuscularSnake/muscular_snake.py +++ b/examples/MuscularSnake/muscular_snake.py @@ -260,8 +260,8 @@ class MuscularSnakeSimulator( ) = get_connection_vector_straight_straight_rod( rod_one, rod_two, - (muscle_start_connection_index[idx],muscle_end_connection_index[idx]), - (0, rod_two.n_elems) + (muscle_start_connection_index[idx], muscle_end_connection_index[idx]), + (0, rod_two.n_elems), ) straight_straight_rod_connection_list.append( [ diff --git a/examples/MuscularSnake/post_processing.py b/examples/MuscularSnake/post_processing.py index 8de47eb14..44b36eb2b 100644 --- a/examples/MuscularSnake/post_processing.py +++ b/examples/MuscularSnake/post_processing.py @@ -12,16 +12,16 @@ def plot_video_with_surface( - rods_history: Sequence[Dict], - video_name="video.mp4", - fps=60, - step=1, - vis2D=True, - **kwargs, + rods_history: Sequence[Dict], + video_name="video.mp4", + fps=60, + step=1, + vis2D=True, + **kwargs, ): plt.rcParams.update({"font.size": 22}) - folder_name = kwargs.get("folder_name","") + folder_name = kwargs.get("folder_name", "") # 2d case import matplotlib.animation as animation @@ -115,7 +115,7 @@ def plot_video_with_surface( ax.add_artist(sphere_artists[sphere_idx]) # ax.set_aspect("equal") - video_name_3D = folder_name+"3D_" + video_name + video_name_3D = folder_name + "3D_" + video_name with writer.saving(fig, video_name_3D, dpi): with plt.style.context("seaborn-whitegrid"): @@ -127,7 +127,7 @@ def plot_video_with_surface( ) if not inst_position.shape[1] == inst_radius.shape[0]: inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] + inst_position[..., 1:] + inst_position[..., :-1] ) rod_scatters[rod_idx]._offsets3d = ( @@ -206,7 +206,7 @@ def plot_video_with_surface( ax.add_artist(sphere_artists[sphere_idx]) ax.set_aspect("equal") - video_name_2D = folder_name+"2D_xy_" + video_name + video_name_2D = folder_name + "2D_xy_" + video_name with writer.saving(fig, video_name_2D, dpi): with plt.style.context("seaborn-whitegrid"): @@ -218,7 +218,7 @@ def plot_video_with_surface( ) if not inst_position.shape[1] == inst_radius.shape[0]: inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] + inst_position[..., 1:] + inst_position[..., :-1] ) rod_lines[rod_idx].set_xdata(inst_position[0]) @@ -297,7 +297,7 @@ def plot_video_with_surface( ax.add_artist(sphere_artists[sphere_idx]) ax.set_aspect("equal") - video_name_2D = folder_name+"2D_zy_" + video_name + video_name_2D = folder_name + "2D_zy_" + video_name with writer.saving(fig, video_name_2D, dpi): with plt.style.context("seaborn-whitegrid"): @@ -309,7 +309,7 @@ def plot_video_with_surface( ) if not inst_position.shape[1] == inst_radius.shape[0]: inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] + inst_position[..., 1:] + inst_position[..., :-1] ) rod_lines[rod_idx].set_xdata(inst_position[2]) @@ -390,7 +390,7 @@ def plot_video_with_surface( ax.add_artist(sphere_artists[sphere_idx]) ax.set_aspect("equal") - video_name_2D = folder_name+"2D_xz_" + video_name + video_name_2D = folder_name + "2D_xz_" + video_name with writer.saving(fig, video_name_2D, dpi): with plt.style.context("seaborn-whitegrid"): @@ -402,7 +402,7 @@ def plot_video_with_surface( ) if not inst_position.shape[1] == inst_radius.shape[0]: inst_position = 0.5 * ( - inst_position[..., 1:] + inst_position[..., :-1] + inst_position[..., 1:] + inst_position[..., :-1] ) rod_lines[rod_idx].set_xdata(inst_position[0]) @@ -437,10 +437,11 @@ def plot_video_with_surface( # See https://github.com/matplotlib/matplotlib/issues/8560/ plt.close(plt.gcf()) + def plot_snake_velocity( - plot_params: dict, - period, - filename="slithering_snake_velocity.png", + plot_params: dict, + period, + filename="slithering_snake_velocity.png", ): time_per_period = np.array(plot_params["time"]) / period avg_velocity = np.array(plot_params["avg_velocity"]) @@ -469,6 +470,7 @@ def plot_snake_velocity( fig.legend(prop={"size": 20}) fig.savefig(filename) + def compute_projected_velocity(plot_params: dict, period): time_per_period = np.array(plot_params["time"]) / period @@ -494,7 +496,7 @@ def compute_projected_velocity(plot_params: dict, period): center_of_mass[(i + 1) * period_step : (i + 2) * period_step] - center_of_mass[(i + 0) * period_step : (i + 1) * period_step], axis=0, - ) + ) # Average the rod directions over multiple periods and get the direction of the rod. direction_of_rod = np.mean(center_of_mass_averaged_over_one_period, axis=0) @@ -526,4 +528,3 @@ def compute_projected_velocity(plot_params: dict, period): average_velocity_over_simulation[0], average_velocity_over_simulation[1], ) - diff --git a/examples/RestartExample/restart_example.py b/examples/RestartExample/restart_example.py index 816b0285b..e97cbb3ec 100644 --- a/examples/RestartExample/restart_example.py +++ b/examples/RestartExample/restart_example.py @@ -3,9 +3,11 @@ """ import sys + sys.path.append("../../") from elastica import * + class RestartExampleSimulator(BaseSystemCollection, Constraints, Forcing, CallBacks): pass @@ -62,7 +64,7 @@ class RestartExampleSimulator(BaseSystemCollection, Constraints, Forcing, CallBa # After finalize you can load restart file. This step is important for current implementation of restart functions, # it is required to load restart files after the finalize step. -LOAD_FROM_RESTART = False +LOAD_FROM_RESTART = False SAVE_DATA_RESTART = True restart_file_location = "data/" @@ -76,11 +78,16 @@ class RestartExampleSimulator(BaseSystemCollection, Constraints, Forcing, CallBa dt = 0.01 * dl total_steps = int(final_time / dt) -time = integrate(timestepper, restart_example_simulator, final_time, total_steps, restart_time=restart_time) +time = integrate( + timestepper, + restart_example_simulator, + final_time, + total_steps, + restart_time=restart_time, +) # Save all the systems appended on the simulator class. Since in this example have only one system, under the # `restart_file_location` directory there is one file called system_0.npz . For each system appended on the simulator # separate system_#.npz file will be created. if SAVE_DATA_RESTART: save_state(restart_example_simulator, restart_file_location, time, True) - diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py index 139bbe8f3..8f71addb7 100644 --- a/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py +++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_inclined_validation.py @@ -35,7 +35,9 @@ class InclinedRodRodContact( shear_modulus = E / (poisson_ratio + 1.0) # Rod orientations -start = np.zeros(3,) +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)]) @@ -56,7 +58,7 @@ class InclinedRodRodContact( nu, E, poisson_ratio, - shear_modulus = shear_modulus, + shear_modulus=shear_modulus, ) rod_one.velocity_collection[:] += 0.05 * -normal.reshape(3, 1) @@ -94,9 +96,7 @@ class InclinedRodRodContact( # Add call backs class RodCallBack(CallBackBaseClass): - """ - - """ + """ """ def __init__(self, step_skip: int, callback_params: dict): CallBackBaseClass.__init__(self) @@ -130,7 +130,9 @@ def make_callback(self, system, time, current_step: int): ) # 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, + RodCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_rod1, ) post_processing_dict_rod2 = defaultdict( @@ -138,7 +140,9 @@ def make_callback(self, system, time, current_step: int): ) # 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, + RodCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_rod2, ) inclined_rod_rod_contact_sim.finalize() diff --git a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py index 7ba4b4487..09b18ca3d 100644 --- a/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py +++ b/examples/RodContactCase/RodRodContact/rod_rod_contact_parallel_validation.py @@ -35,7 +35,9 @@ class ParallelRodRodContact( shear_modulus = E / (poisson_ratio + 1.0) # Rod orientations -start = np.zeros(3,) +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)]) @@ -55,7 +57,7 @@ class ParallelRodRodContact( density, nu, E, - shear_modulus = shear_modulus, + shear_modulus=shear_modulus, ) rod_one.velocity_collection[:] += 0.05 * -normal.reshape(3, 1) @@ -89,9 +91,7 @@ class ParallelRodRodContact( # Add call backs class RodCallBack(CallBackBaseClass): - """ - - """ + """ """ def __init__(self, step_skip: int, callback_params: dict): CallBackBaseClass.__init__(self) @@ -125,7 +125,9 @@ def make_callback(self, system, time, current_step: int): ) # 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, + RodCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_rod1, ) post_processing_dict_rod2 = defaultdict( @@ -133,7 +135,9 @@ def make_callback(self, system, time, current_step: int): ) # 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, + RodCallBack, + step_skip=step_skip, + callback_params=post_processing_dict_rod2, ) parallel_rod_rod_contact_sim.finalize() From eea4a7b85cd231deee57e3323f94796084fe8487 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 15 May 2022 21:49:15 -0500 Subject: [PATCH 41/47] Changelog for upcoming release update-0.2.3 --- RELEASE.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/RELEASE.md b/RELEASE.md index 5e74ac3ec..c22480825 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -1,3 +1,15 @@ +# Release Note (version 0.2.3) + + +## Developer Note +The major updates are knot theory module added to the Cosserat rod as *mixin*, and muscular snake example is added. + +## Notable Changes +- #70: Knot theory module to compute topological quantities. +- #71: Reorganize rod constructor warning messages and collect messages in log. +- #72: Muscular snake example is added. +--- + # Release Note (version 0.2.2) ## Developer Note From 7d85dea97c5785ee8ee7436568d423298cbf1532 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Sun, 15 May 2022 21:58:52 -0500 Subject: [PATCH 42/47] black version file --- elastica/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elastica/version.py b/elastica/version.py index 826c403a7..ae50f1f69 100644 --- a/elastica/version.py +++ b/elastica/version.py @@ -1 +1 @@ -VERSION = "0.2.3" \ No newline at end of file +VERSION = "0.2.3" From dfdf497849fd0fdc51ffa249566da32c7056fe17 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Tue, 17 May 2022 11:36:15 -0500 Subject: [PATCH 43/47] typo:external_forces documentation --- docs/api/external_forces.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/api/external_forces.rst b/docs/api/external_forces.rst index 3d351abb7..fdb55b398 100644 --- a/docs/api/external_forces.rst +++ b/docs/api/external_forces.rst @@ -8,7 +8,7 @@ External Forces / Interactions Description ----------- -External force and environmental interaction are represented as force/torque bounaary condition at different location. +External force and environmental interaction are represented as force/torque boundary condition at different location. .. rubric:: Available Forcing From ee846228db21bd51e347ffe17c1d79ee7b51b356 Mon Sep 17 00:00:00 2001 From: armantekinalp Date: Tue, 17 May 2022 11:37:07 -0500 Subject: [PATCH 44/47] typo:fix binder tutorial title --- docs/guide/pybind.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/guide/pybind.md b/docs/guide/pybind.md index ca31ca154..54393f0da 100644 --- a/docs/guide/pybind.md +++ b/docs/guide/pybind.md @@ -1,4 +1,4 @@ -# PyBind Tutorials +# Binder Tutorials