• Create Account

## Sequential Impulse Solver and Sliding Contacts

Old topic!

Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

20 replies to this topic

### #1RastaRunna  Members

159
Like
0Likes
Like

Posted 21 May 2014 - 01:01 AM

Hi All,

Been implementing a 2d sequential impulse engine (mostly for the academics of it). The main pipeline of my current engine is:

1) Detect Collisions

2) Resolve % of each penetration

3) Repeat 1 & 2 'n' times

3) Resolve velocities

4) Integrate

5) Repeat

This seems to give fairly descent results. However, no matter what I try when objects stack in vertex-edge configurations, or edge-edge in pyramid-like structure, then the angled objects will push others out of the way regardless of how many objects are in the pile being pushed. Here are some videos showing the results:

Initially I thought there was something wrong with my friction implementation but I've tried increasing friction to no avail. Next, for edge-edge cases I only resolve a single contact, so I've experimented with different contact points also to no avail.

I can only think of two things which may be causing this, either:

1) This is un-avoidable for a sequential impulse solver -- which doesn't feel entirely correct. sliding from sequential solvers should be gradual drifts; not this dramatic but perhaps I am wrong.

2) Steps 1 & 2 in my loop place a large emphasis on resolving penetrations (via position corrections) rather than giving more weight to #3. As I've been reading around online, it seems many don't even do #2 but rather only resolve the velocity by applying the appropriate impulse plus some small portion of the penetration (baumgaurte stabilization?).

As far as my penetration resolution implementation: it resolves the penetrating objects in the direction of the contact normal according to the objects' mass. Each pair is resolved  in the order of greatest separating velocity.

Simple answer seems to be put everything to sleep.... though I feel like that's just masking a stability problem.

I'd appreciate any feedback on the ordering of my main loop, and any suggestions on what may be causing my stacks to push apart / how I might be able to mitigate that? Would a TOI solver fix this, or is simultaneous resolution the only answer?  Anyways, I feel like I have lots of ideas, but no clear direction; so some guidance would be very helpful.

Thanks

Edited by RastaRunna, 21 May 2014 - 01:14 AM.

### #2Dirk Gregorius  Members

2438
Like
0Likes
Like

Posted 21 May 2014 - 11:21 AM

Are you warmstarting the friction? You can get pretty bad friction otherwise.

### #3RastaRunna  Members

159
Like
0Likes
Like

Posted 21 May 2014 - 03:20 PM

I'm not. Any recommended ways for doing this? Or would simply feeding in the previous frame's friction impulse work?

Here's my current friction implementation:

Geometry::Vector2D tangent;
tangent.x = relVel->x - normal->x  * relativeVelOnNormal;
tangent.y = relVel->y - normal->y  * relativeVelOnNormal;
Geometry::Normalize(tangent, tangent);
//tangent.x = -normal->y;
//tangent.y =  normal->x;

flttype planarVel = (relVel->x * tangent.x + relVel->y * tangent.y);
flttype changeInPlanarVel = planarVel;
flttype relAPT = (posRelA->x * tangent.y - posRelA->y * tangent.x);
flttype relBPT = (posRelB->x * tangent.y - posRelB->y * tangent.x);

flttype relACrossTangent = entityAInertiaTensorInv * relAPT;
flttype relBCrossTangent = entityBInertiaTensorInv * relBPT;

flttype angularMassAT = relACrossTangent * relAPT;
flttype angularMassBT = relBCrossTangent * relBPT;
flttype frictionImpulse = changeInPlanarVel / (curContact->totalLinearMass + angularMassAT + angularMassBT);

// clip friction mag to min (impulseMag * frictionCoefficient, frictionMag)
flttype frictionMag;
if (frictionImpulse < curContact->friction * impulseMag)
{
// static friction
frictionMag = -frictionImpulse;
}
else
{
// dynamic friction
frictionMag = -curContact->friction * frictionImpulse;
}

/////////////////////////////
//
//   Update velocities
//
//

