CSC417-Physics Based Animation-Lab-3: Vanilla FEM Simulation

为什么要开这个坑:

其实我原本可以选择$\text{GAMES103}$的。。。但是一方面$\text{GAMES103}$用的$\text{Unity}$引擎做作业，我不太适应。另一方面，比较下来感觉$\text{CSC417}$对物理原理和数学推导的阐述比$\text{GAMES103}$更清楚（至少对我来说如此）。因此用这门课的作业来学习物理模拟。

本篇$\text{blog}$直接从$\text{Lab-3}$开始，因为前面的弹簧质点说实话。。。不难，就是麻烦其实是我第二个实验不知道为什么怎么写都炸，当然比$\text{GAMES101}$那个绳子模拟要难些，但是总体而言新东西不多，除了从拉格朗日视角推出了模拟的微分方程。

这次$\text{Lab}$重点在于引入了有限元这一重要方法，推导过程也较长较复杂。

$\text{0. 基本原理前置}$

这方面前置需要了解的东西是最小作用量原理/变分法，拉格朗日量以及$\text{Euler-Lagrange}$方程。这些在数学物理方法以及课程里都有讲，不再赘述。

当然，我在那些地方涉及到这个的时候也不会赘述

$\text{1. 有限元法推导}$

在这次$\text{Lab}$中，我们面对的是四面体网格，下面我们从头推导出完整的算法。

$\text{1.1. 计算线性插值形状函数}$

我们以四面体为考虑的单元，拿到的信息也是四面体顶点的信息，但是为了后面对整个体积求积分，我们必须先通过四面体顶点插值出四面体内部的点的性质。在这里我们采用简单的线性插值，系数也就是形状函数。

下面的推导全部默认向量为列向量。

我们设一个四面体在未形变时（也就是物质空间）四个顶点的坐标分别是$X_{1...4}$，在物质空间下四面体内部一点$X$可以写成四个顶点的凸组合（下面使用了爱因斯坦求和约定）:

$$X=\phi_i(X)X_i$$

并且满足$\sum_i \phi_i(X)=1$

因此可以列出如下线性方程组：

$$\begin{bmatrix}X_1 & X_2 & X_3 & X_4\\ 1&1&1&1\end{bmatrix}\begin{bmatrix}\phi_1(X)\\ \phi_2(X)\\ \phi_3(X) \\ \phi_4(X)\end{bmatrix}=\begin{bmatrix}X\\1\end{bmatrix}$$

即可解出所有的$\phi_i(X)$，我们不妨设$N(X)=\begin{bmatrix}\phi_1(X)\\ \phi_2(X)\\ \phi_3(X) \\ \phi_4(X)\end{bmatrix}$

或者我们作另一种表示，用来方便后续的计算，只采用$3$个变元$\phi$

$$\begin{bmatrix}X_2 - X_1&X_3-X_1&X_4-X_1 \end{bmatrix}\begin{bmatrix}\phi_2(X)\\\phi_3(X)\\\phi_4(X)\end{bmatrix}=X-X_1$$

而$\phi_1$可以由凸组合得到。这个写法的好处是便于求逆，也方便后续的求导计算。

$\text{1.2. 计算四面体动能}$

实际上我们的四面体是会动起来的，我们当然要用运动的观点看问题考虑形变以后的情况。不妨设物质空间下坐标为$X$的点现在在$x(X)$，四个顶点此时的坐标为$x_1,x_2,x_3,x_4$，根据我们使用的线性插值，有$x = \phi_i(X)x_i$。

接下来我们计算任意时刻四面体的动能，四面体的形状可能变了，但是质量没有变，我们考虑对原始的四面体来分析。

对于原始四面体，我们的假设是其密度，为常数$\rho$。考虑物质空间下的一个小体元$\text{d}\Omega$，其质量为$\rho\text{d}\Omega$，于是可以写出其动能为

$$\text{d}T=\frac{1}{2}\rho v^Tv\text{d}\Omega$$

整个四面体动能就是

$$T=\frac{1}{2}\rho \int_\Omega v^Tv\text{d}\Omega$$

这里$\Omega$是原始的四面体。

