【物理模拟】PBD算法详解

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

原理

PBD的算法主体分为三步:

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

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

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

算法主体

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

如图,先根据外力更新速度(这里只考虑了重力)

然后把旧的位置存到p中

最后利用速度更新位置

(如果有碰撞,也在这里处理)

整个过程粒子无需考虑相互关系。
在这里插入图片描述

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

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

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

二维情况:一个弹簧

我们只考虑一个弹簧的约束
在这里插入图片描述

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

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

因此,弹簧要回到原长。

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

C ( x 1 , x 2 ) = ∣ x 1 − x 2 ∣ − L 0 C(x1, x2) = |x1 - x2| - L0 C(x1,x2)=x1x2∣L0
其中x1, x2分别是粒子1和粒子2的位置。L0是原长。

如何让系统满足约束呢?

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

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

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

如何让其趋于0?

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

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

于是我们就找C的梯度。

∇ C 1 = ( x 1 − x 2 ) / ∣ x 1 − x 2 ∣ ∇ C 2 = ( x 2 − x 1 ) / ∣ x 2 − x 1 ∣ \nabla C1 = (x1 - x2) / |x1 - x2| \nabla C2 = (x2 - x1) / |x2 - x1| C1=(x1x2)/∣x1x2∣∇C2=(x2x1)/∣x2x1∣
我们无需那样严谨地去推导数学公式,然后给出C梯度的表达式。我们直接从物理意义上来理解C梯度。那就是让C函数上涨最快的方向。在这个例子中,也就是x1与x2相对的方向。因此我们给出上面的公式。而且我们这里是归一化了的,只求出一个方向。
在这里插入图片描述
那么大小是多少呢?我们给出如下公式
在这里插入图片描述
大小是由系数lambda和质量倒数w决定的。其中lambda的正式名称叫做拉格朗日乘数。

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

三维情况:一个四面体

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

其约束误差就是
C = ( V − V 0 ) C = (V - V0) C=(VV0)

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

在这里插入图片描述

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

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

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

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

因此
g r a d C 4 = ( x 2 − x 1 ) × ( x 3 − x 1 ) gradC4 = (x2 - x1) \times (x3 - x1) gradC4=(x2x1)×(x3x1)
这就是粒子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

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)     

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)

postSolve

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

代码

拷贝自tenMinutePhysics


在这里插入图片描述

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

  • 2
    点赞
  • 3
    收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:技术黑板 设计师:CSDN官方博客 返回首页
评论

打赏作者

beidou111

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值