Verlet Physics Project

BlitzMax Forums/BlitzMax Beginners Area/Verlet Physics Project

Chroma(Posted 2011) [#1]
I've started doing research on verlet integration on this forum and on google etc. I'm going to take a stab at a very simple BMax verlet module, MaxVerlet...that will eventually find a home with Monkey.

If anyone wants to chime in that would be great. I've been a big supporter of Euler physics method for many many years. But since I've been checking into verlets, I see that it is MUCH more simple to get fantasic results.

So far I've gathered that the verlet system is comprised of Points which are Linked together. Some people call the points things such as particles, atoms, vertices etc etc...but here I'll refer to the point type as VPoint (V for verlet), which is the basic building block of a verlet object. Then of course the Points have to be connected together. I'll call this type VLink (also known as constraints, sticks etc etc).

So the initial VPoint type would look like this:
Type VPoint
   Field x#, y#
   Field oldx#, oldy#
End Type


And then VLink type would initially look like this:
Type VLink
   Field pointA:VPoint
   Field pointB:VPoint
End Type


The basic update for the VPoints looks like this:
Local tempx# = x
Local tempy# = y
x :+ x - oldx
y :+ y - oldy
oldx = tempx
oldy = tempy


If you decide to post here, please explain as much as you can about anything you add. This is here in the Beginner section for anyone who wants to learn this and create a light robust verlet system in the process.

Last edited 2011


Wiebo(Posted 2011) [#2]
Hi Chroma,
Good to see this happening, I've been using verlet integration for lots of things and while being simple, it indeed gives great results.

Something that is very handy to have is the ability to add weight to a VPoint, so some points can pull harder than other points.
So if you add a weight field to the VPoint class (it will be Class, not Type, in Monkey) and let the VLink update method take this into consideration by clever wizardry then that would be cool. The default value of the weight should then be 1.0 in the VPoint.New() method.

Also, it would be good to have a UpdateVerletLinks() function to update all links in one call. And maybe points as well?

Last edited 2011


Chroma(Posted 2011) [#3]
Cool thanks Wiebo. I've seen some verlet demos that are amazing.

So the New method you're talking about. I've never done it, I just use a Create function. Is that correct below? Makes sense that a point would need a mass. Now to research more about how to implement it.


VPoint: soon to be Class! Added the mass# field, New method and Update method.
Type VPoint
   Field x#, y#
   Field oldx#, oldy#
   Field mass#

   Method New()
      Self.mass = 1.0
   End Method

   Method Update()
      Local tempx# = Self.x
      Local tempy# = Self.y
      Self.x :+ Self.x - Self.oldx
      Self.y :+ Self.y - Self.oldy
      Self.oldx = tempx
      Self.oldy = tempy
   End Method
End Type



VLink: added length and elasticity field. Definitely need to know the what the correct length between pointA and pointB should be. And be able to adjust how rigid the springiness is. And put in a Create function.
Type VLink
   Field pointA:VPoint
   Field pointB:VPoint
   Field length#
   Field elasticity#

   Function Create:VLink(a:VPoint, b:VPoint, elasticity#)
      Local link:VLink = New VLink
      link.pointA = a
      link.pointB = b
      Local distx# = a.x - b.x
      Local disty# = a.y - b.y
      link.length = Sqr(distx * distx + disty * disty)
      link.elasticity = elasticity
      Return link
   End Function
End Type


Last edited 2011


Wiebo(Posted 2011) [#4]
I use the New() method for constructor stuff I always want executed, because extensions of the type will always run the super.new() as well. Using a Create function for setting defaults works, but you might forget it some day =] I always do. hehe. I recall reading that the New() method on Monkey will also accept passed values (and defaults for that) so the create Method might not even be needed later on!

The mass, length and elasticity (good call) fields might also need a getter and setter later on.

Last edited 2011


GW(Posted 2011) [#5]
Here is some code by Andreas from back in the day. Its very clean and clear.
I just converted it to Bmax.

'------------------------------------------
'// SIMPLE VERLET: by Andreas Blixt 2004 //
'------------------------------------------ 


Global PHY_GRAVITY#		= 100.0 ' Gravity decides how much the particles are pulled towards the ground.
Global PHY_TIMESTEP#	= 0.05 ' The timestep is how much the physics engine will progress per PHY_Update.
Global PHY_ITERATIONS%	= 5 ' More iterations gives more accuracy While less iterations give more speed.

' Constraints are the lines between the particles that keep them in a certain length from eachother.
Type PHY_Constraint
	Field p1:PHY_Particle
	Field p2:PHY_Particle
	Field length#
	Global List:TList
	
	Method New()
		If Not List Then list = CreateList()
		List.addlast(Self)
	End Method
End Type

' Particles are points that move based on their old positions, gravity And mass.
Type PHY_Particle
	Field x#,y#
	Field ox#,oy#
	Field mass#
	Global List:TList
	
	
	
	Method New()
		If Not List Then list = CreateList()
		List.addlast(Self)
	End Method

End Type

' This Function will make sure that all constraints Try To keep their length. It will be better the higher PHY_ITERATIONS is.
Function PHY_SatisfyConstraints()
	Local c:PHY_Constraint,p1:PHY_Particle,p2:PHY_Particle
	Local m1#,m2#,dx#,dy#,deltalength#,diff#
	Local i%

	' Do everything PHY_ITERATIONS times.
	For i = 1 To PHY_ITERATIONS
		For c = EachIn PHY_Constraint.list
			p1 = c.p1
			p2 = c.p2

			' Get difference in distance between the two particles.
			dx = p2.x - p1.x
			dy = p2.y - p1.y

			' Get the masses of the particles.
			m1 = -p1.mass
			m2 = -p2.mass

			' Get the length between the particles.
			deltalength = Sqr(dx * dx + dy * dy)
			If deltalength <> 0.0 ' Avoid division by 0
				' Get a factor of how much To move the particles from eachother.
				diff = (deltalength - c.length) / (deltalength * (m1 + m2))

				' Move the particles.
				p1.x = p1.x + dx * m1 * diff
				p1.y = p1.y + dy * m1 * diff
				p2.x = p2.x - dx * m2 * diff
				p2.y = p2.y - dy * m2 * diff
			EndIf
		Next
	Next
End Function

' This Function will move the particles.
Function PHY_Verlet()
	Local p:PHY_Particle
	Local x#,y#,dx#,dy#

	For p = EachIn PHY_Particle.list
		' Save Current position.
		x = p.x
		y = p.y

		' Get the difference between the old And New positions.
		dx = x - p.ox
		dy = y - p.oy

		' Move the particle based on the difference.
		p.x = p.x + dx
		If p.mass > 0.0 ' Avoid division by zero
			p.y = p.y + dy + PHY_GRAVITY / p.mass * PHY_TIMESTEP * PHY_TIMESTEP
		Else
			p.y = p.y + dy
		EndIf

		' Keep the previous position.
		p.ox = x
		p.oy = y
	Next
End Function










'// EXAMPLE //



AppTitle= "Verlet basics"

Graphics 800,600

' Create all particles (Normally this should be done in a Function To avoid so much creation code.)
Global p1:PHY_Particle = New PHY_Particle
p1.x = 368
p1.y = 268
p1.ox = 368
p1.oy = 268
p1.mass = 1.0

Global p2:PHY_Particle = New PHY_Particle
p2.x = 432
p2.y = 268
p2.ox = 432
p2.oy = 268
p2.mass = 1.0

Global p3:PHY_Particle = New PHY_Particle
p3.x = 368
p3.y = 332
p3.ox = 368
p3.oy = 332
p3.mass = 1.0

Global p4:PHY_Particle = New PHY_Particle
p4.x = 432
p4.y = 332
p4.ox = 432
p4.oy = 332
p4.mass = 1.0

' Create all constraints.
Global c:PHY_Constraint = New PHY_Constraint
c.p1 = p1
c.p2 = p2
c.length = 64

c = New PHY_Constraint
c.p1 = p1
c.p2 = p3
c.length = 64

c= New PHY_Constraint
c.p1 = p2
c.p2 = p4
c.length = 64

c = New PHY_Constraint
c.p1 = p3
c.p2 = p4
c.length = 63

' Diagonal lengths are Sqr(64 * 64 + 64 * 64)
c = New PHY_Constraint
c.p1 = p1
c.p2 = p4
c.length = 90.5097

c = New PHY_Constraint
c.p1 = p2
c.p2 = p3
c.length = 90.5097


While Not KeyHit(key_escape)
	Cls

	' Let the user control one of the particles with the arrow keys.
	If KeyDown(key_up) Then p1.y = p1.y - 5
	If KeyDown(key_down) Then p1.y = p1.y + 5
	If KeyDown(key_left) Then p1.x = p1.x - 5
	If KeyDown(key_right) Then p1.x = p1.x + 5

	' Apply forces To all the particles.
	PHY_Verlet

	' Make sure no particle leaves the screen.
	For Local p:PHY_Particle = EachIn PHY_Particle.list
		If p.x < 10 Then p.x = 10 ElseIf p.x > 790 Then p.x = 790
		If p.y < 10 Then p.y = 10 ElseIf p.y > 590 Then p.y = 590
	Next

	' Enforce lengths of constraints.
	PHY_SatisfyConstraints

	' Draw all the constraints.
	For c:PHY_Constraint = EachIn PHY_Constraint.list
		DrawLine c.p1.x,c.p1.y,c.p2.x,c.p2.y
	Next

	Flip
Wend
End




Chroma(Posted 2011) [#6]
Good stuff GW. Now to study and hack some of that into this! I see that above, the Points are kept in a global point list.


How about make each Verlet "object" a separate type with an array for link and points like so:

Type VEntity
   Field points:VPoint[]
   Field links:VLink[]

   Method AddPoint(p:VPoint)
      Self.points :+ [p]
   End Method

   Method AddLink[a:VLink, b:VLink, elasticity#=1.0)
      Self.links :+ [VLink.Create(a, b, elasticity)]
   End Method
End Type


Not sure what a good number for "eslasticity" is...

And looks like I need to put in a "Satisfy Contraints" or "Resolve" method to make the Links try to keep the right length between the two points in the VLink.

Last edited 2011


Chroma(Posted 2011) [#7]
Would you want iterations per frame or iterations per second? I'd probably go with per second. Because I normally detach logic from render anyhow.

I'll probably have a small demo to post here by the weekend. Already been experimenting but the objects are exploding within a few sec...


GW(Posted 2011) [#8]
the Points are kept in a global point list.

Thats more of an artifact of being written in Blitz3d. Arrays are (nearly) always the better solution.

As long as the number of points and links are kept low, appending to the array works ok, but slows down exponentially. Another method would be to add them to a LL then call some kind of finalize() method after creation that converts the list to an array.

Last edited 2011


Wiebo(Posted 2011) [#9]
I would keep the VLinks in a list or array and not the points, because you are satisfying the constraints and not the points. Also, the links contain the points anyway.

Also, I am not sure there is a need for a VEntity. A list (or array) which can be updated should be enough for basic use. In which scenario do you think a VEntitiy would be handy?

I am all for iterations per second.

The idea of changing a LL into an array after setup sounds interesting. I never had any problems with LL being too slow though, I guess it depends on the kind of app you are making.

Question: does the command "Self.points:+[p]" actually change the array size and adds p to the end?? if yes, awesome, I didn't know that.


Chroma(Posted 2011) [#10]
Ah...didn't think of that. Don't need an array/list for points...only on initial creation, then they can be dumped because they're in the Links. Nice.

I just figured have an "Entity" type would be easier for keeping track of and showing/hiding a specific verlet entity when needed.

Also put in AddForce where you just put in the angle from 0 to 360(359.9999) and the amount of force and bam..it pushes all points the same direction. I find this easier than trying to figure out how much x#,y# force to do manually.

does the command "Self.points:+[p]" actually change the array size and adds p to the end??
Aye it sure does.


EDIT: On second thought...not sure about not keeping Points in their own list. What if you need to find a certain Point quickly? It'd be a bit tedious to go through all the links trying to find a specific point because link[0] might not necessarily have point[0] in it.


Here's what it's ballooned to so far. I've added a few extra tidbits here and there.
' MaxVerlet

Type VEntity
	Field points:VPoint[]
	Field links:VLink[]

	Method AddPoint(x#, y#, mass#=1.0)
		Self.points :+ [VPoint.Create(x, y, mass)]
	End Method

	Method AddLink(a:VPoint, b:VPoint, elasticity#)
		Self.links :+ [VLink.Create(a, b, elasticity)]
	End Method

	' Add force at a certain angle with 0 being straight up
	Method AddForce(angle#, force#)
		Local i, fx#, fy#, p:VPoint		
		Local fx# = Sin(angle) * force
		Local fy# = Cos(angle) * force
		For p = EachIn Self.points
			p.fx = fx
			p.fy = fy
		Next		
	End Method

	Method Update()
		' Still researching...
	End Method

	Method Draw()
		Local l:VLink, p:VPoint
		
		' Draw VLinks first
		SetColor 0,0,255
		For l = EachIn Self.links
			DrawLine l.pointA.x, l.pointA.y, l.pointB.x, l.pointB.y
		Next
		
		' Draw VPoints next so they show over VLinks
		SetColor 0,255,0
		For p = EachIn Self.points
			DrawOval p.x-2, p.y-2, 5, 5
		Next

	End Method
End Type



Type VPoint
	Field x#, y#
	Field oldx#, oldy#
	Field fx#, fy#
	Field mass#

	Function Create:VPoint(x#, y#, mass#=1.0)
		Local p:VPoint = New VPoint
		p.SetPos(x y)
		p.mass = mass
		Return p
	End Function

	Method SetPos(x#, y#)
		Self.x = x
		Self.y = y
		Self.oldx = x
		Self.oldy = y
	End Method

	Method Update()
		Local tempx# = Self.x
		Local tempy# = Self.y
		Self.x :+ (Self.x - Self.oldx) '+ Self.fx * dt * dt
		Self.y :+ (Self.y - Self.oldy) '+ Self.fy * dt * dt
		Self.oldx = tempx
		Self.oldy = tempy
	End Method
End Type


' Vlink - attaches two VPoints at a specific elasticity
Type VLink
	Field pointA:VPoint
	Field pointB:VPoint
	Field length#
	field k#

	Function Create:VLink(a:VPoint, b:VPoint, elasticity#)
		Local l:VLink = New VLink
		l.pointA = a
		l.pointB = b
		l.k = elasticity
		Local dx# = a.x - b.x
		Local dy# = a.y - b.y
		l.length = Sqr(dx * dx + dy * dy)
		Return l
	End Function

	Method Update()
		' Still researching
	End method
End Type


Last edited 2011


Chroma(Posted 2011) [#11]
Shapes are outlined with VPoints, and then VLinked together to form the outside shape and then crisscross VLinks are used to give the shape support. I read somewhere about angular constraints so there wouldn't be a need for the inner VLinks if not wanted. This would give a more "rigid" body. Anyone know how to implement an "angular" constraint?

Also...any information anyone already knows about verlet collision would be very helpful.


Chroma(Posted 2011) [#12]
Well I have a demo going. Will post soon.


Qcat(Posted 2011) [#13]
This looks an interesting project I will be relay interested to take a look when this is done.


Chroma(Posted 2011) [#14]
Oh here's what I've been messing about with. Needs tidying up but it's "working". Yeah, needs a lot of tidying up...


Here's the needed Vec2.bmx file:


And here's MaxVerlet:



Qcat(Posted 2011) [#15]
I have been having a mess around with the code I quite impressed with it. Do you have any plans to add collisions?

I have taken a look at the maths involved and it is well above my head!


Chroma(Posted 2011) [#16]
I'm doing a rewrite of the existing code to make it sure clean and organized, then I'm definitely going to tackle collisions. I think it's be smart to start with simple point / line collision or point in box or circle type stuff. I've been looking at so much online source code from various articles...I'm think I'm getting my head wrapped around this pretty decently.

I've seen quite a few posts here on blitzbasic about this too. It seems for almost every post or project of verlet that everything goes really well until it gets to the collision part...then it just dies off. I'm going to avoid that here hopefully. If you want to do some research on verlet collision and post it here it's be MUCH appreciated.


Chroma(Posted 2011) [#17]
Wow, I've been sitting here floundering with getting the angle between two points. It's working on my graph paper and it's working when I input the numbers and have it print it in bmax.

Atan2 (x1-x2,y1-y2)

Pretty simple. Yeah, unless you forget that in screen coordinates that Y is positive going DOWN instead of negative... Funny..


Function AngleBetweenPoints#(x1#, y1#, x2#, y2#)
   Local dx# = x1 - x2
   Local dy# = y1 - y2
   Return Atan2(dx,-dy)   'Yeah, negative dy....
End Function


Gotta love it.


Qcat(Posted 2011) [#18]
The collisions seem to be the stumbling block for a lot of these projects. I started to work on sum-think based on nates tutorial.

http://blitzbasic.com/Community/posts.php?topic=93912

I did find a few snippets info around the web witch might help.

http://www.codeproject.com/KB/GDI-plus/PolygonCollision.aspx?print=true
http://www.gamedev.net/topic/415538-verlet-object-object-collision/
http://codeflow.org/entries/2010/nov/29/verlet-collision-with-impulse-preservation/


Chroma(Posted 2011) [#19]
Looks like the easiest method to start with would be line to line collision. Well we'll see how this goes.


Chroma(Posted 2011) [#20]
This little collision snippet I found in the archives is actually letting me stack blocks. Still needs a lot of work though!


Function PointInShape(p:TPoint, s:TShape)
	Local in = 0
	Local i, x#, y#, x1#, y1#, x2#, y2#

	x = p.curr.x
	y = p.curr.y
	For i = 0 To s.constraint.length-3
		x1 = s.constraint[i].p1.curr.x
		y1 = s.constraint[i].p1.curr.y
		x2 = s.constraint[i].p2.curr.x
		y2 = s.constraint[i].p2.curr.y
		
		If ((((y1 <= y) And (y < y2)) Or ((y2 <= y) And (y < y1))) And (x < (((x2 - x1) * (y - y1)) / (y2 - y1)) + x1))
			in = 1 - in
		EndIf
	Next
	Return in
End Function



Chroma(Posted 2011) [#21]
Looks like it's easier to have the MaxVerletEngine store the verlets shapes as actual shapes instead of just storing a list of all lines.

And it's good to know what constraints are used for the edges and which are used for supports. Makes it WAY easier when doing collision checks.

As in:
Type TShape
   Field point:TPoint[]
   Field edge:TConstraint[]
   Field support:TConstraint[]
End Type


Last edited 2011


Qcat(Posted 2011) [#22]
It sounds like you are making progress! That code snippet sounds like it is doing the job!


Chroma(Posted 2011) [#23]
Slowly but surely. I was quite happy about the stacking. Now I just need to get it to where the lines collide with points. And that should do it for collision. Once that part is working I'll be looking to optimize it so hopefully some of the veterans here will help out if they have time.

And what I mean by edges is that there are constraints that make the edges of the shape and then there's the inside constraints that give it support. Storing the supports along with the edges is counterproductive. Definitely going to split them up.

And with the basic TShape, you can just make a bunch of custom functions to create the shape you need, and even put in a scale value so you can make it any size you want.

Still cleaning up code. My God it's gets messy fast...


Chroma(Posted 2011) [#24]
I'm looking to break the collision checking up just like Blitz3D does.

You specify the collision group for each object (might go to VShape...seems to make more sense in this case) you make. Then do a MaxVE.Collisions(src_group, dest_group) to make them interact.

So something like:
Local ball:VShape = New VShape
ball.SetCollisionGroup(1)

Local box:VShape = New VShape
box.SetCollisionGroup(2)

MaxVE.Collisions(1, 2)


I want to use TLinks for quick access to objects in lists. I think that will be plenty fast.

Last edited 2011


Qcat(Posted 2011) [#25]
The blitz3d method of collision sounds like a good way to set up the collisions. I will be interested to see this working.

Hopeful the method you are using will horsefly be fast!


Chroma(Posted 2011) [#26]
Ever since Wiebo said that I should only track constraints and not points because the points are in the constraints, I've been wracking my brain to try and figure it out.

And here it is! The cool thing is that like Wiebo suggested...it only stores constraints (of two types). The points are updated first as they should be by looping through the edge constraints (which contain all the points in the shape anyhow) and doing a Step 2. Which skips every other constraint and on the last constraint it only updates the first point in the constraint if the shape contains an odd number of constraints. Whew!

Any easy example would be a box and triangle. A box has 4 points. So you would only update the points in constraint 0 and 2..which would effectively be all points in the shape. For a triangle you would have 3 points. Updating constraint 0 would update points 0 and 1. Then skip to constraint 2 with holds points 2 and 0. We've already updated point 0. So the simple "<" check cures that nicely.

The only catch is that you have to build your VShape in a clockwise fashion and make sure you create the constraints in the same manner..point 0 and point 1 counting upwards to the last point created.

Type VShape
	' con_edge contains all VPoints in a VShape and only constraints that form a VShapes edge
	Field con_edge:VConstraint[]

	' con_supp contains only constraints that give inner support (not used for collision)
	Field con_supp:VConstraint[]

	Method AddConstraint(p1:VPoint, p2:VPoint, mode=0, stiffness#=.5)
		Local c:VConstraint = VConstraint.Create(p1, p2, stiffness)
		Select mode
			Case 1
				Self.con_edge :+ [c]
			Default
				Self.con_supp :+ [c]
		End Select
	End Method

	Method Update()
		Local i

		'Refresh all VPoints first (making sure not to refresh a point twice)
		Local conlen = con_edge.length-1
		For i = 0 to conlen Step 2
			constraint[i].p1.Refresh()
			If i < conlen constraint[i].p2.Refresh()
		Next

		' Resolve edge constraints
		For i = 0 to conlen
			constraint[i].Resolve()
		Next

		' Resolve inner support constraints
		For i = 0 to con_supp.length-1
			con_supp[i].Resolve()
		Next
	End Method
End Type


Last edited 2011