armantekinalp
commented
3 years ago Hi @neelay-j
I would suggest couple of things:
- You have a stiff rod (E=750MPa) so it is expected that to use lower time-steps. Here I would suggest sometimes tuning
nucan help. - Your boundary condition is changing middle node velocity to move it towards target position. This might be tricky and I think it is not well defined boundary condition. Instead I suggest you to assume a spring between current node position and target node position and apply forces to the node. I am attaching an example code for you. Note that for fast convergence I change the material properties so you need to play around with
k,dtetc.
Quick note: PyElastica assumes that rods are incompressible so setting Poisson Ratio some value different than 0.5 will only change the shear modulus of the rod.
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] = [10.0, 0.0, 0.0] # This position is in cm units
return
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
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([1.0, 0.0, 0.0])
normal = np.array([0.0, 1.0, 0.0])
# Below values are in units of centimeters (cm) and grams (g)
base_length = 10 # cm
base_radius = 0.2 # cm
density = 3e-3 # g/cm^3
nu = 0.1 # g/s
E = 7.5e3 # g/(cm * s^2)
poisson_ratio = 2.0
final_time = 2.0
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,
poisson_ratio
)
suture_sim.append(suture_thread)
suture_sim.constrain(suture_thread).using(
SutureTest
)
suture_sim.add_forcing_to(suture_thread).using(MoveMidPointByApplyingForces, k=5000, ramp_up_time=1.0, target_mid_point_position=np.array([5.0, 0.0, 0.8]) )
suture_sim.finalize()
integrate(timestepper, suture_sim, final_time, total_steps)
Hello,
I am using PyElastica to simulate surgical suture thread. To test if my simulation works, I am applying a small velocity to a single node at the middle of the simulated thread until the node reaches a goal point, which is a small distance away from the node’s rest position. Both ends of the simulated suture are held in place. However, often when I run the simulation, the PyElastica rod’s velocities and internal forces become incredibly large very quickly and generate NaN values that destroy the simulation. The suture thread I’m simulating is mostly made of collagen and has the following properties:
Length: 10 cm Radius: 0.2 mm Density: 3e-3 g/cm^3 Young’s Modulus: 750 MPa Poisson Ratio: 2
I ran multiple simulations using different units (g and mm, g and dm, g and cm) and timesteps (1e-4 dl, 1e-6 dl, 1e-8 dl), and most became unstable within the first few iterations. However, in the following simulation, the rod seems to be stable. Unfortunately, this simulation requires an incredibly small timestep (1e-8 dl), resulting in the simulation taking several dozen hours to complete 2 seconds of simulation time.
What is the best way to make the simulation more stable while keeping its timestep at a reasonable size? Ideally, I want the simulation to complete in at most a few minutes. Thank you for your help!
Sincerely, Neelay Joglekar