Part 2
void OpDiffuse(int b, float []cell, float []cell0, float diff, float dt, int iter){
float a;
if (b==1) a = dt * diff * (G2 - 2) * (G3 - 2);
else if(b==2) a = dt * diff * (G1 - 2) * (G3 - 2);
else if(b==3) a = dt * diff * (G1 - 2) * (G2 - 2);
else a = dt * diff * (G1 - 2) * (G2 - 2);
LinearSolver(b, cell, cell0, a, 1+6*a, iter);
}
void OpProject(float []velocX, float []velocY, float []velocZ, float []p, float []div, int iter){
int [] ic = new int [6];
for (int k=1; k<G3-1; k++){
for (int j=1; j<G2-1; j++){
for (int i=1; i<G1-1; i++){
if (solid[IndexCell(i,j,k)]==0){
ic[0] = IndexCell(i+1, j, k);
ic[1] = IndexCell(i-1, j, k);
ic[2] = IndexCell(i, j+1, k);
ic[3] = IndexCell(i, j-1, k);
ic[4] = IndexCell(i, j, k+1);
ic[5] = IndexCell(i, j, k-1);
for (int cc = 0; cc<5; cc++){
if (solid[ic[cc]]==1) ic[cc] = IndexCell(i,j,k);
}
div[IndexCell(i,j,k)] = -0.5f*(
velocX[ic[0]] -
velocX[ic[1]] +
velocY[ic[2]] -
velocY[ic[3]] +
velocZ[ic[4]] -
velocZ[ic[5]]
) /G1;
p[IndexCell(i,j,k)] = 0;
}
}
}
}
Boundary(0, div);
Boundary(0, p);
LinearSolver(0,p,div,1,6,iter);
for (int k=1; k<G3-1; k++){
for(int j=1; j<G2-1; j++){
for(int i=1; i<G1-1; i++){
if (solid[IndexCell(i,j,k)]==0){
ic[0] = IndexCell(i+1, j, k);
ic[1] = IndexCell(i-1, j, k);
ic[2] = IndexCell(i, j+1, k);
ic[3] = IndexCell(i, j-1, k);
ic[4] = IndexCell(i, j, k+1);
ic[5] = IndexCell(i, j, k-1);
for (int cc = 0; cc<5; cc++){
if (solid[ic[cc]]==1) ic[cc] = IndexCell(i,j,k);
}
velocX[IndexCell(i,j,k)] -= 0.5f * (p[ic[0]] -
p[ic[1]]) * G1;
velocY[IndexCell(i,j,k)] -= 0.5f * (p[ic[2]] -
p[ic[3]]) * G2;
velocZ[IndexCell(i,j,k)] -= 0.5f * (p[ic[4]] -
p[ic[5]]) * G3;
}
}
}
}
Boundary(1,velocX);
Boundary(2,velocY);
Boundary(3,velocZ);
}
void OpAdvect(int b, float []d, float []d0, float []velocX, float []velocY, float []velocZ, float dt){
float i0, i1, j0, j1, k0, k1;
float dtx = dt*(G1-2);
float dty = dt*(G2-2);
float dtz = dt*(G3-2);
float s0,s1,t0,t1,u0,u1;
float tmp1,tmp2,tmp3,x,y,z;
float G1float = G1;
float G2float = G2;
float G3float = G3;
float ifloat, jfloat, kfloat;
int i,j,k;
for(k=1, kfloat=1f; k<G3-1; k++, kfloat++){
for(j=1, jfloat=1f; j<G2-1; j++, jfloat++){
for(i=1, ifloat=1f; i<G1-1; i++, ifloat++){
if(solid[IndexCell(i,j,k)]==0){
tmp1=dtx*velocX[IndexCell(i,j,k)];
tmp2=dty*velocY[IndexCell(i,j,k)];
tmp3=dtz*velocZ[IndexCell(i,j,k)];
x=ifloat-tmp1;
y=jfloat-tmp2;
z=kfloat-tmp3;
if(x<0.5f) x=0.5f;
if(x>G1float-1.5f) x=G1float-1.5f;
i0=floor(x);
i1=i0+1.0f;
if(y<0.5f) y=0.5f;
if(y>G2float-1.5f) y=G2float-1.5f;
j0=floor(y);
j1=j0+1.0f;
if(z<0.5f) z=0.5f;
if(z>G3float-1.5f) z=G3float-1.5f;
k0=floor(z);
k1=k0+1.0f;
s1=x-i0;
s0=1.0f-s1;
t1=y-j0;
t0=1.0f-t1;
u1=z-k0;
u0=1.0f-u1;
int i0i = (int)i0;
int i1i = (int)i1;
int j0i = (int)j0;
int j1i = (int)j1;
int k0i = (int)k0;
int k1i = (int)k1;
d[IndexCell(i,j,k)] =
s0 * (t0 * (u0 * d0[IndexCell(i0i,j0i,k0i)]
+u1 * d0[IndexCell(i0i, j0i, k1i)])
+(t1 * (u0 * d0[IndexCell(i0i, j1i, k0i)]
+u1 * d0[IndexCell(i0i, j1i, k1i)])))
+ s1 * (t0 * (u0 * d0[IndexCell(i1i,j0i,k0i)]
+u1 * d0[IndexCell(i1i, j0i, k1i)])
+(t1 * (u0 * d0[IndexCell(i1i, j1i, k0i)]
+u1 * d0[IndexCell(i1i, j1i, k1i)])));
}
}
}
}
Boundary(b,d);
}
void Boundary(int b, float []cell){
for(int j=1; j<G2-1; j++){
for(int i=1; i<G1-1; i++){
cell[IndexCell(i,j,0)] = b == 3 ? -cell[IndexCell(i,j,1)] : cell[IndexCell(i,j,1)];
cell[IndexCell(i,j,G3-1)] = b == 3 ? -cell[IndexCell(i,j,G3-2)] : cell[IndexCell(i,j,G3-2)];
}
}
for(int k=1; k<G3-1; k++){
for(int i=1; i<G1-1; i++){
cell[IndexCell(i,0,k)] = b == 2 ? -cell[IndexCell(i,1,k)] : cell[IndexCell(i,1,k)];
cell[IndexCell(i,G2-1,k)] = b == 2 ? -cell[IndexCell(i,G2-2,k)] : cell[IndexCell(i,G2-2,k)];
}
}
for(int k=1; k<G3-1; k++){
for(int j=1; j<G2-1; j++){
cell[IndexCell(0,j,k)] = b == 1 ? -cell[IndexCell(1,j,k)] : cell[IndexCell(1,j,k)];
cell[IndexCell(G1-1,j,k)] = b == 1 ? -cell[IndexCell(G1-2,j,k)] : cell[IndexCell(G1-2,j,k)];
}
}
//first bottom row, first corner
cell[IndexCell(0,0,0)] = 0.33f * (cell[IndexCell(1,0,0)] +
cell[IndexCell(0,1,0)] +
cell[IndexCell(0,0,1)]);
//first bottom row, second corner
cell[IndexCell(0,G2-1,0)] = 0.33f * (cell[IndexCell(1,G2-1,0)] +
cell[IndexCell(0,G2-2,0)] +
cell[IndexCell(0,G2-1,1)]);
//first top row, first corner
cell[IndexCell(0,0,G3-1)] = 0.33f * (cell[IndexCell(1,0,G3-1)] +
cell[IndexCell(0,1,G3-1)] +
cell[IndexCell(0,0,G3-2)]);
//first top row, second corner
cell[IndexCell(0,G2-1,G3-1)]= 0.33f * (cell[IndexCell(1,G2-1,G3-1)] +
cell[IndexCell(0,G2-2,G3-1)] +
cell[IndexCell(0,G2-1,G3-2)]);
//last bottom row, first corner
cell[IndexCell(G1-1,0,0)] = 0.33f * (cell[IndexCell(G1-2,0,0)] +
cell[IndexCell(G1-1,1,0)] +
cell[IndexCell(G1-1,0,1)]);
//last bottom row, second corner
cell[IndexCell(G1-1,G2-1,0)]= 0.33f * (cell[IndexCell(G1-2,G1-1,0)] +
cell[IndexCell(G1-1,G2-2,0)] +
cell[IndexCell(G1-1,G2-1,1)]);
//last top row, first corner
cell[IndexCell(G1-1,0,G3-1)]= 0.33f * (cell[IndexCell(G1-2,0,G3-1)] +
cell[IndexCell(G1-1,1,G3-1)] +
cell[IndexCell(G1-1,0,G3-2)]);
//last top row, second corner
cell[IndexCell(G1-1,G2-1,G3-1)]=0.33f * (cell[IndexCell(G1-2,G2-1,G3-1)] +
cell[IndexCell(G1-1,G2-2,G3-1)] +
cell[IndexCell(G1-1,G2-1,G3-2)]);
BoundaryInternal(b, cell);
}
void BoundaryInternal(int b, float [] cell){
int index;
for (int m=0; m<solid.length; m++){
if (solid[m]==1){
switch(b)
{
case 1:
if (solidx[m]>0){
index =IndexCell(solidx[m]-1, solidy[m], solidz[m]);
if (solid[index]==0) cell[index] = -cell[index];
}
if (solidx[m]<G1){
index=IndexCell(solidx[m]+1, solidy[m], solidz[m]);
if (solid[index]==0) cell[index] = -cell[index];
}
break;
case 2:
if (solidy[m] > 0){
index = IndexCell(solidx[m], solidy[m]-1, solidz[m]);
if (solid[index]==0) cell[index] = -cell[index];
}
if(solidy[m] < G2){
index = IndexCell(solidx[m], solidy[m]+1, solidz[m]);
if (solid[index]==0) cell[index] = -cell[index];
}
break;
case 3:
if (solidz[m] > 0){
index = IndexCell(solidx[m],solidy[m],solidz[m]-1);
if (solid[index]==0) cell[index] = -cell[index];
}
if (solidz[m] < G3){
index = IndexCell(solidx[m],solidy[m],solidz[m]+1);
if (solid[index]==0) cell[index] = -cell[index];
}
break;
}
cell[IndexCell(solidx[m], solidy[m], solidz[m])] = 0f;
}
}
}
void LinearSolver(int b, float []cell, float []cell0, float a, float c, int iter){
int [] ic = new int [6];
float cRecip =1.0f / c;
for (int m=0; m<iter; m++){
for(int k=1; k<G3-1; k++){
for(int j=1; j<G2-1; j++){
for(int i=1; i<G1-1; i++){
if (solid[IndexCell(i,j,k)] == 0){
ic[0] = IndexCell(i+1, j, k);
ic[1] = IndexCell(i-1, j, k);
ic[2] = IndexCell(i, j+1, k);
ic[3] = IndexCell(i, j-1, k);
ic[4] = IndexCell(i, j, k+1);
ic[5] = IndexCell(i, j, k-1);
for (int cc = 0; cc<5; cc++){
if (solid[ic[cc]]==1) ic[cc] = IndexCell(i,j,k);
}
cell[IndexCell(i,j,k)] =
(cell0[IndexCell(i,j,k)] + a*
(cell[ic[0]] +
cell[ic[1]] +
cell[ic[2]] +
cell[ic[3]] +
cell[ic[4]] +
cell[ic[5]]
)) * cRecip;
}
}
}
}
Boundary(b, cell);
}
}
int IndexCell(int x, int y, int z){
return (x) + (y* (G1)) + (z* (G2)*(G1));
}
}
public void gui(){
currCamMatrix = new PMatrix3D(g3.camera);
camera();
pushStyle();
strokeWeight(2);
stroke(0);
fill(0);
text("FAST FLUID DYNAMICS", width/2-80, 30);
font = createFont("Calibri Bold Italic",10,true);
textFont(font);
text("V.005 // 16th Aug. 2012", width/2-80, 40);
text("based on Jos Stam's Stable Fluids", width/2-80, 50);
text("& Mike Ash's Fluid Simulation for Dummies", width/2-80, 60);
fill(0,0,255,120);
noStroke();
rect(0,0,width,90);
if (guiON == true){
if (mouseY < 90) cam.setActive(false);
else cam.setActive(true);
}
cP5.draw();
g3.camera = currCamMatrix;
popStyle();
}
public void guiLoad(int XYZ){
//cP5.setColorLabel(0);
cP5.addSlider("Viscosity", 0, 0.01f, 10, 10, 180, 7);
cP5.addSlider("Diffusion", 0, 0.01f, 10, 20, 180, 7);
cP5.addSlider("Dt", 0, 2, 10, 30, 180, 7);
cP5.addSlider("DensityIn", 0, 5, 10, 40, 180, 7);
cP5.addSlider("VelocityIn", 0, 50, 10, 50, 180, 7);
switch (XYZ){
case 0:
cP5.addSlider("SectionXYZ", 0, FFD_GridSizeX-1, 10,60,180,7);
cP5.addToggle("SectionX", true, width-250, 30, 20, 7);
cP5.addToggle("SectionY", false, width-250, 50, 20, 7);
cP5.addToggle("SectionZ", false, width-250, 70, 20, 7);
FFD_btnVal6 = 0;
FFD_intXYZ = 0;
break;
case 1:
cP5.addSlider("SectionXYZ", 0, FFD_GridSizeY-1, 10,60,180,7);
cP5.addToggle("SectionX", false, width-250, 30, 20, 7);
cP5.addToggle("SectionY", true, width-250, 50, 20, 7);
cP5.addToggle("SectionZ", false, width-250, 70, 20, 7);
FFD_btnVal6 = 0;
FFD_intXYZ = 1;
break;
case 2:
cP5.addSlider("SectionXYZ", 0, FFD_GridSizeZ-1, 10,60,180,7);
cP5.addToggle("SectionX", false, width-250, 30, 20, 7);
cP5.addToggle("SectionY", false, width-250, 50, 20, 7);
cP5.addToggle("SectionZ", true, width-250, 70, 20, 7);
FFD_btnVal6 = 0;
FFD_intXYZ = 2;
break;
}
if (FFD_GridSizeY > FFD_GridSizeX) cP5.addSlider("IntBoundary", 1, FFD_GridSizeY-1, 10,70,180,7);
else cP5.addSlider("IntBoundary", 1, FFD_GridSizeX-1, 10,70,180,7);
cP5.addToggle("Obstacle", false, width-250, 10, 20, 7);
cP5.addToggle("VelocityArrows", false, width-170, 30, 20, 7);
cP5.addToggle("DensityCubes", false, width-90, 30, 20, 7);
cP5.setAutoDraw(false);
}
public void Viscosity(float val) {
FFD_btnVal1 = val;
}
public void Diffusion(float val) {
FFD_btnVal2 = val;
}
public void Dt(float val) {
FFD_btnVal3 = val;
}
public void DensityIn(float val) {
FFD_btnVal4 = val;
}
public void VelocityIn(float val) {
FFD_btnVal5 = val;
}
public void SectionXYZ(int val){
FFD_btnVal6 = val;
}
public void IntBoundary(int val){
FFD_btnVal7 = val;
}
public void VelocityArrows(boolean theFlag) {
if(theFlag==true) {
FFD_tglVal1 = true;
} else {
FFD_tglVal1 = false;
}
}
public void DensityCubes(boolean theFlag) {
if(theFlag==true) {
FFD_tglVal2 = true;
} else {
FFD_tglVal2 = false;
}
}
public void Obstacle(boolean theFlag) {
if(theFlag==true) {
FFD_tglVal3 = true;
} else {
FFD_tglVal3 = false;
}
}
public void SectionX(boolean theFlag) {
if(theFlag==true) {
FFD_tglVal4 = true;
FFD_tglVal5 = false;
FFD_tglVal6 = false;
guiLoad(0);
} else {
FFD_tglVal4 = false;
FFD_tglVal5 = true;
FFD_tglVal6 = true;
}
}
public void SectionY(boolean theFlag) {
if(theFlag==true) {
FFD_tglVal5 = true;
FFD_tglVal4 = false;
FFD_tglVal6 = false;
guiLoad(1);
} else {
FFD_tglVal5 = false;
FFD_tglVal4 = true;
FFD_tglVal6 = true;
}
}
public void SectionZ(boolean theFlag) {
if(theFlag==true) {
FFD_tglVal6 = true;
FFD_tglVal5 = false;
FFD_tglVal4 = false;
guiLoad(2);
} else {
FFD_tglVal6 = false;
FFD_tglVal5 = true;
FFD_tglVal4 = true;
}
}