Processing Fluid Motion

gistfile1.pde

import java.util.Arrays;

import processing.core.PApplet;
import processing.core.PConstants;

public class Fluid extends PApplet {
    class Physics extends Thread {
	@Override
	public void run() {
	    while (true) {
		for (int i = 0; i < n + 2; i++) {
		    Arrays.fill(u_prev[i], 0);
		    Arrays.fill(v_prev[i], 0);
		    Arrays.fill(dens_prev[i], 0);
		}

		if (mouseX > 1 && mouseX <= n && mouseY > 1 && mouseY <= n) {
		    if (mousePressed && mouseButton == PConstants.LEFT) {
			u[mouseY][mouseX] = force * (mouseX - oMouseX);
			v[mouseY][mouseX] = force * (oMouseY - mouseY);
		    }

		    if (mousePressed && mouseButton == PConstants.RIGHT) {
			dens[mouseY][mouseX] = source;
		    }
		}
		oMouseX = mouseX;
		oMouseY = mouseY;

		vel_step(visc, dt);
		dens_step(diff, dt);
	    }
	}
    }

    private static final long serialVersionUID = 1L;

    public final int n = 270;
    private int gridSize = n + 2;

    private float u[][] = new float[gridSize][gridSize];
    private float v[][] = new float[gridSize][gridSize];
    private float u_prev[][] = new float[gridSize][gridSize];
    private float v_prev[][] = new float[gridSize][gridSize];
    private float dens[][] = new float[gridSize][gridSize];
    private float dens_prev[][] = new float[gridSize][gridSize];

    private float dt;

    private float diff;

    private float visc;

    private float force;

    private float source;

    private int oMouseY;

    private int oMouseX;

    private void add_source(float array[][], float sources[][], float deltaTime) {
	for (int i = 0; i < gridSize; i++) {
	    for (int j = 0; j < gridSize; j++) {
		array[i][j] += deltaTime * sources[i][j];
	    }
	}
    }

    private void diffuse(int b, float array[][], float arrayPrev[][],
	    float diff, float deltaTime) {
	float a = deltaTime * diff * n * n;

	float div = 1 + 4 * a;

	for (int k = 0; k < 20; k++) {
	    for (int i = 1; i <= n; i++) {
		for (int j = 1; j <= n; j++) {
		    array[i][j] = (arrayPrev[i][j] + a
			    * (array[i - 1][j] + array[i + 1][j]
				    + array[i][j - 1] + array[i][j + 1]))
			    / div;
		}
	    }
	    set_bnd(b, array);
	}
    }

    private void set_bnd(int b, float array[][]) {
	if (b == 1) {
	    for (int i = 1; i <= n; i++) {
		array[0][i] = -array[1][i];
		array[n + 1][i] = -array[n][i];
		array[i][0] = array[i][1];
		array[i][n + 1] = array[i][n];
	    }
	} else if (b == 2) {
	    for (int i = 1; i <= n; i++) {
		array[0][i] = array[1][i];
		array[n + 1][i] = array[n][i];
		array[i][0] = -array[i][1];
		array[i][n + 1] = -array[i][n];
	    }
	} else {
	    for (int i = 1; i <= n; i++) {
		array[0][i] = array[1][i];
		array[n + 1][i] = array[n][i];
		array[i][0] = array[i][1];
		array[i][n + 1] = array[i][n];
	    }
	}
	array[0][0] = 0.5f * (array[1][0] + array[0][1]);
	array[0][n + 1] = 0.5f * (array[1][n + 1] + array[0][n]);
	array[n + 1][0] = 0.5f * (array[n][0] + array[n + 1][1]);
	array[n + 1][n + 1] = 0.5f * (array[n][n + 1] + array[n + 1][n]);
    }

