Fast Mass Spring 方法学习与布料仿真

符号约定

• 离散时间节点 $t_1,t_2,...,t_n,...$ 和固定时间间隔 $h=t_n-t_{n-1}$

• 系统节点数目 $m$ ，系统弹簧数目 $s$

• 系统状态 ${\mathbf{q}}_n\in {\mathbb{R}}^{3m}$，描述在 $t_n$ 时刻每一个节点的空间位置

• 质量矩阵 $\mathbf{M}\in {\mathbb{R}}^{3m\times 3m}$，是一个对角矩阵

• 速度 ${\mathbf{v}}_n\in {\mathbb{R}}^{3m}$，描述在 $t_n$ 时刻每一个节点的速度向量

• 外力 $\mathbf{f}({\mathbf{q}}_n):{\mathbb{R}}^{3m}\to {\mathbb{R}}^{3m}$

• 系统能量 $E:{\mathbb{R}}^{3m}\to \mathbb{R}$，有 $\mathbf{f}=-\nabla E$

公式推导

根据牛顿第二定律得到关于系统状态的隐式积分
$$\begin{array}{rl}{\mathbf{q}}_{n+1}& ={\mathbf{q}}_n+h{\mathbf{v}}_{n+1}\\ {\mathbf{v}}_{n+1}& ={\mathbf{v}}_n+h{\mathbf{M}}^{-1}\mathbf{f}({\mathbf{q}}_{n+1})\end{array}\text{\hspace{1em}(1)}$$
我们消去 $\mathbf{v}$ 可以得到
$${\mathbf{q}}_{n+1}-2{\mathbf{q}}_n+{\mathbf{q}}_{n-1}=h^2{\mathbf{M}}^{-1}\mathbf{f}({\mathbf{q}}_{n+1})\text{\hspace{1em}(3)}$$
物理仿真的目的是根据已知条件求解 ${\mathbf{q}}_{n+1}$，因此我们可以简记 $\mathbf{x}:={\mathbf{q}}_{n+1}$，$\mathbf{y}:=2{\mathbf{q}}_n-{\mathbf{q}}_{n-1}$，那么方程就变成了
$$\mathbf{M}(\mathbf{x}-\mathbf{y})=h^2\mathbf{f}(\mathbf{x})\text{\hspace{1em}(4)}$$
式中 $\mathbf{y}$ 可以看成是惯性，如果没有外力，那么下一刻系统状态将会变成 ${\mathbf{q}}_n+({\mathbf{q}}_n-{\mathbf{q}}_{n-1})$。要注意的是我们平时容易把 $\mathbf{y}$ 当成函数来看，但这里方程(4)中只有 $\mathbf{x}$ 和 $\mathbf{f}(\mathbf{x})$ 是未知的，其它都是已知的常量。

构造辅助函数
$$g(\mathbf{x})=\frac{1}{2}(\mathbf{x}-\mathbf{y}{)}^T\mathbf{M}(\mathbf{x}-\mathbf{y})+h^{2}E(\mathbf{x})\text{\hspace{1em}(5)}$$
求解(4)等价于求解(5)的极小值点。

典型的求解非线性函数极值点的方法是牛顿法：
$${\mathbf{x}}_{n+1}={\mathbf{x}}_n-Hes{s}_g({\mathbf{x}}_n{)}^{-1}\nabla g({\mathbf{x}}_n)$$
但是牛顿法每一次迭代中，海森矩阵都会变化，因此每次迭代都要求一次海森矩阵的逆，非常费时间。

接下来就是该算法的精华所在。

