参考：
Matthia Muller的十分钟物理（他就是PBD算法的发明者）
https://matthias-research.github.io/pages/tenMinutePhysics/

原理

PBD的算法主体分为三步：

1. 根据外力更新粒子速度位置，无需考虑粒子间关系
2. 求解约束，使粒子满足粒子间关系
3. 更新粒子的位置，并反向更新速度

由于第3步是先求出粒子位置，再反求速度的。因此它是基于位置的方法，故被称为position based dynamics。

PBD是纯粒子法。在PBD的世界中，只有粒子。其他的一切（三角面、四面体等）都是辅助性的。

1. 根据外力更新粒子速度位置，无需考虑粒子间关系

如图，先根据外力更新速度（这里只考虑了重力）

然后把旧的位置存到p中

最后利用速度更新位置

（如果有碰撞，也在这里处理）

整个过程粒子无需考虑相互关系。

2. 求解约束，使粒子满足粒子间关系

我们分两种情况：一种是最简单的二维情况，一种是三维情况。

我们先从简单的二维情况来：

二维情况：一个弹簧

我们只考虑一个弹簧的约束

这个弹簧只有两个质点和中间的一个边。质点具有位置和质量，边具有原长度和现长度。

随意移动两个粒子，弹簧目前不处于原长。

因此，弹簧要回到原长。

弹簧回到原长，这就是这个系统的约束。

$$C(x1,x2)=|x1-x2|-L0$$
其中x1, x2分别是粒子1和粒子2的位置。L0是原长。

如何让系统满足约束呢？

所谓的满足约束，就是让约束误差等于0。

上面这个式子中的C，就是约束的误差（有时候约束和约束误差这两个词混用）。

通过迭代，让误差趋于0，那就是求解过程。

如何让其趋于0？

那就是求梯度，然后向着梯度减小的方向移动（即所谓的梯度下降思想）。

实际上，从另一个角度来看，梯度就是一维的导数在高维的推广。让导数等于0的位置（即驻点），就是原函数最小化的位置。

于是我们就找C的梯度。

$$\nabla C1=(x1-x2)/|x1-x2|\nabla C2=(x2-x1)/|x2-x1|$$
我们无需那样严谨地去推导数学公式，然后给出C梯度的表达式。我们直接从物理意义上来理解C梯度。那就是让C函数上涨最快的方向。在这个例子中，也就是x1与x2相对的方向。因此我们给出上面的公式。而且我们这里是归一化了的，只求出一个方向。

那么大小是多少呢？我们给出如下公式

大小是由系数lambda和质量倒数w决定的。其中lambda的正式名称叫做拉格朗日乘数。

那么，我们求解出来了gradC，因此就求解出来了粒子间的相互关系。因此，也就知道粒子为了满足相互关系，该向哪里移动。因此给出了dx（如上图）。这个移动的方向，就是最小化约束误差的方向，就是负梯度的方向（因此lambda分母有个负号）。

三维情况：一个四面体

对于三维，我们期望的约束是四面体的体积保持原体积。

其约束误差就是
$$C=(V-V0)$$

而四面体的体积公式是
$$\frac{1}{6}[(x2-x1)\times (x3-x1)]\cdot (x4-x1)$$
先叉乘求出底面积（叉乘得到的是菱形面积，还要除以2），然后再乘以高度（点乘会消除非垂直部分），再乘以1/3(四面体公式中本来的1/3)

而约束误差的梯度是什么呢？我们仍然跳过数学推导，从物理意义上解释。

梯度即函数增长最快的方向，也就是体积增长最快的方向。对于某个粒子来说，哪个方向让体积增长最快呢？

那就是垂直于底面的方向。

而垂直与底面的方向，就是底面的三角形其中两个边叉乘的方向（叉乘满足右手定则）。

因此
$$gradC4=(x2-x1)\times (x3-x1)$$
这就是粒子4所对应约束的梯度。

那么大小如何确定呢？

仍然利用拉格朗日乘数lambda。

值得注意的是：
求解一个三维弹性物体，我们必须同时满足弹簧约束和体积约束。体积约束保证弹性体不发生体积的膨胀和缩小，而弹簧约束保证物体上的每个质点回到原位（也就是最终弹性体恢复原状）。

代码

我们参考的是
https://matthias-research.github.io/pages/tenMinutePhysics/
中第10讲的代码

数据结构：
存储在class SoftBody内。SoftBody可多次实例化以添加多个物体。

simulate函数为算法的主体，它分为三步：

1. preSolve
2. solve
3. postSolve

preSolve

for all particles:
	vel[i] += gravity * dt
	prevPos[i] = pos[i]
    pos[i] += vel[i] *dt
    collision处理：当y<0时把pos挪到0
1
2
3
4
5

solve

又被分为两个部分：

1. solveEdge
2. solveVolume

也就对应上面说的弹簧约束和体积约束。

solveEdge

for i in all edges:
	alpha = C/dt^2
	grad = pos0- pos1
	grad归一化
	s = - (L - L0)/(1/m + alpha)
	pos0[i] += grad * (s/m)
	pos1[i] += grad * (-s/m)     
1
2
3
4
5
6
7

solveVolume

for i in all tets:
	alpha = C/dt^2
	for j in 4:
		temp0 = pos1 - pos0
		temp1 = pos2 - pos0
		grad[j] = (1/6) * tem0.cross(temp1)
	w = sum((1/m) * ||grad[j]||^2 )
	s = -(V - V0) / (w + alpha)
	for j in 4:
		pos[i] = grad[j] * (s/m)
1
2
3
4
5
6
7
8
9
10

postSolve

for i in all particles:
	vel[i] = (pos[i] - prevPos[i])/dt
1
2

代码

拷贝自tenMinutePhysics

---