Fluid Simulation

Game development specific discussions.
Boromir
Posts: 463
Joined: Apr 30, 2015 19:28
Location: Oklahoma,U.S., Earth,Solar System
Contact:

Fluid Simulation

Post by Boromir »

Simple Physics fluid simulation for a platformer I'm working on.
Press "E" to render blue water color.
Move mouse on screen to pull points.

Code: Select all

const k = 0.025
type wpoint
    y as single
    velocity as single
end type
dim points(100) as wpoint
for i as integer=0 to ubound(points)
    points(i).y=240
next

screen 18,32
dim as integer i

do
screenlock
cls

for i=0 to ubound(points)
if i<>ubound(points) then line (i*10,points(i).y)-(i*10+10,points(i+1).y),rgb(0,0,255)
if not multikey(18) then line (i*10-3,points(i).y-3)-(i*10+3,points(i).y+3),rgb(255,0,255),bf
next
if multikey(18) then paint (1,479),rgb(0,0,255),rgb(0,0,255)
screenunlock
    for i=0 to ubound(points)
    dim as single x
    x = 240-(points(i).y)
    points(i).velocity += k * x-0.025*points(i).velocity
    points(i).y += points(i).velocity
next

dim as single l(ubound(points))
dim as single r(ubound(points))

for j as integer = 0 to 8
    for i = 0 to ubound(points)
        if i > 0 then
            l(i) = 0.025 * (points(i).y - points(i - 1).y)
            points(i-1).velocity+=l(i)
        end if
        if i < ubound(points) then
            r(i)=0.025 * (points(i).y - points(i + 1).y)
            points(i+1).velocity += r(i)
        end if
    next
 
    for i = 0 to ubound(points)
        if i > 0 then points(i - 1).y += l(i)
        if i < ubound(points) then points(i + 1).y += r(i)
    next
next

dim as integer x,y
getmouse x,y
for i=0 to ubound(points)
    dim as integer dist,x2,y2,test
    x2=i*10
    y2=points(i).y

    '            test = SQR(abs((x-x2)*(x-x2))+abs((y-y2)*(y-y2)))
    'test=x-x2
    'if y<y2 then test=0
    'if test>=0 and test<=50 then test=50-test
    if x>0 and x>x2-50 and x<=x2 then points(i).y+=x-x2+50
    if x>0 and x>x2 and x<x2+50 then points(i).y-=x-x2-50

next

sleep 50,1
loop until multikey(1)

MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: Fluid Simulation

Post by MrSwiss »

@Boromir,

you might like it, now in fast (optimized for speed):

Code: Select all

' FluidSim.bas -- original by Boromir, optimized for speed by MrSwiss

const As Single k = 0.025
Const           maxPts = 100    ' Integer
Const As ULong  blue = RGB(0,0,255), pink = RGB(255,0,255)

type wpoint
    as Single y, velocity
end Type

dim as wpoint points(maxPts)

' init Type
for i as UInteger = 0 to maxPts
    With points(i)
        .y = 479 / 2
        .velocity = 0.0
    End With
next

' ===== MAIN =====
screen 18, 32

dim as integer  i, j, x, y
Dim as single   l(maxPts), r(maxPts)

Do
    screenlock
    Cls
    for i = 0 to maxPts
        Var x2 = i * 10
        If i <> maxPts then line (x2, points(i).y)-(x2 + 10, points(i + 1).y), blue
        If not multikey(18) then line (x2 - 3, points(i).y - 3)-(x2 + 3, points(i).y + 3), pink, bf
    Next
    if multikey(18) then paint (0,479), blue, blue
    Locate 2, 2, 0 : Print "press [Enter] to EXIT!";
    ScreenUnLock
    
    For i = 0 to maxPts
        With points(i)
            Var x1 = 240 - .y
            .velocity += k * x1 - k * .velocity
            .y += .velocity
        End with
    Next
    
    for j = 0 to 7
        for i = 0 to maxPts
            if i > 0 then
                l(i) = k * (points(i).y - points(i - 1).y)
                points(i - 1).velocity += l(i)
            end if
            if i < maxPts then
                r(i) = k * (points(i).y - points(i + 1).y)
                points(i + 1).velocity += r(i)
            end If
            if i > 0 then points(i - 1).y += l(i)
            if i < maxPts then points(i + 1).y += r(i)
        next
    next
    
    getmouse x, y
    
    for i = 0 to maxPts
        With points(i)
            Var x2 = i * 10, y2 = .y
            If x > 0 AndAlso x > x2 - 50 AndAlso x <= x2 then .y += x - x2 + 50
            if x > 0 AndAlso x > x2 AndAlso x < x2 + 50 then .y -= x - x2 - 50
        End With
    Next
    
    sleep 50,1