根据胡克定律，系统能量（弹性势能）表示为：
$$E(\mathbf{x})=\frac{1}{2}\sum _{i=1}^s{k}_i(||{\mathbf{p}}_{i_1}-{\mathbf{p}}_{i_2}||-r{)}^2\text{\hspace{1em}(6)}$$
其中 $k_i$ 是弹性系数，${\mathbf{p}}_{i_1}$ 和 ${\mathbf{p}}_{i_2}$ 是弹簧的端点坐标，i1,i2∈{1,2,...,m}i_1,i_2\in\{1,2,...,m\}i1​,i2​∈{1,2,...,m}，$r$ 是弹簧原长。但是这个式子是不好处理的，因此论文中引入辅助变量 ${\mathbf{d}}_i$，将弹性势能等价表示为：
$$E(\mathbf{x})=\underset{\mathbf{d}\in U}{\min }\frac{1}{2}\sum _{i=1}^s{k}_i||{\mathbf{p}}_{i_1}-{\mathbf{p}}_{i_2}-{\mathbf{d}}_i|{|}^2\text{\hspace{1em}(7)}$$
其中 U={(d1,...,ds)∈R3s:∣∣di∣∣=ri}U=\{(\boldsymbol{d}_1,...,\boldsymbol{d}_s)\in\mathbb{R}^{3s}:||\boldsymbol{d}_i|| = r_i \}U={(d1​,...,ds​)∈R3s:∣∣di​∣∣=ri​} （论文里这里有一点纰漏）。

(6)和(7)等价的原因是 $(||{\mathbf{p}}_1-{\mathbf{p}}_2||-r{)}^2=\underset{||\mathbf{d}||=r}{\min }||{\mathbf{p}}_1-{\mathbf{p}}_2-\mathbf{d}|{|}^2$ ，这个很好理解，一个点距离球面最近的点，就是这个点与球心连线（的延长线）和球面的交点。变换成 (7) 式之后，论文中直接华丽地把公式改写成了矩阵形式，我这里做一点点简单推导：
$$\begin{array}{rl}||{\mathbf{p}}_1-{\mathbf{p}}_2-\mathbf{d}|{|}^2& =||{\mathbf{p}}_1|{|}^2+||{\mathbf{p}}_2|{|}^2+||d|{|}^2-2{\mathbf{p}}_1\cdot {\mathbf{p}}_2-2{\mathbf{p}}_1\cdot \mathbf{d}+2{\mathbf{p}}_2\cdot \mathbf{d}\\ & =||{\mathbf{p}}_1-{\mathbf{p}}_2|{|}^2-2({\mathbf{p}}_1-{\mathbf{p}}_2)\cdot \mathbf{d}+r^2\end{array}$$
由于是优化问题，常数项 $r^2$ 可以丢弃，因此式中剩下一个很像二次函数的东西，于是就可以得到论文中的华丽变换（论文中这个公式其实不是严格等号，差了一个 $r^2$ 常数项，加上 min 会好一些）：
$$\underset{\mathbf{d}\in U}{\min }\frac{1}{2}\sum _{i=1}^s{k}_i||{\mathbf{p}}_{i_1}-{\mathbf{p}}_{i_2}-{\mathbf{d}}_i|{|}^2=\underset{\mathbf{d}\in U}{\min }\frac{1}{2}{\mathbf{x}}^T\mathbf{L}\mathbf{x}-{\mathbf{x}}^T\mathbf{J}\mathbf{d}\text{\hspace{1em}(8)}$$
其中 $\mathbf{x}=({\mathbf{p}}_1,...,{\mathbf{p}}_m)$，矩阵 $\mathbf{L}\in {\mathbb{R}}^{3m\times 3m},\mathbf{J}\in {\mathbb{R}}^{3m\times 3s}$ 定义如下：
$$\begin{gathered} \mathbf{L}=(\sum _{i=1}^s{k}_i{\mathbf{A}}_i{\mathbf{A}}_i^T)⨂{\mathbf{I}}_3 \\ \mathbf{J}=(\sum _{i=1}^s{k}_i{\mathbf{A}}_i{\mathbf{S}}_i^T)⨂{\mathbf{I}}_3\text{\hspace{1em}(9)} \end{gathered}$$
其中 ${\mathbf{A}}_i\in {\mathbb{R}}^m$ 是关于弹簧端点索引 $i_1,i_2$ 的向量，有 $A_{i,i_1}=1,A_{i,i2}=-1$，其它向量元素都是 0 ；而 ${\mathbf{S}}_i\in {\mathbb{R}}^s$ 是关于弹簧索引 $i$ 的向量，有 $S_{i,i}=1$，其它向量元素都是 0 ；$⨂$ 是克罗内克积。

