Create a pre-curved rod

[emersonkt]

Hello PyElastica,

I am a student and I am trying to apply the library 'PyElastica' in my project for a rod with a pre-curvature. But I could only create a straight rod. Could you please advise me if there is a way to make a pre-curved rod?

Thank you in advance for your help and congratulations for this beautiful work !

Émerson

[armantekinalp]

Hi @emersonkt

There are couple of things you need to do before creating a curved rod:

1. You need to create node positions as an array of (3,n_elem+1) which contains the x,y,z coordinates (row wise) of the nodes. Note that number of nodes is number of element +1 .
2. You need to compute tangents vector. Since you have positions vector already (in step 1), procedure is straight forward.

tangents = position[…,1:] - position[…,:-1]

Tangents are unit vectors so we should normalize them by dividing its length. In PyElastica we have a build in function to compute lengths of vectors. Import the following;

from elastica._linalg import _batch_norm

tangents /= _batch_norm(tangents)

3. The last array you need is for directors. Shape of directors is (3,3,n_elem). We call the rows of director as d1, d2, and d3. You have already computed the d3 in step 2 (d3=tangents). Now you only need to compute d1 or d2. For example if you know d1 then d2 is cross product of d3 and d1.

At this point you know position and directors of the rod. Next step is to create the rod and send position and directors as kwargs;

shearable_rod = CosseratRod.straight_rod (n_elem, 
                                          start, 
                                          direction, 
                                          normal, 
                                          base_length, 
                                          base_radius, 
                                          density, 
                                          nu, 
                                          youngs_modulus, 
                                          poisson_ratio=0.5, 
                                          position=position, 
                                          directors=directors)

Since rod is not straight and curved initially we should compute the rest curvature and strain. However PyElastica has already computed curvature and strain for you when shearable_rod created. You just need to update rest kappa and rest sigma as follows.

shearable_rod.rest_kappa[:] = shearable_rod.kappa[:]
shearable_rod.rest_sigma[:] = shearable_rod.sigma[:]

Now you have a curved rod.

[emersonkt]

Thank you very much ! It worked well and I got the desired pre-curved rod.

In the other hand, when a force is defined in the end tip of the robot, I got a movement just in the last element of the rod as the attached image. Why are the elements not connected ?

Here below I send you the code as I typed it.

Code in Python

# Create rod
n_elem = 100                           # number of elements
start = np.array([0.0, 0.0, 0.0])      # Starting position of first node in rod
direction = np.array([0.0, 0, 1.0])    # Direction the rod extends
normal = np.array([0.0, 1.0, 0.0])     # normal vector of rod
base_length = 3                        # original length of rod (m)
base_radius = 25e-2                    # original radius of rod (m)
density = 1e3                          # density of rod (kg/m^3)
nu = 1e-1                              # Energy dissipation of rod
E = 1e6                                # Elastic Modulus (Pa)
poisson_ratio = 0.495                  # Poisson Ratio

# Definition of the nodes position
position = np.zeros((n_elem+1,3))
for i in range(0, n_elem+1):
  position[i][0] = (i**2)/100000
  position[i][1] = 0
  position[i][2] = (base_length/n_elem)*i

position = np.transpose(position)


# Definition of the tangent of the rod
tangents = np.zeros((3,n_elem))

for i in range(0, 3):
  tangents[i] = position[i,1:] - position[i,:-1]

tangents /= _batch_norm(tangents)
d3 = tangents


# Definition of the other axis
d2 = np.zeros((n_elem,3))
d1 = np.zeros((n_elem,3))
for i in range(0,n_elem):
  d2[i][0] = 0
  d2[i][1] = 1
  d2[i][2] = 0
  d1[i] = np.cross(d2[i],np.transpose(tangents)[i])

d2 = np.transpose(d2)
d1 = np.transpose(d1)


# Putting all direction in the director matrix
directors = np.zeros((3,3,n_elem))
directors[0] = d1
directors[1] = d2
directors[2] = d3

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