// linear vel
entityA->velocity.x -= (impulseX + frictionMag * tangent.x) * entityAInvMass;
entityA->velocity.y -= (impulseY + frictionMag * tangent.y) * entityAInvMass;

entityB->velocity.x += (impulseX + frictionMag * tangent.x) * entityBInvMass;
entityB->velocity.y += (impulseY + frictionMag * tangent.y) * entityBInvMass;

// angular vel
entityA->angularVelocity -= (relACrossNormal * impulseMag + relACrossTangent * frictionMag);
entityB->angularVelocity += (relBCrossNormal * impulseMag + relBCrossTangent * frictionMag);


Edited by RastaRunna, 22 May 2014 - 12:37 AM.

### #4Dirk Gregorius  Members

2438
Like
4Likes
Like

Posted 21 May 2014 - 05:19 PM

First I would make sure that your tangent doesn't go to zero when the relative velocity vanishes. What is the impulseMag? Is this the incremental or the accumulated impulse? You need to clamp against the accumulated impulse.

I would clamp the friction impulse like this (assuming impulseMag > 0) :

frictionImpulse = clamp( frictionImpulse, -curContact->friction * impulseMag, curContact->friction * impulseMag );

### #5RastaRunna  Members

159
Like
0Likes
Like

Posted 22 May 2014 - 12:26 AM

Clamping the friction as you mentioned seems to have helped a bit, though at first there is definitely this odd jarring to right of the pyramid before it falls. Once it falls though, the objects do settle down pretty quickly without pushing apart as much as before, though that could just be due to the configuration of how the blocks lie once they fall:

impulseMag is supposed to be the total impulse required to cause the appropriate change in velocity. I assume it is the accumulated value because I solve for this impulse in one go rather than in several small iterations (per time step). Or does accumulated mean something else? Here is how I compute impulseMag:

flttype relativeVelOnNormal = (relVel->x * normal->x + relVel->y * normal->y);

if (relativeVelOnNormal > 0.0)
continue;

// get applied velocity from previous step along normal (i.e. gravity)
flttype velFromAcc = (entityB->appliedVel.x - entityA->appliedVel.x) * normal->x +
(entityB->appliedVel.y - entityA->appliedVel.y) * normal->y;

// determine desired change in velocity (w/out applied vel from prev step)
flttype changeInVel = -velFromAcc + (-1.0 - curContact->coefficientOfRestitution) * (relativeVelOnNormal - velFromAcc);

// determine impulse required to cause desired change in vel
flttype relAPN = (posRelA->x * normal->y - posRelA->y * normal->x);
flttype relBPN = (posRelB->x * normal->y - posRelB->y * normal->x);

flttype relACrossNormal = entityAInertiaTensorInv * relAPN;
flttype relBCrossNormal = entityBInertiaTensorInv * relBPN;
flttype angularMassA = relACrossNormal * relAPN;
flttype angularMassB = relBCrossNormal * relBPN;

flttype totalLinearMass = entityAInvMass + entityBInvMass;
flttype impulseMag = changeInVel / (totalLinearMass + angularMassA + angularMassB);

// keep impulse positive
impulseMag = impulseMag < 0.0 ? 0.0 : impulseMag;


Edited by RastaRunna, 22 May 2014 - 12:30 AM.

### #6Dirk Gregorius  Members

2438
Like
0Likes
Like

Posted 22 May 2014 - 09:50 AM

1) Detect Collisions

2) Resolve % of each penetration

3) Repeat 1 & 2 'n' times

You don't re-detect collisions. You iterate only on the current contact points. I recommend looking at Box2D Lite to see how this works. Erin also published a bunch of papers that should help you.

Edited by Dirk Gregorius, 22 May 2014 - 03:19 PM.

### #7RastaRunna  Members

159
Like
0Likes
Like

Posted 22 May 2014 - 12:00 PM

Isn't it possible though that by resolving one contact pair, new contacts or collisions are created?

