Hello PyElastica, I am a student and I am trying to apply the librar

Hi <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="

Here below the forcing class: <div class="snippet-clipboard-content notranslate po

Hi <a class="user-mention notranslate" data-hovercard-type="user" data-hovercard-url="

Create a pre-curved rod about pyelastica HOT 12 CLOSED

gazzolalab commented on September 24, 2024

Create a pre-curved rod

from pyelastica.

Comments (12)

armantekinalp commented on September 24, 2024

Hi @emersonkt

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

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 .
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)

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.

from pyelastica.

emersonkt commented on September 24, 2024

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)

from pyelastica.

armantekinalp commented on September 24, 2024

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?

from pyelastica.

emersonkt commented on September 24, 2024

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]

from pyelastica.

armantekinalp commented on September 24, 2024

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.

from pyelastica.

emersonkt commented on September 24, 2024

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 ?

from pyelastica.

armantekinalp commented on September 24, 2024

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?

from pyelastica.

emersonkt commented on September 24, 2024

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!

from pyelastica.

MehtaMeet54 commented on September 24, 2024

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?

from pyelastica.

armantekinalp commented on September 24, 2024

Hi @MehtaMeet54 ,

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

from pyelastica.

MehtaMeet54 commented on September 24, 2024

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()

from pyelastica.

skim0119 commented on September 24, 2024

@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

from pyelastica.

Create a pre-curved rod about pyelastica HOT 12 CLOSED

Comments (12)

Code I tested

Result (python 3.8.0, pyelastica 0.2.1)

Related Issues (20)

Recommend Projects

React

Vue.js

Typescript

TensorFlow

Django

Laravel

D3

Recommend Topics

javascript

web

server

Machine learning

Visualization

Game

Recommend Org

Facebook

Microsoft

Google

Alibaba

D3

Tencent