// A Processing Implementation of Jos Stam's ``Real-Time Fluid Dynamics for Games'' fluid simulation method // The code is modified from Stam's excellent C code // // Right-Click to add density // Left-Click and drag to add velocity // Density is visualised. // Light particles are added to help visualise to velocity field. // // BP160208 // visualisation final int WIDTH = 200; final int HEIGHT = 200; // simulation variables // N = dimensions of simulation grid // DT = delta time // DIFF = diffusion coefficient // VISC = viscosity // FORCE = force at which velocity is applied via right mouse drag // SOURCE = density strength applied via left mouse click // SOURCE_RADIUS = radius of density application final int N = 40; final float DT = 0.1; final float DIFF = 0.001; final float VISC = 0; final float FORCE = 5; final float SOURCE = 5; final int SOURCE_RADIUS = 3; // Particles are used to visualise velocity field // They are not necessary for simulation // PART_RADIUS = visual radius of particle // JITTER = amount of brownian motion applied to particles (in pixels) // WIND_EFFECT = a large value indicates that the particles are very light and sensitive to the velocity field final int NUM_PARTICLES = 400; final float PART_RADIUS = 1.5; final int JITTER = (int)(PART_RADIUS); final float WIND_EFFECT = 500; // INVARIANTS / SYSTEM VARS final float DW = 1.0*WIDTH / N; final float DH = 1.0*HEIGHT / N; int omx, omy; final int SIZE = (N+2)*(N+2); float[] u = new float [SIZE]; float[] u_prev = new float [SIZE]; float[] v = new float [SIZE]; float[] v_prev = new float [SIZE]; float[] dens = new float [SIZE]; float[] dens_prev = new float [SIZE]; int IX(int i, int j){ return i+(N+2)*j; } void add_source ( int N, float[] x, float[] s, float dt ) { int i, size=(N+2)*(N+2); for ( i=0 ; iN+0.5) x=N+ 0.5; i0=(int)x; i1=i0+1; if (y<0.5) y=0.5; if (y>N+0.5) y=N+ 0.5; j0=(int)y; j1=j0+1; s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1; d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)]+t1*d0[IX(i0,j1)])+ s1*(t0*d0[IX(i1,j0)]+t1*d0[IX(i1,j1)]); } } set_bnd ( N, b, d ); } void dens_step ( int N, float[] x, float[] x0, float[] u, float[] v, float diff, float dt ) { add_source ( N, x, x0, dt ); float[] tmp = x; x = x0; x0 = tmp; // SWAP ( x0, x ); diffuse ( N, 0, x, x0, diff, dt ); tmp = x; x = x0; x0 = tmp; // SWAP ( x0, x ); advect ( N, 0, x, x0, u, v, dt ); } void vel_step ( int N, float[] u, float[] v, float[] u0, float[] v0, float visc, float dt ) { add_source ( N, u, u0, dt ); add_source ( N, v, v0, dt ); float[] tmp = u; u = u0; u0 = tmp; diffuse ( N, 1, u, u0, visc, dt ); tmp = v; v = v0; v0 = tmp; diffuse ( N, 2, v, v0, visc, dt ); project ( N, u, v, u0, v0 ); tmp = u; u = u0; u0 = tmp; tmp = v; v = v0; v0 = tmp; advect ( N, 1, u, u0, u0, v0, dt ); advect ( N, 2, v, v0, u0, v0, dt ); project ( N, u, v, u0, v0 ); } void project ( int N, float[] u, float[] v, float[] p, float[] div ) { int i, j, k; float h; h = 1.0/N; for ( i=1 ; i<=N ; i++ ) { for ( j=1 ; j<=N ; j++ ) { div[IX(i,j)] = -0.5*h*(u[IX(i+1,j)]-u[IX(i-1,j)]+ v[IX(i,j+1)]-v[IX(i,j-1)]); p[IX(i,j)] = 0; } } set_bnd ( N, 0, div ); set_bnd ( N, 0, p ); for ( k=0 ; k<20 ; k++ ) { for ( i=1 ; i<=N ; i++ ) { for ( j=1 ; j<=N ; j++ ) { p[IX(i,j)] = (div[IX(i,j)]+p[IX(i-1,j)]+p[IX(i+1,j)]+ p[IX(i,j-1)]+p[IX(i,j+1)])/4; } } set_bnd ( N, 0, p ); } for ( i=1 ; i<=N ; i++ ) { for ( j=1 ; j<=N ; j++ ) { u[IX(i,j)] -= 0.5*(p[IX(i+1,j)]-p[IX(i-1,j)])/h; v[IX(i,j)] -= 0.5*(p[IX(i,j+1)]-p[IX(i,j-1)])/h; } } set_bnd ( N, 1, u ); set_bnd ( N, 2, v ); } void set_bnd ( int N, int b, float[] x ) { for (int i=1 ; i<=N ; i++ ) { x[IX(0,i)]=(b==1)?-x[IX(1,i)]:x[IX(1,i)]; x[IX(N+1,i)]=(b==1)?-x[IX(N,i)]:x[IX(N,i)]; x[IX(i,0)]=(b==2)?-x[IX(i,1)]:x[IX(i,1)]; x[IX(i,N+1)]=(b==2)?-x[IX(i,N)]:x[IX(i,N)]; } x[IX(0 ,0 )] = 0.5*(x[IX(1,0 )]+x[IX(0 ,1)]); x[IX(0 ,N+1)] = 0.5*(x[IX(1,N+1)]+x[IX(0 ,N )]); x[IX(N+1,0 )] = 0.5*(x[IX(N,0 )]+x[IX(N+1,1)]); x[IX(N+1,N+1)] = 0.5*(x[IX(N,N+1)]+x[IX(N+1,N )]); } void draw_dens(int n, float[] dens) { for(int i=0;iN || j<1 || j>N ) return; if ( mouseButton == LEFT ) { u[IX(i,j)] = FORCE * (mx-omx); v[IX(i,j)] = FORCE * (my-omy); } if ( mouseButton == RIGHT ) { set_dens(dens,i,j,SOURCE_RADIUS,SOURCE); } omx = mx; omy = my; } void set_dens(float[] dens, int i, int j, int r, float am) { for(int ii=-r;ii <= r;ii++) { for(int jj=-r;jj<=r;jj++) { int ni = i+ii, nj = j+jj; if (ni < 0) ni += N; else if (ni >= N) ni -= N; if (nj < 0) nj += N; else if (nj >= N) nj -= N; dens[IX(ni,nj)] = am; } } } // PARTICLE float[] px = new float[NUM_PARTICLES]; float[] py = new float[NUM_PARTICLES]; void part_setup() { for(int i=0;i= WIDTH) dip1 = 1000; else dip1 = dens[IX(ip1,j)]; if (im1 < 0) dim1 = 1000; else dim1 = dens[IX(im1,j)]; if (jp1 >= HEIGHT) djp1 = 1000; else dip1 = dens[IX(i,jp1)]; if (jm1 < 0) djm1 = 1000; else djm1 = dens[IX(i,jm1)]; if (dipl < dm1) { x += dt*(dip1 - di); } */ float vx = u[IX(i,j)], vy = v[IX(i,j)]; x += vx*DT*WIND_EFFECT; y += vy*DT*WIND_EFFECT; x += JITTER - random(2*JITTER); y += JITTER - random(2*JITTER); if (x < 0) x = 0; if (x >= WIDTH) x = WIDTH; if (y < 0) y = 0; if (y >= HEIGHT) y = HEIGHT; px[p] = x; py[p] = y; } } void draw_particle(int i) { ellipse(px[i],py[i],PART_RADIUS,PART_RADIUS); } void part_draw() { fill(100,250,50); for(int i=0;i