### #8Dirk Gregorius  Members

2438
Like
0Likes
Like

Posted 22 May 2014 - 03:19 PM

No, you only change velocities. New contacts would be only created if you also update the position/orientation after each impulse. That's not how it works

### #9RastaRunna  Members

159
Like
0Likes
Like

Posted 22 May 2014 - 05:53 PM

Ah, okay that makes sense. I was thinking in this manner: fix overlap (via position corrections) then resolve velocity (via vel corrections or impulse). But from what you're saying it sounds like I should not fix overlap directly via position corrections, but only through impulses. Is that correct?

### #10Dirk Gregorius  Members

2438
Like
3Likes
Like

Posted 23 May 2014 - 09:41 AM

Yes! You first solve the velocities. The projection steo give you the 'best' velocities to advance to the positions. Then you correct positions. I would not re-run the collision detection though, but keep contact points local and measure the penetration along the normal.

Again, I recommend having look at Box2D and Box2D Lite.

### #11RastaRunna  Members

159
Like
0Likes
Like

Posted 23 May 2014 - 10:16 AM

Okay, I'll give that a try and post my results. Thanks for your help, I do appreciate it. Sorry if I keep beating a dead horse -- Im trying to really understand what's going on more than just copy code. Any intuition on why this gives better results (besides being correct)?

### #12RastaRunna  Members

159
Like
0Likes
Like

Posted 24 May 2014 - 12:46 AM

So came across slide #54+ of Erin's "Understanding Constraints" presentation titled "Why not work with positions directly?" which kind-of answers my question (though it seems to be addressing "position-only" solvers rather than position and vel solvers):

I think it's sounding like both solvers could be accurate, but overkill since a velocity solver should do just fine in most cases according to the presentation.

### #13Dirk Gregorius  Members

2438
Like
4Likes
Like

Posted 24 May 2014 - 10:50 AM

There are two things with position solvers.

1) Friction and motors are usually modeled as velocity constraints. It is tricky to do this on the position level

2) If you solve on the position level you need to solve a non-linear system. Solving on the velocity level is a linear problem. Also constraint equations and Jacobians are functions of the positions and orientations. So you need to re-evaluate them while you iterate over them. This is *pretty* costly. On the velocity level they remain constant and can be pre-computed. For particle systems it is a great choice though. Check all the Position-Based-Dynamics (PBD) work by Mueller and the original presentation by Jacobsen using Verlet integration.

If you really want to understand what is going on, here are a couple of links:

http://www8.cs.umu.se/kurser/5DV058/VT09/lectures/old/Lecture11_Claude2006.pdf

http://www8.cs.umu.se/kurser/5DV058/VT09/lectures/old/Lecture12a.pdf

http://www8.cs.umu.se/kurser/5DV058/VT09/lectures/old/Lecture12b.pdf

http://www8.cs.umu.se/kurser/5DV058/VT09/lectures/spooknotes.pdf

Everything from Erin Catto here:

These are good as well:

http://www.digitalrune.com/portals/0/documents/a%20unified%20framework%20for%20rigid%20body%20dynamics.pdf

Finally you can read also everything from David Baraff and Brian Mirtich. Just google for those names

If you want to write your own physics engine I recommend starting from Box2D Lite and then port it slowly over to 3D. Once you have the basic understanding you can improve broadphase, collision detection, contact persistence, etc. based on you personal preferences. You can use the real Box2D as guidance then

Edit: Here are links to the Mueller and Jacobsen papers

http://www.matthiasmueller.info/

For a formal background look into technical mechanics books or online lectures.

Happy studying!

### #14RastaRunna  Members

159
Like
0Likes
Like

Posted 24 May 2014 - 01:03 PM

Wow, this is great info, and will definitely keep me be busy for a while . Thank you for taking the time to list all of these.

I've had an interest in 2D physics simulations for the past 12 years; started in high school building really small orbit and spring simulations. About 5 years ago I decided I want to build an engine so started collecting / reading books in my spare time (so very slowly). My first go at it I built a basic engine w/spatial grids for broadphase, and SAT for narrow but then got stuck on resolution so essentially copied the parts of Box2D I didn't understand. Feeling unsatisfied that I didn't really understand the resolution math and wanting to do something in C++ I started over back in November and built an engine that does GJK / MPR for narrow (saving broadphase for after I get resolution working).

But all of this time I've always had this notion of having to solve penetration through position corrections. I think my previous mindset came from "Game Physics Engine Development" by Ian Millington. For instance Pg. 122 he says "We could combine interpenetration resolution with collision resolution. A better solution, in practice, is to separate the two into distinct phases. First we resolve collisions... Second we resolve the interpenetrations."  And I think I just always assumed that "interpenetration resolution" meant fixing pos... Which doesn't make any sense looking back because the book only ever talks about impulses.... Ahhh!! All this time....... (bangs head on desk).

Anyways, an epiphany for a slow guy. Thanks again for your help.

### #15RastaRunna  Members

159
Like
0Likes
Like

Posted 30 May 2014 - 01:29 AM

Okay, just to update with some results...

So I went back and revamped my resolution system; gutted out the penetration resolution via direct position corrections and switched to impulse only adjustments. This along with warmstarting my impulses (both normal and tangential) made a rather dramatic difference to both the stability and performance of my system.

Here are some videos demonstrating the improved results (sorry for the poor compression quality... but I think you'll get the point. Also I made the blocks smaller so I could fit more on screen):

Here's one with 50 rows (slow w/out any broadphase implementation, but stable!!):

http://youtu.be/zNgasUe-IKg

Another with 30 rows and projectile box w/higher mass:

http://youtu.be/hR2_JyLQtMY

Also, I went back to the book I had been following for so long and realized that it indeed does focus on doing what I was doing wrong: direct position corrections. So I'm not as slow / crazy as I was beginning to think; just misled.   Direct position corrections is so much slower due to the constant need to re-run the detection system (as Dirk already mentioned).

I was actually very surprised how easy it was to apply penetration corrections via impulses. Simply convert the contact's overlap by the change in time to obtain a velocity and add it to your desired change in velocity when computing your resolution impulse. You'll want to only take a fraction of it otherwise you end up adding energy to the system resulting in jitters. Also you'll want to clamp overlap to be positive (which I'm not showing here):

contact.penetrationAsVel = amntOverlapResolve * contact.overlap * dtInv;
impulseMag = (changeInVel + contact.penetrationAsVel) * contact.effectiveNormalMassInv;


Thanks Dirk for the help.

### #16Dirk Gregorius  Members

2438
Like
0Likes
Like

Posted 31 May 2014 - 10:40 AM

Nice work! I am glad you figured it out!

### #17h4tt3n  Members

1917
Like
0Likes
Like

Posted 04 June 2014 - 12:38 PM

Your videos look very impressive! Nice calm stacks of boxes. Would you perhaps list a few of the ressources ( websites, books, articles ) you found the most helpful in getting these results? I'm building a physics engine in c++ and am at the moment finishing the heavily optimized spring part and Keplerian orbit part, and my next big task will be implementing collision detection and resolution. Cheers Mike.

### #18RastaRunna  Members

159
Like
2Likes
Like

Posted 11 June 2014 - 10:01 PM

Thanks!

For collision detection my engine implements a mix between GJK and the Xenocollide algorithm with Minkowski Portal Refinement (MPR) to get contact points. You can find the algorithm described in detail in Game Programming Gems 7 section 2.5. A simple 2D example is given on the authors website. Other special-case detection I follow Real-Time Collision Detection by Ericson modifying to my specific needs.

For resolution I used a mixture of the following resources to understand the equations of motion and implement a sequential impulse solver:

• Game Physics Engine Development by Ian Millington -- while a good resource for general physics engine development, I do not recommend following this approach if you want a stable system for stacking more than 3 or 4 boxes high. The latest edition has an addendum for 2D physics if that is your interest.

