# 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.

## 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 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 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 on other sites
Why not use Erin Catto's excellent Box2d library?

http://www.box2d.org/

Edit: BTW, a "physician" is a doctor of medicine, "physicist" is the word you want.

##### 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 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 on other sites
A->m_vVelocity += J * (m1 / m2);B->m_vVelocity -= J * (m2 / m);

supposed to be m1/m?

##### Share on other sites
Yep, but that not fix crazy spinning stuff...

##### 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 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 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 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 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 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 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 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 on other sites
inline void LimitAngle(float& fAngle){	if ((fAngle < 0) || (fAngle > 6.28))	{		fAngle = fmod(fAngle, 6.28);	}}

##### 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.

## Create an account

Register a new account

• ### Forum Statistics

• Total Topics
628711
• Total Posts
2984342

• 23
• 11
• 10
• 13
• 14