已知四个顶点的速度$v_{1...4}$，我们使用线性插值计算四面体中每个点的速度，于是有

$$v(X)=\begin{bmatrix}v_1 & v_2 & v_3 & v_4\end{bmatrix}\begin{bmatrix}\phi_1(X)\\ \phi_2(X)\\ \phi_3(X) \\ \phi_4(X)\end{bmatrix}=N^T(X)V$$

其中$N = \begin{bmatrix}\phi_1(X)\\ \phi_2(X)\\ \phi_3(X) \\ \phi_4(X)\end{bmatrix}$，于是代入前面的积分就得到

$$T=\frac{1}{2}\rho \int_\Omega V^TN(X) N^T(X)V\text{d}\Omega=\frac{1}{2}\rho V^T\int_\Omega N(X) N^T(X)\text{d}\Omega V$$

我们能把$V$提出来，是因为$V$是个跟四面体内部点无关的量（$V$只和四个顶点的速度相关）。

$\rho$是参数，所以实际上我们只需要算$\int_\Omega N(X) N^T(X)\text{d}\Omega$即可，这就是爆算而已：因为已知四个顶点的坐标，根据点的坐标$X$即可算出$N(X)$，进而可以算出积分，这种事情属于体力劳动，一般用符号计算软件解决。

$\text{1.3. 计算四面体势能}$

为了计算势能，我们引入形变梯度（可以参考连续介质力学笔记）以及应变能密度，首先考虑如何计算形变梯度$F$。

回忆$F$的定义：

$$F=\frac{\partial x}{\partial X}$$

也就是说我们要计算$\frac{\partial x}{\partial X}$,尴尬的是每次到这里我都会短路傻掉忘记怎么推。

我们已知

$$x(X) = x_1 + [-\phi_2(X)-\phi_3(X)-\phi_4(X)]x_1+\phi_2(X)x_2 + \phi_3(X)x_3+ \phi_4(X)x_4$$

$$=x_1 + \begin{bmatrix}x_1&x_2&x_3&x_4\end{bmatrix}\begin{bmatrix}-\phi_2-\phi_3-\phi_4\\\phi_2\\\phi_3\\\phi_4\end{bmatrix}$$

$$=x_1+ \begin{bmatrix}x_1&x_2&x_3&x_4\end{bmatrix}\begin{bmatrix}-1 -1 -1\\I\end{bmatrix}\begin{bmatrix}\phi_2\\\phi_3\\\phi_4\end{bmatrix}$$

$$=x_1+ \begin{bmatrix}x_1&x_2&x_3&x_4\end{bmatrix}\begin{bmatrix}-1 -1 -1\\I\end{bmatrix}\begin{bmatrix}X_2 - X_1&X_3-X_1&X_4-X_1 \end{bmatrix}^{-1}(X-X_1)$$

这样求导就很方便，直接可得

$$F= \begin{bmatrix}x_1&x_2&x_3&x_4\end{bmatrix}\begin{bmatrix}-1 -1 -1\\I\end{bmatrix}\begin{bmatrix}X_2 - X_1&X_3-X_1&X_4-X_1 \end{bmatrix}^{-1}$$

如果用四个$\phi$变元这里会很麻烦。这里好算的原因个人认为是由于自由度是$3$，用$3$个变量就刚好能算而且不多余。

我们可以发现这个形变梯度有个很好的性质，它和$X$无关，而且是线性映射，所以我们设 $D=\begin{bmatrix}-1 -1 -1\\I\end{bmatrix}\begin{bmatrix}X_2 - X_1&X_3-X_1&X_4-X_1 \end{bmatrix}^{-1}$，在任何时刻，用$D$乘以$x_{1...4}$立刻就可以得到形变梯度，这个矩阵在后面十分有用。

对于应变能密度，我们可以自己推导，也可以选择已有的模型。当然，这并不是这个作业的重难点，因此选用对弹性体模拟非常适合的$\text{neo-hookean}$模型。于是有

$$V=\int_\Omega\psi(F)\text{d}\Omega$$

$F$我们前面算出来了，而且很幸运是个和$X$无关的量，也就是说形变张量在四面体内处处相同。所以我们要计算这个积分值，直接用四面体体积乘上应变能密度即可。

