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

Started by
24 comments, last by BadProgrammingGuide 5 months, 2 weeks ago

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.

Advertisement

@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++;
}

This topic is closed to new replies.

Advertisement