Verlet physics tutorial (only the basics)
Blitz3D Forums/Blitz3D Tutorials/Verlet physics tutorial (only the basics)
| ||
| Hi 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 b3d code Type pointmass Field x# Field y# Field ox# Field oy# End Type Function updatepointmasses() For p.pointmass = Each pointmass 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. Graphics 640,480,0,2 SetBuffer BackBuffer() createpointmass(100,100,5,1) While Not KeyDown(1) Cls updatepointmasses() Flip Wend End Type pointmass Field x# Field y# Field ox# Field oy# End Type Function updatepointmasses() Color 255,255,255 For p.pointmass = Each pointmass Oval p\x-4,p\y-4,8,8,0 ;draws an oval at the verlet's position Local dx# = p\x - p\ox ;gets the velocity Local dy# = p\y - p\oy p\ox = p\x ;updates ox and oy p\oy = p\y p\x = p\x + dx# ;adds velocity to get new x and new y 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 Return p.pointmass End Function 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 = Each constraint 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. Graphics 640,480,0,2 SetBuffer BackBuffer() Global CONST_ITERATIONS = 5 p1.pointmass = createpointmass(100,100,1,0) p2.pointmass = createpointmass(100,150,0,0) c.constraint = createconstraint(p1,p2) While Not KeyDown(1) Cls updatepointmasses() updateconstraints() drawconstraints() Flip Wend End Type pointmass Field x# Field y# Field ox# Field oy# End Type Function updatepointmasses() Color 255,255,255 For p.pointmass = Each pointmass Oval p\x-4,p\y-4,8,8,0 ;draws an oval at the verlet's position Local dx# = p\x - p\ox ;gets the velocity Local dy# = p\y - p\oy p\ox = p\x ;updates ox and oy p\oy = p\y p\x = p\x + dx# ;adds velocity to get new x and new y 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 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 = Each constraint 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.pointmass = p1.pointmass c\p2.pointmass = p2.pointmass c\length# = Sqr((c\p1\x-c\p2\x)^2 + (c\p1\y-c\p2\y)^2) Return c.constraint End Function Function drawconstraints() For c.constraint = Each constraint Line c\p1\x,c\p1\y,c\p2\x,c\p2\y Next End Function 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 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: Now the only thing that needs to be done is a collision engine... the hard part (for me at least) |
| ||
| Impressive. Helpful bit of reference material, this. |
| ||
| thanks. glad to know I helped someone understand it more :) |
| ||
| Wow, that's really simple Here's this, I added a triangle and movement controls just to see what would happen. |
| ||
yeah it was even a learning experience for me, the creator because I realized some stuff that I was doing that had slowed my other physics engines down in the past... so here's how to give them a little starting velocity:Global p1.pointmass = createpointmass(400,100,-2,0) Global p2.pointmass = createpointmass(400,150,0,0) Global p3.pointmass = createpointmass(450,150,0,0) Global p4.pointmass = createpointmass(450,100,0,0) Global t1.pointmass = createpointmass(200,50,0,3) Global t2.pointmass = createpointmass(180,100,0,0) Global t3.pointmass = createpointmass(220,100,0,0) put this before the loop the whole concept is really a lot simpler than people make it out to be at first. When I have time, I'll convert this to blitz max for the bmax tutorial section for those that have left blitz 3d |
   |