I am glad it works, but personally I would be suspicious and assume I might potentially miss a transform somewhere.

I'm almost sure this code works fine. The rigid body behaves quite predictably and very similar to simulation I made in Unity.

I could share my code if anyone wanted to play with it.

I'd also be happy if anyone could look again at the original code, which works, and its other variant that I think should work as well but for some reason doesn't.

Original code:

        oriMat = quaternion.ToMatrix(self.ori)
        inertia_world = self.inertia_local * oriMat
        matrix.Invert(inertia_world)
        angAcc = vector.Transform(self.totalTorque, inertia_world)

To recap: oriMat is effectively local-to-world transform. self.totalTorque is torque in world space. So this code transforms local inertia to world inertia and uses its inverse along with world torque to compute world angular acceleration.

Modified version:

        oriMat = quaternion.ToMatrix(self.ori)
        self.totalTorque = vector.Transform(self.totalTorque, matrix.Inverted(oriMat))
        angAcc = vector.Transform(self.totalTorque, matrix.Inverted(self.inertia_local))
        angAcc = vector.Transform(angAcc, oriMat)

Here I wanted to do the torque-inertia transform in local space. To do so, I transform self.totalTorque to rigid body's local space using inverse of oriMat. Then I mul that by inverted local inertia to get angular acceleration in local space. Finally, I transform angular acceleration to world space by transforming it with oriMat.

For some reason this variant doesn't yield correct behaviour of my rigid body and I really can't tell why. Anyone?

I was pointed out by a friend that in the modified version I should skip the last line. So the coude should look this:

        oriMat = quaternion.ToMatrix(self.ori)
        self.totalTorque = vector.Transform(self.totalTorque, matrix.Inverted(oriMat))
        angAcc = vector.Transform(self.totalTorque, matrix.Inverted(self.inertia_local))

And indeed that works. I wrote down on paper all these transforms and this code gives the exact same series of transforms as the original one.

Original:

angAcc = torque * inertia^-1 = torque_world * (inertia_local * localToWorld)^-1 = torque_world * localToWorld^-1 * inertia_local^-1

Clearly this gives angAcc in world-space what I need.

The modified new:

angAcc = torque * inertia^-1 = (torque_world * localToWorld^-1) * inertia_local^-1 = torque_world * localToWorld^-1 * inertia_local^-1

This I would expect to give me angAcc in local-space (which further needs to be mul'led by localToWorld to be in world space) but math clearly shows this is the same as the previous formula.

Anyone sees what is wrong with my understanding? :)

Just bumping up the thread. It's my recent post that I would love to know the answer for.

The mistake is that inertiaWorld != localToWorld * inertiaLocal

The reason is that inertia is a tensor and not a matrix and that they transform differently. Hence the physics definition of a tensor: "A tensor is what transforms like a tensor".

You have to build your example the other way around:

localAngularAcc = localInertia^-1 * localTorque

We can, of course, compute the body space torque like:

localTorque = localToWorld^T * worldTorque

Now we substitute equation (2) into equation (1) which yields:

localAngularAcc = localInertia^-1 * ( localToWorld^T * worldTorque )

Then we transform the local angular acceleration to world space

worldAngularAcc = localToWorld * localAngularAcc = localToWorld * localInertia^-1 * ( localToWorld^T * worldTorque )

Finally we define the world space inertia as:

worldSpaceInertia^-1 = localToWorld * localInertia^-1 * localToWorld^T

This yields the pleasing result:

worldAngularAcc = worldSpaceInertia^-1 * worldTorque 

The key point is really that inertia is a tensor and not a matrix. Note how all my operations are conceptually applied to a vector and not to the tensor. I find that confusing as well, but I am not aware of a better way to explain it! Hopefully the example makes clear how this works for the inertia tensor and how your math above is wrong missing a transformation.

Maybe you should think of it more as a change of basis. You transform the world torque into the local body space. Then you apply the local inertia tensor. And finally you transform the the result back to world space.

HTH,

-Dirk

I found this link. It explains all different transformations.

http://www.continuummechanics.org/cm/coordxforms.html

Sorry for such a massive delay but I completely forgot about that thread...

Okay, I've gone through your equations and the link you provided but it's still not right.

