particle based interactive fluid

[h4tt3n]

Hello folks,

Some time ago I posted a simple smoothed particle hydrodynamics example, and recently I've made it somewhat better. This code sample uses grid based collision detection and handles up to around 10.000 particles depending on your system. Ideas, suggestions and critique are welcome :-)

Cheers,
Mike

Code: Select all

''	rough smoothed particle hydrodynamics (SPH) draft
''	Michael "h4tt3n" Nissen, july 5th 2010
''	
''	todo:
''	predefined cell neighbors
''	better kernel (less bouncy liquid)
''	better broad phase collision detection
''	better mouse interaction
''	better code structure

''	includes
#Include Once "GL/gl.bi"

''	constants
Const as single		dt						= 0.01											''	timestep
Const as single		half_dt				= 0.5*dt										''	half timestep
const as single		r 						= 10												''	interaction radius
const as single		r_2 					= 2*r												''	two radius (for grid)
const as single		r_sqrd 				= r*r												''	radius squared
Const as single		k			 				= 400												''	global "spring stiffnes"
Const as single		k_near			 	= 3000											''	local "spring stiffnes"
Const as single		viscosity			= 2													''	viscosity
Const as single		rest_dens 		= 6													''	rest density
const as single		grav_acc			= 400												''	gravity
Const as single		pi						= 4*atn(1)									''	pi
Const as single		radius				= 2													''	particle distance on startup
const as integer  block_width		= 60												''	<--- change to reduce / increase num masses
const as integer  block_height	= 100												''	
Const as integer  num_masses		= block_width*block_height	''	total number of masses
const as integer  neighbors			= 48												''	max number of neighbors
Const as integer 	screen_wid 		= 800												''	screen width
Const as integer 	screen_hgt 		= 400												''	screen height
Const As Single 	wall_friction	= 1													''	wall friction (1 = full, 0 = off)
Const As Single 	interaction_Radius	= 8										''	mouse interaction radius

Const As Integer cell_dia						= r
Const As Integer cell_max_masses		= 24
const as integer num_cells_row			= screen_wid/cell_dia
const As integer num_cells_col			= screen_hgt/cell_dia

Const As Integer num_cells					= num_cells_row*num_cells_col

Const As Single cell_wid						= screen_wid/num_Cells_row
Const As single cell_hgt						= screen_hgt/num_Cells_Col
const As single half_cell_wid				= cell_wid/2

''	types
Type CellType
	As Integer Num_Masses, Mass(1 To cell_max_masses)
End Type

type vectortype
	as single x, y
end Type

Type NeighborType
	As Integer a
	As Single Dist
	As vectortype norm_dist
End Type

type masstype
	As Integer num_neighbors
	As NeighborType Neighbor(1 To neighbors)
	as vectortype frc, vel, psn
	as Single density, angvel, torque
end Type

''	dim variables
Dim as integer 	center_x 	= screen_wid\2
Dim as integer 	center_y	= screen_hgt\2

dim as integer mouse_x, mouse_y, mouse_x_old, mouse_y_old, mouse_vel_x, mouse_vel_y

dim Shared As masstype mass(1 to num_masses)
dim Shared As celltype cell(1 to num_cells_row, 1 To num_cells_col)

''	narrow phase collision detection. Make particles neighbors if distance < r
Sub CheckIfNeighbor(ByVal a_ As Integer, ByVal b_ As Integer)
	
	If mass(a_).num_neighbors = neighbors Then Exit sub
	If mass(b_).num_neighbors = neighbors Then Exit Sub
	
	Dim As vectortype dist
	
	dist.x = mass(b_).psn.x-mass(a_).psn.x: If Abs(dist.x) > r Then Exit Sub
	dist.y = mass(b_).psn.y-mass(a_).psn.y: If Abs(dist.y) > r Then Exit Sub
	
	Dim As Single dist_sqrd = dist.x*dist.x+dist.y*dist.y
	
	If dist_sqrd > r_sqrd Or dist_sqrd = 0 Then Exit Sub
	
	Dim As Single distance = Sqr(dist_sqrd)
	
	Dim As vectortype norm_dist = (dist.x/distance, dist.y/distance)
	
	Dim As Single density = 1 - distance / r
	
	'' create neighbors and save data
	With mass(a_)
		.num_neighbors += 1
		.Neighbor(.num_neighbors).a = b_
		.Neighbor(.num_neighbors).dist = distance
		.Neighbor(.num_neighbors).norm_dist.x = norm_dist.x
		.Neighbor(.num_neighbors).norm_dist.y = norm_dist.y
		.Density += density
	End With
	
	With mass(b_)
		.num_neighbors += 1
		.Neighbor(.num_neighbors).a = a_
		.Neighbor(.num_neighbors).dist = distance
		.Neighbor(.num_neighbors).norm_dist.x = -norm_dist.x
		.Neighbor(.num_neighbors).norm_dist.y = -norm_dist.y
		.Density += density
	End With
	
End Sub

Randomize Timer

Dim As Integer n = 1

''	place particles
For i As Integer = 1 To block_width
	For j As Integer = 1 To block_height
		With mass(n)
			.psn.x = i * radius * 2
			.psn.y = screen_hgt - j * radius * 2
			
			.psn.x += (Rnd-Rnd)
			.psn.y += (Rnd-Rnd)
			
			.vel.x += (Rnd-Rnd)
			.vel.y += (Rnd-Rnd)

		End With
		n += 1
	Next
Next

''	create opengl screen
ScreenRes screen_wid, screen_hgt, 32,, 2
glMatrixMode(GL_PROJECTION)  
glLoadIdentity()
glViewport(0, 0, screen_wid, screen_hgt)
glOrtho(0, screen_wid, screen_hgt, 0, 0, 32)
glMatrixMode(GL_MODELVIEW)
glLoadIdentity()
glClearColor (0.6, 0.7, 1.0, 0.0)
'glenable(GL_ALPHA)
glEnable(GL_BLEND)
glblendfunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
'glblendfunc(GL_DST_COLOR , GL_SRC_COLOR)

getmouse(mouse_x, mouse_y)

