How to simulate water in 2D physics with codea?

[littfass]

Looking for bestpractice advices on how to simulate water for interaction with contaos physical bodies.
Eg. Material point method for fluid simulation
Would like to implement the following with Lua:
http://www.youtube.com/watch?v=IQCPuVOOq2o
This was originally implemented in java by http://grantkot.com/MPM/Liquid.html

Tried to port this script from processing to Lua, but failed dramatically.
Here the work from Golan Levin that I tried to port.
https://github.com/golanlevin/MPM-Fluid

The main porting work is done, but there are some caveats in the logic.
If one is interested to help i’ll try to post the code here (is there some way to share projects via this forum?)
Should be useful for all of us, so maybe we can put our brains together to solve this.

[littfass]

Edited main post

[ruilov]

having that interacting with gravity would be pretty cool. I wonder, can you just use the physics api and lots of tiny bodies to get a good effect?

[littfass]

Sorry, I guess not… Too slow. It might need some kind of particles logic handling a hole grid of bodies but not sure how to achieve this within the physics engine.

Yeah, interacting with gravity should rock!

[Bri_G]

Hi @lttlfas,

I had trouble porting code from processing and I think it was because it was using a library. When the library was removed on processing I got the same result as with my own code in Codea.

Try posting your code - there may be something obvious, but check the library situation first. If you could point us to the processing code that may help too.

Bri_G

[Bri_G]

Hi @lttlfas,

Ignore my last comment found the processing code.

Bri_G

[Connorbot999]

On your code the fps is?

[littfass]

@Bri_G
Didnt use any other library. But I had to fake the map() and millis() methods from processing.
Maybe thats the point do following functions make sense in Lua, are they equivalent to their respective processing methods?

function millis()
    return ElapsedTime * 1000
end

function map(value,low1,high1,low2,high2)
  if (value < low1) or (value > high1) then
      return nil
  end
  local srcMax = high1 - low1
  local dstMax = high2 - low2
  local adjValue = value - low1
  return (adjValue * dstMax / srcMax) + low2
end

@Connorbot999
Difficult to say, the particles grid (table) is called out of bounds after half a second and throws an error

How would you check for framerate? Is that the right way to get it:

--setup
time_start=ElapsedTime
--draw
FPS = frames / ( ElapsedTime - time_start)
frames=frames+1

[littfass]

Wow, will have to have a look on how to pretty print code on the forums, sorry
@Bri_G
Here’s the code for now, but i guess it’s useless in this form. Still looking for how to print the code in a more readable form.

[littfass]

–Here the code: init settings and grid

MPM = class()

function MPM:init()

    self.settings = {}
    self.settings[0]=2.5 --"Density" 2.5
    self.settings[1]=0.5 --"Stiffness" 0.5
    self.settings[2]=3.0--"Bulk Viscosity" 3.0
    self.settings[3]=1 --"Elasticity" 1.0
    self.settings[4]=1 --"Viscosity" 1.0
    self.settings[5]=1 --"Yield Rate" 1.0
    self.settings[6]=0.02 --"Gravity" //0.02 to 0.004
    self.settings[7]=1 --"Smoothing" 1.0
        
    self.particles = {}
    self.particlesX=0
    self.particlesY=0
    self.nParticles=0

    self.grid={} -- with [][]
    self.activeNodeArray = {}
    self.nActiveNodes=0

    self.gsizeX=0
    self.gsizeY=0

    self.prevMousePressed = false
    self.prevMouseX = 0
    self.prevMouseY = 0
    self.mode = -1
    self.SC = 10 --scaling
    self.timer=0
    self.np =10
    self.steps=0
    self:setup()
end

function MPM:setup()
  self.width=WIDTH
  self.height=HEIGHT

  self.np =15
  self:initLiquid(1+math.floor(self.width/self.SC), 1+math.floor(self.height/self.SC), self.np,self.np)
  self:initGridAndParticles()
end

function MPM:draw()
  background(35, 45, 56, 255)
  self:simulate()
  self:render()
end

