算法描述

物体的动力学可由 $N$ 个顶点和 $M$ 个关于顶点位置的约束进行描述．第 $i$ 个顶点的质量、位置、速度分别记作 $m_i,\mathbf{x}_i,\mathbf{v}_i$．第 $j$ 个约束可能为等式或不等式，记作

$$C_{j}(\mathbf{x}_{1}, \cdots, \mathbf{x}_{N})=0\quad \text{或}\quad C_{j}(\mathbf{x}_{1}, \cdots, \mathbf{x}_{N})\ge0.$$

我们用符号 $\succ$ 来表示 $=$ 或 $\ge$，并且令 $\mathbf{x}=(\mathbf{x}_{1}, \cdots, \mathbf{x}_{N})$，则第 $j$ 个约束可简记为 $C_{j}(\mathbf{x})\succ 0$．

基于位置的动力学（Position Based Dynamics）简称为PBD．该算法会先计算出每个顶点的初始预测位置 $\mathbf{p}_i(i=1,2,\cdots,N)$，然后尝试移动 $\mathbf{p}_i(i=1,2,\cdots,N)$ 使其满足所有约束

$$C_j\succ 0(j=1,2,\cdots,M).$$

PBD的具体算法可由以下伪代码表示．

图1 - PBD 伪代码

求解器

在PBD伪代码中，第8至第10行是求解器的迭代求解．这里的函数

$$\mathrm{ projectConstraints }\left(C_{1}, \cdots, C_{M + M_{\mathrm{Coll }}}, \mathbf{p}_{1}, \ldots, \mathbf{p}_{N}\right)$$

是方程/不等式组