''	main loop
Do
	
	'' save old mouse state
	mouse_x_old = mouse_x
	mouse_y_old = mouse_y
	
	''	get new mouse state
	getmouse(mouse_x, mouse_y)
	
	'' calc mouse velocity
	mouse_vel_x = (mouse_x-mouse_x_old)/dt
	mouse_vel_y = (mouse_y-mouse_y_old)/dt
	
	''	mouse interaction
	for i as integer = 1 to num_masses
		with mass(i)	
			
			Dim As vectortype dst
			
			dst.x = .psn.x - mouse_x:	if dst.x > interaction_radius then continue for
			dst.y = .psn.y - mouse_y:	if dst.y > interaction_radius then continue for
			
			if dst.x*dst.x+dst.y*dst.y > interaction_Radius*interaction_Radius then continue for
			
			if mouse_vel_x then .vel.x += (mouse_vel_x-.vel.x) * 0.16
			if mouse_vel_y then .vel.y += (mouse_vel_y-.vel.y) * 0.16
			
		end with	
	Next	
	
	''	reset particles
	for i as integer = 1 to num_masses
		With mass(i)
			.density = 0
			.num_neighbors = 0
		End With
	Next
	
	''	reset cells
	For i As Integer = 1 To num_cells_row
		For j As Integer = 1 To num_cells_col
			cell(i, j).Num_Masses = 0
		Next
	Next
	
	''	update grid (broad phase collision detection)
	For i as integer = 1 to num_masses
		Dim As Integer cell_row = (mass(i).Psn.y\cell_hgt)+1
		Dim as integer cell_col = (mass(i).Psn.x\cell_wid)+1
		With cell(cell_col, cell_row)
			If .num_masses = cell_max_masses Then Continue For
			.num_masses += 1
			.Mass(.num_masses) = i
		End With
	Next
	
	''	narrow phase collision detection
	For i As Integer = 1 To num_cells_row
		For j As Integer = 1 To num_cells_col
			If cell(i, j).num_masses = 0 Then Continue For
			For k_ As Integer = 1 To cell(i, j).num_masses
				Dim As Integer a_ = cell(i, j).mass(k_)
				
				If mass(a_).num_neighbors = neighbors Then Continue For
				
				''	cell(x, y) - self 
				For l_ As Integer = k_+1 To cell(i, j).num_masses
					CheckIfNeighbor(a_, cell(i, j).mass(l_))
				Next
				
				''	cell(x, y) - cell(x+1, y)
				If i < num_cells_row Then
					For l_ As Integer = 1 To cell(i+1, j).num_masses
						CheckIfNeighbor(a_, cell(i+1, j).mass(l_))
					Next
				End If
				
				If j = num_cells_col Then Continue For
				
				''	cell(x, y) - cell(x, y+1)
				For l_ As Integer = 1 To cell(i, j+1).num_masses
					CheckIfNeighbor(a_, cell(i, j+1).mass(l_))
				Next
				
				''	cell(x, y) - cell(x+1, y+1)
				If i < num_cells_row Then
					For l_ As Integer = 1 To cell(i+1, j+1).num_masses
						CheckIfNeighbor(a_, cell(i+1, j+1).mass(l_))
					Next
				End If
				
				''	cell(x, y) - cell(x-1, y+1)
				If i > 1 Then
					For l_ As Integer = 1 To cell(i-1, j+1).num_masses
						CheckIfNeighbor(a_, cell(i-1, j+1).mass(l_))
					Next
				End If
				
			Next
		Next
	Next
		
	''	set pressure
	For i as integer = 1 to num_masses
		
		If mass(i).num_neighbors = 0 Then Continue For

		''	mass i global pressure
		Dim As Single pressure = -k * (mass(i).density - rest_dens)
		
		For j as integer = 1 to mass(i).num_neighbors
			
			Dim As Single gradient = 1 - mass(i).neighbor(j).dist / r
				
			Dim As Integer l = mass(i).neighbor(j).a
			
			''	mass i + l global pressure
			Dim As Single press = (pressure - k * (mass(l).density - rest_dens) )
			
			''	mass i + l local pressure
			press -= k_near * gradient
			
			''	viscosity
			Dim As vectortype vel = (mass(i).vel.x-mass(l).vel.x, mass(i).vel.y-mass(l).vel.y)
				
			Dim As Single proj_vel = vel.x*mass(i).neighbor(j).norm_dist.x+vel.y*mass(i).neighbor(j).norm_dist.y
			
			press -=  viscosity * proj_vel * gradient
			
			press *= gradient
				
			mass(i).frc.x += press * mass(i).neighbor(j).norm_dist.x
			mass(i).frc.y += press * mass(i).neighbor(j).norm_dist.y
			
			mass(l).frc.x -= press * mass(i).neighbor(j).norm_dist.x
			mass(l).frc.y -= press * mass(i).neighbor(j).norm_dist.y
			
		Next
		
	Next
	
	''	integrate (since particle mass = 1, force and acceleration are identical)
	for i as integer = 1 to num_masses
		with mass(i)
			
			.frc.y += grav_acc
			
			.vel.x += .frc.x * half_dt
			.vel.y += .frc.y * half_dt
			
			.angvel += .torque * half_dt
			
			.psn.x += .vel.x * dt + .frc.x * dt * half_dt
			.psn.y += .vel.y * dt + .frc.y * dt * half_dt
			
			.vel.x += .frc.x * half_dt
			.vel.y += .frc.y * half_dt
			
			.angvel += .torque * half_dt
			
			.frc.x = 0
			.frc.y = 0
			
			.torque = 0
			
		end with
	Next
	
	''	keep particles inside screen
	For i as integer = 1 to num_masses
		with mass(i)
			
			if .psn.x < radius then .psn.x = radius: .vel.x = -.vel.x*wall_friction: .vel.y *= wall_friction
			if .psn.x > screen_wid - radius then .psn.x = screen_wid -radius: .vel.x = -.vel.x*wall_friction: .vel.y *= wall_friction
			
			If .psn.y < radius then .psn.y = radius: .vel.y = -.vel.y*wall_friction: .vel.x *= wall_friction
			If .psn.y > screen_hgt - radius Then .psn.y = screen_hgt -radius: .vel.y = -.vel.y*wall_friction: .vel.x *= wall_friction
			
		end with
	Next
	
	''	draw to screen
	glClear(GL_DEPTH_BUFFER_BIT Or GL_COLOR_BUFFER_BIT)
	
	glLoadIdentity()
	
	glenable(GL_POINT_SMOOTH)
	glpointsize(r)
	glcolor4f(0.0f, 0.0f, 0.0f, 0.25f)
	
	GLbegin(GL_POINTS)
		for i as integer = 1 to num_masses
			glvertex2f(mass(i).psn.x, mass(i).psn.y)
		Next
	GLend()
	
	gldisable(GL_POINT_SMOOTH)
	glpointsize(1)
	glcolor4f(1.0f, 1.0f, 0.0f, 0.75f)
	
	glbegin(GL_POINTS)
		for i as integer = 1 to num_masses
			glvertex2f(mass(i).psn.x, mass(i).psn.y)
		Next
	GLend()
	
	GLflush()
	
	Flip()
	
	sleep 1, 1
	
Loop Until MultiKey(1)

end

[rolliebollocks]

Cool.

[pestery]

Cool, you can see the shock waves moving through the fluid. Ha, it goes berserk when you move the mouse onto the screen from beneath the fluid. The code is really good too, pretty well labeled and easy to follow, although I had to read through the viscosity and integration bits a couple of times before I got it.

How hard would it be to make it 3d? You could model the flow across a height-map maybe. Could be useful for simulating a dam breaking and washing away stuff in a game for example.

[h4tt3n]

... I had to read through the viscosity and integration bits a couple of times before I got it.

How hard would it be to make it 3d? You could model the flow across a height-map maybe. Could be useful for simulating a dam breaking and washing away stuff in a game for example.

Viscosity is basically a velocity based spring damping (see my tutorial for details!) scaled by distance. The integration algorithm is just plane velocity verlet, which imho is the best.

Making this 3D would be fairly simple, since it's based entirely on linear forces between particles. All you need to do is add the z component to all vectors. Opengl will do all the messy graphics for you.

Feel free to tamper, experiment and re-post funny discoveries!

Cheers,
Mike

[rolliebollocks]

ceiling circles!

Code: Select all

''        rough smoothed particle hydrodynamics (SPH) draft
''        Michael "h4tt3n" Nissen, july 5th 2010
''        
''        todo:
''        predefined cell neighbors
''        better kernel (less bouncy liquid)
''        better broad phase collision detection
''        better mouse interaction
''        better code structure

''        includes
#Include Once "GL/gl.bi"

''        constants
Const As Single                dt                                                = 0.01                                                                                        ''        timestep
Const As Single                half_dt                                = 0.5*dt                                                                                ''        half timestep
Const As Single                r                                                 = 10                                                                                                ''        interaction radius
Const As Single                r_2                                         = 2*r                                                                                                ''        two radius (for grid)
Const As Single                r_sqrd                                 = r*r                                                                                                ''        radius squared
Const As Single                k                                                        = 400                                                                                                ''        global "spring stiffnes"
Const As Single                k_near                                = 3000                                                                                        ''        local "spring stiffnes"
Const As Single                viscosity                        = 2                                                                                                        ''        viscosity
Const As Single                rest_dens                 = 6                                                                                                        ''        rest density
Const As Single                grav_acc                        = 400                                                                                                ''        gravity
Const As Single                pi                                                = 4*Atn(1)                                                                        ''        pi
Const As Single                radius                                = 2                                                                                                        ''        particle distance on startup
Const As Integer  block_width                = 60                                                                                                ''        <--- change to reduce / increase num masses
Const As Integer  block_height        = 100                                                                                                ''        
Const As Integer  num_masses                = block_width*block_height+100        ''        total number of masses
Const As Integer  neighbors                        = 48                                                                                                ''        max number of neighbors
Const As Integer         screen_wid                 = 800                                                                                                ''        screen width
Const As Integer         screen_hgt                 = 400                                                                                                ''        screen height
Const As Single         wall_friction        = 1                                                                                                        ''        wall friction (1 = full, 0 = off)
Const As Single         interaction_Radius        = 8                                                                                ''        mouse interaction radius

Const As Integer cell_dia                                                = r
Const As Integer cell_max_masses                = 24
Const As Integer num_cells_row                        = screen_wid/cell_dia
Const As Integer num_cells_col                        = screen_hgt/cell_dia

Const As Integer num_cells                                        = num_cells_row*num_cells_col

Const As Single cell_wid                                                = screen_wid/num_Cells_row
Const As Single cell_hgt                                                = screen_hgt/num_Cells_Col
Const As Single half_cell_wid                                = cell_wid/2