function MPM:initLiquid(ingsizeX, ingsizeY, inparticlesX, inparticlesY)
  self.gsizeX     = ingsizeX
  self.gsizeY     = ingsizeY
  self.particlesX = inparticlesX
  self.particlesY = inparticlesY

  self.prevMousePressed = false
  self.prevMouseX = 0
  self.prevMouseY = 0
end

function MPM:initGridAndParticles()
  self.grid = {} --new Node[gsizeX][gsizeY]
  for i=0,self.gsizeX do
    self.grid[i] = {}
    for j=0,self.gsizeY do
      self.grid[i][j] = Node()
    end
  end

  self.activeNodeArray = {} --new Node[gsizeX * gsizeY];
  self.nActiveNodes = 0


  self.nParticles = self.particlesX * self.particlesY
  self.particles = {} --new Particle [nParticles]; 
  local np = 0
  local fi = self.gsizeX / (self.particlesX+5)
  for i=0,self.particlesX do
    for j=0,self.particlesY do
      local px = (i+1)*fi
      local py = map(j, 0,self.particlesY-1, 2,self.gsizeY/3)
      local p = Particle(px, py, 0.0, 0.0)
      self.particles[np] = p
      np = np + 1
    end
  end
end

function MPM:render()
  noSmooth()
  --stroke(0, 255, 239, 255)
  --beginShape(POINTS)
fill(0, 211, 255, 255)
  for ip=0,self.nParticles do
    local p = self.particles[ip]
    ellipse(self.SC*p.x, self.SC*p.y,10,10) --vertex()
  end
  --endShape()
end

[littfass]

– simulate function part 1

function MPM:simulate ()
        
  local mouseX=CurrentTouch.x
  local mouseY=CurrentTouch.y
  local drag = false
  local mdx = 0.0
  local mdy = 0.0
  local mx = mouseX/self.SC
  local my = mouseY/self.SC

  local phi=0
  local mousePressed=(not CurrentTouch.state == ENDED)

  if mousePressed and self.prevMousePressed then
    drag = true
    mdx = (mouseX - self.prevMouseX)/self.SC
    mdy = (mouseY - self.prevMouseY)/self.SC
  end

  self.prevMousePressed = mousePressed
  self.prevMouseX = mouseX
  self.prevMouseY = mouseY

  for i=0,self.gsizeX do
    for j=0,self.gsizeY do
      self.grid[i][j]:clear()
    end
  end
  self.nActiveNodes = 0

  --Particles pass 1
  local pcxTmp=0
  local pcyTmp=0
  for ip=0,self.nParticles do

    local p = self.particles[ip]

    pcxTmp = math.floor(p.x - 0.5)
    pcxTmp = math.min(self.gsizeX - 4, math.max(0, pcxTmp))
    local pcx = pcxTmp
    p.cx = pcxTmp

    pcyTmp = math.floor(p.y - 0.5)
    pcyTmp = math.min(self.gsizeY - 4, math.max(0, pcyTmp))
    local pcy = pcyTmp
    p.cy = pcyTmp

    local px = p.px
    local py = p.py
    local gx = p.gx
    local gy = p.gy
    local pu = p.u
    local pv = p.v

    local x = p.cx - p.x
    px[0] = 0.5 * x * x + 1.5 * x + 1.125
    gx[0] = x + 1.5
    x = x + 1
    px[1] = -x * x + 0.75
    gx[1] = -2 * x
    x = x + 1
    px[2] = (0.5 * x * x - 1.5 * x) + 1.125
    gx[2] = x - 1.5

    local y = p.cy - p.y
    py[0] = 0.5 * y * y + 1.5 * y + 1.125
    gy[0] = y + 1.5
    y = y + 1
    py[1] = -y * y + 0.75
    gy[1] = -2 * y
    y = y + 1
    py[2] = (0.5 * y * y - 1.5 * y) + 1.125
    gy[2] = y - 1.5

    for i = 0,3 do
      local nrow = self.grid[pcx + i] -- potential for array index out of bounds here if simulation explodes.

      local pxi = px[i]
      local gxi = gx[i]

      for j = 0,3 do
        local n = nrow[pcy + j] -- potential for array index out of bounds here if simulation explodes.

        if not n.active then
          n.active = true
          self.activeNodeArray[self.nActiveNodes] = n
          self.nActiveNodes = self.nActiveNodes + 1
        end
        phi   = pxi * py[j]
        n.m  = n.m + phi
        n.gx = n.gx +  gxi * py[j]
        n.gy = n.gy +  pxi * gy[j]
        n.u = n.u +  phi * pu
        n.v = n.v +  phi * pv
      end
    end
  end


  for ni=0,self.nActiveNodes do
    local n = self.activeNodeArray[ni]
    if n and (n.m > 0) then
      n.u = n.u / n.m
      n.v = n.u / n.m
    end
  end

