ball-ball collisions (untested)

6 files changed, 311 insertions, 21 deletions

diff --git a/Makefile b/Makefile
index 1c7b23a..49b04ef 100644
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,7 @@ CC = gcc
 CFLAGS = -std=c99 -Wall -pedantic -Wno-deprecated-declarations
 LDFLAGS = -lGLEW -lGL -lX11 -lGLU -lOpenGL -lOpenCL -lglut -lGLX
 
-SRC = balls.c sysfatal.c geo.c rand.c
+SRC = balls.c sysfatal.c geo.c rand.c partition.c
 OBJ = ${SRC:.c=.o}
 
 balls: ${OBJ}
diff --git a/balls.c b/balls.c
index bf81517..9e2522a 100644
--- a/balls.c
+++ b/balls.c
@@ -129,8 +129,13 @@ initCL(void) {
 	#endif
 
 	/* Get GPU device. */
-	if (clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL) < 0)
-		sysfatal("No GPUs available.\n");
+	err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);
+	if (err < 0) {
+		fprintf(stderr, "No GPU available, falling back to CPU.\n");
+		err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_CPU, 1, &device, NULL);
+		if (err < 0)
+			sysfatal("No CL devices available.\n");
+	}
 
 	/* Create context. */
 	context = clCreateContext(properties, 1, &device, NULL, NULL, &err);
diff --git a/balls.cl b/balls.cl
index e3450db..05a223b 100644
--- a/balls.cl
+++ b/balls.cl
@@ -1,24 +1,15 @@
 #define G_FACTOR 3600.0
 #define G (9.81 / G_FACTOR)
+#define DENSITY 1500.0
 
