Written by David Moore
on March 07, 2018

Phys 141 - Sample Quadtree Code

Here’s some very functional Quadtree code. I don’t really want to call it “minimal”, because it includes some bonuses that make the code nicer. Files included are:

Vector2.h. This allows me to write position+=velocity*dt; where position and velocity both have type Vector2.
ImageUtil.h and ImageUtil.cpp. These allow me to create and save bitmaps from code.
Physics.h and Physics.cpp. These two files contain all of the quadtree construction and the entirety of the Barnes-Hut algorithm! Physics.cpp is heavily commented and is still only 108 lines.
main.cpp. This puts everything together.

Download Quadtree.zip Here. You will have to compile it with g++. (g++ -o quadtree <all .cpp files> should work)

Here’s a 10 million particle run, taking about 12 hours. (60 seconds per force calculation, 8 seconds to build the quadtree, 2 seconds to delete it, negligible time to timestep and render the image).

Here is the commented code which builds the quadtree (addParticle) and traverses it (calculateForce)

 
//The magic function that builds the quadtree 
void QuadNode::addPointMass(PointMass *arg)
{
	if(val==nullptr) { //if this node does not store a single particle
		if(subnodes[0][0]==nullptr)//no current mass & no subnodes
		{
			val=arg; //store the mass at this node
		} else { //else we have subnodes: we should pass the mass along so that the subnodes can deal with the particle.
			//update COM and total M.
			masspos+=arg->position*arg->mass;
			mass+=arg->mass;
			Vector2 mid=(tl+br)*.5;
			int x=arg->position.x>mid.x; //  (int)true==1, (int)false==0
			int y=arg->position.y<mid.y;
			subnodes[y][x]->addPointMass(arg);
		}
	}
	else {
		//Else this node stores a single particle and no subnodes. We must create four subnodes and decide which subnode to add 
		//each particle (arg and val) to.

		//If the particles are right next to each other, just give up. In a scientific quality algorithm this might be bad.
		if((fabs(val->position.x-arg->position.x )<0.0001) && (fabs(val->position.y - arg->position.y)<0.0001))
			return;
		PointMass *valtemp=val;
		val=nullptr;
		createSubNodes();
		masspos=Vector2(0,0);
		mass=0;
		addPointMass(valtemp);
		addPointMass(arg);
	}
}

//Add to the PointMass's acceleration vector.
//This is the force calculation with the algorithm of Barnes and Hut.
void QuadNode::calculateForces(PointMass *arg, const double& G) const
{
	if (subnodes[0][0] != nullptr)
	{
		// we get here iff this node has subnodes
		Vector2 pt = masspos / mass;  //center of mass position
		Vector2 diff = pt - arg->position; //vector to the center of mass
		double d = diff.x*diff.x + diff.y*diff.y;
		register double diffw = (tl.x - br.x);

		if (diffw*diffw >= d * .7) //if the size of the node divided by the distance is large, then recurse
		{
			for (int n = 0; n<2; n++)
				for (int m = 0; m<2; m++)
					subnodes[m][n]->calculateForces(arg, G);
		}
		else //else just calculate the force to the center of mass
		{
			Vector2 diff = pt - arg->position;
			double d = diff.x*diff.x + diff.y*diff.y;
			double force = G * arg->mass*mass / d;
			d = sqrt(d);
			if (d>0.000001)
				arg->gravforce += diff * force / d;
		}
	}
	else if (val != nullptr)
	{
		//we get here iff this node does not have subnodes and it does have a particle.
		Vector2 diff = val->position - arg->position;
		double d = diff.x*diff.x + diff.y*diff.y;
		double force = G * arg->mass*val->mass / d;
		d = sqrt(d);
		if (d>0.5)
			arg->gravforce += diff * force / d;
	}
}

← → Top