[littfass]

– simulate function part 2

-- Particles pass 2
  local densitySetting = self.settings[0]
  local stiffness      = self.settings[1]
  local bulkViscosity  = self.settings[2]
  local elasticity     = self.settings[3]
  local viscosity      = self.settings[4]
  local yieldRate      = self.settings[5]

  local T = millis()/1500
  self.timer=T
  for ip=0,self.nParticles do
    local p = self.particles[ip]
    local px = p.px
    local py = p.py
    local gx = p.gx
    local gy = p.gy
    local pcy = p.cy
    local pcx = p.cx

    local dudx = 0.0
    local dudy = 0.0
    local dvdx = 0.0
    local dvdy = 0.0

    local gxi=0
    local pxi=0
    local gxf=0
    local gyf=0
    for i = 0,3 do
      local nrow = self.grid[pcx + i]
      gxi = gx[i]
      pxi = px[i]

      local n0 = nrow[pcy]
      gxf = gxi * py[0]
      gyf = pxi * gy[0]
      dudx = dudx +  n0.u * gxf
      dudy = dudy +  n0.u * gyf
      dvdx = dvdx +  n0.v * gxf
      dvdy = dvdy +  n0.v * gyf

      local n1 = nrow[pcy+1]
      gxf = gxi * py[1]
      gyf = pxi * gy[1]
      dudx = dudx +  n1.u * gxf
      dudy = dudy +  n1.u * gyf
      dvdx = dvdx +  n1.v * gxf
      dvdy = dvdy +  n1.v * gyf

      local n2 = nrow[pcy+2]
      gxf = gxi * py[2]
      gyf = pxi * gy[2]
      dudx = dudx +  n2.u * gxf
      dudy = dudy +  n2.u * gyf
      dvdx = dvdx +  n2.v * gxf
      dvdy = dvdy +  n2.v * gyf
    end

    local w1  = dudy - dvdx
    local wT0 = w1 * p.T01
    local wT1 = 0.5 * w1 * (p.T00 - p.T11)
    local D00 = dudx
    local D01 = 0.5 * (dudy + dvdx)
    local D11 = dvdy
    local trace = 0.5 * (D00 + D11)
    D00 = D00-trace
    D11 = D11-trace

    p.T00 = p.T00 + (-wT0 + D00) - yieldRate * p.T00
    p.T01 = p.T01 + ( wT1 + D01) - yieldRate * p.T01
    p.T11 = p.T11 + ( wT0 + D11) - yieldRate * p.T11
    local norma = p.T00 * p.T00 + 2.0 * p.T01 * p.T01 + p.T11 * p.T11
    if norma > 5.0 then
       p.T00 = 0.0
       p.T01 = 0.0
       p.T11 = 0.0
    end

    local cx = math.abs(math.floor(p.x))
   -- if cx>self.gsizeX then cx=5 end
    local cy = math.abs(math.floor(p.y))
  --  if cy>self.gsizeY then cy=5 end
  --  cx=50
 --   cy=50
   -- test3=self.gsizeX.." / "..self.gsizeY.." / "..cx.." / "..cy
    
   --test="uups "..cx.." - "..cy.." timer:"..self.timer.." steps:"..self.steps
    
    local cxi = cx + 1
    local cyi = cy + 1
    
    local n00  = self.grid[cx][cy]
    local n01  = self.grid[cx][cyi]
    local n10  = self.grid[cxi][cy]
    local n11  = self.grid[cxi][cyi]

    local p00 = n00.m
    local x00 = n00.gx
    local y00 = n00.gy
    local p01 = n01.m
    local x01 = n01.gx
    local y01 = n01.gy
    local p10 = n10.m
    local x10 = n10.gx
    local y10 = n10.gy
    local p11 = n11.m
    local x11 = n11.gx
    local y11 = n11.gy

    local pdx = p10 - p00
    local pdy = p01 - p00
    local C20 =  3.0 * pdx - x10 - 2.0 * x00
    local C02 =  3.0 * pdy - y01 - 2.0 * y00
    local C30 = -2.0 * pdx + x10 + x00
    local C03 = -2.0 * pdy + y01 + y00
    local csum1 = p00 + y00 + C02 + C03
    local csum2 = p00 + x00 + C20 + C30
    local C21 =   3.0 * p11 - 2.0 * x01 - x11 - 3.0 * csum1  - C20
    local C31 = (-2.0 * p11 +        x01 + x11 + 2.0 * csum1) - C30
    local C12 =   3.0 * p11 - 2.0 * y10 - y11 - 3.0 * csum2  - C02
    local C13 = (-2.0 * p11 +        y10 + y11 + 2.0 * csum2) - C03
    local C11 = x01 - C13 - C12 - x00

    local u  = p.x - cx*1.00 --float cx
    local u2 = u * u
    local u3 = u * u2
    local v  = p.y - cy*1.00 -- float cy
    local v2 = v * v
    local v3 = v * v2
    local density = p00 + x00 * u + y00 * v + C20 * u2 + C02 * v2 + C30 * u3 + C03 * v3 + C21 * u2 * v  + C31 * u3 * v + C12 * u  * v2 + C13 * u  * v3 + C11 * u  * v

    local DS = densitySetting
    --[[ spatially varying density field:
    local DS = densitySetting * ( math.pow(  (p.y / self.gsizeY), 3.0) )
    ]]--
    
    --if not DS then DS=2.5 end
    
    local pressure = (stiffness / math.max(1.0, DS)) * (density - DS)

    if pressure > 2.0 then
      pressure = 2.0
    end
    
    local fx = 0.0
    local fy = 0.0
    if p.x < 3.0 then
      fx = fx +  3.0 - p.x
    elseif p.x > (self.gsizeX - 4) then
      fx = fx +  (self.gsizeX - 4) - p.x
    end

    if p.y < 3 then
      fy = fy +  3 - p.y
    elseif p.y > (self.gsizeY - 4) then
      fy = fy +  (self.gsizeY - 4) - p.y
    end

    trace = trace * self.settings[1]
    local T00 = elasticity * p.T00 + viscosity * D00 + pressure + bulkViscosity * trace
    local T01 = elasticity * p.T01 + viscosity * D01
    local T11 = elasticity * p.T11 + viscosity * D11 + pressure + bulkViscosity * trace

    local dx=nil
    local dy=nil
    for i=0,3 do
      local nrow = self.grid[pcx + i]
      local ppxi = px[i]
      local pgxi = gx[i]

      local n0 = nrow[pcy]
      phi = ppxi * py[0]
      dx  = pgxi * py[0]
      dy  = ppxi * gy[0]
      n0.ax = n0.ax +  fx * phi -(dx * T00 + dy * T01)
      n0.ay = n0.ay + fy * phi -(dx * T01 + dy * T11)

      local n1 = nrow[pcy + 1]
      phi = ppxi * py[1]
      dx  = pgxi * py[1]
      dy  = ppxi * gy[1]
      n1.ax = n1.ax + fx * phi -(dx * T00 + dy * T01)
      n1.ay = n1.ay + fy * phi -(dx * T01 + dy * T11)

      local n2 = nrow[pcy + 2]
      phi = ppxi * py[2]
      dx  = pgxi * py[2]
      dy  = ppxi * gy[2]
      n2.ax = n2.ax + fx * phi -(dx * T00 + dy * T01)
      n2.ay = n2.ay + fy * phi -(dx * T01 + dy * T11)
    end
  end


  for ni=0,self.nActiveNodes do
    local n = self.activeNodeArray[ni]
    if n and (n.m > 0.0) then
      n.ax = n.ax/n.m
      n.ay = n.ay/n.m
      n.u = 0.0
      n.v = 0.0
    end
  end
