need help with 2d circle collision physics :(

Author

Message

Neuro Fuzzy

17

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 6th Aug 2010 08:38

Link

I've been struggling with this for a while. I dunno if any of you have noticed my n-body sim code anywhere, but I've discovered the collision method I was using is more bullcrap than I thought it was (I knew it was innacurate but not that innacurate), so I've been struggling to find a decent collision algorithm that will preserve momentum and rotational momentum over many steps. I have the physics pretty good, so that the total momentum of the system only changes in the range of 10^-5 or 10^-6 per step... but the rotational velocity of bodies keeps increasing

Meaning, if you have, say, 20 bodies colliding, at first they behave like you would expect, and form an approximate circle... but then - seemingly out of nowhere - the 20 bodies start spinning in one direction, faster and faster until the centripetal force rips the bodies apart, all from nowhere.
Here's the source code:

#include "math.h"
#include <vector>
#include "DarkGDK.h"
#include "D3DFunc.h"
using namespace std;
class circPhy;

class vector2
{
    public:
        float x;
        float y;
        vector2();
		vector2(const float& a,const float& b);
		const vector2 operator+	(const vector2& b) const;	//vector addition
		vector2& operator+=		(const vector2& b);
		const vector2 operator-	(const vector2& b) const;	//vector subtraction
		vector2& operator-=		(const vector2& b);
		const vector2 operator*	(const float& b) const;	//vector scalar multiplication
		vector2& operator*=		(const float& b);
		const vector2 operator/	(const float& b) const;	//vector scalar division
		vector2& operator/=		(const float& b);
		const float operator*	(const vector2& b) const;	//dot product
		void normalize();
		void normalize(const float& dist);
		vector2 normalized();
		vector2 normalized(const float& dist);
		float length();
};
static const vector2 operator* (const float& b, const vector2& a);

class verletPoint
{   
    public:   
        vector2 pos;
        vector2 vel;
        vector2 acc;
		vector2 lAcc;
        verletPoint();
        verletPoint(const vector2& position,const vector2& velocity,const vector2& acceleration);
        verletPoint(const vector2& position,const vector2& velocity);
        verletPoint(const vector2& position);
        void integrate(const float& timestep);
};
class circPhy
{
    public:
		static float G;
		static float COR;
        float M;
        float R;
		bool updateable;
        verletPoint center;
        circPhy(){}
        circPhy(const float& center_x,const float& center_y,const float& radius,const float& mass, const bool& moveable);
        void update(const float& timestep);
        void coll(circPhy* colObj);
		void addGrav(const vector2& pos, const float& mass);
		void zeroAcc();
};

void DarkGDK(void)
{
	dbSetDisplayMode(640,480,32);
	int number_of_bodies = 1;
	int lastclick=0;
	d3dInit();
	dbSyncOn();
	dbRandomize(1);

vector<circPhy> r;
	int moveobj=0;
	for(int n=0;n<number_of_bodies;n++)
		r.push_back(circPhy(dbRnd(dbScreenWidth()),dbRnd(dbScreenHeight()),20, 10,true));
	while ( LoopGDK ( ) )
	{
		dbCLS();
		dbInk(dbRGB (0, 0, 255),dbRGB (0, 0, 255));
		vector2 sum=vector2(0.0,0.0);
		if(!dbSpaceKey())
		{
			for(int n=0;n<number_of_bodies;n++)
			{
				for(int i=0;i<number_of_bodies;i++)
					if(i!=n)
						r[n].addGrav(r[i].center.pos,r[i].M);
			}
			for(int n=0;n<number_of_bodies;n++)
			{
				for(int i=0;i<number_of_bodies;i++)
					if(i!=n)
						r[n].coll(&r[i]);
			}			
			for(int n=0;n<number_of_bodies;n++)
			{
				sum+=r[n].center.vel*r[n].M;
				r[n].update(2.0f);
			}
		}
		if(dbMouseClick() && lastclick!=dbMouseClick())
		{
			if(dbMouseClick()==1)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(20)+10,dbRnd(10)+2,true));
				number_of_bodies++;
			}
			if(dbMouseClick()==2)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(5)+5,1,false));
				number_of_bodies++;
			}
			if(dbMouseClick()==4)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(100)+5,100,false));
				number_of_bodies++;
				moveobj=r.size()-1;
			}
		}
		if(dbMouseClick()==4)
		{
			r[moveobj].center.pos=vector2(dbMouseX(),dbMouseY());
		}
		if(lastclick==4 && dbMouseClick()!=4)
		{
			r.erase(r.begin()+moveobj);
			number_of_bodies--;
		}
		lastclick=dbMouseClick();
		for(int n=0;n<number_of_bodies;n++)
			d3dCircle(r[n].center.pos.x,r[n].center.pos.y,r[n].R,0);
		dbText(0,0,"number of bodies:");
		dbText(0,16,dbStr(number_of_bodies));
		dbText(0,32,"framerate (capped @ 60):");
		dbText(0,48,dbStr(dbScreenFPS()));
		dbText(0,64,dbStr((float)sum.x));
		dbText(0,80,dbStr((float)sum.y));
		dbSync();
	}
	return;
}
//class vector2
	vector2::vector2()
	{
		x = y = 0.0f;
	}
	vector2::vector2(const float& a,const float& b)
	{
		x=a;
		y=b;
	}
	const vector2 vector2::operator+ (const vector2& b) const
	{
		vector2 temp;
		temp.x=x+b.x;
		temp.y=y+b.y;
		return temp;
	}
	vector2& vector2::operator+= (const vector2& b)
	{
		if(this!=&b)
		{
			x=x+b.x;
			y=y+b.y;
		}
		return *this;
	}
	const vector2 vector2::operator-(const vector2& b) const
	{
		vector2 temp;
		temp.x=x-b.x;
		temp.y=y-b.y;
		return temp;
	}
	vector2& vector2::operator-=(const vector2& b)
	{
		if(this!=&b)
		{
			x=x-b.x;
			y=y-b.y;
		}
		return *this;
	}
	const float vector2::operator* (const vector2& b) const
	{
		return (x*b.x+y*b.y);
	}
	const vector2 vector2::operator* (const float& b) const
	{
		vector2 temp;
		temp.x=x*b;
		temp.y=y*b;
		return temp;
	}
	vector2& vector2::operator*= (const float& b)
	{
		x=x*b;
		y=y*b;
		return *this;
	}
	const vector2 vector2::operator/ (const float& b) const
	{
		vector2 temp=vector2(0.0f,0.0f);
		if(b>0)
		{
			temp.x=x/b;
			temp.y=y/b;
		}
		return temp;
	}
	vector2& vector2::operator/= (const float& b)
	{
		if(b>0)
		{
			x=x/b;
			y=y/b;
		}
		return *this;
	}
	void vector2::normalize()
	{
		float dist = sqrt(x*x+y*y);
		x/=dist;
		y/=dist;
	}
	void vector2::normalize(const float& dist)
	{
		x/=dist;
		y/=dist;
	}
	vector2 vector2::normalized()
	{
		float dist = sqrt(x*x+y*y);
		return vector2(x/dist,y/dist);
	}
	vector2 vector2::normalized(const float& dist)
	{
		return vector2(x/dist,y/dist);
	}
	float vector2::length()
	{
		float n=x*x+y*y;
		return sqrt(n>0?n:0);
	}
