Sign in to follow this  

Contact solver, physics goes crazy :D

This topic is 3587 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

I try to fastly implemet some sort of ball physics for 2d case. And because I'm not so Enstein in physics, I don't know how exactly I suppouse to solve that when I also have angular velocity... Maybe someone who has write that know where I have the bug? Code listed below... [Edited by - Lolmen on February 14, 2008 9:59:43 AM]

Share this post


Link to post
Share on other sites
it's often a sign problem. Try without the last part of the collision response first and just the overlap test, then with the dynamic response, without rotational components, then with rotational components.

Share this post


Link to post
Share on other sites
The separation works just fine, same for rotational and not.
Dynamic test without rotational components it's same, and balls doesn't doing something like staple or glue (usually sign problem), they just explodes from each other like have tons of energy. (impulses problem)

Because I'm not physician and try it out by slides from Erin Catto, I'm don't usderstand much and not see how else I can avoid problem.

Share this post


Link to post
Share on other sites
Using box2d doesn't make me independent and it's not exactly what I'm need.
btw it takes time to learn all stuff to manage objects or writing wrappers.
It's also I think not too smart to use library just for 10 balls physics :D

Share this post


Link to post
Share on other sites
So finally I rewrite it to make more flexibly and clear, but still dunno what potencially can be wrong here...



void ResolveCollision( CCircle *A, CCircle *B, contact_t info )
{
// compute contact in local space
Vector ra = info.position - A->GetOrigin(); // circles have one contact point
Vector rb = info.position - B->GetOrigin();

// compute velocities
Vector va = A->GetVelocity() + ra.Perp() * A->m_flAngularVelocity;
Vector vb = B->GetVelocity() + rb.Perp() * B->m_flAngularVelocity;

// relative velocity at contact
Vector dv = vb - va;

Vector N = info.normal; // collision normal

float vn = dv.Dot(N);

float m1 = 1.0f / A->GetMass();
float m2 = 1.0f / B->GetMass();
float m = (m1 + m2);

// separation
// overlaps? I'm using sweep vs static test
if( vn > 0.0f || info.entire_range ) // both boolean test mark overlapping
{
Vector D = N * info.separation; // use penetration depth
D *= 0.5f; // penetration relaxing? :D
Vector D0, D1;
// should I count rotations and velocities here?
if (m1 > 0.0f)
{
D0 = D * (m1 / m);
A->SetOrigin( A->GetOrigin() + D0 );
}
if (m2 > 0.0f)
{
D1 = D * (m2 / m);
B->SetOrigin( B->GetOrigin() - D1 );
}
}

float friction_s = A->GetCOF() * B->GetCOF();
float elasticity = A->GetCOR() * B->GetCOR();

Vector J = ( ( elasticity * vn ) * N ) * friction_s;

A->m_vVelocity += J * (m1/m);
B->m_vVelocity -= J * (m2/m);

A->m_flAngularVelocity += Cross(ra, J) / A->m_flInertia;
B->m_flAngularVelocity -= Cross(rb, J) / B->m_flInertia;

}








It's works and I think it works fun, not sure how much it's correct, but I'm doesn't seen anything what annoys me, but something still unnatural.

[Edited by - Lolmen on February 14, 2008 10:11:53 AM]

Share this post


Link to post
Share on other sites
One guy told me that the problem is in computing the "Normal Impulse"
and probably because of absence "Tangential Impulse".
But, how can I calculate and apply them right?

Oliver, maybe you can help me?

Thanks for every one in advance.

Share this post


Link to post
Share on other sites
there were too many errors, so should be easier for me to directly fix the code...



void ResolveIntersection( CCircle *A, CCircle *B, contact_t info )
{
// we have penetration
if(info.separation > 0.0f)
{
float ma = 1.0f / A->GetMass();
float mb = 1.0f / B->GetMass();
float m = (ma + mb);

// EDIT: for consistency with relative velocity,
// normal should be pointing from A centre towards B centre.
// so that in both cases, A is the reference body, B is the acting body.
// you should have to swap the normal around.
Vector PA = A->GetOrigin();
Vector PB = B->GetOrigin()
assert(info.normal.dot(PB - PA) > 0.0f); // should always be true

float relaxation = 0.5f;
Vector D = (info.normal * info.separation) * relaxation; // spearation vector

// separation vectors
Vector DA = D * (ma / m);
Vector DB = D * (mb / m);

// separate
A->SetOrigin( PA - DA );
B->SetOrigin( PB + DB );
}
}

