GAMES103课程笔记05-Physics-Based Cloth Simulation
这个系列是GAMES103：基于物理的计算机动画入门(GAMES 103: Intro to Physics Based Animation)的同步课程笔记。本课程是计算机动画的入门课程，着重介绍各种基于物理的动画仿真模拟技术。本节课主要介绍弹簧质点系统与布料模拟。
Spring Mass System
Ideal Spring
在正式介绍弹簧质点系统(spring mass system)前我们首先来回忆下弹簧的力学性质。根据胡克定律一维理想弹簧的能量和弹力可以定义为：
\[E(x) = \frac{1}{2} k (x - L)^2\] \[f(x) = -\frac{d E}{d x} = -k (x - L)\]其中\(k\)为弹簧的刚度，\(L\)为弹簧的原长。
而对于三维空间中的弹簧系统的内力，我们同样可以通过对系统能量求导来计算：
\[E (\mathbf{x}) = \frac{k}{2} (\Vert \mathbf{x}_i - \mathbf{x}_j \Vert - L)^2\] \[\mathbf{f}_i (\mathbf{x}) = - \nabla_i E = -k (\Vert \mathbf{x}_i - \mathbf{x}_j \Vert - L) \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert}\] \[\mathbf{f}_j (\mathbf{x}) = - \nabla_j E = -k (\Vert \mathbf{x}_j - \mathbf{x}_i \Vert - L) \frac{\mathbf{x}_j - \mathbf{x}_i}{\Vert \mathbf{x}_j - \mathbf{x}_i \Vert}\]
而当系统包含多个弹簧时，系统总能量等于单个弹簧的能量之和，每个点受到的弹力等于与它相连的所有弹簧内力之和：
\[E = \sum_j \frac{k}{2} (\Vert \mathbf{x}_i - \mathbf{x}_j \Vert - L_j)^2\] \[\mathbf{f}_i (\mathbf{x}) = - \nabla_i E = \sum_j -k (\Vert \mathbf{x}_i - \mathbf{x}_j \Vert - L_j) \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert}\]
Representation
使用弹簧质点系统进行布料仿真时我们需要将布料转化为弹簧构成的网络。除了要在规则的节点之间添加水平和竖直方向的弹簧外我们还需要对角线方向(diagonal)以及弯曲(bending)弹簧来抵抗布料在不同方向的拉伸和翻折，一些结构化的弹簧网络可参考下图：
除了结构化弹簧网格外我们也可以直接使用模型的三角网格来作为弹簧系统。在这种情况下弹簧网络往往是非结构化的，同时还要在网格的内部边上添加对角的弯曲弹簧以抵抗弯曲：
接下来的问题是如何来表示和存储弹簧网络。最直接的想法是使用列表来存储模型的顶点和三角形，然后每个三角形的三条边分别对应三根弹簧。但这种做法的缺陷在于三角网格内部的边会被重复计算，而且没有办法直接表示三角网格内的弯曲弹簧。
因此我们需要对三角形的列表进行一些处理：首先将每个三角形用一个三元组来表示，其中前两位表示三角形的两个顶点最后一位为三角形的编号；然后我们对所有的三元组进行排序，这样同一条边对应的不同三角形会排在相邻的位置上；最后把重复的边删除掉就可以得到网格上的所有弹簧。而在删除重复边的过程中我们已经得到了网格内部的边，我们可以利用这些重复的边来构造内部的弯曲弹簧。
Explicit Integration
有了弹簧网络后我们就可以利用质点动力学的方法来实现布料的仿真，整个系统的计算过程如下：
按照上式更新质点速度和位置的方法称为显式积分(explicit integration)。显式积分的主要缺陷在于在积分过程中系统的能量可能会发散，这种现象也称为数值不稳定性(numerical instability)。要提高系统的稳定性通常可以考虑缩减时间步长\(\Delta t\)，但这样做会降低仿真的效率。因此现代计算机模拟中一般不会使用显式积分的方式。
Implicit Integration
和显式积分相比隐式积分(implicit integration)具有更好的数值稳定性。在Rigid Body Dynamics一节中我们已经介绍过隐式积分的基本概念，对应的积分公式为：
\[\mathbf{v}^{[1]} = \mathbf{v}^{[0]} + \Delta t \mathbf{M}^{-1} \mathbf{f}^{[1]}\] \[\mathbf{x}^{[1]} = \mathbf{x}^{[0]} + \Delta t \mathbf{v}^{[1]}\]通过交换计算顺序，上式也可以等价地表示为：
\[\mathbf{x}^{[1]} = \mathbf{x}^{[0]} + \Delta t \mathbf{v}^{[0]} + \Delta t^2 \mathbf{M}^{-1} \mathbf{f}^{[1]}\] \[\mathbf{v}^{[1]} = \frac{\mathbf{x}^{[1]} - \mathbf{x}^{[0]}}{\Delta t}\]Newton-Raphson Method
假设外力\(\mathbf{f}\)只与质点的位置有关，则隐式积分等价于求解方程：
\[\mathbf{x}^{[1]} = \mathbf{x}^{[0]} + \Delta t \mathbf{v}^{[0]} + \Delta t^2 \mathbf{M}^{-1} \mathbf{f} (\mathbf{x}^{[1]})\]进一步还可以证明它还等价于求解如下定义的优化问题：
\[\mathbf{x}^{[1]} = \arg \min \frac{1}{2 \Delta t^2} \Vert \mathbf{x} - \mathbf{x}^{[0]} - \Delta t \mathbf{v}^{[0]} \Vert_\mathbf{M} + E(\mathbf{x})\] \[\Vert \mathbf{x} \Vert_\mathbf{M} = \mathbf{x}^T \mathbf{M} \mathbf{x}\]证明方法很简单，只需要对优化函数求导并让导数取0即可：
\[\begin{aligned} \frac{d F}{d \mathbf{x}} &= \frac{d}{d \mathbf{x}} \frac{1}{2 \Delta t^2} \Vert \mathbf{x} - \mathbf{x}^{[0]} - \Delta t \mathbf{v}^{[0]} \Vert_\mathbf{M} + E(\mathbf{x}) \\ &= \frac{1}{\Delta t^2} \mathbf{M} (\mathbf{x} - \mathbf{x}^{[0]} - \Delta t \mathbf{v}^{[0]}) - \mathbf{f} (\mathbf{x}) \\ &= \mathbf{0} \end{aligned}\] \[\mathbf{x} - \mathbf{x}^{[0]} - \Delta t \mathbf{v}^{[0]} - \Delta t^2 \mathbf{M}^{-1} \mathbf{f} (\mathbf{x}) = \mathbf{0}\] \[\mathbf{x} = \mathbf{x}^{[0]} + \Delta t \mathbf{v}^{[0]} + \Delta t^2 \mathbf{M}^{-1} \mathbf{f} (\mathbf{x})\]因此我们可以将隐式积分转换为求解一个非线性优化问题。求解非线性优化问题的基本方法是Newton-Raphson迭代(Newton-Raphson method)，也称为Newton’s method，它的基本思想是对优化目标的梯度进行线性化来计算增量\(\Delta \mathbf{x}\)，然后不断迭代直到收敛：
回到我们的弹簧系统，我们需要计算出\(F(\mathbf{x})\)的二阶导数：
\[\frac{d^2 F}{d \mathbf{x}^2} = \frac{1}{\Delta t^2} \mathbf{M} + \mathbf{H}\]其中\(\mathbf{H}\)为弹簧系统的Hessian矩阵：
\[\mathbf{H} (\mathbf{x}) = \sum_{e = \{ i, j \}} \begin{bmatrix} \frac{\partial^2 E}{\partial \mathbf{x}_i^2} & \frac{\partial^2 E}{\partial \mathbf{x}_i \partial \mathbf{x}_j} \\ \frac{\partial^2 E}{\partial \mathbf{x}_i \partial \mathbf{x}_j} & \frac{\partial^2 E}{\partial \mathbf{x}_j^2} \end{bmatrix} = \sum_{e = \{ i, j \}} \begin{bmatrix} \mathbf{H}_e & -\mathbf{H}_e \\ -\mathbf{H}_e & \mathbf{H}_e \end{bmatrix}\]其中，
\[\mathbf{H}_e = k \frac{\mathbf{x}_{ij} \mathbf{x}_{ij}^T}{\Vert \mathbf{x}_{ij} \Vert^2} + k \bigg(1 - \frac{L}{\Vert \mathbf{x}_{ij} \Vert} \bigg) \bigg(\mathbf{I} - \frac{\mathbf{x}_{ij} \mathbf{x}_{ij}^T}{\Vert \mathbf{x}_{ij} \Vert^2} \bigg)\] \[\mathbf{x}_{ij} = \mathbf{x}_i - \mathbf{x}_j\]我们把\(F(\mathbf{x})\)的梯度和Hessian矩阵带入迭代公式就得到了整个弹簧质点系统的状态更新方法：
这样我们只需要通过求解线性方程组就能够实现整个仿真过程。不过需要注意的是Newton-Raphson迭代不能保证一定能收敛到最小值，实际上它只能保证收敛到局部极值，要判断收敛到极大值还是极小值还需要借助\(F(\mathbf{x})\)的Hessian矩阵：当Hessian矩阵为正定矩阵时为极小值，负定时为极大值：
显然我们希望系统的Hessian矩阵永远是正定的，这样我们总能找到最小值点。那么整个优化函数的Hessian矩阵究竟如何呢？我们首先来看一眼\(\mathbf{H}_e\)矩阵的形式：\(\mathbf{H}_e\)矩阵可以看成两个矩阵的和，其中第一项是一个半正定矩阵而第二项在弹簧处于压缩的状态时是一个半负定的矩阵：
这种形式就导致整个Hessian矩阵在弹簧处于压缩时可能是非正定的：
当Hessian矩阵出现非正定的情况时整个系统可能存在多个极小值，从力学角度上来看这表示当弹簧出现压缩时可能会出现多个能量极小值位置，它们都是系统的可行解：
如果要保证Hessian矩阵在求解时是正定的，一般可以考虑减小时间步长\(\Delta t\)或者在出现压缩时丢掉\(\mathbf{H}_e\)矩阵的第二项。
Jacobi Method
最后我们来考虑如何求解线性方程组\(\mathbf{A} \Delta \mathbf{x} = \mathbf{b}\)，最基本的解法是Jacobi迭代(Jacobi method)：
除此之外也可以使用一些直接解法或是其它的迭代解法来进行计算：
Bending and Locking Issues
The Bending Spring Issue
使用弹簧质点系统进行仿真时的一个主要缺陷在于弯曲(bending)。假设有如下图所示的弹簧质点系统，其中蓝色线为弯曲弹簧。当\(\mathbf{x}_1\)和\(\mathbf{x}_2\)绕\(\mathbf{x}_0\)进行弯曲时，弯曲弹簧的长度不会发生明显变化，因此几乎不会产生抵抗弯曲的弹力。这种现象说明当弹簧质点系统位于一个平面上时，整个系统几乎无法抵抗出平面方向上的弯曲，而显然这与现实生活中的纸张、布料等材质是不符的。
因此我们需要额外添加抵抗出平面弯曲的内力。一种简单的做法是假设内力与两个三角形平面的夹角有关，两个面的法向夹角越大抵抗弯曲的内力就越大，这种方法称为二面角模型(dihedral angle model)。在二面角模型中每个二面角对应4个节点上的内力，这些内力都是二面角\(\theta\)的函数\(\mathbf{f}_i = f(\theta) \mathbf{u}_i\)，而且它们在几何上需要满足约束：
- \(\mathbf{u}_1\)和\(\mathbf{u}_2\)的方向需要与所在三角形的法向相同；
- \(\mathbf{u}_3\)和\(\mathbf{u}_4\)不会使弹簧伸长，因此\(\mathbf{u}_4 - \mathbf{u}_3\)需要垂直于所在边，即\((\mathbf{u}_4 - \mathbf{u}_3) \perp \mathbf{n}_1, \mathbf{n}_2\)；
- 系统的内力为\(\mathbf{0}\)，即\(\mathbf{u}_1 + \mathbf{u}_2 + \mathbf{u}_3 + \mathbf{u}_4 = \mathbf{0}\)，因此\(\mathbf{u}_3\)和\(\mathbf{u}_4\)必为\(\mathbf{n}_1\)和\(\mathbf{n}_2\)的线性组合。
求解这个系统得到每个顶点上内力的方向和大小如下：
二面角模型通过计算系统内力来实现仿真，比较适合显式积分来求解。相对地，二面角模型很难计算能量项以及导数，因此基本不会使用隐式积分来求解。
除了二面角模型外还可以使用二次弯曲模型(quadratic bending model)来处理。此时我们假设系统的初始状态为平面，而且在仿真过程中系统基本会发生拉伸，这样可以定义系统的能量为：
由于能量项具有二次的形式，我们可以方便地计算出系统内力和Hessian矩阵：
\[\mathbf{f} (\mathbf{x}) = -\nabla E(\mathbf{x}) = -\mathbf{Q} \begin{bmatrix} \mathbf{x}_0 \\ \mathbf{x}_1 \\ \mathbf{x}_3 \\ \mathbf{x}_4 \end{bmatrix}\] \[\mathbf{H} (\mathbf{x}) = \frac{\partial^2 E}{\partial \mathbf{x}^2} = \mathbf{Q}\]显然二次弯曲模型非常适合隐式积分，但需要注意当系统存在不能忽略的拉伸或是初始状态不是平面时需要对能量进行修正。
The Locking Issue
布料的力学性质之一是它具有很强的抗拉刚度但它的抗弯刚度却很小，在日常生活中我们可以轻松地弯曲衣物但对它进行拉伸却非常困难。使用弹簧质点系统来模拟布料的这种力学行为容易出现锁死(locking)的问题：假设有如下图所示的弹簧质点系统，当弹簧刚度很大时我们仍然可以沿左图虚线方向来进行翻折，但不能沿右图的虚线方向来翻折，即系统的弯曲被锁死了。
锁死现象的本质是我们的弹簧系统有过多的约束，这导致了进行仿真时系统的自由度不够。对于三角形网格(流形)，根据欧拉公式有：
\[\vert E \vert = 3 \times \vert V \vert - 3 - \vert E_\text{boundary} \vert\]使用弹簧质点系统进行仿真时的变量数为\(3 \vert V \vert\)，边构成的约束为\(\vert E \vert\)，因此整个系统的自由度为：
\[3 \vert V \vert - \vert E \vert = 3 + \vert E_\text{boundary} \vert\]这说明使用弹簧质点系统进行仿真时的自由度是不够的，因此在弹簧刚度比较大或是网格的分辨率比较低的情况下需要额外注意。