Energy in a spring mass system
Author
Hi,
I'm putting together my first spring mass system, based on this article:
http://www.gamasutra.com/features/20011005/oliveira_01.htm
My code is basically the same as they have, although I'm only doing calculations in 2D. I'm playing about with a string and a cloth, and they both have problems. The main control over the systems (other than the physical arrangement of the springs) is tweaking values for k and the damping factor. For some values of these, things seem pretty stable, but those aren't the values I want to use ;) I want k to be somewhere around 0.9 (ie not very springy) and the damping to be about 0.9 as well (not too much damping, but enough to let the system settle down relatively easily). The problem is when I do that, as soon as I start trying to drag points around the whole system very quickly spirals into chaos. The accumulated forces grow exponentially and spread throughout the system and the whole thing falls apart.
Unless I've missed something about these systems, the values I want to use should be within reason. Other combinations of values seem to lead to the same problem (although some don't), so I'm guessing the problem is in the modelling of the system somewhere. The individual force calculations seem (to my laymanish eyes) reasonable enough, but when they're accumulated for each point they get scary. I was under the impression that in the less-than-optimal example code that a lot of these forces should be negated - i.e P1 exerts a force of whatever size on P0, at the other end of a spring; but then you check P1, and P0 exerts a similar sized force back, so overall the system stays stable.
I'm kinda working in the dark here though, with not much previous experience of anything other than very basic physics coding. Does anyone have any ideas what might be causing the problem?
Author
Okay, I've had a hunt around for information, and it looks like the problem might well be some shoddy euler integration in the sample code (although let me know if it's something else which could be the problem). Joy of joys, that means I've got to put together something more complex - what would be a recommended method? I've found some stuff on RK4 and Verlet integration, and I gather there might be some others, but having never done this before I don't know which is best to use, how best to use it, and if there are any potential caveats I need to be aware of. Haylp!
indeed its due to the sucky nature of explicit integration.
explicit integration only evaluates forces at the beginning of a timestep, and assumes them to be constant troughout the timestep. but in case of a stiff spring, this isnt anywhere near correct. a stiff spring will oscilate during the timestep, producing a net force wildly different than extrapolating the derivative at the beginning of the timestep, hence a large error appears, which unfortunatly also magnifies itself.
other explicit integration methods might be able to cope better with this problem, but the fundamentals of it remain the same. youre not going to get masive improvements by choosing another explicit scheme.
if you have the time you should look into implicit integration. it can handle arbitrary stiffness without problems. it might seem a little complex, but its not that hard really. its just that there seems to be a complete lack of decent resources explaining it..
explicit integration only evaluates forces at the beginning of a timestep, and assumes them to be constant troughout the timestep. but in case of a stiff spring, this isnt anywhere near correct. a stiff spring will oscilate during the timestep, producing a net force wildly different than extrapolating the derivative at the beginning of the timestep, hence a large error appears, which unfortunatly also magnifies itself.
other explicit integration methods might be able to cope better with this problem, but the fundamentals of it remain the same. youre not going to get masive improvements by choosing another explicit scheme.
if you have the time you should look into implicit integration. it can handle arbitrary stiffness without problems. it might seem a little complex, but its not that hard really. its just that there seems to be a complete lack of decent resources explaining it..
The velocity-less Verlet method is extremely stable. You can implement a stable string/rope/cloth system very quickly and easily. See:
http://www.gamasutra.com/resource_guide/20030121/jacobson_01.shtml
http://www.gamasutra.com/resource_guide/20030121/jacobson_01.shtml
Author
Yup, that Gamasutra article was one of the bits and pieces I found which talked about other methods of integration. Out of interest, how does Verlet cope with things which ARE springy? I know I've been asking about strings and cloths but it would be nice to still keep my options open in terms of using my system for springy things as well. The article says
"The constraints are not guaranteed to be satisfied after one iteration only, but because of the Verlet integration scheme, the system will quickly converge to the correct state over some frames. In fact, using only one iteration and approximating the square root removes the stiffness that appears otherwise when the sticks are perfectly stiff."
Which kinda implies that springiness comes about as a result of deliberately rough calculations as opposed to explicitly having control over k. Is that right?
"The constraints are not guaranteed to be satisfied after one iteration only, but because of the Verlet integration scheme, the system will quickly converge to the correct state over some frames. In fact, using only one iteration and approximating the square root removes the stiffness that appears otherwise when the sticks are perfectly stiff."
Which kinda implies that springiness comes about as a result of deliberately rough calculations as opposed to explicitly having control over k. Is that right?
Yes.
Although you could, I suppose, implement the constraints in a more springy manner, but you might loose the robustness of verlet that way...
Although you could, I suppose, implement the constraints in a more springy manner, but you might loose the robustness of verlet that way...
Quote:Original post by ElectroDruid
"The constraints are not guaranteed to be satisfied after one iteration only, but because of the Verlet integration scheme, the system will quickly converge to the correct state over some frames. In fact, using only one iteration and approximating the square root removes the stiffness that appears otherwise when the sticks are perfectly stiff."
Which kinda implies that springiness comes about as a result of deliberately rough calculations as opposed to explicitly having control over k. Is that right?
Not exactly. When the Verlet constraints are relaxed (e.g. over several frames), the visual result is a soft system, which will appear highly damped and not very spring-like.
You can mix both methods (rigid Verlet and spring), switch methods on the fly, and perhaps even blend (lerp) between methods on the fly. You could model a constraint as a spring for small displacements, then switch (or blend) to the rigid Verlet constraint as the limit is approached (for blending/lerp, you would compute both methods, resulting in new velocities and positions, then lerp between the results. The blend amount, alpha, can be linear or non-linear).
Verlet:
v = xc - xo + a*dt*dt
xo = xc
xc += v
Euler-Verlet:
v += a*dt*dt
x += v
To switch Euler-Verlet(ev_) to Verlet(v_):
v_v = ev_v
v_xc = ev_x
v_xo = v_xc - v_v
To switch Verlet(v_) to Euler-Verlet(ev_):
ev_x = v_xc
ev_v = v_xc - v_xo
or, if you've been computing and storing velocity:
ev_v = v_v
To keep it simple, you can just use Verlet and compute and apply spring forces when the contraint is a spring (apply acceleration but don't directly modify position), otherwise when the constraint is rigid, update the position only (where the new velocity is derived from xc-xo).
This brings to mind a simplification to allow arbitrary blending between rigid and spring constraints:
1. Use the Verlet integration method.
2. During rigid constraint update, if alpha (blend factor) is 0, don't update position. If alpha is .5, update position by 1/2 the amount, if alpha is 1.0, update position the full amount. In general, update position by displacement*alpha.
3. During spring update, set springAlpha = 1 - alpha. Scale the resulting force or acceleration by springAlpha.
This should allow you to blend between rigid and spring constraints by varying alpha. Alpha can be controlled in many ways (a constant factor, linear, non-linear (e.g. curve, smooth function), lookup table, etc.). You could have the constraint behave like a spring then switch (blend) to rigid at the limits (when the spring would have gone unstable).
[Edited by - John Schultz on April 18, 2005 2:08:09 PM]
im not too sure about that whole blending thing. you never know how thats going to mess with stability (unless youre a much better mathematician than me). but its worth a try.
anyway, i got my implicit integrator running just minutes ago, and im never going back. handles all kinds of springness without quircks and its unconditionally stable. all in all its only 200-some lines of code, and in hindsight, it wasnt all that complicated.
unfortunatly i dont have any sources to refer you to. well, i do have them, but they all suck. maybe ill write an article on implict integration that doesnt assume you already know the subject, because they are sorely lacking.
anyway, i got my implicit integrator running just minutes ago, and im never going back. handles all kinds of springness without quircks and its unconditionally stable. all in all its only 200-some lines of code, and in hindsight, it wasnt all that complicated.
unfortunatly i dont have any sources to refer you to. well, i do have them, but they all suck. maybe ill write an article on implict integration that doesnt assume you already know the subject, because they are sorely lacking.
Author
The blending idea sounds interesting. However, right now I'm struggling a bit with just trying to create the solid Verlet cloth thing. My original code was based on this:
http://www.gamasutra.com/features/20011005/oliveira_01.htm
My cloth mesh is generated procedurally, and I've tried two forms of it:
In general, using Euler with the second mesh is too saggy and unstable - there doesn't seem to be enough structure to the cloth. Using Verlet seems to work roughly (but not quite) the same either way.
The code is a bit of a mess right now, but I'll try to post some tidied up snippets to give you an idea of what I've been trying.
When I'm using Euler, AccumulateForces() adds up all the forces for each node, as in the Oliviera article I mentioned about. When I'm using Verlet, AccumulateForces() sets all forces to 0 - the nodes are moved directly in SatisfyConstraints(). By a similar token, SatisfyConstraints() has no effect when I'm using Euler. AnchorCheck() is as it is in Oliviera - my cloth is attached to the top of the screen by the top row of nodes. When I click and drag a node with the mouse, that also becomes an anchor until I release the mouse button.
Spot the work-in-progress code! ;) Obviously in Euler node_curr->force has just been calculated for each node. In Verlet, only the gravity should apply. I've noticed that gravity has far less (if any) effect in Verlet, probably because right after this function I'm setting positions outright, and not using the velocity vector at all (right now I'm storing pos, oldpos and vel for every node, to allow for me switching between different integration methods). I'm not so worried about gravity right now though, it should be fixable.
You can see the commented out Euler integration as well as the Verlet there. There is no damping applied to the Verlet.
There are a few bits of weirdness here I suppose. Partly the instability caused by not normalising the direction vector to move in before scaling it by diff - whether this is an oversight in Jacobson's description of the Verlet stuff or a hack on my part to "fix" some other problem I'm not sure. Also, the directions I'm moving the nodes in is the opposite to in Jacobson's pseudocode, but again, doing things my way gives me something which is a bit like a cloth (ish), and the other way just breaks. The other is NUM_ITERATIONS. I expected to only need to iterate once (as in the code above), but I get a very floppy rubbery sheet. Setting the value up to 10 improves it somewhat, but it's still probably more rubbery than using Euler at low enough values of k to keep things stable. If I set NUM_ITERATIONS much higher than 10, the cloth shakes itself apart.
I'm obviously doing something wrong, but I don't know what. Although the behaviour of the cloth is slightly different using the different techniques, they both suffer the same basic problem that when I drag a point around, the cloth takes a long time to move other nearby points to follow it. So the point of the cloth I've moved stretches a long way, and there's a second or so, during which time the rest of the cloth is largely unaffected before slowly starting to move towards the cursor. Setting dt to 1 (ie just removing the term from the integration) in Euler makes the cloth much more responsive, but also much more unstable, meaning I have to use very low values of k. Doing the same for Verlet makes the cloth hugely soft and unresponsive.
http://www.gamasutra.com/features/20011005/oliveira_01.htm
My cloth mesh is generated procedurally, and I've tried two forms of it:
A--B--- (In which each pair of points is connected twice, e.g there is|\/|\/| a spring between A and B, and another one between B and A)|/\|/\|-------A--B--- (In which each pair of points is connected only once, ie there is| /| /| only one spring connecting A and B)|/ |/ |-------In general, using Euler with the second mesh is too saggy and unstable - there doesn't seem to be enough structure to the cloth. Using Verlet seems to work roughly (but not quite) the same either way.
The code is a bit of a mess right now, but I'll try to post some tidied up snippets to give you an idea of what I've been trying.
// Called once per framevoid ComputeForces (SPRING* spring){ AccumulateForces (spring); ApplyForces (spring); SatisfyConstraints(spring); AnchorCheck (spring);}When I'm using Euler, AccumulateForces() adds up all the forces for each node, as in the Oliviera article I mentioned about. When I'm using Verlet, AccumulateForces() sets all forces to 0 - the nodes are moved directly in SatisfyConstraints(). By a similar token, SatisfyConstraints() has no effect when I'm using Euler. AnchorCheck() is as it is in Oliviera - my cloth is attached to the top of the screen by the top row of nodes. When I click and drag a node with the mouse, that also becomes an anchor until I release the mouse button.
void ApplyForces(SPRING* spring){ //apply resultant forces to increase vel of each spring for (S32 i = 0; i < spring->node_count; ++i) { SPRING_NODE *node_curr = &spring->node; Vector2f acc; acc.x = (node_curr->force.x + gSpringGravity.x); acc.y = (node_curr->force.y + gSpringGravity.y); double dt = gTiming->GetDelta(); // Euler integration // node_curr->vel += acc * dt; // node_curr->pos += node_curr->vel * dt; //energy loss for velocity// node_curr->vel *= spring->energy_loss; // Verlet integration node_curr->pos += node_curr->pos - node_curr->oldpos + (acc * dt * dt); node_curr->oldpos = node_curr->pos; }}Spot the work-in-progress code! ;) Obviously in Euler node_curr->force has just been calculated for each node. In Verlet, only the gravity should apply. I've noticed that gravity has far less (if any) effect in Verlet, probably because right after this function I'm setting positions outright, and not using the velocity vector at all (right now I'm storing pos, oldpos and vel for every node, to allow for me switching between different integration methods). I'm not so worried about gravity right now though, it should be fixable.
You can see the commented out Euler integration as well as the Verlet there. There is no damping applied to the Verlet.
// pseudocoded a bit to make it easier to readvoid SatisfyConstraints (SPRING* spring){ for (/*every node in the mesh*/) { SPRING_NODE *node_curr = &/*current node*/ Vector2f *me = &/*get my position*/ for (/*every neighbour to this node*/) { Vector2f *myneighbour = /*get this neighbour's position*/ int NUM_ITERATIONS = 1; for (int interations = 0; interations < NUM_ITERATIONS; ++interations) { //get the distance Vector2f delta; delta.x = myneighbour->x - me->x; delta.y = myneighbour->y - me->y; F32 distance = delta.Length(); // Normalise the distance vector to just give a direction: // Although Jacobson says use the delta vector itself in // recalculating the position, this shakes my cloth apart in a // matter of seconds Vector2f force_dir = delta; force_dir /= distance; F32 invmass1 = 1.0f; // Masses of 1 at the moment F32 invmass2 = 1.0f; F32 restlength = /*get at rest distance between these nodes*/; F32 diff = (distance - restlength) / ((distance) * (invmass1 + invmass2)); if(diff) { me->x += invmass1 * force_dir.x * diff; me->y += invmass1 * force_dir.y * diff; myneighbour->x -= invmass1 * force_dir.x * diff; myneighbour->y -= invmass1 * force_dir.y * diff; } } } }}There are a few bits of weirdness here I suppose. Partly the instability caused by not normalising the direction vector to move in before scaling it by diff - whether this is an oversight in Jacobson's description of the Verlet stuff or a hack on my part to "fix" some other problem I'm not sure. Also, the directions I'm moving the nodes in is the opposite to in Jacobson's pseudocode, but again, doing things my way gives me something which is a bit like a cloth (ish), and the other way just breaks. The other is NUM_ITERATIONS. I expected to only need to iterate once (as in the code above), but I get a very floppy rubbery sheet. Setting the value up to 10 improves it somewhat, but it's still probably more rubbery than using Euler at low enough values of k to keep things stable. If I set NUM_ITERATIONS much higher than 10, the cloth shakes itself apart.
I'm obviously doing something wrong, but I don't know what. Although the behaviour of the cloth is slightly different using the different techniques, they both suffer the same basic problem that when I drag a point around, the cloth takes a long time to move other nearby points to follow it. So the point of the cloth I've moved stretches a long way, and there's a second or so, during which time the rest of the cloth is largely unaffected before slowly starting to move towards the cursor. Setting dt to 1 (ie just removing the term from the integration) in Euler makes the cloth much more responsive, but also much more unstable, meaning I have to use very low values of k. Doing the same for Verlet makes the cloth hugely soft and unresponsive.
ElectroDruid,
If gravity is not having an effect, something is broken. Make sure your particles work by themselves before working on springs/constraints. With gravity enabled, the particles should fall correctly. After the particles are working, test just one spring/constraint, etc.
Compare the following Verlet constraint update to your implementation:
Eelco,
Implicit integrators are excellent for their stability and generality. How did you deal with energy loss and geometric error (non-symplectic-ness)? Any memory or performance issues for large systems? What method did you use to solve the system of equations?
If gravity is not having an effect, something is broken. Make sure your particles work by themselves before working on springs/constraints. With gravity enabled, the particles should fall correctly. After the particles are working, test just one spring/constraint, etc.
Compare the following Verlet constraint update to your implementation:
void DistanceConstraint::applyConstraints(void) { PhysicalParticle & pa = physicsSystem->physicalParticles[pIndexA]; PhysicalParticle & pb = physicsSystem->physicalParticles[pIndexB]; Vec3 delta = pa.cur - pb.cur;#ifndef USE_NEWTON_RAPHSON_CONSTRAINT_ADJUST flt deltaLength = delta.length(); if (deltaLength > 0.f) { flt diff = (deltaLength - restLength) / (deltaLength * (pa.invMass+pb.invMass)); Vec3 deltaDiff = delta*diff; pa.cur -= pa.invMass*deltaDiff; pb.cur += pb.invMass*deltaDiff; } // if#else if (pa.mass == pb.mass) { delta *= restLength2 / (delta.lengthSqr() + restLength2) - .5f; // Newton-Raphson sqrt approximation with initial guess restLength (1st order Taylor expansion). pa.cur += delta; pb.cur -= delta; } else { flt scaleTotal = 2.f / (pa.invMass+pb.invMass); delta *= restLength2 / (delta.lengthSqr() + restLength2) - .5f; // Newton-Raphson sqrt approximation with initial guess restLength (1st order Taylor expansion). delta *= scaleTotal; pa.cur += delta * pa.invMass; pb.cur -= delta * pb.invMass; } // if#endif} // applyConstraintsEelco,
Implicit integrators are excellent for their stability and generality. How did you deal with energy loss and geometric error (non-symplectic-ness)? Any memory or performance issues for large systems? What method did you use to solve the system of equations?
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
