算法描述

物体的动力学可由 \(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} \]