技術標籤：圖形學演算法圖形學流體力學Fluid

效果圖

實現基於非常非常經典的 Paper ——Jos Stam 的Real-Time Fluid Dynamics for Games

程式碼和註釋

float visc = 0.0f;
float diff = 0.0f;

int N = 250;
int size = (N + 2) * (N + 2);
float[] u = new float[size]; // 速度, 水平/豎直
float[] v = new float[size];
float[] u_prev = new float[size];
float[] v_prev = new float[size];
float[] dens = new float[size]; // 密度
float[] dens_prev = new float[size];

float dt = 0.4f; // 時間間隔
float densityAdd = 50.0f;
float forceX = 0.0f;
float forceY = -7.0f;
float circleRadius = 1.0f;

void setup()
{
    size(252, 252);
}

// reference one dimensional array element
/// @note 二維轉一維，便於索引
int IX(int i, int j)
{
    return i + (N + 2) * j;
}

/// @note 新增密度
void addSource(float[] x, float s[], float dt)
{
    for (int index = 0; index < size; index++)
    {
        x[index] += dt * s[index];
    }
}

// diffusion solver
// possible diffusion at rate diff
/// @note 密度擴散（4 鄰域）
/// 由於是稀疏矩陣所以採用 Gauss-Seidel relaxation（鬆弛）法 求解
void diffuse(int N, int b, float[] x, float[] x0, float diff, float dt)
{
    int i, j, k;
    float a = dt * diff * N * N;
    for (k = 0; k < 20; k++) ///< Gauss-Seidel relaxation
    {
        for (i = 1; i <= N; i++)
        {
            for (j = 1; j <= N; j++)
            {
                x[IX(i, j)] = (x0[IX(i, j)] + a * (x[IX(i - 1, j)] + x[IX(i + 1, j)] + x[IX(i, j - 1)] + x[IX(i, j + 1)])) / (1 + 4 * a);
            }
        }
        setBounds(N, b, x);
    }
}

/// @note 平流
/// 當作一個粒子集合進行建模求解移動的密度值
/// 注意是沿這速度場逆向回溯求解，通過上一個時間片的密度值進行插值求解
void advect(int N, int b, float[] d, float[] d0, float[] u, float v[], float dt)
{
    int i, j, i0, j0, i1, j1;
    float x, y, s0, t0, s1, t1;
    float dt0 = dt * N;

    for (i = 1; i <= N; i++)
    {
        for (j = 1; j <= N; j++)
        {
            /// @note 柑橘速度場得到上一個時間片的位置
            x = i - dt0 * u[IX(i, j)];
            y = j - dt0 * v[IX(i, j)];

            if (x < 0.5)
            {
                x = 0.5;
            }
            if (x > N + 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;

            /// @note 插值的權重
            s1 = x - i0;
            s0 = 1 - s1;
            t1 = y - j0;
            t0 = 1 - t1;

            /// @note 四鄰域的插值
            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)]);
        }
    }

    setBounds(N, b, d);
}

void densityStep(int N, float[] x, float[] x0, float[] u, float[] v, float diff, float dt)
{
    float[] tmp;
    addSource(x, x0, dt);
    tmp = x0;
    x0 = x;
    x = tmp; // Swap x0, x
    diffuse(N, 0, x, x0, diff, dt);
    tmp = x0;
    x0 = x;
    x = tmp; // Swap Back x0, x
    advect(N, 0, x, x0, u, v, dt);
}

/// @note 投射算符（強制速度的質量不變性，mass conserving，通過 Hodge 分界實現 ）
/// 每個速度場都是質量保留場和梯度場的和
/// 由於是稀疏矩陣所以採用 Gauss-Seidel relaxation（鬆弛）法 求解 泊松方程
void project(int N, float[] u, float[] v, float[] p, float[] div)
{
    int i, j, k;
    float h = 1.0f / 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;
        }
    }

    // 邊界處理
    setBounds(N, 0, div);
    setBounds(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.0f;
            }
        }
        setBounds(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;
        }
    }

    setBounds(N, 1, u);
    setBounds(N, 2, v);
}

/// @note 速度的更新
/// 來源於三個因素 —— 加力、粘性散度、自平流
void velocityStep(int N, float[] u, float[] v, float[] u0, float[] v0, float visc, float dt)
{
    float[] tmp;

    addSource(u, u0, dt);
    addSource(v, v0, dt);

    tmp = u0;
    u0 = u;
    u = tmp; // Swap u0 u
    diffuse(N, 1, u, u0, visc, dt);

    tmp = v0;
    v0 = v;
    v = tmp; // Swap v0 v
    diffuse(N, 2, v, v0, visc, dt);
    project(N, u, v, u0, v0);

    tmp = u0;
    u0 = u;
    u = tmp; // Swap Back u0 u
    tmp = v0;
    v0 = v;
    v = tmp; // Swap Back v0 v
    advect(N, 1, u, u0, u0, v0, dt);
    advect(N, 2, v, v0, u0, v0, dt);
    project(N, u, v, u0, v0);
}

void setBounds(int N, int b, float[] x)
{
    /// @note 邊界上速度的反向
    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)];
    }

    /// @note 四個角處的插值
    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)]);
}

/// @note 跟新速度和密度
void updateSimulation(float dt)
{
    velocityStep(N, u, v, u_prev, v_prev, visc, dt);
    densityStep(N, dens, dens_prev, u, v, diff, dt);
}

void drawDensity()
{
    background(0);
    loadPixels();

    float density;
    for (int index = 0; index < size; index++)
    {
        density = dens[index] * 255.0f;
        pixels[index] = color(density, density, density);
    }
    updatePixels();
}

void drawCircle(float dt, float val, float radius)
{
    // create B&W mask and apply density where pixels > 0.0f
    background(0);
    noStroke();
    ellipse(mouseX, mouseY, radius, radius);
    loadPixels();

    float pixelVal;
    for (int index = 0; index < size; index++)
    {
        pixelVal = red(pixels[index]);

        if (pixelVal > 0.0f)
        {
            dens_prev[index] = val;
            u_prev[index] += dt * forceX; ///< 加力
            v_prev[index] += dt * forceY;
        }
    }
}

void draw()
{
    for (int index = 0; index < size; index++)
    {
        dens_prev[index] = 0;
        u_prev[index] = 0;
        v_prev[index] = 0;
    }

    if (mousePressed == true)
    {
        if (mouseButton == LEFT)
        {
            drawCircle(dt, densityAdd, circleRadius);
        }
    }

    updateSimulation(dt);
    drawDensity();
}
 

擴充套件延申

《GPU Gem3》(Chapter 38. Fast Fluid Dynamics Simulation on the GPU)

Simulation and Visualization of a 3D Fluid

流體模擬FluidSimulation：流體模擬基礎 —— 很詳盡的 Navier Stokes 方程的數學推導，很棒棒