$$V= \psi(F)\times \text{vol}(\Omega)$$

$\text{1.4. 矩阵展开为向量}$

我们计算了单个四面体的这些那些，现在可以把他们组装起来成为一个四面体网格系统。

能量是可加的，因此可以计算总动能

$$T=\frac{1}{2}\rho \sum\limits_{i=1}^nV_i^T\int_{\Omega_i} N_i(X) N_i^T(X)\text{d}\Omega V_i$$

到这里事情不太好做了，因为我们要处理一堆矩阵的乘法，如果用选择矩阵来表示$V$，由于$V$的列数和广义速度$\dot{q}$的列数不一致，我们就不能这么做。实际上更麻烦的地方还在后面，等到对矩阵求导的时候求出来一堆高阶张量就不好收拾了。

这里的核心在于我们重新排列矩阵的元素，或者说将矩阵展平成向量。其实对于质点弹簧系统我们也做过这样的事情，只不过那个时候每个考虑的单元是$2$个顶点，这里是$4$个。

考虑重新定义$V$和$N$:

$$V=\begin{bmatrix}v_1\\v_2\\v_3\\v_4\end{bmatrix},N(X)=\begin{bmatrix}\phi_1(X)I&\phi_2(X)I&\phi_3(X)I&\phi_4(X)I\end{bmatrix}$$

同样的，我们表示四个顶点的$x$的时候，也把它们展成一个向量。我们这么做，前面的式子还是满足的，只不过有些矩阵的形状和内容要相应地修改一下。

比如说求微分用的矩阵$D$，原本$D$是需要乘上$x_{1...4}$拼成的矩阵，现在需要乘上$x_{1...4}$拼成的大向量，当然，我们输出的结果也要是一个大向量。

这个倒没啥难的，把$D$矩阵的元素合理分配一下就行。我们的$F$也展开成一个$9$维向量，因此我们实际需要的矩阵是一个$9\times 12$的矩阵

$$\frac{\partial x}{\partial X}= \begin{bmatrix} D_{11}&0&0&D_{21}&0&0&D_{31}&0&0&D_{41}&0&0\\ D_{12}&0&0&D_{22}&0&0&D_{32}&0&0&D_{42}&0&0\\ D_{13}&0&0&D_{23}&0&0&D_{33}&0&0&D_{43}&0&0\\ 0&D_{11}&0&0&D_{21}&0&0&D_{31}&0&0&D_{41}&0\\ 0&D_{12}&0&0&D_{22}&0&0&D_{32}&0&0&D_{42}&0\\ 0&D_{13}&0&0&D_{23}&0&0&D_{33}&0&0&D_{43}&0\\ 0&0&D_{11}&0&0&D_{21}&0&0&D_{31}&0&0&D_{41}\\ 0&0&D_{12}&0&0&D_{22}&0&0&D_{32}&0&0&D_{42}\\ 0&0&D_{13}&0&0&D_{23}&0&0&D_{33}&0&0&D_{43}\\ \end{bmatrix} \begin{bmatrix}x_1\\x_2\\x_3\\x_4\end{bmatrix}$$

这玩意儿写出来感觉也没啥技巧。。。暴力展开也没啥问题，毕竟我们已经把$x_{1...4}$组织成列向量了，去考虑矩阵每一个行向量和这个列向量的点积会比较方便。

为了方便后续的使用我们不妨称这个矩阵为$B_i$，这时来组装就很方便。

$\text{1.5. 组装}$

前面得到了式子：

$$T=\frac{1}{2}\rho \sum\limits_{i=1}^nV_i^T\int_{\Omega_i} N_i(X) N_i^T(X)\text{d}\Omega V_i$$

按照套路，现在再来用选择矩阵的技术处理，首先记住$V$是个$12$维列向量。

我们的选择矩阵是一个$12\times 3n$的矩阵，它可以看成从上往下$4$个$3\times 3n$的子矩阵拼成，每一个大概长成这个样子：

$$\begin{bmatrix}O&O&O&...&I&O&...&O\\\end{bmatrix}$$

