Skip to content
Browse files

fixed cpu python sph

  • Loading branch information...
1 parent 667f38f commit 6124e8849bfedc0d8df886cd83ddb6440b441f1a @enjalot committed Sep 9, 2011
View
171 teach/pycpu/forces.py
@@ -0,0 +1,171 @@
+import pygame
+from pygame.locals import *
+
+from kernels import *
+from vector import Vec
+
+#from timing import print_timing
+from timing import Timing
+timings = Timing()
+
+#@print_timing
+@timings
+def density_update(sphp, particles):
+ #brute force
+ for pi in particles:
+ pi.dens = 0.
+ for pj in particles:
+ r = pi.pos - pj.pos
+ #print r
+ if mag(r) > pi.h: continue
+ #pi.dens += pj.mass*Wpoly6(pi.h, r)
+ pi.dens += pj.mass * sphp.kernels.poly6(r)
+
+#@print_timing
+@timings
+def force_update(sphp, particles):
+ #brute force
+ rho0 = sphp.rho0
+ K = sphp.K
+
+ for pi in particles:
+ pi.force = Vec([0.,0.])
+ pi.xsph = Vec([0.,0.])
+ di = pi.dens
+ #print "di", di
+ Pi = K*(di - rho0)
+ #print "Pi", Pi
+ #print "pi.h", pi.h
+ for pj in particles:
+ if pj == pi:
+ continue
+ r = pi.pos - pj.pos
+ #temporary optimization until we do efficient neighbor search
+ if mag(r) > pi.h: continue
+ #print "r", r
+
+ dj = pj.dens
+ Pj = K*(dj - rho0)
+
+ #print "dWspiky", dWspiky(pi.h, r)
+ #kern = .5 * (Pi + Pj) * dWspiky(pi.h, r)
+ kern = .5 * (Pi + Pj) * sphp.kernels.dspiky(r)
+ f = r*kern
+ #does not account for variable mass
+ f *= pi.mass / (di * dj)
+
+ #print "force", f
+ pi.force += f
+
+ #XSPH
+ #float4 xsph = (2.f * sphp->mass * Wijpol6 * (velj-veli)/(di.x+dj.x));
+ xsph = pj.veleval - pi.veleval
+ xsph *= 2. * pi.mass * sphp.kernels.poly6(r) / (di + dj)
+ pi.xsph += xsph
+
+
+ #print "pi.force", pi.force
+
+
+
+#@timings
+def calcRepulsionForce(normal, vel, sphp, distance):
+ repulsion_force = (sphp.boundary_stiffness * distance - sphp.boundary_dampening * np.dot(normal, vel))*normal;
+ return repulsion_force;
+
+#@timings
+def calcFrictionForce(v, f, normal, friction_kinetic, friction_static_limit):
+ pass
+
+#@print_timing
+@timings
+def collision_wall(sphp, domain, particles):
+
+ dmin = domain.bnd_min * sphp.sim_scale
+ dmax = domain.bnd_max * sphp.sim_scale
+ bnd_dist = sphp.boundary_distance
+ #float diff = params->boundary_distance - (p.z - gp->bnd_min.z);
+ #if (diff > params->EPSILON)
+ #r_f += calculateRepulsionForce(normal, v, params->boundary_stiffness, params->boundary_dampening, diff);
+ #f_f += calculateFrictionForce(v, f, normal, friction_kinetic, friction_static_limit);
+
+
+ for pi in particles:
+ f = Vec([0.,0.])
+ p = pi.pos
+ #X Min
+ diff = bnd_dist - (p.x - dmin.x)
+ if diff > 0:
+ normal = Vec([1.,0.])
+ f += calcRepulsionForce(normal, pi.vel, sphp, diff)
+ #calcFrictionForce(pi.v, pi.f, normal, friction_kinetic, friction_static_limit)
+ #X Max
+ diff = bnd_dist - (dmax.x - p.x)
+ if diff > 0:
+ normal = Vec([-1.,0.])
+ f += calcRepulsionForce(normal, pi.vel, sphp, diff)
+ #calcFrictionForce(pi.v, pi.f, normal, friction_kinetic, friction_static_limit)
+ #Y Min
+ diff = bnd_dist - (p.y - dmin.y)
+ if diff > 0:
+ normal = Vec([0.,1.])
+ f += calcRepulsionForce(normal, pi.vel, sphp, diff)
+ #calcFrictionForce(pi.v, pi.f, normal, friction_kinetic, friction_static_limit)
+ #Y Max
+ diff = bnd_dist - (dmax.y - p.y)
+ if diff > 0:
+ normal = Vec([0.,-1.])
+ f += calcRepulsionForce(normal, pi.vel, sphp, diff)
+ #calcFrictionForce(pi.v, pi.f, normal, friction_kinetic, friction_static_limit)
+
+ #print "collision force", f
+ pi.force += f
+
+
+
+#@print_timing
+@timings
+def euler_update(sphp, particles):
+ dt = .001
+
+ for pi in particles:
+ f = pi.force
+
+ f.y += -9.8
+
+ speed = mag(f);
+ #print "speed", speed
+ if speed > sphp.velocity_limit:
+ f *= sphp.velocity_limit/speed;
+
+
+ #print "force", f
+ pi.vel += f * dt
+ pi.pos += pi.vel * dt
+
+#@print_timing
+@timings
+def leapfrog_update(sphp, particles):
+ dt = .001
+
+ #print "LEAPFROG++++++++++++++++++++++"
+ for pi in particles:
+ f = pi.force
+ #print "f", f
+
+ f.y += -9.8
+
+ speed = mag(f);
+ #print "speed", speed
+ if speed > sphp.velocity_limit:
+ f *= sphp.velocity_limit/speed;
+
+ #print "force", f
+ vnext = pi.vel + f * dt
+ vnext += sphp.xsph_factor * pi.xsph
+ pi.pos += vnext * dt
+
+ veval = .5 * (pi.vel + vnext)
+ pi.vel = vnext
+ pi.veleval = veval
+
View
103 teach/pycpu/glutil.py
@@ -0,0 +1,103 @@
+#from OpenGL.GL import *
+import OpenGL.GL as gl
+from OpenGL.GLU import *
+from OpenGL.raw.GL.VERSION.GL_1_5 import glBufferData as rawGlBufferData
+from OpenGL.arrays import ArrayDatatype as ADT
+
+from vector import Vec
+
+
+class VBO:
+#hacking this together because on Mac I can't use the PyOpenGL vbo class
+#mixed with PyOpenCL for some reason
+ def __init__(self, data):
+ self.data = data
+ self.vbo_id = gl.glGenBuffers(1)
+ gl.glBindBuffer(gl.GL_ARRAY_BUFFER, self.vbo_id)
+ rawGlBufferData(gl.GL_ARRAY_BUFFER, ADT.arrayByteCount(data), ADT.voidDataPointer(data), gl.GL_DYNAMIC_DRAW)
+ gl.glBindBuffer(gl.GL_ARRAY_BUFFER, 0)
+
+ def bind(self):
+ gl.glBindBuffer(gl.GL_ARRAY_BUFFER, self.vbo_id)
+
+
+
+def init((width, height)):
+
+ #glEnable(GL_DEPTH_TEST)
+ gl.glEnable(gl.GL_NORMALIZE)
+ gl.glShadeModel(gl.GL_SMOOTH)
+
+
+ gl.glMatrixMode(gl.GL_PROJECTION)
+ gl.glLoadIdentity()
+ gluPerspective(60.0, width/float(height), 1, 1000.0)
+ #glEnable(GL_DEPTH_TEST)
+ gl.glMatrixMode(gl.GL_MODELVIEW)
+
+
+
+
+def lights():
+ gl.glEnable(gl.GL_LIGHTING)
+ gl.glEnable(gl.GL_COLOR_MATERIAL)
+
+ light_position = [10., 10., 200., 0.]
+ light_ambient = [.2, .2, .2, 1.]
+ light_diffuse = [.6, .6, .6, 1.]
+ light_specular = [2., 2., 2., 0.]
+ gl.glLightfv(gl.GL_LIGHT0, gl.GL_POSITION, light_position)
+ gl.glLightfv(gl.GL_LIGHT0, gl.GL_AMBIENT, light_ambient)
+ gl.glLightfv(gl.GL_LIGHT0, gl.GL_DIFFUSE, light_diffuse)
+ gl.glLightfv(gl.GL_LIGHT0, gl.GL_SPECULAR, light_specular)
+ gl.glEnable(gl.GL_LIGHT0)
+
+"""
+ mat_ambient = [.2, .2, 1.0, 1.0]
+ mat_diffuse = [.2, .8, 1.0, 1.0]
+ mat_specular = [1.0, 1.0, 1.0, 1.0]
+ high_shininess = 3.
+
+ mat_ambient_back = [.5, .2, .2, 1.0]
+ mat_diffuse_back = [1.0, .2, .2, 1.0]
+
+ glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient);
+ glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
+ glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
+ glMaterialfv(GL_FRONT, GL_SHININESS, high_shininess);
+
+ glMaterialfv(GL_BACK, GL_AMBIENT, mat_ambient_back);
+ glMaterialfv(GL_BACK, GL_DIFFUSE, mat_diffuse_back);
+ glMaterialfv(GL_BACK, GL_SPECULAR, mat_specular);
+ glMaterialfv(GL_BACK, GL_SHININESS, high_shininess);
+"""
+
+
+def draw_line(v1, v2):
+ gl.glBegin(gl.GL_LINES)
+ gl.glVertex3f(v1.x, v1.y, v1.z)
+ gl.glVertex3f(v2.x, v2.y, v2.z)
+ gl.glEnd()
+
+
+def draw_axes():
+ #X Axis
+ gl.glColor3f(1,0,0) #red
+ v1 = Vec([0,0,0])
+ v2 = Vec([1,0,0])
+ draw_line(v1, v2)
+
+ #Y Axis
+ gl.glColor3f(0,1,0) #green
+ v1 = Vec([0,0,0])
+ v2 = Vec([0,1,0])
+ draw_line(v1, v2)
+
+ #Z Axis
+ gl.glColor3f(0,0,1) #blue
+ v1 = Vec([0,0,0])
+ v2 = Vec([0,0,1])
+ draw_line(v1, v2)
+
+
+
View
118 teach/pycpu/hash.py
@@ -0,0 +1,118 @@
+import math
+from vector import Vec
+
+class Domain(object):
+ def __init__(self, bnd_min, bnd_max):#, surface):
+
+ #the boundary we want particles to stay within
+ self.bnd_min = bnd_min
+ self.bnd_max = bnd_max
+
+ #get the volume of the domain
+ self.width = bnd_max.x - bnd_min.x
+ self.height = bnd_max.y - bnd_min.y
+ self.depth = bnd_max.z - bnd_min.z
+ if self.depth == 0.: self.depth = 1.
+
+ self.V = self.width * self.height * self.depth
+
+ #self.surface = surface
+
+
+ def setup(self, cell_size):
+ #we create 2 cells of padding around the bounds
+ s2 = 2.*cell_size;
+ self.min = self.bnd_min - Vec([s2, s2, s2])
+ self.max = self.bnd_max + Vec([s2, s2, s2])
+
+ self.size = self.max - self.min
+ self.res = Vec([ math.ceil(self.size.x / cell_size),
+ math.ceil(self.size.y / cell_size),
+ math.ceil(self.size.z / cell_size) ])
+
+ self.size = self.res * cell_size
+ self.max = self.min + self.size
+
+ self.delta = Vec([ self.res.x / self.size.x,
+ self.res.y / self.size.y,
+ self.res.z / self.size.z ])
+
+ self.nb_cells = int(self.res.x * self.res.y * self.res.z)
+
+
+ def calc_cell(self, v):
+ """Calculate the grid cell from a vertex"""
+ vv = (v - self.min)*self.delta
+
+ ii = Vec([0,0,0])
+ ii.x = int(vv.x)
+ ii.y = int(vv.y)
+ ii.z = int(vv.z)
+ return ii
+
+ def calc_hash(self, p):
+ """Calculate the grid hash from a grid position"""
+ #tODO: wrapEdges boolean
+ gx = p.x
+ gy = p.y
+ gz = p.z
+ return (gz*self.res.y + gy) * self.res.x + gx;
+
+
+ def draw(self):
+ """draw the lines of the grid"""
+ pass
+
+ def __str__(self):
+ s = "min: %s\n" % self.min
+ s += "max: %s\n" % self.max
+ s += "bnd min: %s\n" % self.bnd_min
+ s += "bnd max: %s\n" % self.bnd_max
+ s += "size: %s\n" % self.size
+ s += "res: %s\n" % self.res
+ s += "delta: %s\n" % self.delta
+ s += "number of cells: %s\n" % self.nb_cells
+ return s
+
+ def make_struct(self, scale_factor):
+
+ print("scale factor", scale_factor)
+ size = self.size * scale_factor
+ min = self.min * scale_factor
+ max = self.max * scale_factor
+ bnd_min = self.bnd_min * scale_factor
+ bnd_max = self.bnd_max * scale_factor
+ res = self.res
+ delta = self.delta / scale_factor
+
+ import struct
+ gpstruct = struct.pack('ffff'+
+ 'ffff'+
+ 'ffff'+
+ 'ffff'+
+ 'ffff'+
+ 'ffff'+
+ 'ffff'+
+ 'i',
+ size.x, size.y, size.z, 0.,
+ min.x, min.y, min.z, 0.,
+ max.x, max.y, max.z, 0.,
+ bnd_min.x, bnd_min.y, bnd_min.z, 0.,
+ bnd_max.x, bnd_max.y, bnd_max.z, 0.,
+ res.x, res.y, res.z, 0.,
+ delta.x, delta.y, delta.z, 0.,
+ self.nb_cells
+ )
+ return gpstruct
+
+
+if __name__ == "__main__":
+
+ #doing some testing on hashing out of bounds
+ dmin = Vec([0,0,0])
+ dmax = Vec([5,5,5])
+ domain = Domain(dmin, dmax)
+
+
+
+
View
70 teach/pycpu/initialize.py
@@ -0,0 +1,70 @@
+import pygame
+from pygame.locals import *
+import numpy as np
+
+from vector import Vec
+from hash import Domain
+from kernels import Kernel
+
+
+
+import numpy as np
+
+def addRect_old(num, pmin, pmax, sphp):
+ #Create a rectangle with at most num particles in it. The size of the return
+ #vector will be the actual number of particles used to fill the rectangle
+ print "**** addRect ****"
+ print "rest dist:", sphp.rest_distance
+ print "sim_scale:", sphp.sim_scale
+ spacing = 1.0 * sphp.rest_distance / sphp.sim_scale;
+ print "spacing", spacing
+
+ xmin = pmin.x# * scale
+ xmax = pmax.x# * scale
+ ymin = pmin.y# * scale
+ ymax = pmax.y# * scale
+
+ print "min, max", xmin, xmax, ymin, ymax
+ rvec = [];
+ i=0;
+ for y in np.arange(ymin, ymax, spacing):
+ for x in np.arange(xmin, xmax, spacing):
+ if i >= num: break
+ print "x, y", x, y
+ rvec += [ Vec([x,y]) * sphp.sim_scale];
+ #rvec += [[x, y, 0., 1.]]
+ i+=1;
+ print "%d particles added" % i
+ #rvecnp = np.array(rvec, dtype=np.float32)
+ #return rvecnp;
+ return rvec
+
+
+
+def init_particles(num, sphp, domain, surface):
+ particles = []
+ p1 = Vec([.5, 2.]) * sphp.sim_scale
+ particles += [ Particle(p1, sphp, [0,0,255], surface) ]
+
+ pmin = Vec([.5, .5])
+ pmax = Vec([2., 3.])
+ ps = addRect_old(num, pmin, pmax, sphp)
+ for p in ps:
+ particles += [ Particle(p, sphp, [255,0,0], surface) ]
+
+ """
+ p2 = Vec([400., 400.]) * sphp.sim_scale
+ p3 = Vec([415., 400.]) * sphp.sim_scale
+ p4 = Vec([400., 415.]) * sphp.sim_scale
+ p5 = Vec([415., 415.]) * sphp.sim_scale
+ p6 = Vec([430., 430.]) * sphp.sim_scale
+ particles += [ Particle(p1, sphp, [255,0,0], surface) ]
+ particles += [ Particle(p2, sphp, [0,0,255], surface) ]
+ particles += [ Particle(p3, sphp, [0,205,0], surface) ]
+ particles += [ Particle(p4, sphp, [0,205,205], surface) ]
+ particles += [ Particle(p5, sphp, [0,205,205], surface) ]
+ particles += [ Particle(p6, sphp, [0,205,205], surface) ]
+ """
+ return particles
+
+
View
113 teach/pycpu/kernels.py
@@ -0,0 +1,113 @@
+from vector import Vec
+
+from math import *
+import numpy as np
+from numpy import pi
+
+def mag(x):
+ return sqrt(np.dot(x, x))
+def mag2(x):
+ return np.dot(x,x)
+
+
+class Kernel(object):
+ def __init__(self, h):
+ self.h = h
+ h6 = h**6
+ h9 = h6*h*h*h
+ self.coeffs = {}
+ self.coeffs["poly6"] = 315./(64.*pi*h9)
+ self.coeffs["dspiky"] = -45./(pi*h6)
+ self.coeffs["ddvisc"] = 45./(pi*h6)
+
+ def poly6(self, r):
+ mr2 = mag2(r)
+ if mr2 < self.h*self.h:
+ return self.coeffs["poly6"] * (self.h**2 - mr2)**3
+ else:
+ return 0.
+
+ def dspiky(self, r):
+ magr = mag(r)
+ if magr == 0:
+ magr = 1E-6
+
+ hr2 = self.h - magr
+
+ if magr < self.h:
+ return self.coeffs["dspiky"] * hr2 * hr2 / magr
+ return 0
+
+ def ddvisc(self, r):
+ magr = mag(r)
+ if magr == 0:
+ magr = 1E-6
+
+ hr = self.h - magr
+
+ if magr < self.h:
+ return self.coeffs["ddvisc"] * hr
+ return 0
+
+
+
+def Wpoly6(h, r):
+ coeff = 315./(64.*pi*h**9)
+ #magr = abs(r)
+ magr = mag(r)
+
+ if magr < h:
+ return coeff*(h**2 - magr**2)**3
+ return 0
+
+def Wspiky(h, r):
+ coeff = 15./(pi*h**6)
+ #magr = abs(r)
+ magr = mag(r)
+
+ if magr < h:
+ return coeff*(h-magr)**3
+ return 0
+
+def dWspiky(h, r):
+ """ still need to multiply this quantity by r """
+ #magr = abs(r)
+ magr = mag(r)
+ #print "magr", magr
+ if magr == 0:
+ magr = 1E-6
+
+ den = magr * pi * h**6
+ coeff = -45./den
+ hr2 = h - magr
+
+ if magr < h:
+ return coeff*hr2*hr2
+ return 0
+
+def main():
+ import numpy as np
+ import pylab
+
+ h = 1
+ k = Kernel(h)
+
+ X = np.linspace(-h, h, 100)
+ Y = []
+ dY = []
+ for x in X:
+ #Y += [Wpoly6(h, [0,x])]
+ Y += [Wspiky(h,[0, x])]
+ #dY += [dWspiky(h, [0,x])]
+ #dY += [k.dspiky([0,x])]
+ dY += [k.ddvisc([0,x])]
+
+ #pylab.plot(X,Y)
+ pylab.plot(X,dY)
+ pylab.show()
+
+
+if __name__ == "__main__":
+ main()
+
+
View
168 teach/pycpu/main.py
@@ -0,0 +1,168 @@
+import pygame
+from pygame.locals import *
+
+from forces import *
+from vector import Vec
+import sph
+#import magnet
+from particle import Particle
+from hash import Domain
+
+dt = .001
+
+
+def fromscreen(p, surface):
+ #v.x
+ p.y = surface.get_height() - p.y
+ return p
+
+
+def toscreen(p, surface, screen_scale):
+ translate = Vec([0,0])
+ p.x = translate.x + p.x*screen_scale
+ p.y = surface.get_height() - (translate.y + p.y*screen_scale)
+ return p
+
+
+@timings
+def draw_particles(ps):
+ ps.reverse()
+ i = 0
+ for p in ps:
+ if i == len(ps)-1:
+ p.draw(True)
+ else:
+ p.draw()
+ i += 1
+ ps.reverse()
+
+
+def main():
+ pygame.init()
+ screen = pygame.display.set_mode((800, 800))
+ pygame.display.set_caption('SPH Forces')
+
+ background = pygame.Surface(screen.get_size())
+ background = background.convert()
+ background.fill((250, 250, 250))
+
+ clock = pygame.time.Clock()
+
+ max_num = 2**12 #4096
+ #max_num = 2**10 #1024
+ #max_num = 2**8 #256
+ #max_num = 2**7 #128
+
+ dmin = Vec([0,0,0])
+ dmax = Vec([50,50,50])
+ domain = Domain(dmin, dmax)#, screen)
+ system = sph.SPH(max_num, domain)
+ #system = magnet.Magnet(max_num, domain)
+
+ #particles = sph.init_particles(5, system, domain, screen)
+ particles = []
+ p1 = Vec([25, 25]) * system.sim_scale
+ np = Vec([.8, .8])
+ particles += [ Particle(p1, np, system, [0,128,128], screen) ]
+ particles[0].lock = True
+
+ p = Vec([25 + 8, 25]) * system.sim_scale
+ particles += [ Particle(p, np, system, [128,128,0], screen) ]
+ #"""
+ p = Vec([25, 25 + 8]) * system.sim_scale
+ particles += [ Particle(p, np, system, [128,128,0], screen) ]
+ p = Vec([25 - 8, 25]) * system.sim_scale
+ particles += [ Particle(p, np, system, [128,128,0], screen) ]
+ p = Vec([25, 25 - 8]) * system.sim_scale
+ particles += [ Particle(p, np, system, [128,128,0], screen) ]
+ p = Vec([25 + 8, 25 + 8]) * system.sim_scale
+ particles += [ Particle(p, np, system, [128,128,0], screen) ]
+ #"""
+
+
+
+ print "p0.pos:", particles[0].pos
+
+
+ mouse_down = False
+ pause = True
+ pi = 1
+ tpi = pi
+ while 1:
+ tpi = pi
+ tcol = particles[pi].col[:]
+ particles[pi].col = [0, 200, 0]
+ clock.tick(60)
+ key = pygame.key.get_pressed()
+ for event in pygame.event.get():
+ if event.type == QUIT or key[K_ESCAPE] or key[K_q]:
+ print "quit!"
+ return
+ elif key[K_t]:
+ print timings
+
+ elif key[K_0]:
+ pi = 0
+ elif key[K_1]:
+ pi = 1
+ elif key[K_2]:
+ pi = 2
+ elif key[K_3]:
+ pi = 3
+ elif key[K_4]:
+ pi = 4
+ #elif key[K_5]:
+ # pi = 5
+
+ elif key[K_p]:
+ pause = not pause
+
+ elif event.type == MOUSEBUTTONDOWN:
+ if pygame.mouse.get_pressed()[2]:
+ v = Vec([event.pos[0], event.pos[1]])
+ v = fromscreen(v, screen)
+ p = Particle([0,0], [0,0], system, [128,128,0], screen)
+ p.move(v)
+ particles += [ p ]
+ break
+ tcol = particles[pi].col[:]
+ mouse_down = True
+ elif event.type == MOUSEMOTION:
+ if(mouse_down):
+ v = Vec([event.pos[0], event.pos[1]])
+ v = fromscreen(v, screen)
+ print v
+ particles[pi].move(v)
+ particles[pi].vel = Vec([0.,0.])
+ particles[pi].force = Vec([0.,0.])
+ particles[pi].angular = Vec([0.,0.])
+ elif event.type == MOUSEBUTTONUP:
+ particles[pi].vel = Vec([0.,0.])
+ particles[pi].force = Vec([0.,0.])
+ particles[pi].angular = Vec([0.,0.])
+ mouse_down = False
+
+
+ screen.blit(background, (0, 0))
+
+ density_update(system, particles)
+ if not pause:
+ for i in range(10):
+ #magnet_update(system, particles)
+ #density_update(system, particles)
+ force_update(system, particles)
+
+ #contact_update(system, particles)
+ collision_wall(system, domain, particles)
+ #euler_update(system, particles, dt)
+ leapfrog_update(system, particles)
+
+ draw_particles(particles)
+ pygame.display.flip()
+
+ particles[tpi].col = tcol[:]
+
+
+
+if __name__ == "__main__":
+ main()
View
71 teach/pycpu/particle.py
@@ -0,0 +1,71 @@
+import pygame
+from pygame.locals import *
+import numpy as np
+
+from vector import Vec
+from hash import Domain
+from kernels import Kernel
+
+
+def toscreen(p, surface, screen_scale):
+ translate = Vec([0,0])
+ p.x = translate.x + p.x*screen_scale
+ p.y = surface.get_height() - (translate.y + p.y*screen_scale)
+ return p
+
+class Particle:
+ def __init__(self, pos, np, system, color, surface):
+ #physics stuff
+ #position
+ self.pos = pos
+ #north pole (vector with origin at position)
+ #self.np = np
+ #south pole (vector with origin at position)
+ #self.sp = -1 * np
+
+ #self.r = system.radius
+ self.h = system.smoothing_radius
+ self.scale = system.sim_scale
+ self.mass = system.mass
+ self.dens = system.rho0
+
+ self.force = Vec([0., 0.])
+ self.vel = Vec([0.,0.])
+ self.veleval = Vec([0.,0.])
+
+ #lock a particle in place (force updates don't affect it)
+ self.lock = False
+
+ #pygame stuff
+ self.col = color
+ self.surface = surface
+ self.screen_scale = self.surface.get_width() / system.domain.width
+
+ def move(self, pos):
+ self.pos = pos * self.scale / self.screen_scale
+
+ #print "dens", self.dens
+
+ def draw(self, show_dense = False):
+ #draw circle representing particle smoothing radius
+
+ dp = toscreen(self.pos / self.scale, self.surface, self.screen_scale)
+ dp = [int(dp.x), int(dp.y)]
+ r = self.screen_scale * self.h / self.scale
+ #print dp, r
+ #pygame.draw.circle(self.surface, self.col, dp, r, 1)
+ #pygame.draw.circle(self.surface, self.col, dp, int(r), 1)
+ #draw filled circle representing particle size
+ pygame.draw.circle(self.surface, self.col, dp, int(self.dens / 40.), 0)
+
+ #dnp = toscreen(self.pos + self.np / self.scale, self.surface, self.screen_scale)
+ #dsp = toscreen(self.pos + self.sp / self.scale, self.surface, self.screen_scale)
+ #pygame.draw.circle(self.surface, [200,0,0], dnp, 5)
+ #pygame.draw.circle(self.surface, [0,0,200], dsp, 5)
+
+
+ #TODO draw force vector (make optional)
+ #vec = [self.x - f[0]*fdraw/fscale, self.y - f[1]*fdraw/fscale]
+ #pygame.draw.line(self.surface, pj.col, self.pos, vec)
+
+
View
110 teach/pycpu/sph.py
@@ -0,0 +1,110 @@
+import pygame
+from pygame.locals import *
+import numpy as np
+
+from vector import Vec
+from hash import Domain
+from kernels import Kernel
+
+
+class SPH:
+ def __init__(self, max_num, domain):
+ self.max_num = max_num
+
+ rho0 = 1000. #rest density [ kg/m^3 ]
+ VF = .0262144 #simulation volume [ m^3 ]
+ VP = VF / max_num #particle volume [ m^3 ]
+ m = rho0 * VP #particle mass [ kg ]
+ re = (VP)**(1/3.) #particle radius [ m ]
+ #re = (VP)**(1/2.) #particle radius [ m ]
+ print "re, m, VP", re, m, VP
+ rest_distance = .87 * re #rest distance between particles [ m ]
+
+ smoothing_radius = 2.0 * rest_distance #smoothing radius for SPH Kernels
+ boundary_distance = .5 * rest_distance #for calculating collision with boundary
+
+ #the ratio between the particle radius in simulation space and world space
+ print "VF", VF
+ print "domain.V: ", domain.V
+ print "VF/domain.V", VF/domain.V
+ print "scale calc", (VF/domain.V)**(1/3.)
+ #print "scale calc", (VF/domain.V)**(1/2.)
+ sim_scale = (VF / domain.V)**(1/3.) #[m^3 / world m^3 ]
+ #sim_scale = (VF / domain.V)**(1/2.) #[m^2 / world m^2 ]
+
+ self.rho0 = rho0
+ self.VF = VF
+ self.mass = m
+ self.VP = VP
+ self.re = re
+ self.rest_distance = rest_distance
+ self.smoothing_radius = smoothing_radius
+ self.boundary_distance = boundary_distance
+ self.sim_scale = sim_scale
+
+ print "====================================================="
+ print "particle mass:", self.mass
+ print "Fluid Volume VF:", self.VF
+ print "simulation scale:", self.sim_scale
+ print "smoothing radius:", self.smoothing_radius
+ print "rest distance:", self.rest_distance
+ print "====================================================="
+
+ #Other parameters
+ self.K = 15. #Gas constant
+ self.boundary_stiffness = 20000.
+ self.boundary_dampening = 256.
+ #friction
+ self.friction_coef = 0.
+ self.restitution_coef = 0.
+ #not used yet
+ self.shear = 0.
+ self.attraction = 0.
+ self.spring = 0.
+
+ self.velocity_limit = 600.
+ self.xsph_factor = .05
+
+ self.viscosity = .01
+ self.gravity = -9.8
+
+ self.EPSILON = 10E-6
+
+ #Domain
+ self.domain = domain
+ self.domain.setup(self.smoothing_radius / self.sim_scale)
+
+ #Kernels
+ self.kernels = Kernel(self.smoothing_radius)
+
+ def make_struct(self, num):
+ import struct
+ sphstruct = struct.pack('ffff'+
+ 'ffff'+
+ 'ffff'+
+ 'ffff'+
+ 'ffff'+
+ 'ffff'+
+ 'ffff'+
+ 'iiii',
+ self.mass, self.rest_distance, self.smoothing_radius, self.sim_scale,
+ self.boundary_stiffness, self.boundary_dampening, self.boundary_distance, self.K,
+ self.viscosity, self.velocity_limit, self.xsph_factor, self.gravity,
+ self.friction_coef, self.restitution_coef, self.shear, self.attraction,
+ self.spring, self.EPSILON, np.pi, self.kernels.coeffs["poly6"],
+ 0., 0., 0., self.kernels.coeffs["dspiky"], #some kernels not used so coeffs
+ 0., 0., 0., self.kernels.coeffs["ddvisc"], #not calculated right now
+ num, 0, 0, self.max_num #nb_vars and coice not used
+ )
+ return sphstruct
+
+
+
+
+if __name__ == "__main__":
+ dmin = Vec([0.,0.,0.])
+ dmax = Vec([500.,500.,1.])
+ domain = Domain(dmin, dmax)
+ system = SPH(2**14, domain) #16384
+ #system = SPH(2**12, domain) #4096
+
View
63 teach/pycpu/test_hash.py
@@ -0,0 +1,63 @@
+from vector import Vec
+from hash import Domain
+import sph
+
+
+def print_hash(domain, v):
+ print "position", v
+ gp = domain.calc_cell(v)
+ print "grid cell", gp
+ print "grid hash", domain.calc_hash(gp)
+
+
+if __name__ == "__main__":
+
+ #doing some testing on hashing out of bounds
+ dmin = Vec([0.,0.,0.])
+ dmax = Vec([5.,5.,5.])
+ domain = Domain(dmin, dmax, None)
+
+ max_num = 8192
+ system = sph.SPH(max_num, domain)
+
+ dmin_s = dmin * system.sim_scale
+ dmax_s = dmax * system.sim_scale
+ domain_s = Domain(dmin_s, dmax_s, None)
+ domain_s.setup(system.smoothing_radius)
+
+ #print domain_s
+ print domain
+
+ print "=================="
+
+ v = Vec([1., 1., 1.])
+ print_hash(domain, v)
+
+ v = Vec([4.9, 4.9, 4.9])
+ print_hash(domain, v)
+
+ v = Vec([4.9999, 4.9999, 4.9999])
+ print_hash(domain, v)
+
+ print "out of positive bounds"
+ v = Vec([6., 6., 6.])
+ print_hash(domain, v)
+
+ v = Vec([6.5, 6.5, 6.5])
+ print_hash(domain, v)
+
+ v = Vec([7., 7., 7.])
+ print_hash(domain, v)
+
+ print "out of negative bounds"
+ #these are not accurate because hash should be unsigned int
+ v = Vec([-1., -1., -1.])
+ print_hash(domain, v)
+
+ v = Vec([-2., -2., -2.])
+ print_hash(domain, v)
+
+
+
+
+
View
126 teach/pycpu/timing.py
@@ -0,0 +1,126 @@
+
+import time
+from time import clock as clk
+
+def print_timing(func):
+ def wrapper(*arg):
+ t1 = time.time()
+ res = func(*arg)
+ t2 = time.time()
+ #print "%s took %f ms with args: %s" % (func.func_name, (t2-t1)*1000.0, arg)
+ print "%s took %f ms" % (func.func_name, (t2-t1)*1000.0)
+ return res
+ return wrapper
+
+#----------------------------------------------------------------------
+
+def timing_collector(out):
+ while True:
+ (name, t) = (yield)
+ if name in out:
+ out[name] += [t]
+ else:
+ out[name] = [t]
+
+
+class Timing(object):
+ def __init__(self):#,cor):
+ self.timings = {}
+ self.col = self.__collector()
+ self.col.next() #coroutine syntax
+
+ def __collector(self):
+ while True:
+ (name, t) = (yield) #coroutine syntax
+ if name in self.timings:
+ self.timings[name]["timings"] += [t]
+ self.timings[name]["count"] += 1
+ self.timings[name]["total"] += t
+ else:
+ self.timings[name] = {} #if this entry doesn't exist yet
+ self.timings[name]["timings"] = [t]
+ self.timings[name]["count"] = 1
+ self.timings[name]["total"] = t
+
+ def __call__(self, func):
+ """Turn the object into a decorator"""
+ def wrapper(*arg, **kwargs):
+ t1 = time.time() #start time
+ res = func(*arg, **kwargs) #call the originating function
+ t2 = time.time() #stop time
+ t = (t2-t1)*1000.0 #time in milliseconds
+ data = (func.__name__, t)
+ self.col.send(data) #collect the data
+ return res
+ return wrapper
+
+ def __str__(self):
+ s = "Timings:\n"
+ print dir(self)
+ keys = self.timings.keys()
+ for key in sorted(keys):
+ s += "%s | " % key
+ ts = self.timings[key]["timings"]
+ count = self.timings[key]["count"]
+ total = self.timings[key]["total"]
+ s += "average: %s | total: %s | count: %s\n" % (total / count, total, count)
+ return "%s" % s
+
+#----------------------------------------------------------------------
+
+
+class Timer:
+ def tic(self):
+ self.time = clk()
+
+ def toc(self, msg=""):
+ self.time = clk() - self.time
+ if msg:
+ print "(%s), time: %f sec" % (msg, 1000.*self.time)
+ else:
+ print "time: %f sec" % (1000.*self.time)
+
+#----------------------------------------------------------------------
+
+
+
+
+if __name__ == "__main__":
+
+
+ timings = Timing()
+ #@print_timing
+
+ @timings
+ def add(x,y):
+ for i in range(10000):
+ c = x + y
+ return c
+
+ #@print_timing
+ @timings
+ def multiply(x,y):
+ for i in range(10000):
+ c = x * y
+ return c
+
+
+ for i in range(100):
+ add(3.,4.)
+ multiply(3., 4.)
+
+ print timings
+
+ """
+ t = Timer()
+ x = 3.
+
+ t.tic()
+ add(5.,7.)
+ t.toc()
+
+ t.tic()
+ for i in range(1000):
+ x = log(2.+cos(x))
+ t.toc()
+ """
View
103 teach/pycpu/vector.py
@@ -0,0 +1,103 @@
+
+
+#Vector class
+#Author - Ian Johnson | enjalot@gmail.com
+#http://docs.scipy.org/doc/numpy/user/basics.subclassing.html#slightly-more-realistic-example-attribute-added-to-existing-array
+
+import numpy as np
+import math
+
+class Vec(np.ndarray):
+ props = ['x', 'y', 'z', 'w']
+
+ def __new__(cls, input_array):
+ # Input array is an already formed ndarray instance
+ # We first cast to be our class type
+ obj = np.asarray(input_array).view(cls)
+ # add the new attribute to the created instance
+ if len(obj) < 2 or len(obj) > 4:
+ #not a 2,3 or 4 element vector!
+ return None
+ return obj
+
+ def __array_finalize__(self, obj):
+ #this gets called after object creation by numpy.ndarray
+ #print "array finalize", obj
+ if obj is None: return
+ #we set x, y etc properties by iterating through ndarray
+ for i in range(len(obj)):
+ setattr(self, Vec.props[i], obj[i])
+
+ def __array_wrap__(self, out_arr, context=None):
+ #this gets called after numpy functions are called on the array
+ #out_arr is the output (resulting) array
+ #oa = Vec(out_arr)
+ for i in range(len(out_arr)):
+ #setattr(oa, Vec.props[i], oa[i])
+ setattr(out_arr, Vec.props[i], out_arr[i])
+ #out_arr = oa
+ return np.ndarray.__array_wrap__(self, out_arr, context)
+ #return np.ndarray.__array_wrap__(self, oa, context)
+
+ def __repr__(self):
+ desc="""Vec2(data=%(data)s,""" # x=%(x)s, y=%(y)s)"""
+ for i in range(len(self)):
+ desc += " %s=%s," % (Vec.props[i], getattr(self, Vec.props[i]))
+ desc = desc[:-1] #cut off last comma
+ return desc % {'data': str(self) }
+
+
+ def __setitem__(self, ind, val):
+ #print "SETITEM", self.shape, ind, val
+ if self.shape == (): #dirty hack for dot product
+ return
+ self.__dict__[Vec.props[ind]] = val
+ return np.ndarray.__setitem__(self, ind, val)
+
+ def __setattr__(self, item, val):
+ self[Vec.props.index(item)] = val
+
+
+ ###
+ #math utilities
+ ###
+
+def normalize(u):
+ return u / (math.sqrt(np.dot(u, u)))
+
+
+if __name__ == "__main__":
+ v1 = Vec(np.arange(2))
+ v2 = Vec(np.arange(2))
+ print "v1:", repr(v1)
+ print "v2:", repr(v2)
+ v3 = v1 + v2
+ print "v3:", repr(v3)
+ v1[0] = 5
+ print "v1.x, v1[0]:", v1.x, v1[0]
+ v1.y = 6
+ print "v1.y, v1[1]:", v1.y, v1[1]
+ print "v1:", repr(v1)
+
+ print "====================="
+ print "v1.shape", v1.shape
+ print "multiply(v1,v1):", np.multiply(v1,v1)
+ print "dot(v1,v1):", np.dot(v1, v1)
+
+ print "3D ================"
+ va = Vec([0,0,0])
+ print "va: ", repr(va)
+ va.x = 2
+ va.y = 9
+ va.z = 5
+ print "va: ", repr(va)
+ vb = Vec([1,1,1])
+ vc = va + vb
+ print "v3: ", repr(vc)
+ print "dot(va, vb):", np.dot(va, vb)
+
+
+
+
+
+

0 comments on commit 6124e88

Please sign in to comment.
Something went wrong with that request. Please try again.