static const vector2 operator* (const float& b,const vector2& a)
{
	return a*b;
}

verletPoint::verletPoint()
{
	vector2 pos (0.0f,0.0f);
	vector2 vel (0.0f,0.0f);
	vector2 acc (0.0f,0.0f);
	vector2 lAcc(0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position, const vector2 & velocity, const vector2 & acceleration)
{
	pos=position;
	vel=velocity;
	acc=acceleration;
	vector2 lAcc(0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position,const vector2 & velocity)
{
	pos=position;
	vel=velocity;
	vector2 acc (0.0f,0.0f);
	vector2 lAcc(0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position)
{
	pos=position;
	vector2 vel (0.0f,0.0f);
	vector2 acc (0.0f,0.0f);
	vector2 lAcc(0.0f,0.0f);
}
void verletPoint::integrate(const float& timestep)
{
	//	vel+=acc*timestep;
	//pos+=vel*timestep;

pos+=vel*timestep+lAcc*timestep*timestep/2;	//P(t+dt)~=P(t)+V(t)*dt+A(t)*((dt)^2)/2
	vel+=(lAcc+acc)*timestep/2;					//V(t+dt)~=V(t)+(A(t)+A(t+dt))*dt/2
	lAcc=acc;
}

float circPhy::G=2.0f;
	float circPhy::COR=0.0f;
circPhy::circPhy(const float& center_x,const float& center_y,const float& radius,const float& mass, const bool& moveable)
{
	M=mass;
	R=radius;
	updateable=moveable;
	center=verletPoint(vector2(center_x,center_y));
}
void circPhy::update(const float& timestep)
{
	if(updateable)
		center.integrate(timestep);
	zeroAcc();
}
void circPhy::coll(circPhy* colObj)
{
	float d=0.0f;
	vector2 x=center.pos-colObj->center.pos;
	d=x.length();
	x.normalize(d);
	float r1=R+colObj->R;
	if(r1>d)
	{
		float v1x=center.vel*x;
		float v2x=colObj->center.vel*x;
		center.vel=center.vel-(v1x-((M*v1x+colObj->M*v2x)/(M+colObj->M)))*x;
		colObj->center.vel=colObj->center.vel-(v2x-((M*v1x+colObj->M*v2x)/(M+colObj->M)))*x;
		center.pos+=(r1-d)/2*x;
		colObj->center.pos-=(r1-d)/2*x;
	}
}
void circPhy::addGrav(const vector2& pos, const float& mass)
{
	vector2 d=pos-center.pos;
	float n=d.length();
	center.acc+=d*(G*mass/(n*n*n));
}
void circPhy::zeroAcc()
{
	center.acc=vector2(0.0f,0.0f);
}

+ Code Snippet

#include "math.h"
#include <vector>
#include "DarkGDK.h"
#include "D3DFunc.h"
using namespace std;
class circPhy;

class vector2
{
    public:
        float x;
        float y;
        vector2();
		vector2(const float& a,const float& b);
		const vector2 operator+	(const vector2& b) const;	//vector addition
		vector2& operator+=		(const vector2& b);
		const vector2 operator-	(const vector2& b) const;	//vector subtraction
		vector2& operator-=		(const vector2& b);
		const vector2 operator*	(const float& b) const;	//vector scalar multiplication
		vector2& operator*=		(const float& b);
		const vector2 operator/	(const float& b) const;	//vector scalar division
		vector2& operator/=		(const float& b);
		const float operator*	(const vector2& b) const;	//dot product
		void normalize();
		void normalize(const float& dist);
		vector2 normalized();
		vector2 normalized(const float& dist);
		float length();
};
static const vector2 operator* (const float& b, const vector2& a);

class verletPoint
{   
    public:   
        vector2 pos;
        vector2 vel;
        vector2 acc;
		vector2 lAcc;
        verletPoint();
        verletPoint(const vector2& position,const vector2& velocity,const vector2& acceleration);
        verletPoint(const vector2& position,const vector2& velocity);
        verletPoint(const vector2& position);
        void integrate(const float& timestep);
};
class circPhy
{
    public:
		static float G;
		static float COR;
        float M;
        float R;
		bool updateable;
        verletPoint center;
        circPhy(){}
        circPhy(const float& center_x,const float& center_y,const float& radius,const float& mass, const bool& moveable);
        void update(const float& timestep);
        void coll(circPhy* colObj);
		void addGrav(const vector2& pos, const float& mass);
		void zeroAcc();
};



void DarkGDK(void)
{
	dbSetDisplayMode(640,480,32);
	int number_of_bodies = 1;
	int lastclick=0;
	d3dInit();
	dbSyncOn();
	dbRandomize(1);

	vector<circPhy> r;
	int moveobj=0;
	for(int n=0;n<number_of_bodies;n++)
		r.push_back(circPhy(dbRnd(dbScreenWidth()),dbRnd(dbScreenHeight()),20, 10,true));
	while ( LoopGDK ( ) )
	{
		dbCLS();
		dbInk(dbRGB (0, 0, 255),dbRGB (0, 0, 255));
		vector2 sum=vector2(0.0,0.0);
		if(!dbSpaceKey())
		{
			for(int n=0;n<number_of_bodies;n++)
			{
				for(int i=0;i<number_of_bodies;i++)
					if(i!=n)
						r[n].addGrav(r[i].center.pos,r[i].M);
			}
			for(int n=0;n<number_of_bodies;n++)
			{
				for(int i=0;i<number_of_bodies;i++)
					if(i!=n)
						r[n].coll(&r[i]);
			}			
			for(int n=0;n<number_of_bodies;n++)
			{
				sum+=r[n].center.vel*r[n].M;
				r[n].update(2.0f);
			}
		}
		if(dbMouseClick() && lastclick!=dbMouseClick())
		{
			if(dbMouseClick()==1)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(20)+10,dbRnd(10)+2,true));
				number_of_bodies++;
			}
			if(dbMouseClick()==2)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(5)+5,1,false));
				number_of_bodies++;
			}
			if(dbMouseClick()==4)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(100)+5,100,false));
				number_of_bodies++;
				moveobj=r.size()-1;
			}
		}
		if(dbMouseClick()==4)
		{
			r[moveobj].center.pos=vector2(dbMouseX(),dbMouseY());
		}
		if(lastclick==4 && dbMouseClick()!=4)
		{
			r.erase(r.begin()+moveobj);
			number_of_bodies--;
		}
		lastclick=dbMouseClick();
		for(int n=0;n<number_of_bodies;n++)
			d3dCircle(r[n].center.pos.x,r[n].center.pos.y,r[n].R,0);
		dbText(0,0,"number of bodies:");
		dbText(0,16,dbStr(number_of_bodies));
		dbText(0,32,"framerate (capped @ 60):");
		dbText(0,48,dbStr(dbScreenFPS()));
		dbText(0,64,dbStr((float)sum.x));
		dbText(0,80,dbStr((float)sum.y));
		dbSync();
	}
	return;
}
//class vector2
	vector2::vector2()
	{
		x = y = 0.0f;
	}
	vector2::vector2(const float& a,const float& b)
	{
		x=a;
		y=b;
	}
	const vector2 vector2::operator+ (const vector2& b) const
	{
		vector2 temp;
		temp.x=x+b.x;
		temp.y=y+b.y;
		return temp;
	}
	vector2& vector2::operator+= (const vector2& b)
	{
		if(this!=&b)
		{
			x=x+b.x;
			y=y+b.y;
		}
		return *this;
	}
	const vector2 vector2::operator-(const vector2& b) const
	{
		vector2 temp;
		temp.x=x-b.x;
		temp.y=y-b.y;
		return temp;
	}
	vector2& vector2::operator-=(const vector2& b)
	{
		if(this!=&b)
		{
			x=x-b.x;
			y=y-b.y;
		}
		return *this;
	}
	const float vector2::operator* (const vector2& b) const
	{
		return (x*b.x+y*b.y);
	}
	const vector2 vector2::operator* (const float& b) const
	{
		vector2 temp;
		temp.x=x*b;
		temp.y=y*b;
		return temp;
	}
	vector2& vector2::operator*= (const float& b)
	{
		x=x*b;
		y=y*b;
		return *this;
	}
	const vector2 vector2::operator/ (const float& b) const
	{
		vector2 temp=vector2(0.0f,0.0f);
		if(b>0)
		{
			temp.x=x/b;
			temp.y=y/b;
		}
		return temp;
	}
	vector2& vector2::operator/= (const float& b)
	{
		if(b>0)
		{
			x=x/b;
			y=y/b;
		}
		return *this;
	}
	void vector2::normalize()
	{
		float dist = sqrt(x*x+y*y);
		x/=dist;
		y/=dist;
	}
	void vector2::normalize(const float& dist)
	{
		x/=dist;
		y/=dist;
	}
	vector2 vector2::normalized()
	{
		float dist = sqrt(x*x+y*y);
		return vector2(x/dist,y/dist);
	}
	vector2 vector2::normalized(const float& dist)
	{
		return vector2(x/dist,y/dist);
	}
	float vector2::length()
	{
		float n=x*x+y*y;
		return sqrt(n>0?n:0);
	}