''        types
Type CellType
        As Integer Num_Masses, Mass(1 To cell_max_masses)
End Type

Type vectortype
        As Single x, y
End Type

Type NeighborType
        As Integer a
        As Single Dist
        As vectortype norm_dist
End Type

Type masstype
        As Integer num_neighbors
        As NeighborType Neighbor(1 To neighbors)
        As vectortype frc, vel, psn
        As Single density, angvel, torque
End Type

''        dim variables
Dim As Integer         center_x         = screen_wid\2
Dim As Integer         center_y        = screen_hgt\2

Dim As Integer mouse_x, mouse_y, mouse_x_old, mouse_y_old, mouse_vel_x, mouse_vel_y

Dim Shared As masstype mass(1 To num_masses)
Dim Shared As celltype cell(1 To num_cells_row, 1 To num_cells_col)

''        narrow phase collision detection. Make particles neighbors if distance < r
Sub CheckIfNeighbor(Byval a_ As Integer, Byval b_ As Integer)
        
        If mass(a_).num_neighbors = neighbors Then Exit Sub
        If mass(b_).num_neighbors = neighbors Then Exit Sub
        
        Dim As vectortype dist
        
        dist.x = mass(b_).psn.x-mass(a_).psn.x: If Abs(dist.x) > r Then Exit Sub
        dist.y = mass(b_).psn.y-mass(a_).psn.y: If Abs(dist.y) > r Then Exit Sub
        
        Dim As Single dist_sqrd = dist.x*dist.x+dist.y*dist.y
        
        If dist_sqrd > r_sqrd Or dist_sqrd = 0 Then Exit Sub
        
        Dim As Single distance = Sqr(dist_sqrd)
        
        Dim As vectortype norm_dist = (dist.x/distance, dist.y/distance)
        
        Dim As Single density = 1 - distance / r
        
        '' create neighbors and save data
        With mass(a_)
                .num_neighbors += 1
                .Neighbor(.num_neighbors).a = b_
                .Neighbor(.num_neighbors).dist = distance
                .Neighbor(.num_neighbors).norm_dist.x = norm_dist.x
                .Neighbor(.num_neighbors).norm_dist.y = norm_dist.y
                .Density += density
        End With
        
        With mass(b_)
                .num_neighbors += 1
                .Neighbor(.num_neighbors).a = a_
                .Neighbor(.num_neighbors).dist = distance
                .Neighbor(.num_neighbors).norm_dist.x = -norm_dist.x
                .Neighbor(.num_neighbors).norm_dist.y = -norm_dist.y
                .Density += density
        End With
        
End Sub

Randomize Timer

Dim As Integer n = 1

''        place particles
For i As Integer = 1 To block_width
        For j As Integer = 1 To block_height
                With mass(n)
                        .psn.x = i * radius * 2
                        .psn.y = screen_hgt - j * radius * 2
                        
                        .psn.x += (Rnd-Rnd)
                        .psn.y += (Rnd-Rnd)
                        
                        .vel.x += (Rnd-Rnd)
                        .vel.y += (Rnd-Rnd)

                End With
                n += 1
        Next
Next

''        create opengl screen
ScreenRes screen_wid, screen_hgt, 32,, 2
glMatrixMode(GL_PROJECTION)  
glLoadIdentity()
glViewport(0, 0, screen_wid, screen_hgt)
glOrtho(0, screen_wid, screen_hgt, 0, 0, 32)
glMatrixMode(GL_MODELVIEW)
glLoadIdentity()
glClearColor (0.6, 0.7, 1.0, 0.0)
'glenable(GL_ALPHA)
glEnable(GL_BLEND)
glblendfunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
'glblendfunc(GL_DST_COLOR , GL_SRC_COLOR)

getmouse(mouse_x, mouse_y)

''        main loop
Do
        
        '' save old mouse state
        mouse_x_old = mouse_x
        mouse_y_old = mouse_y
        
        ''        get new mouse state
        getmouse(mouse_x, mouse_y)
        
        '' calc mouse velocity
        mouse_vel_x = (mouse_x-mouse_x_old)/dt
        mouse_vel_y = (mouse_y-mouse_y_old)/dt
        
        ''        mouse interaction
        For i As Integer = 1 To num_masses
                With mass(i)        
                        
                        Dim As vectortype dst
                        
                        dst.x = .psn.x - mouse_x:        If dst.x > interaction_radius Then Continue For
                        dst.y = .psn.y - mouse_y:        If dst.y > interaction_radius Then Continue For
                        
                        If dst.x*dst.x+dst.y*dst.y > interaction_Radius*interaction_Radius Then Continue For
                        
                        If mouse_vel_x Then .vel.x += (mouse_vel_x-.vel.x) * 0.16
                        If mouse_vel_y Then .vel.y += (mouse_vel_y-.vel.y) * 0.16
                        
                End With        
        Next        
        
        ''        reset particles
        For i As Integer = 1 To num_masses
                With mass(i)
                        .density = 0
                        .num_neighbors = 0
                End With
        Next
        
        ''        reset cells
        For i As Integer = 1 To num_cells_row
                For j As Integer = 1 To num_cells_col
                        cell(i, j).Num_Masses = 0
                Next
        Next
        
        ''        update grid (broad phase collision detection)
        For i As Integer = 1 To num_masses
                Dim As Integer cell_row = (mass(i).Psn.y\cell_hgt)+1
                Dim As Integer cell_col = (mass(i).Psn.x\cell_wid)+1
                With cell(cell_col, cell_row)
                        If .num_masses = cell_max_masses Then Continue For
                        .num_masses += 1
                        .Mass(.num_masses) = i
                End With
        Next
        
        ''        narrow phase collision detection
        For i As Integer = 1 To num_cells_row
                For j As Integer = 1 To num_cells_col
                        If cell(i, j).num_masses = 0 Then Continue For
                        For k_ As Integer = 1 To cell(i, j).num_masses
                                Dim As Integer a_ = cell(i, j).mass(k_)
                                
                                If mass(a_).num_neighbors = neighbors Then Continue For
                                
                                ''        cell(x, y) - self 
                                For l_ As Integer = k_+1 To cell(i, j).num_masses
                                        CheckIfNeighbor(a_, cell(i, j).mass(l_))
                                Next
                                
                                ''        cell(x, y) - cell(x+1, y)
                                If i < num_cells_row Then
                                        For l_ As Integer = 1 To cell(i+1, j).num_masses
                                                CheckIfNeighbor(a_, cell(i+1, j).mass(l_))
                                        Next
                                End If
                                
                                If j = num_cells_col Then Continue For
                                
                                ''        cell(x, y) - cell(x, y+1)
                                For l_ As Integer = 1 To cell(i, j+1).num_masses
                                        CheckIfNeighbor(a_, cell(i, j+1).mass(l_))
                                Next
                                
                                ''        cell(x, y) - cell(x+1, y+1)
                                If i < num_cells_row Then
                                        For l_ As Integer = 1 To cell(i+1, j+1).num_masses
                                                CheckIfNeighbor(a_, cell(i+1, j+1).mass(l_))
                                        Next
                                End If
                                
                                ''        cell(x, y) - cell(x-1, y+1)
                                If i > 1 Then
                                        For l_ As Integer = 1 To cell(i-1, j+1).num_masses
                                                CheckIfNeighbor(a_, cell(i-1, j+1).mass(l_))
                                        Next
                                End If
                                
                        Next
                Next
        Next
                
        ''        set pressure
        For i As Integer = 1 To num_masses
                
                If mass(i).num_neighbors = 0 Then Continue For

                ''        mass i global pressure
                Dim As Single pressure = -k * (mass(i).density - rest_dens)
                
                For j As Integer = 1 To mass(i).num_neighbors
                        
                        Dim As Single gradient = 1 - mass(i).neighbor(j).dist / r
                                
                        Dim As Integer l = mass(i).neighbor(j).a
                        
                        ''        mass i + l global pressure
                        Dim As Single press = (pressure - k * (mass(l).density - rest_dens) )
                        
                        ''        mass i + l local pressure
                        press -= k_near * gradient
                        
                        ''        viscosity
                        Dim As vectortype vel = (mass(i).vel.x-mass(l).vel.x, mass(i).vel.y-mass(l).vel.y)
                                
                        Dim As Single proj_vel = vel.x*mass(i).neighbor(j).norm_dist.x+vel.y*mass(i).neighbor(j).norm_dist.y
                        
                        press -=  viscosity * proj_vel * gradient
                        
                        press *= gradient
                                
                        mass(i).frc.x += press * mass(i).neighbor(j).norm_dist.x
                        mass(i).frc.y += press * mass(i).neighbor(j).norm_dist.y
                        
                        mass(l).frc.x -= press * mass(i).neighbor(j).norm_dist.x
                        mass(l).frc.y -= press * mass(i).neighbor(j).norm_dist.y
                        
                Next
                
        Next
        
        ''        integrate (since particle mass = 1, force and acceleration are identical)
        For i As Integer = 1 To num_masses
                With mass(i)
                        
                        .frc.y -= grav_acc
                        
                        .vel.x += .frc.x * half_dt * cos(.vel.x) - sin(.vel.y)
                        .vel.y += .frc.y * half_dt * cos(.vel.y) + sin(.vel.x)
                        
                        .angvel -= .torque * half_dt
                        
                        .psn.x += .vel.x * dt + .frc.x * dt * half_dt
                        .psn.y += .vel.y * dt + .frc.y * dt * half_dt
                        
                        .vel.x += .frc.x * half_dt
                        .vel.y += .frc.y * half_dt
                        
                        .angvel += .torque * half_dt
                        
                        .frc.x = 0
                        .frc.y = 0
                        
                        .torque = 0
                        
                End With
        Next
        
        ''        keep particles inside screen
        For i As Integer = 1 To num_masses
                With mass(i)
                        
                        If .psn.x < radius Then .psn.x = radius: .vel.x = -.vel.x*wall_friction: .vel.y *= wall_friction
                        If .psn.x > screen_wid - radius Then .psn.x = screen_wid -radius: .vel.x = -.vel.x*wall_friction: .vel.y *= wall_friction
                        
                        If .psn.y < radius Then .psn.y = radius: .vel.y = -.vel.y*wall_friction: .vel.x *= wall_friction
                        If .psn.y > screen_hgt - radius Then .psn.y = screen_hgt -radius: .vel.y = -.vel.y*wall_friction: .vel.x *= wall_friction
                        
                End With
        Next
        
        ''        draw to screen
        glClear(GL_DEPTH_BUFFER_BIT Or GL_COLOR_BUFFER_BIT)
        
        glLoadIdentity()
        
        glenable(GL_POINT_SMOOTH)
        glpointsize(r)
        glcolor4f(0.0f, 0.0f, 0.0f, 0.25f)
        
        GLbegin(GL_POINTS)
                For i As Integer = 1 To num_masses
                        glvertex2f(mass(i).psn.x, mass(i).psn.y)
                Next
        GLend()
        
        gldisable(GL_POINT_SMOOTH)
        glpointsize(1)
        glcolor4f(1.0f, 1.0f, 0.0f, 0.75f)
        
        glbegin(GL_POINTS)
                For i As Integer = 1 To num_masses
                        glvertex2f(mass(i).psn.x, mass(i).psn.y)
                Next
        GLend()
        
        GLflush()
        
        Flip()
        
        Sleep 1, 1
        