    private void advect(int b, float densityNew[][], float densityPrev[][],
	    float velocityU[][], float velocityV[][], float deltaTime) {
	float dt0 = deltaTime * n;
	float x, y, s0, s1, t0, t1;
	int i0, i1, j0, j1;
	for (int i = 1; i <= n; i++) {
	    for (int j = 1; j <= n; j++) {
		x = i - dt0 * velocityU[i][j];
		y = j - dt0 * velocityV[i][j];

		if (x < 0.5) {
		    x = 0.5f;
		}
		if (x > n + 0.5) {
		    x = n + 0.5f;
		}

		i0 = (int) x;
		i1 = i0 + 1;

		if (y < 0.5) {
		    y = 0.5f;
		}
		if (y > n + 0.5) {
		    y = n + 0.5f;
		}

		j0 = (int) y;
		j1 = j0 + 1;
		s1 = x - i0;
		s0 = 1 - s1;
		t1 = y - j0;
		t0 = 1 - t1;
		densityNew[i][j] = s0
			* (t0 * densityPrev[i0][j0] + t1 * densityPrev[i0][j1])
			+ s1
			* (t0 * densityPrev[i1][j0] + t1 * densityPrev[i1][j1]);
	    }
	}
	set_bnd(b, densityNew);
    }

    void dens_step(float diff, float deltaTime) {
	float temp[][];
	synchronized (dens) {
	    add_source(dens, dens_prev, deltaTime);
	    temp = dens_prev;
	    dens_prev = dens;
	    dens = temp;
	    diffuse(0, dens, dens_prev, diff, deltaTime);
	    temp = dens_prev;
	    dens_prev = dens;
	    dens = temp;
	    advect(0, dens, dens_prev, u, v, deltaTime);
	}
    }

    void vel_step(float viscosity, float deltaTime) {
	float temp[][];
	add_source(u, u_prev, deltaTime);
	add_source(v, v_prev, deltaTime);
	temp = u_prev;
	u_prev = u;
	u = temp;
	diffuse(1, u, u_prev, viscosity, deltaTime);
	temp = v_prev;
	v_prev = v;
	v = temp;
	diffuse(2, v, v_prev, viscosity, deltaTime);
	project(u, v, u_prev, v_prev);
	temp = u_prev;
	u_prev = u;
	u = temp;
	temp = v_prev;
	v_prev = v;
	v = temp;
	advect(1, u, u_prev, u_prev, v_prev, deltaTime);
	advect(2, v, v_prev, u_prev, v_prev, deltaTime);
	project(u, v, u_prev, v_prev);
    }

    private void project(float velocityU[][], float velocityV[][], float p[][],
	    float div[][]) {
	float h = 1.0f / n;
	for (int i = 1; i <= n; i++) {
	    for (int j = 1; j <= n; j++) {
		div[i][j] = -0.5f
			* h
			* (velocityU[i + 1][j] - velocityU[i - 1][j]
				+ velocityV[i][j + 1] - velocityV[i][j - 1]);
		p[i][j] = 0;
	    }
	}
	set_bnd(0, div);
	set_bnd(0, p);
	for (int k = 0; k < 20; k++) {
	    for (int i = 1; i <= n; i++) {
		for (int j = 1; j <= n; j++) {
		    p[i][j] = (div[i][j] + p[i - 1][j] + p[i + 1][j]
			    + p[i][j - 1] + p[i][j + 1]) / 4;
		}
	    }
	    set_bnd(0, p);
	}
	for (int i = 1; i <= n; i++) {
	    for (int j = 1; j <= n; j++) {
		velocityU[i][j] -= 0.5f * (p[i + 1][j] - p[i - 1][j]) / h;
		velocityV[i][j] -= 0.5f * (p[i][j + 1] - p[i][j - 1]) / h;
	    }
	}
	set_bnd(1, velocityU);
	set_bnd(2, velocityV);
    }

    @Override
    public void setup() {
	size(n, n, P2D);

	dt = 0.1f;
	diff = 0.0f;
	visc = 0.0f;
	force = 5.0f;
	source = 100.0f;
	colorMode(RGB, 1.f);

	new Physics().start();
	noStroke();
	frameRate(9999);
    }

    @Override
    public void draw() {
	background(0);

	loadPixels();
	synchronized (dens) {
	    for (int i = 1; i < n; i++) {
		for (int j = 1; j < n; j++) {
		    pixels[j + n * i] = color(dens[i][j], 0, 0);
		}
	    }
	}
	updatePixels();
    }

    public static void main(String[] args) {
	PApplet.main(new String[] { "Fluid" });
    }
}