It seems you''re trying to implement implicit integration for rigid body simulation. The thing is, that this is not very useful, because the forces rigid bodies get are rarely dependent on the other rigid bodies. Therefore, as there''s very little interaction between the bodies, the matrix ends up really sparse.

The backward Euler method is useful on systems where we have a great number of objects (usually particles) which are connected somehow (for example springs). Cloth simulation is a great example of this. The matrix will still end up being sparse (on the average 8 nodes on each row).

You derived your matrix correctly. You''re correct, that kind of "does the quantity we''re derivating depend on the quantity with which we''re derivating with respect to?" is an intuitive way to understand the sparsity of the generated matrices.

I don''t know what you''re thinking, but the implicit Euler method is definitely NOT to be applied on the bodies one by one. The whole idea of the method is that it supports GLOBAL interactions between the objects. Therefore we need to construct one, huge but luckily sparse matrix.

- Mikko Kauppila

Well, my idea was to understand the whole process for a single rigid object and then extend it for more, but I wasn''t aware that in this case you need to take into account all interacting objects.

I planned to use this method for single rigid bodies, multiple rigid bodies linked by joints (articulated), cloth and maybe other things I find I can do with this method.

The reason why I decided to explore this method is because my current integration methods fail under certain circumstances, especially under the influence of springs.

I looked briefly into LCP based methods and read another of Baraff''s papers ("Linear Time Dynamics using Lagrange Multipliers") and it just went flying straight above my head so I thought I would tackle something else.

Ok, I''m not going to give up on this one though.

Let''s take the example of cloth and consider three particles then (I assume this could be a minimum number so I can understand the process without making it too complicated).

So because we are considering particles then our state vector is reduced as there is no more orientation component:

	X(t) =	[ x(t) ]	position		[ P(t) ]	linear momentum 

The derivative of the state vector is

	.	X(t) =	[ v(t) ]	linear velocity		[ F(t) ]	force 

If we consider three particles such as particle A is linked to B, and B to C, but A is not linked to C

	A----B----C 

We can put all the state vectors into a joint state vector and the same for the derivative vector

				.	X(t) =	[ A.x(t) ]	X(t) =	[ A.v(t) ]		[ A.P(t) ]		[ A.F(t) ]		[ B.x(t) ]		[ B.v(t) ]		[ B.P(t) ]		[ B.F(t) ]		[ C.x(t) ]		[ C.v(t) ]		[ C.P(t) ]		[ C.F(t) ] 

The Big Matrix (which I''ll call BigM - was A previously) will now be like a grid where all the particles cross.

I''m just thinking out loud here so I may be completely wrong in which case I''d appreciate if you could correct me.

BigM = [	AA	AB	AC	][	BA	BB	BC	][	CA	CB	CC	] 

where the particles are not interacting then those components will be zero.

BigM = [	AA	AB	0	][	BA	BB	BC	][	0	CB	CC	] 

For the diagonal components then we have the previous A matrix from the previous post.
For example the AA matrix would be:

AA =[	1 	d( A.v(t) )	][	-	-----------	][	h	d( A.P(t) )	][				][ 		    1		][	0	    -		][		    h		] 

For the components where the particles interact then we should have something like this.

AB = [   1 	d( A.v(t) )	d( A.v(t) )	][   - -	-----------	-----------	][   h	d( B.x(t) )	d( B.P(t) )	][					][	d( A.F(t) )	1   d( A.F(t) )	][	-----------	- - -----------	][	d( B.x(t) )	h   d( B.P(t) )	] 

So at this point we could think "is the velocity of A dependant on the position of B?" and maybe some of these would be quite obvious but maybe that depends on the type of relationship that the particles have.
I remember reading that there is different types of springs for the cloth models and maybe depending on the type then this matrix could be different (just guessing here).


Whether I''m right or not I will still need to calculate expressions like:

d( A.v(t) )----------- = ?d( A.P(t) )d( A.v(t) )----------- = ?d( B.P(t) ) 

I have a feeling that these may have to be implemented as a case by case, because as I said above they may be different for the different types of spring (I think).
There must be documentation that explains this, or maybe Baraff''s paper explains it and I just didn''t get it. I''ll have to read it again.

Thanks again Mikko.
lusinho

>There must be documentation that explains this,
>or maybe Baraff''s paper explains it and I just
>didn''t get it. I''ll have to read it again.

Baraff''s paper doesn''t explain construction of the matrix very well (quite badly to speak the truth).

>X(t) =
>[ x(t) ]position
>[ P(t) ]linear momentum

It''s of course possible to derive the equations this way by constructing a state vector which is a combination of positions and linear momentums. I haven''t implemented it myself. It''s slower and more complex. We get off MUCH more easily if we state in the original implicit Euler formulation:

deltaX = (v+deltaV)*timeStep
deltaV = f(x+deltaX,y+deltaY)

From this we only solve deltaV. deltaX is then trivial to compute. Computing deltaV only results in a simpler and smaller matrix. This is explained really well in Baraff''s paper (Large Steps [98] and the one from "physically based modeling principles and practice".

>where the particles are not interacting then those components >will be zero.
>
>BigM =
>[AA AB 00]
>[BA BB BC]
>[00 CB CC]

This is not just a bit but completely right! You''ve understood the method this far perfecly. In cloth simulation the final matrix is a sum of three matrices, A=M-h*dFdV-h^2*dFdX (they all have the same sparsity patterns, so summing them is easy). Let''s concentrate on the dFdX matrix (force derived wrt. position) which is the only "compulsory" matrix (dFdV is not compulsory if you have no damping forces in the first place, which I strongly recommend). In this case, AA in BigM means "force exerted on A derived wrt. position of A" and so on.

>d( A.v(t) )
>----------- = ?
>d( A.P(t) )

The actual derivatives are the hardest ones to compute. By solving only deltaV all we need to compute is derivatives of forces wrt. positions (and with damping forces forces wrt. velocities). External forces (mouse interactions) don''t affect the derivatives. Only if the force depends on the positions/velocities of the particles do they affect the derivatives. The only forces of this kind are the spring forces (I personally handle air friction by multiplying velocities with 0.999 at the end of the step

So in the end, the only problem left is to derivate the spring forces. The formulas are simple, and it''s shown in this paper: Kwang-Jin Choi, Hyeong-Seok Ko: Stable but responsive cloth. ACM Transactions on Graphics 21(3): 604-611 (2002)

- Mikko Kauppila