Interactive Simulation of Elastic Deformable Materials论文中的公式纠错

1.简介

布料仿真中我们通常将布料离散化为三角网格或四边形网格，并用弹簧质点模型对其建模，最终得到如下的动力学方程：

$$\begin{equation} \begin{cases} \dot{x}=v\\ M\dot{v}=F_{external} + F_{internal} \end{cases} \end{equation}$$
根据内力的表示 $F_{internal}$的表示不同，可分别到得到（1）基于惩罚力（penalty force）的动力学方程：
$$\begin{equation} \begin{cases} \dot{x}=v\\ M\dot{v}=F_{external} - kJ^Tc(x) \end{cases} \end{equation}$$
（2）基于约束力（constraint force）的动力学方程：
$$\begin{equation} \begin{cases} \dot{x}=v\\ M\dot{v}=F_{external} + J^T\lambda\\ c(x)=0 \end{cases} \end{equation}$$
其中， $J=\nabla c(x)$， $c(x)$是约束，例如形如 $c(x_1,x_2)=||x_1-x_2||-l_0=0$的距离约束。

　这篇文章通过令约束力（constraint force）与惩罚力（penalty force）相等，得到如下动力学方程（称之为 soft constrained dynamics）：

$$\begin{equation} \begin{cases} \dot{x}=v\\ M\dot{v}=F_{ext} + J^T\lambda\\ \alpha\lambda(x,t)=-c(x) - \beta \dot{c}(x) \end{cases} \end{equation}$$
并用类似于Symplectic Euler的方法进行求解。

2.求$\lambda$的区间均值

用$\lambda$在$[t,t+\Delta t]$（实际上应该是$[t_n,t_n+\Delta t]$，所以论文里的积分上下限也写错了，下面的是更正后的）的均值$\tilde{\lambda}$代替$\lambda$,$\tilde{\lambda}=\Delta t^{-1}\int_{t_n}^{t_n+\Delta t}\alpha^{-1}(c+\beta\dot{c})dt$
将$c(x)$在$t_n$处做一阶泰勒展开，得到，$c(x)=c_n+\dot{c}(t-t_n)$，则，

$$\tilde{\lambda}=\Delta t^{-1}\alpha^{-1}\int_{t_n}^{t_n+\Delta t}[c_n+\dot{c}|_n(t-t_n)+\beta \dot{c}]dt$$
$\dot{c}=J^TV$是 $c$对时间求导，令$\dot{c}=\dot{c}|_n$得，
$$\tilde{\lambda}=\Delta t^{-1}\alpha^{-1}\Big(\Delta tc_n+\frac{1}{2}\Delta t^2\dot{c}|_n+\beta \Delta t \dot{c}|_n\Big)=\alpha^{-1}\Big(c_n+(\frac{1}{2}\Delta t+\beta )J_n^TV_{n+1}\Big)$$
其中， $\dot{c}|_n=J_n^TV_{n+1}$（ 我也不知道为什么不是$V_n$）令， $\gamma=(\frac{1}{2}\Delta t+\beta )^{-1}$，则，
$$\alpha \tilde{\lambda}=c_n+\gamma^{-1}J_n^TV_{n+1}$$
最终得到如下方程，
$$\begin{equation} {\left[ \begin{array}{cc} M & -J_n^T\\ J_n& \gamma\alpha\Delta t^{-1} \end{array} \right ]} {\left[ \begin{array}{c} V_{n+1} \\ \Delta t\tilde{\lambda} \end{array} \right ]}= {\left[ \begin{array}{c} MV_n+\Delta tF_n \\ -\gamma C_n \end{array} \right ]} \end{equation}$$

　实际中，我们会固定一些点，这一般通过令固定点的质量无穷大实现，而计算机中无法表示数学意义上的无穷大，因此使用质量的逆等于0来等价表示。因此通常我们都会让公式只出现质量的逆矩阵$M^{-1}$（由于$M$是对角矩阵，求逆耗时很少，而且是once and for all的工作），因此上式改写成如下形式：
　

$$\begin{equation} {\left[ \begin{array}{cc} I & -M^{-1}J_n^T\\ J_n& \gamma\alpha\Delta t^{-1} \end{array} \right ]} {\left[ \begin{array}{c} V_{n+1} \\ \Delta t\tilde{\lambda} \end{array} \right ]}= {\left[ \begin{array}{c} V_n+\Delta tM^{-1}F_n \\ -\gamma C_n \end{array} \right ]} \end{equation}$$
注意，在我的实际操作中，求解这个矩阵只能用LU分解，Eigen库的Conjugate Gradient和SimplicialLDLT都会导致结果不正确。

3.Schur complement

根据上面得到的正确系数矩阵，可得到其舒尔补（Schur complement）为，

$$S_n=J_nM^{-1}J_n^T-\Delta t^{-1}\gamma\alpha$$
$$S_n\tilde{\lambda}=-\Delta t^{-1}J_nV_n-J_nM^{-1}F_n-\gamma\Delta t^{-1}C_n$$
其中， $F_n$表示非约束力之外的力。最终得到速度和位置的更新公式：
$$\begin{equation} \begin{cases} x_{n+1}=x_n+\Delta t v_{n+1}\\ v_{n+1}=v_n+hM^{-1}F_n+hM^{-1}J_n^T\tilde{\lambda} \end{cases} \end{equation}$$

如果使用原来错误的公式，仿真一块两点固定的布料时，布料的右侧会出现ghost force：不正常的被向前拉。而使用纠正后的公式表现正常。如下图所示：

4.总结

1. 本文纠正了论文Interactive Simulation of Elastic Deformable Materials公式推导的错误，给出了写程序时应该使用的求解方法建议。
2. 有关基于约束的动力学（constrained dynamics）的推导和求解，如果有人感兴趣可以在下面留言，有机会我会介绍。
3. 本文所讲的论文（06年发表）提出了soft constraint（这篇论文里称之为regularized constraint），15年，16年都有论文重新讲述soft constraint的求解，以及如何应用到PBD（position based dynamics）上。同样的，感兴趣的可以在下面留言，有机会我会介绍。
4. 在实际求解中我发现对于不同的稀疏系数矩阵，需要使用不同的求解方法（我用的是Eigen），否则得不到解，具体原因之后会详细给出。
5. 截图使用的仿真环境为GitHub上PBD原作者公布的源代码:https://github.com/InteractiveComputerGraphics/PositionBasedDynamics.