RealBasic: isolated project 3 -> Liquid physics


This is a very rudimentary physics simulation of liquid water. It cannot show al lot of particles in plain RealBasic. I suspect if you rewrite it for GDI+ or even OpenGL you could use thousands of particles.

It was based on the work of Grant Kot (who made several demos for HTML5) and jgittins/quinbd’s version for Android.

It’s Grant Kots implementation of the Material Point Method (MPM). Wiki: http://en.wikipedia.org/wiki/Material_Point_Method

Grant Kot:

For interpolation, I use the quadratic B-spline presented here: Analysis and Reduction of Quadrature Errors in the Material Point Method. Please note that there is an error in equation 17, for the cubic B-spline. The middle two equations should end with 2/3. I don’t use the cubic spline though, because I’m interested in real-time simulation and I feel that the quadratic spline is the best balance of speed and quality.

Instead of integrating the density over time (which is what most of the MPM papers do), I do a density summation every frame. Because this is not dependent on previously calculated values of density, there is no accumulated error. To minimize grid artifacts, I use the cubic interpolation method presented here: A Semi-Lagrangian CIP Fluid Solver without Dimensional Splitting.

I had to read it also several times before it made any sense. 🙂

If you want to get deeper into the mechanism read:
Analysis and Reduction of Quadrature Errors in the Material Point Method: http://www.gorgeousapps.com/steffen_08_IJNME_mpm_preprint.pdf

and
A Semi-Lagrangian CIP Fluid Solver without Dimensional Splitting: http://www.gorgeousapps.com/PA-08-09-19.pdf