Loop Until MultiKey(1)

End
 

[Voltage]

Absolutely love it h4tt3n. I'm a sucker for graphics code, and specifically fluid dynamics at the moment.

I've only breifly looked at the code, but a couple of things that may help:

- Better broadphase collision detection (I see it on the TODO)
- Do you need the SQR? Or can you just test on the squared result. <-- maybe answered already in the code which I haven't studied.
- I'd run a profiler over it to see where optimizations would be best focused.

I hope to have a play around with this a bit later. If I come up with anything interesting I'll post it here.

[h4tt3n]

@ Rolliebollocks - Too easy! :-7

@Voltage -

A better broad phase will definitely improve this. It just looks better with more particles. My next step is to implement the "sweep and prune" or "linear linked list" methods, which should be able to handle somewhat larger number of particles.

I've tried implementing a sqr() less version with *some* success, but I haven't worked all the math out yet - not even sure you can. It wasn't much faster btw, so I don't even think it's worth the trouble.

I'm looking forward to see what you figure out.

Cheers,
Mike

[rolliebollocks]

I wanted to see if it would create ripples or something! Plus I got polygons to fill and I didn't want to get distracted.

[kinem]

Very cool

[kiyotewolf]

O.o How do you come up with this stuff..?


Ohhh yeah.. physics..



~Kiyote!

[rdc]

Hello folks,

Some time ago I posted a simple smoothed particle hydrodynamics example, and recently I've made it somewhat better. This code sample uses grid based collision detection and handles up to around 10.000 particles depending on your system. Ideas, suggestions and critique are welcome :-)

Cheers,
Mike

I just have one word for you: Geek! :) Totally awesome Mike. Very smooth execution on this. I too love this kind of stuff.

[mrToad]

Awesome :D
Wish I had the mind to make this kind of thing.

[h4tt3n]

Thanks guys, I'm glad you like it. Here's a sph primer if you're intersted:

http://www8.cs.umu.se/kurser/5DV058/VT1 ... cture8.pdf

You'll notice that my impementation is fairly simple (if not naive) compared to this :-)

Cheers,
Mike

[h4tt3n]

Hello folks,

Okay, after 1½ years I had another look at my old SPH experiments and tuned them up a bit. Here's what I got so far. It's still not perfect, there are a few ugly hacks I need to iron our, but it's more realistic and handles a *lot* more particles than my previous samples. This is partly because I got it running without calling the sqr() function. If you have a decent machine, it'll run ahem... fluently with around 20.000 particles.

First, the vector library. Save as "vec2.bi"

Code: Select all

''*******************************************************************************
''		
''		Freebasic 2d floating point and integer vector library
''		version 0.6b, march 2011, Michael "h4tt3n" Nissen, jernmager@yahoo.dk
''		Integer vectors have been added for screen and mouse operations.   
''		
''		function syntax:
''		
''   	(return type) (function name) (argument type (, ...))
''		
''	 	floating point vector function list:
''		
''   	vector absolute          	(vector)          				- absolute value
''   	vector normal            	(vector)          				- normal vector
''   	vector normalised        	(vector)          				- normalised vector
''		vector normalisednormal		(vector)									-	normalised normal vector
''   	scalar magnitude         	(vector)          				- magnitude
''   	scalar magnitudesquared		(vector)          				- magnitude squared
''   	scalar distance          	(vector, vector)  				- vector distance
''   	scalar distancesquared   	(vector, vector)  				- vector distance squared
''   	scalar dot               	(vector, vector)  				- dot product
''   	scalar dotnormal         	(vector, vector)  				- normal dot product
''   	vector project           	(vector, vector)					-	vector projection
''		vector component					(vector, vector)					-	vector component
''		vector randomise					(scalar)									-	randomise in range +/- value
''		vector lerp								(vector, vector, scalar)	- linear interpolation
''
''		integer vector function list:
''		
''		
''	 	function useage, member and non-member style:
''
''		vector_a.function(vector_b),	function(vector_a, vector_b)
''
''		
''*******************************************************************************

''	
Type float As Single
Type integ As Integer

''  vec2f 2d float vector structure
Type vec2f
	
  ''  vec2f variables
  As float x, y
	
  ''  vec2f constructor declarations
	Declare Constructor ()
  Declare Constructor (ByVal x As float, ByVal y As float)
  Declare Constructor (ByRef vec As vec2f)
	
  ''  vec2f compound arithmetic member operator declarations
  Declare Operator += (ByRef rhs As vec2f)
  Declare Operator -= (ByRef rhs As vec2f)
  Declare Operator *= (ByRef rhs As vec2f)
  Declare Operator *= (ByRef rhs As float)
  Declare Operator /= (ByRef rhs As float)
  Declare Operator Let (ByRef rhs As vec2f)
	
	''	vec2f member function declarations
  Declare Function absolute() As vec2f
  Declare Function normal() As vec2f 
  Declare Function normalised() As vec2f
  Declare Function normalisednormal() As vec2f
  Declare Function magnitude() As float
  Declare Function magnitudesquared() As float
  Declare Function distance(ByRef rhs As vec2f) As float
  Declare Function distancesquared(ByRef rhs As vec2f) As float
  Declare Function dot(ByRef rhs As vec2f) As float
  Declare Function dotnormal(ByRef rhs As vec2f) As float
  Declare Function project(ByRef rhs As vec2f) As vec2f
  Declare Function component(ByRef rhs As vec2f) As vec2f
  Declare Function randomise(ByVal rhs As float) As vec2f
  Declare Function lerp(ByRef rhs As vec2f, ByVal i As float) As vec2f
  