# Compute the rest curvature and strain
rod.rest_kappa[:] = rod.kappa[:]
rod.rest_sigma[:] = rod.sigma[:]


# Add rod to SystemSimulator
simulation_teste.append(rod)

[armantekinalp]

Hi @emersonkt,

Can you also provide how you provide your forcing class? Also did you see any numerical instabilities in your simulation i.e. NaN in position, velocity arrays?

[emersonkt]

Here below the forcing class:

#Define 1x3 array of the applied forces
origin_force = np.array([0.0, 0.0, 0.0])
end_force = np.array([-0.03, 0.0, 0.0]) 
final_time = 10
ramp_up_time=final_time / 2

simulation_teste.add_forcing_to(rod).using(
    EndpointForces,                 # Traction BC being applied
    origin_force,                   # Force vector applied at first node
    end_force,                      # Force vector applied at last node
    ramp_up_time                    # Ramp up time 
)

Looking at the position and velocity arrays, I got the following results. So maybe the velocities are good! However, just on the last element of the x axis I am getting an ununderstandable result.

Result from "print(rod.position_collection[0,:])" :
[ 0.00000000e+00  1.00000000e-05  4.00000000e-05  9.00000000e-05
  1.60000000e-04  2.50000000e-04  3.60000000e-04  4.90000000e-04
  6.40000000e-04  8.10000000e-04  1.00000000e-03  1.21000000e-03
  1.44000000e-03  1.69000000e-03  1.96000000e-03  2.25000000e-03
  2.56000000e-03  2.89000000e-03  3.24000000e-03  3.61000000e-03
  4.00000000e-03  4.41000000e-03  4.84000000e-03  5.29000000e-03
  5.76000000e-03  6.25000000e-03  6.76000000e-03  7.29000000e-03
  7.84000000e-03  8.41000000e-03  9.00000000e-03  9.61000000e-03
  1.02400000e-02  1.08900000e-02  1.15600000e-02  1.22500000e-02
  1.29600000e-02  1.36900000e-02  1.44400000e-02  1.52100000e-02
  1.60000000e-02  1.68100000e-02  1.76400000e-02  1.84900000e-02
  1.93600000e-02  2.02500000e-02  2.11600000e-02  2.20900000e-02
  2.30400000e-02  2.40100000e-02  2.50000000e-02  2.60100000e-02
  2.70400000e-02  2.80900000e-02  2.91600000e-02  3.02500000e-02
  3.13600000e-02  3.24900000e-02  3.36400000e-02  3.48100000e-02
  3.60000000e-02  3.72100000e-02  3.84400000e-02  3.96900000e-02
  4.09600000e-02  4.22500000e-02  4.35600000e-02  4.48900000e-02
  4.62400000e-02  4.76100000e-02  4.90000000e-02  5.04100000e-02
  5.18400000e-02  5.32900000e-02  5.47600000e-02  5.62500000e-02
  5.77600000e-02  5.92900000e-02  6.08400000e-02  6.24100000e-02
  6.40000000e-02  6.56100000e-02  6.72400000e-02  6.88900000e-02
  7.05600000e-02  7.22500000e-02  7.39600000e-02  7.56900000e-02
  7.74400000e-02  7.92100000e-02  8.10000000e-02  8.28100000e-02
  8.46400000e-02  8.64900000e-02  8.83600000e-02  9.02500000e-02
  9.21600000e-02  9.40900000e-02  9.60397896e-02  9.92031162e-02
 -4.90488926e-01]

