Skip to content

Instantly share code, notes, and snippets.

@seibe
Created April 27, 2015 06:44
Show Gist options
  • Save seibe/c3e68ccf5c910cfe9e42 to your computer and use it in GitHub Desktop.
Save seibe/c3e68ccf5c910cfe9e42 to your computer and use it in GitHub Desktop.
Stable Fluid for C# ()
/**
* Jos Stam's Stable Fluids (for C#)
**/
public static class FluidMath
{
public static void VelStep(int N, ref float[] u, ref float[] v, ref float[] u0, ref float[] v0, float visc, float dt)
{
AddSource(N, ref u, ref u0, dt);
AddSource(N, ref v, ref v0, dt);
Swap(ref u0, ref u);
Diffuse(N, 1, ref u, ref u0, visc, dt);
Swap(ref v0, ref v);
Diffuse(N, 2, ref v, ref v0, visc, dt);
Project(N, ref u, ref v, ref u0, ref v0);
Swap(ref u0, ref u);
Swap(ref v0, ref v);
Advect(N, 1, ref u, ref u0, ref u0, ref v0, dt);
Advect(N, 2, ref v, ref v0, ref u0, ref v0, dt);
Project(N, ref u, ref v, ref u0, ref v0);
}
public static void DensStep(int N, ref float[] x, ref float[] x0, ref float[] u, ref float[] v, float diff, float dt)
{
AddSource(N, ref x, ref x0, dt);
Swap(ref x0, ref x);
Diffuse(N, 0, ref x, ref x0, diff, dt);
Swap(ref x0, ref x);
Advect(N, 0, ref x, ref x0, ref u, ref v, dt);
}
public static int IX(int N, int i, int j)
{
return i + (N + 2) * j;
}
private static void Swap(ref float[] x, ref float[] x0)
{
float[] tmp = x0;
x0 = x;
x = tmp;
}
private static void AddSource(int N, ref float[] x, ref float[] s, float dt)
{
int size = (N + 2) * (N + 2);
for (int i = 0; i < size; i++) x[i] += dt * s[i];
}
private static void SetBnd(int N, int b, ref float[] x)
{
for (int i = 1; i <= N; i++)
{
x[IX(N, 0 , i)] = b == 1 ? -x[IX(N, 1, i)] : x[IX(N, 1, i)];
x[IX(N, N + 1, i)] = b == 1 ? -x[IX(N, N, i)] : x[IX(N, N, i)];
x[IX(N, i , 0)] = b == 2 ? -x[IX(N, i, 1)] : x[IX(N, i, 1)];
x[IX(N, i, N + 1)] = b == 2 ? -x[IX(N, i, N)] : x[IX(N, i, N)];
}
x[IX(N, 0 , 0)] = 0.5f * (x[IX(N, 1, 0)] + x[IX(N, 0, 1)]);
x[IX(N, 0 , N + 1)] = 0.5f * (x[IX(N, 1, N + 1)] + x[IX(N, 0, N)]);
x[IX(N, N + 1 , 0)] = 0.5f * (x[IX(N, N, 0)] + x[IX(N, N + 1, 1)]);
x[IX(N, N + 1, N + 1)] = 0.5f * (x[IX(N, N, N + 1)] + x[IX(N, N + 1, N)]);
}
private static void LinSolve(int N, int b, ref float[] x, ref float[] x0, float a, float c)
{
for (int k = 0; k < 20; k++)
{
for (int i = 1; i <= N; i++)
{
for (int j = 1; j <= N; j++)
{
x[IX(N, i, j)] = (x0[IX(N, i, j)] + a * (x[IX(N, i - 1, j)] + x[IX(N, i + 1, j)] + x[IX(N, i, j - 1)] + x[IX(N, i, j + 1)])) / c;
}
}
SetBnd(N, b, ref x);
}
}
private static void Diffuse(int N, int b, ref float[] x, ref float[] x0, float diff, float dt)
{
float a = dt * diff * N * N;
LinSolve(N, b, ref x, ref x0, a, 1 + 4 * a);
}
private static void Advect(int N, int b, ref float[] d, ref float[] d0, ref float[] u, ref float[] v, float dt)
{
int i0, j0, i1, j1;
float x, y, s0, t0, s1, t1, dt0;
dt0 = dt * N;
for (int i = 1; i <= N; i++)
{
for (int j = 1; j <= N; j++)
{
x = i - dt0 * u[IX(N, i, j)]; y = j - dt0 * v[IX(N, i, j)];
if (x < 0.5f) x = 0.5f; if (x > N + 0.5f) x = N + 0.5f; i0 = (int)x; i1 = i0 + 1;
if (y < 0.5f) y = 0.5f; if (y > N + 0.5f) y = N + 0.5f; j0 = (int)y; j1 = j0 + 1;
s1 = x - i0; s0 = 1 - s1; t1 = y - j0; t0 = 1 - t1;
d[IX(N, i, j)] = s0 * (t0 * d0[IX(N, i0, j0)] + t1 * d0[IX(N, i0, j1)]) + s1 * (t0 * d0[IX(N, i1, j0)] + t1 * d0[IX(N, i1, j1)]);
}
}
SetBnd(N, b, ref d);
}
private static void Project(int N, ref float[] u, ref float[] v, ref float[] p, ref float[] div)
{
int i, j;
for (i = 1; i <= N; i++)
{
for (j = 1; j <= N; j++)
{
div[IX(N, i, j)] = -0.5f * (u[IX(N, i + 1, j)] - u[IX(N, i - 1, j)] + v[IX(N, i, j + 1)] - v[IX(N, i, j - 1)]) / N;
p[IX(N, i, j)] = 0;
}
}
SetBnd(N, 0, ref div);
SetBnd(N, 0, ref p);
LinSolve(N, 0, ref p, ref div, 1, 4);
for (i = 1; i <= N; i++)
{
for (j = 1; j <= N; j++)
{
u[IX(N, i, j)] -= 0.5f * N * (p[IX(N, i + 1, j)] - p[IX(N, i - 1, j)]);
v[IX(N, i, j)] -= 0.5f * N * (p[IX(N, i, j + 1)] - p[IX(N, i, j - 1)]);
}
}
SetBnd(N, 1, ref u);
SetBnd(N, 2, ref v);
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment