emscripten
Downloading...
Resize canvas Lock/hide mouse pointer    
This is an SPH/lagrangian fluid simulation. First a local density is computed for each particle based on its denstiance to its nearby neighbors. A pressure is then computed for each particle based on its local density: densities over the target density are positive and those under the target negative. Particles with a positive pressure are pushed away from their neighbors (and objects with a negative density are drawn to their neighbors) with a force proportional to the magnitude of their pressure.

Following the approach from Clavet et al. particles experience two pressures, one is the normal SPH pressure that tries to achieve rest density, the other is based on a special "near" density that rapidly decreases with distance. This "near" pressure always pushes particles away from each other, which results in more coherent structures with less clustering. 

 import "Vec2.plai"
import "Draw.plai"

width = 800
aspect_ratio = 5/8
init_window("2D Rigid Body Dynamics",width,width*aspect_ratio)

k_smooth_radius = 0.035
k_stiff = 150.0
k_stiffN = 1000.0
k_rest_density = 0.2
grab_radius = 0.08

struct Particle
  pos, oldPos, vel
  press, dens
  pressN, densN
  grabbed

def Particle(x,y)
  Particle{Vec2(x,y), Vec2(x,y), Vec2(0,0), 0.0, 0.0, 0.0, 0.0, false}

particles = [:Particle]
r = 35
n = 40
numParticles = 0
for i in 0...r
  for j in 0...n
    particles << Particle(-1+2*j/n+0.1*i/r,aspect_ratio*(0.5 - 1.5*i/r))
    numParticles += 1

struct Pair
  p1, p2      #particles in the pair
  q, q2, q3   #q is (1-d/h), q2 is q^2, q3 is q^3

def Pair(p1, p2, q)
  Pair{p1, p2, q, q*q, q*q*q}

def simulateSPH(particles,numParticles,k_smooth_radius,k_rest_density,k_stiff,k_stiffN,mousePos,grab_radius,aspect_ratio,dt)
  #Integrate velocity & gravity to update positions
  for p in particles
    p.vel = (p.pos - p.oldPos)/dt   #Compute vel from positions updated via pressure forces
    p.vel += Vec2(0,-5)*dt          #Integrate velocity based on gravity
    
    if p.pos.y < -aspect_ratio      #Collision response
      p.pos.y = -aspect_ratio
      p.vel.y *= -0.3
    if p.pos.x < -1.0
      p.pos.x = -1.0
      p.vel.x *= -0.3
    if p.pos.x > 1.0
      p.pos.x = 1.0
      p.vel.x *= -0.3

    p.vel += 10*dt*((mousePos - p.pos)/grab_radius - p.vel) if p.grabbed #Move grabbed particles towards mouse
    
    p.oldPos = p.pos
    p.pos += p.vel*dt              #Integrate position based on velocity
    p.dens = 0.0
    p.densN = 0.0

  #Find all neighboring particles [TODO: Slow for large SPH systems!]
  pairs = [:Pair]
  for i in 0...numParticles
    for j in 0...numParticles
      dist = particles[i].pos.distanceTo(particles[j].pos)
      if dist < k_smooth_radius and i < j
        q = 1 - (dist/k_smooth_radius)
        pairs << Pair(ref particles[i], ref particles[j], q)

  #Accumulate per-particle density
  for p in pairs
    p.p1.dens += p.q2
    p.p2.dens += p.q2
    p.p1.densN += p.q3
    p.p2.densN += p.q3

  #Computer per-particle pressure: stiffness*(density - rest_density)
  for p in particles
    p.press = k_stiff*(p.dens - k_rest_density)
    p.pressN = k_stiffN*(p.densN)
    if p.press > 30 p.press = 30      #maximum pressure
    if p.pressN > 300 p.pressN = 300  #maximum near pressure

  for pair in pairs
    a = pair.p1
    b = pair.p2
    total_pressure = (a.press + b.press) * pair.q + (a.pressN + b.pressN) * pair.q2
    displace = total_pressure * (dt**2)
    a.pos += displace * (a.pos - b.pos).normalized
    b.pos += displace * (b.pos - a.pos).normalized

last_frame_start = now()
while window_is_open()
  dt = duration(last_frame_start)/1000 #ms to s
  last_frame_start = now()

  #Note which particles to grab/relase based on mouse
  mouse = getMouse()
  mousePos = Vec2(mouse.x,mouse.y)
  if mouse.clicked
    for p in particles p.grabbed = true if mousePos.distanceTo(p.pos) < grab_radius
  if mouse.released
    p.grabbed = false for p in particles 


  #Substep simulation for increased numerical stability
  for 0..10
    sim_dt = dt/10
    sim_dt = 0.003 if sim_dt > 0.003 #enforce maximum timestep
    simulateSPH(particles,numParticles,k_smooth_radius,k_rest_density,k_stiff,k_stiffN,mousePos,grab_radius,aspect_ratio,sim_dt)

  #Draw the background & SPH particles
  color(0,0,0)
  drawBox(0,0,2,2*aspect_ratio,0)
  for p in particles
    q = (p.press)/30                     #normalize pressure to 1
    color(0.7-q*0.5,0.8-q*0.4,1.0-q*0.2) #color particles by pressure (white->blue)
    drawCircle(p.pos.x,p.pos.y,0.012)    #x,y,radius

  drawFrame()