--test=self.settings[6]

[littfass]

– simulate function part 3

  -- Particles pass 3
  local gravity =self.settings[6]
  for ip=0,self.nParticles do
    local p = self.particles[ip]

    local px = p.px -- with []
    local py = p.py -- with []
    local pcy = p.cy
    local pcx = p.cx
    for i=0,3 do
      local nrow = self.grid[pcx + i]
      local ppxi = px[i]

      local n0 = nrow[pcy]
      phi = ppxi * py[0]
      p.u = p.u +  phi * n0.ax
      p.v = p.v +  phi * n0.ay

      local n1 = nrow[pcy + 1]
      phi = ppxi * py[1]
      p.u = p.u + phi * n1.ax
      p.v = p.v + phi * n1.ay

      local n2 = nrow[pcy + 2]
      phi = ppxi * py[2]
      p.u = p.u + phi * n2.ax
      p.v = p.v + phi * n2.ay
    end

    p.v = p.v + gravity
    if drag then
      local vx = abs(p.x - mx)
      local vy = abs(p.y - my)
      if (vx < 10) and (vy < 10) then
        local weight = (1.0 - vx / 10) * (1.0 - vy / 10)
        p.u = p.u + weight * (mdx - p.u)
        p.v = p.v + weight * (mdy - p.v)
      end
    end

    local xf = p.x + p.u
    local yf = p.y + p.v
    if xf < 2.0 then
      p.u = p.u + (2.0 - xf) + math.random()*0.01
    elseif (xf > (self.gsizeX - 3.0)) then
      p.u = p.u + (self.gsizeX - 3.0) - xf - math.random()*0.01
    end
    if (yf < 2.0) then
      p.v = p.v + (2.0 - yf) + math.random()*0.01
    elseif (yf > (self.gsizeY - 3.0)) then
      p.v = p.v + (self.gsizeY - 3.0) - yf - math.random()*0.01
    end


    local pu = p.u
    local pv = p.v
    for i=0,3 do
      local nrow = self.grid[pcx + i]
      local ppxi = px[i]

      local n0 = nrow[pcy]
      phi = ppxi * py[0]
      n0.u = n0.u + phi * pu
      n0.v = n0.v + phi * pv

      local n1 = nrow[pcy + 1]
      phi = ppxi * py[1]
      n1.u = n1.u + phi * pu
      n1.v = n1.v + phi * pv

      local n2 = nrow[pcy + 2]
      phi = ppxi * py[2]
      n2.u = n2.u + phi * pu
      n2.v = n2.v + phi * pv
    end
  end

  for ni=0,self.nActiveNodes do
    local n = self.activeNodeArray[ni]
    if n and (n.m > 0.0) then
      n.u = n.u/n.m
      n.v = n.v/n.m
    end
  end

  --Particles pass 4
  local gu=nil
  local gv=nil
  local smoothing = 1.00--self.settings[7]
  for ip=0,self.nParticles do
    local p = self.particles[ip]


    gu = 0.0
    gv = 0.0

    local px = p.px
    local py = p.py
    local pcy = p.cy
    local pcx = p.cx
    for i=0,3 do
      local nrow = self.grid[pcx + i]
      local ppxi = px[i]

      local n0 = nrow[pcy]
      phi = ppxi * py[0]
      gu = gu + phi * n0.u
      gv = gv + phi * n0.v

      local n1 = nrow[pcy + 1]
      phi = ppxi * py[1]
      gu = gu + phi * n1.u
      gv = gv + phi * n1.v

      local n2 = nrow[pcy + 2]
      phi = ppxi * py[2]
      gu = gu + phi * n2.u
      gv = gv + phi * n2.v
    end
    

    
    p.gu = gu
    p.gv = gv
    p.x = p.x + gu
    p.y = p.y + gv
    
    p.u = p.u + smoothing * (gu - p.u)
    p.v = p.v + smoothing * (gv - p.v)
    self.steps = self.steps + 1
