Please help me debug this iterative impulse solver (infinite spinning, random movement) (Video examples and code)

It looks to me to be a problem with the location of the center of mass, which causes the collision to behave in a non-realistic way. I would double check that the center of mass is at the correct position, and that any COM→contact point vectors are in the right coordinate system. It may also be an issue with the inertia tensor, if you don't account for the COM correctly using parallel axis theorem.

@Aressera Thank you for the pointers! I'll experiment tomorrow! :-)

@Aressera

The red square in the middle is the center of mass, and it follows the object around. I'll probably check the vector from the contact point to the center. The vector should be in world space, right?

Since the contact points seem to be correctly computed and the center of mass seems to be correctly computed (for ease it's the center of the object) i decided to check out whether the math behind the WorldTransformMatrix was incorrect. But it wasn't.

The way I rotate the inertiatensor by the WorldTransformMatrix doesn't seem to be where it's wrong either… Hmm.. But it does really seem that rotations/torque are being calculated badly.. So far I haven't found the culprit

I'm a bit at a loss. I discovered an error in the computation of the angular correction in the position/penetration resolution part of the solver, but it didn't fix the endless impulse buildup.

(frictionless mode)

I've gone through the code. Last thing I'll try is to update the relative contact position (to COM) between the penetration solver and the velocity solver

My guess is that the bug is in the collision resolution.

On the last video, the box spins along its long axis. Integration respecting inertia looks right.
But when a corner hits the ground, even without friction, the impact should counteract the spinning. But it looks like the opposite: The impact accelerates the spin.
Looks like a wrong sign, wrong order of cross products, or something like that.

Maybe it helps to compare with anther resource beside the book. I had used this back the day: https://chrishecker.com/Rigid_Body_Dynamics​

@JoeJ I think you're right.

To simplify the whole thing i chose a cube instead, and set the inverse inertia tensor diagonals to that of a sphere (0.5 & 0.5 & 0.5).

And without friction the cube still behaves weirdly and accumlutes some velocity

I'm getting more and more certain, that's it's a combination of the torque part of the impulse calculation in relation to perhaps the iterative solving method.

Notice how the box rotates without any forces to apply, but only to a certain degree, then it stops. It's only processing two contactpoints out of the four before the solver thinks there are no changes to the impulses for the last two contacts.

That last rotation part seems to be related to the position correction.

I think I'll try the more simple method in the Game Physics Cookbook and see what results i get.

I think specific test/debug scene setups can help.
Currently you test like playing a game, using random human input. That's needed ofc., but once a problem shows up, it can be worth to make a setup so the problem shows each run, and always in the exact same way. When you change the code, you get specific feedback.

In this case, i would place the box axis aligned a bit over the floor, with some initial angular velocity along the x axis. Then the box collides each time with the floor the same way and you don't need to push it manually.

@JoeJ Yeah, I'll do prebuilt testing.

I managed to swap the whole collision resolution to a more simple version (From the Physics cookbook).

The system is prebuilt to handle several contact points, which Ians doesn't handle out of the box.

As you can see the exact same error seems to be playing out. It's like when calculating the torque or current velocity (including rotation) the system can't seem to recognize from an objective point of view which way the object is already spinning….

@JoeJ I accidentally set up the automated test. Now the object just drops into the world and starts spinning.

@Aressera @joej Alright, here's an update. So I've been wondering for a while if my contact normal is correct. I've always been returning the “correct” contact normal with respect what I see on the screen. The problem may be that the algorithms I use are for multiple objects, where usually the normal points towards object “B”, not “A”, which it actually does in my implementation. So I tried flipping the contactnormal in the Physics Cookbook stuff. Also i decided that the relative vector between the contactpoint and the center of mass might need to be rotated. Just by looking at the videos it seemed that everything rotated in one infinite loop. And I mean, the current velocity from which the correctional impulse stems is related to that vector (the relativecontactposition). So I decided to transform it into worldspace rotation instead of objectspace.

Here's the result of using the velocity correction from the cookbook and the position correction from Ian

Seems like you guys were on to something.. It seems like it seems more reasonable now

So when I don't transform the relative contactposition, but i draw a line from the contactpoint and to the contactpoint + the vector, it looks like this.

So it does seem to point in the direction of the center

When I negate they look like this:

Lastly this is negated and rotated by the world transform matrix

BadProgrammingGuide said:
When I negate they look like this:

This one looks right to me, because the vectors show the arm from the contact to the com. Although if you would draw them from the com, the negation would not be needed. And i assume arm from the com is likely the convention most people use.

The last involving world transform matrix seems clearly wrong.

@JoeJ Currently when I just use the relative contact position based on the COM (both in world space) I end up with endless rotation - doesn't seem like that normal change mattered

I've created the standard test. Drop at 45 degree angle over the x-axis. Friction really makes things weird. Code below.

while (currentiteration < 8) {
   for (int cidx = 0; cidx < contactDataFrame.size(); cidx++) {
    //apply impulse, frictionless
    {
     GVECTOR relativeContactPoint = Math::VectorSubtract(contactDataFrame[cidx].contactpoint, mPhysicsOBJData[0].position);
     GVECTOR closingVelocity = Math::VectorNegate(Math::VectorAdd(mPhysicsOBJData[0].velocity, Math::Vector3Cross(mPhysicsOBJData[0].rotation, relativeContactPoint)));

     float dp = Math::Vector3Dot(closingVelocity, Math::VectorNegate(contactDataFrame[cidx].contactnormal)).m128_f32[0];

     if (dp > 0.0f) {
      continue;
     }

     float numerator = (-(1.0f + restitution)) * dp;
     
     GMATRIX tensor = transformedTensors[0];;  
     GVECTOR torque = Math::Vector3Cross(Math::Vec3MultiplyMatrixRHS(Math::Vector3Cross(relativeContactPoint, Math::VectorNegate(contactDataFrame[cidx].contactnormal)), tensor), relativeContactPoint);

     float denominator = mPhysicsOBJData[0].inverseMass + Math::Vector3Dot(Math::VectorNegate(contactDataFrame[cidx].contactnormal), torque).m128_f32[0];
     if (denominator == 0.0f) {
      PRINT_N("ERROR - denominator J == 0");
      continue;
     }

     float j = numerator / denominator;
     j /= contactDataFrame[cidx].contactcount;

     GVECTOR impulse = Math::VecMultiplyScalar(Math::VectorNegate(contactDataFrame[cidx].contactnormal), j);
   
     mPhysicsOBJData[0].velocity = Math::VectorSubtract(mPhysicsOBJData[0].velocity, Math::VecMultiplyScalar(impulse, mPhysicsOBJData[0].inverseMass));
     mPhysicsOBJData[0].rotation = Math::VectorSubtract(mPhysicsOBJData[0].rotation, Math::Vec3MultiplyMatrixRHS(Math::Vector3Cross(relativeContactPoint, impulse), tensor));

     GVECTOR tangentNormal = Math::VectorSubtract(closingVelocity, Math::VecMultiplyScalar(Math::VectorNegate(contactDataFrame[cidx].contactnormal), Math::Vector3Dot(closingVelocity, Math::VectorNegate(contactDataFrame[cidx].contactnormal)).m128_f32[0]));

     if (sqrt(Math::Vector3Dot(tangentNormal, tangentNormal).m128_f32[0]) == 0.0f) {
      PRINT_N("ERROR - NO TAGENT VECTOR");
      continue;
     }
     tangentNormal = Math::Vec3Normalize(tangentNormal);

     numerator = -Math::Vector3Dot(closingVelocity, tangentNormal).m128_f32[0];
     denominator = mPhysicsOBJData[0].inverseMass + Math::Vector3Dot(tangentNormal, Math::Vector3Cross(Math::Vec3MultiplyMatrixRHS(Math::Vector3Cross(relativeContactPoint, tangentNormal), tensor), relativeContactPoint)).m128_f32[0];
     if (denominator == 0.0f) {
      PRINT_N("ERROR - DENOMINATOR TANGET == 0");
      continue;
     }
     float jt = numerator / denominator;
     jt /= contactDataFrame[cidx].contactcount;

     if (jt == 0.0f) {
      PRINT_N("ERROR - NO JT IMPULSE");
      continue;
     }

     float friction = 3.75f;
     if (jt > j * friction) {
      jt = j * friction;

     }
     else if (jt < -j * friction) {
      jt = -j * friction;
     }

     GVECTOR tagentImpulse = Math::VecMultiplyScalar(tangentNormal, jt);
     mPhysicsOBJData[0].velocity = Math::VectorSubtract(mPhysicsOBJData[0].velocity, tagentImpulse);
     mPhysicsOBJData[0].rotation = Math::VectorSubtract(mPhysicsOBJData[0].rotation, Math::Vec3MultiplyMatrixRHS(Math::Vector3Cross(relativeContactPoint, tagentImpulse), tensor));
    }
   }
   currentiteration++;
}