1 files changed, 153 insertions, 0 deletions
diff --git a/kernel.cl b/kernel.cl
new file mode 100644
index 0000000..8bbe28a
--- /dev/null
+++ b/kernel.cl
@@ -0,0 +1,153 @@
+#define RADIUS      $radius
+#define N_PARTICLES $n_particles
+#define DELTA_T     $delta_t
+
+typedef float  scalar_t;
+typedef float2 vec_t;
+
+#define SCALAR_MAX FLT_MAX
+
+__constant vec_t wall_normals[4] = {
+	(vec_t)(-1.,0.),
+	(vec_t)( 1.,0.),
+	(vec_t)(0., 1.),
+	(vec_t)(0.,-1.),
+};
+
+__constant scalar_t wall_loc[4] = {
+	-RADIUS, 1.-RADIUS, 1.-RADIUS, -RADIUS
+};
+
+scalar_t solve_wall_collision(vec_t n, scalar_t loc, vec_t p, vec_t v) {
+	if (dot(n,v) > 0.) {
+		vec_t wall_v = dot(n,v) * n;
+		return (loc - dot(p,n)) / dot(wall_v,n);
+	} else {
+		return SCALAR_MAX;
+	}
+}
+
+scalar_t pos_min_or_infty(scalar_t t0, scalar_t t1) {
+	if (t0 >= 0.) {
+		if (t1 >= 0.) {
+			return min(min(t0, t1), SCALAR_MAX);
+		} else {
+			return min(t0, SCALAR_MAX);
+		}
+	} else {
+		if (t1 >= 0.) {
+			return min(t1, SCALAR_MAX);
+		} else {
+			return SCALAR_MAX;
+		}
+	}
+}
+
+scalar_t solve_collision(vec_t p, vec_t v, vec_t p_, vec_t v_) {
+	scalar_t a = dot(v-v_, v-v_);
+	scalar_t b  = 2.*dot(p-p_, v-v_);
+	scalar_t c = dot(p-p_, p-p_) - 4.*RADIUS*RADIUS;
+	scalar_t d = b*b - 4.*a*c;
+
+	if (d >= 0.) {
+		scalar_t t0 = (-b + sqrt(d))/(2.*a);
+		scalar_t t1 = (-b - sqrt(d))/(2.*a);
+
+		return pos_min_or_infty(t0, t1);
+	} else {
+		return SCALAR_MAX;
+	}
+}
+
+__kernel void evolve(__global vec_t* pos_a,
+                     __global vec_t* vel_a,
+                     __global vec_t* pos_b,
+                     __global vec_t* vel_b,
+                     __global volatile unsigned int* last_collide)
+{
+	unsigned int i = get_global_id(0);
+
+	vec_t p = pos_a[i];
+	vec_t v = vel_a[i];
+
+	unsigned int jParticle = N_PARTICLES;
+	scalar_t min2intersect = SCALAR_MAX;
+
+	for (unsigned int iParticle=0; iParticle < N_PARTICLES; ++iParticle) {
+		if (iParticle != i && !(last_collide[i] == iParticle && last_collide[iParticle] == i)) {
+			vec_t p_ = pos_a[iParticle];
+			vec_t v_ = vel_a[iParticle];
+
+			scalar_t time2intersect = solve_collision(p, v, p_, v_);
+
+			if (time2intersect < min2intersect) {
+				min2intersect = time2intersect;
+				jParticle = iParticle;
+			}
+		}
+	}
+
+	unsigned int jWall = N_PARTICLES;
+	scalar_t min2wall = SCALAR_MAX;
+
+	for (unsigned int iWall=0; iWall < 4; ++iWall) {
+		scalar_t time2wall = solve_wall_collision(wall_normals[iWall], wall_loc[iWall], p, v);
+		if (time2wall < min2wall) {
+			min2wall = time2wall;
+			jWall = iWall;
+		}
+	}
+
+	if (min2intersect < DELTA_T) {
+		if (min2wall < min2intersect) {
+			p += min2wall * v;
+			v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall];
+			p += (DELTA_T - min2wall) * v;
+			last_collide[i] = N_PARTICLES;
+
+			pos_b[i] = p;
+			vel_b[i] = v;
+		} else {
+			vec_t p_ = pos_a[jParticle];
+			vec_t v_ = vel_a[jParticle];
+
+			p  += min2intersect * v;
+			p_ += min2intersect * v_;
+
+			vec_t omega = normalize(p - p_);
+
+			while (dot(omega,omega) < 1.0) {
+				omega *= 1.001;
+			}
+
+			v  -= dot(vel_a[i] - vel_a[jParticle], omega) * omega;
+			v_ -= dot(vel_a[jParticle] - vel_a[i], omega) * omega;
+
+			if (i < jParticle) {
+				p  += (DELTA_T - min2intersect) * v;
+				p_ += (DELTA_T - min2intersect) * v_;
+
+				pos_b[i] = p;
+				vel_b[i] = v;
+
+				pos_b[jParticle] = p_;
+				vel_b[jParticle] = v_;
+
+				last_collide[i] = jParticle;
+				last_collide[jParticle] = i;
+			}
+		}
+	} else {
+		if (min2wall < DELTA_T) {
+			p += min2wall * v;
+			v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall];
+			p += (DELTA_T - min2wall) * v;
+			last_collide[i] = N_PARTICLES;
+		} else {
+			p += DELTA_T * v;
+		}
+
+		pos_b[i] = p;
+		vel_b[i] = v;
+	}
+}