[ODE] Damping

[Joe Ante]

>> Furthermore applying velocity directly is less robust than applying a
>> damping force directly.
> 
> Yeah, I understand that in many cases setting velocities directly is a bit
> taboo, but in this case, all it would mean is it would become a function
> of mass (which it really shouldn't be). I mean, I could do F=MA, but then
> it would just turn around and do V=V+(dT*F/M), and seeing how worried
> people are about the speed impact of clamping a floating point value, I
> think I will try and save the CPU from factor then de-factoring mass. I
> don't think we will be losing any robustness.

Since you will have to use the powf function to calculate the velocity scale
it will be slower. Saying that there are no completely deterministic physics
engines out there is untrue and certainly not a reason to make a
deterministic physics system be undeterministic. In this case its not only
not deterministic but deterministically leads to different simulations
depending on the delta time you use in the simulation.
 
Setting the force is more robust since the constraint system would evaluate
after the damping. If you set the velocity directly the constraint system
has no chance of saying hey this value is not good for me, ill ignore it.

Here is code that robustly adds a drag force. It works correctly also for
rotational drag (The code example given by Jon doesn't work correctly for
non-spherical inertia tensors)
And most importantly it clamps the force value to a maximum value, which is
the force which will make the body stop within one physics frame. This is
extremely important for stability. (Using the velocity appraoch it would
obviously not be necessary but that has the other drawbacks outlined above)

The code is a bit intermixed with our own Vector3 library but I am sure you
can figure it out.


        float invDeltaTime = 1.0F / deltaTime;
        const dReal *vel = dBodyGetLinearVel(body);
        const dReal *w = dBodyGetAngularVel(body);

        // Add linear drag
        Vector3f dragVec, maxDragVec;
        float maxDragScale = -bodyComponent.GetMass () * invDeltaTime;

        // Drag = 0 now means no drag. compensates for mass. Independent of
delta time
        float drag = -bodyComponent.GetDrag () * bodyComponent.GetMass ();
        dragVec = Vector3f (drag * vel[0], drag * vel[1], drag * vel[2]);
        maxDragVec = Vector3f (vel[0] * maxDragScale, vel[1] * maxDragScale,
vel[2] * maxDragScale);
        if (SqrMagnitude (dragVec) > SqrMagnitude (maxDragVec))
            dragVec = maxDragVec;

        dBodyAddForce (body, dragVec.x, dragVec.y, dragVec.z);

        // Rotational drag has to be multiplied by inertia tensor
        drag = -bodyComponent.GetRotationDrag ();
        dVector3 temp;
        dMULTIPLY1_331 (temp, dBodyGetMassFast (body)->I, w);
        
        dragVec = Vector3f (drag * temp[0], drag * temp[1], drag * temp[2]);
        maxDragVec = Vector3f (temp[0] * -invDeltaTime, temp[1] *
-invDeltaTime, temp[2] * -invDeltaTime);
        if (SqrMagnitude (dragVec) > SqrMagnitude (maxDragVec))
            dragVec = maxDragVec;
        
        dBodyAddTorque (body, dragVec.x, dragVec.y, dragVec.z);

Joachim Ante