I am no expert when it comes to physics programming (for that see Christer Ericson’s blog) but like lots of games programmers I have written a few physics demos using Verlet integration and point masses.
So that’s nothing new, but a while back I came across a really nice approach for iterative constraint solving (which is the common approach in Verlet engines). Usually you would formulate your constraint function and solve it directly, for example the classic distance constraint would look like:
error = restlength - |A-B|
In the case of the distance constraint, fixing this error is straight forward, you simply adjust A and B along their axis by error/2 and -error/2 respectively.
For more complicated constraints things can become difficult, because although it can be easy to calculate the error it is not always obvious how to correct for it.
Generally what you want to do is move along the derivative of your constraint function, ie: minimize the error function using Newton’s method. This works pretty fine and provided you do enough iterations will converge quite quickly. However finding the derivative function itself can be problematic, let’s take the example of an angular constraint function:
Given three points and a rest angle (the angle between the two line segments) this function will the return the current error:
error = restangle - acos(dot(B-A/|B-A|, C-B/|C-B|))
Now you have the error, how do you move the three points to correct for it? Well you need to calculate the partial derivative of the error function for each component of A, B and C. To do this symbolically is a pain in the ass and quickly gets messy, I tried plugging it into Mathematica and got a huge mess before eventually doing it by hand which was very tedious given my rusty calculus.
But here’s another way, if you can’t calculate the derivative directly you can at least approximate it using finite differences. This turns out to be so sweet because all of a sudden you only have to write your error function and you’re done. If that’s not clear here’s an example to hopefully elucidate things:
First of all, calculate the current error, I’ll stick with the angular constraint example which takes three points.
currentError = Error(A, B, C)
Now calculate the finite difference for each component of the three points by evaluating the error at some small step away from the current position, here h is our step size and A,B and C are two-vectors but the same thing applies to other dimensions:
X = {h, 0}
Y = {0, h}
derivA.x = (Error(A + X, B, C) - currentError)/h
derivA.y = (Error(A + Y, B, C) - currentError)/h
derivB.x = (Error(A, B + X, C) - currentError)/h
derivB.y = (Error(A, B + Y, C) - currentError)/h
derivC.x = (Error(A, B, C + X) - currentError)/h
derivC.y = (Error(A, B, C + Y) - currentError)/h
(Note the short hand notation for the differentiation of a scalar function with respect to a vector is the del operator)
Now you just need to scale according to the current error and relative masses of the points and apply each derivative to each point respectively.
So there you have it, you’ve solved the constraint with no more information than a way of measuring the current error. Because this is less direct it is also generally more costly but can be a great for prototyping new or more complex constraints.
For more information see the following:
http://www.beosil.com/download/PositionBasedDynamics_VRIPHYS06.pdf
http://www.bulletphysics.com/Bullet/wordpress/bullet
http://www.gamasutra.com/resource_guide/20030121/jacobson_01.shtml
