Mercurial > lcfOS

comparison python/other/bouncing_cube_opencl.py @ 290:7b38782ed496

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
File moves
author Windel Bouwman
date Sun, 24 Nov 2013 11:24:15 +0100
parents python/bouncing_cube_opencl.py@5a965e9664f2
children
comparison
equal deleted inserted replaced
289:bd2593de3ff8 290:7b38782ed496
1 from PyQt4.QtGui import *
2 from PyQt4.QtCore import *
3 from PyQt4.QtOpenGL import QGLWidget
4 from OpenGL.GL import *
5 from OpenGL.GLU import gluPerspective
6 import sys
7 from random import random
8 from math import pi, cos, sin, fabs, sqrt
9 from numpy import mat, array, ones, zeros, eye
10 from numpy.linalg import norm
11 import numpy as np
12 import time
13 import scipy.integrate
14 import pyopencl
15
16 """
17 Test script that lets a dice bounce.
18 Converted from 20-sim equations into python code.
19
20 20-sim website:
21 http://www.20sim.com
22
23 """
24 def drawCube(w):
25 glBegin(GL_QUADS) # Start Drawing The Cube
26 glColor3f(0.0,1.0,0.0) # Set The Color To Blue
27 glVertex3f( w, w,-w) # Top Right Of The Quad (Top)
28 glVertex3f(-w, w,-w) # Top Left Of The Quad (Top)
29 glVertex3f(-w, w, w) # Bottom Left Of The Quad (Top)
30 glVertex3f( w, w, w) # Bottom Right Of The Quad (Top)
31
32 glColor3f(1.0,0.5,0.0) # Set The Color To Orange
33 glVertex3f( w,-w, w) # Top Right Of The Quad (Bottom)
34 glVertex3f(-w,-w, w) # Top Left Of The Quad (Bottom)
35 glVertex3f(-w,-w,-w) # Bottom Left Of The Quad (Bottom)
36 glVertex3f( w,-w,-w) # Bottom Right Of The Quad (Bottom)
37
38 glColor3f(1.0,0.0,0.0) # Set The Color To Red
39 glVertex3f( w, w, w) # Top Right Of The Quad (Front)
40 glVertex3f(-w, w, w) # Top Left Of The Quad (Front)
41 glVertex3f(-w,-w, w) # Bottom Left Of The Quad (Front)
42 glVertex3f( w,-w, w) # Bottom Right Of The Quad (Front)
43
44 glColor3f(1.0,1.0,0.0) # Set The Color To Yellow
45 glVertex3f( w,-w,-w) # Bottom Left Of The Quad (Back)
46 glVertex3f(-w,-w,-w) # Bottom Right Of The Quad (Back)
47 glVertex3f(-w, w,-w) # Top Right Of The Quad (Back)
48 glVertex3f( w, w,-w) # Top Left Of The Quad (Back)
49
50 glColor3f(0.0,0.0,1.0) # Set The Color To Blue
51 glVertex3f(-w, w, w) # Top Right Of The Quad (Left)
52 glVertex3f(-w, w,-w) # Top Left Of The Quad (Left)
53 glVertex3f(-w,-w,-w) # Bottom Left Of The Quad (Left)
54 glVertex3f(-w,-w, w) # Bottom Right Of The Quad (Left)
55
56 glColor3f(1.0,0.0,1.0) # Set The Color To Violet
57 glVertex3f( w, w,-w) # Top Right Of The Quad (Right)
58 glVertex3f( w, w, w) # Top Left Of The Quad (Right)
59 glVertex3f( w,-w, w) # Bottom Left Of The Quad (Right)
60 glVertex3f( w,-w,-w) # Bottom Right Of The Quad (Right)
61 glEnd() # Done Drawing The Quad
62
63 def drawFloor(w, h):
64 glBegin(GL_QUADS) # Start Drawing The Cube
65
66 glColor3f(1.0,0.5,0.0) # Set The Color To Orange
67 glVertex3f( w,-w,h)# Top Right Of The Quad (Bottom)
68 glVertex3f(-w,-w,h)# Top Left Of The Quad (Bottom)
69 glVertex3f(-w,w,h)# Bottom Left Of The Quad (Bottom)
70 glVertex3f( w,w,h)# Bottom Right Of The Quad (Bottom)
71 glEnd() # Done Drawing The Quad
72
73 def drawAxis():
74 glLineWidth(0.5)
75 glBegin(GL_LINES)
76 glColor3f(1.0, 0.0, 0.0)
77 glVertex3f(0,0,0)
78 glVertex3f(1,0,0)
79 glColor3f(0.0, 1.0, 0.0)
80 glVertex3f(0,0,0)
81 glVertex3f(0,1,0)
82 glColor3f(0.0, 0.0, 1.0)
83 glVertex3f(0,0,0)
84 glVertex3f(0,0,1)
85 glEnd()
86
87
88 def cross(A, B):
89 a = A.A1
90 b = B.A1
91 return mat(np.cross(a, b)).T
92
93 def skew(X):
94 Y = mat(zeros( (3, 3) ))
95 a,b,c = X.A1
96 Y[0,1] = -c
97 Y[0,2] = b
98 Y[1,0] = c
99 Y[1,2] = -a
100 Y[2,0] = -b
101 Y[2,1] = a
102 return Y
103
104 def adjoint(T):
105 W = T[0:3, 0]
106 V = T[3:6, 0]
107 a = mat(zeros( (6,6) ) )
108 a[0:3, 0:3] = skew(W)
109 a[3:6, 0:3] = skew(V)
110 a[3:6, 3:6] = skew(W)
111 return a
112
113 def Adjoint(H):
114 R = H[0:3, 0:3]
115 P = H[0:3, 3]
116 a = mat(zeros( (6,6) ) )
117 a[0:3, 0:3] = R
118 a[3:6, 3:6] = R
119 a[3:6, 0:3] = skew(P) * R
120 return a
121
122 def quatToR(q):
123 x, y, z, w = q.A1
124 r = mat(eye(3))
125 r[0,0] = 1 - (2*y**2+2*z**2)
126 r[0,1] = 2*x*y+2*z*w
127 r[0,2] = 2*x*z - 2*y*w
128 r[1,0] = 2*x*y-2*z*w
129 r[1,1] = 1 - (2*x**2 + 2*z**2)
130 r[1,2] = 2*y*z + 2*x*w
131 r[2,0] = 2*x*z+2*y*w
132 r[2,1] = 2*y*z - 2*x*w
133 r[2,2] = 1 - (2*x**2+2*y**2)
134 return r
135
136 def rotateAbout(axis, angle):
137 ax, ay, az = (axis/norm(axis)).A1
138 qx = ax*sin(angle/2.0)
139 qy = ay*sin(angle/2.0)
140 qz = az*sin(angle/2.0)
141 qw = cos(angle/2.0)
142 q = mat(array([qx,qy,qz,qw])).T
143 return q
144
145 def normalizeQuaternion(quat):
146 x,y,z,w = quat.A1
147 magnitude = sqrt(x*x + y*y + z*z + w*w)
148 x = x / magnitude
149 y = y / magnitude
150 z = z / magnitude
151 w = w / magnitude
152 quat[0, 0] = x
153 quat[1, 0] = y
154 quat[2, 0] = z
155 quat[3, 0] = w
156 return quat
157
158 def VTo4x4(V):
159 v1, v2, v3 = V.A1
160 return mat(array( \
161 [[0.0, -v3, v2, -v1], \
162 [ v3, 0.0, -v1, -v2], \
163 [-v2, v1, 0.0, -v3], \
164 [v1, v2, v3, 0.0] ]))
165
166 def homogeneous(R,p):
167 H = mat(eye(4))
168 H[0:3, 0:3] = R
169 H[0:3, 3] = p
170 return H
171
172 def rateOfChange(states, thetime, parameters):
173 quat = states[0:4, 0] # Orientation (4)
174 pos = states[4:7, 0] # Position (3)
175 P = states[7:13, 0] # Momentum (6)
176 massI, gravity = parameters
177 # Rigid body parts:
178 # Forward Kinematic chain:
179 H = homogeneous(quatToR(quat), pos) # Forward kinematics
180
181 AdjX2 = mat(eye(6)) # The connectionpoint in the real world
182 adjAdjX2_1 = adjoint(AdjX2[0:6,0])
183 adjAdjX2_2 = adjoint(AdjX2[0:6,1])
184 adjAdjX2_3 = adjoint(AdjX2[0:6,2])
185 adjAdjX2_4 = adjoint(AdjX2[0:6,3])
186 adjAdjX2_5 = adjoint(AdjX2[0:6,4])
187 adjAdjX2_6 = adjoint(AdjX2[0:6,5])
188 AdjInv2 = Adjoint(H.I)
189 M2 = AdjInv2.T * (massI * AdjInv2) # Transfor mass to base frame
190 MassMatrix = M2
191
192 wrenchGrav2 = mat( zeros((1,6)) )
193 wrenchGrav2[0, 0:3] = -cross(gravity, pos).T
194 wrenchGrav2[0, 3:6] = gravity.T
195
196 Bk = mat( zeros( (6,6) ))
197 Bk[0:3, 0:3] = skew(P[0:3, 0])
198 Bk[3:6, 0:3] = skew(P[3:6, 0])
199 Bk[0:3, 3:6] = skew(P[3:6, 0])
200
201 # TODO: do this a cholesky:
202 v = np.linalg.solve(MassMatrix, P) # Matrix inverse like thingy !
203
204 T2_00 = v # Calculate the relative twist!
205 TM2 = T2_00.T * M2
206 twistExternal = T2_00 # Twist van het blokje
207 TMSum2 = TM2
208
209 PDotBodies = mat( zeros( (6,1)) )
210 PDotBodies[0,0] = TMSum2 * (adjAdjX2_1 * T2_00)
211 PDotBodies[1,0] = TMSum2 * (adjAdjX2_2 * T2_00)
212 PDotBodies[2,0] = TMSum2 * (adjAdjX2_3 * T2_00)
213 PDotBodies[3,0] = TMSum2 * (adjAdjX2_4 * T2_00)
214 PDotBodies[4,0] = TMSum2 * (adjAdjX2_5 * T2_00)
215 PDotBodies[5,0] = TMSum2 * (adjAdjX2_6 * T2_00)
216
217 PDot = -PDotBodies - Bk * v
218 PDot += wrenchGrav2.T
219
220 ##### Contact wrench part:
221 HB_W = H # Is H-matrix van het blokje
222 WrenchB = mat(zeros( (1,6) ))
223 for px in [-0.5, 0.5]:
224 for py in [-0.5, 0.5]:
225 for pz in [-0.5, 0.5]:
226 HB1_B = homogeneous(mat(eye(3)), mat([px,py,pz]).T)
227 HB1_W = HB_W * HB1_B
228 HW1_W = homogeneous(mat(eye(3)), HB1_W[0:3,3])
229 HW_W1 = HW1_W.I
230 HB_W1 = HW_W1 * HB_W
231
232 AdjHB_W1 = Adjoint(HB_W1)
233 TB_W1_W1 = AdjHB_W1 * twistExternal
234 z = HB1_W[2, 3]
235 vx, vy, vz = TB_W1_W1[3:6, 0].A1
236 if z < 0:
237 # Contact forces:
238 Fx = -50.0*vx
239 Fy = -50.0*vy
240 Fz = -z*50000.0
241 else:
242 Fx = 0.0
243 Fy = 0.0
244 Fz = 0.0
245 # TODO: reflect impulse
246 WrenchW1 = mat([0,0,0,0,0,Fz])
247 # Transform it back:
248 WrenchB += (AdjHB_W1.T * WrenchW1.T).T
249 ##### End of contact wrench
250
251 PDot += (WrenchB * AdjInv2).T
252
253 # Position and orientation rates:
254 QOmega = VTo4x4(v[0:3, 0])
255 quatDot = 0.5 * QOmega * quat
256 vel = v[3:6, 0]
257 posDot = skew(v[0:3]) * pos + vel
258 # The rate vector:
259 rates = mat(zeros( (13,1) ))
260 rates[0:4, 0] = quatDot
261 rates[4:7, 0] = posDot
262 rates[7:13, 0] = PDot
263 return rates
264
265 def fWrapper(y, t, parameters):
266 y = mat(y).T
267 dy = rateOfChange(y, t, parameters)
268 return dy.T.A1
269
270 def SCIPY(endtime, dt, state, parameters):
271 times = np.arange(0.0, endtime, dt)
272 y0 = state.T.A1
273 res = scipy.integrate.odeint(fWrapper, y0, times, args=(parameters,))
274 states = []
275 res = res.T
276 r,c = res.shape
277 for ci in range(0,c):
278 states.append(mat(res[:,ci]).T)
279 return states
280
281 def RK4(endtime, dt, state, parameters):
282 t = 0.0
283 states = []
284 while t < endtime:
285 newstate = mat (zeros( state.shape )) # Create a new object
286
287 #### Runge Kutta integration:
288 k1 = rateOfChange(state, t, parameters)
289 k2 = rateOfChange(state + 0.5*dt*k1, t, parameters)
290 k3 = rateOfChange(state + 0.5*dt*k1, t, parameters)
291 k4 = rateOfChange(state + dt*k3, t, parameters)
292 newstate = state + (dt/6.0)*(k1+2*k2+2*k3+k4)
293
294 # Normalize quat:
295 newstate[0:4, 0] = normalizeQuaternion(newstate[0:4, 0])
296 states.append(newstate)
297
298 state = newstate
299
300 t += dt
301 print state[6,0], t, ' ', (t/endtime)*100.0, '%'
302 return states
303
304
305 def simulate(endtime, dt):
306 PInitial = mat( zeros((6,1)) )
307 posInitial = mat(array([-1.2, -1.3, 2.8])).T
308 quatInitial = rotateAbout(mat(array([0.2, 1.0, 0.4])).T, 0.5)
309 # Parameters:
310 gravity = mat( array([0,0,-9.81]) ).T
311 massI = mat(eye(6)) * 1.01 # Mass matrix
312 parameters = (massI, gravity)
313
314 # The state vector!
315 state = mat(zeros( (13,1) ))
316 state[0:4, 0] = quatInitial
317 state[4:7, 0] = posInitial
318 state[7:13, 0] = PInitial
319
320 return SCIPY(endtime, dt, state, parameters)
321
322 class W(QGLWidget):
323 time = 0.0
324 index = 0
325 def __init__(self, states, dt, parent=None):
326 super(W, self).__init__(parent)
327 self.firstRun = True
328 self.savePNGS = False
329 self.dt = dt
330 self.resize(500,500)
331 self.states = states
332 self.UP()
333 t = QTimer(self)
334 t.timeout.connect(self.UP)
335 t.start(self.dt*1000.0)
336
337 def UP(self):
338 self.time += self.dt
339 if self.states:
340 state = self.states[self.index]
341 self.Dicequat = state[0:4, 0]
342 self.Dicepos = state[4:7, 0]
343 self.index += 1
344 if self.index == len(self.states):
345 self.index = 0
346 self.firstRun = False
347 # Paint:
348 self.update()
349 if self.firstRun:
350 # Create png images for the movie:
351 if self.savePNGS:
352 pm = self.renderPixmap()
353 pm.save('image'+str(self.index)+'.png')
354
355 def initializeGL(self):
356 glClearColor(0.0, 0.5, 0.0, 1.0)
357 glEnable(GL_DEPTH_TEST)
358 glDepthFunc(GL_LESS)
359 glShadeModel(GL_SMOOTH)
360 def resizeGL(self,w,h):
361 glViewport(0, 0, w, h)
362 glMatrixMode(GL_PROJECTION)
363 glLoadIdentity()
364 gluPerspective(45.0, float(w)/float(h), 0.1, 100.0)
365 glMatrixMode(GL_MODELVIEW)
366 def paintGL(self):
367 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) # Clear buffers
368 glLoadIdentity() # Reset The View
369
370 glLoadIdentity()
371 glTranslatef(0.0,-2.0,-10.0) # Move Left And Into The Screen
372 glRotatef(-90.0, 1.0, 0.0, 0.0)
373 drawFloor(2.0, 0.0)
374 drawAxis()
375 self.renderText(1.0, 0.0, 0.0, 'X')
376 self.renderText(0.0, 1.0, 0.0, 'Y')
377 self.renderText(0.0, 0.0, 1.0, 'Z')
378
379 self.renderText(0.0,0.0,1.2,str(self.time))
380
381 x,y,z = self.Dicepos.A1
382 R = quatToR(self.Dicequat)
383
384 glTranslatef(x, y, z)
385 # Trick to rotate the openGL matrix:
386 r = R.A1
387 rotR = (r[0], r[3], r[6], 0.0,
388 r[1], r[4], r[7], 0.0,
389 r[2], r[5], r[8], 0.0,
390 0.0, 0.0, 0.0, 1.0)
391 glMultMatrixd(rotR)
392
393 drawCube(0.6)
394
395 et = 20.0
396 dt = 0.04
397 print 'starting integration... endtime =', et, ' stepsize =', dt
398 t0 = time.time()
399 states = simulate(et, dt)
400 t1 = time.time()
401 print 'That was heavy, it took me ', t1-t0, ' seconds!'
402 app = QApplication(sys.argv)
403 w = W(states, dt)
404 w.show()
405 sys.exit(app.exec_())
406