I've finally got my mini engine to work (with some vibrations but oh well), but I really need to add angle constraints to my ragdolls and I just can't figure out how they fit in the contact response calculations.

The main issue is that it is impossible to restrain rotation by only applying force/impule to one point. Currently I'm solving for impulse that would give the two bodies the same velocity at the contact point. If I split up the velocity and angular momentum response, then I end up with too many unknowns to solve for. And that's not even touching interpenetration handling. Also, if the constraints are given in euler angles, is it possible to express them in terms of quaternions to simplify the calculations somehow?

There seems to be an awful lack of materials on this topic. 99% of the online sources cover the very basic stuff or refer to off-the-shelf engines, and I can't just randomly buy every book and hope that it contains what I'm looking for (okay reading contents helps, but last time I tried that I ended up being wrong).

Also I have an issue that is possibly related to the original problem - if I swing a flail, I expect the ball to follow a large arc. However, in my simulation it just drags after the handle. I suppose forcing the chain to remain within a certain angle from the straight direction would fix it somewhat, but I might be missing something crucial. The flail handle is controlled by predefined animations, so I have to simulate the rest of the flail ragdoll based on its position/orientation. Currently I'm just taking the velocity of the tip (position change since the last frame) and plugging it into the response equations. However, it would seem that in order to get the result I need, the acceleration (or orientation change) would have to be factored in as well somehow, but none of my equations use accelerations anyway so I don't know where to start.

Okay I think I kind of figured it out - we can separate the impulse response into the linear and angular parts, and the two constraints are (1) conservation of linear momentum (P1+P2=0), and (2) conservation of angular momentum (using contact point as reference, L1 - c1 x P1 + L2 - c2 x P2 = 0). The joint itself requires that relative velocity at contact point is zero (which involves impulse response matrices).

Next, it seems that if we have 3 angle constraints, which correspond to the 3 components of relative angular velocity in the contact space. For active components (those that are in the extreme position), we allow the corresponding component of L1 - c1 x P1 to be non-zero, and instead require the resulting relative angular velocity to have a zero in that component.

Still not sure what to do with position (interpenetration) resolution, and how to express the constraints (originally in euler angles) in terms of quats to check the angle conditions. I'd imagine it should translate into something involving quat dot products?

The are two problems here. First we need to build a spherical joint (ball-socket) and then we apply a conical limit. Do you know how to model the spherical joint?

Then when you have the joint constraints the solver simply iterates over both joints and contacts. E.g.:

for ( n iterations )

SolveJoints()

SolveContacts()

If by spherical joint you mean locking a point on two rigid bodies together, then yeah I got it to work already. Are you saying I should model the position and angle limits separately? I was thinking that it is actually pretty easy to implement a generic 6-degree constraint based on what I've laid out above, which can model pretty much every kind of joint I might ever need. I also throw all collision contacts and joint contacts together (its in Javascript so I don't even need a base class).

By the way, a little off topic, is there any merit to setting a resolution factor for every iteration so that it actually has a chance to converge in tricky situations?

You can use such a factor. This is called Successive Over/Under Relaxation. I made no good experience with this.

The classical approach to model joints is:

1) Define the position constraint

2) Build time derivative of position constraint to identify velocity constraint

3) From velocity constraint identify Jacobian by inspection

Don't try to build the gradients of the position constraint directly!

For the spherical joint with a conical limit you need 4 constraints. 3 to model the ball socket and 1 for the limit. So the question is how to model the conical limit. In order to define angular constraints between two rigid bodies you need to define a local joint frame on each body. Then you can inspect the relative frame to compute angles From there you need to define a position constraint to model the joint limit. I recommend not using Euler angles. The Havok documentation has a good explanation of this.

All I describe here has been discussed many times before. I recommend looking at the presentations from the GDC physics tutorial and Erin Catto's presentations in particular. They can be found here:

www.box2d.org

Thanks for the link;

What do you mean about position constraint gradients though? Position change per unit impulse? That's exactly what I'm doing right now and it seems to work pretty well, what's wrong with that? (I mean the matrix that looks like -r_x * I^(-1) * r_x + 1/m). I'm coming up with a 6x6 jacobian for a generic joint.

