Constraint Solver problem

[Bilbert]

I have been working on fairly simple physics simulation with a constraint solver. I've ran into a bit of a snag. I received great help here with the Jacobian matrix and I hope someone can point me in the right direction once again. The paper "Iterative Dynamic with Temporal Coherence" by Erin Catto was great an I'm essentially trying to implement the same system. After writing all the code, it seems to work but I believe the results are wrong. I have simplified the problem to the basics - simulating a single body acting as a pendulum hanging from the origin with a length of 1.0. The object has a mass of 1.0. I'm solving the equation:

Code: Select all

J * M^-1 * J^t  * lamba = -J * (1/dt*V + M^-1 * Fext)

With mass = 1.0 and the body starting at rest, that eliminates the M^-1 and 1/dt*V. I'm using a 6x1 vector to represent Fext; the first 3 values represent force and the last 3 values twist. This is a distance constraint so twist should remain unaffected. Also, I'm keeping the z value 0. My constraint 6x3 Jacobian is:

Code: Select all

J = [ 1 0 0  0 -z  y 
      0 1 0  z  0 -x
      0 0 1 -y  x  0 ]

I got this from "Iterative Dynamics" and it seemed to agree with a Game programming gem article by Russ Smith. I implemented a Conjugate Gradient algorithm described in "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain" by Jonathan Richard Shewchuk. Some test LCP problems confirmed it working. With all the pieces together and running, the initial case simulation I got:

Code: Select all

lamba = [ 0 4.9 0 0 0 0 ]

and

Code: Select all

J^t*lamba = [ 0 4.9 0 0 0 -4.9 ]

If i read this correct, that would apply a force upward with a rotation around the Z-axis. Wouldn't a distance constraint force try to pull the object back toward the origin? That would be a force parallel to the X-axis not Y. I think i have all my pieces in line. If someone could please confirm my method and math concepts are correct, then it may be a bug in code or a problem between monitor and chair. Sorry for the long post. Lengthy problem require lengthy descriptions. Any help or suggestions are very welcome.

Thanks,
Bil

[Dirk Gregorius]

For a pendulum I would use the following constraint (assuming the anchor of the pendulum is at the origin):

C = |x| - L
dC/dt = x/|x| * v

J = ( x/|x| 0 )


It seems that you have formulated a ball-socket joint with an implicit fixed reference body. This will work as well, but will result in different forces. You used:

C = x + r
dC/dt = v + cross( omega, r )

J = ( I -skew( r ) )


This difference is that for the later the orientation of the swinging body will change. The first problem is simple since it is a scalar equation. For your other problem try a direct 3x3 matrix solver to validate the result of the CG solver or test if it is an inverse of the effective mass (J*W*JT).

[Bilbert]

Yes, I intentionally use a ball-socket joint. My first milestone is to simulate a chain with links of different mass connected with ball-socket joints. My thinking was that limiting the solution to only 1 mass would make it easier to verify results and then build from there. I did test my CG algorithm on 2x2 matrices with success but having pieced it together from multiple articles I'm not confident with it. I couldn't determine, however, if i was trying to solve the wrong equations or was incorrectly solving them. From your response it would seem my equations are correct. I shall look back into my CG solver. Maybe a Gauss-Seidel algorithm would be easier to implement?

Thanks for the help,
Bil

[Bilbert]

So after much diving into LCP, my constraint program is working...kinda. When the mass is initially placed above or below the origin it holds the constraint correctly. Placed to the left or right it seems to "try" to hold the constraint but then drifts apart. I think this is from numerical drift. Many papers describe a bias to be added to the system to correct for this. Is this simply a spring force pulling it back into correct position?

Bil

[Erin Catto]

You can think of Baumgarte stabilization as a spring + damper. However, if used properly, a Baumgarte stabilization can be stiffer than a spring while being stable.

[Bilbert]

Thank you Erin. I looked into Baumgarte stabilization more and your explanation makes sense.