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