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!

May 25th, 2012 at 11:35

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

October 9th, 2012 at 03:01

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

October 9th, 2012 at 13:51

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

October 9th, 2012 at 15:44

That will be great! My email is wensi.isnew@gmail.com. Would you also explain the uses of grid.gx/grid.gy and how they are computed? Thank you!