From 1ae2b13b97f10bb7d241399189679e262d6ffef9 Mon Sep 17 00:00:00 2001 From: Andrew-Tao Date: Tue, 17 Dec 2024 21:02:37 -0600 Subject: [PATCH] update non-linear examples --- ...antilever_conservative_distributed_load.py | 268 ++++++++++++ ...tilever_distrubuted_load_postprecessing.py | 396 +++++++++++++++++ ...ilever_nonconservative_distributed_load.py | 288 +++++++++++++ .../Tumbling_Unconstrained_Rod.py | 213 +++++++++ ...mbling_Unconstrained_Rod_postprocessing.py | 405 ++++++++++++++++++ 5 files changed, 1570 insertions(+) create mode 100644 examples/CantileverDistributedLoad/cantilever_conservative_distributed_load.py create mode 100644 examples/CantileverDistributedLoad/cantilever_distrubuted_load_postprecessing.py create mode 100644 examples/CantileverDistributedLoad/cantilever_nonconservative_distributed_load.py create mode 100644 examples/TumblingUnconstrainedRod/Tumbling_Unconstrained_Rod.py create mode 100644 examples/TumblingUnconstrainedRod/Tumbling_Unconstrained_Rod_postprocessing.py diff --git a/examples/CantileverDistributedLoad/cantilever_conservative_distributed_load.py b/examples/CantileverDistributedLoad/cantilever_conservative_distributed_load.py new file mode 100644 index 00000000..7386a6e9 --- /dev/null +++ b/examples/CantileverDistributedLoad/cantilever_conservative_distributed_load.py @@ -0,0 +1,268 @@ +# from charset_normalizer.legacy import ResultDict +from matplotlib import pyplot as plt +import numpy as np +import elastica as ea +from cantilever_distrubuted_load_postprecessing import ( + plot_video_with_surface, + Find_Tip_Position, + adjust_square_cross_section, +) + + +def Conservative_Force_Simulator(load, Animation=False): + class StretchingBeamSimulator( + ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping, ea.CallBacks + ): + pass + + stretch_sim = StretchingBeamSimulator() + final_time = 10 + + # Options + PLOT_FIGURE = True + SAVE_FIGURE = False + SAVE_RESULTS = False + # setting up test params + n_elem = 100 + start = np.zeros((3,)) + direction = np.array([1.0, 0.0, 0.0]) + normal = np.array([0.0, 1.0, 0.0]) + base_length = 0.5 + base_radius = 0.01 / ( + np.pi ** (1 / 2) + ) # The Cross-sectional area is 1e-4(we assume its equivalent to a square cross-sectional surface with same area) + base_area = np.pi * base_radius**2 + density = 1000 # nomilized with conservative case F=15 + youngs_modulus = 1.2e7 + dl = base_length / n_elem + dt = 0.1 * dl / 50 + I = (0.01**4) / 12 + end_force_x = (youngs_modulus * I * load) / (density * base_area * (base_length**3)) + # For shear modulus of 1e4, nu is 99! + poisson_ratio = 0.0 + shear_modulus = youngs_modulus / (2 * (poisson_ratio + 1.0)) + + rendering_fps = 30 + + stretchable_rod = ea.CosseratRod.straight_rod( + n_elem, + start, + direction, + normal, + base_length, + base_radius, + density, + youngs_modulus=youngs_modulus, + shear_modulus=shear_modulus, + ) + + adjust_section = adjust_square_cross_section( + n_elem, + direction, + normal, + base_length, + base_radius, + density, + youngs_modulus=youngs_modulus, + shear_modulus=shear_modulus, + rod_origin_position=start, + ring_rod_flag=False, + ) + + stretchable_rod.mass_second_moment_of_inertia = adjust_section[0] + stretchable_rod.inv_mass_second_moment_of_inertia = adjust_section[1] + stretchable_rod.bend_matrix = adjust_section[2] + + stretch_sim.append(stretchable_rod) + stretch_sim.constrain(stretchable_rod).using( + ea.OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,) + ) + + Conservative_Load = np.array([0.0, -end_force_x, 0.0]) + + stretch_sim.add_forcing_to(stretchable_rod).using( + ea.GravityForces, acc_gravity=Conservative_Load + ) + + # add damping + + damping_constant = 0.1 + + stretch_sim.dampen(stretchable_rod).using( + ea.AnalyticalLinearDamper, + damping_constant=damping_constant, + time_step=dt, + ) + + # Add call backs + class AxialStretchingCallBack(ea.CallBackBaseClass): + def __init__(self, step_skip: int, callback_params: dict): + ea.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() + ) + self.callback_params["velocity_magnitude"].append( + ( + stretchable_rod.velocity_collection[-1][0] ** 2 + + stretchable_rod.velocity_collection[-1][1] ** 2 + + stretchable_rod.velocity_collection[-1][2] ** 2 + ) + ** 0.5 + ) + + recorded_history = ea.defaultdict(list) + stretch_sim.collect_diagnostics(stretchable_rod).using( + AxialStretchingCallBack, step_skip=200, callback_params=recorded_history + ) + + stretch_sim.finalize() + timestepper = ea.PositionVerlet() + # timestepper = PEFRL() + + total_steps = int(final_time / dt) + print(stretch_sim) + print("Total steps", total_steps) + ea.integrate(timestepper, stretch_sim, final_time, total_steps) + + relative_tip_position = np.zeros( + 2, + ) + relative_tip_position[0] = ( + Find_Tip_Position(stretchable_rod, n_elem)[0] / base_length + ) + relative_tip_position[1] = ( + -Find_Tip_Position(stretchable_rod, n_elem)[1] / base_length + ) + + print(relative_tip_position) + + if Animation: + plot_video_with_surface( + [recorded_history], + video_name="cantilever_conservative_distributed_load.mp4", + fps=rendering_fps, + step=1, + # The following parameters are optional + x_limits=(-0.0, 0.5), # Set bounds on x-axis + y_limits=(-0.5, 0.0), # Set bounds on y-axis + z_limits=(-0.0, 0.5), # Set bounds on z-axis + dpi=100, # Set the quality of the image + vis3D=True, # Turn on 3D visualization + vis2D=False, # Turn on projected (2D) visualization + ) + + relative_tip_position = np.zeros( + 2, + ) + relative_tip_position[0] = ( + Find_Tip_Position(stretchable_rod, n_elem)[0] / base_length + ) + relative_tip_position[1] = ( + -Find_Tip_Position(stretchable_rod, n_elem)[1] / base_length + ) + + print(relative_tip_position) + return relative_tip_position + + +Conservative_Force_Simulator(15, Animation=True) + +x_tip_experiment = [] +y_tip_experiment = [] +x_tip_paper = [ + 0.9912, + 0.9309, + 0.8455, + 0.7613, + 0.6874, + 0.6249, + 0.5724, + 0.5281, + 0.4906, + 0.4584, + 0.4306, + 0.4064, + 0.3851, +] +y_tip_paper = [ + 0.1241, + 0.3411, + 0.4976, + 0.6031, + 0.6745, + 0.7243, + 0.7603, + 0.7871, + 0.8077, + 0.8239, + 0.8370, + 0.8478, + 0.8568, +] +load_on_rod = np.arange(1, 26, 2) +for i in load_on_rod: + x_tip_experiment.append(Conservative_Force_Simulator(i)[0]) + y_tip_experiment.append(Conservative_Force_Simulator(i)[1]) + + +plt.plot( + load_on_rod, + x_tip_paper, + color="black", + marker="*", + linestyle="--", + label="Theoretical_x", +) +plt.plot( + load_on_rod, + y_tip_paper, + color="black", + marker="*", + linestyle=":", + label="Theoretical_y", +) +plt.scatter( + load_on_rod, + x_tip_experiment, + color="blue", + marker="s", + linestyle="None", + label="x_tip/L", +) +plt.scatter( + load_on_rod, + y_tip_experiment, + color="red", + marker="s", + linestyle="None", + label="y_tip/L", +) + +plt.title("Conservative-Load Elastica Simulation Results") +# Title +plt.xlabel("Load") # X-axis label +plt.ylabel("x_tip/L and y_tip/L") # Y-axis label +plt.grid() +plt.legend() # Optional: Add a grid +plt.show() # Display the plot diff --git a/examples/CantileverDistributedLoad/cantilever_distrubuted_load_postprecessing.py b/examples/CantileverDistributedLoad/cantilever_distrubuted_load_postprecessing.py new file mode 100644 index 00000000..210bbfe1 --- /dev/null +++ b/examples/CantileverDistributedLoad/cantilever_distrubuted_load_postprecessing.py @@ -0,0 +1,396 @@ +import numpy as np +from matplotlib import pyplot as plt +from matplotlib import cm +from tqdm import tqdm +from typing import Dict, Sequence +from typing import Optional +import logging +from numpy.testing import assert_allclose +from elastica.utils import MaxDimension, Tolerance +from elastica._linalg import _batch_cross, _batch_norm, _batch_dot +from elastica.rod.factory_function import ( + _assert_dim, + _position_validity_checker, + _directors_validity_checker, + _position_validity_checker_ring_rod, +) + + +def Find_Tip_Position(rod, n_elem): + x_tip = rod.position_collection[0][n_elem] + y_tip = rod.position_collection[1][n_elem] + z_tip = rod.position_collection[2][n_elem] + + return x_tip, y_tip, z_tip + + +def adjust_square_cross_section( + n_elements, + direction, + normal, + base_length, + base_radius, + density, + youngs_modulus: float, + *, + rod_origin_position: np.ndarray, + ring_rod_flag: bool, + shear_modulus: Optional[float] = None, + position: Optional[np.ndarray] = None, + directors: Optional[np.ndarray] = None, + **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." + ) + + # sanity checks here + assert n_elements > 2 if ring_rod_flag else n_elements > 1 + assert base_length > Tolerance.atol() + assert np.sqrt(np.dot(normal, normal)) > Tolerance.atol() + assert np.sqrt(np.dot(direction, direction)) > Tolerance.atol() + + # define the number of nodes and voronoi elements based on if rod is + # straight and open or closed and ring shaped + n_nodes = n_elements if ring_rod_flag else n_elements + 1 + n_voronoi_elements = n_elements if ring_rod_flag else n_elements - 1 + + # check if position is given. + if position is None: # Generate straight and uniform rod + # Set the position array + position = np.zeros((MaxDimension.value(), n_nodes)) + if not ring_rod_flag: # i.e. a straight open rod + + start = rod_origin_position + 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) + else: # i.e a closed ring rod + ring_center_position = rod_origin_position + binormal = np.cross(direction, normal) + for i in range(n_elements): + position[..., i] = ( + base_length + / (2 * np.pi) + * ( + np.cos(2 * np.pi / n_elements * i) * binormal + + np.sin(2 * np.pi / n_elements * i) * direction + ) + ) + ring_center_position + _position_validity_checker_ring_rod( + position, ring_center_position, n_elements + ) + + # Compute rest lengths and tangents + position_for_difference = ( + np.hstack((position, position[..., 0].reshape(3, 1))) + if ring_rod_flag + else position + ) + position_diff = position_for_difference[..., 1:] - position_for_difference[..., :-1] + rest_lengths = _batch_norm(position_diff) + tangents = position_diff / rest_lengths + normal /= np.linalg.norm(normal) + + if directors is None: # Generate straight uniform rod + print("Directors not specified") + # 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 + # directors will be wrong!! + assert_allclose( + _batch_dot(normal_collection, tangents), + 0, + atol=Tolerance.atol(), + err_msg=(" Rod normal and tangent are not perpendicular to each other!"), + ) + 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_dim(radius_temp, 2, "radius") + radius[:] = radius_temp + # Check if the elements of radius are greater than tolerance + 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_dim(density_temp, 2, "density") + density_array[:] = density_temp + # Check if the elements of density are greater than tolerance + assert np.all( + density_array > Tolerance.atol() + ), " Density has to be greater than 0." + + # Second moment of inertia + + side_length = np.zeros(n_elements) + side_length.fill(0.01) + + A0 = np.pi * radius * radius + + I0_1 = (side_length**4) / 12 + I0_2 = (side_length**4) / 12 + I0_3 = I0_2 * 2 + + I0 = np.array([I0_1, I0_2, I0_3]).transpose() + + # Mass second moment of inertia for disk cross-section + mass_second_moment_of_inertia = np.zeros( + (MaxDimension.value(), MaxDimension.value(), n_elements), np.float64 + ) + + mass_second_moment_of_inertia_temp = np.einsum( + "ij,i->ij", I0, density * rest_lengths + ) + + for i in range(n_elements): + np.fill_diagonal( + mass_second_moment_of_inertia[..., i], + mass_second_moment_of_inertia_temp[i, :], + ) + # 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." + log.warning(message) + + # Inverse of second moment of inertia + inv_mass_second_moment_of_inertia = np.zeros( + (MaxDimension.value(), MaxDimension.value(), n_elements) + ) + for i in range(n_elements): + # Check rank of mass moment of inertia matrix to see if it is invertible + assert ( + np.linalg.matrix_rank(mass_second_moment_of_inertia[..., i]) + == MaxDimension.value() + ) + inv_mass_second_moment_of_inertia[..., i] = np.linalg.inv( + mass_second_moment_of_inertia[..., i] + ) + + # Shear/Stretch matrix + if not shear_modulus: + 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 + # "On Timoshenko's correction for shear in vibrating beams" by Kaneko, 1975 + alpha_c = 27.0 / 28.0 + shear_matrix = np.zeros( + (MaxDimension.value(), MaxDimension.value(), n_elements), np.float64 + ) + for i in range(n_elements): + np.fill_diagonal( + shear_matrix[..., i], + [ + alpha_c * shear_modulus * A0[i], + alpha_c * shear_modulus * A0[i], + youngs_modulus * A0[i], + ], + ) + + # Bend/Twist matrix + bend_matrix = np.zeros( + (MaxDimension.value(), MaxDimension.value(), n_voronoi_elements + 1), np.float64 + ) + for i in range(n_elements): + np.fill_diagonal( + bend_matrix[..., i], + [ + youngs_modulus * I0_1[i], + youngs_modulus * I0_2[i], + shear_modulus * I0_3[i], + ], + ) + if ring_rod_flag: # wrap around the value in the last element + bend_matrix[..., -1] = bend_matrix[..., 0] + 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 + rest_lengths_temp_for_voronoi = ( + np.hstack((rest_lengths, rest_lengths[0])) if ring_rod_flag else rest_lengths + ) + bend_matrix = ( + bend_matrix[..., 1:] * rest_lengths_temp_for_voronoi[1:] + + bend_matrix[..., :-1] * rest_lengths_temp_for_voronoi[0:-1] + ) / (rest_lengths_temp_for_voronoi[1:] + rest_lengths_temp_for_voronoi[:-1]) + + return ( + mass_second_moment_of_inertia, + inv_mass_second_moment_of_inertia, + bend_matrix, + ) + + +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) + + ax.view_init(elev=90, azim=-90) + + 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-v0_8-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()) diff --git a/examples/CantileverDistributedLoad/cantilever_nonconservative_distributed_load.py b/examples/CantileverDistributedLoad/cantilever_nonconservative_distributed_load.py new file mode 100644 index 00000000..5fc566d2 --- /dev/null +++ b/examples/CantileverDistributedLoad/cantilever_nonconservative_distributed_load.py @@ -0,0 +1,288 @@ +from matplotlib import pyplot as plt +import numpy as np +import elastica as ea + + +from cantilever_distrubuted_load_postprecessing import ( + plot_video_with_surface, + Find_Tip_Position, + adjust_square_cross_section, +) + + +def cantilever_subjected_to_a_nonconservative_load( + n_elem, + base_length, + base_radius, + youngs_modulus, + dimentionless_varible, + Animation=False, + PLOT_FIGURE_Equilibrium=False, +): + class StretchingBeamSimulator( + ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping, ea.CallBacks + ): + pass + + stretch_sim = StretchingBeamSimulator() + + stretchable_rod = ea.CosseratRod.straight_rod( + n_elem, + start, + direction, + normal, + base_length, + base_radius, + density, + youngs_modulus=youngs_modulus, + shear_modulus=shear_modulus, + ) + + adjust_section = adjust_square_cross_section( + n_elem, + direction, + normal, + base_length, + base_radius, + density, + youngs_modulus=youngs_modulus, + shear_modulus=shear_modulus, + rod_origin_position=start, + ring_rod_flag=False, + ) + + stretchable_rod.mass_second_moment_of_inertia = adjust_section[0] + stretchable_rod.inv_mass_second_moment_of_inertia = adjust_section[1] + stretchable_rod.bend_matrix = adjust_section[2] + + stretch_sim.append(stretchable_rod) + + stretch_sim.constrain(stretchable_rod).using( + ea.OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,) + ) + + load = (youngs_modulus * I * dimentionless_varible) / ( + density * base_area * (base_length**3) + ) + + stretch_sim.add_forcing_to(stretchable_rod).using( + ea.external_forces.NonconserativeForce, load + ) + + # add damping + dl = base_length / n_elem + dt = 0.1 * dl / 50 + damping_constant = 0.2 + + stretch_sim.dampen(stretchable_rod).using( + ea.AnalyticalLinearDamper, + damping_constant=damping_constant, + time_step=dt, + ) + + # Add call backs + class AxialStretchingCallBack(ea.CallBackBaseClass): + """ + Tracks the velocity norms of the rod + """ + + def __init__(self, step_skip: int, callback_params: dict): + ea.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() + ) + self.callback_params["velocity_magnitude"].append( + ( + stretchable_rod.velocity_collection[-1][0] ** 2 + + stretchable_rod.velocity_collection[-1][1] ** 2 + + stretchable_rod.velocity_collection[-1][2] ** 2 + ) + ** 0.5 + ) + + recorded_history = ea.defaultdict(list) + + stretch_sim.collect_diagnostics(stretchable_rod).using( + AxialStretchingCallBack, step_skip=200, callback_params=recorded_history + ) + + stretch_sim.finalize() + timestepper = ea.PositionVerlet() + # timestepper = PEFRL() + + total_steps = int(final_time / dt) + print(stretch_sim) + print("Total steps", total_steps) + ea.integrate(timestepper, stretch_sim, final_time, total_steps) + + if PLOT_FIGURE_Equilibrium: + + plt.plot( + recorded_history["time"], + recorded_history["velocity_magnitude"], + lw=1.0, + label="velocity_magnitude", + ) + plt.xlabel("Time (s)") + plt.ylabel("Speed (m/s)") + plt.title("Tip Position and Speed vs. Time") + plt.legend() + plt.grid() + plt.show() + + rendering_fps = 30 + if Animation: + plot_video_with_surface( + [recorded_history], + video_name="cantilever_Non-conservative_distributed_load.mp4", + fps=rendering_fps, + step=1, + # The following parameters are optional + x_limits=(-0.0, 0.5), # Set bounds on x-axis + y_limits=(-0.5, 0.0), # Set bounds on y-axis + z_limits=(-0.0, 0.5), # Set bounds on z-axis + dpi=100, # Set the quality of the image + vis3D=True, # Turn on 3D visualization + vis2D=False, # Turn on projected (2D) visualization + ) + + pos = stretchable_rod.position_collection.view() + + tip_position = Find_Tip_Position(stretchable_rod, n_elem) + relative_tip_position = np.zeros((2,)) + + relative_tip_position[0] = tip_position[0] / base_length + relative_tip_position[1] = -tip_position[1] / base_length + + print(relative_tip_position) + return relative_tip_position + + +final_time = 10 +# setting up test params +n_elem = 100 +start = np.zeros((3,)) +direction = np.array([1.0, 0.0, 0.0]) +normal = np.array([0.0, 1.0, 0.0]) +base_length = 0.5 +base_radius = 0.01 / (np.pi ** (1 / 2)) +base_area = np.pi * base_radius**2 +density = 1000 +dimentionless_varible = 15 +youngs_modulus = 1.2e7 +# For shear modulus of 1e4, nu is 99! +poisson_ratio = 0 +shear_modulus = youngs_modulus / (2 * (poisson_ratio + 1.0)) +I = (0.01**4) / 12 + + +cantilever_subjected_to_a_nonconservative_load( + n_elem, base_length, base_radius, youngs_modulus, -15, True, False +) + +x_tip_paper = [ + 0.9910, + 0.9213, + 0.7935, + 0.6265, + 0.4422, + 0.2596, + 0.0705, + -0.0585, + -0.1677, + -0.2585, + -0.3258, + -0.3731, + -0.4042, +] +y_tip_paper = [ + 0.1251, + 0.3622, + 0.5628, + 0.7121, + 0.8049, + 0.8451, + 0.8420, + 0.8071, + 0.7513, + 0.6840, + 0.6123, + 0.5411, + 0.4734, +] +x_tip_experiment = [] +y_tip_experiment = [] +load_on_rod = np.arange(1, 26, 2) +for i in load_on_rod: + x_tip_experiment.append( + cantilever_subjected_to_a_nonconservative_load( + n_elem, base_length, base_radius, youngs_modulus, i, False, False + )[0] + ) + y_tip_experiment.append( + -cantilever_subjected_to_a_nonconservative_load( + n_elem, base_length, base_radius, youngs_modulus, i, False, False + )[1] + ) + +plt.plot( + load_on_rod, + x_tip_paper, + color="black", + marker="*", + linestyle="--", + label="Theoretical_x", +) +plt.plot( + load_on_rod, + y_tip_paper, + color="black", + marker="*", + linestyle=":", + label="Theoretical_y", +) +plt.scatter( + load_on_rod, + x_tip_experiment, + color="blue", + marker="s", + linestyle="None", + label="x_tip/L", +) +plt.scatter( + load_on_rod, + y_tip_experiment, + color="red", + marker="s", + linestyle="None", + label="y_tip/L", +) + +plt.title("Non-Conservative-Load Elastica Simulation Results") +# Title +plt.xlabel("Load") # X-axis label +plt.ylabel("x_tip/L and y_tip/L") # Y-axis label +plt.grid() +plt.legend() # Optional: Add a grid +plt.show() # Display the plot diff --git a/examples/TumblingUnconstrainedRod/Tumbling_Unconstrained_Rod.py b/examples/TumblingUnconstrainedRod/Tumbling_Unconstrained_Rod.py new file mode 100644 index 00000000..27c215cf --- /dev/null +++ b/examples/TumblingUnconstrainedRod/Tumbling_Unconstrained_Rod.py @@ -0,0 +1,213 @@ +import elastica as ea + +import numpy as np +from elastica.timestepper.symplectic_steppers import PositionVerlet +from elastica.timestepper import integrate +from elastica.external_forces import ( + EndPointTorque, + EndpointForces_with_time_factor, + EndPointTorque_with_time_factor, +) +from elastica.external_forces import UniformTorques +from Tumbling_Unconstrained_Rod_postprocessing import ( + plot_video_with_surface, + adjust_parameter, + lamda_t_function, +) +from matplotlib import pyplot as plt + +n_elem = 256 +start = np.array([0.0, 0.0, 8.0]) +end = np.array([6.0, 0.0, 0.0]) +direction = np.array([0.6, 0.0, -0.8]) +normal = np.array([0.0, 1.0, 0.0]) +base_length = 10 + +side_length = 0.01 + + +base_radius = side_length / (np.pi ** (1 / 2)) + + +density = 1e4 +youngs_modulus = 1e7 +poisson_ratio = 0 +shear_modulus = youngs_modulus / (poisson_ratio + 1.0) + + +class NonConstrainRodSimulator( + ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping, ea.CallBacks +): + pass + + +sim = NonConstrainRodSimulator() + +rod = ea.CosseratRod.straight_rod( + n_elem, + start, + direction, + normal, + base_length, + base_radius, + density, + youngs_modulus=youngs_modulus, +) + +adjust_section = adjust_parameter( + n_elem, + direction, + normal, + base_length, + side_length, + base_radius, + density, + youngs_modulus=youngs_modulus, + rod_origin_position=start, + ring_rod_flag=False, +) + +rod.mass_second_moment_of_inertia = adjust_section[0] +rod.inv_mass_second_moment_of_inertia = adjust_section[1] +rod.bend_matrix = adjust_section[2] + +print("mass_second_moment_of_inertia=", rod.mass_second_moment_of_inertia) +print("bend_matrix=", rod.bend_matrix) + +sim.append(rod) + +dl = base_length / n_elem +dt = 0.01 * dl / 1 + +origin_force = np.array([0.0, 0.0, 0.0]) +end_force = np.array([20.0, 0.0, 0.0]) + +sim.add_forcing_to(rod).using( + EndpointForces_with_time_factor, origin_force, end_force, lamda_t_function +) + + +sim.add_forcing_to(rod).using( + EndPointTorque_with_time_factor, + 1, + lamda_t_function, + direction=np.array([0.0, 200.0, -100.0]), +) + +sim.dampen(rod).using( + ea.AnalyticalLinearDamper, + damping_constant=0.0, + time_step=dt, +) + +print("Forces added to the rod") + +final_time = 20 +total_steps = int(final_time / dt) + +print("Total steps to take", total_steps) + +rendering_fps = 30 +step_skip = int(1.0 / (rendering_fps * dt)) + + +class AxialStretchingCallBack(ea.CallBackBaseClass): + def __init__(self, step_skip: int, callback_params: dict): + ea.CallBackBaseClass.__init__(self) + self.every = step_skip + self.callback_params = callback_params + + def make_callback(self, system, time, current_step: int): + if current_step % self.every == 0: + self.callback_params["time"].append(time) + self.callback_params["step"].append(current_step) + self.callback_params["position"].append(system.position_collection.copy()) + self.callback_params["radius"].append(system.radius.copy()) + self.callback_params["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() + ) + + +recorded_history = ea.defaultdict(list) +sim.collect_diagnostics(rod).using( + AxialStretchingCallBack, step_skip=step_skip, callback_params=recorded_history +) + +sim.finalize() +print("System finalized") + + +timestepper = PositionVerlet() +integrate(timestepper, sim, final_time, total_steps) + +time_analytic = [0.0, 2.0, 3.0, 3.8, 4.4, 5.0, 5.5, 5.8, 6.1, 6.5, 7.0] +mass_center_analytic = [ + [ + 3.0, + 4.0667, + 6.5667, + 9.7304, + 12.5288, + 15.500, + 18.000, + 19.500, + 21.00, + 23.00, + 25.500, + ], + [4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], +] +plt.plot( + time_analytic, + mass_center_analytic[0], + marker="*", + color="black", + label="x_analytic", +) +plt.plot( + time_analytic, + mass_center_analytic[1], + marker="*", + color="black", + label="y_analytic", +) +plt.plot( + time_analytic, + mass_center_analytic[2], + marker="*", + color="black", + label="z_analytic", +) + + +mass_center = np.array(recorded_history["center_of_mass"]) + +plt.plot(recorded_history["time"][0:240], mass_center[:, 0][0:240], label="x") +plt.plot(recorded_history["time"][0:240], mass_center[:, 1][0:240], label="y") +plt.plot(recorded_history["time"][0:240], mass_center[:, 2][0:240], label="z") + +plt.xlabel("Time/(second)") # X-axis label +plt.ylabel("Center of mass") # Y-axis label +plt.grid() +plt.legend() # Optional: Add a grid +plt.show() + + +plot_video_with_surface( + [recorded_history], + video_name="Tumbling_Unconstrained_Rod.mp4", + fps=rendering_fps, + step=1, + # The following parameters are optional + x_limits=(0, 200), # Set bounds on x-axis + y_limits=(-4, 4), # Set bounds on y-axis + z_limits=(0.0, 8), # Set bounds on z-axis + dpi=100, # Set the quality of the image + vis3D=True, # Turn on 3D visualization + vis2D=False, # Turn on projected (2D) visualization +) diff --git a/examples/TumblingUnconstrainedRod/Tumbling_Unconstrained_Rod_postprocessing.py b/examples/TumblingUnconstrainedRod/Tumbling_Unconstrained_Rod_postprocessing.py new file mode 100644 index 00000000..7988d26d --- /dev/null +++ b/examples/TumblingUnconstrainedRod/Tumbling_Unconstrained_Rod_postprocessing.py @@ -0,0 +1,405 @@ +from typing import Optional +import logging + +from numpy.testing import assert_allclose +from elastica.utils import MaxDimension, Tolerance +from elastica._linalg import _batch_cross, _batch_norm, _batch_dot +from elastica.rod.factory_function import ( + _assert_dim, + _position_validity_checker, + _directors_validity_checker, + _position_validity_checker_ring_rod, +) +import numpy as np +from matplotlib import pyplot as plt + +from matplotlib import cm +from tqdm import tqdm + +from typing import Dict, Sequence + + +def lamda_t_function(time): + if time < 2.5: + factor = time * (1 / 2.5) + elif time > 2.5 and time < 5.0: + factor = -time * (1 / 2.5) + 2 + else: + factor = 0 + + return factor + + +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) + + ax.view_init(elev=0, azim=0) + + time_idx = 0 + rod_lines = [None for _ in range(n_visualized_rods)] + rod_com_lines = [None for _ in range(n_visualized_rods)] + rod_scatters = [None for _ in range(n_visualized_rods)] + + for rod_idx in range(n_visualized_rods): + inst_position, inst_radius = rod_history_unpacker(rod_idx, time_idx) + 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-v0_8-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( + 1000 * 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()) + + +def adjust_parameter( + n_elements, + direction, + normal, + base_length, + length, + base_radius, + density, + youngs_modulus: float, + *, + rod_origin_position: np.ndarray, + ring_rod_flag: bool, + shear_modulus: Optional[float] = None, + position: Optional[np.ndarray] = None, + directors: Optional[np.ndarray] = None, + rest_sigma: Optional[np.ndarray] = None, + rest_kappa: Optional[np.ndarray] = None, + **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." + ) + + # sanity checks here + assert n_elements > 2 if ring_rod_flag else n_elements > 1 + assert base_length > Tolerance.atol() + assert np.sqrt(np.dot(normal, normal)) > Tolerance.atol() + assert np.sqrt(np.dot(direction, direction)) > Tolerance.atol() + + # define the number of nodes and voronoi elements based on if rod is + # straight and open or closed and ring shaped + n_nodes = n_elements if ring_rod_flag else n_elements + 1 + n_voronoi_elements = n_elements if ring_rod_flag else n_elements - 1 + + # check if position is given. + if position is None: # Generate straight and uniform rod + # Set the position array + position = np.zeros((MaxDimension.value(), n_nodes)) + if not ring_rod_flag: # i.e. a straight open rod + + start = rod_origin_position + 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) + else: # i.e a closed ring rod + ring_center_position = rod_origin_position + binormal = np.cross(direction, normal) + for i in range(n_elements): + position[..., i] = ( + base_length + / (2 * np.pi) + * ( + np.cos(2 * np.pi / n_elements * i) * binormal + + np.sin(2 * np.pi / n_elements * i) * direction + ) + ) + ring_center_position + _position_validity_checker_ring_rod( + position, ring_center_position, n_elements + ) + + # Compute rest lengths and tangents + position_for_difference = ( + np.hstack((position, position[..., 0].reshape(3, 1))) + if ring_rod_flag + else position + ) + position_diff = position_for_difference[..., 1:] - position_for_difference[..., :-1] + rest_lengths = _batch_norm(position_diff) + tangents = position_diff / rest_lengths + normal /= np.linalg.norm(normal) + + if directors is None: # Generate straight uniform rod + print("Directors not specified") + # 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 + # directors will be wrong!! + assert_allclose( + _batch_dot(normal_collection, tangents), + 0, + atol=Tolerance.atol(), + err_msg=(" Rod normal and tangent are not perpendicular to each other!"), + ) + 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_dim(radius_temp, 2, "radius") + radius[:] = radius_temp + # Check if the elements of radius are greater than tolerance + 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_dim(density_temp, 2, "density") + density_array[:] = density_temp + # Check if the elements of density are greater than tolerance + assert np.all( + density_array > Tolerance.atol() + ), " Density has to be greater than 0." + + # Second moment of inertia + + side_length = np.zeros(n_elements) + side_length.fill(length) + + A0 = np.pi * radius * radius + + I0_1 = ((side_length**4) / 12) * (1200000) + I0_2 = ((side_length**4) / 12) * (1200000) + I0_3 = I0_2 * 2 + + I0 = np.array([I0_1, I0_2, I0_3]).transpose() + + # Mass second moment of inertia for disk cross-section + mass_second_moment_of_inertia = np.zeros( + (MaxDimension.value(), MaxDimension.value(), n_elements), np.float64 + ) + + mass_second_moment_of_inertia_temp = np.einsum( + "ij,i->ij", I0, density * rest_lengths + ) + + for i in range(n_elements): + np.fill_diagonal( + mass_second_moment_of_inertia[..., i], + mass_second_moment_of_inertia_temp[i, :], + ) + # 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." + log.warning(message) + + # Inverse of second moment of inertia + inv_mass_second_moment_of_inertia = np.zeros( + (MaxDimension.value(), MaxDimension.value(), n_elements) + ) + for i in range(n_elements): + # Check rank of mass moment of inertia matrix to see if it is invertible + assert ( + np.linalg.matrix_rank(mass_second_moment_of_inertia[..., i]) + == MaxDimension.value() + ) + inv_mass_second_moment_of_inertia[..., i] = np.linalg.inv( + mass_second_moment_of_inertia[..., i] + ) + + # Shear/Stretch matrix + if not shear_modulus: + 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 + # "On Timoshenko's correction for shear in vibrating beams" by Kaneko, 1975 + alpha_c = 27.0 / 28.0 + shear_matrix = np.zeros( + (MaxDimension.value(), MaxDimension.value(), n_elements), np.float64 + ) + for i in range(n_elements): + np.fill_diagonal( + shear_matrix[..., i], + [ + alpha_c * shear_modulus * A0[i], + alpha_c * shear_modulus * A0[i], + youngs_modulus * A0[i] * 100, + ], + ) + + # Bend/Twist matrix + bend_matrix = np.zeros( + (MaxDimension.value(), MaxDimension.value(), n_voronoi_elements + 1), np.float64 + ) + for i in range(n_elements): + np.fill_diagonal( + bend_matrix[..., i], + [ + youngs_modulus * I0_1[i] * (1 / 20), + youngs_modulus * I0_2[i] * (1 / 20), + shear_modulus * I0_3[i], + ], + ) + if ring_rod_flag: # wrap around the value in the last element + bend_matrix[..., -1] = bend_matrix[..., 0] + 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 + rest_lengths_temp_for_voronoi = ( + np.hstack((rest_lengths, rest_lengths[0])) if ring_rod_flag else rest_lengths + ) + bend_matrix = ( + bend_matrix[..., 1:] * rest_lengths_temp_for_voronoi[1:] + + bend_matrix[..., :-1] * rest_lengths_temp_for_voronoi[0:-1] + ) / (rest_lengths_temp_for_voronoi[1:] + rest_lengths_temp_for_voronoi[:-1]) + + return ( + mass_second_moment_of_inertia, + inv_mass_second_moment_of_inertia, + bend_matrix, + )