Originally I had planned to develop that a bit further before presenting it here but now I spend weeks chasing bugs (some even imaginary) to make that work, that I decide to leave that here as one option how to do 3D particle collision in case somebody is looking for an example.
There are two versions of the code, which are slightly different. The first version does event identification for every particle after each event was resolved. The second version avoids this and only does this for the particles involved in the last event. The second runs faster with the current setup (12 fps instead of 9 fps on my old chromebook) but from time to time there is an overlap between two particles, that happens. It is rare and resolved in the next round but it shows something is not working perfectly there.
Thanks to @scudly for the help and pointing me in the right direction for some stuff.
Slower, more robust version:
//License:CC BY-NC-SA 4.0
//https://creativecommons.org/licenses/by-nc-sa/4.0/
Particle[] particle;
Box box;
Quadtree qt;
int qtcapacity = 4;
int boxsize = 400;
float angle = 45.0;
//variable parameters
float speed = 10.0; //average particle speed
float speedsd = 0.2; // standard deviation speed
float psize = 10.0; //particle size
float psizesd = 0.2; //standard deviation of size distribution
int pnumber = 200; //particle number
float maxvel = 0.0; //maximum velocity (used for quadtree)
ArrayList<Event> eventlist; //Array of priority queue
float restime = 0.0; //residual time in this round; 0 to pass first particle resolution
float maxsize = 0.0; //maximum size for particles
void setup() {
size(700, 700, P3D);
particle = new Particle[pnumber];
for (int i = 0; i < pnumber; i++) {
PVector loc, vel;
loc = new PVector();
float pisize = -1;
while (pisize < psize*(1-psizesd*3) || pisize > psize*(1+psizesd*3) || pisize <=0) {
pisize = (randomGaussian()*psizesd+1)+psize;
}
float pdist = 0;
float doubler = 1;
while (pdist < doubler) {
loc = new PVector (map(random(1), 0, 1, -boxsize/2+pisize, boxsize/2-pisize), map(random(1), 0, 1, -boxsize/2+pisize, boxsize/2-pisize), map(random(1), 0, 1, -boxsize/2+pisize, boxsize/2-pisize));
pdist = boxsize;
for (int j=0; j < i; j++) {
doubler = particle[j].size + pisize;
float pdist2 = PVector.dist(loc, particle[j].location);
if (pdist2 < pdist) {
pdist = pdist2;
}
}
}
vel = PVector.random3D();
float velmagrd = -1;
while (velmagrd < speed*(1-speedsd*3) || velmagrd > speed*(1+speedsd*3) || velmagrd <=0) {
velmagrd=(randomGaussian()+1)*speed;
}
float velimag = velmagrd*psize / pisize; //larger particles are slower
vel.setMag(velimag);
particle[i] = new Particle(loc, vel, pisize, i);
if (vel.mag() > maxvel) {
maxvel =vel.mag();
}
if (pisize > maxsize) {
maxsize =pisize;
}
}
}
void draw() {
background(0);
lights();
translate(350, 350, 0);
rotateY(angle);
noFill();
stroke(33, 203, 255);
strokeWeight(1);
box(boxsize);
//***reset of variables
box = new Box(0, 0, 0, boxsize);
maxvel = 0;
qt = new Quadtree(qtcapacity, box);
for (Particle p : particle) {
//resolve particle if time is left
if (restime > 0) {
PVector step = p.velocity.copy();
step.mult(restime);
p.location.add(step);
}
//reset particles and quadtree
if (p.velocity.mag() > maxvel) {
maxvel =p.velocity.mag();
}
qt.insert(p); //insertion into quadtree
p.display(); //draw particles
}
restime = 1; //changed here, so that in first runthrough particles are not moved above
//first collision testing for all particles and population of the event
eventlist = new ArrayList<Event>();
eventidentification(particle);
//collision resolution until all events are resolved
while (eventlist.size() > 0 ) {
//get lowest time event from eventlist
Event neid = extractMax(); //next event id
//collision resolution for particle and wall collisions
if (neid.type == "particle") {
collresolution(neid.aid, neid.bid, neid.etime);
} else {
wallresolution(neid.aid, neid.etime, neid.type);
}
//update of all other particles to the same time and reset of particles and quadtree
qt = new Quadtree(qtcapacity, box);
maxvel = 0;
for (Particle p : particle) {
if (p != particle[neid.aid] && p != particle[neid.bid]) {
PVector ostep = p.velocity.copy();
ostep.mult(neid.etime);
p.location.add(ostep);
}
if (p.velocity.mag() > maxvel) {
maxvel =p.velocity.mag();
}
qt.insert(p); //insertion into quadtree
}
// updates the residual time in the round, then new event selection
restime -= neid.etime;
eventlist = new ArrayList<Event>();
eventidentification(particle);
}
}
void mousePressed() {
angle += 0.1;
}
class Box {
float rx;
float ry;
float rz;
float rsize;
Box(float x_, float y_, float z_, float size_) {
rx = x_;
ry = y_;
rz = z_;
rsize = size_;
}
boolean contains (Particle p) {
boolean bcheck = false;
if (p.location.x >= rx - rsize/2 &&
p.location.x < rx +rsize/2 &&
p.location.y >= ry -rsize/2 &&
p.location.y < ry +rsize/2 &&
p.location.z >= rz - rsize/2 &&
p.location.z < rz +rsize/2) {
bcheck = true;
}
return bcheck;
}
boolean intersects(Box range) {
boolean icheck = false;
if (range.rx - range.rsize/2 > box.rx + box.rsize/2 ||
range.rx + range.rsize/2 < box.rx - box.rsize/2 ||
range.ry - range.rsize/2 > box.ry + box.rsize/2 ||
range.ry + range.rsize/2 < box.ry -box.rsize/2 ||
range.rz - range.rsize/2 > box.rz + box.rsize/2 ||
range.rz + range.rsize/2 < box.rz -box.rsize/2) {
icheck = false;
} else {
icheck = true;
}
return icheck;
}
}
//this tab contains all collision detection and resolution functions
//function to detect collision in selected particle array and return list of events in an ArrayList
void eventidentification(Particle[] ep_) {
Particle[] eparticle = new Particle[ep_.length];
eparticle = ep_;
for (Particle p : eparticle) {
//walltest enters wallevents; event addition in function to allow definition of wall
walltest(p);
ArrayList<Particle> pfound = new ArrayList<Particle>();
float qtboxsize = 2*maxvel + 2*maxsize + 2*p.size + 2*p.velocity.mag(); //location of other particles can maximum be one step + one radius of both particles on each side
Box range = new Box(p.location.x, p.location.y, p.location.z, qtboxsize);
qt.query(range, pfound);
float ctime = 0;
for (Particle i : pfound) {
if (p.pid < i.pid) { //this eliminates double counting for the first event location
ctime = Colltest(p, i);
if (ctime >= 0 && ctime <= restime) {
Event x = new Event(i.pid, p.pid, ctime, "particle");
pqueue_insert(x);
}
}
}
}
}
void walltest(Particle p) {
float wtimex=restime+1;
float wtimey=restime+1;
float wtimez=restime+1;
PVector ploc = PVector.add(p.location, p.velocity);
if (ploc.x+p.size >= boxsize/2) {
wtimex = (boxsize/2 - p.location.x-p.size) / p.velocity.x;
} else if (ploc.x-p.size <= -boxsize/2 ) {
wtimex = (-boxsize/2 - p.location.x+p.size) / p.velocity.x;
}
if (ploc.y+p.size >= boxsize/2) {
wtimey = (boxsize/2 - p.location.y-p.size) / p.velocity.y;
} else if (ploc.y-p.size <= -boxsize/2 ) {
wtimey = (-boxsize/2 - p.location.y+p.size) / p.velocity.y;
}
if (ploc.z+p.size >= boxsize/2) {
wtimez = (boxsize/2 - p.location.z-p.size) / p.velocity.z;
} else if (ploc.z-p.size <= -boxsize/2 ) {
wtimez = (-boxsize/2 - p.location.z+p.size) / p.velocity.z;
}
if ((wtimex <= wtimey && wtimex <= wtimez) && wtimex <= restime) {
Event a = new Event(p.pid, p.pid, wtimex, "wallx");
pqueue_insert(a);
} else if ((wtimey < wtimex && wtimey <= wtimez) && wtimey <= restime) {
Event a = new Event(p.pid, p.pid, wtimey, "wally");
pqueue_insert(a);
} else if (wtimez <= restime) {
Event a = new Event(p.pid, p.pid, wtimez, "wallz");
pqueue_insert(a);
}
}
void wallresolution(int id, float wtime, String ttype) {
Particle p = particle[id];
String type = ttype;
PVector step = p.velocity.copy();
step.mult(wtime);
p.location.add(step);
if (type == "wallx") {
p.velocity.x *= -1;
} else if (type == "wally") {
p.velocity.y *= -1;
} else {
p.velocity.z *= -1;
}
}
//collision test function returning time of collision of both particles, if any
float Colltest (Particle pa, Particle pb) {
PVector A0 = pa.location.copy();
PVector va = pa.velocity.copy();
float ra = pa.size;
PVector B0 = pb.location.copy();
PVector vb = pb.velocity.copy();
float rb = pb.size;
float u0=restime+1; //normalized time of first collision
float u1; //mormalized time of second collision
PVector AB = PVector.sub(B0, A0);
PVector vab = PVector.sub(vb, va); //relative velocity (in normalized time)
float rab = ra+rb;
float a = vab.dot(vab); //u*u coefficient
float b = 2*vab.dot(AB); //u coefficient
float c = AB.dot(AB) - rab*rab; //constant term
//check if they're currently overlapping (this should never happen)
if ( AB.dot(AB) <= rab*rab ) {
u0 = 0;
u1 = 0;
//collision test would be considered true
}
// Check if they hit each other (u0, u1, real numbers) or not (u0,u1 complex numbers)
float q = b*b - 4*a*c;
if ( q >= 0 ) {
float sq = sqrt(q);
float d = 1 / (2*a);
u0 = ( -b + sq ) * d;
u1 = ( -b - sq ) * d;
if ( u0 > u1 ) {
//swap the values of u0 and u1
float temp = u0;
u0=u1;
u1 = temp;
//collision is true, lowest time is reported
} else {
//complex roots, that means collision is false
}
}
return u0;
}
//collision resolution of two particles, updating particle location, velocity for this particle
void collresolution(int aid, int bid, float u0) {
PVector A0 = particle[aid].location;
PVector B0 = particle[bid].location;
PVector va = particle[aid].velocity;
PVector vb = particle[bid].velocity;
float ra = particle[aid].size;
float rb = particle[bid].size;
PVector da = new PVector();
da = PVector.mult(va, u0);
PVector A1 = new PVector();
A1 = PVector.add(A0, da);
PVector db = new PVector();
db = PVector.mult(vb, u0);
PVector B1 = new PVector();
B1 = PVector.add(B0, db);
float ma = 4/3 * PI *pow(ra, 3); //density equals 1
float mb = 4/3 * PI *pow(rb, 3); //density equals 1
PVector vecx = PVector.sub(A1, B1);
vecx.normalize();
PVector vecv1 = va.copy();
float x1 = PVector.dot(vecx, vecv1);
PVector vecv1x = new PVector();
PVector.mult(vecx, x1, vecv1x);
PVector vecv1y = PVector.sub(vecv1, vecv1x);
PVector vecx2 = vecx.mult(-1);
PVector vecv2 = vb.copy();
float x2 = vecx2.dot(vecv2);
PVector vecv2x = new PVector();
PVector.mult(vecx2, x2, vecv2x);
PVector vecv2y = PVector.sub(vecv2, vecv2x);
float consta1 = (ma - mb)/(ma+mb);
float consta2 = 2* mb / (ma+mb);
float constb1 = 2 * ma / (ma+mb);
float constb2 = (mb - ma)/(ma+mb);
PVector vanew = new PVector();
PVector.mult(vecv1x, consta1, vanew);
PVector inta1 = new PVector();
PVector.mult(vecv2x, consta2, inta1);
vanew.add(inta1);
vanew.add(vecv1y);
PVector vbnew = new PVector();
PVector.mult(vecv1x, constb1, vbnew);
PVector intb1 = new PVector();
PVector.mult(vecv2x, constb2, intb1);
vbnew.add(intb1);
vbnew.add(vecv2y);
//update of new particle parameters
particle[aid].location = A1.copy();
particle[bid].location = B1.copy();
particle[aid].velocity = vanew.copy();
particle[bid].velocity = vbnew.copy();
}
class Event {
int aid; //ID of particle A
int bid; //ID of particle B
float etime; //time of event in total timeline of round
String type = "particle"; //defines if event is collision with particle or wallx, wally, wallz
Event(int aid_, int bid_, float etime_, String type_) {
aid=aid_;
bid=bid_;
etime=etime_;
type = type_;
}
}
class Particle {
PVector location;
PVector velocity;
float size;
int pr = 248; //rvalue particle color
int pg = 255; //gvalue particle color
int pb = 38; //bvalue particle color
int pid = 0; //particle id
Particle(PVector location_, PVector velocity_, float size_, int pid_) {
location = location_;
velocity = velocity_;
size = size_;
pid = pid_;
}
void display() {
pushMatrix();
translate(location.x, location.y, location.z);
noStroke();
fill(pr, pg, pb);
sphere(size);
popMatrix();
}
}
// Function to insert a new element in the Binary Heap
void pqueue_insert(Event e) {
eventlist.add(e);
// Shift Up to maintain heap property
shiftUp(eventlist.size()-1);
}
// Function to extract the element with maximum priority
Event extractMax() {
Event result = eventlist.get(0);
// Replace the value at the root with the last leaf
eventlist.set(0, eventlist.get(eventlist.size()-1));
//remove entry
eventlist.remove(eventlist.size()-1);
// Shift down the replaced element to maintain the heap property
shiftDown(0);
return result;
}
// Function to get value of the current maximum element
Event getMax() {
return eventlist.get(0);
}
// Function to remove the element located at given index
void pqueue_remove(int i) {
eventlist.get(i).etime = -1;
// Shift the node to the root of the heap
shiftUp(i);
// Extract the node
extractMax();
}
// Function to return the index of the parent node of a given node
int parent(int i) {
return (i - 1) / 2;
}
// Function to return the index of the left child of the given node
int leftChild(int i) {
return ((2 * i) + 1);
}
// Function to return the index of the right child of the given node
int rightChild(int i) {
return ((2 * i) + 2);
}
// Function to shift up the node in order to maintain the heap property
void shiftUp(int i) {
while (i > 0 &&
(eventlist.get(parent(i)).etime > eventlist.get(i).etime ||
(eventlist.get(parent(i)).etime == eventlist.get(i).etime && eventlist.get(parent(i)).type != "particle")) ) {
// Swap parent and current node
swap(parent(i), i);
// Update i to parent of i
i = parent(i);
}
}
// Function to shift down the node in order to maintain the heap property
void shiftDown(int i) {
int maxIndex = i;
// Left Child
int l = leftChild(i);
if (l <= eventlist.size()-1 &&
(eventlist.get(l).etime < eventlist.get(maxIndex).etime ||
(eventlist.get(l).etime == eventlist.get(maxIndex).etime && eventlist.get(maxIndex).type != "particle")) ) {
maxIndex = l;
}
// Right Child
int r = rightChild(i);
if (r <= eventlist.size()-1 &&
(eventlist.get(r).etime < eventlist.get(maxIndex).etime ||
(eventlist.get(r).etime == eventlist.get(maxIndex).etime && eventlist.get(maxIndex).type != "particle")) ) {
maxIndex = r;
}
// If i not same as maxIndex
if (i != maxIndex) {
swap(i, maxIndex);
shiftDown(maxIndex);
}
}
void swap(int i, int j) {
Event temp= eventlist.get(i);
eventlist.set(i, eventlist.get(j));
eventlist.set(j, temp);
}
// This code is contributed by 29AjayKumar
//https://www.geeksforgeeks.org/priority-queue-using-binary-heap/
class Quadtree {
ArrayList<Particle> pmembers = new ArrayList<Particle>();
int qtcapacity;
Box box = new Box(0, 0, 0, 0);
boolean divided = false;
Quadtree northwestfront;
Quadtree northeastfront;
Quadtree southwestfront;
Quadtree southeastfront;
Quadtree northwestback;
Quadtree northeastback;
Quadtree southwestback;
Quadtree southeastback;
Quadtree(int cap_, Box box_) {
qtcapacity = cap_;
box = box_;
}
void subdivide() {
Box nwf = new Box(box.rx-box.rsize/4, box.ry-box.rsize/4, box.rz+box.rsize/4, box.rsize/2);
northwestfront = new Quadtree(qtcapacity, nwf);
Box nef = new Box(box.rx+box.rsize/4, box.ry-box.rsize/4, box.rz+box.rsize/4, box.rsize/2);
northeastfront = new Quadtree(qtcapacity, nef);
Box swf = new Box(box.rx-box.rsize/4, box.ry+box.rsize/4, box.rz+box.rsize/4, box.rsize/2);
southwestfront = new Quadtree(qtcapacity, swf);
Box sef = new Box(box.rx+box.rsize/4, box.ry+box.rsize/4, box.rz+box.rsize/4, box.rsize/2);
southeastfront = new Quadtree(qtcapacity, sef);
Box nwb = new Box(box.rx-box.rsize/4, box.ry-box.rsize/4, box.rz-box.rsize/4, box.rsize/2);
northwestback = new Quadtree(qtcapacity, nwb);
Box neb = new Box(box.rx+box.rsize/4, box.ry-box.rsize/4, box.rz-box.rsize/4, box.rsize/2);
northeastback = new Quadtree(qtcapacity, neb);
Box swb = new Box(box.rx-box.rsize/4, box.ry+box.rsize/4, box.rz-box.rsize/4, box.rsize/2);
southwestback = new Quadtree(qtcapacity, swb);
Box seb = new Box(box.rx+box.rsize/4, box.ry+box.rsize/4, box.rz-box.rsize/4, box.rsize/2);
southeastback = new Quadtree(qtcapacity, seb);
}
void insert(Particle particle_) {
if (box.contains(particle_) == false) {
return;
}
if (pmembers.size() <= qtcapacity) {
pmembers.add(particle_);
} else {
if (divided == false) {
subdivide();
divided = true;
}
if (divided == true) {
northwestfront.insert(particle_);
northeastfront.insert(particle_);
southwestfront.insert(particle_);
southeastfront.insert(particle_);
northwestback.insert(particle_);
northeastback.insert(particle_);
southwestback.insert(particle_);
southeastback.insert(particle_);
}
}
}
ArrayList<Particle> query(Box range, ArrayList<Particle> pfound ) {
if (box.intersects(range) == false) {
} else {
for (Particle i : pmembers) {
if (range.contains(i)) {
pfound.add(i);
}
}
if (divided == true) {
northwestfront.query(range, pfound);
northeastfront.query(range, pfound);
southwestfront.query(range, pfound);
southeastfront.query(range, pfound);
northwestback.query(range, pfound);
northeastback.query(range, pfound);
southwestback.query(range, pfound);
southeastback.query(range, pfound);
}
}
return pfound;
}
}