End Type

''  vec2f unary arithmetic non-member operator declarations
Declare Operator - (ByRef rhs As vec2f) As vec2f

''  vec2f binary arithmetic non-member operator declarations
Declare Operator + (ByRef lhs As vec2f, ByRef rhs As vec2f) As vec2f
Declare Operator - (ByRef lhs As vec2f, ByRef rhs As vec2f) As vec2f
Declare Operator * (ByVal lhs As float, ByRef rhs As vec2f) As vec2f
Declare Operator * (ByRef lhs As vec2f, ByVal rhs As float) As vec2f
Declare Operator / (ByRef lhs As vec2f, ByVal rhs As float) As vec2f

''  vec2f non-member function declarations
Declare Function absolute (ByRef lhs As vec2f) As vec2f
Declare Function normal (ByRef lhs As vec2f) As vec2f
Declare Function normalised (ByRef lhs As vec2f) As vec2f
Declare Function normalisednormal(ByRef lhs As vec2f) As vec2f
Declare Function magnitude (ByRef lhs As vec2f) As float
Declare Function magnitudesquared (ByRef lhs As vec2f) As float
Declare Function distance (ByVal lhs As vec2f, ByRef rhs As vec2f) As float
Declare Function distancesquared (ByVal lhs As vec2f, ByRef rhs As vec2f) As float
Declare Function dot (ByVal lhs As vec2f, ByRef rhs As vec2f) As float
Declare Function dotnormal (ByVal lhs As vec2f, ByRef rhs As vec2f) As float
Declare Function project (ByVal lhs As vec2f, ByRef rhs As vec2f) As vec2f
Declare Function component(ByVal lhs As vec2f, ByRef rhs As vec2f) As vec2f
Declare Function trigonometry(ByRef lhs As float) As vec2f
Declare Function randomise(ByVal lhs As float) As vec2f
Declare Function lerp(ByVal lhs As vec2f, ByRef rhs As vec2f, ByVal i As float) As vec2f

''  vec2f constructors
Constructor vec2f(): This.x = 0.0: This.y = 0.0: End Constructor
Constructor vec2f(ByVal x As float, ByVal y As float): This.x = x: This.y = y: End Constructor
Constructor vec2f(ByRef vec As vec2f): This.x = vec.x: This.y = vec.y: End Constructor

''  vec2f compound arithmetic member operators
Operator vec2f.+= (ByRef rhs As vec2f): This.x += rhs.x: This.y += rhs.y: End Operator
Operator vec2f.-= (ByRef rhs As vec2f): This.x -= rhs.x: This.y -= rhs.y: End Operator
Operator vec2f.*= (ByRef rhs As vec2f): This.x *= rhs.x: This.y *= rhs.y: End Operator
Operator vec2f.*= (ByRef rhs As float): This.x *= rhs: This.y *= rhs: End Operator
Operator vec2f./= (ByRef rhs As float): This.x /= rhs: This.y /= rhs: End Operator
Operator vec2f.let (ByRef rhs As vec2f): This.x = rhs.x: This.y = rhs.y: End Operator

''  vec2f member functions
Function vec2f.absolute() As vec2f: Return vec2f(Abs(This.x), Abs(This.y)): End Function
Function vec2f.normal() As vec2f: Return vec2f(This.y, -This.x): End Function 
'Function vec2f.normalised() As vec2f: If This.x <> 0 And This.y <> 0 Then Return This/This.magnitude(): Else Return vec2f(0, 0): End If: End Function
Function vec2f.normalised() As vec2f: If This.x And This.y = 0 Then Return vec2f(): Else Return This/This.magnitude(): End If: End Function
Function vec2f.normalisednormal() As vec2f: Return This.normal()/magnitude(): End Function
Function vec2f.magnitude() As float: Return Sqr(This.magnitudesquared()): End Function
Function vec2f.magnitudesquared() As float: Return This.dot(This): End Function
Function vec2f.distance(ByRef rhs As vec2f) As float: Return Sqr(This.distancesquared(rhs)): End Function
Function vec2f.distancesquared(ByRef rhs As vec2f) As float: Return (This-rhs).dot(This-rhs): End Function
Function vec2f.dot(ByRef rhs As vec2f) As float: Return (This.x*rhs.x+This.y*rhs.y): End Function
Function vec2f.dotnormal(ByRef rhs As vec2f) As float: Return This.dot(rhs.normal()): End Function
Function vec2f.project(ByRef rhs As vec2f) As vec2f: Return (This.dot(rhs)/This.magnitudesquared())*rhs: End Function
Function vec2f.component(ByRef rhs As vec2f) As vec2f: Return (This.dot(rhs)/rhs.magnitudesquared)*rhs: End Function
Function vec2f.randomise(ByVal rhs As float) As vec2f: Return vec2f((Rnd-Rnd)*rhs, (Rnd-Rnd)*rhs): End Function
Function vec2f.lerp(ByRef rhs As vec2f, ByVal i As float) As vec2f: Return This+(rhs-This)*i: End Function

''  vec2f unary arithmetic non-member operators
Operator - (ByRef rhs As vec2f) As vec2f: Return vec2f(-rhs.x, -rhs.y): End Operator

''  vec2f binary arithmetic non-member operators
Operator + (ByRef lhs As vec2f, ByRef rhs As vec2f) As vec2f: Return vec2f(lhs.x+rhs.x, lhs.y+rhs.y): End Operator
Operator - (ByRef lhs As vec2f, ByRef rhs As vec2f) As vec2f: Return vec2f(lhs.x-rhs.x, lhs.y-rhs.y): End Operator
Operator * (ByVal lhs As float, ByRef rhs As vec2f) As vec2f: Return vec2f(lhs*rhs.x, lhs*rhs.y): End Operator
Operator * (ByRef lhs As vec2f, ByVal rhs As float) As vec2f: Return vec2f(lhs.x*rhs, lhs.y*rhs): End Operator
Operator / (ByRef lhs As vec2f, ByVal rhs As float) As vec2f: Return vec2f(lhs.x/rhs, lhs.y/rhs): End Operator

''  vec2f non-member functions
Function absolute (ByRef lhs As vec2f) As vec2f: Return lhs.absolute(): End Function
Function normal (ByRef lhs As vec2f) As vec2f: Return lhs.normal(): End Function
Function normalised (ByRef lhs As vec2f) As vec2f: Return lhs.normalised(): End Function
Function normalisednormal(ByRef lhs As vec2f) As vec2f: Return lhs.normalisednormal(): End Function
Function magnitude (ByRef lhs As vec2f) As float: Return lhs.magnitude(): End Function
Function magnitudesquared (ByRef lhs As vec2f) As float: Return lhs.magnitudesquared(): End Function
Function distance (ByVal lhs As vec2f, ByRef rhs As vec2f) As float: Return lhs.distance(rhs): End Function
Function distancesquared (ByVal lhs As vec2f, ByRef rhs As vec2f) As float: Return lhs.distancesquared(rhs): End Function
Function dot (ByVal lhs As vec2f, ByRef rhs As vec2f) As float: Return lhs.dot(rhs): End Function
Function dotnormal (ByVal lhs As vec2f, ByRef rhs As vec2f) As float: Return lhs.dotnormal(rhs): End Function
Function project (ByVal lhs As vec2f, ByRef rhs As vec2f) As vec2f: Return lhs.project(rhs): End Function
Function component(ByVal lhs As vec2f, ByRef rhs As vec2f) As vec2f: Return lhs.component(rhs): End Function
Function trigonometry(ByRef lhs As float) As vec2f: Return vec2f(Cos(lhs), Sin(lhs)): End Function
Function randomise(ByVal lhs As float) As vec2f: Return vec2f((Rnd-Rnd)*lhs, (Rnd-Rnd)*lhs): End Function
Function lerp(ByVal lhs As vec2f, ByRef rhs As vec2f, ByVal i As float) As vec2f: Return lhs.lerp(rhs, i): End Function

''  vec2i 2d integer vector structure
Type vec2i
	
  ''  vec2i variables
  As integ x, y
	
  ''  vec2i constructor declarations
	Declare Constructor ()
  Declare Constructor (ByVal x As integ, ByVal y As integ)
  Declare Constructor (ByRef vec As vec2i)
	
  ''  vec2i compound arithmetic member operator declarations
  Declare Operator += (ByRef rhs As vec2i)
  Declare Operator -= (ByRef rhs As vec2i)
  Declare Operator *= (ByRef rhs As vec2i)
  Declare Operator *= (ByVal rhs As integ)
  Declare Operator \= (ByVal rhs As integ)
  Declare Operator Let (ByRef rhs As vec2i)
  
  ''  vec2i member function declarations
  