Result from "print(rod.velocity_collection[0,:])" :
[ 0.00000000e+000  0.00000000e+000  0.00000000e+000  0.00000000e+000
  0.00000000e+000  0.00000000e+000  0.00000000e+000  0.00000000e+000
  0.00000000e+000  0.00000000e+000  0.00000000e+000  0.00000000e+000
  0.00000000e+000  0.00000000e+000  0.00000000e+000  0.00000000e+000
  0.00000000e+000  0.00000000e+000  0.00000000e+000  0.00000000e+000
  0.00000000e+000  0.00000000e+000  0.00000000e+000  0.00000000e+000
  0.00000000e+000  0.00000000e+000  0.00000000e+000  0.00000000e+000
 -4.05163029e-317  2.35993845e-312 -1.35594067e-307  7.68381801e-303
 -4.29364406e-298  2.36538276e-293 -1.28444227e-288  6.87344240e-284
 -3.62398684e-279  1.88215235e-274 -9.62673614e-270  4.84793391e-265
 -2.40315398e-260  1.17231509e-255 -5.62641767e-251  2.65599781e-246
 -1.23285198e-241  5.62544491e-237 -2.52253046e-232  1.11126054e-227
 -4.80791535e-223  2.04228156e-218 -8.51419087e-214  3.48245509e-209
 -1.39695194e-204  5.49369171e-200 -2.11719240e-195  7.99261337e-191
 -2.95434589e-186  1.06876320e-181 -3.78219757e-177  1.30868553e-172
 -4.42517822e-168  1.46149105e-163 -4.71178896e-159  1.48197925e-154
 -4.54458778e-150  1.35787114e-145 -3.95034482e-141  1.11817162e-136
 -3.07713262e-132  8.22615315e-128 -2.13446211e-123  5.37065783e-119
 -1.30917101e-114  3.08852990e-110 -7.04405030e-106  1.55133868e-101
 -3.29510908e-097  6.74127363e-093 -1.32652337e-088  2.50691511e-084
 -4.54280247e-080  7.88007109e-076 -1.30609515e-071  2.06456485e-067
 -3.10612618e-063  4.43851817e-059 -6.01107295e-055  7.69877028e-051
 -9.30520518e-047  1.05924907e-042 -1.13359001e-038  1.13873112e-034
 -1.07222265e-030  9.44860328e-027 -7.77317739e-023  5.94246588e-019
 -4.18562158e-015  2.67733251e-011 -1.52084988e-007  7.40153882e-004
 -1.50973127e-001]

[armantekinalp]

Hi @emersonkt,

I realized there was a bug in the public version of PyElastica. Latest PR #14 should fix the issue.

Can you pull from master and try again?

Let me know if you have any problems.

[emersonkt]

Hi @armantekinalp,

Thank you very much for the support, however I am still having the same problem. My project is being done via Google Colab, so I guess the library is pulled on each build, right ?

[armantekinalp]

Hi @emersonkt ,

I am not familiar with how Google Colab pulls the library. If it is doing pip install pyelastica then I don't think you will see the updates, because we haven't done a new package release for pip.

Can you check if you can clone pyelastica from master to your Google Colab directory?

[emersonkt]

I get it. Your comment led me to this website https://medium.com/@ashwindesilva/how-to-use-google-colaboratory-to-clone-a-github-repository-e07cf8d3d22b which explains how to use git with google colab (if anyone else is interested). Now it is perfect !

Thank you very much for your excellent work!

[MehtaMeet54]

Hi @armantekinalp ,
I tried to create a pre-curved rod but am facing issues with the director array. It is coming up as a zero 3rd order vector instead of what I am providing as input for the directors while creating the rod model. The position array of the rod is as desired, and it is just the director array that is inconsistent.
How can I fix this?

[armantekinalp]

Hi @MehtaMeet54 ,

Can you please provide the version of PyElastica you are using and also how you are initializing the rod?

[MehtaMeet54]

VERSION = "0.1.0.post2"
Code:

n_elem = 100
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])
base_length = 1.0
base_radius = 0.05
base_area = np.pi * base_radius ** 2
density = 8000
nu = 0.2
E = 1e11
poisson_ratio=0.3
position = np.zeros((n_elem+1,3))
for i in range(0, n_elem+1):
  position[i][0] = (i/n_elem)**2/base_length #y
  position[i][1] = 0 #z
  position[i][2] = (base_length/n_elem)*i #x
