Sign in to follow this  
Sirisian

Impulse 2D physics response

Recommended Posts

I'm trying to convert the information in this: http://www.cs.cmu.edu/~baraff/sigcourse/notesd2.pdf into code. Basically page 18 has been confusing me. It's about the equation for the impulse J on the top of page 17. I can't seem to figure out how to convert the idea into code. "n * ((Iinv * (ra x n)) x ra)" if you know that for x and y that are vectors that * is a dot product and x is the cross product. Iinv is just an inverted matrix, ra and n are vectors. I can't seem to understand how to turn this into code since the vectors are 2D. "n * (Iinv * Cross(ra, n) x ra" is what I've come up with since I can't figure out what to do next. The cross of 2 2D vectors is just a scalar right? I've been told it's the z component. So vector * (matrix * scalar) x vector. I'm completely lost as to how to interpret this for 2D. I also can't seem to figure out what: c->a->P += force; c->b->P -= force; c->a->L += ra x force; c->b->L -= rb x force; it's never explained what P and L are. If anyone can shed some light on some of these things it would be helpful. Also with the inertia tensor I was told for 2D you'd have the (1,1) to be the max Y extent squared times mass and then (2,2) to be the max X extent squared time mass. Is this right?
            Vector2 cb1ContactPointVel = new Vector2((float)(cb1.Velocity.X - (cp.X - cb1.Position.X) * cb1.RotationRate * Math.Sin(cb1.RotationRate * t_ + cb1.Rotation) - (cp.Y - cb1.Position.Y) * cb1.RotationRate * Math.Cos(cb1.RotationRate * t_ + cb1.Rotation)),
                                                     (float)(cb1.Velocity.Y - (cp.Y - cb1.Position.Y) * cb1.RotationRate * Math.Sin(cb1.RotationRate * t_ + cb1.Rotation) + (cp.X - cb1.Position.X) * cb1.RotationRate * Math.Cos(cb1.RotationRate * t_ + cb1.Rotation)));
            Vector2 cb2ContactPointVel = new Vector2((float)(cb2.Velocity.X - (cp.X - cb2.Position.X) * cb2.RotationRate * Math.Sin(cb2.RotationRate * t_ + cb2.Rotation) - (cp.Y - cb2.Position.Y) * cb2.RotationRate * Math.Cos(cb2.RotationRate * t_ + cb2.Rotation)),
                                                     (float)(cb2.Velocity.Y - (cp.Y - cb2.Position.Y) * cb2.RotationRate * Math.Sin(cb2.RotationRate * t_ + cb2.Rotation) + (cp.X - cb2.Position.X) * cb2.RotationRate * Math.Cos(cb2.RotationRate * t_ + cb2.Rotation)));


            Vector2 vRel = Vector2.Dot(normal, cb1ContactPointVel - cb2ContactPointVel);
            const float THRESHOLD = 0.001f;
            if (vrel > THRESHOLD) /* moving away */
                return FALSE;
            if (vrel > -THRESHOLD) /* resting contact */
                return FALSE;
            else /* vrel < -THRESHOLD */
                return TRUE;

            Vector2 r1 = cp - (cb1.Position + cb1.Velocity * t_);
            Vector2 r2 = cp - (cb2.Position + cb2.Velocity * t_);
            
            Matrix inertiaTensor1 = new Matrix();
            inertiaTensor1.M11 = cb1.GetMaxYextentSquared() * cb1.Mass; //distanceSquared1Y
            inertiaTensor1.M22 = cb1.GetMaxXextentSquared() * cb1.Mass; //distanceSquared1X
            inertiaTensor1.M33 = 1;
            inertiaTensor1.Invert();

            Matrix inertiaTensor2 = new Matrix();
            inertiaTensor2.M11 = cb2.GetMaxYextentSquared() * cb2.Mass; //distanceSquared2Y
            inertiaTensor2.M22 = cb2.GetMaxXextentSquared() * cb2.Mass; //distanceSquared2X
            inertiaTensor2.M33 = 1;
            inertiaTensor2.Invert();

            float term1 = 1 / cb1.Mass;
            float term2 = 1 / cb2.Mass;
            float term3 = n * ((c->a->Iinv * (ra x normal)) x ra);
            float term4 = n * ((c->b->Iinv * (rb x normal)) x rb);

            float J = (-(1 + e) / vRel) / (term1 + term2 + term3 + term4);
            Vector2 force = J * normal;
            
            //Not sure what P and L are
            c->a->P += force;
            c->b->P -= force;
            c->a->L += Vector2.Ccw(r1, force);
            c->b->L -= Vector2.Ccw(r1, force);
            /* recompute auxiliary variables */
            cb1.Velocity = c->a->P / cb1.Mass;
            cb2.Velocity = c->b->P / cb2.Mass;

            c->a->omega = c->a->Iinv * c->a->L;
            c->b->omega = c->b->Iinv * c->b->L;

Share this post


Link to post
Share on other sites
Quote:
I can't seem to understand how to turn this into code since the vectors are 2D.


Use 3D vectors, then.

Quote:
it's never explained what P and L are.


P is the linear momentum
L is the angular momentum

P = m * V
L = r x P

The cross product will always be perpendicular to the two vectors you're multiplying. If you're using two dimensions, then that vector will not be in your defined space. You can still determine the magnitude by standard rules, however.

|AxB| = AxBy - AyBx

The problem with this, is that you can't really use this number in a convenient way to update your quantities. L will stick out of your plane, and not necessarily in the same direction as your axis of rotation. You'll lose a ton of information.

The inertia tensor, assuming all parts have a 0 z component is:
[my2 -mxy 0]
[-mxy mx2 0]
[0 0 m(x2+y2) 0]

Noticeably, this is still in 3D.

Share this post


Link to post
Share on other sites
Well, I suppose you don't have to do a full conversion, as long as you remember that the value of L represents the z-component of that vector. The inertia tensor, in this case, can really be reduced to the a33 value, because any axis of rotation is going to also have only a z-component.

Share this post


Link to post
Share on other sites
P looks like the net force on the body and L the net torque.

But I'll just throw my hat in the ring... It's based on Chris Hecker's articles, just another spin.

http://members.gamedev.net/oliii/polygon.rar



// apply collision impusles at contacts
for(int i = 0; i < m_contacts.m_count; i ++)
{
Vector pa = m_contacts.m_contact[i].m_position[0];
Vector pb = m_contacts.m_contact[i].m_position[1];
Vector ra = pa - a->m_position;
Vector rb = pb - b->m_position;
Vector va = a->m_velocity + ra.perp() * a->m_angvelocity;
Vector vb = b->m_velocity + rb.perp() * b->m_angvelocity;

float vn = ((va - vb) * m_ncoll);
if(vn > 0.0f) continue;

float ta = (ra ^ m_ncoll) * (ra ^ m_ncoll) * a->m_invinertia;
float tb = (rb ^ m_ncoll) * (rb ^ m_ncoll) * b->m_invinertia;
float m = a->m_invmass + b->m_invmass;
float denom = m + ta + tb;
float j = -(1.0f + CoR) * vn / denom;

Vector impulse = m_ncoll * j;

a->m_velocity += impulse * a->m_invmass;
b->m_velocity -= impulse * b->m_invmass;

a->m_angvelocity += (ra ^ impulse) * a->m_invinertia;
b->m_angvelocity -= (rb ^ impulse) * b->m_invinertia;
}




Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this