[ODE] ODE Source code vs. Baraff

[david@csworkbench.com]

Yeah, my last attempt at a solution tried this exact approach... it
calculates the inertia matrices for each body then calculates A, rhs,
lamda (with dSolveLCP), applies the forces, and MoveAndRotate's the
bodies, all inside the iteration loop.  I tried moving A and the inertia
matrices out, but as those do change a little bit after moving the bodies,
the inaccuracy overpowered the speed increase by forcing me to raise the
number of iterations to keep things stable.  I was able to take the
Jacobian calculations out, though, which actually improved the results.

Joakim:  That method uses very little space, since it only stores one 6x6
max block of the system matrix that used to take number of DOF's removed
from the entire system (m) ^2 space.  It still stores the Jacobian
vectors, though, but that's just an mx8 matrix, and worth it.  If you
wanted a really memory miser algorithm, you could recalculate J for each
joint at each iteration, making the integration step actually end up
taking O(1) space instead of O(n)....  It would be simple to make that a
compiler switch if you're interested.  On the other hand, the algorithm
I'm working on now takes a max of 2 * m * 8, so we're still talking O(n)
space, just twice as much.  You'll see great memory advantages from any of
the solutions I've tried for sufficiently large systems.

I'm just trying to figure out if there's a way to keep all that out of the
iteration loop and limit it to actually finding a solution for one full
step as quickly as possible.  If there's a way to use only the forces
added and optimize them iteratively, it would blow away the
integrate-every-iteration approach for speed, only double the memory
required (by precalculating and storing A for each joint), and possibly be
more stable/accurate due to actually attempting to converge on a correct
system instead of just hoping that the error introduced is small enough to
be taken out next timestep. It also gives you a chance to leave the loop
early if last iteration added a small enough amount of force.

David

> (S4) Integration
>
> In fact, the simple solution we have outlined in past e-mails is
> quite similar to your scheme. In fact, the same ODE source code gives
> clues to formulating the "iterative" algorithm.  Using  dSolveLCP
> simplifies a lot and also (because already performs some sort of
>  error reduction) limits the number of iterations until convergence.
>
> Step S1 (PrepareLCP): The code for computing A and RHS for a single
> joint can be  adapted from the existing ODE source code for computing  A
> and RHS with (multiple) joints.  Note that the single joint A matrix
>  is equivalent to computing the diagonal blocks of the multiple
> joints A matrix.
>
> Steps S2 (SolveLCP) and S3 (Acumulate): Don't need any adjustment at
> all.  Simply pass the right parameters for the single joint scenario.
>
> Step S4 (Integration): Don't forget to integrate using MoveAndRotate
>
> Also, I found useful to precompute the inertia matrices first (probably
> along with the A matrices)
>
> Sergi