Yeah I guess I didn't mean euler angles per se; I suppose taking the relative frame quat and turning it into a scaled vector would be the thing I'm looking for, assuming I have ranges like -pi/4..pi/4 as my constraints?

The Jacobian should be 6 x 12. Then the effective mass Me = ( J * M^-1 * JT ) ^-1 will be 6x6 (e.g. 6 x 12 * 12 x 12 * 12 x 6) for a generic joint.

These are all well understood and solved problems. If you want to come up with something on your own this fine, but I am not sure if I can help you with this. Sorry!

I've looked at Bullet's code, they do specify the angular limits in euler angles.

Let's say I've built the 6x12 Jacobian. Now I have several options on how to solve it:

1. Solve every row individually, as if they were separate constraints, and clamp the lambdas accordingly. Seems good, but might not be very stable if the iteration count isn't high enough? I think bullet uses this approach.

2. Solve it as a system, then clamp the result. Probably incorrect, as I'm basically treating LCP as a linear system?

3. Solve equality constraints as a system, then solve inequalities individually. This should work better, right? I.e. if the position limits are usually min=0 max=0, then I could just solve it as a 3x3 matrix, just like a regular point-to-point joint. That seems like a lot of 'if' statements though, plus I don't feel like inverting a 6x6 matrix.

For each joint we need to define a body-local joint frame. Here is a sketch that hopefully illustrates the idea:

Of course you want the two local anchor point to be coincident using the point-to-point constraint. From what I understand you already figured that out, so I skip this for now, but we can talk about later if you like or need another example of what I will show now below.

Now we need to model the conical limit. A very simple way to do this would be using the following constraint:

C = u2 * u1 > cos( theta / 2 )

Here theta is the (full) cone angle and * denotes the dot product. From there we follow the basic procedure and build the velocity constraint:

dC/dt = du2/dt * u1 + u2 * du1/dt

The change of a fixed axis in a moving frame is simply omega x axis. Hence

dC/dt = ( omega2 x u2 ) * u1 + u2 * ( omega1 x u1 )

Using the scalar triple product identity and arranging terms should then yield:

dC/dt = ( u2 x u1 ) * ( omega2 - omega1 )

From the velocity constraint we can now identify the Jacobian by inspection:

J = ( 0^t -(u2 x u1)^T 0^T (u2 x u1)^T )

I used the ^T here to make clear that this is a 1 x 12 matrix.

So in the solver we just check the position constraint. If it is satisfied we can ignore this constraint in this frame. Otherwise we need to solve:

J * W * JT * lambda = -J * v - beta * C / dt

This approach has a couple of advantages. Using the time derivative and utilizing that dC/dt = del_C/del_x * dx/dt + del_C/del_t is much easier than trying to build the gradients directly. Though a good geometrical understanding helps a lot to simplify constraints and Jacobians. You also have the stabilization term naturally handy and don't need to guess it like when you would start from the velocity constraint. BTW, the fancy formula is called the total derivative and a summary can be found here: https://en.wikipedia.org/wiki/Total_derivative

A little implementation gotcha. In your implementation you would store the constraint axis local in each body frame and then u1 = R1 * u1_local and u2 = R2 * u2_local to identify the current joint frames. Then you basically simply inspect the relative orientation.

This yields a 4x4 system when solved together with the point to point constraint. You have a couple of options now:

1) Solve each scalar constraint individually (this is essentially what the ODE dQuicksolve does)
2) You solve the point to point constraint as 3x3 block constraint and the limit as a scalar constraint. This is the common way these days from my experience
3) You solve the point to point constraint together with conical scalar constraint. This yields a 4x4 MLCP which can be solved with direct enumeration

Erin does the 3) in Box2D. This yields the best results as it boost converges since the point to point constraint and angular limit don't fight each other. In 3D I found this very expensive so I use 2) personally in Rubikon.

Since in an iterative solver the order in which we solve constraints matters the remaining question is if I solve the point to point constraint before the limit or vice versa. I solve the point to point constraint always last to emphasize it. I found that stretch at anchor points is more visually noticeable than a violated constraint.

PS:
Here is an example for Direct Enumaration from Murphy's webbook:
http://ioe.engin.umich.edu/people/fac/books/murty/linear_complementarity_webbook/kat1.pdf