$$\left\{ \begin{array}{l} C_{1}\left(\mathbf{p}_{1}, \cdots, \mathbf{p}_{N}\right)\succ 0, \\ C_{2}\left(\mathbf{p}_{1}, \cdots, \mathbf{p}_{N}\right)\succ 0, \\ \quad\ldots\\ C_{M + M_{\mathrm{Coll }}}\left(\mathbf{p}_{1}, \cdots, \mathbf{p}_{N}\right)\succ 0, \\ \end{array} \right.$$

的一个非线性Gauss-Seidel求解器．求解器 $\mathrm{ projectConstraints}$ 会对约束 $C_{1},C_{2},\cdots,C_{M + M_{\mathrm{Coll }}}$ 依次进行一阶线性化求解．每次求解完 $C_{j}$，求解器会当即更新 $\mathbf{p}_{1},\mathbf{p}_{2},\cdots,\mathbf{p}_{N}$ 的值，使之影响对 $C_{j+1}$ 的求解．下面我们将具体介绍 $C_{j}$ 的求解方法．

如果约束是一个等式，我们希望找到一个 $\Delta \mathbf{p}$ 使得 $C_{j}(\mathbf{p}+\Delta \mathbf{p})=0$，其中 $\mathbf{p}=\left(\mathbf{p}_{1}, \ldots, \mathbf{p}_{N}\right)$．对其进行一阶线性化可得

$$C_{j}(\mathbf{p}+\Delta \mathbf{p})\approx C_{j}(\mathbf{p})+\nabla C_{j}(\mathbf{p}) \cdot \Delta \mathbf{p} = 0$$

上述方程的一般会有无穷个解．如果我们将解 $\Delta \mathbf{p}=\left(\Delta\mathbf{p}_{1},\Delta\mathbf{p}_{2},\cdots,\Delta\mathbf{p}_{N}\right)$ 限制在

$$\left(w_1\nabla_\mathbf{p_1} C_{j}(\mathbf{p}),w_2\nabla_\mathbf{p_2} C_{j}(\mathbf{p}),\cdots,w_N\nabla_\mathbf{p_N} C_{j}(\mathbf{p})\right)$$

的方向上，即令

$$\begin{equation} \Delta \mathbf{p}_i=\lambda w_i\nabla_{\mathbf{p}_i} C_{j}(\mathbf{p}),\quad i=1,2,\cdots,N, \end{equation}$$

我们就可以求出一个特解

$$\Delta \mathbf{p}_i=-\frac{w_iC_j\left(\mathbf{p}\right)}{\sum_{k=1}^N w_{k}\left|\nabla_{\mathbf{p}_{j}} C_j\left(\mathbf{p}\right)\right|^{2}}\nabla_{\mathbf{p}_i} C_{j}(\mathbf{p}),\quad i=1,2,\cdots,N.$$

$C_{j}$ 的求解方法，就是用 $\mathbf{p}+\Delta \mathbf{p}$ 来更新原来的 $\mathbf{p}$．因此，$\mathrm{projectConstraints}$ 函数可用伪代码表示成

图2 - $\mathrm{projectConstraints}$ 求解器伪代码

特定约束

拉伸

拉伸约束为

$$C(\mathbf{p}_1,\mathbf{p}_2) = \left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert - d,$$

其中 $d$ 代表初始长度。计算梯度，可得

$$\begin{aligned} \nabla_{\mathbf{p}_1} C(\mathbf{p}_1,\mathbf{p}_2)&=\nabla_{\mathbf{p}_1} \left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert\\ &=\nabla_{\mathbf{p}_1} \left(\left(\mathbf{p}_{1x}-\mathbf{p}_{2x}\right)^2+\left(\mathbf{p}_{1y}-\mathbf{p}_{2y}\right)^2+\left(\mathbf{p}_{1z}-\mathbf{p}_{2z}\right)^2\right)^{\frac{1}{2}}\\ &=\frac{1}{2\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert}2(\mathbf{p}_1-\mathbf{p}_2)\\ &=\frac{\mathbf{p}_1-\mathbf{p}_2}{\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert}, \end{aligned}$$ $$\begin{aligned} \nabla_{\mathbf{p}_2} C(\mathbf{p}_1,\mathbf{p}_2)&=\nabla_{\mathbf{p}_2} \left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert\\ &=\nabla_{\mathbf{p}_2} \left(\left(\mathbf{p}_{1x}-\mathbf{p}_{2x}\right)^2+\left(\mathbf{p}_{1y}-\mathbf{p}_{2y}\right)^2+\left(\mathbf{p}_{1z}-\mathbf{p}_{2z}\right)^2\right)^{\frac{1}{2}}\\ &=\frac{1}{2\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert}2(\mathbf{p}_2-\mathbf{p}_1)\\ &=\frac{\mathbf{p}_2-\mathbf{p}_1}{\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert}. \end{aligned}$$

于是

$$\begin{aligned} \lambda&=\frac{\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert - d}{w_1+w_2},\\ \Delta \mathbf{p}_1&=-\frac{\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert - d}{w_1+w_2}w_1\frac{\mathbf{p}_1-\mathbf{p}_2}{\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert}=\frac{w_1}{w_1+w_2}\left(1-\frac{d}{\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert}\right)(\mathbf{p}_2-\mathbf{p}_1),\\ \Delta \mathbf{p}_2&=-\frac{\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert - d}{w_1+w_2}w_2\frac{\mathbf{p}_2-\mathbf{p}_1}{\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert}=\frac{w_2}{w_1+w_2}\left(1-\frac{d}{\left\Vert \mathbf{p}_1-\mathbf{p}_2\right\Vert}\right)(\mathbf{p}_1-\mathbf{p}_2). \end{aligned}$$

SDF碰撞

SDF碰撞约束为

$$C(\mathbf{p}_i) = \mathrm{SDF}(\mathbf{p}_i).$$

矫正量

$$\Delta\mathbf{p}_i = - \frac{\mathrm{SDF}(\mathbf{p}_i)}{w_i\left\Vert\nabla_{\mathbf{p}_i}\mathrm{SDF}(\mathbf{p}_i)\right\Vert^2}w_i\nabla_{\mathbf{p}_i}\mathrm{SDF}(\mathbf{p}_i)= -\frac{\mathrm{SDF}(\mathbf{p}_i)}{\left\Vert\nabla_{\mathbf{p}_i}\mathrm{SDF}(\mathbf{p}_i)\right\Vert^2}\nabla_{\mathbf{p}_i}\mathrm{SDF}(\mathbf{p}_i).$$

例子：球面SDF碰撞

给定半径为 $r$ 的球面，它的SDF碰撞约束为

$$C(\mathbf{p}_i) = \left\Vert\mathbf{p}_i\right\Vert-r.$$

矫正量

$$\Delta\mathbf{p}_i = -\frac{\left\Vert\mathbf{p}_i\right\Vert-r}{\left\Vert\mathbf{p}_i\right\Vert^2/\left\Vert\mathbf{p}_i\right\Vert^2}\frac{\mathbf{p}_i}{\left\Vert\mathbf{p}_i\right\Vert}=\left(\frac{r}{\left\Vert\mathbf{p}_i\right\Vert}-1\right)\mathbf{p}_i.$$

数值微分的计算

我们可以使用四面体技巧对 SDF 数值微分的计算进行优化。设

$$\begin{aligned} \mathbf{k}_0&=(1,-1,-1)\\ \mathbf{k}_1&=(-1,-1,1)\\ \mathbf{k}_2&=(-1,1,-1)\\ \mathbf{k}_3&=(1,1,1)\\. \end{aligned}$$

考虑求和

$$\begin{aligned} m &=\sum_{i=0}^3 \mathbf{k}_{i}\frac{f\left(\mathbf{p}+h \mathbf{k}_{i}\right)}{\sqrt{3}h} \\ &=\sum_{i=0}^3 \mathbf{k}_{i} \frac{f\left(\mathbf{p}+h \mathbf{k}_{i}\right)-f(\mathbf{p})}{\sqrt{3}h}\\ &\approx\sum_{i=0}^3\mathbf{k}_{i} \nabla_{\widehat{\mathbf{k}}_{i}}f(\mathbf{p})\\ &=\sum_{i=0}^3 \begin{pmatrix} k_{ix}\\ k_{iy}\\ k_{iz}\\ \end{pmatrix} \left(\widehat{k}_{ix},\widehat{k}_{iy},\widehat{k}_{iz}\right) \begin{pmatrix} \nabla_{x}f(\mathbf{p})\\ \nabla_{y}f(\mathbf{p})\\ \nabla_{z}f(\mathbf{p}) \end{pmatrix}\\ &=\frac{1}{\sqrt{3}}\sum_{i=0}^3 \begin{pmatrix} k_{ix}^2&k_{ix}k_{iy}&k_{ix}k_{iz}\\ k_{ix}k_{iy}&k_{iy}^2&k_{iy}k_{iz}\\ k_{ix}k_{iz}&k_{iy}k_{iz}&k_{iz} \end{pmatrix} \begin{pmatrix} \nabla_{x}f(\mathbf{p})\\ \nabla_{y}f(\mathbf{p})\\ \nabla_{z}f(\mathbf{p})\\ \end{pmatrix}\\ &=\frac{4}{\sqrt{3}}\nabla f(\mathbf{p}) \end{aligned}.$$

由此可以得到

$$\nabla f(\mathbf{p})\approx\sum_{i=0}^3 \mathbf{k}_{i}\frac{f\left(\mathbf{p}+h \mathbf{k}_{i}\right)}{4h}$$