Applying a constant torque in one direction about pyelastica HOT 12 CLOSED

class EndPointTorque(NoForces):
    def __init__(self,torque,direction):
        self.torque= (torque * direction).reshape(3, 1) # defined in global frame, shape(3,1)
    def apply_torques(self, system, time: np.float = 0.0):
        system.external_torques[..., -1] += system.director_collection[..., -1] @ self.torque

import numpy as np
# Import Wrappers 
from elastica.wrappers import BaseSystemCollection, Connections, Constraints, Forcing, CallBacks
# Import Cosserat Rod Class
from elastica.rod.cosserat_rod import CosseratRod
# Import Boundary Condition Classes
from elastica.boundary_conditions import OneEndFixedRod, FreeRod
from elastica.external_forces import NoForces, EndpointForces, UniformTorques, UniformForces
# Import Timestepping Functions
from elastica.timestepper.symplectic_steppers import PositionVerlet
from elastica.timestepper import integrate
from elastica.callback_functions import CallBackBaseClass
from collections import defaultdict
# from elastica import endpointtorque


class TimoshenkoBeamSimulator(BaseSystemCollection, Connections, Constraints, Forcing, CallBacks):
    pass

timoshenko_sim = TimoshenkoBeamSimulator()

# setting up test params
n_elem = 100

density = 1000
nu = 0.1
E = 1e6
# For shear modulus of 1e4, nu is 99!
poisson_ratio = 99


start = np.zeros((3,))
direction = np.array([0.0, 0.0, 1.0])
normal = np.array([0.0, 1.0, 0.0])
base_length = 3.0
base_radius = 0.25
base_area = np.pi * base_radius ** 2

shearable_rod = CosseratRod.straight_rod(
    n_elem,
    start,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    nu,
    E,
    poisson_ratio,
)

timoshenko_sim.append(shearable_rod)

timoshenko_sim.constrain(shearable_rod).using(
    OneEndFixedRod, 
    constrained_position_idx=(0,), 
    constrained_director_idx=(0,)
)
print('One end of the rod is now fixed in place')

origin_force = np.array([0.0, 0.0, 0.0])
end_force = np.array([-10.0, 0.0, 0.0])
ramp_up_time = 5.0

timoshenko_sim.add_forcing_to(shearable_rod).using(
    EndpointForces, 
    origin_force, 
    end_force, 
    ramp_up_time=ramp_up_time
)

class EndPointTorque(NoForces):
    def __init__(self,torque,direction):
        self.torque= torque * direction # defined in global frame, shape(3,1)
    def apply_torques(self, system, time: np.float = 0.0):
        system.external_torques[..., -1] += system.director_collection[..., -1] @ self.torque

direction1 = np.array([1.0, 0.0, 0.0])
toque_end = -10
timoshenko_sim.add_forcing_to(shearable_rod).using(
EndPointTorque,torque=toque_end,direction=direction1)

print('Forces added to the rod')

class TimCallBack(CallBackBaseClass):
    """
    Call back function for continuum snake
    """

    def __init__(self, step_skip: int, callback_params: dict):
        CallBackBaseClass.__init__(self)
        self.every = step_skip
        self.callback_params = callback_params

    def make_callback(self, system, time, current_step: int):

        if current_step % self.every == 0:

            self.callback_params["time"].append(time)
            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()
            )
            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()
            )

            return

pp_list = defaultdict(list)
timoshenko_sim.collect_diagnostics(shearable_rod).using(
    TimCallBack, step_skip=200, callback_params=pp_list
)

print('Callback function added to the simulator')

timoshenko_sim.finalize()
print('System finalized')

final_time = 20.0
dl = base_length / n_elem
dt = 0.01 * dl
total_steps = int(final_time / dt)
print("Total steps to take", total_steps)

timestepper = PositionVerlet()
integrate(timestepper, timoshenko_sim, final_time, total_steps)

def analytical_result(arg_rod, arg_end_force, shearing=True, n_elem=500):
    base_length = np.sum(arg_rod.rest_lengths)
    arg_s = np.linspace(0.0, base_length, n_elem)
    if type(arg_end_force) is np.ndarray:
        acting_force = arg_end_force[np.nonzero(arg_end_force)]
    else:
        acting_force = arg_end_force
    acting_force = np.abs(acting_force)
    linear_prefactor    = -acting_force / arg_rod.shear_matrix[0, 0, 0]
    quadratic_prefactor = -acting_force / 2.0 * np.sum(arg_rod.rest_lengths / arg_rod.bend_matrix[0, 0, 0])
    cubic_prefactor     = (acting_force / 6.0) / arg_rod.bend_matrix[0, 0, 0] 
    if shearing:    
        return arg_s, arg_s*linear_prefactor + arg_s**2*quadratic_prefactor + arg_s**3*cubic_prefactor
    else:
        return arg_s, arg_s**2 * quadratic_prefactor + arg_s**3 * cubic_prefactor

def plot_timoshenko(shearable_rod, end_force):
    import matplotlib.pyplot as plt
    analytical_shearable_positon = analytical_result(shearable_rod, end_force, shearing=True)
   

    fig = plt.figure(figsize=(5, 4), frameon=True, dpi=150)
    ax = fig.add_subplot(111)
    ax.grid(b=True, which="major", color="grey", linestyle="-", linewidth = 0.25)

    ax.plot(analytical_shearable_positon[0],   analytical_shearable_positon[1],  "k--", label="Timoshenko")
  
    ax.plot(shearable_rod.position_collection[2, :], 
            shearable_rod.position_collection[0, :], 
            "b-", label="n="+str(shearable_rod.n_elems))
    

    
    ax.legend(prop={"size": 12})
    ax.set_ylabel('Y Position (m)', fontsize = 12)
    ax.set_xlabel('X Position (m)', fontsize = 12)
    plt.show()

plot_timoshenko(shearable_rod, end_force)
    
import matplotlib.pyplot as plt
from IPython import display

evolve_for_time = 20.0
update_interval = 1.0e-1

# update the plot every 1 second

from IPython.display import Video

def plot_video( plot_params: dict, video_name="video.mp4", margin=0.2, fps=15):  
    from matplotlib import pyplot as plt
    import matplotlib.animation as manimation

    positions_over_time = np.array(plot_params["position"])

    print("creating video -- this can take a few minutes")
    FFMpegWriter = manimation.writers["ffmpeg"]
    metadata = dict(title="Movie Test", artist="Matplotlib", comment="Movie support!")
    writer = FFMpegWriter(fps=fps, metadata=metadata)
    fig = plt.figure()
    plt.axis("equal")
    with writer.saving(fig, video_name, dpi=100):
        for time in range(1, len(plot_params["time"])):
            
            x = positions_over_time[time][2]
         
            y = positions_over_time[time][0]
            fig.clf()
            plt.plot(x, y, "-", linewidth=3)
            
            plt.xlim([0 - margin, 3 + margin])
            plt.ylim([-0.05 - margin, 0.01 + margin])
            writer.grab_frame()
    plt.close(fig)

filename_video = "continuum100.mp4"
plot_video(pp_list, video_name=filename_video, margin=0.02, fps=25)
    
Video("continuum100.mp4")

print(shearable_rod.kappa[2])