loop until MultiKey(&h1C)   ' press [Enter] to quit
' ===== END-MAIN =====
' ----- EOF -----
Boromir
Posts: 463
Joined: Apr 30, 2015 19:28
Location: Oklahoma,U.S., Earth,Solar System
Contact:

Re: Fluid Simulation

Post by Boromir »

Nice. I think I learned a few things going over what you did.;)

Stress test with 1280 points

Code: Select all

' FluidSim.bas -- original by Boromir, optimized for speed by MrSwiss
dim as any ptr water,temp
const As Single k = 0.025
const as single damp = 0.145
Const           maxPts = 1280    ' Integer
Const As ULong  blue = RGB(0,100,255), pink = RGB(255,0,255)
const pdist=1
const scrw=1280
const scrh=1024

type wpoint
    as Single y, velocity
end Type

dim as wpoint points(maxPts)

' init Type
for i as UInteger = 0 to maxPts
    With points(i)
        .y = 479 / 2
        .velocity = 0.0
    End With
next

' ===== MAIN =====
screenres scrw,scrh, 32
water=imagecreate(scrw,scrh)
temp=imagecreate(scrw,scrh)
for i as uinteger=1 to scrh
    line water,(0,i)-(scrw,i),rgb(0,0,(scrh-i+100)/5)
next
dim as integer  i, j, x, y
Dim as single   l(maxPts), r(maxPts)

Do
    screenlock
    Cls
    if multikey(18) then put temp,(0,0),water,alpha
    for i = 0 to maxPts
        Var x2 = i * pdist
        If i <> maxPts then line temp,(x2, points(i).y)-(x2 + pdist, points(i + 1).y), blue
        If not multikey(18) then line (x2 - 3, points(i).y - 3)-(x2 + 3, points(i).y + 3), pink, bf
    Next
    if multikey(18) then paint temp,(0,0), blue, blue
    if multikey(18) then put(0,0),temp,alpha
    Locate 2, 2, 0 : Print "press [ESC] to EXIT!";
    ScreenUnLock
   
    For i = 0 to maxPts
        With points(i)
            Var x1 = 240 - .y
            .velocity += k * x1 - damp * .velocity
            .y += .velocity
        End with
    Next
   
    for j = 0 to 7
        for i = 0 to maxPts
            if i > 0 then
                l(i) = k * (points(i).y - points(i - 1).y)
                points(i - 1).velocity += l(i)
            end if
            if i < maxPts then
                r(i) = k * (points(i).y - points(i + 1).y)
                points(i + 1).velocity += r(i)
            end If
            if i > 0 then points(i - 1).y += l(i)
            if i < maxPts then points(i + 1).y += r(i)
        next
    next
   
    getmouse x, y
   
    for i = 0 to maxPts
        With points(i)
            Var x2 = i * pdist, y2 = .y
            If x > 0 AndAlso x > x2 - 50 AndAlso x <= x2 then .y += x - x2 + 50
            if x > 0 AndAlso x > x2 AndAlso x < x2 + 50 then .y -= x - x2 - 50
        End With
    Next
   
    sleep 50,1
loop until MultiKey(1)   ' press [ESC] to quit
imagedestroy(water)
' ===== END-MAIN =====
' ----- EOF -----
Last edited by Boromir on Mar 10, 2017 15:13, edited 1 time in total.
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: Fluid Simulation

Post by MrSwiss »

Boromir, you are not really testing the same thing ... (comparing Apples with Bananas, so to speak).
Because: using Images & Put, is ruining the speed (Put is slow).
Boromir
Posts: 463
Joined: Apr 30, 2015 19:28
Location: Oklahoma,U.S., Earth,Solar System
Contact:

Re: Fluid Simulation

Post by Boromir »

MrSwiss wrote:Boromir, you are not really testing the same thing ... (comparing Apples with Bananas, so to speak).
Because: using Images & Put, is ruining the speed (Put is slow).
You have to push "e" to render. I was just testing if there was a decrease in operation speed with 1280 points without doing rendering.
The implementation of "put" was for a better visual experience when it comes time to merging into my platformer.
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: Fluid Simulation

Post by MrSwiss »

Boromir wrote:I was just testing if there was a decrease in operation speed with 1280 points without doing rendering.
Not quite true - because, besides increasing the Window's size (array's), you've at the same time,
increased the number of calculated points by a factor of 10!

Code: Select all

