How to simulate the gravity of cable being applied forces? about pyelastica HOT 11 CLOSED

import numpy as np

from elastica.wrappers import BaseSystemCollection, Connections, Constraints, Forcing, CallBacks
from elastica.rod.cosserat_rod import CosseratRod
from elastica.boundary_conditions import FreeRod
from elastica.timestepper.symplectic_steppers import PositionVerlet
from elastica.timestepper import integrate
from elastica import *

# Make simulator
class SutureSimulator(BaseSystemCollection, Connections, Constraints, Forcing, CallBacks):
    pass

# Make constraint to apply velocity to rod midpoint while both ends are held in place
class SutureTest(FreeRod):
    def __init__(self):
        pass

    def constrain_values(self, system, time: np.float64 = 0.0):
        # Anchor both rod ends
        system.position_collection[..., 0] = np.zeros(3)
        # system.position_collection[..., -1] = [0, 0.0, 5.0] # This position is in cm units

    def constrain_rates(self, system, time: np.float64 = 0.0):
        # # Calculate displacement from goal position (offset from midpoint by 8 mm in y direction)
        # displacement = (
        #     np.array([5.0, 0.0, 0.8]) # Goal point, in cm units
        #     - system.position_collection[..., n_elem//2] # Rod midpoint
        # )
        # # Move at low speed in direction of displacement
        # # Stop movement if control node is close enought to goal position
        # if (np.linalg.norm(displacement) > 1e-6):
        #     system.velocity_collection[..., n_elem//2] = displacement / 2 # arbitrarily scale displacement to reduce speed
        # else:
        #     system.velocity_collection[..., n_elem//2] *= 0

        # Anchor both rod ends
        system.velocity_collection[..., 0] *= 0
        # system.velocity_collection[..., -1] *= 0


class MoveMidPointByApplyingForces(NoForces):

    def __init__(self, k, ramp_up_time, target_mid_point_position):
        self.k = k
        self.ramp_up_time = ramp_up_time
        self.target_mid_point_position = target_mid_point_position

    def apply_forces(self, system, time: np.float = 0.0):
        factor = min(time/self.ramp_up_time, 1.0)

        displacement = self.target_mid_point_position -system.position_collection[...,n_elem//2]

        force = self.k * factor * displacement

        system.external_forces[...,n_elem//2] += force



class ContinuumSnakeCallBack(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["position"].append(
                system.position_collection.copy()
            )

            return

suture_sim = SutureSimulator()

#initialize all relevant rod and time variables
n_elem = 50
start = np.array([0.0, 0.0, 0.0])
direction = np.array([0.0, 0.0,1.0])
normal = np.array([0.0, 1.0, 0.0])
# Below values are in units of centimeters (cm) and grams (g)
base_length = 5 # cm
base_radius = 0.2 # cm
density = 3e-3 # g/cm^3
nu = 0.1 # g/s
E = 7.5e3 # g/(cm * s^2)
youngs_module = 0.03


final_time = 1
dl = base_length / n_elem
dt = 1e-4 * dl
total_steps = int(final_time / dt)

timestepper = PositionVerlet()

suture_thread = CosseratRod.straight_rod(
    n_elem,
    start,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    nu,
    E,
    youngs_module
)
suture_sim.append(suture_thread)

suture_sim.constrain(suture_thread).using(
    SutureTest
)

suture_sim.add_forcing_to(suture_thread).using(MoveMidPointByApplyingForces, k=8000, ramp_up_time=1.0, target_mid_point_position=np.array([0, 0.8, 2.5]))
gravitational_acc = -9.80665
suture_sim.add_forcing_to(suture_thread).using(
        GravityForces, acc_gravity=np.array([0.0, gravitational_acc, 0.0])
    )


pp_list = defaultdict(list)
suture_sim.collect_diagnostics(suture_thread).using(ContinuumSnakeCallBack, step_skip=200, callback_params=pp_list
    )


suture_sim.finalize()
integrate(timestepper, suture_sim, final_time, total_steps)


SAVE_RESULTS = True

if SAVE_RESULTS:
    import pickle

    filename = "suture_thread.dat"
    file = open(filename, "wb")
    pickle.dump(pp_list, file)
    file.close()`

class CustomForce(NoForces):

    def __init__(self, k, ramp_up_time, node):
        self.k = k
        self.node = node
        self.ramp_up_time = ramp_up_time

    def apply_forces(self, system, time: np.float = 0.0):
        factor = min(time/self.ramp_up_time, 1.0)

        if time > 0.5:
            force = self.k * factor * np.array([0,0,1])
        else:
            force = self.k * factor * np.array([1,0,0])
        system.external_forces[...,self.node] += force

suture_thread = CosseratRod.straight_rod(
    n_elem,
    start,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    nu,
    E=E_stretch,
    shear_modulus=shear_modulus,
)
I = np.pi/4 * base_radius**4
suture_thread.bend_matrix[0,0,:] = I * E_bending
suture_thread.bend_matrix[1,1,:] = I * E_bending
suture_sim.append(suture_thread)