Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# seibe/FluidMath.cs

Created Apr 27, 2015
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); } }
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.