static const vector2 operator* (const float& b,const vector2& a)
{
	return a*b;
}


	
verletPoint::verletPoint()
{
	vector2 pos (0.0f,0.0f);
	vector2 vel (0.0f,0.0f);
	vector2 acc (0.0f,0.0f);
	vector2 lAcc(0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position, const vector2 & velocity, const vector2 & acceleration)
{
	pos=position;
	vel=velocity;
	acc=acceleration;
	vector2 lAcc(0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position,const vector2 & velocity)
{
	pos=position;
	vel=velocity;
	vector2 acc (0.0f,0.0f);
	vector2 lAcc(0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position)
{
	pos=position;
	vector2 vel (0.0f,0.0f);
	vector2 acc (0.0f,0.0f);
	vector2 lAcc(0.0f,0.0f);
}
void verletPoint::integrate(const float& timestep)
{
	//	vel+=acc*timestep;
	//pos+=vel*timestep;


	pos+=vel*timestep+lAcc*timestep*timestep/2;	//P(t+dt)~=P(t)+V(t)*dt+A(t)*((dt)^2)/2
	vel+=(lAcc+acc)*timestep/2;					//V(t+dt)~=V(t)+(A(t)+A(t+dt))*dt/2
	lAcc=acc;
}

	float circPhy::G=2.0f;
	float circPhy::COR=0.0f;
circPhy::circPhy(const float& center_x,const float& center_y,const float& radius,const float& mass, const bool& moveable)
{
	M=mass;
	R=radius;
	updateable=moveable;
	center=verletPoint(vector2(center_x,center_y));
}
void circPhy::update(const float& timestep)
{
	if(updateable)
		center.integrate(timestep);
	zeroAcc();
}
void circPhy::coll(circPhy* colObj)
{
	float d=0.0f;
	vector2 x=center.pos-colObj->center.pos;
	d=x.length();
	x.normalize(d);
	float r1=R+colObj->R;
	if(r1>d)
	{
		float v1x=center.vel*x;
		float v2x=colObj->center.vel*x;
		center.vel=center.vel-(v1x-((M*v1x+colObj->M*v2x)/(M+colObj->M)))*x;
		colObj->center.vel=colObj->center.vel-(v2x-((M*v1x+colObj->M*v2x)/(M+colObj->M)))*x;
		center.pos+=(r1-d)/2*x;
		colObj->center.pos-=(r1-d)/2*x;
	}
}
void circPhy::addGrav(const vector2& pos, const float& mass)
{
	vector2 d=pos-center.pos;
	float n=d.length();
	center.acc+=d*(G*mass/(n*n*n));
}
void circPhy::zeroAcc()
{
	center.acc=vector2(0.0f,0.0f);
}

If you run it (d3d_circle is used, but you could replace that with whatever command you want to draw a dot, circle, or box; the effect is the same), you can see that the total momentum (x and y values are the second to last and last values printed on the top left of the screen) stays close to zero. If you have two bodies collide (click to create one), you can see them behave as expected, even with preserving momentum when they collide. But... as you add more and more, you can see that it all spins out of control!!!
What am i doing wrong? First I find the force on every point, then I update collision on every point, then I update velocity and move every point. I've tried switching the order of those to no avail, and the problem doesn't just occur with objects of differing radius/mass/density, though it happens faster with smaller/denser objects.

I could imagine that when two objects aren't moving/accelerating towards each other in exactly opposite directions, they would slide and this would add to rotational velocity... but i've no clue how to solve this.

Is't life, I ask, is't even prudence, to bore thyself and bore thy students?

Back to top

Profile PM Email

Neuro Fuzzy

17

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 7th Aug 2010 02:49

Link

I figure I should post some more info on what I have for the physics so far. The following function

+ Code Snippet

void circPhy::addGrav(const vector2& pos, const float& mass)
{
	vector2 d=pos-center.pos;
	float n=d.length();
	center.acc+=d*(G*mass/(n*n*n));
}

adds gravity based on a vector and mass fed in. Explaining the math of it, this is basically newtons law of universal gravitation:
F=G*m1*m2/d^2
in this case, F=G*m1*m2/|d1|^2

x| being the length of vector x)
however, this is two dimensional, so we want a vector with length F in the direction of d (the difference between the two positions)
direction of d = d/|d|
so F in the direction of d=G*m1*m2/|d1|^2 * d/|d|=d*G*m1*m2/|d1|^3
to find the acceleration of body 1, divide by the mass of m1. This cancels out, so...
acc=d*G*m2/|d1|^3
since this is applied to every object, the law every object has an equal and opposite reaction is obeyed (unless you want it not to be).

this function:

+ Code Snippet

void circPhy::coll(circPhy* colObj)
{
	float d=0.0f;
	vector2 x=center.pos-colObj->center.pos;
	d=x.length();
	x.normalize(d);
	float r1=R+colObj->R;
	if(r1>d)
	{
		float v1x=center.vel*x;
		float v2x=colObj->center.vel*x;
		center.vel=center.vel-(v1x-((M*v1x+colObj->M*v2x)/(M+colObj->M)))*x;
		colObj->center.vel=colObj->center.vel-(v2x-((M*v1x+colObj->M*v2x)/(M+colObj->M)))*x;
		center.pos+=(r1-d)/2*x;
		colObj->center.pos-=(r1-d)/2*x;
	}
}

handles the velocity changes in a collision between two bodies. Basically, you have an inelastic collision that conserves momentum. An inelastic collision is one that leaves both bodies with the same velocity after the collision, and we also say momentum (velocity times mass) is conserved (momentum before collision=momentum after collision), that means:
v1*m1+v2*m2=(m1+m2)*vf (vf is the final velocity)
or
vf=(v1*m1+v2*m2)/(m1+m2)
That's the case in one dimension. Soo... to expand this to two dimensions, we first zero the velocity in the direction of x (again, x is the direction of object 2 to object 1). So
v1=v1-(v1*x)x
(v1*x is the dot product, but the other x is scalar multiplication)
The collision can be seen in one dimension, if we just apply it in the direction of x. However, v1 and v2 are also vectors. We want their magnitude in direction x, so we just take their dot product.
1d delta v1 = (m1*v1*x+m2*v2*x)/(m1+m2)
2d delta v1 = (m1*v1*x+m2*v2*x)/(m1+m2)x
so
v1=v1-(v1*x)x+(m1*v1*x+m2*v2*x)/(m1+m2)x
or
v1=v1-(v1*x-(m1*v1*x+m2*v2*x)/(m1+m2))x
which is all whats going on in this function.

Now... I'm pretty sure this is where I run into trouble. This:
center.pos+=(r1-d)/2*x;
colObj->center.pos-=(r1-d)/2*x;
makes sure the two objects aren't colliding. The sum of the radii is r1, and d is the distance between the two, so if they are colliding (and if this code is being run, they are), then r1-d is the amount they intersect eachother by. Moving each body half the colliding distance apart makes sure they don't intersect (again, multiplying by x ensures its in the direction of x).
If I take out this code, three bodies colliding works fine, but everything intersects everything!

I'm not really sure what I need to do... Is my momentum preserving velocity collision thingymajig going too far? Do I need to do something with normal forces or changing acceleration? I need to have code that makes sure the bodies don't intersect, because otherwise errors could build up, and you could have bodies in almost the exact same spot... which would be bad.

Is't life, I ask, is't even prudence, to bore thyself and bore thy students?

Back to top

Profile PM Email

Dr Tank

15

Years of Service

User Offline

Joined: 1st Apr 2009

Location: Southampton, UK

Posted: 9th Aug 2010 01:51 Edited at: 9th Aug 2010 01:52

Link

I think your problem may be quite simple.

In the stage where you separate overlapping objects, just changing their positions will affect the angular momentum of the system, unless the particles are colliding along a common axis, like the simple 1D problem.

Imagine an ice skater spinning around. She draws her arms in, (like gravity pulling your circles together), and spins around faster. Then your collision system say's "bam your arms are back out again", but she's still spinning around faster. Then her arms come in again and she spins faster and this keeps repeating.

I imagine to fix this, using a centre of mass frame (like
you use for collisions anyway), you could try preserving the angular momentum of this two body system. Maybe the rotational energy. My mechanics is a bit rusty. Basically when you move two circles apart, the component of their velocities about the centre of mass should decrease.

Back to top

Profile PM

Neuro Fuzzy

17

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 9th Aug 2010 09:41 Edited at: 9th Aug 2010 21:17

Link

Just looked into that, and I don't think that's the problem. If I change the order things update in, I can make sure that the gravity pulling the circles together isn't calculated when they are intersecting each other. When that happens, the rotation problem doesn't seem to occur, but another problem is noticeable [edit] (when I also disable collision :X) [/edit]: the objects keep moving towards each other. For two bodies colliding, this happens (each picture is another five seconds or so).

So... they keep moving towards each other when they're colliding, just very slowly. A similar thing happens with three bodies colliding.
The interesting case is with four bodies colliding like this:

the two bodies that aren't colliding move towards each other faster than the ones that are. This means there's a net force acting on the circles... when in reality it should be at a stable equilibrium.

... I think I need to update the normal force of the objects? The thing is I'm not sure of how to do that with objects accelerating and colliding like this.

Is't life, I ask, is't even prudence, to bore thyself and bore thy students?

Back to top

Profile PM Email

Diggsey

18

Years of Service

User Offline

Joined: 24th Apr 2006

Location: On this web page.

Posted: 9th Aug 2010 15:09 Edited at: 9th Aug 2010 15:12

Link

If you can calculate the time that two circle will collide at, you can find the first collision, step all of the movements up to that point, calculate the collision, and then find the next soonest collision, etc. until you have exhausted the time-step.

Unless you can perform the collision calculation at the exact point in time when it happens you will not get an accurate result.

There is also the problem that acceleration isn't constant across the time-step, making the simulation inaccurate to start with.

[b]

Back to top

Profile PM Email Website

Neuro Fuzzy

17

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 9th Aug 2010 21:35 Edited at: 9th Aug 2010 22:54

Link

Quote: "There is also the problem that acceleration isn't constant across the time-step, making the simulation inaccurate to start with.
"

yeah... my dad started asking me numerical analysis questions about this >.>

Quote: "If you can calculate the time that two circle will collide at, you can find the first collision, step all of the movements up to that point, calculate the collision, and then find the next soonest collision, etc. until you have exhausted the time-step."

That sounds... scary... and I don't think I can afford for collision to be that computationally expensive.

I'm pretty sure that my problem is that forces don't cancel out... so bodies are trying to accelerate towards each other, and the velocities are partially canceling out, and the position check is handling the rest (which is bad). After a lot of thought, I think my problems can be solved with the statement: "If two objects will be colliding next time step, their acceleration in direction x will be the same" where x is the vector from one object to another.

But... I think I'll have to write a collision solver for this...

[edit]
yeah, I'm pretty sure I do need to write a collision solver for this. I changed the collision routine so that for each two particles colliding, the acceleration would be the same in the direction they were colliding (with no prediction of whether they would be colliding next time step). This works fine for a couple bodies colliding, but it only ensures the acceleration is the same for the last two bodies checked. The more you iterate this check before updating it, the more accurate your result will be, but that's a lot slower than using a collision solver.

Is't life, I ask, is't even prudence, to bore thyself and bore thy students?

Back to top

Profile PM Email

Sorry your browser is not supported!

Dark GDK / need help with 2d circle collision physics :(