End Type

''  vec2i unary arithmetic non-member operator declarations
Declare Operator - (ByRef rhs As vec2i) As  vec2i

''  vec2i binary arithmetic non-member operator declarations
Declare Operator + (ByRef lhs As vec2i, ByRef rhs As vec2i) As vec2i
Declare Operator - (ByRef lhs As vec2i, ByRef rhs As vec2i) As vec2i
Declare Operator * (ByVal lhs As float, ByRef rhs As vec2i) As vec2i
Declare Operator * (ByRef lhs As vec2i, ByVal rhs As float) As vec2i
Declare Operator \ (ByRef lhs As vec2i, ByVal rhs As float) As vec2i

''  vec2i non-member function declarations

''  vec2i constructors
Constructor vec2i(): This.x = 0: This.y = 0: End Constructor
Constructor vec2i(ByVal x As integ, ByVal y As integ): This.x = x: This.y = y: End Constructor
Constructor vec2i(ByRef vec As vec2i): This.x = vec.x: This.y = vec.y: End Constructor

''  vec2i compound arithmetic member operators
Operator vec2i.+= (ByRef rhs As vec2i): This.x += rhs.x: This.y += rhs.y: End Operator
Operator vec2i.-= (ByRef rhs As vec2i): This.x -= rhs.x: This.y -= rhs.y: End Operator
Operator vec2i.*= (ByRef rhs As vec2i): This.x *= rhs.x: This.y *= rhs.y: End Operator
Operator vec2i.*= (ByVal rhs As integ): This.x *= rhs: This.y *= rhs: End Operator
Operator vec2i.\= (ByVal rhs As integ): This.x \= rhs: This.y \= rhs: End Operator
Operator vec2i.let (ByRef rhs As vec2i): This.x = rhs.x: This.y = rhs.y: End Operator

''  vec2i member functions

''  vec2i unary arithmetic non-member operators
Operator - (ByRef rhs As vec2i) As vec2i: Return vec2i(-rhs.x, -rhs.y): End Operator

''  vec2i binary arithmetic non-member operators
Operator + (ByRef lhs As vec2i, ByRef rhs As vec2i) As vec2i: Return vec2i(lhs.x+rhs.x, lhs.y+rhs.y): End Operator
Operator - (ByRef lhs As vec2i, ByRef rhs As vec2i) As vec2i: Return vec2i(lhs.x-rhs.x, lhs.y-rhs.y): End Operator
Operator * (ByVal lhs As integ, ByRef rhs As vec2i) As vec2i: Return vec2i(lhs*rhs.x, lhs*rhs.y): End Operator
Operator * (ByRef lhs As vec2i, ByVal rhs As integ) As vec2i: Return vec2i(lhs.x*rhs, lhs.y*rhs): End Operator
Operator \ (ByRef lhs As vec2i, ByVal rhs As integ) As vec2i: Return vec2i(lhs.x\rhs, lhs.y\rhs): End Operator

''  vec2i non-member functions

''  shared binary arithmetic non-member operator declarations
Declare Operator + (ByRef lhs As vec2f, ByRef rhs As vec2i) As vec2f
Declare Operator + (ByRef lhs As vec2i, ByRef rhs As vec2f) As vec2f
Declare Operator - (ByRef lhs As vec2f, ByRef rhs As vec2i) As vec2f
Declare Operator - (ByRef lhs As vec2i, ByRef rhs As vec2f) As vec2f

''  shared non-member function declarations

''  shared binary arithmetic non-member operators
Operator + (ByRef lhs As vec2f, ByRef rhs As vec2i) As vec2f: Return vec2f(lhs.x+rhs.x, lhs.y+rhs.y): End Operator
Operator + (ByRef lhs As vec2i, ByRef rhs As vec2f) As vec2f: Return vec2f(lhs.x+rhs.x, lhs.y+rhs.y): End Operator
Operator - (ByRef lhs As vec2f, ByRef rhs As vec2i) As vec2f: Return vec2f(lhs.x-rhs.x, lhs.y-rhs.y): End Operator
Operator - (ByRef lhs As vec2i, ByRef rhs As vec2f) As vec2f: Return vec2f(lhs.x-rhs.x, lhs.y-rhs.y): End Operator

''  shared non-member functions

And here the sph example. Have fun :-)

Code: Select all

''*******************************************************************************
''  Smoothed Particle Hydrodynamics (SPH) fluid Simulation
''  
''  Version 8, december 6th 2011
''  Mike "h4tt3n", micha3l_niss3n@yahoo.dk
''
''  Description:
''
''  
''  Controls:
''
''  press left mouse button to stir fluid
''  roll mouse wheel up / down to create current eddies
''  press escape key to exit
''      
''  Reference:
''  
''  http://image.diku.dk/projects/media/kelager.06.pdf
''  http://www8.cs.umu.se/kurser/5DV058/VT10/lectures/Lecture8.pdf
''  http://www.cs.clemson.edu/~bpelfre/sph_tutorial.pdf
''  http://graphics.stanford.edu/projects/lgl/papers/apkg-aspf-sig07/apkg-aspf-sig07.pdf
''
''  Todo-list:
''  
''  
''
''      
''*******************************************************************************

''   Includes
#Include Once "vec2.bi"
#Include Once "GLFW.bi"

''   global constants
Const As Integer  TRUE                = -1                           ''   boolean true
Const As Integer  FALSE               = Not TRUE                     ''   boolean false
Const As float  dt                  = 0.01                         ''   
Const As float  inv_dt              = 1 / dt                       ''   
Const As float  half_dt             = 0.5 * dt                     ''   
Const As float  dt_Sqrd             = dt * dt                      ''   
Const As float  G                   = 400                          ''   
Const As float  zero                = 4                            ''   min float value
Const As float  pi                  = 4*Atn(1)                     ''   pi
Const As float  r                   = 10                           ''   particle interaction radius
Const As float  inv_r               = 1 / r                        ''   particle interaction radius
Const As float  r_sqrd              = r * r                        ''  radius squared
Const As float  inv_r_sqrd          = 1 / r_sqrd                   ''  inverse radius squared
Const As float  k                   = 400000                       ''   global stiffnes coefficient
Const As float  k_near              = 400                         ''   local stiffnes coefficient
Const As float  mju                 = 4000                         ''   viscosity coefficient
Const As float  rest_density        = 4                            ''	rest density
Const As Integer  block_width         = 120                          ''                                                                                     ''        <--- change to reduce / increase num masses
Const As Integer  block_height        = 100                          ''   
Const As Integer  num_particles       = block_width*block_height     ''   number of particles
Const As Integer  max_neighbors       = 32                           ''   max number of neighbors per particles
Const As Integer  max_num_neighbors   = num_particles*max_neighbors  ''   global max number of neighbors

Const As Integer  screen_wid          = 900                          ''   screen width
Const As Integer  screen_hgt          = 600                          ''   screen height

Const As Integer  cell_dia            = 0.75 * r                     ''   cell diameter
Const As Integer  max_cell_particles  = 16                           ''   max number of particles held by one cell
Const As Integer  num_cells_row       = screen_wid/cell_dia          ''   number of cells per row / on x axis
Const As Integer  num_cells_col       = screen_hgt/cell_dia          ''   number of cells per column / on y axis
Const As Integer  num_cells           = num_cells_row*num_cells_col  ''	  global number of cells
Const As float  cell_wid            = screen_wid/num_cells_row     ''  cell width
Const As float  cell_hgt            = screen_hgt/num_cells_col     ''  cell height

Const As Integer  border              = r                            ''
Const As Integer  r_plus_border       = r + border                   ''
Const As float  border_k            = 3600                         ''
Const As float  border_d            = 0                            ''

Const As float  interaction_radius  = 32                           '' mouse interaction radius
Const As float  int_radius_sqrd     = interaction_radius^2         ''

''
Type ParticleType
   As Vec2f psn
   As Vec2f psn_old
   As Vec2f vel
   As Vec2f frc
   As float density
   As float pressure
   As float volume
   As float mass
   As Integer num_neighbors
End Type

''   
Type CellType
   As Integer Num_Particles
   As Integer Particle(1 To max_cell_particles)
   As Integer NumCellNeighbors
   As CellType Ptr CellNeighbor(1 To 4)
