I am trying to implement joint constraints for connected rods. I was looking at the derivation for the particular Jacobian on this site:

http://myselph.de/gamePhysics/equalityConstraints.html

They chose the constraint to be C = (Pa - Pb)*(Pa - Pb), which makes the jacobian [-d -r1Xd d r2Xd]

d=current distance between constraint point on each object

r=lever arm from center of mass of each object to its constraint point

They also have a javascript demo at the bottom which appears to work. My question is, why does this work. Looking at the math, it looks like it should not.

Because:

When the constraint is valid, Pa - Pb should equal 0. Then the jacobian is 0, so JV is 0 as well, therefore lambda is 0 whenever the constraint itself is valid, even if the derivative of the constraint is invalid. It won't attempt to create any corrective velocity until Pa and Pb are no longer equal.

This is different from a pendulum constraint where Pa - Pb is supposed to be non zero, allowing lambda to be nonzero even when the constraint position is valid.

Is there something I am missing, or is this formula for joint constraint wrong? It looks like the only reason their demo works is because the B term fixes it on the next frame.

I can't find any other decent explanations of how this is supposed to work other than the classic pendulum or penetration constraints.

Your observation is absolutely correct. Choosing a constraint function such that the Jacobian is not defined if the constraint error is zero (C = 0) is a very poor choice. For exactly the reasons you describe. Usually a distance constraint is formulated like C = |pb - pa| - L. If L is reasonable large the Jacobian will be fine. Also note that the correct constraint to hold two points coincident is C = pb - pa. The later is 2d/3d constraint (hence the bold notation).

The reason is that the 1d distance constraint is not only a poor choice because it is undefined when C = 0. Let's assume we actually had a small displacement and get a valid Jacobian. The constraint can only handle forces along the offset vector of the two connected points. If you apply a force orthogonal to that direction the constraint will not handle this as expected and the objects will drift apart until the constraint will handle this in the next frame.

It looks nearly correct. I think you just have the signs and dimension a bit mixed up. Here is how to derive the constraint.

// Position constraint

C = pb - pa

// Velocity constraint (is the time derivative of the position constraint)

dC/dt = pb/dt - pa/dt

= vb + wb x rb - va - wa x ra

= vb - rb x wb - va + ra x wa

= vb - skew( rb ) * wb - va + skew( ra ) * wa

The velocity constraint has a special form:

dC/dt = del_C / del_x * dx /dt + del_C / del_t = J * v + del_C / del_t

If the velocity constraint doesn't depend explicitly on time (del_C / del_t = 0) - which it doesn't in our case - this simply becomes.

dC/dt = J * v

So finally we can find the Jacobian by inspection:

J = [ -1 skew( ra ) 1 -skew( rb ) ]

There are a couple of things here. v the the linear velocity at the center of mass and w (read omega) is the angular velocity. r are the offset vectors (lever arms) from the center of mass to the constrained points and skew( r ) is the cross product matrix such that r x w = skew( r ) * w. Also note that '1' is not a vector, but the identity matrix. The Jacobian is a 3 x 12 matrix. If you build the effective mass matrix Me^-1 = J * W * JT the result will be a 3x3 matrix. You can solve the constraint as 3D block using e.g. sequential impulses. Also for clarification del_ is the *partial* derivative. If you are interested in the math check here: https://en.wikipedia.org/wiki/Total_derivative

There were plenty of great talks about this at the GDC Physics tutorial over the last years. In particular I would check out the talks by Erin Catto and Richard Tongue. You can find all presentations here: http://box2d.org/downloads/

Code for solving the velocity constraint as a 3D block using sequential impulses.

HTH