第四次作业就是流体模拟了，作业中给了若干的实现步骤，以及一些模板代码。

首先第一步，在update函数的开头，加载水面mesh的高度，然后在update的结束时，把计算后的高度更新到mesh中。这个很简单：

void Update () 
{
    Mesh mesh = GetComponent<MeshFilter> ().mesh;
    Vector3[] X    = mesh.vertices;
    float[,] new_h = new float[size, size];
    float[,] h     = new float[size, size];

    //TODO: Load X.y into h.
    for(int i = 0; i < size; i++)
    {
        for(int j = 0; j < size; j++)
        {
            h[i, j] = X[i * size + j].y;
        }
    }

    //TODO: Store h back into X.y and recalculate normal.
    for(int i = 0; i < size; i++)
    {
        for(int j = 0; j < size; j++)
        {
            X[i * size + j].y = h[i, j];
        }
    }
    mesh.vertices = X;
    mesh.RecalculateNormals();
}

第二步，当玩家按下r键时，需要在水面随机落下一个水滴，那么该位置的水面高度需要增加，为了避免水面溢出，还需要把该位置附近的水面高度减小，这样保证整体的水量不变。

if (Input.GetKeyDown ("r")) 
{
    //TODO: Add random water.
    int ri = Random.Range(0, size);
    int rj = Random.Range(0, size);
    float rh = Random.Range(0.1f, 1.0f);
    h[ri, rj] += rh;

    int neighbors = 0;
    if(ri > 0) neighbors++;
    if(ri < size - 1) neighbors++;
    if(rj > 0) neighbors++;
    if(rj < size - 1) neighbors++;

    rh /= neighbors;
    if(ri > 0) h[ri - 1, rj] -= rh;
    if(ri < size - 1) h[ri + 1, rj] -= rh;
    if(rj > 0) h[ri, rj - 1] -= rh;
    if(rj < size - 1) h[ri, rj + 1] -= rh;
}

这个时候不停按r可以看到出现很多小坑，坑里还有个小尖峰：

下一步，就要把这个变化传递出去，作业里要求使用Neumann Boundaries算法来更新，具体可以参考PPT第25页：

void Shallow_Wave(float[,] old_h, float[,] h, float [,] new_h)
{		
    for(int i = 0; i < size; i++)
    {
        for(int j = 0; j < size; j++)
        {
            new_h[i, j] = h[i, j] + damping * (h[i, j] - old_h[i, j]);
            if(i > 0) new_h[i, j] += rate * (h[i - 1, j] - h[i, j]);
            if(i < size - 1) new_h[i, j] += rate * (h[i + 1, j] - h[i, j]);
            if(j > 0) new_h[i, j] += rate * (h[i, j - 1] - h[i, j]);
            if(j < size - 1) new_h[i, j] += rate * (h[i, j + 1] - h[i, j]);
        }
    }

    for(int i = 0; i < size; i++)
    {
        for(int j = 0; j < size; j++)
        {
            old_h[i, j] = h[i, j];
            h[i, j] = new_h[i, j];
        }
    }
}

此时按r键就可以看到涟漪了：

那么接下来，就是要加上水面和方块的交互了。先要找到水面与方块相交的区域，然后计算出它们的low_h，这里计算相交可以使用Unity内置的Bounds.IntersectRay，每个水面格子从底部发射一条射线，然后判断是否和方块的包围盒相交，相交的距离就是low_h。

var go = GameObject.Find(name);
var renderer = go.GetComponent<Renderer>();
var bounds = renderer.bounds;
var pos = go.transform.position;
int grid_x = (int)(pos.x * 10 + size * 0.5f);
int grid_z = (int)(pos.z * 10 + size * 0.5f);
int li = grid_x - 6;
int ui = grid_x + 6;
int lj = grid_z - 6;
int uj = grid_z + 6;

for(int i = li; i <= ui; i++)
{
    for(int j = lj; j <= uj; j++)
    {
        if (i >= 0 && j >= 0 && i < size && j < size)
        {
            float dist = 99999;
            bounds.IntersectRay(new Ray(new Vector3(i * 0.1f - size * 0.05f, 0, j * 0.1f - size * 0.05f),
                Vector3.up), out dist);
            low_h[i, j] = dist;
        }
    }
}

然后需要计算vh，参考PPT第31页的流程：

for(int i = 0; i < size; i++)
{
    for(int j = 0; j < size; j++)
    {
        if (low_h[i, j] <= new_h[i, j])
        {
            b[i, j] = (new_h[i, j] - low_h[i, j]) / rate;
            cg_mask[i, j] = true;
        }
        else
        {
            b[i, j] = 0;
            cg_mask[i, j] = false;
            vh[i, j] = 0;
        }
    }
}

Conjugate_Gradient(cg_mask, b, vh, li, ui, lj, uj);

得到vh之后，就可以拿来迭代，计算出最终的new_h了：

for(int i = 0; i < size; i++)
{
    for(int j = 0; j < size; j++)
    {
        if (cg_mask[i, j])
        {
            vh[i, j] *= gamma;
        }
    }
}

for(int i = 0; i < size; i++)
{
    for(int j = 0; j < size; j++)
    {
        if(i > 0) new_h[i, j] += rate * (vh[i - 1, j] - vh[i, j]);
        if(i < size - 1) new_h[i, j] += rate * (vh[i + 1, j] - vh[i, j]);
        if(j > 0) new_h[i, j] += rate * (vh[i, j - 1] - vh[i, j]);
        if(j < size - 1) new_h[i, j] += rate * (vh[i, j + 1] - vh[i, j]);
    }
}

最终的效果如下：