也就是说，除了该四面体的四个顶点对应的位置是三阶单位阵，其他的位置都是零矩阵。

然后事情就好说了，$V=S_i\dot{q}$，代入前面的求和式

$$T=\frac{1}{2}\rho\dot{q}^T \sum\limits_{i=1}^nS_i^T\int_{\Omega_i} N_i(X) N_i^T(X)\text{d}\Omega S_i\dot{q}$$

我们成功把$\dot{q}$从里提取了出来，中间求和式里有了选择矩阵，意义就十分明确了，将积分得到的矩阵组装成一整个大矩阵就行，按照定义，组装成的大矩阵就是质量矩阵$M$。

也就是说：

$$M=\sum\limits_{i=1}^nS_i^T\int_{\Omega_i} N_i(X) N_i^T(X)\text{d}\Omega S_i$$

计算$T$正是为了得到这个$M$。

然后是势能：

$$V=\sum\limits_{i=1}^n \psi(F_i) \times \text{vol}(\Omega_i)$$

为了数值求解微分方程，我们现在需要做的事情是求出势能对$q$的全微分（也就是负的广义力）。这里实际上会随$q$变化而变化的量只有$F$（因为$\text{vol}(\Omega_i)$是未形变时的体积），于是

$$\frac{\partial V}{\partial q}=\sum\limits_{i=1}^n\text{vol}(\Omega_i)\frac{\partial \psi(F_i)}{\partial F_i} \frac{\partial F_i}{\partial q}$$

$$=\sum\limits_{i=1}^n\text{vol}(\Omega_i)\frac{\partial \psi(F_i)}{\partial F_i} \frac{\partial F_i}{\partial q}$$

$\frac{\partial \psi(F_i)}{\partial F_i}$我们可以通过符号计算软件求出公式，因此我们实际上只用化简$\frac{\partial F_i}{\partial q}$

$$\frac{\partial F_i}{\partial q}=\frac{\partial (B_iS_iq)}{\partial q}=S_i^TB_i^T$$

所以

$$\frac{\partial V}{\partial q} = \sum\limits_{i=1}^n\text{vol}(\Omega_i)S_i^TB_i^T\frac{\partial \psi(F_i)}{\partial F_i}$$

可以看到$S_i$又出现了，所以这还是个组装过程。

我们还需要求出刚度矩阵，也就是$V$对$q$的二阶导

$$\frac{\partial^2 V}{\partial q^2}=\sum\limits_{i=1}^n\text{vol}(\Omega_i)S_i^TB_i^T\frac{\partial^2 \psi(F_i)}{\partial F_i^2}\frac{\partial F_i}{\partial q}$$

$$=\sum\limits_{i=1}^n\text{vol}(\Omega_i)S_i^TB_i^T\frac{\partial^2 \psi(F_i)}{\partial F_i^2}B_iS_i$$

这样，必须的量我们就推——完——了！

$\text{2. 解微分方程}$

$\text{2.1. 转化为优化问题}$

使用有限元，我们把该求的量都求出来了，现在该考虑考虑怎么解微分方程了。

微分方程是

$$M\ddot{q}=-\frac{\partial V}{\partial q}$$

当然，为了方便，我们可以记$f(q)=-\frac{\partial V}{\partial q}$

我们用白金之星隐式欧拉法来数值求解这个方程的话，写出来式子就是

$$M\dot{q}_{t+1}=M\dot{q}_t+f(q_{t+1})\Delta t$$

$$q_{t+1} =q_t + \dot{q}_{t+1}\Delta t$$

把下式代入上式，有

$$M\dot{q}_{t+1}=M\dot{q}_t+f(q_t+\dot{q}_{t+1}\Delta t)\Delta t$$

这是关于$\dot{q}_{t+1}$ 的一元方程，只需要求出它的根就好了，为了好写，设 $v=\dot{q}_{t+1}$。

这个方程不可以直接解算，下面我们将使用一种重要的手法来求它的数值解。

方程求零点问题往往可以转化为最优化问题，变成最优化问题我们就有很多好用的方法去求解，所以我们考虑对这个方程积分，来得到一个原函数，这样我们求出原函数的极值点，就相当于求出了方程的零点。

