作业二其实就是要使用隐式积分和PBD两种方式来实现布料求解。难度相对于作业一来说要简单一些，在文档中基本把步骤都写清楚了。主要逻辑首先参考Lecture 05 PPT的第18页：

然后我们按照文档的步骤一步一步地来。注意0号顶点和20号顶点是不参与更新的，它们相当于就是定死了位置。第一步是初始化，这个很简单：

for (int i = 0; i < V.Length; i++)
{
    if (!Skip_Update(i))
    {
        V[i] *= damping;
    }
}
for(int i = 0; i < X.Length; i++)
{
    if (!Skip_Update(i))
    {
        X_hat[i] = X[i] + V[i] * t;
        X[i] = X_hat[i];
    }
}

第二步就是计算梯度。这里的梯度就是PPT里的
$$\nabla F({\text{x}}^{(k)})=\frac{1}{\Delta t^2}\text{M}({\text{x}}^{(k)}-{\text{x}}^{[0]}-\Delta t^2{\text{v}}^{[0]})-\text{f}({\text{x}}^{(k)})$$
这里的f就是重力和弹簧间的弹力。重力是个常量好办，弹力的计算需要遍历所有的边，找到边的两个顶点，分别进行计算，参考PPT的第11页：

void Get_Gradient(Vector3[] X, Vector3[] X_hat, float t, Vector3[] G)
{
    //Momentum and Gravity.
    for (int i = 0; i < G.Length; i++)
    {
        if (!Skip_Update(i))
        {
            G[i] = mass * (X[i] - X_hat[i]) / (t * t) - mass * gravity;
        }
    }

    //Spring Force.
    for (int e = 0; e < L.Length; e++)
    {
        int i = E[e * 2 + 0];
        int j = E[e * 2 + 1];
        float x = Vector3.Distance(X[i], X[j]);
        Vector3 f = spring_k * (1 - L[e] / x) * (X[i] - X[j]);

        if (!Skip_Update(i))
        {
            G[i] += f;
        }

        if (!Skip_Update(j))
        {
            G[j] -= f;
        }
    }
}

然后需要计算：
$$\frac{{\partial }^{2}F({\text{x}}^{(k)})}{\partial {\text{x}}^2}=\frac{1}{\Delta t^2}\text{M}+\text{H}({\text{x}}^{(k)})$$
作业文档中给了近似求解的方法，我们就不用算$\text{H}$了。我们使用Chebyshev加速牛顿法迭代，可以参考PPT中第26页：

作业中是固定的32次迭代次数，这里的break判断就可以省略掉；另外，x的更新直接使用作业里的公式即可：

for(int k=0; k<32; k++)
{
    Get_Gradient(X, X_hat, t, G);

    if(k == 0)
    {
        omega = 1.0f;
    }
    else if(k == 1)
    {
        omega = 2.0f / (2.0f - rho * rho);
    }
    else
    {
        omega = 4.0f / (4.0f - rho * rho * omega);
    }
    
    //Update X by gradient.
    for(int i = 0; i < X.Length; i++)
    {
        if (!Skip_Update(i))
        {
            Vector3 old = X[i];
            X[i] = omega * (X[i] - 1 / (mass / (t * t) + 4 * spring_k) * G[i]) + (1 - omega) * last_X[i];
            last_X[i] = old;
        }
    }
}

迭代完别忘记更新下V，这里要使用+=，因为一开始算$\tilde{x}$的时候加过v了：

for(int i = 0; i < X.Length; i++)
{
    if (!Skip_Update(i))
    {
        V[i] += (X[i] - X_hat[i]) / t;
    }
}

最后的碰撞检测很简单，算一下点到球心的距离，如果小于半径就说明发生碰撞：

void Collision_Handling()
{
    Mesh mesh = GetComponent<MeshFilter> ().mesh;
    Vector3[] X = mesh.vertices;

    //Handle colllision.
    Vector3 c = sphere.transform.position;
    for(int i = 0; i < X.Length; i++)
    {
        if(!Skip_Update(i))
        {
            float d = Vector3.Distance(X[i], c);
            if (d < r)
            {
                V[i] = V[i] + 1 / t * (c + r * (X[i] - c) / d - X[i]);
                X[i] = c + r * (X[i] - c) / d;
            }
        }
    }

    mesh.vertices = X;
}

最后效果如下：

然后我们再看下PBD的实现。第一步就是让每个顶点自由更新：

for(int i=0; i<X.Length; i++)
{
    if(i==0 || i==20)	continue;
    //Initial Setup
    //...
    V[i] *= damping;
    V[i] += gravity * t;
    X[i] += V[i] * t;
}

接下来，通过若干次迭代施加约束，这里可以参考Lecture 06 PPT的第10页：

作业里$\alpha$取的0.2：

void Strain_Limiting()
{
    Mesh mesh = GetComponent<MeshFilter> ().mesh;
    Vector3[] vertices = mesh.vertices;

    //Apply PBD here.
    //...
    Vector3[] sum_X = new Vector3[vertices.Length];
    int[] sum_N = new int[vertices.Length];

    for (int e = 0; e < L.Length; e++)
    {
        int i = E[e * 2 + 0];
        int j = E[e * 2 + 1];
        float x = Vector3.Distance(vertices[i], vertices[j]);
        Vector3 f = L[e] * (vertices[i] - vertices[j]) / x;

        sum_X[i] += 0.5f * (vertices[i] + vertices[j] + f);
        sum_N[i]++;
        sum_X[j] += 0.5f * (vertices[i] + vertices[j] - f);
        sum_N[j]++;
    }

    for(int i = 0; i < V.Length; i++)
    {
        if(i == 0 || i == 20)	continue;
        Vector3 f = (0.2f * vertices[i] + sum_X[i]) / (0.2f + sum_N[i]);
        V[i] += 1 / t * (f - vertices[i]);
        vertices[i] = f;
    }

    mesh.vertices = vertices;
}

最后效果如下：

如果你觉得我的文章有帮助，欢迎关注我的微信公众号 我是真的想做游戏啊