Anyhow, the result of it was this piece of code I translated to RealBasic:

  Dim drag As Boolean = False
  Dim mdx As Single
  Dim mdy As Single
  
  Dim i, j as Integer
  Dim cxi As Integer
  Dim cyj As Integer
  Dim cyi as Integer
  Dim n As ABNode
  
  Dim phi As Single
  Dim dx, dy As Single
  Dim x, y As Single
  
  Dim cx As Integer
  Dim cy As Integer
  
  Dim p00 As Single
  Dim x00 As Single
  Dim y00 As Single
  Dim p01 As Single
  Dim x01 As Single 
  Dim y01 As Single 
  Dim p10 As Single 
  Dim x10 As Single 
  Dim y10 As Single 
  Dim p11 As Single 
  Dim x11 As Single 
  Dim y11 As Single 
  
  Dim pdx As Single 
  Dim pdy As Single 
  Dim C20 As Single 
  Dim C02 As Single 
  Dim C30 As Single 
  Dim C03 As Single 
  Dim csum1 As Single 
  Dim csum2 As Single 
  Dim C21 As Single 
  Dim C31 As Single 
  Dim C12 As Single 
  Dim C13 As Single 
  Dim C11 As Single 
  
  Dim u As Single 
  Dim u2 As Single 
  Dim u3 As Single 
  Dim v As Single 
  Dim v2 As Single 
  Dim v3 As Single 
  Dim density As Single 
  
  Dim gx,gy as Single
  Dim mu, mv As Single 
  Dim gu, gv As Single 
  
  Dim pressure As Single 
  
  Dim fx As Single
  Dim fy As Single
  
  Dim p as ABParticle
  
  Dim vx As Single
  Dim vy As Single
  Dim weight As Single
  
  If (pressed) And (pressedprev) Then
    drag = True
    mdx = 0.25 * (mx - mxprev)
    mdy = 0.25 * (my - myprev)
  End If
  
  pressedprev = pressed
  mxprev = mx
  myprev = my
  
  For Each n  In active
    n.m = 0
    n.d = 0
    n.gx = 0
    n.gy = 0
    n.u = 0
    n.v = 0
    n.ax = 0
    n.ay = 0
    n.active = False
  Next
  Redim active(-1)
  
  For Each p  In Particles
    p.cx = p.x - 0.5
    p.cy = p.y - 0.5
    
    x = p.cx - p.x
    p.px(0) = (0.5 * x * x + 1.5 * x + 1.125)
    p.gx(0) = (x + 1.5)
    x = x + 1.0
    p.px(1) = (-x * x + 0.75)
    p.gx(1) = (-2.0 * x)
    x = x + 1.0
    p.px(2) = (0.5 * x * x - 1.5 * x + 1.125)
    p.gx(2) = (x - 1.5)
    
    y = p.cy - p.y
    p.py(0) = (0.5 * y * y + 1.5 * y + 1.125)
    p.gy(0) = (y + 1.5)
    y = y + 1.0
    p.py(1) = (-y * y + 0.75)
    p.gy(1) = (-2.0 * y)
    y = y + 1.0
    p.py(2) = (0.5 * y * y - 1.5 * y + 1.125)
    p.gy(2) = (y - 1.5)
    
    For i  = 0 To 2
      For j  = 0 To 2
        cxi = p.cx + i
        cyj = p.cy + j
        n = grid(cxi,cyj)
        If Not n.active Then
          n.active = True
          active.Append n
        End If
        phi = p.px(i) * p.py(j)
        n.m = n.m + phi * p.mat.m
        n.d = n.d + phi
        dx = p.gx(i) * p.py(j)
        dy = p.px(i) * p.gy(j)
        n.gx = n.gx + dx
        n.gy = n.gy + dy
      Next
    Next
  Next
  
  For Each p In Particles
    cx = p.x
    cy = p.y
    cxi  = cx + 1
    cyi  = cy + 1
    
    p00 =grid(cx,cy).d
    x00 =grid(cx,cy).gx
    y00 =grid(cx,cy).gy
    p01 =grid(cx,cyi).d
    x01 =grid(cx,cyi).gx
    y01 =grid(cx,cyi).gy
    p10 =grid(cxi,cy).d
    x10 =grid(cxi,cy).gx
    y10 =grid(cxi,cy).gy
    p11 =grid(cxi,cyi).d
    x11 =grid(cxi,cyi).gx
    y11 =grid(cxi,cyi).gy
    
    pdx = p10 - p00
    pdy = p01 - p00
    C20 = 3.0 * pdx - x10 - 2.0 * x00
    C02 = 3.0 * pdy - y01 - 2.0 * y00
    C30 = -2.0 * pdx + x10 + x00
    C03 = -2.0 * pdy + y01 + y00
    csum1 = p00 + y00 + C02 + C03
    csum2 = p00 + x00 + C20 + C30
    C21 = 3.0 * p11 - 2.0 * x01 - x11 - 3.0 * csum1 - C20
    C31 = -2.0 * p11 + x01 + x11 + 2.0 * csum1 - C30
    C12 = 3.0 * p11 - 2.0 * y10 - y11 - 3.0 * csum2 - C02
    C13 = -2.0 * p11 + y10 + y11 + 2.0 * csum2 - C03
    C11 = x01 - C13 - C12 - x00
    
    u = p.x - cx
    u2 = u * u
    u3 = u * u2
    v = p.y - cy
    v2 = v * v
    v3 = v * v2
    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
    
    pressure  = density - 1.0
    If pressure > 2.0 Then
      pressure = 2.0
    End If
    
    fx = 0.0
    fy = 0.0
    
    If p.x < 4.0 Then
      fx = fx + p.mat.m * (4.0 - p.x)
    ElseIf p.x > gsizeX - 5 Then
      fx = fx + p.mat.m * (gsizeX - 5 - p.x)
    End If
    If p.y < 4.0 Then
      fy = fy + p.mat.m * (4.0 - p.y)
    ElseIf p.y > gsizeY - 5 Then
      fy = fy + p.mat.m * (gsizeY - 5 - p.y)
    End If
    
    If drag Then
      vx = Abs(p.x - 0.25 * mx)
      vy = Abs(p.y - 0.25 * my)
      If (vx < 10.0) And (vy < 10.0) Then
        weight = p.mat.m * (1.0 - vx / 10.0) * (1.0 - vy / 10.0)
        fx = fx + weight * (mdx - p.u)
        fy = fy + weight * (mdy - p.v)
      End If
    End If
    
    For i  = 0 To 2
      For j  = 0 To 2
        n = grid((p.cx + i),(p.cy + j))
        phi  = p.px(i) * p.py(j)
        gx  = p.gx(i) * p.py(j)
        gy  = p.px(i) * p.gy(j)
        
        n.ax = n.ax + -(gx * pressure) + fx * phi
        n.ay = n.ay + -(gy * pressure) + fy * phi
      Next
    Next
  Next
  
  For Each n  In Active
    If n.m > 0.0 Then
      n.ax = n.ax / n.m
      n.ay = n.ay / n.m
      n.ay = n.ay + 0.03
    End If
  Next
  For Each p in Particles
    For i  = 0 To 2
      For j  = 0 To 2
        n = grid((p.cx + i),(p.cy + j))
        phi = p.px(i) * p.py(j)
        p.u = p.u + phi * n.ax
        p.v = p.v + phi * n.ay
      Next
    Next
    mu = p.mat.m * p.u
    mv = p.mat.m * p.v
    For i  = 0 To 2
      For j  = 0 To 2
        n = grid((p.cx + i),(p.cy + j))
        phi  = p.px(i) * p.py(j)
        n.u = n.u + phi * mu
        n.v = n.v + phi * mv
      Next
    Next
  Next
  For Each n  In Active
    If n.m > 0.0 Then
      n.u = n.u / n.m
      n.v = n.v / n.m
    End If
  Next
  For Each p  In Particles
    gu = 0.0
    gv = 0.0
    For i  = 0 To 2
      For j  = 0 To 2
        n = grid((p.cx + i),(p.cy + j))
        phi = p.px(i) * p.py(j)
        gu = gu + phi * n.u
        gv = gv + phi * n.v
      Next
    Next
    p.x = p.x + gu
    p.y = p.y + gv
    p.u = p.u + 1.0 * (gu - p.u)
    p.v = p.v + 1.0 * (gv - p.v)
    If p.x < 1.0 Then
      p.x = (1.0 + rnd() * 0.01)
      p.u = 0.0
    ElseIf p.x > gsizeX - 2 Then
      p.x = (gsizeX - 2 - rnd() * 0.01)
      p.u = 0.0
    End If
    If p.y < 1.0 Then
      p.y = (1.0 + rnd() * 0.01)
      p.v = 0.0
    ElseIf p.y > gsizeY - 2 Then
      p.y = (gsizeY - 2 - rnd() * 0.01)
      p.v = 0.0
    End If
  Next 