void ResolveCollision(CCircle *A, CCircle *B, contact_t info)
{
// intersection relaxation
ResolveIntersection(A, B, info);

// compute contact in local space
// EDIT: contact points should be one for each body, especially when you have intersection.
Vector ra = info.positionA - A->GetOrigin(); // circles have one contact point
Vector rb = info.positionB - B->GetOrigin();

// compute velocities
Vector va = A->GetVelocity() + ra.Perp() * A->m_flAngularVelocity;
Vector vb = B->GetVelocity() + rb.Perp() * B->m_flAngularVelocity;

Vector v = vb - va; // relative velocity at contact
Vector n = info.normal; // collision normal

// impact speed
float vn = v.Dot(n);

// objects are seperating. do not compute a collision and friction impulse.
if (vn > 0.0f)
return;

float ma = 1.0f / A->GetMass();
float mb = 1.0f / B->GetMass();
float m = (ma + mb);
float ia = A->m_flInertia;
float ib = B->m_flInertia;

// EDIT : more accurate representation fo the combined friction, although not accurate enough.
float elasticity = (A->GetCOR() + B->GetCOR()) * 0.5f;

// calculate collision impulse
float ta = ra.cross(info.normal); // 'a' term
float tb = rb.cross(info.normal); // 'b' term
float c = vn / (m + (ta*ta)*ia + (tb*tb)*ib); // final impulse

// collision impulse
Vector C = (c * elasticity) * n;

// friction impulse (0 by default).
// EDIT : compute friction impulse separately
Vector F(0, 0);

if(1)
{
// average friction
float friction = (A->GetCOF() + B->GetCOF()) * 0.3f;

// direction for friction.
n = info.normal.Perp();
if(vab.dot(n) > 0.0f) n *= -1.0f; // flick around if wrong direction.

float vn = v.Dot(n);
float ta = ra.cross(n);
float tb = rb.cross(n);
float f = vn / (m + (ta*ta)*ia + (tb*tb)*ib);

// friction impulse
F = (f * friction) * n;
}

// combined impulse.
Vector J = C + F;

// EDIT : error with the unit consistency
A->m_vVelocity -= J * ma;
B->m_vVelocity += J * mb;

// EDIT : error with the unit consistency
A->m_flAngularVelocity -= Cross(ra, J) * ia;
B->m_flAngularVelocity += Cross(rb, J) * ib;
}




try that. The signs should be consistent now, as well as the units and values. For sign checks, paper and pen usually best.

I've marked with EDIT tags important things.

Share this post


Link to post
Share on other sites
Well, thanks a lot!
I will modify couple of things, because of normal is pointing in other side, and second contact point, well now I understand why even for spheres they sould be two... And thats work defenitely better except that Friction impulse, not sure how
much friction impulse depends on two contact points except one, and I fix that contacts and normal problems slightly later (just can't think much), but just for sure, make sure that it's as good as rest part, because if you places computations of "friction" in within brackets then I suspect that is have some reason for that. :)

Also very big thanks that you undestand what I'm talking about, and even don't scary of reading my kinda code :D

Now after couple of really crazy nights fighting with these poor ball physics simulation I finally can go take some sleep.

Share this post


Link to post
Share on other sites
I put the friction in 'brackets' as it's not a very accurate model, and you can already get good results without using the friction at all, but the ball wont rotate on impact. It's the friction that adds spin to the ball motion.

In fact, it should probably be in another function like ResolveFriction(A, B, collision_info).

the code would be a bit clearer, and you can use alternative friction models.

The friction model I use is basically like a collison impulse, but set at 90 degrees from the collision normal, and going opposite the relative velocity (to reduce the energy on impact). I could even it's completely 'made up', I'm not sure if it bare any reality at all, but it works (or should if I got my signs right) to some degree.

Share this post


Link to post
Share on other sites
So the results I achieve, is really unstability, when I put it all together, I mean improvements in Intersection test, and to make it work more stably, I modify maths like this, and at this moment it's work stable. But I have no clue why it suppose to be like that one, I modify inertia and signs and also couple of things.