$$Mv-M\dot{q}_t-f(q_t+v\Delta t)\Delta t =0$$

对$v$积分得到

$$\frac{1}{2}(v-\dot{q}_t)M(v-\dot{q}_t)+V(q_t+v\Delta t)$$

这里其实做了一步近似，就是把$f$外面的$\Delta t$看成$\frac{\partial q}{\partial v}$然后运用链式法则，倒也没啥毛病。

所以说现在我们需要求的就是

$$v^{\ast} = \argmin\limits_v [\frac{1}{2}(v-\dot{q}_t)M(v-\dot{q}_t)+V(q_t+v\Delta t)]$$

不妨记上面这个优化目标为$E$。

$\text{2.2. Newton's Method}$

如果不能够数值求解，那么几乎肯定只能慢慢迭代，我们首先选定一个方向$d$，然后沿着这个方向走一小步，不断循环直到我们发现结果几乎不再变化就可以停止了，我们就求出了近似的数值解。

如何选择这个方向呢？我们可以考虑优化目标下降最快的方向，也就是

$$d = \argmin\limits_{\tilde{d}} E(v_i + \tilde{d})$$

这个意思很明确，解释一下记号，假设我们现在处于第$i+1$轮迭代，第$i$轮迭代后我们处于$v_i$。

但是实际上我们这是把自己绕进了另一个坑。。。另一个不好直接求解的优化问题，这就只能靠近似了，进行二阶展开。

$$d = \argmin [E(v_i) + g_i^T\tilde{d} + \frac{1}{2}\tilde{d}^TH_i\tilde{d}]=\argmin [g_i^T\tilde{d} + \frac{1}{2}\tilde{d}^TH_i\tilde{d}]$$

这个要优化起来总算容易些了，先不考虑怎么算梯度和算$\text{Hessian}$矩阵，求个$d$总归是没问题的。

$$d = -H_i^{-1}g_i$$

这里其实已经能看出一些牛顿迭代法的味道了。

那么一个问题就是$H_i$和$g_i$应该怎么算，回想一下$E$是我们积分求出来的。。。所以$g_i$自然就是原来那个东西，也就是

$$g_i = M(v_i-\dot{q}_i)-f(q_t+v_i\Delta t)\Delta t$$

进一步求导得$\text{Hessian}$矩阵

$$H_i = M - \frac{\partial f}{\partial q}(q_t+v_i\Delta t)\Delta t^2=M+\frac{\partial^2 V}{\partial q^2}(q_t+v_i\Delta t)\Delta t^2$$

下一步是选择一个合适的步长$\alpha$，然后更新$v_{i+1}=v_i+\alpha d$，选择这个$\alpha$我们还是可以变成一个优化问题

我们现在的目标是

$$\alpha=\argmin\limits_{\alpha^\ast} E(v_i+\alpha^\ast d)$$

这可以用一点简单的线搜索方法。

1. 选择一个初始猜测值$\alpha_{max}$

2. 如果$E(v_i+\alpha d)\le E(v_i)+cg_i^Td$或者$\alpha$已经足够小

3. $\alpha = p\alpha$

当然，这里$p, c$以及我们认为$\alpha$“足够小”都是参数也就是说最终还是逃不过调参的命运。

$\text{3. 用稀疏网格控制密集网格}$

最后一步，如果我们的模型网格非常细密，我们不可能继续用上面的方法，这样我们会遇到巨量的四面体，算不了。

所幸，四面体网格的方法在这方面体现了优势。我们可以对一个相对稀疏的网格进行物理模拟，然后运用线性插值的方法，把细密的网格运动给插值出来。

自然而然，我们会发现这个插值系数就是第一步计算出来的线性插值形状函数，所以我们只需要把这些所有的形状函数拼成一个超级巨大的稀疏矩阵，乘以我们计算出来的稀疏网格的广义坐标即可。

$\text{4. Coding}$

理解了上面的推导，写代码就几乎没什么实现上的难度其实上一个实验我也这么自信，但是不知道为何怎么写怎么炸。

真的有些麻烦的主要是那些需要符号计算软件进行求微分/积分的部分不过我懒得用符号计算软件所以直接借鉴别人的。