Keep points in local space and transform them to world. It's not really slow in practice. Obviously this is slower than baumgarte, but the solver should already be very fast, and this memory should already be in the cache. Plus, only a couple iterations of NGS are necessary for pretty good results.

The method I'm using is to reconstruct it from the contact points, based off of what kind of shapes I'm dealing with and what kind of feature they came from. If there was a plane hitting a point, I'd store the plane in local space of one shape, and the point in the other. Then the plane can be rotated into world space, and the point transformed.

Obviously for spheres the normal would just constantly be the vector between sphere centers.

I keep the world normal constant and just re-evaluate the distance. E.g. d = ( p2 - p1 ) * n + d0. This is the easiest of all variations and works great in practice.

Obviously you can go arbitrary fancy. If you rotate the normal you very likely need to take that also in the constraint into account:

C = ( x2 + r2 - x1 - r1 ) * n

dC/dt = ( v2 + w2 x r - v1 - w1 x r1 ) * n + ( x2 + r2 - x1 - r1 ) * dn/dt

dn/dt = omega x n. Obviously you will get a slightly different Jacobian.

Usually the bold term is zero since either x2 + r2 - x1 - r1 is zero or it is parallel with the normal resulting in a zero cross product. Depending on how you setup the non-penetration constraint. During position projection both might not be true anymore necessarily.

Obviously you can go arbitrary fancy. If you rotate the normal you very likely need to take that also in the constraint into account:

No, velocity solving is a linear problem and during iterations you are projecting velocities onto the 'fixed' constraint manifold. Position projection is a non-linear problem (hence the N for Newton in NGS). The Jacobian is a function of the position. I guess you could try to anticipate the change in the position due to the velocity change, but I am not sure if this would makes sense.

Oh, I might have misunderstood your question. In the case of the velocity iteration we still have the following:

dC/dt = ( v2 + w2 x r - v1 - w1 x r1 ) * n + ( x2 + r2 - x1 - r1 ) * dn/dt

But now either x2 + r2 - x1 - r1 = 0 if you setup the local contact points from global ones. Or in the case where the local points are disjoint (I think you asked about this earlier) you compute the normal from the difference vector anyway in which case the cross product ( x2 + r2 - x1 - r1 ) x n would vanish since the vectors are parallel.

This might not matter in practice at all. You have a similar problem with prismatic joints. If you ignore the right term you get problems if the bodies are sliding apart quite a bit. Penetrations are usually small so it might not matter. Formally this would be the correct way though.