Position Based Dynamics简介

发表于 2021-07-20 分类于 CG

前言

最近一直在研究关于仿真模拟的PBD[2006]方法，整了一堆论文看。下一步的目标是找一个框架把这些方法实现出来看看效果，甚至进一步试试有没有改进的空间，在此之前，先对之前的工作做一个小小的总结吧。

PBD基本方法简述

作为一种仿真模拟的一种新兴方法，PBD的目标与其他方法一致，即通过某个时刻每个粒子的位置，结合其他一些速度或者相互作用的信息，计算推断出下一个时刻每个粒子的位置，如此即可形成物体完整的运动轨迹。一般最容易想到的方式就是利用牛顿第二定律通过每个粒子受到的合力求出加速度，然后依次更新速度与位置信息。这样的方法虽然快捷，但是非常容易产生发散的情况。

一种最为常见的处理办法是将这种显式格式变成隐式的以提高稳定性，即隐式欧拉方法： \[ \begin{cases} x_{n+1}=x_n+hv_{n+1}\\ v_{n+1}=v_n+hM^{-1}(f_{int}(x_{n+1})+f_{ext}) \end{cases} \] 其中\(h,M\)分别为时间步长与质量矩阵。可以看到，它与显式方法的不同就是更新速度时要用新的位置信息来计算内力，这也就意味着我们无法和之前一样简单地直接递推得到新的位置，而要每一步求解一个维数为粒子数的一个超大型非线性方程组，时间成本会比较高。

PBD方法(Position based dynamics)顾名思义，是一种基于位置的动力学模拟方式。它放弃了一部分的物理准确性，来追求处理的效率性。事实上如今，在一些对于模拟精度要求不是那么高的场景下，PBD已然成为了一种非常常用的模拟框架。

PBD步骤

PBD方法的基本思路就是首先通过当前的速度与外力信息直接推算出每个粒子大致落在的位置，再对满足约束的空间进行投影，如此便得到了新的位置。举个例子：

当前有两个不受外力的质点，它们之间的作用力看成一根弹簧，已知质点的位置与速度如图，假设开始时弹簧处于原长状态。

那么首先直接通过单位时间步长的速度估算出新的位置，但是我们发现弹簧被拉长了，实际作用中由于弹簧拉力的存在，这样的运动是不会发生的，此时我们要对质点的位置做一些调整。

在弹簧的例子中，我们就是尽量把质子位置往它们距离恰为弹簧原长的空间上投影，如下图，接着再更新速度。这样，在时间步长小的时候，产生的效果就有些像是真实的运动了。

接下来要处理的问题就是如何将顶点位置投影到符合约束的空间中去。

约束投影

约束可以用关于位置信息的一个等式或不等式来表示，即 \[ C(p_1,\cdots,p_n)=0~or~\ge 0 \] 例如弹簧的约束就可以用 \[ C(p_1,p_2)=|p_1-p_2|-d=0 \] 来表示。而碰撞约束这种，保持一个点始终在某一个面一侧的约束就是典型的不等式约束。

先考虑最简单的等式约束以及质量均相等的情形。约束如何满足呢？由于\(C\)关于刚体变换是不变的，则\(\nabla_p C\)就与刚体变换组成的空间垂直，从而沿此方向变化可以保证动量守恒。具体变化多少？令\(C(p+\Delta p)\approx C(p)+\nabla_p C\cdot \Delta p=0,\Delta p=\lambda \nabla_p C\)，于是可求得\(\lambda=-\frac{C(p)}{|\nabla_p C(p)|^2}\)即可。

考虑质量，只需关于逆质量加权即可 \[ \Delta p_i=-s\omega_i \Delta_{p_i}C(p_1,\cdots,p_n) \]

\[ s=\frac{C(p_1,\cdots.p_n)}{\sum_j\omega_j |\nabla_{p_j}C(p_1,\cdots,p_n)|^2} \]

为了表示约束的强弱（类似于弹簧的弹性系数的概念），简单的处理办法只需要在投影时\(\Delta p\)前额外再乘一个系数即可。

但是，约束总是很多的。如果整体求解这样的约化后的线性系统，又将遇到许多难题。论文中提出的就是最简单易行的Gauss-Seidel式的迭代求解，即每个约束依次投影，重复总过程若干遍即可得到近似解。

若干优化

这样的原始算法显然是有许多优化空间的。下面简述几种：

解耦系数关系

其实在当前算法下，物理性质的表现与约束顺序、约束投影迭代次数、时间步长等基本参数都有很强的关联性，这与我们的期望显然是不符的。有一些算法在一定程度上缓解了这些问题：

HPBD

Hierarchical Position Based Dynamics[2008]采用了一种类似于多分辨率分析的思路，将原始网格变为分层级的网络结构。首先移动框架再逐步调整细节，就可以缓解约束求解顺序造成的收敛慢的问题。

XPBD

Extended Position Based Dynamics[2016]从牛顿定律出发，结合力与势能的关系，通过对于Lagrange乘子的同时优化，采用拟牛顿法求解可得新的投影公式 \[ \Delta \lambda_j =\frac{-C_j(x_i)-\tilde{\alpha_j}\lambda_{ij}}{\nabla C_jM^{-1}\nabla C_j^T+\tilde{\alpha_j}} \]

\[ \Delta x=M^{-1}\nabla C(x_i)^T\Delta\lambda \]

据称，可以避免时间步长严重影响物理性质，等实现出来可以尝试一下。

最优化思路改善结果+加速

上面提到Lagrange乘子，我当时就在想这个PBD方法是不是可以用优化的思路进行改善，然后就看到了这些论文~~，失去了发论文的机会啊（不是）~~。

总之就是PBD的更一般的方法Projective Dynamics与Implicit Euler方法的结合[2014]，通过Local/Global求解器可以得到更整体优化的方法。于此之外类似的思路进一步有Chebyshev迭代[2015]与拟牛顿迭代[2016]的方法，可以使之收敛得更快。

GPU并行

GS迭代需要同时读写一些顶点位置信息，对于并行是不安全的。一种方法是采用Jacobi式的迭代代替，最终对位移取均值，或者使用超松弛等方式加速。

另一种方式是尝试让不同的约束处理不同时读写相同顶点的坐标。对此只要先将约束分为没有重叠顶点的小组，每个组内并行即可。

增强物理性

目前为止，我们处理弹性体的方式都是将其视为弹簧，通过调整弹簧弹性系数来模拟物体实际属性，但是这种方法显然有很强的局限性。比如，我们很难去表现各向异性的物体性质。而且我们实际上对于弹性力学已经有了越来越深的研究，一种非常自然的想法就是用能量来表示约束（例如之前的约束为弹性势能为零）。这就又有基于四面体网格[2014]与基于SPH[2021]两种想法，通过位置表示变形梯度，再进一步表示能量密度即可。

结尾

纸上谈兵先写这么多吧，如果后面实现出了效果再更。