position = np.transpose(position)
tangents = np.zeros((3,n_elem))

for i in range(0, 3):
  tangents[i] = position[i,1:] - position[i,:-1]
tangents /= _batch_norm(tangents)

d3 = tangents
d2 = np.zeros((n_elem,3))
d1 = np.zeros((n_elem,3))

for i in range(0,n_elem):
  d2[i][0] = 0
  d2[i][1] = 1
  d2[i][2] = 0
  d1[i] = np.cross(d2[i],np.transpose(tangents)[i])

d2 = np.transpose(d2)
d1 = np.transpose(d1)

directors = np.zeros((3,3,n_elem))
directors[0] = d1
directors[1] = d2
directors[2] = d3

shearable_rod = CosseratRod.straight_rod(n_elem,start,direction,normal,base_length,base_radius,density,nu,E,poisson_ratio,position=position,directors=directors)
shearable_rod.rest_kappa[:] = shearable_rod.kappa[:]
shearable_rod.rest_sigma[:] = shearable_rod.sigma[:]
dynamic_update_sim.append(shearable_rod)
dynamic_update_sim.constrain(shearable_rod).using(OneEndFixedRod, constrained_position_idx=(0,), constrained_director_idx=(0,))
force = 100
dynamic_update_sim.add_forcing_to(shearable_rod).using(UniformForces, force, np.array([-1.0,0.0,0.0]))

dynamic_update_sim.finalize()

[skim0119]

@MehtaMeet54 Hi

Thank you for providing the code. I checked the director is not initialized properly in the version 0.1.0.post2.

There was a bug in initializing rod with custom directors. The bug is fixed in 0.1.0.post3 with 8af5172. Please test your code with updated pyelastica pip install --upgrade pyelastica. I tested on both version 0.2.1 and 0.1.0.post3.

Code I tested

Click to expand

import numpy as np
from elastica import *

from elastica._linalg import _batch_norm

n_elem = 100
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])
base_length = 1.0
base_radius = 0.05
base_area = np.pi * base_radius ** 2
density = 8000
nu = 0.2
E = 1e11
poisson_ratio=0.3
position = np.zeros((n_elem+1,3))
for i in range(0, n_elem+1):
  position[i][0] = (i/n_elem)**2/base_length #y
  position[i][1] = 0 #z
  position[i][2] = (base_length/n_elem)*i #x
position = np.transpose(position)
tangents = np.zeros((3,n_elem))

for i in range(0, 3):
  tangents[i] = position[i,1:] - position[i,:-1]
tangents /= _batch_norm(tangents)

d3 = tangents
d2 = np.zeros((n_elem,3))
d1 = np.zeros((n_elem,3))

for i in range(0,n_elem):
  d2[i][0] = 0
  d2[i][1] = 1
  d2[i][2] = 0
  d1[i] = np.cross(d2[i],np.transpose(tangents)[i])

d2 = np.transpose(d2)
d1 = np.transpose(d1)

directors = np.zeros((3,3,n_elem))
directors[0] = d1
directors[1] = d2
directors[2] = d3

shearable_rod = CosseratRod.straight_rod(n_elem,start,direction,normal,base_length,base_radius,density,nu,E,poisson_ratio,position=position,directors=directors)
shearable_rod.rest_kappa[:] = shearable_rod.kappa[:]
shearable_rod.rest_sigma[:] = shearable_rod.sigma[:]

print(f"{np.allclose(shearable_rod.position_collection, position)=}")
print(f"{np.allclose(shearable_rod.director_collection, directors)=}")

Result (python 3.8.0, pyelastica 0.2.1)

Click to expand

np.allclose(shearable_rod.position_collection, position)=True
np.allclose(shearable_rod.director_collection, directors)=True