This models the agglomeration of particles moving around in a box. The stickiness of the particles can be changed. It defines, how likely particles stick together when they meet. Particles can de-agglomerate as well, when a particle hits them. The probability for this is (1-stickiness), however the number of neighbors a particle has, reduces the probability to de-attach from the group exponentially.
One should be able to see the differences of diffusion-controlled agglomeration (stickiness high) and reaction controlled agglomeration (stickiness low).
A video can be found on Reddit.
It is also possible to make the walls sticky. If the whole surface is sticky, the particles cover the wall very quickly, if one confines the stickiness to a smaller area, structures of particles grow into the box from this area.
//License:CC BY-NC-SA 4.0
//https://creativecommons.org/licenses/by-nc-sa/4.0/
//smoother slider
Particle[] particle;
ParticleGroup[] pgroup;
Box box;
Quadtree qt;
PGraphics slider;
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
float stickiness = 0.3; //defines, how likely two particle agglomerate, 0 never, 1 always
float volpercent = 0; //volume percent of particles in box
float avgrg = 0; //average radius of gyration of all groups
float avgrh = 0; //average hydrodymanic radius
int groupcount = 0;
float slidery=map(stickiness, 0, 1, 437, 87); //stickiness slider y-coordinate
void setup() {
size(900, 700, P3D);
particle = new Particle[pnumber];
pgroup = new ParticleGroup[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; //requires restriction to at least gaussian/2 to avoid glitches
}
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);
pgroup[i] = new ParticleGroup();
particle[i] = new Particle(loc, vel, pisize, i, pgroup[i]);
pgroup[i].grouplist.add(particle[i]);
if (vel.mag() > maxvel) {
maxvel =vel.mag();
}
if (pisize > maxsize) {
maxsize =pisize;
}
volpercent += particle[i].pmass;
}
volpercent = volpercent / (boxsize*boxsize*boxsize);
//color definition for groups
for (int i=0; i< pgroup.length; i++) {
float ca = i;
float cb = i-pnumber/3;
float cc = i-2*pnumber/3;
if (i <= pnumber/3) {
pgroup[i].gr = 255;
pgroup[i].gg = 125;
pgroup[i].gb = int(map(abs(ca), 0, pnumber/3, 0, 255));
} else if (i <= 2*pnumber/3) {
pgroup[i].gr = 125;
pgroup[i].gg = int(map(abs(cb), 0, pnumber/3, 0, 255));
pgroup[i].gb = 255;
} else {
pgroup[i].gr = int(map(abs(cc), 0, pnumber/3, 0, 255));
pgroup[i].gg = 255;
pgroup[i].gb = 125;
}
}
//image of slider for particle stickiness adjustment
slider = createGraphics(86, 26);
slider.beginDraw();
slider.colorMode(RGB, 255, 255, 255);
slider.fill(255);
slider.strokeWeight(6);
slider.stroke(190);
slider.rect(3, 3, 80, 20);
slider.stroke(125);
slider.strokeWeight(2);
slider.line(43, 8, 43, 18);
slider.line(33, 8, 33, 18);
slider.line(53, 8, 53, 18);
slider.endDraw();
}
void draw() {
background(0);
perspective();
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.ppgroup.groupvel.copy();
step.mult(restime);
p.location.add(step);
}
//reset particles and quadtree
qt.insert(p); //insertion into quadtree
p.display(); //draw particles
}
//update of group parameters
avgrg=0;
avgrh=0;
groupcount=0;
for (ParticleGroup g : pgroup) {
g.pgupdate();
if (g.groupvel.mag() > maxvel) {
maxvel =g.groupvel.mag();
}
if (g.grouplist.size()>0) {
groupcount++;
avgrg += g.Rg;
avgrh += g.Rh;
}
}
avgrg /= groupcount;
avgrg /= groupcount;
//***code for testing overlapping particles. There are two types of overlapping particles:
//1. particles within the same group, which due a rounding error during collsion resolution are not exactly 1.0 * (p.size+b.size) distanced, but only 0.9999999x
//2. sometimes two particles show up which have 0.93xxxx or 0.95xxxx of (p.size+b.size) distance and are not in the same group, these pairs always create an event in the next round which resolves the issue
/*
println(frameCount);
for (Particle p : particle){
for (Particle b: particle){
if (p.pid < b.pid){
PVector pnew = p.ppgroup.groupvel.copy();
pnew.mult(restime);
pnew.add(p.location);
PVector bnew = b.ppgroup.groupvel.copy();
bnew.mult(restime);
bnew.add(b.location);
float distance = PVector.dist(pnew, bnew);
boolean samegroup = false;
if (p.ppgroup == b.ppgroup){
samegroup = true;
}
if (distance < p.size+b.size){
println(p.pid, b.pid, distance/(p.size+b.size),samegroup);
}
}
}
}
*/
//*** end of test code
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>();
for (Particle p : particle) {
eventidentification(p);
}
//collision resolution until all events are resolved
while (eventlist.size() > 0 ) {
//get lowest time event from eventlist
Event neid = extractMax(); //next event id
//store information about group of particles because they might get changed when deagglomerating
ParticleGroup neidoldga = particle[neid.aid].ppgroup;
ParticleGroup neidoldgb = particle[neid.bid].ppgroup;
//println(frameCount,"Next Event", neid.aid, neid.bid, neid.etime, neid.type);
//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] && p.ppgroup != neidoldga && p.ppgroup != neidoldgb) {
PVector ostep = p.ppgroup.groupvel.copy();
ostep.mult(neid.etime);
p.location.add(ostep);
}
qt.insert(p); //insertion into quadtree
}
//update of all group parameters
avgrg=0;
avgrh=0;
groupcount=0;
for (ParticleGroup g : pgroup) {
g.pgupdate();
if (g.groupvel.mag() > maxvel) {
maxvel =g.groupvel.mag();
}
if (g.grouplist.size()>1) {
groupcount++;
avgrg += g.Rg;
avgrh += g.Rh;
}
}
avgrg /= groupcount;
avgrh /= groupcount;
// updates the residual time in the round, then new event selection
restime -= neid.etime;
eventlist = new ArrayList<Event>();
for (Particle p : particle) {
eventidentification(p);
}
}
datascreen();
}
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 p = ep_;
//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.ppgroup.groupvel.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;
p.pneighbors = new ArrayList<Particle>();
for (Particle i : pfound) {
float distance = PVector.dist(p.location, i.location);
if (distance < 1.02*(p.size+i.size) && p.ppgroup == i.ppgroup) {
p.pneighbors.add(i);
}
if (p.pid < i.pid && p.ppgroup != i.ppgroup) { //this eliminates double counting for the first event location and excludes particles of the same group
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.ppgroup.groupvel);
if (ploc.x+p.size >= boxsize/2) {
wtimex = (boxsize/2 - p.location.x-p.size) / p.ppgroup.groupvel.x;
} else if (ploc.x-p.size <= -boxsize/2 ) {
wtimex = (-boxsize/2 - p.location.x+p.size) / p.ppgroup.groupvel.x;
}
if (ploc.y+p.size >= boxsize/2) {
wtimey = (boxsize/2 - p.location.y-p.size) / p.ppgroup.groupvel.y;
} else if (ploc.y-p.size <= -boxsize/2 ) {
wtimey = (-boxsize/2 - p.location.y+p.size) / p.ppgroup.groupvel.y;
}
if (ploc.z+p.size >= boxsize/2) {
wtimez = (boxsize/2 - p.location.z-p.size) / p.ppgroup.groupvel.z;
} else if (ploc.z-p.size <= -boxsize/2 ) {
wtimez = (-boxsize/2 - p.location.z+p.size) / p.ppgroup.groupvel.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.ppgroup.groupvel.copy();
step.mult(wtime);
for (Particle b : p.ppgroup.grouplist) {
b.location.add(step);
if (type == "wallx") {
b.velocity.x *= -1;
} else if (type == "wally") {
b.velocity.y *= -1;
} else if (type == "wallz") {
b.velocity.z *= -1;
}
}
if (p.ppgroup.grouplist.size()>1) { //deagglomeration only, when more than one particle in the group
float deaggtest = random(1);
if (deaggtest <= pow((1-stickiness), p.pneighbors.size())) { //deagglomeration less likely the more neighbors a particle is attached to
deagglomeration(p);
}
}
}
//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.ppgroup.groupvel.copy();
float ra = pa.size;
PVector B0 = pb.location.copy();
PVector vb = pb.ppgroup.groupvel.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 and resttime for events 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].ppgroup.groupvel;
PVector vb = particle[bid].ppgroup.groupvel;
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 = particle[aid].ppgroup.groupmass;
float mb = particle[bid].ppgroup.groupmass;
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
PVector atrans = PVector.sub(vanew, particle[aid].ppgroup.groupvel);
PVector btrans = PVector.sub(vbnew, particle[bid].ppgroup.groupvel);
for (Particle a : particle[aid].ppgroup.grouplist) {
a.location.add(da);
a.velocity.add(atrans);
}
for (Particle b : particle[bid].ppgroup.grouplist) {
b.location.add(db);
b.velocity.add(btrans);
}
float sticktest = random(1);
float deaggtest = random(1);
//particles get added to the largest group, when they stick
if (sticktest <= stickiness) {
ParticleGroup oldg;
ParticleGroup nextg;
if (particle[aid].ppgroup.grouplist.size() >= particle[bid].ppgroup.grouplist.size()) {
oldg = particle[bid].ppgroup;
nextg = particle[aid].ppgroup;
} else {
oldg = particle[aid].ppgroup;
nextg = particle[bid].ppgroup;
}
for (int i = oldg.grouplist.size()-1; i >= 0; i--) {
oldg.grouplist.get(i).ppgroup = nextg;
nextg.grouplist.add(oldg.grouplist.get(i));
oldg.grouplist.remove(i);
}
} else {
if (particle[aid].ppgroup.grouplist.size()>1) { //deagglomeration only, when group has more than one particle
if (deaggtest <= pow((1-stickiness), particle[aid].pneighbors.size())) { //chance for deagglomeration lower the more neighbors a particle is attached to
deagglomeration(particle[aid]);
}
}
if (particle[bid].ppgroup.grouplist.size()>1) {
if (deaggtest <= pow((1-stickiness), particle[bid].pneighbors.size())) {
deagglomeration(particle[bid]);
}
}
}
}
void deagglomeration(Particle p_) {
Particle p = p_;
if (p.ppgroup.grouplist.size()>1) {
for (ParticleGroup g : pgroup) {
boolean groupfound = false;
if (g.grouplist.size() == 0) {
int index = p.ppgroup.grouplist.indexOf(p);
g.grouplist.add(p);
p.ppgroup.grouplist.remove(index);
p.ppgroup = g; //has to happen at the end, so that all other actions are carried out correctly
groupfound = true;
}
if (groupfound) {
break;
}
}
}
}
void datascreen() {
hint(DISABLE_DEPTH_TEST); //required to make sure datascreen is not in 3D environment
camera();
ortho();
float volpv = float(int(volpercent*10000))/100;
String volps = str(volpv);
String rhs = str(int(avgrh));
String gc = str(groupcount);
//this section allows the slider button to move smoothly
if (mousePressed && mouseX > 730 && mouseX < 820) {
if (slidery <= 437 && slidery >= 87) {
int my = mouseY;
if (my > 437) {
my=437;
} else if (my < 87) {
my = 87;
}
slidery += (my-slidery)*0.3;
}
}
textSize(24);
fill(125);
textAlign(CENTER);
text("Vol%", 650, 500, 125, 50);
text("Avg. Rh,agg", 750, 500, 125, 50);
text("No. of aggregates", 650, 600, 250, 50);
textSize(28);
text(volps, 650, 550, 125, 50);
text(rhs, 750, 550, 125, 50);
text(gc, 650, 650, 250, 50);
text("Particle", 700, 20, 150, 50);
text("Stickiness", 700, 50, 150, 50);
fill(255);
rect(750, 100, 50, 350);
fill (60);
rect(750, 100, 10, 350);
rect(790, 100, 10, 350);
image(slider, 732, slidery);
hint(ENABLE_DEPTH_TEST);
}
void keyPressed() {
if (key == 't' || key == 'T') {
angle += 0.1;
}
}
void mouseReleased() {
stickiness = map(slidery, 437, 87, 0, 1);
}
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;
float pmass;
int pr = 248; //rvalue particle color
int pg = 255; //gvalue particle color
int pb = 38; //bvalue particle color
int pid = 0; //particle id
ParticleGroup ppgroup; //group the particle belongs to
ArrayList <Particle> pneighbors; //neighbors close to particle, required for deagglomeration
Particle(PVector location_, PVector velocity_, float size_, int pid_, ParticleGroup ppgroup_) {
location = location_;
velocity = velocity_;
size = size_;
pid = pid_;
ppgroup = ppgroup_;
pmass = 4/3 * PI *pow(size, 3);
pneighbors = new ArrayList<Particle>();
}
void display() {
pushMatrix();
translate(location.x, location.y, location.z);
noStroke();
if (ppgroup.grouplist.size()>1) { //group colors are used, if more than one particle are in group
fill(ppgroup.gr, ppgroup.gg, ppgroup.gb);
} else {
fill(pr, pg, pb);
}
sphere(size);
popMatrix();
}
}
class ParticleGroup {
float groupmass;
PVector com = new PVector(); //center of gravity location
float Rg; //radius of gyration of group
float Rh; //hydrodynamic radius
PVector groupvel = new PVector();
ArrayList <Particle> grouplist;
int gr = 248;
int gg = 255;
int gb=38;
ParticleGroup() {
grouplist = new ArrayList<Particle>();
}
void pgupdate() {
//mass update
groupmass =0;
for (Particle m : grouplist) {
groupmass += m.pmass;
}
//Center of mass and group velocity update
PVector comcalc = new PVector();
PVector groupvelcalc = new PVector();
float invgmass = 1 / groupmass;
for (Particle m : grouplist) {
PVector ploc = m.location.copy();
ploc.mult(m.pmass);
comcalc.add(ploc);
PVector pmvel = m.velocity.copy();
pmvel.mult(m.pmass);
groupvelcalc.add(pmvel);
}
comcalc.mult(invgmass);
com = comcalc.copy();
groupvelcalc.mult(invgmass);
groupvel = groupvelcalc.copy();
//Update of Rg and Rh
float rgcalc = 0;
float rhcalc = 0;
int mostdist = 0;
for (Particle m : grouplist) {
float rdist = PVector.dist(m.location, com);
rgcalc += m.pmass * rdist * rdist;
if (rdist > rhcalc) {
rhcalc = rdist;
mostdist = m.pid;
}
}
rhcalc +=particle[mostdist].size;
Rh=rhcalc;
rgcalc *= invgmass;
Rg = sqrt(rgcalc);
}
}
// Function to insert a new element in the Binary Heap
void pqueue_insert(Event e) {
//hsize = hsize + 1;
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;
}
}