Incorrect rotation for 'thin/long' objects

[Norbo]

This is a cross post from the Math and Physics forum at gamedev.net (http://www.gamedev.net/community/forums ... _id=470954) in an effort to get some more exposure, but here goes:

In my physics engine, objects which are not of approximately the same height/length/width (cubelike) develop inaccurate angular velocities.

At first, I had assumed I had made an error in my numerical method for finding the inertia tensor and continued working within the bounds of cubelike objects, planning to fix it at a later stage. Now, as I add support objects whose inertia tensors can be evaluated symbolically (spheres, cones, cylinders), I observe the same behavior despite my confidence in the calculation.

For testing purposes, I focused on the cylinder case. I calculate the body inertia tensor in the same manner as listed at http://scienceworld.wolfram.com/physics ... inder.html, with a 'minor' difference. As I am developing specifically for the XNA Framework, I am using their Matrix class which is exclusively 4x4. For my implementation, I simply used the 4x4 identity matrix as a base, setting the M11, M22, and M33 fields to the appropriate values (as shown in the link).

The body inertia tensor inverse is then calculated based on that, and each frame's inertia tensor is calculated using:

Code: Select all

myInertiaTensorInverse = rotationMatrix * myBodyInertiaTensorInverse * Matrix.Transpose(rotationMatrix);

Video:
(removed)

I filmed a test set up with two cylinders of different sizes and the same mass under two different timescales (1 and .25). In each, the iteration count for the iterative solver is very high (100) so the impulse calculation converges accurately. Additionally, the larger, more cubelike cylinder stabilizes quickly after rolling around a little bit, while the smaller cylinder jitters noticably at timescale = 1 as it bounces and attempts to flip around. In both timescales, significant rotational error can be observed after the cylinder falls off the platform. As an extra note, the larger cylinder and other similarly cube-like objects experience no such noticable behavior upon falling off the platform.

What is the most likely culprit for the wildly accelerating deathspin after falling from the platform, and the related hyperactivity of the small cylinder on the platform? So far, there's a fairly short list of options I'm currently looking at:

* The inertia tensor representation. As I am not 100% sure, I would like to verify that my 4x4 inertia tensor matrix interpretation is valid.

* My current method of moving the simulation forward is a very simple euler integration implementation; as some of the jitter was removed with a lower timestep (more updates per movement), I thought perhaps that changing to an RK4 or other higher level integration method would help alleviate the issue, though I'm not sure it can fully explain the extreme case of spinning after falling off the platform.

* Another option would be that I have misused the inertia tensor somewhere else in the implementation, outside of contact impulse application (given that it spins inappropriately even when not in contact with any object). This narrows it down to essentially a couple of statements:
-As listed above, the inertia tensor inverse update in the forces/velocities update method:

Code: Select all

myInertiaTensorInverse = rotationMatrix * myBodyInertiaTensorInverse * Matrix.Transpose(rotationMatrix);

-And also in the forces/velocities update method:

Code: Select all

angularVelocity = Vector3.Transform(myAngularMomentum, myInertiaTensorInverse);

* I don't currently have separate bias velocities/impulses in position correction, but the spinning problem can be observed after leaving the platform leading me to believe that the primary culprit is elsewhere. This may, however, contribute to some of the jitter visible by the small cylinder in the first simulation of the video.

* I use a quaternion and matrix method for storing rotations, similar to David Baraff's approach on pages 21 and 22 in the paper at http://www.cs.cmu.edu/~baraff/sigcourse/notesd1.pdf, as follows:
First, I take into account forces to generate the new velocities.

Code: Select all

//Angular Momentum:
myAngularMomentum = myAngularMomentum + myTorque * wait;
//Tensor:
myInertiaTensorInverse = Matrix.Transpose(rotationMatrix) * myBodyInertiaTensorInverse * rotationMatrix;
//Angular Velocity:
angularVelocity = Vector3.Transform(myAngularMomentum, myInertiaTensorInverse);

Second, any constraints/contacts are resolved, and finally the rotation is updated:

Code: Select all

rotation += ((new Quaternion(angularVelocity * wait, 0)) * .5f) * rotation;
rotation = Quaternion.Normalize(rotation);
rotationMatrix = Matrix.CreateFromQuaternion(rotation);[/list]

I'd greatly appreciate any input about anything I'm doing wrong or if there are other possibilities I've overlooked.


*Edit:
Upon creating a more controlled environment for demonstration purposes, I filmed the following video:

(removed)
Left to right, the momentums provided were:
(1,0,0) (0,1,0) (0,0,1)

Having observed the extremely rapid spin of the thin cylinder under the last angular momentum setting and the equal angular velocity between the first two, it would seem likely that I've reversed the inertia tensor calculation's axes. This could be due to XNA/DirectX's interpretation of Y as 'up', and somewhere in my construction of the tensor I could have used an inconsistent reference.

This is very promising and I'll investigate after I wake up. Perhaps my troubles have come to a premature end

[melax]

As you mentioned, first ensure that you get the little bugs sorted out (y up vs z up etc).

Let me help you with the next issue you may face next...
It looks like you are conserving angular momentum - which is the correct thing to do. Some physics engines dont bother and just keep the spin constant.
After you get all the kinks out of your code, try spinning the cylinder about another arbitrary axis. not just (1,0,0), (0,1,0) and (0,0,1).
When you do this I suspect you will notice that the cylinder will pick up speed (angular velocity) and reorient itself so that its long axis becomes aligned with the axis of rotation. The angular momentum will remain the same.

To solve that, you can use runge-kutta integration just for the orientation update based on the angular momentum of the rigidbody after you have done all that dynamics solving sequential impulses etc etc. Baraff's notes show how to do this for then entire rigidbody state, but that isn't compatable with the way a lot of physics engines update their velocities. Just use the orientation update from those references.

In code this might look something like:

Code: Select all

Quaternion DiffQ(const Quaternion &orientation,const float3x3 &tensorinv,const float3 &angular)
{
	Quaternion s_orientation_normed = normalize(orientation);
	float3x3 s_orientation_matrix = s_orientation_normed.getmatrix();
	float3x3 Iinv = Transpose(s_orientation_matrix) * tensorinv * s_orientation_matrix;
	float3 halfspin = (Iinv * angular) * 0.5f;
	return Quaternion(halfspin.x,halfspin.y,halfspin.z,0)*s_orientation_normed;
}

Quaternion rkupdateq(const Quaternion &s,const float3x3 &tensorinv,const float3 &angular,float dt)
{
	Quaternion d1 = DiffQ(s ,tensorinv,angular);
	Quaternion s2(s+d1*(dt/2));

	Quaternion d2 = DiffQ(s2,tensorinv,angular);
	Quaternion s3(s+d2*(dt/2));

	Quaternion d3 = DiffQ(s3,tensorinv,angular);
	Quaternion s4(s+d3*(dt));

	Quaternion d4 = DiffQ(s4,tensorinv,angular);
	return normalize( s + d1*(dt/6) + d2*(dt/3) + d3*(dt/3) + d4*(dt/6)) ;  
}
...
  // in the update code when its time to compute the orientation for the next frame...

          rb->orientation_next  = rkupdateq(rb->orientation,rb->tensorinv,rb->rotation,DeltaT);

The code may not be optimal, but I think its basically correct.
The supplied inertia tensor rb->tensorinv is not rotated - its based on the object being aligned with the world axis.