作业三的内容主要就是实现一下FVM。我们按照文档中的步骤，第一步就是去独立地更新mesh的速度和位置，在初始化每个顶点的受力时，需要考虑到重力的影响。

for(int i=0 ;i<number; i++)
{//TODO: Add gravity to Force.Force[i] = new Vector3(0, -9.8f * mass, 0);
}for(int i=0; i<number; i++)
{//TODO: Update X and V here.V[i] += Force[i] * dt / mass;V[i] *= damp;X[i] += V[i]*dt;
}

接下来要考虑mesh的碰撞处理。在作业给的场景中，就只有一个简单的floor。

floor的世界坐标是(0, -3, 0)，那么这里就直接简化碰撞处理了：

//TODO: (Particle) collision with floor.
if (X[i].y < -3f)
{V[i].y += (-3f - X[i].y) / dt;X[i].y = -3f;
}

下一步，我们要在Start方法中预先计算好常量矩阵${\text{D}}_m$以及它的逆。根据PPT的第25页，有
$${\text{D}}_m=[{\text{X}}_{10}{\text{X}}_{20}{\text{X}}_{30}]$$
每个四面体都有这样的一个矩阵：

//TODO: Need to allocate and assign inv_Dm
inv_Dm = new Matrix4x4[tet_number];
for(int i = 0; i < tet_number; i++)
{Matrix4x4 Dm = Build_Edge_Matrix(i);inv_Dm[i] = Dm.inverse;
}Matrix4x4 Build_Edge_Matrix(int tet)
{Matrix4x4 ret=Matrix4x4.zero;//TODO: Need to build edge matrix here.Vector4 e1 = X[Tet[tet * 4 + 1]] - X[Tet[tet * 4 + 0]];Vector4 e2 = X[Tet[tet * 4 + 2]] - X[Tet[tet * 4 + 0]];Vector4 e3 = X[Tet[tet * 4 + 3]] - X[Tet[tet * 4 + 0]];ret.SetColumn(0, e1);ret.SetColumn(1, e2);ret.SetColumn(2, e3);ret.SetColumn(3, new Vector4(0, 0, 0, 1));return ret;
}

接下来，计算Deformation gradient F和Green strain G，计算的方法也在PPT第25页：
$$\text{F}=[{\text{x}}_{10}{\text{x}}_{20}{\text{x}}_{30}]{\text{D}}_m^{-1}$$

$$\text{G}=\frac{1}{2}({\text{F}}^T\text{F}-\text{I})$$

//TODO: Deformation Gradient
Matrix4x4 m = Build_Edge_Matrix(tet);
Matrix4x4 F = m * inv_Dm[tet];
//TODO: Green Strain
Matrix4x4 G = Matrix_Mul(0.5f, Matrix_Sub(F.transpose * F, Matrix4x4.identity));

作业里使用的是the second Piola-Kirchho stress，也给出了相应的公式：
$$\text{S}=2s_1\text{G}+s_{0}tr(\text{G})\text{I}$$

//TODO: Second PK Stress
Matrix4x4 S = Matrix_Add(Matrix_Mul(2.0f * stiffness_1, G),Matrix_Mul(stiffness_0 * Matrix_Trace(G), Matrix4x4.identity));

那么根据ppt，就可以计算出四面体上每个顶点的力了：
$$\text{P}=\text{F}\text{S}$$

$$[{\text{f}}_1{\text{f}}_2{\text{f}}_3]=-\frac{1}{6det({\text{D}}_m^{-1})}\text{P}{\text{D}}_m^{-\text{T}}$$

$${\text{f}}_0=-{\text{f}}_1-{\text{f}}_2-{\text{f}}_3$$

//TODO: Elastic Force
Matrix4x4 fm = Matrix_Mul(-1.0f / (6.0f * inv_Dm[tet].determinant),F * S * inv_Dm[tet].transpose);
f1 = fm.GetColumn(0);
f2 = fm.GetColumn(1);
f3 = fm.GetColumn(2);Force[Tet[tet * 4 + 0]] -= f1 + f2 + f3;
Force[Tet[tet * 4 + 1]] += f1;
Force[Tet[tet * 4 + 2]] += f2;
Force[Tet[tet * 4 + 3]] += f3;