You're saying (and the said link) that to transform tensor you go like transform * tensor * transform^T. But my code doesn't do it. It simply does transform * tensor and it works (in world-space). Look again at the code that computes world-space angular acc:

        localToWorld = quaternion.ToMatrix(self.ori)
        inertia_world = self.inertia_local * localToWorld
        angAcc_world = vector.Transform(self.totalTorque, matrix.Inverted(inertia_world))

Expanded and written in math notation it would be like:

angAcc_world =
self.totalTorque * inertia_world^-1 =
self.totalTorque * (self.inertia_local * localToWorld)^-1 =
self.totalTorque * localToWorld^-1 * self.inertia_local^-1

Again, localToWorld is only used once, not twice. And this code really works, he simulation behaves in predictable manner, as opposed to any other approach I try.

I don't really know where this:

Finally we define the world space inertia as:

worldSpaceInertia^-1 = localToWorld * localInertia^-1 * localToWorld^T

is taking place in my code.

One thing interesting here is that the last form of this equations could be written like this:

angAcc_world =
self.totalTorque * localToWorld^-1 * self.inertia_local^-1 =
(self.totalTorque * localToWorld^-1) * self.inertia_local^-1 =
totalTorque_local * self.inertia_local^-1 =
angAcc_local

Now that is darn confusing. Does it really mean that angAcc in world-space is the same as angAcc in local-space?

Every math and physics book you will find will tell you that this is how it works, E.g. here:

https://en.wikipedia.org/wiki/Moment_of_inertia#The_inertia_matrix_in_different_reference_frames

Just because something looks correct, doesn't actually proof it is correct. Secondly, you essentially want to convince me here that you are right and the whole physics world is wrong (even though I have shown you the mistake several times backing it up with external resources). My suggestion to approach your problem would be to apply the correct transformation and now that you can rule that out as an error source find your real bug. Or maybe there is no bug and you have a wrong expectation what the outcome of your experiment is.

And obviously:

angAcc_world == angAcc_local would only be correct iff local_to_world is the identity matrix. You have essentially proven yourself that your are wrong.

If you do Transform * Tensor the results will be *hugely* incorrect. Imagine if the body is almost symmetric (Tensor is almost a unit matrix times mass) - you will effectively multiply your acceleration by the transform matrix. So if you rotate the body 180 degrees around the X axis, all new torques will be completely invertex around that axis.

That is a nice example! Let me quickly extend on it. Say we actually have a sphere inertia which is invariant under rotation.

I_world = R * I_local * RT

The sphere inertia I_local = 2/5 * mass * radius^2 * identity_matrix. Let's write this as I_local = s * E (where s = 2/5 * mass * radius^2 and E = identity matrix). Then

I_world = R * ( s * E ) * RT = s (R * E * RT) = s (R * RT) = s * E

Which is the correct result because of the symmetry. You formula breaks apart here.

It's not that I'm stubborn and decided to proof that the whole world is wrong and I'm right. But I re-checked my simulation nubmer of times and know it behaves right. And I discovered why...

I thought that I would proof myself right by implementing the simulaton in Unity so that it's easier to compare with Unity's built-in rigidbody implementation (that one is for sure working right). I typed out the same simulation code I did for my Python app. And... in Unity it didn't work right. So I used the actual formula you're suggesting for transforming the inertia tensor. And it worked! Now that was a discovery... It pushed me to write it in yet another language, this time C++ (I thougt maybe I had screwed and/or misunderstood how objects in Python worked). And my simulation worked identically to the one in Python.

Since I'm stubborn to keep insisting on using row-major matrices order (and Unity uses column-major) I took a quick look at the function that adds vector to quaternion (to alter orientation with new angular velocity). Here it comes:

    quat operator + (const vec3& v) const
    {
        quat temp = quat(0.0f, v.x, v.y, v.z);
        temp = temp * *this;
        temp.w *= 0.5f;
        temp.x *= 0.5f;
        temp.y *= 0.5f;
        temp.z *= 0.5f;
        
        return quat(w + temp.w, x + temp.x, y + temp.y, z + temp.z);
    }

I took a blind shot and swapped the multiplication order in this line:

        temp = temp * *this;

to this:

        temp = *this * temp;

and surprisingly it worked!

To quickly recap:

When I use my original add-vector-to-quaternion function and inertia_local*local_to_world the simulation works.

When I use the new add-vector-to-quaternion (with quaternions mul swapped) and the "valid" inertia transformation the simulation works.

Why the wrong ordering of quaternion multiplication could compensate for the lack of the localToWorld^T matrix?