-
Notifications
You must be signed in to change notification settings - Fork 20
Illustrative ex for viz #11
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: master
Are you sure you want to change the base?
Changes from 1 commit
06398c4
11a0059
d5f2d44
a61e602
e69d5cf
89db91f
b10d0e5
b1a923c
aa6a338
17f8f0b
5a79f75
134bf88
f769203
31e7b65
7c348e9
9c10bf6
bcbca28
f0aa71b
3c82531
8441d29
ae2a57b
e60fe75
f950c97
cb5ba46
8314218
d11bd84
5ef35c4
733f2ca
85b715e
9fcf6bb
4ea2488
0fec414
18ad2b2
cde87c8
53bddd4
896a29f
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
pydy_viz
- Loading branch information
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,221 @@ | ||
#This program represents a hypothetical situation for a complete workflow | ||
#for simulating a three link pendulum with links as rigid bodies | ||
from sympy import symbols,sympify | ||
from sympy.physics.mechanics import * | ||
|
||
#Number of links = 3 | ||
N_links = 3 | ||
|
||
#Number of masses = 3 | ||
N_bobs = 3 | ||
|
||
#Defining Dynamic Symbols ................ | ||
|
||
#Generalized coordinates(angular) ... | ||
|
||
alpha = dynamicsymbols('alpha1 alpha2 alpha3') | ||
beta = dynamicsymbols('beta1 beta2 beta3') | ||
|
||
#Generalized speeds(angular) ... | ||
alphad = dynamicsymbols('alpha1 alpha2 alpha3',1) | ||
betad = dynamicsymbols('beta1 beta2 beta3',1) | ||
|
||
#Mass of each bob: | ||
m = symbols('m:'+str(N_bobs)) | ||
|
||
|
||
#Length mass and radii of each link(assuming as rods) .. | ||
l = symbols('l:' + str(N_links)) | ||
M = symbols('M:' + str(N_links)) | ||
radii = symbols('radii:' + str(N_links)) | ||
|
||
#For storing Inertia for each link : | ||
Ixx = symbols('Ixx:'+str(N_links)) | ||
Iyy = symbols('Iyy:'+str(N_links)) | ||
|
||
#gravity and time .... | ||
g, t = symbols('g t') | ||
|
||
|
||
#Now defining an Inertial ReferenceFrame First .... | ||
|
||
I = ReferenceFrame('I') | ||
|
||
#And some other frames ... | ||
|
||
A = ReferenceFrame('A') | ||
A.orient(I,'Body',[alpha[0],beta[0],0],'ZXY') | ||
B = ReferenceFrame('B') | ||
B.orient(I,'Body',[alpha[1],beta[1],0],'ZXY') | ||
C = ReferenceFrame('C') | ||
C.orient(I,'Body',[alpha[2],beta[2],0],'ZXY') | ||
|
||
|
||
#Setting angular velocities of new frames ... | ||
A.set_ang_vel(I, alphad[0] * I.z + betad[0] * I.x) | ||
B.set_ang_vel(I, alphad[1] * I.z + betad[1] * I.x) | ||
C.set_ang_vel(I, alphad[2] * I.z + betad[2] * I.x) | ||
|
||
|
||
|
||
# An Origin point, with velocity = 0 | ||
O = Point('O') | ||
O.set_vel(I,0) | ||
|
||
#Three more points, for masses .. | ||
P1 = O.locatenew('P1', l[0] * A.y) | ||
P2 = O.locatenew('P2', l[1] * B.y) | ||
P3 = O.locatenew('P3', l[2] * C.y) | ||
|
||
#Setting velocities of points with v2pt theory ... | ||
P1.v2pt_theory(O, I, A) | ||
P2.v2pt_theory(P1, I, B) | ||
P3.v2pt_theory(P2, I, C) | ||
points = [P1,P2,P3] | ||
|
||
Pa1 = Particle('Pa1', points[0], m[0]) | ||
Pa2 = Particle('Pa2', points[1], m[1]) | ||
Pa3 = Particle('Pa3', points[2], m[2]) | ||
particles = [Pa1,Pa2,Pa3] | ||
|
||
|
||
|
||
#defining points for links(RigidBodies) | ||
#Assuming CoM as l/2 ... | ||
P_link1 = O.locatenew('P_link1', l[0]/2 * A.y) | ||
P_link2 = O.locatenew('P_link1', l[1]/2 * B.y) | ||
P_link3 = O.locatenew('P_link1', l[2]/2 * C.y) | ||
|
||
#setting velocities of these points with v2pt theory ... | ||
P_link1.v2pt_theory(O, I, A) | ||
P_link2.v2pt_theory(P_link1, I, B) | ||
P_link3.v2pt_theory(P_link2, I, C) | ||
|
||
points_rigid_body = [P_link1,P_link2,P_link3] | ||
|
||
|
||
#defining inertia tensors for links | ||
|
||
inertia_link1 = inertia(A,Ixx[0],Iyy[0],0) | ||
inertia_link2 = inertia(B,Ixx[1],Iyy[1],0) | ||
inertia_link3 = inertia(C,Ixx[2],Iyy[2],0) | ||
|
||
#Defining links as Rigid bodies ... | ||
|
||
link1 = RigidBody('link1', P_link1, A, M[0], (inertia_link1, O)) | ||
link2 = RigidBody('link2', P_link2, B, M[1], (inertia_link2, P1)) | ||
link3 = RigidBody('link3', P_link3, C, M[2], (inertia_link3, P2)) | ||
links = [link1,link2,link3] | ||
|
||
|
||
#Defining a basic shape for links .. | ||
rod1 = Cylinder(length=l[0],radii=radii[0]) | ||
rod2 = Cylinder(length=l[1],radii=radii[1]) | ||
rod3 = Cylinder(length=l[2],radii=radii[2]) | ||
|
||
link1.shape(rod1) | ||
link2.shape(rod2) | ||
link3.shape(rod3) | ||
|
||
#Applying forces on all particles , and adding all forces in a list.. | ||
forces = [] | ||
for particle in particles: | ||
|
||
mass = particle.get_mass() | ||
point = particle.get_point() | ||
forces.append((point, -mass * g * I.y) ) | ||
|
||
#Applying forces on rigidbodies .. | ||
for link in links: | ||
mass = link.get_mass() | ||
point = link.get_masscenter() | ||
forces.append((point, -mass * g * I.y) ) | ||
|
||
kinetic_differentials = [] | ||
for i in range(0,N_bobs): | ||
kinetic_differentials.append(alphad[i] - alpha[i]) | ||
kinetic_differentials.append(betad[i] - beta[i]) | ||
|
||
#Adding particles and links in the same system ... | ||
total_system = [] | ||
for particle in particles: | ||
total_system.append(particle) | ||
|
||
for link in links: | ||
total_system.append(link) | ||
|
||
q = [] | ||
for angle in alpha: | ||
q.append(angle) | ||
for angle in beta: | ||
q.append(angle) | ||
print q | ||
u = [] | ||
|
||
for vel in alphad: | ||
u.append(vel) | ||
for vel in betad: | ||
u.append(vel) | ||
|
||
print u | ||
kane = KanesMethod(I, q_ind=q, u_ind=u, kd_eqs=kinetic_differentials) | ||
fr, frstar = kane.kanes_equations(forces, total_system) | ||
|
||
print fr | ||
|
||
#Now we have symbolic equations of motion. .. | ||
# we integrate them numerically. .. | ||
|
||
params = [g ,l1,l2,l3,m1,m2,m3,M1,M2,M3] | ||
|
||
param_vals = [9.8 ,1.0,1.0,1.0,2,2,2,5,5,5] | ||
|
||
right_hand_side = code_generator(kane,params) | ||
|
||
#setting initial conditions .. | ||
init_conditions = [radians(45),radians(45),radians(30),\ | ||
radians(30),radians(15),radians(15),\ | ||
0, 0, 0,\ | ||
0, 0, 0] | ||
t = [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0] | ||
|
||
numerical_vals = odeint(right_hand_side,init_conditions,t) | ||
|
||
#Now for each t, we have numerical vals of coordinates .. | ||
#Now we set up a visualization frame, | ||
|
||
frame1 = VisualizationFrame('frame1',I,O) | ||
|
||
frame1.add_rigidbodies(links) | ||
|
||
frame1.add_particles(particles) | ||
|
||
param_vals_for_viz = {'g':9.8 ,'l1':1.0,'l2':1.0,'l3':1.0,'m1':2,'m2':2,'m3':2,'M1':5,'M1':5,'M1':5] | ||
|
||
json = frame1.generate_json(initial_conditions,q) | ||
#Here we can replace initial_conditions with the conditions at any | ||
#specific time interval .... | ||
|
||
##This line does following things ... | ||
##calls RigidBody.transform_matrix() on all rigidbodies, | ||
#so that we have info of rigidbodies(CoM,rotation,translation) | ||
# w.r.t VisualizationFrame('I' in this case) | ||
##Even if they are defined in any other frame .... | ||
##calls point.set_pos() for all particles with arg | ||
##as Point in VisualizationFrame(O here) .. | ||
##With these and init_conditions, we have a starting point of visualization .. | ||
|
||
scene = Scene() | ||
Scene.view(json) # Just visualize/view the system at initial_conditions, | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. lower-case scene There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. yes it should be scene ... my bad .... On 6/23/13, Christopher Dembia [email protected] wrote:
|
||
#w.r.t VIsualizationFrame(I in this case) | ||
# No simulation here. .. | ||
|
||
Scene.simulate(json,numerical_vals) # modify the input json, | ||
#Add the values at different times from numerical_vals, | ||
#To the json, which is then passed to | ||
#javascript | ||
# | ||
#(alpha1,alpha2,alpha3,beta1,beta2,beta3) at time=t | ||
|
||
|
||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,157 @@ | ||
#For this one we assume RigidBody links | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Why is this file here? I though you have this in another pull request. Do not include it in both PR's. |
||
from sympy import symbols,sympify | ||
from sympy.physics.mechanics import * | ||
|
||
#Number of links = 3 | ||
N_links = 3 | ||
|
||
#Number of masses = 3 | ||
N_bobs = 3 | ||
|
||
#Defining Dynamic Symbols ................ | ||
|
||
#Generalized coordinates(angular) ... | ||
|
||
alpha = dynamicsymbols('alpha1 alpha2 alpha3') | ||
beta = dynamicsymbols('beta1 beta2 beta3') | ||
|
||
#Generalized speeds(angular) ... | ||
alphad = dynamicsymbols('alpha1 alpha2 alpha3',1) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. These are not technically the generalized speeds. The generalized speeds should be defined as: omega = alpha' You need to define some generalized speeds, not just the derivatives of the coordinates. so: omega1, omega2, omega3 = dynamicsymbols('omega1 omega2 omega3') And similarly for beta. |
||
betad = dynamicsymbols('beta1 beta2 beta3',1) | ||
|
||
#Mass of each bob: | ||
m = symbols('m:'+str(N_bobs)) | ||
|
||
|
||
#Length and mass of each link .. | ||
l = symbols('l:' + str(N_links)) | ||
M = symbols('M:' + str(N_links)) | ||
#For storing Inertia for each link : | ||
Ixx = symbols('Ixx:'+str(N_links)) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. PEP8: spaces around operators |
||
Iyy = symbols('Iyy:'+str(N_links)) | ||
|
||
#gravity and time .... | ||
g, t = symbols('g t') | ||
|
||
|
||
#Now defining an Inertial ReferenceFrame First .... | ||
|
||
I = ReferenceFrame('I') | ||
|
||
#And some other frames ... | ||
|
||
A = ReferenceFrame('A') | ||
A.orient(I,'Body',[alpha[0],beta[0],0],'ZXY') | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. PEP8: space after comma |
||
B = ReferenceFrame('B') | ||
B.orient(I,'Body',[alpha[1],beta[1],0],'ZXY') | ||
C = ReferenceFrame('C') | ||
C.orient(I,'Body',[alpha[2],beta[2],0],'ZXY') | ||
|
||
|
||
#Setting angular velocities of new frames ... | ||
A.set_ang_vel(I, alphad[0] * I.z + betad[0] * I.x) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. these should be A.set_ang_vel(I, omega1 * I.z + gamma1 * I.x) |
||
B.set_ang_vel(I, alphad[1] * I.z + betad[1] * I.x) | ||
C.set_ang_vel(I, alphad[2] * I.z + betad[2] * I.x) | ||
|
||
|
||
|
||
# An Origin point, with velocity = 0 | ||
O = Point('O') | ||
O.set_vel(I,0) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. PEP8, space after comma |
||
|
||
#Three more points, for masses .. | ||
P1 = O.locatenew('P1', l[0] * A.y) | ||
P2 = O.locatenew('P2', l[1] * B.y) | ||
P3 = O.locatenew('P3', l[2] * C.y) | ||
|
||
#Setting velocities of points with v2pt theory ... | ||
P1.v2pt_theory(O, I, A) | ||
P2.v2pt_theory(P1, I, B) | ||
P3.v2pt_theory(P2, I, C) | ||
points = [P1,P2,P3] | ||
|
||
Pa1 = Particle('Pa1', points[0], m[0]) | ||
Pa2 = Particle('Pa2', points[1], m[1]) | ||
Pa3 = Particle('Pa3', points[2], m[2]) | ||
particles = [Pa1,Pa2,Pa3] | ||
|
||
|
||
|
||
#defining points for links(RigidBodies) | ||
#Assuming CoM as l/2 ... | ||
P_link1 = O.locatenew('P_link1', l[0]/2 * A.y) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. PEP8: Space around operators. |
||
P_link2 = O.locatenew('P_link1', l[1]/2 * B.y) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Shouldn't this be: P_link2 = P1.locatenew('P_link1', l[1] / 2 * B.y) |
||
P_link3 = O.locatenew('P_link1', l[2]/2 * C.y) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This should be: P_link3 = P2.locatenew('P_link1', l[2] / 2 * C.y) |
||
#setting velocities of these points with v2pt theory ... | ||
P_link1.v2pt_theory(O, I, A) | ||
P_link2.v2pt_theory(P_link1, I, B) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. P_link1 and P_link2 do not exist on the same rigid body so this is incorrect. It should be: P_link2.v2pt_theory(P1, I, B) There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. did you meant do not exist on same referenceframe? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. rigid body and reference frame are synonymous in this case. The point should be the one that is at the joints. |
||
P_link3.v2pt_theory(P_link2, I, C) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Same as above. |
||
|
||
points_rigid_body = [P_link1,P_link2,P_link3] | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. PEP8: Spaces after commas. |
||
|
||
|
||
#defining inertia tensors for links | ||
|
||
inertia_link1 = inertia(A,Ixx[0],Iyy[0],0) | ||
inertia_link2 = inertia(B,Ixx[1],Iyy[1],0) | ||
inertia_link3 = inertia(C,Ixx[2],Iyy[2],0) | ||
|
||
#Defining links as Rigid bodies ... | ||
|
||
link1 = RigidBody('link1', P_link1, A, M[0], (inertia_link1, O)) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Why do you define your inertia about the origin? Isn't the better way to define your inertia about the center of mass of each link. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I defined on the rotation axis of each rigidbody, There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The "clean" way to define your inertia tensor is about the body's reference frame and the center of mass. (at least in this case). Other times it may be useful to define inertia about other points, depending on the problem. |
||
link2 = RigidBody('link2', P_link2, B, M[1], (inertia_link2, P1)) | ||
link3 = RigidBody('link3', P_link3, C, M[2], (inertia_link3, P2)) | ||
links = [link1,link2,link3] | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. PEP8 |
||
|
||
|
||
#Applying forces on all particles , and adding all forces in a list.. | ||
forces = [] | ||
for particle in particles: | ||
|
||
mass = particle.get_mass() | ||
point = particle.get_point() | ||
forces.append((point, -mass * g * I.y) ) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. A one liner: forces = [(p.get_point(), -p.get_mass() * g * I.y) for p in particles] |
||
|
||
#Applying forces on rigidbodies .. | ||
for link in links: | ||
mass = link.get_mass() | ||
point = link.get_masscenter() | ||
forces.append((point, -mass * g * I.y) ) | ||
kinetic_differentials = [] | ||
for i in range(0,N_bobs): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
|
||
kinetic_differentials.append(alphad[i] - alpha[i]) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. this should be kinetic_differentials.append(omega1 - alphad[i]) omega = alpha' There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. And it really should be |
||
kinetic_differentials.append(betad[i] - beta[i]) | ||
|
||
#Adding particles and links in the same system ... | ||
total_system = [] | ||
for particle in particles: | ||
total_system.append(particle) | ||
|
||
for link in links: | ||
total_system.append(link) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This can be a one liner:
|
||
|
||
q = [] | ||
for angle in alpha: | ||
q.append(angle) | ||
for angle in beta: | ||
q.append(angle) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. q = alpha + beta There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. you may have to do: q = list(q) + list(beta) |
||
print q | ||
u = [] | ||
|
||
for vel in alphad: | ||
u.append(vel) | ||
for vel in betad: | ||
u.append(vel) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. u = omega + delta |
||
|
||
print u | ||
kane = KanesMethod(I, q_ind=q, u_ind=u, kd_eqs=kinetic_differentials) | ||
fr, frstar = kane.kanes_equations(forces, total_system) | ||
|
||
print fr | ||
|
||
|
||
|
||
|
||
|
||
|
||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -148,10 +148,3 @@ | |
fr, frstar = kane.kanes_equations(forces, total_system) | ||
|
||
print fr | ||
|
||
|
||
|
||
|
||
|
||
|
||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is a picky comment, but differentials are different from differential equations. Also are these kinematic or kinetic? The word kinetic is sometimes used to mean "dynamic" (that is, coming from F=ma).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
these are kinematic equations.
On 6/23/13, Christopher Dembia [email protected] wrote:
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
yes, the correct term is "kinematic differential equations"