End Type

''
Type NeighborType
   As Integer a
   As Integer b
   As float distance_squared
   As vec2f n_over_distance
End Type

''
Type ScreenType
   
   Declare Function CreateScreen(w As Integer, h As Integer, r As Integer, g As Integer, _
      b As Integer, a As Integer, d As Integer, s As Integer, m As Integer) As Integer
   Declare Function DeleteScreen() As Integer
   
   As Integer wid
   As Integer Hgt
   As Integer Bpp
   
End Type

''
Type MouseType
	
   As vec2i Psn
   As vec2i Psn_Old
   As vec2f vel
   As Integer Lbt
   As Integer Lbt_Old
   As Integer Rbt
   As Integer Rbt_Old
   As Integer Whl
   As Integer Whl_Old
   
   Declare sub Update()
   
End Type


Dim Shared As ScreenType Scr
Dim Shared As MouseType Mouse

Declare Sub run_simulation()
Declare Sub initiate_simulation()
Declare Sub broad_phase_collision()
Declare Sub narrow_phase_collision()
Declare Sub reset_cells()
Declare Sub reset_particles()
Declare Sub calculate_pressure()
Declare Sub calculate_internal_force()
Declare Sub screen_border_force()
Declare Sub Integrate()
Declare Sub draw_particles()

Declare Function Check_if_neighbor(ByVal a As Integer, ByVal b As Integer) As Integer

Dim Shared As CellType cell(1 To num_cells_row, 1 To num_cells_col)
Dim Shared As ParticleType particle(1 To num_particles)
Dim Shared As NeighborType Neighbor(1 To max_num_neighbors)

Dim Shared As Integer num_neighbors

''   run simulation
initiate_simulation()
Scr.CreateScreen(screen_wid, screen_hgt, 8, 8, 8, 8, 0, 0, GLFW_window)
run_simulation()
Scr.DeleteScreen

End


Sub run_simulation()
   
   ''	main program loop
   Do
      
      reset_particles()
      reset_cells()
      
      broad_phase_collision()
      narrow_phase_collision()
      
      calculate_pressure()
      calculate_internal_force()
      screen_border_force()
      mouse.Update()
      
      Integrate()
      
      draw_particles()
      
      'glfwsleep(0.01)
      
   Loop While (Not glfwgetkey(GLFW_key_esc) And glfwgetwindowparam(GLFW_opened))
   
End Sub

Sub initiate_simulation()
	
	''	dereference cell neighbor pointers
	For x As Integer = 1 To num_cells_row
      For y As Integer = 1 To num_cells_col
         
			With Cell(x, y)
				
				.numCellNeighbors = 0
				
				''	cell(x, y) - cell(x+1, y)
				If x < num_cells_row Then
					.numCellNeighbors += 1
					.CellNeighbor(.numCellNeighbors) = @Cell(x+1, y)
				EndIf
				
				''	cell(x, y) - cell(x, y+1)
				If y < num_cells_col Then
					.numCellNeighbors += 1
					.CellNeighbor(.numCellNeighbors) = @Cell(x, y+1)
				EndIf
				'
				''	cell(x, y) - cell(x+1, y+1)
				If x < num_cells_row And y < num_cells_col Then
					.numCellNeighbors += 1
					.CellNeighbor(.numCellNeighbors) = @Cell(x+1, y+1)
				EndIf
				
				''	cell(x, y) - cell(x-1, y+1)
				If x > 1 And y < num_cells_col Then
					.numCellNeighbors += 1
					.CellNeighbor(.numCellNeighbors) = @Cell(x-1, y+1)
				EndIf
				
			End With
         
      Next
   Next
	
   
   Randomize Timer
   
   Dim As Integer n = 0

   ''   place particles
   For x As Integer = 1 To block_width
      For y As Integer = 1 To block_height
         
         n += 1
         
         With Particle(n)
            
            .psn.x = border + r * 0.5 + x * r * 0.5
            .psn.y = screen_hgt - border - y * r * 0.5
            
            .psn.x += (Rnd-Rnd) * r * 0.1
            .psn.y += (Rnd-Rnd) * r * 0.1
            
            .vel.x += (Rnd-Rnd) * 10
            .vel.y += (Rnd-Rnd) * 10
            
            .psn_old.x = .psn.x' + (Rnd-Rnd)*2
            .psn_old.y = .psn.y' + (Rnd-Rnd)*2
            
         End With
      Next
   Next
   
End Sub

Sub broad_phase_collision()
     
   ''	
   For i As Integer = 1 To num_particles
      
      Dim As Integer cell_row = Particle(i).Psn.x \ cell_wid + 1: If (cell_row < 1) Or (cell_row > num_cells_row) Then Continue For
      Dim As Integer cell_col = Particle(i).Psn.y \ cell_hgt + 1: If (cell_col < 1) Or (cell_col > num_cells_col) Then Continue For
      
      With cell(cell_row, cell_col)
         
         If .num_particles = max_cell_particles Then Continue For
         .num_particles += 1
         .particle(.num_particles) = i
         
      End With
      
   Next
   
End Sub

Sub narrow_phase_collision()
   
   ''   reset number of neighbors
   num_neighbors = 0
   
   ''   loop through all cells and their neighbors to find particles within interaction radius
   For x As Integer = 1 To num_cells_row
      For y As Integer = 1 To num_cells_col
      	
      	With Cell(x, y)
      		
      		If .num_particles = 0 Then Continue For
      		
      		'' cell(x, y) - self
      		For i As Integer = 1 To .num_particles-1
            	
	            For l_ As Integer = i+1 To .num_particles
	            	
	               Check_if_neighbor(.particle(i), .particle(l_))
	               
	            Next
   				
      		Next
      		
      		'' cell(x, y) - neighboring cells
   			For i As Integer = 1 To .num_particles
   				
   				For j As Integer = 1 To .numCellNeighbors
   					
   					For k_ As Integer = 1 To .CellNeighbor(j)->num_particles
   						
   						Check_if_neighbor(.particle(i), .CellNeighbor(j)->Particle(k_))
   						
   					Next
   					
   				Next
   				
   			Next
   			
      	End With
      	
      Next
   Next
   
End Sub

Sub reset_cells()
	
   For x As Integer = 1 To num_cells_row
      For y As Integer = 1 To num_cells_col
         cell(x, y).Num_particles = 0
      Next
   Next
   
End Sub

Sub reset_particles()
	
   For i As Integer = 1 To num_particles
      With particle(i)
         .density       = 1
         .num_neighbors = 0
         .frc           = Vec2f(0.0f, G)
      End With
   Next
   
End Sub

Sub calculate_pressure()
	
   For i as Integer = 1 to num_particles
      With particle(i)
         .pressure = k * (.density - rest_density)
      End With
   Next
	
End Sub

Sub calculate_internal_force()
   
   For i As Integer = 1 To num_neighbors
         
      With neighbor(i)
         
         Dim as float F_pressure
         
         Dim as float F_viscosity
         
         Dim as float vel, W
         
         If .distance_squared > zero Then	
         
	         ''   velocity vector projected onto distance vector
	         vel = (particle(.a).vel - particle(.b).vel).dot(.n_over_distance)
	         
	         ''   viscosity force
	         F_viscosity = mju * (vel / (particle(.a).density * particle(.b).density))
	         
	         ''   pressure force
	         'F_pressure = (0.5 * (particle(.a).pressure + particle(.b).pressure)) / _
	         '   (particle(.a).density + particle(.b).density) 
	         
	         ''   pressure force
	         F_pressure = (particle(.a).pressure + particle(.b).pressure) / _
	            (particle(.a).density * particle(.a).density + particle(.b).density * particle(.b).density) 
	     		
         End If
         
         F_pressure += 0.5 * (particle(.a).density + particle(.b).density) * k_near
     		 
         ''   force vector
         Dim as Vec2f F = .n_over_distance * (F_pressure + F_viscosity)
         
         ''   apply opposite equal forces
         particle(.a).frc -= F
         particle(.b).frc += F
         
      End With
      
   Next
   
End Sub