-float
-min(float a, float b) {
-	if (a < b)
-		return a;
-	return b;
-}
-
-float
-max(float a, float b) {
-	if (a > b)
-		return a;
-	return b;
-}
-
-float
-clamp(float x, float lo, float hi) {
-	return min(hi, max(x, lo));
-}
+float mass(float radius);
+int isCollision(float2 p1, float r1, float2 p2, float r2);
+void setPosition(float2 *p1, float r1, float2 *p2, float r2);
+float2 reaction(float2 p1, float2 v1, float m1, float2 p2, float2 v2, float m2);
+float2 unitNorm(float2 v);
+float fdot(float2 a, float2 b);
+float len(float2 v);
+float volume(float radius);
 
 __kernel void
 move(__global float2 *positions, __global float2 *velocities) {
@@ -63,6 +54,44 @@ collideWalls(__global float2 *positions, __global float2 *velocities, __global f
 }
 
 __kernel void
+collideBalls(
+	__global size_t *ballIndices,
+	__global float2 *positions,
+	__global float2 *velocities,
+	__global float *radii
+) {
+	size_t id, i1, i2;
+	float2 p1, p2, v1, v2;
+	float r1, r2, m1, m2;
+
+	id = get_global_id(0);
+	i1 = ballIndices[id];
+	i2 = ballIndices[id+1];
+
+	p1 = positions[i1];
+	p2 = positions[i2];
+	v1 = velocities[i1];
+	v2 = velocities[i2];
+	r1 = radii[i1];
+	r2 = radii[i2];
+	m1 = mass(r1);
+	m2 = mass(r2);
+
+	if (!isCollision(p1, r1, p2, r2))
+		return;
+
+	setPosition(&p1, r1, &p2, r2);
+
+	v1 = reaction(p1, v1, m1, p2, v2, m2);
+	v2 = reaction(p2, v2, m2, p1, v1, m1);
+
+	positions[i1] = p1;
+	positions[i2] = p2;
+	velocities[i1] = v1;
+	velocities[i2] = v2;
+}
+
+__kernel void
 genVertices(__global float2 *positions, __global float *radii, __global float2 *vertices) {
 	size_t ball, nsegs;
 	float2 center;
@@ -80,3 +109,62 @@ genVertices(__global float2 *positions, __global float *radii, __global float2 *
 
 	vertices[ball*get_local_size(0)] = center;
 }
+
+float
+mass(float radius) {
+	return volume(radius) * DENSITY;
+}
+
+/* Return true if the two balls are colliding. */
+int
+isCollision(float2 p1, float r1, float2 p2, float r2) {
+	float2 dist;
+	float rhs;
+
+	dist = p1 - p2;
+	rhs = r1 + r2;
+	return (dist.x*dist.x + dist.y*dist.y) <= rhs*rhs;
+}
+
+/* Set the positions of two balls at the moment of collision. */
+void
+setPosition(float2 *p1, float r1, float2 *p2, float r2) {
+	float2 mid, n;
+
+	mid = (*p1 + *p2) / 2.0f;
+	n = unitNorm(*p2 - *p1);
+	*p1 = mid - n*r1;
+	*p2 = mid + n*r2;
+}
+
+/* Return the velocity of ball 1 after colliding with ball 2. */
+float2
+reaction(float2 p1, float2 v1, float m1, float2 p2, float2 v2, float m2) {
+	float mrat, coef;
+	float2 dist;
+
+	mrat = 2.0 * m2 / (m1 + m2);
+	dist = p1 - p2;
+	coef = fdot(v1-v2, dist) / (len(dist)*len(dist));
+	return v1 - dist*mrat*coef;
+}
+
+float2
+unitNorm(float2 v) {
+	return v / len(v);
+}
+
+float
+fdot(float2 a, float2 b) {
+	return a.x*b.x + a.y*b.y;
+}
+
+float
+len(float2 v) {
+	return sqrt(v.x*v.x + v.y*v.y);
+}
+
+float
+volume(float radius) {
+	return 4.0 * M_PI_F * radius*radius*radius / 3.0;
+}
diff --git a/balls.h b/balls.h
index 051c185..a767db0 100644
--- a/balls.h
+++ b/balls.h
@@ -4,6 +4,22 @@ typedef struct {
 	float2 min, max;
 } Rectangle;
 
+/*
+ * A partition of the set of all possible collisions between pairs of balls.
+ * Collisions within a cell of the partition can run concurrently.  Cells must
+ * run sequentially.
+*/
+typedef struct {
+	struct cell {
+		size_t (*ballIndices)[2]; /* Array of pairs of ball indices. */
+		size_t size; /* Length of ballIndices. */
+	} *cells; /* Array of cells. */
+	size_t size; /* Length of cell array. */
+} Partition;
+
+Partition partitionCollisions(size_t nBalls);
+void freePartition(Partition part);
+
 int isCollision(float2 p1, float r1, float2 p2, float r2);
 Rectangle insetRect(Rectangle r, float n);
 
diff --git a/geo.c b/geo.c
index 84c73cf..6178bf6 100644
--- a/geo.c
+++ b/geo.c
@@ -1,3 +1,5 @@
+#include <stddef.h>
+
 #include "balls.h"
 
 int
diff --git a/partition.c b/partition.c
new file mode 100644
index 0000000..aaac220
--- /dev/null
+++ b/partition.c
@@ -0,0 +1,179 @@
+#include <stdlib.h>
+
+#include "sysfatal.h"
+#include "balls.h"
+
+/*
+ * A graph on the set of balls. Nodes are ball indices and edges are potential
+ * collisions between two balls.
+ */
+typedef struct {
+	size_t (*edges)[2]; /* Array of edges. Each edge endpoint is a ball index. */
+	size_t nEdges; /* Length of edge array. */
+} Graph;
+
+static Partition newPartition(void);
+static Graph completeGraph(size_t nBalls);
+static Graph matching(Graph g);
+static Graph addEdge(Graph g, size_t edge[2]);
+static Graph removeEdges(Graph g, Graph edges);
+static Graph removeEdge(Graph g, size_t edge[2]);
+static Partition addCell(Partition part, Graph cell);
+static Graph newGraph(void);
+static void freeGraph(Graph g);
+
+/*
+ * Partition the set of all possible collisions between pairs of balls into
+ * chunks that can be computed concurrently. Collisions within a cell of the
+ * partition can run concurrently. Cells must run sequentially.
+ */
+Partition
+partitionCollisions(size_t nBalls) {
+	Partition part;
+	Graph graph, cell;
+
+	part = newPartition();
+	graph = completeGraph(nBalls);
+	while (graph.nEdges > 0) {
+		cell = matching(graph);
+		graph = removeEdges(graph, cell);
+		part = addCell(part, cell);
+	}
+
+	freeGraph(graph);
+
+	return part;
+}
+
+void
+freePartition(Partition part) {
+	while (part.size-- > 0)
+		free(part.cells[part.size].ballIndices);
+	free(part.cells);
+}
+
+/* Allocate an empty partition. Partition should be freed by the caller after use. */
+static Partition
+newPartition(void) {
+	Partition part;
+
+	if ((part.cells = malloc(0)) == NULL)
+		sysfatal("Failed to allocate partition.\n");
+	part.size = 0;
+
+	return part;
+}
+
+/* Create a complete graph on the set of balls. */
+Graph
+completeGraph(size_t nBalls) {
+	Graph g;
+	size_t e, i, j;
+
+	g.nEdges = nBalls*(nBalls-1)/2;
+	if ((g.edges = malloc(g.nEdges * 2*sizeof(size_t))) == NULL)
+		sysfatal("Failed to allocate graph.\n");
+
+	e = 0;
+	for (i = 0; i < nBalls; i++) {
+		for (j = i+1; j < nBalls; j++) {
+			g.edges[e][0] = i;
+			g.edges[e++][1] = j;
+		}
+	}
+
+	return g;
+}
+
+/*
+ * Find and return an independent edge set in g. Two edges are independent if
+ * they don't share any vertices.
+ */
+Graph
+matching(Graph g) {
+	Graph match;
+	size_t e, m;
+
+	match = newGraph();
+	for (e = 0; e < g.nEdges; e++) {
+		for (m = 0; m < match.nEdges; m++)
+			if (g.edges[e][0] == match.edges[m][0]
+				|| g.edges[e][0] == match.edges[m][1]
+				|| g.edges[e][1] == match.edges[m][0]
+				|| g.edges[e][1] == match.edges[m][1]
+			)
+				break;
+		if (m == match.nEdges) /* No shared vertices. */
+			match = addEdge(match, g.edges[e]);
+	}
+
+	return match;
+}
+
+Graph
+addEdge(Graph g, size_t edge[2]) {
+	if ((g.edges = realloc(g.edges, (g.nEdges+1)*2*sizeof(size_t))) == NULL)
+		sysfatal("Failed to add edge to graph.\n");
+	g.edges[g.nEdges][0] = edge[0];
+	g.edges[g.nEdges][1] = edge[1];
+	g.nEdges++;
+	return g;
+}
+
+static Graph
+removeEdges(Graph g, Graph edges) {
+	size_t i;
+
+	for (i = 0; i < edges.nEdges; i++)
+		g = removeEdge(g, edges.edges[i]);
+	return g;
+}
+
+Graph
+removeEdge(Graph g, size_t edge[2]) {
+	size_t i, j;
+
+	for (i = 0; i < g.nEdges; i++) {
+		if (g.edges[i][0] == edge[0] && g.edges[i][1] == edge[1]) {
+			/* Shift left. */
+			for (j = i+1; j < g.nEdges; j++) {
+				g.edges[j-1][0] = g.edges[j][0];
+				g.edges[j-1][1] = g.edges[j][1];
+			}
+
+			/* Shrink array. */
+			if ((g.edges = realloc(g.edges, (g.nEdges-1)*2*sizeof(size_t))) == NULL)
+				sysfatal("Failed to reallocate graph edge array.\n");
+			g.nEdges--;
+
+			break;
+		}
+	}
+
+	return g;
+}
+
+Partition
+addCell(Partition part, Graph cell) {
+	if ((part.cells = realloc(part.cells, (part.size+1)*2*sizeof(size_t))) == NULL)
+		sysfatal("Failed to grow partition array.\n");
+	part.cells[part.size].ballIndices = cell.edges;
+	part.cells[part.size].size = cell.nEdges;
+	part.size++;
+	return part;
+}
+
+static Graph
+newGraph(void) {
+	Graph g;
+
+	if ((g.edges = malloc(0)) == NULL)
+		sysfatal("Failed to allocate graph.\n");
+	g.nEdges = 0;
+	return g;
+}
+
+static void
+freeGraph(Graph g) {
+	free(g.edges);
+}