Edited main post

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!

Hi @lttlfas,

Ignore my last comment found the processing code.

Bri_G

:)

@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

--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

-- 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]

-- 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

Is not running the error is

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

Pausing playback

@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

Hi All,

I suspect this ma be a bug, but possibly only with physics. The local definitions should not be available to any of the other code , yet it runs for several cycles.

Bri_G

:)

Hi @littfass, I'm a newbie with Codea/Lua but I like a lot your project, because I want to use Codea for trying to simulate real world problems (if possible ) :-D With regard to your "quick test", I can't explain (for the moment) the reason why it slows down when increasing the particle number, but I think one reason could be the way the fluid particles are "created": at the screen-center and, due to interactions between them, they displace each other simultanously while they are created. I'm not sure if you've fixed that already, but I think, an improvement in the "creation" procedure could help... but I don't know it exactly :/ Greetings from Mexico,

Victor

@littfass, just a tiny remark: the "creation" procedure in your "quick test" creates 101 particles (instead of, e.g., 100); I mean, your initial index should be 1 instead of 0 in order to create your maximum number of items.

Edit: I've trying to figure out what could be wrong when playing with more than 200 particles, focusing for the moment on your quick test and with my scarce knowledge of Codea/Lua. It seems to me that the "creation" function isn't the main poblem, which begins to appear when tilting the ipad... It might be that the problem is the drawing ability of the device (computation-speed required for the drawing loop in your code for every frame)... :-? ~O) Anyway, I attach here a short video from a similar simulation that I made: the "flow" of perfectly elastic "fluid particles" inside a box, "applying" a momentum in x-direction to them in order to start with an initial velocity. After that, they interact with a "solid body" (e.g. a cylinder) and the walls of the box.

I'll attach my code here later, after "cleaning it", for the use of beginners like me. Saludos desde México, Victor ~O)

If you are intersted, maybe you can do operations with the image, a wave ripple must be easy thought, you only have to change the colors info of an image using set and get with an image data structure, Check this simple JavaScript http://badassjs.com/post/932217000/water-ripple-effect-in-canvas