--    test2=p.x.." - "..p.y.." steps:"..self.steps
    
 --       if p.x<0 then p.x=0 end
  --  if p.x>WIDTH/self.SC-4 then p.x=WIDTH/self.SC-5 end
 --   if p.y<0 then p.y=0 end
 --   if p.y>HEIGHT/self.SC-4 then p.y=HEIGHT/self.SC-4 end
  end
  
end

[littfass]

– The particle and node classes

Particle = class()

function Particle:init(inx, iny, inu, inv)
    if not inx then inx=0 end
    if not iny then iny=0 end
    if not inu then inu=0 end
    if not inv then inv=0 end
    
    self.x = inx
    self.y = iny
    self.u = inu
    self.v = inv
    
    self.cx = 0
    self.cy = 0
    
    self.gu=0
    self.gv=0
    self.T00=0
    self.T01=0
    self.T11=0

    self.px={0,0,0}
    self.py={0,0,0}
    self.gx={0,0,0}
    self.gy={0,0,0}
end



Node = class()

function Node:init()
    self.m  = 0
    self.d  = 0
    self.gx = 0
    self.gy = 0
    self.u  = 0
    self.v  = 0
    self.ax = 0
    self.ay = 0
    self.active = false
end

function Node:clear()
    self.m = 0
    self.d = 0
    self.gx = 0
    self.gy = 0
    self.u = 0
    self.v = 0
    self.ax = 0
    self.ay = 0
    self.active = false
