This is in reference to "Stable Constrained Dynamics": https://hal.inria.fr/hal-01157835/document

Equation (18) / (21)

I'm having trouble understanding how to build this K "geometric stiffness" term.

K = (∂J^T / ∂x) λ

Where J is the constraints jacobian and λ is the constraint force magnitudes.

What I do know - based on its usage in (21), K should be a (3n x 3n) matrix in 3D where n is the number of particles lets say.

What I'm confused about - the jacobian J is a (C x 3n) matrix where C is the number of constraints. λ is (C x 1). This doesn't seem to work out in terms of the matrix dimensions...

What am I missing here?

If I consider only a single constraint, then it does appear to work out in terms of the dimensions - I end up with λ being a scalar and K ultimately being (3n x 3n). However, that leads to the question of how to then build K such that it contains all of the individual constraint K's (one K for each constraint I guess)?

The Jacobian is the matrix of transposed gradients for each constraint function. The partial derivative of the gradient is the Hessian matrix. In the multidimensional case you have a Hessian for each constraint. So dJ/dx would be a tensor of Hessians (array of matrices) I guess. 

In the one dimensional case K = H * lambda. Extrapolating from here I would guess in the n dimensional case it is the linear combination of Hessians. E.g. K = H1 * lambda1 + H2 * lambda2 + ... + H_n * lambda_n

This is a good summary of the calculus used here:

https://www.value-at-risk.net/functions/

Yeah that sounds right. 

I emailed the author and asked for clarification on this. He confirmed that the final K matrix is indeed the linear combination of hessians. I've implemented it this way into my physics engine and it seems to be working as described in the paper.

I found matlab to be of great help when deriving K. For anyone else that gets stuck, here is an example of computing K for a simple length constraint between points "A" and "B" with rest length "l":

syms A1 A2 B1 B2 lambda l
A = [A1; A2];
B = [B1; B2];
phi = norm(A - B) - l;

J = jacobian(phi, [A1, A2, B1, B2]);
K_til = jacobian(transpose(J), [A1, A2, B1, B2]) * lambda;
pretty(K_til)

The same logic should work for any type of constraint.

Cool, thanks for confirming this!