The full source code and demo can be downloaded from: http://www.gorgeousapps.com/ABLiquid.zip

Until next time!

Advertisements

About Alwaysbusy

My name is Alain Bailleul and I'm the Senior Software Architect/Engineer at One-Two. I like to experiment with new technologies, Computer Vision and A.I. My projects are programmed in B4X , Xojo, C#, java, HTML, CSS and JavaScript. View all posts by Alwaysbusy

4 responses to “RealBasic: isolated project 3 -> Liquid physics

  • Alwaysbusy

    Tip: Use your mouse to play with the ‘water’ 🙂

  • Vincent Wen

    Hey Alain,

    I’m trying to do something similar in C. I’m using Stam’s Stable Fluid solver for the fluid dynamics. I wanted to use the USCIP method for advection but can’t get it working. Would you please explain how you did it? I haven’t programmed in Basic before so I can’t fully understand your code.

    Thank you,
    Vincent

    • Alwaysbusy

      Hi Vincent,

      This is a rather old project so I’m not very familiar with the code anymore. I’ve looked a little into the work of Stam’s Stable Fluid solver and it looks very interesting! I think it goes a little further than my old project. The code in the article is also only the simulate function which is only part of the full system. I’ll have a look if I can find the original java version and if you want I’ll email it to you. Let me know.

      Cheers,
      Alain

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: