I have an impulse-based position constraint for contacts along the lines of this paper by Bender & Schmitt.

I'm reaching a case where the loop doesn't converge for the tangent (after being stable for about 20 seconds or 700 frames). The normal reaches target after about 4 iterations, but the tangent will only get so close before getting pushed further away. Continuing to run the loop, the normal will begin to oscillate over a negative and positive range, centered on zero, and the tangent will oscillate in some positive range but never reach the target.

In the paper, the tangent unit vector `t` is calculated by dropping the distance vector `d` between the contact points at the end of the timestep, onto the plane of the contact normal unit vector `n` at the end of the timestep: `t = normalize (d - n*dot(d,n))`. I noticed in Catto that he uses 2 tangent unit vectors, though that paper he is calculating a velocity constraint, not a position constraint. But with the tangent calculated above, the 2nd tangent given by `cross (n,t)` is going to have a distance of 0.0 so no impulse will be applied in that direction.

Possibilities:

1) calculate tangents based on relative velocities instead of distances

2) switch to a velocity constraints

3) add more points to the contact and solve pairs iteratively (currently it's only one pair of contact points)

I'd say do option 1. The tangent vector is for applying contact friction, so surely it has to be calculated from the relative velocities of the contact points? I had a read through the paper and it seems to me that the tangent vector is described in the middle of page 4, as this.

That's quite a nice, readable paper :-)

Is this 2D or 3D?

By the papers I've read so far I've never seen one satisfying velocity dependent friction constraints in a position-based dynamics (if you guys have any I'll be happy to read).

Anyhow, if you switch to a velocity constraint based approach in 3D then you also want to make sure you still have a tangent basis (or bases if you're using more than 3 linearized directions) even when if the tangential speed vanishes. One way to accomplish this is as follows.

tangents[0] = vRel - vRelN * normal; //vRel = relative velocity at the contact point, vRelN = Dot(vRel, normal)
float tangentLenSq = LenSq(tangents[0]);
if (tangentLenSq > epsilon * epsilon) {
    float invLen = 1.0f / Sqrt(tangentLenSq);
    tangents[0] *= invLen;
    tangents[1] = Cross(normal, tangents[0]);
}
else {
    ComputeBasis(normal, tangents[0], tangents[1]);
}

Of course if the normal doesn't change very much then you will be fine just with a basis. The rest is the same as in Catto (e.g. solving it as a motor joint).

For games I recommend 2) and 3) and 1), respectively.

I'd say do option 1. The tangent vector is for applying contact friction, so surely it has to be calculated from the relative velocities of the contact points? I had a read through the paper and it seems to me that the tangent vector is described in the middle of page 4, as this.

That's quite a nice, readable paper :-)

That is the tangent used for collision response, the tangent for contacts is described on page 6.

So it turned out I was looking at the Catto "Iterative Dynamics" paper, which is a bit more advanced than this scenario, and actually what I had been using before was Catto's GDC2006 powerpoint, which I had previously compared to Bender and Schmitt and found to be essentially the same, although Catto goes into picking points on boxes and some other tricks such as a 'bias' velocity.

Then I would say it's time to add more contact points, but I have general meshes, not only boxes. Does anyone have any pointers for finding contact points on a general surface? I'm guessing the brute force approach is just searching along the surface normal and comparing dot products.

Forget about bias velocities and the Baumgarte scheme, since they artificially add more energy to the system.

I think Erin wanted to emphasize the performance of these algorithms in the old slides, which is quite cheap.

If you really want an stable simulation, then use the Full-Non Linear Gauss Seidel, implemented in Box2D.

Contact creation is hard to get it right. You need a combination of frame coherence, stable manifolds, numerical robustness, and so...

At the very end, you will need SAT (for penetrating shapes) and GJK (for overlapping shapes, most of them are symmetric) in your physics engine in order for you to create correct manifolds. So, the adventure can become quite large, if you choose to support various shapes.

For creating contact manifolds between convex shapes, look at Dirk's Gregorius presentations on the Downloads section at box2d.org. You'll find there what you need.

Thanks for the pointers, I was going to modify my post to mention that I do use GJK already for distance and penetration calculation, but it does things pretty minimally so it usually returns a point/face as nearest features in a resting scenario.

I assume you're using EPA or another algorithm for the penetration calculation. I made no good experience with this in the past. Specially for quadratic shapes and degenerate cases. Here, we recommend you to use SAT instead. But keep your GJK. You will definately need it in your engine later.

We get this type of question a lot. For removing duplicate replies, I recommend you to read this thread

http://www.gamedev.net/topic/667499-3d-sat-problem/

Oh, I forgot mentioning that I have an actual contact creation implementation using the SAT if you need some reference code:

https://github.com/irlanrobson/Bounce/blob/master/Src/Collision/b3SAT.cpp

This code actually is not fully optimized, and I haven't release the new version of this engine with new features, but it is fundamentally what should be done to get stable manifolds!