const pdist=1    ' was 10 before
This naturally decreases the speed considerably ... (you can't have it, and eat it, at the same time).

You'll have to find a acceptable balance, for your specific needs.
Boromir
Posts: 463
Joined: Apr 30, 2015 19:28
Location: Oklahoma,U.S., Earth,Solar System
Contact:

Re: Fluid Simulation

Post by Boromir »

The change in "pdist" was so I could see all the points on the screen at once. I don't see how it effects performance?
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: Fluid Simulation

Post by MrSwiss »

Unmodified code can be used up to a ScreenRes(1000, 480, 32) because,
with a spacing of 10, we're covering the whole distance with 101 Points.
(we might have to leave some of them outside of view-port, for quality)

Keep in mind that maxPts is sizing 3 arrays ( array of Type, l() and r() ).
The way I understand the code:
  • arrays l() and r() are for smoothing, from point to point.
    Therefore, calculating every point (without spacing) makes no sense.
It's even more than factor 10, it's more like 12.8 ^ 3 (3 array's) ~ factor 20.
(array's (of Type) grow exponential (by members), that gets sometimes forgotten)
Boromir
Posts: 463
Joined: Apr 30, 2015 19:28
Location: Oklahoma,U.S., Earth,Solar System
Contact:

Re: Fluid Simulation

Post by Boromir »

MrSwiss wrote:Unmodified code can be used up to a ScreenRes(1000, 480, 32) because,
with a spacing of 10, we're covering the whole distance with 101 Points.
(we might have to leave some of them outside of view-port, for quality)

Keep in mind that maxPts is sizing 3 arrays ( array of Type, l() and r() ).
The way I understand the code:
  • arrays l() and r() are for smoothing, from point to point.
    Therefore, calculating every point (without spacing) makes no sense.
It's even more than factor 10, it's more like 12.8 ^ 3 (3 array's) ~ factor 20.
(array's (of Type) grow exponential (by members), that gets sometimes forgotten)
Yes, but in the end it is maxPts that is causing the performance decrease.
Changing pdist alone doesn't have any impact on performance.
Iv'e found setting pdist to 10 has the best effect, but shortening increases the amount of points that can fit on the screen.
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: Fluid Simulation

Post by MrSwiss »

Iv'e found setting pdist to 10 has the best effect
Yep, correct.
but shortening increases the amount of points that can fit on the screen.
The amount of points should more or less match: points * pdist = scr_width (some more left/right maybe).

But what I'm trying to get across is:
  • it's better to increase the array's only slightly (to cover 1280) - and -
    leave the spacing at 10
I'd try, maxPts = 128 + 16 (left/right 8 outside view-port) ... maxPts = 144.

[edit]code updated, 2017-03-11[/edit]

Code: Select all

' FluidSim.bas -- original by Boromir, optimized for speed by MrSwiss

Const As Single k = 0.025
Const           maxPts = 144, pdist = 10    ' Integer
' colors
Const As ULong  blue = &hFF0000FF, pink = &hFFFF00FF, br_blue = &hFF007FFF, white = &hFFFFFFFF
' screen
Const wid = 1280, hei = 800, hhei = hei Shr 1   ' hhei = hei \ 2

Type wpoint
    As Single y, velocity
End Type

Dim As wpoint points(maxPts)

' init Type
For i As UInteger = 0 To maxPts
    With points(i)
        .y = hhei
        .velocity = 0.0
    End With
Next

' ===== MAIN =====
ScreenRes(wid, hei, 32)

Dim As Integer  i, j, x, y
Dim As Single   l(maxPts), r(maxPts)

Do
    ScreenLock
    Cls
    For i = 0 To maxPts
        Var x2 = i * pdist
        If i <> maxPts Then Line (x2, points(i).y)-(x2 + 10, points(i + 1).y), blue
        If not MultiKey(&h12) Then Line (x2 - 3, points(i).y - 3)-(x2 + 3, points(i).y + 3), pink, bf
    Next
    If MultiKey(18) Then Paint (0, hei - 1), blue, blue : Paint (0, 0), br_blue, blue
    Draw String (8, 8), "press [Enter] to EXIT!", white
    ScreenUnLock

    For i = 0 To maxPts
        With points(i)
            Var x1 = hhei - .y
            .velocity += k * x1 - k * .velocity
            .y += .velocity
        End with
    Next

    For j = 0 To pdist - 1
        For i = 0 To maxPts
            If i > 0 Then
                l(i) = k * (points(i).y - points(i - 1).y)
                points(i - 1).velocity += l(i)
            End If
            If i < maxPts Then
                r(i) = k * (points(i).y - points(i + 1).y)
                points(i + 1).velocity += r(i)
            End If
            If i > 0 Then points(i - 1).y += l(i)
            If i < maxPts Then points(i + 1).y += r(i)
        Next
    Next

    GetMouse(x, y)

    For i = 0 To maxPts
        With points(i)
            Var x2 = i * pdist, y2 = .y
            If x > 0 AndAlso x > x2 - 50 AndAlso x <= x2 Then .y += x - x2 + 50
            If x > 0 AndAlso x > x2 AndAlso x < x2 + 50 Then .y -= x - x2 - 50
        End With
    Next
   
    Sleep(50, 1)
Loop Until MultiKey(&h1C)   ' press [Enter] to quit
' ===== END-MAIN ===== ' ----- EOF -----
Boromir
Posts: 463
Joined: Apr 30, 2015 19:28
Location: Oklahoma,U.S., Earth,Solar System
Contact:

Re: Fluid Simulation

Post by Boromir »

I noticed you used andalso instead of and. I've seen this done before but I don't understand how it is faster. When should I use andalso instead of and.
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: Fluid Simulation

Post by MrSwiss »

I noticed you used andalso instead of and. I've seen this done before but I don't understand how it is faster.
Binary Operators are in FB:
  • Not (unary, also called inverter, bit level)
    And, Or, Xor etc.
    all work on a bit-wise level (also, every instance of statement is always evaluated)
Boolean Operators are in FB (short cut capability):
  • AndAlso
    OrElse
    boolean evaluation (TRUE/FALSE), the first *wrong* statement finishes evaluation immediately.
This is the reason, they are also called *short cut Operators*.
NOTE: new since FBC ver. 1.04.0 (together with Boolean variables introduced)
Since in most cases boolean evaluation is working, that is: simple TRUE/FALSE checks, they are faster
than binary evaluation, that also always checks all given conditonals.
When should I use andalso instead of and.
Binary AND is mainly used (since AndAlso) to *mask* a value, e.g. x And &h01 (isolates the lowest bit
of x for a later 1/0 check).
AndAlso can be used for most other tasks (where traditionally, And was used).
h4tt3n
Posts: 698
Joined: Oct 22, 2005 21:12
Location: Denmark

Re: Fluid Simulation

Post by h4tt3n »

This thread reminded me of an old project I had lying around...

This body of water has a constant surface area. The 4th order Runge-Kutta integrator is totally overkill and can be replaced with simple symplectic Euler.

Cheers, Mike

Code: Select all

'******************************************************************************'
'
'   Soft body "water" simulation. Michael "h4tt3n" Nissen, march 2008
'   (Press esc to quit)
'
'******************************************************************************'

Const Pi                =     4*Atn(1)  ''  pi (better not change ^^)
Const dt                =     0.01      ''  timestep, delta time
Const half_dt           =     dt/2      ''  one half timestep (for integration)
Const sixth_dt          =     dt/6      ''  one sixth timestep (for integration)
Const Num_Blobs         =     1         ''  number of blobs
Const Num_Masses        =     96        ''  number of vertices in blob
Const Point_Mass        =     2         ''  mass of each vertice
Const Kp                =     0.5       ''  pressure constant
const area_factor       =     1.0       ''  change rest area to in- or deflate
Const Ks                =     2000      ''  linear spring stiffnes
Const Kld               =     10        ''  linear spring damping
Const Kps               =     2000      ''  mouse-body spring stiffnes
Const Kpd               =     2         ''  mouse-body spring damping
Const Grav_Acc          =     0         ''  gravitational acceleration
Const Scrn_Wid          =     1200      ''  screen width
Const Scrn_Hgt          =     500       ''  screen height
Const Pick_mindist      =     30        ''  pick up mass within distance
Const Prod_Radius       =     50        ''  prod ball radius
Const water_level			=		scrn_hgt * 0.75

Randomize Timer

''  define types
Type Vector_Type
  As Double X, Y
End Type

Type Mass_Type
  As Double Mass
  As Vector_Type Frc, Acc, Vel, Pos
  As Vector_Type cur_vel, cur_pos
  As Vector_Type K1_vel, K1_Pos, K2_Vel, K2_Pos, K3_Vel, K3_Pos, K4_Vel, K4_Pos
End Type 

Type Spring_Type
  As Integer a, b
  As Double Length, Rest_Length
  As Vector_Type Lng_Hat, Vel
End Type

Type Blob_Type
  As Mass_Type Mass(1 to Num_Masses)
  As Spring_Type Spring(1 to Num_Masses)
  As Uinteger col
  As Double Total_Mass, Area, Rest_area, Pressure, Kp
  As Vector_Type Frc, Acc, Vel, Pos
End Type

Declare Sub Init_Blobs()
Declare Sub Integrate_RK4()
Declare Sub Accelerate_Masses()
Declare Sub Collision_Detection()
Declare Sub Mouse_Interaction()
Declare Sub Draw_Blobs()

Dim Shared as Blob_Type Blob(1 to Num_Blobs)
Dim Shared As Vector_Type Dist, Vel, closest_point
Dim Shared As Double Force, Damping, Distance, temp_dist, frc_a, frc_b, _
  sin_angle, cos_angle, torque, dmp_a, dmp_b, blob_area
Dim Shared As Integer i, j, k, l, mouse_x, mouse_y, mouse_btn, pick_state, Picked_Mass, Picked_Blob


'*******************************************************************************


Init_Blobs()

ScreenRes Scrn_Wid, Scrn_hgt, 16
Color RGB(0, 0, 0), RGB(255, 255, 255)

''  main program loop
Do
  

   Collision_Detection()
   Integrate_RK4()
   Draw_Blobs()
  	Sleep 1, 1
  
Loop Until Multikey(1)

End


'*******************************************************************************

Sub Draw_Blobs()
  
  ''  draw to screen
  ScreenLock
	  Cls 
	  For i = Lbound(Blob) to Ubound(Blob)
	    With Blob(i)
	      For j = Lbound(.Spring) To Ubound(.Spring)
	        With .Spring(j)
	          Line _
	            (Blob(i).Mass(.a).Pos.X, Blob(i).Mass(.a).Pos.Y)-_
	            (Blob(i).Mass(.b).Pos.X, Blob(i).Mass(.b).Pos.Y), 0
	        End With
	      Next
	      Paint (1, scrn_hgt-1), .col, 0
	      If Pick_State = 1 Then
	        Line _
	        (Mouse_x, Mouse_Y)-_
	        (Blob(Picked_Blob).Mass(Picked_Mass).Pos.X, Blob(Picked_Blob).Mass(Picked_Mass).Pos.Y), _
	        Rgb(0, 0, 0)
	      End If
	      If Mouse_Btn = 2 Then
	        Circle(Mouse_X, Mouse_Y), Prod_Radius, 0,,,1, f
	      End If
	    End With
	  Next
	  
	  Locate 1, 1: Print Using "Area (pixels)  : ######.#"; Blob_Area
	  Locate 2, 1: Print Using "Area (percent) :    ###.###"; (blob(1).rest_area/blob_area)*100
	  
  ScreenUnlock
  
End Sub

Sub Init_Blobs()
  
  ''  initiate blobs
  For i = Lbound(Blob) To Ubound(Blob)
    With Blob(i)
      .Kp = Kp
      .col = RGB(64, 128, 255)
      .Pos.X = 0
      .Pos.Y = 0
      
     ''	set mass positions
     With .Mass(1)
       .pos.x = scrn_wid
       .pos.y = water_level
     End With
     With .Mass(2)
       .pos.x = scrn_wid
       .pos.y = scrn_hgt
     End With
     With .Mass(3)
       .pos.x = 0
       .pos.y = scrn_hgt
     End With
     With .Mass(4)
       .pos.x = 0
       .pos.y = water_level
     End With
     
     For j = 5 To Ubound(.Mass)
        With .Mass(j)
        	.pos.x = (scrn_wid/(ubound(blob(i).mass)-3))*(j-4)
        	.pos.y = water_level'+(Rnd-Rnd)*12
        End With
     Next
     
     ''	find sum of masses
     For j = Lbound(.Mass) To Ubound(.Mass)
        With .Mass(j)
          .Mass = Point_Mass
          Blob(i).Total_Mass += .Mass
        End With
      next
      
      ''	set springs
      For j = 1 To Ubound(.Mass)
        With .Spring(j)
          k = j+1
          If k > Num_Masses Then k -= Num_Masses
          .a = j
          .b = k
          .Rest_Length = scrn_wid/(Ubound(blob(i).Mass))-6
        End With
      Next
      
      
      ''  ideal area or "rest area"
      For j = Lbound(.Mass) To Ubound(.Mass)
        With .Spring(j)
          Blob(i).Rest_Area += Blob(i).Mass(.a).Pos.X*Blob(i).Mass(.b).Pos.Y-_
                               Blob(i).Mass(.a).Pos.Y*Blob(i).Mass(.b).Pos.X
        End With
      Next
      .Rest_Area *= 0.5
      .Rest_Area *= area_factor
      
    End With
  Next
  
End Sub

Sub Mouse_Interaction()
  
  Getmouse Mouse_X, Mouse_Y,, Mouse_Btn
  
  ''  pick up blob with leftmouse
  If Mouse_Btn = 1 Then
    If Pick_State = 0 Then
      Temp_Dist = Pick_MinDist
      For i = Lbound(Blob) To Ubound(Blob)
        With Blob(i)
          For j = Lbound(.mass) To Ubound(.mass) 
            With .Mass(j)
              Dist.X = .Pos.X-Mouse_X
              Dist.Y = .Pos.Y-Mouse_Y
              Distance = Sqr(Dist.X*Dist.X+Dist.Y*Dist.Y)
              If Distance < Temp_Dist Then 
                Temp_Dist = Distance
                Picked_Blob = i
                Picked_Mass = j
              End If
            End With
          Next
        End With
      Next
      If Picked_Blob <> -1 Then
        Pick_State = 1
      End If
    Else
      With Blob(Picked_Blob).Mass(Picked_Mass)
        Dist.X = .Pos.X-Mouse_X
        Dist.Y = .Pos.Y-Mouse_Y
        Distance = Sqr(Dist.X*Dist.X+Dist.Y*Dist.Y)
        Force = -(Kps*Distance)-Kpd*((Dist.X*.Vel.X+Dist.Y*.Vel.Y)/Distance)
        .Frc.X += Force*(Dist.X/Distance)
        .Frc.Y += Force*(Dist.Y/Distance)
      End With
    End If
  Else
    Pick_State = 0
    Picked_Mass = -1
    Picked_Blob = -1
  End If
  
  ''  prod blobs with right mouse
  If Mouse_Btn = 2 Then
    Pick_State = 0
    Picked_Mass = -1
    Picked_Blob = -1
    For i = Lbound(Blob) To Ubound(Blob) 
      With Blob(i)
        For j = Lbound(.Mass) To Ubound(.Mass)
          Dist.X = .mass(j).Pos.X-Mouse_X
          Dist.Y = .mass(j).Pos.Y-Mouse_Y
          Distance = Sqr(Dist.X*Dist.X+Dist.Y*Dist.Y)
          If Distance < Prod_Radius Then 
            With .Mass(j)
              Force = -16000*(Distance-Prod_Radius)-4*((Dist.X*.Vel.X+Dist.Y*.Vel.Y)/Distance)
              .Frc.X += Force*(Dist.X/Distance)
              .Frc.Y += Force*(Dist.Y/Distance)
            End With
          End If
        Next
      End With
    Next
  End If
End Sub

Sub Accelerate_Masses()
  
  For i = Lbound(Blob) To Ubound(Blob)
    With Blob(i)
      
      ''  blob (center of mass) position and velocity
      .Pos.X = 0
      .Pos.y = 0
      For j = Lbound(.Mass) to Ubound(.Mass)
        With .Mass(j)
          Blob(i).Pos.X += .Mass*.Pos.X
          Blob(i).Pos.Y += .Mass*.Pos.Y
        End With
      Next
      .Pos.X /= .Total_Mass
      .Pos.y /= .Total_Mass
      
      ''  apply spring force
      For j = 4 to Ubound(.Spring)
        With .Spring(j)
          Dist.X = Blob(i).Mass(.b).Pos.X-Blob(i).Mass(.a).Pos.X
          Dist.Y = Blob(i).Mass(.b).Pos.Y-Blob(i).Mass(.a).Pos.Y
          .Length = Sqr(Dist.X*Dist.X+Dist.Y*Dist.Y)
          .Lng_Hat.X = Dist.X/.Length
          .Lng_Hat.Y = Dist.Y/.Length
          .Vel.X = Blob(i).Mass(.b).Vel.X-Blob(i).Mass(.a).Vel.X
          .Vel.Y = Blob(i).Mass(.b).Vel.Y-Blob(i).Mass(.a).Vel.Y
          Force = -Ks*(.Length-.Rest_Length)-Kld*(.Lng_Hat.X*.Vel.X+.Lng_Hat.Y*.Vel.Y)
          Blob(i).Mass(.a).Frc.X -= Force*.Lng_Hat.X
          Blob(i).Mass(.a).Frc.Y -= Force*.Lng_Hat.Y
          Blob(i).Mass(.b).Frc.X += Force*.Lng_Hat.X
          Blob(i).Mass(.b).Frc.Y += Force*.Lng_Hat.Y
        End With
      Next
      
      ''  blob area
      Blob_Area = .Area
      .Area = 0
      For j = Lbound(.Spring) To Ubound(.Spring)
        With .Spring(j)
          Blob(i).Area += Blob(i).Mass(.a).Pos.X*Blob(i).Mass(.b).Pos.Y-_
                          Blob(i).Mass(.b).Pos.X*Blob(i).Mass(.a).Pos.Y
        End With
      Next
      .Area *= 0.5
      
      ''  blob pressure
      .Pressure = -.Kp*(.Area-.Rest_Area)
      
      ''  apply pressure force
      For j = 4 To Ubound(.Spring)
        With .Spring(j)
          Force = .Length*Blob(i).Pressure
          Blob(i).Mass(.a).Frc.X += Force*.Lng_Hat.Y
          Blob(i).Mass(.a).Frc.Y += Force*-.Lng_Hat.X
          Blob(i).Mass(.b).Frc.X += Force*.Lng_Hat.Y
          Blob(i).Mass(.b).Frc.Y += Force*-.Lng_Hat.X
        End With
      Next
      
    End With
  Next
  
  ''	accelerate masses
  For i = Lbound(Blob) To Ubound(Blob)
    With Blob(i)
      For j = 1 to 4
        With .Mass(j)
          .Acc.Y = .Frc.Y/.mass
          '.Acc.Y += Grav_Acc
          .frc.y = 0
        End With
      Next
    End With
  Next
  For i = Lbound(Blob) To Ubound(Blob)
    With Blob(i)
      For j = 5 to Ubound(.Mass)
        With .Mass(j)
          .Acc.X = .Frc.X/.mass
          .Acc.Y = .Frc.Y/.mass
          .Acc.Y += Grav_Acc
          .frc.x = 0
          .frc.y = 0
        End With
      Next
    End With
  Next
  
End Sub

Sub Collision_Detection()
  
  ''  blob - screen edge collision detection
  For i = Lbound(Blob) To Ubound(Blob)
    With Blob(i)
      For j = 5 to Ubound(.Mass)
        With .Mass(j)
          ''  right
          If .Pos.x > Scrn_Wid-1 Then
            .Frc.x -= 64000*(.Pos.x-(Scrn_wid-1))
          End If
          ''  left
          If .Pos.x < 1 Then
            .Frc.x -= 64000*(.Pos.x-1)
          End If
          ''  bottom
          If .Pos.Y > Scrn_Hgt-1 Then
            .Frc.y -= 64000*(.Pos.Y-(Scrn_Hgt-1))
          End If
          ''  top
          If .Pos.Y < 1 Then
            .Frc.y -= 64000*(.Pos.Y-1)
          End If
        End With
      Next
    End With
  Next
  
End Sub


Sub Integrate_RK4()
  
  Mouse_Interaction()
  Accelerate_Masses()
  
  ''  save current position and velocity
  ''  K1 step
  For i = Lbound(Blob) To Ubound(Blob)
    With Blob(i)
      For j = Lbound(.Mass) to Ubound(.Mass)
        With .Mass(j)
          .cur_pos.x = .pos.x
          .cur_pos.y = .pos.y
          .cur_vel.x = .vel.x
          .cur_vel.y = .vel.y
          .k1_pos.x = .vel.x
          .k1_pos.y = .vel.y
          .k1_vel.x = .acc.x
          .k1_vel.y = .acc.y
          .pos.x = .cur_pos.x + .k1_pos.x*half_dt
          .pos.y = .cur_pos.y + .k1_pos.y*half_dt
          .vel.x = .cur_vel.x + .k1_vel.x*half_dt
          .vel.y = .cur_vel.y + .k1_vel.y*half_dt
        End With
      Next
    End With
  Next
  
  Mouse_Interaction()
  Accelerate_Masses()
  
  ''  K2 step
  For i = Lbound(Blob) To Ubound(Blob)
    With Blob(i)
      For j = Lbound(.Mass) to Ubound(.Mass)
        With .Mass(j)
          .k2_pos.x = .vel.x
          .k2_pos.y = .vel.y
          .k2_vel.x = .acc.x
          .k2_vel.y = .acc.y
          .pos.x = .cur_pos.x + .k2_pos.x*half_dt
          .pos.y = .cur_pos.y + .k2_pos.y*half_dt
          .vel.x = .cur_vel.x + .k2_vel.x*half_dt
          .vel.y = .cur_vel.y + .k2_vel.y*half_dt
        End With
      Next
    End With
  Next
  
  Mouse_Interaction()
  Accelerate_Masses()
  
  ''  K3 step
  For i = Lbound(Blob) To Ubound(Blob)
    With Blob(i)
      For j = Lbound(.Mass) to Ubound(.Mass)
        With .Mass(j)
          .k3_pos.x = .vel.x
          .k3_pos.y = .vel.y
          .k3_vel.x = .acc.x
          .k3_vel.y = .acc.y
          .pos.x = .cur_pos.x + .k3_pos.x*dt
          .pos.y = .cur_pos.y + .k3_pos.y*dt
          .vel.x = .cur_vel.x + .k3_vel.x*dt
          .vel.y = .cur_vel.y + .k3_vel.y*dt
        End With
      Next
    End With
  Next
  
  Mouse_Interaction()
  Accelerate_Masses()
  
  ''  K4 step
  ''  Find new position and velocity
  For i = Lbound(Blob) To Ubound(Blob)
    With Blob(i)
      For j = Lbound(.Mass) to Ubound(.Mass)
        With .Mass(j)
          .k4_pos.x = .vel.x
          .k4_pos.y = .vel.y
          .k4_vel.x = .acc.x
          .k4_vel.y = .acc.y
          .pos.x = .cur_pos.x + (.k1_pos.x + (.k2_pos.x + .k3_pos.x)*2 + .k4_pos.x) * sixth_dt
          .pos.y = .cur_pos.y + (.k1_pos.y + (.k2_pos.y + .k3_pos.y)*2 + .k4_pos.y) * sixth_dt
          .vel.x = .cur_vel.x + (.k1_vel.x + (.k2_vel.x + .k3_vel.x)*2 + .k4_vel.x) * sixth_dt
          .vel.y = .cur_vel.y + (.k1_vel.y + (.k2_vel.y + .k3_vel.y)*2 + .k4_vel.y) * sixth_dt
        End With
      Next
    End With
  Next

End Sub
Boromir
Posts: 463
Joined: Apr 30, 2015 19:28
Location: Oklahoma,U.S., Earth,Solar System
Contact:

Re: Fluid Simulation

Post by Boromir »

That's one nice bit of physics. ;D
I used Euler integration in my water demo. I am new to anything beyond basic physics so I was looking to make a very simple and understandable bit of code to introduce myself into this awesome new world.

My program has 3 tweakable factors to obtain different liquid effects. Jello,molasses,water,etc...
k=spring stiffness factor
damp=dampening factor
spread=wave spreading factor

Code: Select all

' FluidSim.bas -- original by Ezekiel Gutierrez, optimized for speed by MrSwiss

dim shared as integer mx,my,click
Const As ULong  blue = RGBA(0,100,255,0), pink = RGB(255,0,255)

const scrw=1280
const scrh=1024

type wpoint
    as Single y, velocity
end Type

type liquid
    dim as any ptr water,temp
    as Single k = 0.025
    as single damp = 0.025
    as single spread=0.075
    maxpts as integer
    pdist  as integer
    reDim as single l(1)
    redim as single r(1)
    redim points(1) as wpoint
    declare sub init(var1 as integer,var2 as integer)
    declare sub move
    declare sub drawr
end type

screenres scrw,scrh, 32,,1
' init
sub liquid.init(var1 as integer,var2 as integer)
maxpts=var1
pdist=var2
reDim as single l(maxPts), r(maxPts)
redim points(maxpts) as wpoint

for i as UInteger = 0 to maxPts
    With points(i)
        .y = 140
        .velocity = 0.0
    End With
next
water=imagecreate(scrw,scrh)
temp=imagecreate(scrw,scrh)
for i as uinteger=1 to scrh
    dim as integer a=(i+100)/2
    if a>230 then a=230
    line water,(0,i)-(scrw,i),rgba(0,0,(scrh-i+200)/5,a)
next
end sub

' ===== MAIN =====
sub liquid.drawr
    if not multikey(18) then put temp,(0,0),water,alpha
    for i as integer = 0 to maxPts
        Var x2 = i * pdist
        If i <> maxPts then line temp,(x2, points(i).y)-(x2 + pdist, points(i + 1).y), blue
        If multikey(18) then line (x2 - 3, points(i).y - 3+600)-(x2 + 3, points(i).y + 3+600), pink, bf
    Next
    if not multikey(18) then paint temp,(0,0), blue, blue
    if not multikey(18) then put(0,0+600),temp,(0,0)-(scrw,424),alpha
    dim as any ptr mask=imagecreate(1280,1024,rgba(0,0,0,0))
    if 0+600+1024<1024 then line mask,(0,0+600+1024)-(scrw,scrh),rgba(0,0,190/5,230),bf:put (0,0),mask,alpha
    imagedestroy(mask)
end sub

sub liquid.move
    For i as integer = 0 to maxPts
        With points(i)
            Var x1 = 140 - .y
            .velocity += k * x1 - damp * .velocity
            .y += .velocity
        End with
    Next
   
    for j as integer = 0 to 7
        for i as integer= 0 to maxPts
            if i > 0 then
                l(i) = spread * (points(i).y - points(i - 1).y)
                points(i - 1).velocity += l(i)
            end if
            if i < maxPts then
                r(i) = spread * (points(i).y - points(i + 1).y)
                points(i + 1).velocity += r(i)
            end If
            if i > 0 then points(i - 1).y += l(i)
            if i < maxPts then points(i + 1).y += r(i)
        next
    next
   
    getmouse mx,my,, click
    if click=1 then
    for i as integer = 0 to maxPts
        With points(i)
            Var x2 = i * pdist, y2 = .y
            If mx > 0 AndAlso mx > x2 - 50 AndAlso mx <= x2 then .y += mx - x2 + 50
            if mx > 0 AndAlso mx > x2 AndAlso mx < x2 + 50 then .y -= mx - x2 - 50
        End With
    Next
    end if
end sub


dim as liquid water1
water1.init(130,10)


do
    water1.move
    
    
    screenlock
    cls
    water1.drawr
    screenunlock
    
    sleep 50,1
loop until multikey(1)
' ===== END-MAIN =====


imagedestroy(water1.water)
imagedestroy(water1.temp)


MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: Fluid Simulation

Post by MrSwiss »

ReDim in the used Type is not needed (use *(Any)* for dynamic array's):

Code: Select all

Type liquid
    Dim As Any Ptr water, temp
    As Integer maxpts, pdist
    As Single  k = 0.025, damp = 0.025, spread = 0.075
    As Single  l(Any), r(Any)
    As wpoint  points(Any)
    Declare Sub init(ByVal var1 As Integer, ByVal var2 As Integer)
    Declare Sub move()	' indicate: no parameters
    Declare Sub drawr()
End Type
Post Reply