void ResolveIntersection( CCircle *A, CCircle *B, contact_t info )
{
// we have penetration
if( info.entire_range )
{
float ma = 1.0f / A->GetMass();
float mb = 1.0f / B->GetMass();
float m = (ma + mb);

// EDIT: for consistency with relative velocity,
// normal should be pointing from A centre towards B centre.
// so that in both cases, A is the reference body, B is the acting body.
// you should have to swap the normal around.
Vector PA = A->GetOrigin();
Vector PB = B->GetOrigin();
// no more need assertion here!
float relaxation = 0.5f;
Vector D = (info.normal * info.separation) * relaxation; // sepearation vector

// separation vectors
Vector DA = D * (ma / m);
Vector DB = D * (mb / m);

// separate
A->SetOrigin( PA - DA );
B->SetOrigin( PB + DB );
}
}
void ResolveCollision(CCircle *A, CCircle *B, contact_t info)
{
// intersection relaxation
ResolveIntersection(A, B, info);

// compute contact in local space
// EDIT: contact points should be one for each body, especially when you have intersection.
Vector ra = info.positionA - A->GetOrigin();
Vector rb = info.positionB - B->GetOrigin();

// compute velocities
Vector va = A->GetVelocity() + ra.Perp() * A->m_flAngularVelocity;
Vector vb = B->GetVelocity() + rb.Perp() * B->m_flAngularVelocity;

Vector v = vb - va; // relative velocity at contact
Vector n = info.normal; // collision normal

// impact speed
float vn = v.Dot(n);

// objects are separating. do not compute a collision and friction impulse.
if (vn > 0.0f)
return;

float ma = 1.0f / A->GetMass();
float mb = 1.0f / B->GetMass();
float m = (ma + mb);
// inertia maths changes
float ia = 1.0f / A->m_flInertia;
float ib = 1.0f / B->m_flInertia;
float i = (ia + ib);

// calculate collision impulse
float ta = ra.Cross(n); // 'a' term
float tb = rb.Cross(n); // 'b' term
float c = vn / (m + (ta*ta)*ia + (tb*tb)*ib); // final impulse

// EDIT : more accurate representation for the combined friction, although not accurate enough.
float elasticity = (A->GetCOR() + B->GetCOR()) * 0.5f;
// collision impulse
Vector C = (c * elasticity) * n;

// friction impulse (0 by default).
// EDIT : compute friction impulse separately
Vector F(0, 0);

// if(1)
{
// average friction
float friction = (A->GetCOF() + B->GetCOF()) * 0.3f;

// direction for friction.
n = n.Perp();
if(v.Dot(n) > 0.0f) n = -n; // flick around if wrong direction.
// no need to declaration twice
vn = v.Dot(n);
ta = ra.Cross(n);
tb = rb.Cross(n);
float f = vn / (m + (ta*ta)*ia + (tb*tb)*ib); // tangential impulse, I suspect... :D

// friction impulse
F = (f * friction) * n;
}

// combined impulse.
Vector J = C + F;

// EDIT : error with the unit consistency
A->m_vVelocity += J * (ma / m); // mass computing changes
B->m_vVelocity -= J * (mb / m); // also signs changes

// EDIT : error with the unit consistency
A->m_flAngularVelocity += Cross(ra, J) * (ia/i); // same changes in inertia
B->m_flAngularVelocity -= Cross(rb, J) * (ib/i);
}





Olivier - any thoughs why it suppose to be like that?
Also I'm curious about friction impulse, you said that it's not so accuratiest, is there exist any information about more accuratiest version? :D

Share this post


Link to post
Share on other sites
ia and ib are inverse inertias so that's correct.

the impulse J should also have a negative sign. try




// collision impulse
Vector C = -(1 + elasticity) * (c) * n;

// friction impulse (0 by default).
// EDIT : compute friction impulse separately
Vector F(0, 0);

if(1)
{
// average friction
float friction = (A->GetCOF() + B->GetCOF()) * 0.3f;

// direction for friction.
n = info.normal.Perp();
if(vab.dot(n) < 0.0f) n *= -1.0f; // flick around if wrong direction.

float vn = v.Dot(n);
float ta = ra.cross(n);
float tb = rb.cross(n);
float f = vn / (m + (ta*ta)*ia + (tb*tb)*ib);

// friction impulse
F = (-friction) * (f) * n;
}

// combined impulse.
Vector J = C + F;

// EDIT : error with the unit consistency
A->m_vVelocity -= J * ma;
B->m_vVelocity += J * mb;

// EDIT : error with the unit consistency
A->m_flAngularVelocity -= Cross(ra, J) * ia;
B->m_flAngularVelocity += Cross(rb, J) * ib;
}


Share this post


Link to post
Share on other sites
No that's ruin it, damn.
I even dunno why, I check all the contacts normals and stuff in intersection test, it's fine. Maybe it's because I forced to use inverted coordinate system?

I mean the Y axis is pointing to ground.... like Vector(0,0) is the top left corner and the (800,600) is the bottom right corner.

Share this post


Link to post
Share on other sites
Just write it down on a piece of paper. The direction of normals, impulses and all should make sense then. It's a problem with signs because the equations should be correct.

Share this post


Link to post
Share on other sites

This topic is 3587 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

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