Sub screen_border_force()
   
   For i As Integer = 1 To num_particles
      With particle(i)
      	
      	If .psn.x < r_plus_border Then 
            .frc.x -= (.psn.x - r_plus_border) * border_k + .vel.x * Border_D
      	ElseIf .psn.x > screen_wid - r_plus_border Then 
            .frc.x -= (.psn.x - (screen_wid - r_plus_border)) * border_k + .vel.x * Border_D
         EndIf
         
         If .psn.y < r_plus_border Then 
            .frc.y -= (.psn.y - r_plus_border) * border_k + .vel.y * Border_D
         ElseIf .psn.y > screen_hgt - r_plus_border Then 
            .frc.y -= (.psn.y - (screen_hgt - r_plus_border)) * border_k + .vel.y * Border_D
         EndIf
         
      End With
   Next
   
End sub

Sub Integrate()
   
   ''   symplectic Euler
   'For i As Integer = 1 To num_particles
   '   With particle(i)
   '      
   '      .vel += .frc * dt
   '      .psn += .vel * dt
   '      .frc = Vec2f(0.0f, G)
   '      
   '   End With
   'Next
   
   ''   midpoint / leapfrog
   For i As Integer = 1 To num_particles
      With particle(i)
         
         .vel += .frc * half_dt
         .psn += .vel * dt + 0.5 *.frc * dt_sqrd
         .vel += .frc * half_dt
         .frc = Vec2f(0.0f, G)
         
      End With
   Next
   
   ''   verlet
   'For i As Integer = 1 To num_particles
   '   With particle(i)
   '      
   '      Dim As vec2f temp = .psn
   '      .psn       += .psn - .psn_old + .frc * dt_Sqrd
   '      .psn_old   = temp
   '      .vel       = (.psn - .psn_old) * inv_dt
   '      
   '   End With
   'Next
   '
End Sub

Sub draw_particles()
   
   glClear(GL_DEPTH_BUFFER_BIT Or GL_COLOR_BUFFER_BIT)
   
   glLoadIdentity()
   
   glBegin(GL_POINTS)
   
      For i As Integer = 1 To num_Particles
      	
      	With Particle(i)
	      	
	      	Dim As float col = (.density-1) * 0.1
	      	
	      	glColor3f(col, 0.0, 1 - col)
	      	
	         glVertex2f(.psn.x, .psn.y)
	         
      	End With
         
      Next
      
   glEnd()
   
   glfwSwapBuffers()
   
End Sub

Function Check_if_neighbor(ByVal a As Integer, ByVal b As Integer) As Integer
   
   If num_neighbors = max_num_neighbors Then Return FALSE
   
   If Particle(a).num_neighbors = max_neighbors Then Return FALSE
   If Particle(b).num_neighbors = max_neighbors Then Return FALSE
   
   Dim As vec2f distance_vector = Particle(b).psn - Particle(a).psn
   
   If Abs(distance_vector.x) > r Then Return FALSE
   If Abs(distance_vector.y) > r Then Return FALSE
   
   Dim As float distance_squared = distance_vector.magnitudesquared()
   
   If (distance_squared > r_sqrd) Then Return FALSE
   
   ''	density kernel W(Xij, r) = ( 1 - Xij^2 / r^2 )^2
   Dim As float density = 1 - distance_squared * inv_r_sqrd: density *= density
   
	'' create neighbor
   num_neighbors += 1
   
   With neighbor(num_neighbors)
      .a = a
      .b = b
      .distance_squared = distance_squared
      
      If distance_squared > 1e-3 Then 
      	.n_over_distance = distance_vector / distance_squared
      Else
      	.n_over_distance = vec2f(0, 0)
      End If
      
   End With
   
   With Particle(a)
      .num_neighbors += 1
      .density += density
   End With
   
   With Particle(b)
      .num_neighbors += 1
      .density += density
   End With
   
End Function

Function ScreenType.CreateScreen(w As Integer, h As Integer, r As Integer, g As Integer, _
   b As Integer, a As Integer, d As Integer, s As Integer, m As Integer) As Integer
   
   ''   this function creates an GLFW opengl screen
   If glfwinit() Then
      
      Dim As GLFWvidmode glfwvm
      glfwGetDesktopMode(@glfwvm)
      
      If w = 0 Or h = 0 Then 
         w = glfwvm.width
         h = glfwvm.height
         m = GLFW_FULLSCREEN
      End If
      
      'glfwopenwindowhint(GLFW_FSAA_SAMPLES, 8)
      glfwOpenWindowHint(GLFW_WINDOW_NO_RESIZE, GL_TRUE)
      
      If glfwopenwindow(w, h, r, g, b, a, d, s, m) Then
         glfwSetWindowPos(glfwvm.width\2 - w\2, glfwvm.height\2 - h\2)
         glfwSetWindowTitle("Mike's square root free SPH example # 8.")
         
         glMatrixMode(GL_PROJECTION)  
         glLoadIdentity()
         glViewport(0, 0, w, h)
         glOrtho(0, w, h, 0, 0, 1)
         glMatrixMode(GL_MODELVIEW)
         glLoadIdentity()
         glEnable(GL_CLEAR)
         glClearColor(0.0, 0.0, 0.0, 0.1)
         glPointSize(3)
         glLineWidth(1)
         glColor3f(0.3, 0.5, 1.0)
         
         Return TRUE
         
      End If
      
      glfwTerminate()
      
   End If
   
   Return FALSE
   
End Function

Function ScreenType.DeleteScreen() As Integer
	
   glfwCloseWindow()
   glfwTerminate()
   Return TRUE
   
End Function

Sub MouseType.Update
   
   ''
   Psn_Old = Psn
   Lbt_Old = Lbt
   Rbt_Old = Rbt
   Whl_Old = Whl
   
   ''
   glfwGetMousePos(@Psn.x, @Psn.y)
   Lbt = glfwgetmousebutton(GLFW_MOUSE_BUTTON_LEFT)
   Rbt = glfwgetmousebutton(GLFW_MOUSE_BUTTON_RIGHT)
   whl = glfwgetmousewheel()
   
   ''
   vel.x = (psn.x-psn_old.x) * inv_dt
   vel.y = (psn.y-psn_old.y) * inv_dt
   
   If (Lbt = GLFW_PRESS) Then
	   
	   For i as Integer = 1 to num_particles
			
		   Dim As vec2f dst = Particle(i).psn - psn
		   
		   If Abs(dst.x) > interaction_radius then continue for
		   If Abs(dst.y) > interaction_radius then continue For
		   
		   Dim As float dst_Sqrd = dst.x*dst.x+dst.y*dst.y
		   
		   If dst_Sqrd > int_radius_sqrd then continue For
		   
		   If vel.x <> 0 Then Particle(i).frc.x -= ( (Particle(i).vel.x - vel.x) * inv_dt ) * ( 1 - dst_Sqrd / int_radius_sqrd ) * 0.15
		   If vel.y <> 0 Then Particle(i).frc.y -= ( (Particle(i).vel.y - vel.y) * inv_dt ) * ( 1 - dst_Sqrd / int_radius_sqrd ) * 0.15
		   
		   
	   Next
	   
   End If
   
   If whl > whl_old Then
   	
   	For i as Integer = 1 to num_particles
		
		   Dim As vec2f dst = Particle(i).psn - psn
		   
		   If Abs(dst.x) < 1 then continue for
		   If Abs(dst.y) < 1 then continue For
		   
		   Dim As float dst_Sqrd = dst.x*dst.x+dst.y*dst.y
		   
		   Particle(i).frc.x +=  dst.y * ( 1 / dst_Sqrd ) * 200000
		   Particle(i).frc.y -=  dst.x * ( 1 / dst_Sqrd ) * 200000
		    
	   Next
   	
   EndIf
   
   If whl < whl_old Then
   	
   	For i as Integer = 1 to num_particles
		
		   Dim As vec2f dst = Particle(i).psn - psn
		   
		   If Abs(dst.x) < 1 then continue for
		   If Abs(dst.y) < 1 then continue For
		   
		   Dim As float dst_Sqrd = dst.x*dst.x+dst.y*dst.y
		   
		   Particle(i).frc.x -=  dst.y * ( 1 / dst_Sqrd ) * 200000
		   Particle(i).frc.y +=  dst.x * ( 1 / dst_Sqrd ) * 200000
		    
	   Next
   	
   EndIf
   
End Sub

Cheers,
Mike

[Voltage]

I seriously love fluid dynamics.

I've read as many articles as I can over the last 12 months.... but I still don't have a stable SPH routine yet.. I'll let you offer your expertise if you like :)

I was getting an error when trying to compile with the latest version of FB... had to add this line to the GLFW.bi include file:

Code: Select all

#define GLFW_WINDOW_NO_RESIZE 131090

More! Oncore!