end

[littfass]

–and finally the main class

function setup()
   displayMode(FULLSCREEN)
    liquid=MPM()
end

function draw()
    liquid:draw()
end

@connorbot, thanks, the 3 tilde kept code formatting

[Connorbot999]

Is not running the error is

error: error: [string "MPM = class()..."]:86: attempt to call global 'map' (a nil value)

Pausing playback

[littfass]

@connorbot
You will have to include the map() and millis() functions from above

Notice, even then, this is a broken implementation. However, changing the settings makes the simulation crash more or less faster

[littfass]

@ruilov
Here is a quick and dirty test with lots of tiny bodies (there is a weird ?bug? though after few seconds, breaking the world bounds)
With 100 bodies it’s ok, with 200 it becomes slow, with 300 it fails

function setup()
    maxitems=100
    drops={}  
    supportedOrientations(LANDSCAPE_LEFT)
    displayMode(FULLSCREEN)
    local x=WIDTH/2
    local y=HEIGHT/2
    for i=0,maxitems do
        table.insert(drops,addbody(x+i,y)) 
    end
    boundingbox()
end

function boundingbox()
    ground = physics.body(POLYGON, vec2(0,1), vec2(0,0), vec2(WIDTH,0), vec2(WIDTH,1))
    ground.type = STATIC
    groundleft = physics.body(POLYGON, vec2(1,0), vec2(0,0), vec2(0,HEIGHT), vec2(1,HEIGHT))
    groundleft.type = STATIC
    groundright = physics.body(POLYGON, vec2(WIDTH - 1,0), vec2(WIDTH,0), vec2(WIDTH,HEIGHT), vec2(WIDTH - 1,HEIGHT))
    groundright.type = STATIC
    ceiling = physics.body(POLYGON, vec2(0,HEIGHT-1), vec2(0,HEIGHT), vec2(WIDTH,HEIGHT), vec2(WIDTH,HEIGHT-1))
     ceiling.type = STATIC
end

function addbody(x,y)
    local r=5
    local circle = physics.body(CIRCLE, r)
    circle.gravityScale = 10+math.random()
    circle.interpolate = true
    circle.x = x
    circle.y = y
    circle.restitution = 0.7
    circle.sleepingAllowed = false   
    return circle
end

function draw()
    local g=vec3(Gravity.x,Gravity.y,Gravity.z)
    physics.gravity(g)
    background(40, 40, 50)
    noStroke()
    for i,body in ipairs(drops) do
        fill(0, 200, 255, 126)
        ellipse(body.x,body.y,5)
    end
end

[littfass]

@Bri_G
Many thanks for the bug fix, updated the above test code.
Maybe the garbage collector is the issue in the MPM class, I’ll have a look later

[Bri_G]

Hi @littfass,

Just tried out your small tester, floor falls away after several cycles. Resolved that by removing local definition on your floor, left, right and ceiling. By that I mean delete the preceding word ‘local’ on each line. Then it runs fine.

Hope that helps.

Bri_G