如果算上外力 ${\mathbf{f}}_{ext}\in {\mathbb{R}}^{3m}$，系统能量可以表示为：
$$E(\mathbf{x})=\underset{\mathbf{d}\in U}{\min }\frac{1}{2}{\mathbf{x}}^T\mathbf{L}\mathbf{x}-{\mathbf{x}}^T\mathbf{J}\mathbf{d}-{\mathbf{x}}^T{\mathbf{f}}_{ext}\text{\hspace{1em}(10)}$$
将(10)代入(5)得到最终的优化问题
$$\underset{\mathbf{x}\in {\mathbb{R}}^{3m},\mathbf{d}\in U}{\min }\frac{1}{2}{\mathbf{x}}^T(\mathbf{M}+h^2\mathbf{L})\mathbf{x}-h^2{\mathbf{x}}^T\mathbf{J}\mathbf{d}-h^2{\mathbf{x}}^T{\mathbf{f}}_{ext}-\frac{1}{2}{\mathbf{x}}^T\mathbf{M}\mathbf{y}-\frac{1}{2}{\mathbf{y}}^T\mathbf{M}\mathbf{x}\text{\hspace{1em}(11)}$$
式中舍弃了一个常数项 $\frac{1}{2}{\mathbf{y}}^T\mathbf{M}\mathbf{y}$ 。

数值计算

论文中描述的求解(11)的方法是块坐标下降法，需要计算若干个迭代，每次迭代有两个步骤：local step 和 global step ，分别说明如下：

1. local step 是固定 $\mathbf{x}$ 计算 $\mathbf{d}$，其实就是把 $\mathbf{d}$ 投影到弹簧原长的位置，具体操作为遍历每个弹簧，找到其两个端点 ${\mathbf{p}}_{i_1}$ 和 ${\mathbf{p}}_{i_2}$，记 ${\mathbf{p}}_{i_{12}}={\mathbf{p}}_{i_1}-{\mathbf{p}}_{i_2}$，那么就更新 ${\mathbf{d}}_i=r_i\frac{{\mathbf{p}}_{i_{12}}}{||{\mathbf{p}}_{i_{12}}||}$ 。

2. global step 是固定 $\mathbf{d}$ 计算 $\mathbf{x}$，也就是求一个凸二次优化问题。$\mathbf{M}+h^2\mathbf{L}$ 是一个对称正定矩阵，求解(11)的最小值可以令其梯度等于 0 ，也就是：
   $$(\mathbf{M}+h^2\mathbf{L})\mathbf{x}=h^2\mathbf{J}\mathbf{d}+\mathbf{M}\mathbf{y}+h^2{\mathbf{f}}_{ext}\text{\hspace{1em}(12)}$$
   因此该步骤就等价于求解一个线性方程组。因为 $\mathbf{M}+h^2\mathbf{L}$ 一直是不变的，所以可以在初始化的时候先做稀疏矩阵 Cholesky 分解（$A=L{L}^T$），可节省大量时间。

如果考虑速度阻尼 $\alpha$，只需令 $\mathbf{y}={\mathbf{q}}_n+\alpha ({\mathbf{q}}_n-{\mathbf{q}}_{n-1})$ 即可。

demo

我基于 fast mass spring 用OpenGL（GLFW）做了布料仿真。布料的配置过程是这样的：

• 正方形布料，边节点数目为 $n$ ，总节点数目为 $n^2$ （$n>1$）；
• 布料宽度为 $w$，相邻节点间隔为 $w/(n-1)$；
• 布料弹簧数目为 $6n^2-10n+2$，其中 structural 类型的弹簧 $2n(n-1)$ 个，shearing类型的弹簧 $2(n-1{)}^2$ 个，bending类型的弹簧 $2n(n-2)$ 个；

设置悬挂节点质量为无穷大，外力为0，得到如下效果：

该算法在 Intel i7-11700K CPU 单核上，3600个节点，2万多个弹簧，每一帧迭代 10次，运行帧率可以超过60帧，非常丝滑~