More verlet stuff
Blitz3D Forums/Blitz3D Programming/More verlet stuff
| ||
As I mentioned in my previous topic, I was going to mess around with learning verlets. I think I've got it, take a look at this little program I've just written:Graphics 1280, 768, 32, 2 SetBuffer BackBuffer() SeedRnd MilliSecs() Const TIMESTEP# = .995 Type dat Field x#, y# Field h#, v# Field oldx#, oldy# Field fric# Field mass# End Type Global p.dat Global temp.dat Type node_dat Field x1#, y1#, x2#, y2# End Type Global node.node_dat For n=0 To 19 p.dat = New dat p\x = Rand(GraphicsWidth()) p\y = Rand(GraphicsHeight()) p\oldx = p\x p\oldy = p\y p\h = Rnd(-1, 1) p\v = Rnd(-1, 1) p\fric = Rnd(.9, 1) p\mass = Rnd(1, 30) Next ; An extra particle for creating a central 'sun', etc. p.dat = New dat p\x = GraphicsWidth() / 2 p\y = GraphicsHeight() / 2 p\oldx = p\x p\oldy = p\y p\h = 0 p\v = 0 p\fric = 0 p\mass = 100 Repeat Cls For p.dat = Each dat p\v = 0 p\h = 0 ; Calculate pull for each particle interaction For temp.dat = Each dat xdiff# = temp\x - p\x ydiff# = temp\y - p\y dist# = Sqr(xdiff * xdiff + ydiff * ydiff) If Abs(dist#) < 100 Then node.node_dat = New node_dat node\x1 = p\x: node\y1 = p\y node\x2 = temp\x: node\y2 = temp\y End If If Abs(dist#) > 0 Then p\h = p\h + (xdiff / dist) / (dist * p\mass / 10) p\v = p\v + (ydiff / dist) / (dist * p\mass / 10) ;p\h = p\h + (xdiff / dist) / p\mass ;p\v = p\v + (ydiff / dist) / p\mass End If Next ; Verlet temp_x# = p\x temp_y# = p\y p\x = p\x + p\x * p\fric - p\oldx * p\fric + p\h * TIMESTEP p\y = p\y + p\y * p\fric - p\oldy * p\fric + p\v * TIMESTEP If p\x < 0 Then p\x = 0 If p\x > GraphicsWidth() Then p\x = GraphicsWidth() If p\y < 0 Then p\y = 0 If p\y > GraphicsHeight() Then p\y = GraphicsHeight() p\oldx = temp_x# p\oldy = temp_y# Next ; Draw node lines Color 32, 32, 64 For node.node_dat = Each node_dat Line node\x1, node\y1, node\x2, node\y2 Delete node Next ; Draw particles For p.dat = Each dat Color 128, 128, 255 radius# = p\mass * .5 Oval p\x - radius, p\y - radius, p\mass, p\mass Next Flip Until KeyHit(1) End It simply simulates gravity around a point and an attempt to put some node stuff in there as well. Have a play around with a values, it's kinda fun. :) EDIT: Note to self: Use more comment lines. ;) |
| ||
Thanks. this is very interesting stuff and something I'll look into when I finish my current project. One thing that's confusing me after a little look is adding a counter for nodes, I'm making about 4 times as many nodes as lines visible. Your not making checks between objects and self though, or duplicate A<>B and B<>A, so I guess that's where the extra lines are coming from, and this also means node resolving is happening twice per particle, right? As an extra note, you don't need Abs(dist#) following a square root(dist#) as a square root is always positive (unless you're dealing with imaginary numbers). |
| ||
You're correct. I'm not doing any checking for duplicate nodes at the moment, the nodes thing was added in afterwards as just an experiment. I'm wondering how I could do it actually. I'm thinking perhaps the node type could store the handles of both the particles, and when creating a new node it checks to see if a node connecting those two particles already exists. If it works it might be an interesting way of making AI follow paths or something similar. The square-root thingy, that's just a bad habit. Ignore it. :D |
| ||
Yup, work by the look of things - take a look:Graphics 1280, 768, 32, 2 SetBuffer BackBuffer() SeedRnd MilliSecs() Const TIMESTEP# = .995 Type dat Field x#, y# Field h#, v# Field oldx#, oldy# Field fric# Field mass# Field rad# End Type Global p.dat Global temp.dat Type node_dat Field x1#, y1#, x2#, y2# Field link[2] End Type Global node.node_dat For n=0 To 29 p.dat = New dat p\x = Rand(GraphicsWidth()) p\y = Rand(GraphicsHeight()) p\oldx = p\x p\oldy = p\y p\h = Rnd(-1, 1) p\v = Rnd(-1, 1) p\fric = 1;Rnd(.9, 1) p\mass = Rnd(5, 80) p\rad = p\mass * .1 If p\rad < 1 Then p\rad = 1 Next Repeat Cls For p.dat = Each dat p\v = 0 p\h = 0 ; Calculate pull for each particle interaction For temp.dat = Each dat xdiff# = temp\x - p\x ydiff# = temp\y - p\y dist# = Sqr(xdiff * xdiff + ydiff * ydiff) If dist# < 100 And p <> temp Then node_duplicated = False For node.node_dat = Each node_dat If node\link[0] = Handle(p.dat) And node\link[1] = Handle(temp.dat) Then node_duplicated = True If node\link[0] = Handle(temp.dat) And node\link[1] = Handle(p.dat) Then node_duplicated = True Next If node_duplicated = False Then node.node_dat = New node_dat node\x1 = p\x: node\y1 = p\y node\x2 = temp\x: node\y2 = temp\y node\link[0] = Handle(p.dat) node\link[1] = Handle(temp.dat) End If End If If dist# > 0 Then ;p\h = p\h + (xdiff / dist) / (dist * p\mass / 10) ;p\v = p\v + (ydiff / dist) / (dist * p\mass / 10) p\h = p\h + (xdiff / dist) / p\mass p\v = p\v + (ydiff / dist) / p\mass End If Next ; Verlet temp_x# = p\x temp_y# = p\y p\x = p\x + p\x * p\fric - p\oldx * p\fric + p\h * TIMESTEP * TIMESTEP p\y = p\y + p\y * p\fric - p\oldy * p\fric + p\v * TIMESTEP * TIMESTEP If p\x < 0 Then p\x = 0 If p\x > GraphicsWidth() Then p\x = GraphicsWidth() If p\y < 0 Then p\y = 0 If p\y > GraphicsHeight() Then p\y = GraphicsHeight() p\oldx = temp_x# p\oldy = temp_y# Next ; Draw node lines Color 32, 32, 64 nodes = 0 For node.node_dat = Each node_dat nodes = nodes + 1 Line node\x1, node\y1, node\x2, node\y2 Delete node Next ; Draw particles For p.dat = Each dat Color 128, 128, 255 radius# = p\rad * .5 Oval p\x - radius, p\y - radius, p\rad, p\rad Next Text 10, 10, "Nodes: " + nodes Flip Until KeyHit(1) End I found that easier than what I thought it would be! |
| ||
I might try and simulate accretion next. That should be pretty easy to do I think, when two particles are at a certain distance and a slow enough velocity so that they cannot escape each others pull they could be considered as 'one' particle. |
| ||
Here's a challenge that I haven't solved with a look at this verlet shennanigans - have a 'rope' stretched across the screen. It should sags a bit due to gravity. And then jiggle one end so that the jiggles move through the length of the rope. By editing out the x component to get a 2D verlet engine in only the y direction, the transmission of the impulse on one end along the length of the rope dies out quick, and I get creepage as the rope drifts, and generally it doesn't work. I also can't see where to specify a range constraint so a point can't move further than n away from it's neighbours. So all in all I'm not sure I get how this works! ;) eg. It seems increasing the value of mass decreases the 'weight' of a particle. |
| ||
Well I've just solved my challenge myself, by returning to the original example posted in the rope thread. Obviously it's doing something (ConstrainMass ?) that your above example isn't, so when I changed your above code I couldn't get the same ropey behaviour. I'm not sure what you're doing above is 'true' verlet-ness because you haven't got a per-particle constraint thing. Aren't you only summing and averaging per particle? I really ought to look harder at the papers, but it'd be nice if they gave more code examples rather than maths equations! |
| ||
Hmm... I've used the formulas given in the Gamasutra article. I even did a little test in Excel to see how the 'Verlet' system I programmed above works different to the simpler Euler system I used in the previous thread. This is how it works as far as I can gather; Using verlets the velocity variable isn't cumulative. It's just a sum of all the forces acting on the particle during that cycle. After that it resets the velocity to 0. Say in the first cycle you set the Y velocity to 0.1; the particle will move 0.1 down the screen. In the next cycle, the velocity is set to 0, but the particle should still move 0.1 units down the screen. This is because it's using the previous coordinates to calculate the next position: the two positions should be exactly 0.1 units apart. The rope example I posted in the other thread I don't think is 'correct' - it works to a degree, but I don't think it's the typical way of doing it. Tomorrow I'm going to have a go at doing the rope program using the 'verlet' system I used above, I have an idea of how I might do it. |
| ||
I also can't see where to specify a range constraint so a point can't move further than n away from it's neighbours. Ah - I did that using unit vectors. I don't know if you know how you use them, but I'll go over it anyway. A unit vector is a vector with a length of exactly 1. Take a vector v(4, 5) for example. That's 4 units up, and 5 units across. Using pythags theorem, you can work out the length to be: vector length = sqr(4 * 4 + 5 * 5) ...Which equals approx. 6.40 units. We want that to be a length of 1, but still point in the same direction. To work out the X & Y coords of the unit vector you use the following: X = X / dist Y = Y / dist This gives you the unit vector. To set the length of the vector, you simple multiply the coords by whatever amount you want it to be. E.g., if your vector needs to be 20 units long, then it's simply: X = (X / dist) * 20 Y = (Y / dist) * 20 When you're calculating each location for the rope joints, you use the previous coordinates to work out where the next joint should be. You make the next joint 'n' units away from the previous point in whatever direction it's currently pointing. I hope that makes some kind of sense! I pretty much taught myself all that over a weekend, so I'm still learning. ;) |
| ||
In the next cycle, the velocity is set to 0, but the particle should still move 0.1 units down the screen. This is because it's using the previous coordinates to calculate the next position: the two positions should be exactly 0.1 units apart. I think you mean acceleration, rather than velocity. It's set to 0 after each cycle because acceleration generally isn't accumulative (for example, gravity, which accelerates objects downward, is constant). Both Verlet and Euler are "integrators" (I think that's the word used...). Euler is an approximate method, while Verlet uses methods derived from the simple calculus problem of integrating the force acting on an object: velocity = (integral of) gravity (dt) velocity = acceleration * t + initial_velocity position = (integral of) velocity (dt) position = (acceleration/2) * t * t + initial_velocity * t + initial_position (Where t = TIMESTEP) It really starts to get fun when you connect particles together with rigid constraints, which form rigid bodys (like described in the Gamasutra article). |