其实这个时候房子已经可以跳动起来了，只是。。。

模拟很不稳定，房子直接飞了。作业中为了避免这一情况，使用了拉普拉斯平滑。我们需要统计每个顶点所邻接的所有其他顶点，计算出一个邻接的平均速度，然后进行加权：

void Laplacian_Smoothing(float blend)
{for (int i = 0; i < number; i++){V_sum[i] = Vector3.zero;V_num[i] = 0;}for(int tet=0; tet<tet_number; tet++){Vector3 sum = V[Tet[tet * 4 + 0]] + V[Tet[tet * 4 + 1]] + V[Tet[tet * 4 + 2]] + V[Tet[tet * 4 + 3]];V_sum[Tet[tet * 4 + 0]] += sum;V_sum[Tet[tet * 4 + 1]] += sum;V_sum[Tet[tet * 4 + 2]] += sum;V_sum[Tet[tet * 4 + 3]] += sum;V_num[Tet[tet * 4 + 0]] += 4;V_num[Tet[tet * 4 + 1]] += 4;V_num[Tet[tet * 4 + 2]] += 4;V_num[Tet[tet * 4 + 3]] += 4;}for(int i = 0; i < number; i++){V[i] = (V[i] + blend * V_sum[i] / V_num[i]) / (1 + blend);}
}

我们取blend为0.1，再看看现在的效果：

基础部分就完成了。下面再看一下附加题部分，其实就是根据PPT的第35页，换一种方法计算P。

首先是利用SVD得到三个矩阵：

Matrix4x4 U = Matrix4x4.zero;
Matrix4x4 S = Matrix4x4.zero;
Matrix4x4 V = Matrix4x4.zero;svd.svd(F, ref U, ref S, ref V);

根据公式：
$$P=\text{U}diag(\frac{\partial W}{\partial {\lambda }_0},\frac{\partial W}{\partial {\lambda }_1},\frac{\partial W}{\partial {\lambda }_2}){\text{V}}^{\text{T}}$$
那么问题就只剩计算这三个偏导数。而W是关于$I_C$和$I{I}_C$的函数（注意，这里PPT给的公式有误）
$$W=\frac{s_0}{8}(I_C-3{)}^2+\frac{s_1}{4}(I{I}_C-2I_C+3)$$
而
$$I_C={\lambda }_0^2+{\lambda }_1^2+{\lambda }_2^2$$

$$I{I}_C={\lambda }_0^4+{\lambda }_1^4+{\lambda }_2^4$$

又有
$$\frac{\partial W}{\partial \lambda }=\frac{\partial W}{\partial I_C}\frac{\partial I_C}{\partial \lambda }+\frac{\partial W}{\partial I{I}_C}\frac{\partial I{I}_C}{\partial \lambda }$$
所以只要分别计算上述式子中右侧的偏导数，就能得到最终的P了：

float lambda0 = S[0, 0];
float lambda1 = S[1, 1];
float lambda2 = S[2, 2];float Ic = lambda0 * lambda0 + lambda1 * lambda1 + lambda2 * lambda2;float dWdIc = 0.25f * stiffness_0 * (Ic - 3f) - 0.5f * stiffness_1;
float dWdIIc = 0.25f * stiffness_1;
float dIcdlambda0 = 2f * lambda0;
float dIcdlambda1 = 2f * lambda1;
float dIcdlambda2 = 2f * lambda2;
float dIIcdlambda0 = 4f * lambda0 * lambda0 * lambda0;
float dIIcdlambda1 = 4f * lambda1 * lambda1 * lambda1;
float dIIcdlambda2 = 4f * lambda2 * lambda2 * lambda2;
float dWd0 = dWdIc * dIcdlambda0 + dWdIIc * dIIcdlambda0;
float dWd1 = dWdIc * dIcdlambda1 + dWdIIc * dIIcdlambda1;
float dWd2 = dWdIc * dIcdlambda2 + dWdIIc * dIIcdlambda2;Matrix4x4 diag = Matrix4x4.zero;
diag[0, 0] = dWd0;
diag[1, 1] = dWd1;
diag[2, 2] = dWd2;
diag[3, 3] = 1;
Matrix4x4 P = U * diag * V.transpose;

最终效果如下，感觉两种计算方式的效果差不多：