Learn Verlet Physics

BlitzMax Forums/BlitzMax Tutorials/Learn Verlet Physics

Nate the Great

(Posted 2009) [#1]

Hi

This is just a translated version of my b3d tutorial: sorry its not in OO but that is because I translated it almost directly from b3d

NOTE: I know this lacks timestep and collision detection etc but the point of this tutorial is to help people that don't get it, get it.

Here is a tutorial I decided to write for anyone who wants to understand verlet integration but really doesn't know how to decode all of those confusing equations on wikipedia.

So where to start? Well verlet integration is based around point masses. These point masses are bound together by constraints that can be rigid or flexible. Thats about it. Now to coding it.

The pointmass Type itself:

Each point mass has an x value and a y value as well as an ox (old x) and an oy (old y) These ox and oy values are used to calculate the velocity of the particle after each frame. So right now our code for the verlet type would look like this in b3d or bmax:

Type pointmass
	Field x#
	Field y#
	Field ox#
	Field oy#
End Type

now this doesn't do much by itself so we must come up with an update function. The update function figures out the velocities of each verlet and moves them one step further. This is easier than it sounds.

The equation of figuring out the velocity couldn't be simpler

dx = x - ox
dy = y - oy

dx is the difference between the ox and the current x

so now we put the current x and y values in our ox and oy fields like this

ox = x
oy = y

next we add our velocity to the current x and y to simulate a time step forward.

x = x + dx
y = y + dy

could it be any simpler? This satisfies newton's first law. Paraphrased: an object in motion stays in motion unless acted on by another force

So here is what our type and update function look like in bmax code

Global pointlist:TList = New TList

Type pointmass
	Field x#
	Field y#
	Field ox#
	Field oy#
End Type

Function updatepointmasses()

	For p:pointmass = EachIn pointlist
		dx# = p.x - p.ox
		dy# = p.y - p.oy
		
		p.ox = p.x
		p.oy = p.y
		
		p.x = p.x + dx
		p.y = p.y + dy
	Next

End Function

now to make a short demo of this, we must add drawing functions and a simple while loop.

Global pointlist:TList = New TList

Type pointmass
	Field x#
	Field y#
	Field ox#
	Field oy#
End Type

Function updatepointmasses()
	
	For p:pointmass = EachIn pointlist
		DrawOval p.x-4,p.y-4,8,8
		dx# = p.x - p.ox
		dy# = p.y - p.oy
		
		p.ox = p.x
		p.oy = p.y
		
		p.x = p.x + dx
		p.y = p.y + dy
	Next

End Function

Function createpointmass:pointmass(x#,y#,vx#,vy#)	' x And y are coords For the verlet. vx And vy are velocity values For the verlet

p:pointmass = New pointmass
p.x = x
p.y = y
p.ox = x-vx	'gives the particle a starting velocity
p.oy = y-vy
pointlist.addlast(p:pointmass)

Return p:pointmass

End Function

Graphics 640,480,0,60

createpointmass(100,100,5,1)

While Not KeyDown(key_escape)
Cls


updatepointmasses()
Flip
Wend
End

so the code above simply makes a particle and gives it a velocity. :-)

Now I will move on to constraints which I view as wooden rods between verlets. So how do constraints work? well constraints are satisfied in a loop by brute force pretty much. All other methods are either too ineffiecient or just not accurate enough. The constraint loop goes after the loop that updates velocity and works like this: if a constraint is too short due to the particles being compressed, then stretch it to make it the right size. If a constraint is too long because the particles were stretched apart then shrink it to make it the right size.

The constraint type will look something like this:

Type constraint
	Field p1:pointmass
	Field p2:pointmass
	Field length#
End Type

Here are the equations we will use for constraints (not optimized but for simplicity instead)

find the distance between the constrained verlets:
dist = sqr((x1-x2)^2 + (y1-y2)^2)

get the differences in the x and y values:
dx = x1-x2
dy = y1-y2

dist is the distance between 2 particles and length is the length of the constraint. diff is the difference between the two

diff = dist-length

diff = diff / length

dx = dx * .5
dy = dy * .5

and finally the constraint is solved with:

x1 = x1 - (diff * dx)
y1 = y1 - (diff * dy)

x2 = x2 + (diff *dx)
y2 = y2 + (diff * dy)

If you really want to know how this works it is best to draw it out on graph paper, and do it step by step.

One final note on constraints before you see the code is that the constraints are updated 5-20 times every loop depending on the accuracy needed for the engine. This is only necessary when there are many constraints but I will put it in anyway.

For now here is the update function we are going to use. this is not for efficiency, it is made so it is easy to understand.

Function updateconstraints()
	
	For cnt = 0 To CONST_ITERATIONS-1	'this is necessary with many constraints To solve them correctly
		For c:constraint = EachIn constraintlist
			dist# = Sqr((c.p1.x-c.p2.x)^2 + (c.p1.y-c.p2.y)^2)	'distance formula
			
			diff# = dist#-c.length#	'shows the margin of error the update loop has created so it can be corrected
			
			dx# = c.p1.x-c.p2.x	'difference between x's and y's
			dy# = c.p1.y-c.p2.y
			
			If c.length > 0 Then	'prevents a divided by 0 error that may occur
				diff = diff / c.length
			Else
				diff = 0
			EndIf
			
			dx = dx * .5
			dy = dy * .5
			
			c.p1.x = c.p1.x - (diff*dx)
			c.p1.y = c.p1.y - (diff*dy)
			
			c.p2.x = c.p2.x + (diff*dx)
			c.p2.y = c.p2.y + (diff*dy)
		Next
	Next
	
End Function

Now if we put this all together, we can make a basic demo to show that constraints work.

Global pointlist:TList = New TList
Global constraintlist:TList = New TList
Global CONST_ITERATIONS = 5

Type pointmass
	Field x#
	Field y#
	Field ox#
	Field oy#
End Type

Function updatepointmasses()
	
	For p:pointmass = EachIn pointlist
		DrawOval p.x-4,p.y-4,8,8
		dx# = p.x - p.ox
		dy# = p.y - p.oy
		
		p.ox = p.x
		p.oy = p.y
		
		p.x = p.x + dx
		p.y = p.y + dy
	Next

End Function

Function createpointmass:pointmass(x#,y#,vx#,vy#)	' x And y are coords For the verlet. vx And vy are velocity values For the verlet

p:pointmass = New pointmass
p.x = x
p.y = y
p.ox = x-vx	'gives the particle a starting velocity
p.oy = y-vy
pointlist.addlast(p:pointmass)

Return p:pointmass

End Function


Type constraint
	Field p1:pointmass
	Field p2:pointmass
	Field length#
End Type


Function updateconstraints()
	
	For cnt = 0 To CONST_ITERATIONS-1	'this is necessary with many constraints To solve them correctly
		For c:constraint = EachIn constraintlist
			dist# = Sqr((c.p1.x-c.p2.x)^2 + (c.p1.y-c.p2.y)^2)	'distance formula
			
			diff# = dist#-c.length#	'shows the margin of error the update loop has created so it can be corrected
			
			dx# = c.p1.x-c.p2.x	'difference between x's and y's
			dy# = c.p1.y-c.p2.y
			
			If c.length > 0 Then	'prevents a divided by 0 error that may occur
				diff = diff / c.length
			Else
				diff = 0
			EndIf
			
			dx = dx * .5
			dy = dy * .5
			
			c.p1.x = c.p1.x - (diff*dx)
			c.p1.y = c.p1.y - (diff*dy)
			
			c.p2.x = c.p2.x + (diff*dx)
			c.p2.y = c.p2.y + (diff*dy)
		Next
	Next
	
End Function

Function createconstraint:constraint(p1:pointmass,p2:pointmass)
	
	c:constraint = New constraint
	c.p1 = p1
	c.p2 = p2
	c.length = Sqr((p1.y-p2.y)^2+(p1.x-p2.x)^2)
	
	constraintlist.addlast(c:constraint)
	Return c:constraint
End Function


Function drawconstraints()
	For c:constraint = EachIn constraintlist
		DrawLine c.p1.x,c.p1.y,c.p2.x,c.p2.y
	Next
End Function


Graphics 640,480,0,60

p1:pointmass = createpointmass(100,100,1,0)
p2:pointmass = createpointmass(100,150,0,0)
c:constraint = createconstraint(p1,p2)

While Not KeyDown(key_escape)
Cls


updatepointmasses()
updateconstraints()
drawconstraints()

Flip
Wend
End

Ok, time for a little demo of a square falling using the same verlet engine in the above post. I just added some verlets and constraints before the loop

Global pointlist:TList = New TList
Global constraintlist:TList = New TList
Global CONST_ITERATIONS = 5

Type pointmass
	Field x#
	Field y#
	Field ox#
	Field oy#
End Type

Function updatepointmasses()
	
	For p:pointmass = EachIn pointlist
		DrawOval p.x-4,p.y-4,8,8
		dx# = p.x - p.ox
		dy# = p.y - p.oy + .1
		
		p.ox = p.x
		p.oy = p.y
		
		p.x = p.x + dx
		p.y = p.y + dy
		
		If p.y > 480 Then
			p.y = 480
			dx = p.x-p.ox
			p.ox = p.x-dx/2
		endif
	Next

End Function

Function createpointmass:pointmass(x#,y#,vx#,vy#)	' x And y are coords For the verlet. vx And vy are velocity values For the verlet

p:pointmass = New pointmass
p.x = x
p.y = y
p.ox = x-vx	'gives the particle a starting velocity
p.oy = y-vy
pointlist.addlast(p:pointmass)

Return p:pointmass

End Function


Type constraint
	Field p1:pointmass
	Field p2:pointmass
	Field length#
End Type


Function updateconstraints()
	
	For cnt = 0 To CONST_ITERATIONS-1	'this is necessary with many constraints To solve them correctly
		For c:constraint = EachIn constraintlist
			dist# = Sqr((c.p1.x-c.p2.x)^2 + (c.p1.y-c.p2.y)^2)	'distance formula
			
			diff# = dist#-c.length#	'shows the margin of error the update loop has created so it can be corrected
			
			dx# = c.p1.x-c.p2.x	'difference between x's and y's
			dy# = c.p1.y-c.p2.y
			
			If c.length > 0 Then	'prevents a divided by 0 error that may occur
				diff = diff / c.length
			Else
				diff = 0
			EndIf
			
			dx = dx * .5
			dy = dy * .5
			
			c.p1.x = c.p1.x - (diff*dx)
			c.p1.y = c.p1.y - (diff*dy)
			
			c.p2.x = c.p2.x + (diff*dx)
			c.p2.y = c.p2.y + (diff*dy)
		Next
	Next
	
End Function

Function createconstraint:constraint(p1:pointmass,p2:pointmass)
	
	c:constraint = New constraint
	c.p1 = p1
	c.p2 = p2
	c.length = Sqr((p1.y-p2.y)^2+(p1.x-p2.x)^2)
	
	constraintlist.addlast(c:constraint)
	Return c:constraint
End Function


Function drawconstraints()
	For c:constraint = EachIn constraintlist
		DrawLine c.p1.x,c.p1.y,c.p2.x,c.p2.y
	Next
End Function


Graphics 640,480,0,60

p1:pointmass = createpointmass(100,100,0,0)
p2:pointmass = createpointmass(100,150,0,0)
p3:pointmass = createpointmass(150,150,0,0)
p4:pointmass = createpointmass(150,100,0,0)

c:constraint = createconstraint(p1,p2)
c1:constraint = createconstraint(p2,p3)
c2:constraint = createconstraint(p3,p4)
c3:constraint = createconstraint(p1,p4)


While Not KeyDown(key_escape)
Cls


updatepointmasses()
updateconstraints()
drawconstraints()

Flip
Wend
End

so run that and see what happens... uh It falls over! so we must support it by adding the following lines to the beginning after all other constraints are created.

c.constraint = createconstraint(p1,p3)
c.constraint = createconstraint(p2,p4)

so the code will look like this:

Global pointlist:TList = New TList
Global constraintlist:TList = New TList
Global CONST_ITERATIONS = 5

Type pointmass
	Field x#
	Field y#
	Field ox#
	Field oy#
End Type

Function updatepointmasses()
	
	For p:pointmass = EachIn pointlist
		DrawOval p.x-4,p.y-4,8,8
		dx# = p.x - p.ox
		dy# = p.y - p.oy + .1
		
		p.ox = p.x
		p.oy = p.y
		
		p.x = p.x + dx
		p.y = p.y + dy
		
		If p.y > 480 Then
			p.y = 480
			dx = p.x-p.ox
			p.ox = p.x-dx/2
		endif
	Next

End Function

Function createpointmass:pointmass(x#,y#,vx#,vy#)	' x And y are coords For the verlet. vx And vy are velocity values For the verlet

p:pointmass = New pointmass
p.x = x
p.y = y
p.ox = x-vx	'gives the particle a starting velocity
p.oy = y-vy
pointlist.addlast(p:pointmass)

Return p:pointmass

End Function


Type constraint
	Field p1:pointmass
	Field p2:pointmass
	Field length#
End Type


Function updateconstraints()
	
	For cnt = 0 To CONST_ITERATIONS-1	'this is necessary with many constraints To solve them correctly
		For c:constraint = EachIn constraintlist
			dist# = Sqr((c.p1.x-c.p2.x)^2 + (c.p1.y-c.p2.y)^2)	'distance formula
			
			diff# = dist#-c.length#	'shows the margin of error the update loop has created so it can be corrected
			
			dx# = c.p1.x-c.p2.x	'difference between x's and y's
			dy# = c.p1.y-c.p2.y
			
			If c.length > 0 Then	'prevents a divided by 0 error that may occur
				diff = diff / c.length
			Else
				diff = 0
			EndIf
			
			dx = dx * .5
			dy = dy * .5
			
			c.p1.x = c.p1.x - (diff*dx)
			c.p1.y = c.p1.y - (diff*dy)
			
			c.p2.x = c.p2.x + (diff*dx)
			c.p2.y = c.p2.y + (diff*dy)
		Next
	Next
	
End Function

Function createconstraint:constraint(p1:pointmass,p2:pointmass)
	
	c:constraint = New constraint
	c.p1 = p1
	c.p2 = p2
	c.length = Sqr((p1.y-p2.y)^2+(p1.x-p2.x)^2)
	
	constraintlist.addlast(c:constraint)
	Return c:constraint
End Function


Function drawconstraints()
	For c:constraint = EachIn constraintlist
		DrawLine c.p1.x,c.p1.y,c.p2.x,c.p2.y
	Next
End Function


Graphics 640,480,0,60

p1:pointmass = createpointmass(100,100,0,0)
p2:pointmass = createpointmass(100,150,0,0)
p3:pointmass = createpointmass(150,150,0,0)
p4:pointmass = createpointmass(150,100,0,0)

c:constraint = createconstraint(p1,p2)
c1:constraint = createconstraint(p2,p3)
c2:constraint = createconstraint(p3,p4)
c3:constraint = createconstraint(p1,p4)

c5:constraint = createconstraint(p2,p4)
c6:constraint = createconstraint(p1,p3)


While Not KeyDown(key_escape)
Cls


updatepointmasses()
updateconstraints()
drawconstraints()

Flip
Wend
End

Now the only thing that needs to be done is a collision engine... the hard part (for me at least)

BlitzSupport

(Posted 2009) [#2]

This looks like good stuff... I'll be sure to read it at some point. Thanks!

ImaginaryHuman

(Posted 2009) [#3]

Hey, this is nice. I was able to follow you all the way through and can now say that I understand how this works! Thanks.

I don't entirely grasp *why* it's necessary to do each of the steps in the constaint processing, but I'm sure it'll become clearer.

How exactly would someone go about collision detection with this? I guess if you just have shapes falling down the screen you could say that if a pointmass's y coordinate goes past the bottom of the screen you fix it to the bottom, but doesn't it need to calculate some kind of rebound force or something?

Also how easy it is to adapt this to using springs instead of fixed constraints?

Nate the Great

(Posted 2009) [#4]

ok nice to know it helped some people! I love teaching.

it's necessary to do the constraint loop many times because when you have many constraints connected to many verlets, one time through wont solve it, it will get closer every loop it does. - hope that clears it up a little?

as for the springs, im not really sure if this is the way you are supposed to do it but there are two things I do (and they work perfectly so I dont see whats wrong with it):

1. I give each constraint an elasticity variable (the value should be 1 or more... any less and it will explode.)
2. in the constraint loop, change this:

c.p1.x = c.p1.x - (diff*dx)
c.p1.y = c.p1.y - (diff*dy)

c.p2.x = c.p2.x + (diff*dx)
c.p2.y = c.p2.y + (diff*dy)

to this:

c.p1.x = c.p1.x - (diff*dx)/c.elasticity
c.p1.y = c.p1.y - (diff*dy)/c.elasticity

c.p2.x = c.p2.x + (diff*dx)/c.elasticity
c.p2.y = c.p2.y + (diff*dy)/c.elasticity

and make sure you set the elasticity to at least 1 where the constraints are created so the system doesnt explode.

Nate the Great

(Posted 2009) [#5]

oops I forgot to answer your question about collision detection

I will post an official tutorial on this sometime later but I don't understand it enough now. Here is what I am experimenting with right now:

loop through all verlets and test them against all other verlets for collision because verlets in my engine have a radius and can collide with eachother.

test for collisions between all constraints and all verlets

to cause them to collide just change the x and y values to a legal value and position the ox and oy values accordingly

for some good code for basing your collisions off of, here is my math.bmx file

'used for line intersect functions

Global Intersection_X#
Global Intersection_Y#
Global Intersection_AB#
Global Intersection_CD#

Function acqrdistance#[](x:Float,y:Float,x2:Float,y2:Float,acqdist:Float)	'x = col point x2 = center of collided circle ditto for y's  acqdist = the distance to acquire between these two points
	
	Local temparay#[]
	Local nx:Float
	Local ny:Float
	Local dx! = x - x2	'doubles because these are the cause of most floating point inacuracy that causes jerking and bouncing around
	Local dy!= y - y2
	Local dist:Float = Sqr(dx*dx + dy*dy)
	
	If dist < acqdist And dist > 0 Then
		
		dx = (dx / dist) * acqdist
		dy = (dy / dist) * acqdist
		nx = x - dx
		ny = y - dy
		
	Else
		nx = x2
		ny = y2
	EndIf
	
	temparay# = [nx,ny]
	
	Return temparay#
	
End Function



Function acqrdistance2#[](x:Float,y:Float,x2:Float,y2:Float,acqdist:Float)	'x = col point x2 = center of collided circle ditto for y's  acqdist = the distance to acquire between these two points
	
	Local temparay#[]
	Local nx:Float
	Local ny:Float
	Local dx! = x - x2	'doubles because these are the cause of most floating point inacuracy that causes jerking and bouncing around
	Local dy!= y - y2
	Local dist:Float = Sqr(dx*dx + dy*dy)
	
	If dist > 0 Then	'took out dist < acqdist
		
		dx = (dx / dist) * acqdist
		dy = (dy / dist) * acqdist
		nx = x + dx
		ny = y + dy
		
	Else
		nx = x2
		ny = y2
	EndIf
	
	temparay# = [nx,ny]
	
	Return temparay#
	
End Function




Function linecircleintersect#[][]( pLineStart#[], pLineEnd#[], pCircleCenter#[], pCircleRadius# )	'many thanks to blitz support for this one!
	
	Local tmpIntersections#[][]
	
	Local p# = pCircleCenter[0], q# = pCircleCenter[1]
	Local m# = (pLineEnd[1]-pLineStart[1])/(pLineEnd[0]-pLineStart[0])
	Local r# = pCircleRadius
	Local t# = pLineEnd[1]- (m*pLineEnd[0])
	Local s# = t-q
	
	Local a# = m*m + 1, b# = (2*m*s) - (2*p), c# = s*s + p*p - (r*r)
	
	Local bsqminfourac# = b*b-4*a*c
	
	If bsqminfourac > 0 Then
		
		bsqminfourac = Sqr(bsqminfourac)
		
		Local x1# = ((-b)+bsqminfourac)/(2*a)
		Local x2# = ((-b)-bsqminfourac)/(2*a)
		
		tmpIntersections = [[x1,(m*x1)+t],[x2,(m*x2)+t]]
		
	ElseIf bsqminfourac = 0 Then
		
		tmpIntersections = [[(-b)/(2*a),(-b*m)/(2*a)+t]]
		
	EndIf
	
	Return tmpIntersections

EndFunction


Function Lines_Intersect:Int(Ax#, Ay#, Bx#, By#, Cx#, Cy#, Dx#, Dy#)	'finds the line intersections... never mind commenting the math.. I dont even understand it anymore
	Local dymcy# = dy - cy
	Local dxmcx# = dx - cx
	
	Local Rn# = (Ay#-Cy#)*(DxmCx#) - (Ax#-Cx#)*(DymCy#)
	Local Rd# = (Bx#-Ax#)*(DymCy#) - (By#-Ay#)*(DxmCx#)
		
	If Rd# = 0 
		
		' Lines are parralel.
		' If Rn# is also 0 Then lines are coincident.  All points intersect. 
		' Otherwise, there is no intersection point.
		Return False
	Else
		' The lines intersect at some point.  Calculate the intersection point.
	Local Sn# = (Ay#-Cy#)*(Bx#-Ax#) - (Ax#-Cx#)*(By#-Ay#)
		Intersection_AB# = Rn# / Rd#
		Intersection_CD# = Sn# / Rd#
		Intersection_X# = Ax# + Intersection_AB#*(Bx#-Ax#)
		Intersection_Y# = Ay# + Intersection_AB#*(By#-Ay#)	
	EndIf
		
	If Intersection_AB#>0 And Intersection_AB#<1 And Intersection_CD#>0 And Intersection_cd#<1
			Return True	
		Else
			Return False
	EndIf


End Function


Function UpdateVerletPolysegCollisions()

Local cnt:Int = 0
Local v:verlet

For cnt = 0 To Verlcount-1
	V:verlet = Verletarray:Verlet[cnt]
	If V.active = True Then
		
		
		checkverletbounds v:verlet
		
		If v.x < activerightbound And v.x > activeleftbound And v.y > activeupbound And v.y < activelowbound Then
			
			If v.ghost = False Then updatepolysegs(v:verlet)
		
		EndIf
		
	EndIf
	
Next

End Function





'|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
'thanks for the following two functions "another JIM"
'Finds out if a point is on the left or right side of a segment
Function LineIntersectOX:Int(lx1:Float, ly1:Float, lx2:Float, ly2:Float, px:Float, py:Float)
	Local newpx:Float
	
	If ((ly1 < py) And (ly2 > py) ..
		Or (ly1 > py) And (ly2 < py)) And (lx1 > px) And (lx2 > px)
		Return True
	End If
	
	Local b:Int = ((ly1 <= py) And (ly2 >= py)) Or ((ly1 >= py) And (ly2 <= py))
	If (b = False)
		Return False
	EndIf
	
	If ly2 = ly1 Then
		Return False
	Else	
		newpx = ((py - ly1) / (ly2 - ly1)) * (lx2 - lx1) + lx1
		If (px < newpx) Return True
	EndIf
	Return False
End Function

'finds out if a point is inside a polygon
Function PointInPoly:Int(px:Float, py:Float, poly:Float[])
	Local i:Int
	Local sum:Int

	For i = 0 To Int(poly.length / 2) - 2
		Local p1x:Float = poly[i * 2]
		Local p1y:Float = poly[i * 2 + 1]
		Local p2x:Float = poly[(i + 1) * 2]
		Local p2y:Float = poly[(i + 1) * 2 + 1]
		
		sum:+LineIntersectOX(p1x, p1y, p2x, p2y, px, py)
	Next
	
	If (poly.length > 4)
		Local p1xd:Float = poly[0]
		Local p1yd:Float = poly[1]
		Local p2xd:Float = poly[poly.length - 2]
		Local p2yd:Float = poly[poly.length - 1]
		
		sum:+LineIntersectOX(p1xd, p1yd, p2xd, p2yd, px, py)
	End If
	
	Return sum Mod 2
End Function
'|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

sorry its not commented yet

ImaginaryHuman

(Posted 2009) [#6]

Great, the springs thing looks really simple to add.

So for collision detection are you saying that you just detect when a point or edge is inside/overlapping one from another object, and move the points until they don't?

Nate the Great

(Posted 2009) [#7]

yeah for the collision detection I do that but it greatly simplifies things if you make each point mass have a radius with only pointmass-pointmass collisions rather than constraint-pointmass collisions. My engine will have both as it seems that many modern engines like box2d, chipmunk, and farseer have poly-poly collison so if I want to compete I must include it too.

You can also use seperating axis theory to give a decent speedup but it unfortunately only works for convex objects.

ImaginaryHuman

(Posted 2009) [#8]

What about if the points are not very close to each other? Wouldn't that make the tip of a box penetrate another?

Nate the Great

(Posted 2009) [#9]

You can always approximate a box shape with no penetration... take a look at the verlet max demo in my sig. but I recommend making verlet-polygon collisions as they are almost essential for any general purpose engine

Fry Crayola

(Posted 2009) [#10]

Wow, Nate, this is good stuff. That sample is instantly impressive.

Are you still planning on producing a collision detection tutorial?

Nate the Great

(Posted 2009) [#11]

thanks

Are you still planning on producing a collision detection tutorial?

yes, when I figure it out! and get it implemented in my own physics engine...

*sigh* its really confusing because there are so many choices of what to do for collision detection, there is no set standard formula for collision detection with verlet physics..

_Skully

(Posted 2009) [#12]

Damn you Nate! Your hooking me into this stuff again! ;)

anyway...

Here is a version with the need for sqr removed... faster :)

Oh, and I rotated the box for more fun lol

Global pointlist:TList = New TList
Global constraintlist:TList = New TList
Global CONST_ITERATIONS = 5

Type pointmass
	Field x#
	Field y#
	Field ox#
	Field oy#
End Type

Function updatepointmasses()
	
	For p:pointmass = EachIn pointlist
		DrawOval p.x-4,p.y-4,8,8
		dx# = p.x - p.ox
		dy# = p.y - p.oy + .1
		
		p.ox = p.x
		p.oy = p.y
		
		p.x = p.x + dx
		p.y = p.y + dy
		
		If p.y > 480 Then
			p.y = 480
			dx = p.x-p.ox
			p.ox = p.x-dx/2
		endif
	Next

End Function

Function createpointmass:pointmass(x#,y#,vx#,vy#)	' x And y are coords For the verlet. vx And vy are velocity values For the verlet

p:pointmass = New pointmass
p.x = x
p.y = y
p.ox = x-vx	'gives the particle a starting velocity
p.oy = y-vy
pointlist.addlast(p:pointmass)

Return p:pointmass

End Function


Type constraint
	Field p1:pointmass
	Field p2:pointmass
	Field length#
End Type


Function updateconstraints()
	
	For cnt = 0 To CONST_ITERATIONS-1	'this is necessary with many constraints To solve them correctly
		For c:constraint = EachIn constraintlist
			dist# = ((c.p1.x-c.p2.x)^2 + (c.p1.y-c.p2.y)^2)	'distance formula
			
			diff# = dist#-c.length#	'shows the margin of error the update loop has created so it can be corrected
			
			dx# = c.p1.x-c.p2.x	'difference between x's and y's
			dy# = c.p1.y-c.p2.y
			
			If c.length > 0 Then	'prevents a divided by 0 error that may occur
				diff = diff / c.length
			Else
				diff = 0
			EndIf
			
			dx = dx * .25
			dy = dy * .25
			
			c.p1.x = c.p1.x - (diff*dx)
			c.p1.y = c.p1.y - (diff*dy)
			
			c.p2.x = c.p2.x + (diff*dx)
			c.p2.y = c.p2.y + (diff*dy)
		Next
	Next
	
End Function

Function createconstraint:constraint(p1:pointmass,p2:pointmass)
	
	c:constraint = New constraint
	c.p1 = p1
	c.p2 = p2
	c.length = ((p1.y-p2.y)^2+(p1.x-p2.x)^2)
	
	constraintlist.addlast(c:constraint)
	Return c:constraint
End Function


Function drawconstraints()
	For c:constraint = EachIn constraintlist
		DrawLine c.p1.x,c.p1.y,c.p2.x,c.p2.y
	Next
End Function


Graphics 640,480,0,60

p1:pointmass = createpointmass(125,100,0,0)
p2:pointmass = createpointmass(150,125,0,0)
p3:pointmass = createpointmass(125,150,0,0)
p4:pointmass = createpointmass(100,125,0,0)

c:constraint = createconstraint(p1,p2)
c1:constraint = createconstraint(p2,p3)
c2:constraint = createconstraint(p3,p4)
c3:constraint = createconstraint(p1,p4)

c5:constraint = createconstraint(p2,p4)
c6:constraint = createconstraint(p1,p3)


While Not KeyDown(key_escape)
Cls


updatepointmasses()
updateconstraints()
drawconstraints()

Flip
Wend
End

_Skully

(Posted 2009) [#13]

This is funner...

Space turns the slow drawline on/off
Left Control adds random vertical velocities to the points
Oh, and I made it a little bouncier and replaced drawoval with an image...

Global pointlist:TList = New TList
Global constraintlist:TList = New TList
Global CONST_ITERATIONS = 5

Type pointmass
	Field x#
	Field y#
	Field ox#
	Field oy#
End Type

Function updatepointmasses()
	
	For p:pointmass = EachIn pointlist
		DrawImage im,p.x,p.y
		dx# = p.x - p.ox
		dy# = p.y - p.oy + .1
		
		p.ox = p.x
		p.oy = p.y
		
		p.x = p.x + dx
		p.y = p.y + dy
		
		If p.y > 480 Then
			p.y = 480
			dx = p.x-p.ox
			p.ox = p.x-dx/2
		endif
	Next

End Function

Function createpointmass:pointmass(x#,y#,vx#,vy#)	' x And y are coords For the verlet. vx And vy are velocity values For the verlet

p:pointmass = New pointmass
p.x = x
p.y = y
p.ox = x-vx	'gives the particle a starting velocity
p.oy = y-vy
pointlist.addlast(p:pointmass)

Return p:pointmass

End Function


Type constraint
	Field p1:pointmass
	Field p2:pointmass
	Field length#
End Type


Function updateconstraints()
	
	For cnt = 0 To CONST_ITERATIONS-1	'this is necessary with many constraints To solve them correctly
		For c:constraint = EachIn constraintlist
			dist# = ((c.p1.x-c.p2.x)^2 + (c.p1.y-c.p2.y)^2)	'distance formula
			
			diff# = dist#-c.length#	'shows the margin of error the update loop has created so it can be corrected
			
			dx# = c.p1.x-c.p2.x	'difference between x's and y's
			dy# = c.p1.y-c.p2.y
			
			If c.length > 0 Then	'prevents a divided by 0 error that may occur
				diff = diff / c.length
			Else
				diff = 0
			EndIf
			
			dx = dx * .15
			dy = dy * .15
			
			c.p1.x = c.p1.x - (diff*dx)
			c.p1.y = c.p1.y - (diff*dy)
			
			c.p2.x = c.p2.x + (diff*dx)
			c.p2.y = c.p2.y + (diff*dy)
		Next
	Next
	
End Function

Function createconstraint:constraint(p1:pointmass,p2:pointmass)
	
	c:constraint = New constraint
	c.p1 = p1
	c.p2 = p2
	c.length = ((p1.y-p2.y)^2+(p1.x-p2.x)^2)
	
	constraintlist.addlast(c:constraint)
	Return c:constraint
End Function


Function drawconstraints()
	For c:constraint = EachIn constraintlist
		DrawLine c.p1.x,c.p1.y,c.p2.x,c.p2.y
	Next
End Function


Graphics 640,480,0,60
AutoMidHandle(True)
Global drawlines=False
Global im:TImage=CreateImage(8,8,1,DYNAMICIMAGE|MASKEDIMAGE)
DrawOval 0,0,8,8
GrabImage(im,0,0)
Cls

For ct=0 To 20
x=Rnd(0,590)
y=Rnd(0,100)
p1:pointmass = createpointmass(25+x,y,0,0)
p2:pointmass = createpointmass(50+x,25+y,0,0)
p3:pointmass = createpointmass(25+x,50+y,0,0)
p4:pointmass = createpointmass(x,25+y,0,0)

c:constraint = createconstraint(p1,p2)
c1:constraint = createconstraint(p2,p3)
c2:constraint = createconstraint(p3,p4)
c3:constraint = createconstraint(p1,p4)

c5:constraint = createconstraint(p2,p4)
c6:constraint = createconstraint(p1,p3)
Next

While Not KeyDown(key_escape)
Cls


updatepointmasses()
updateconstraints()
If drawlines Then drawconstraints()
If KeyHit(KEY_SPACE ) Then drawlines=1-drawlines
If KeyHit(KEY_LCONTROL)
	For p:pointmass=EachIn pointlist
		p.y=p.y-Rnd(0,10)
	Next
End If
Flip
Wend
End

Nate the Great

(Posted 2009) [#14]

Damn you Nate! Your hooking me into this stuff again! ;)

anyway...

Here is a version with the need for sqr removed... faster :)

Oh, and I rotated the box for more fun lol

lol nobody made you read it ;)

and thanks for the examples

edit
uhhhh skully? did you just take the sqr off of there? lol! idk how that works...

edit2: what exactly did you change to get rid of the sqr?

Fry Crayola

(Posted 2009) [#15]

Looks to me like he made the constraint's length variable the square of the actual length, meaning the square root is never required when you're processing the constraints later.

Pretty common trick for speeding things up if you only need to compare values.

_Skully

(Posted 2009) [#16]

exactly... since the comparative values are related by SQR... factor it out!

The other thing I had to do is:

dx = dx * .25
dy = dy * .25

so it didn't self destruct (which is 0.5^2)

Nate the Great

(Posted 2009) [#17]

hmm that just doesnt rest well with me... :/ sqr is not a linear operation like you cant just facor it out... well ill experiment though...

_Skully

(Posted 2009) [#18]

Sqr(4)/Sqr(4)=1
4/4=1

hmmm

If Sqr(9)<Sqr(4) then
If 3<2 then

Nate the Great

(Posted 2009) [#19]

sqr(4)/sqr(2) = irrational number thats not 2 lol
4/2 = 2

sure this works

If Sqr(9)<Sqr(4) then
If 3<2 then

but then you must take the sqr eventually because you are taking a difference in lengths and dividing dx by that dist or something... hmmmm

sqr(25) = distance
4 and 3 are the x and y components of said vector
the following makes said vector length of 1 by dividing the vector x and y by the vector length

4/sqr(25) = 4/5
3/sqr(25) = 3/5

now for without the sqr, your way I belive

4/(25*.25) = 4/6.25 close but not quite right
3/(25*.25) = 3/6.25 again, close but no cigar :/
the bigger the constraint lengths are the worse the error and more unstable it is

_Skully

(Posted 2009) [#20]

hmmm...

Ya maybe I just got lucky there... or perhaps there's a balancing act to be done here between exactness and speed which is quite often done... but there's no denying that it works as it is, but the question is... would that come back to haunt you later?

Zakk	(Posted 2010) [#21]

Demon magic, that's what it is.

Good tutorial :)