• Game Physics by David H. Eberly -- This book is very thorough and gives lots of details on the math for resolution. At an implementation level this book tends to be more focused on block solvers and simultaneous resolution.

Finally, I followed Box2D lite for ordering the integration / detection / resolution pipeline as well as for how to resolve via impulses only (with no direct position corrections), and how to best do warm-starting.

All this said, the two features that added the most to the stability were:

1) Warm starting on normal and tangential impulses

2) Impulses for penetration correction

3) Adjusting the coefficient of restitution (bouncy objects don't lend well to stacking). Box2D lite effectively sets this to 0 along the normal if you look in the code.

Hope this helps!

Edited by RastaRunna, 11 June 2014 - 10:07 PM.

### #19Finalspace  Members

831
Like
0Likes
Like

Posted 14 June 2014 - 02:07 AM

Haha i am in the exactly same position - dreaming since years about making a physics engine. Had no good results so far - therefore had a pause for a while, but got back in the last few years and learned all the math stuff for integration / collisions / response etc. - therefore i understand everything a lot better. Now i am in progress of making the "sequential impulse solver" thingy and had trouble to get it stable - but now after i integrated the position correction as a velocity impulse - it got stable yay. But it seems that my restitution stuff is now off - it creates more energy over time. (A bouncing ball with 1.0 restitution jumps even higher over time).

One thing which bugs is me is the "appliedVelocity" in your example. From the name i thought its the applied acceleration multiplied by delta time for current/last? frame maybe?

Also it would be nice, if you could explain what you do exactly for modifying bounciness - is it something like this in the contact initialization process?:

contact.restitution = Math.max(contact.bodyA.restitution, contact.bodyB.restitution); // Between 0 and 1
math.vec2Sub(relVelocity, contact.bodyB.velocity, contact.bodyA.velocity);
var velAlongNormal = math.vec2Dot(contact.normal, relVelocity);
if (velAlongNormal > this.velocityThreshold) { // threshold are just one pixel
contact.restitution = 0;
}


Thanks in regards,

Final

Edited by Finalspace, 14 June 2014 - 02:14 AM.

C#, Java, JavaScript are nice languages - Very Easy to code with - BUT it shakes your brain when you swap back to C or any other low level language...

### #20RastaRunna  Members

159
Like
0Likes
Like

Posted 15 June 2014 - 11:47 PM

The example above that uses appliedVelocity was part of my less-stable version which I am not using anymore. The idea comes from "Game Physics Engine Development" by Ian Millington (section 16.1.2 of 2nd ed) and is intended to remove the previous frame's applied velocity (primarily acceleration due to gravity) in order to reduce jitteriness. With position correction integrated, my engine was stable enough that the term wasn't needed anymore, so I've since removed it.

If you still want to use it, my understanding is that you determine it's magnitude along the collision normal (appliedVel dot normal), and subtract it from the amount of desired change in velocity for that frame in order to prevent over-compensation.

Regarding restitution, for my needs, I really don't need much to be "bouncy", so I'm setting the change in vel to be only enough to remove the relative velocity w.r.t the normal (i.e. changeInVel = -magnitudeOfRelativeVelOnNormal). But anyways, the added energy may be coming because of the integration of the position correction. Remember position correction only needs to be applied for a single frame; and so when we apply it as an impulse, it will continue to affect the bodies in question over several frames. To account for this, you should experiment with only applying a fraction of this corrective impulse. This combined with warm starting should be enough. If you dig into the full Box2D (rather than Box2D lite) I believe it handles coefficients of restitution more dynamically.

Something else to consider. If you are caching your contacts across multiple frames or have implemented any type of relaxation to solve your velocity constraints, you'll want to make sure that you update the  "velAlongNormal" each time you apply an impulse to either body otherwise you may accumulate error which *could* be causing any extra bounciness.

Old topic!

Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.