我的笔记
PBD初探 https://blog.csdn.net/weixin_43940314/article/details/126065813
XPBD 文献笔记 https://blog.csdn.net/weixin_43940314/article/details/126686064
primal dual PBD 文献笔记 https://blog.csdn.net/weixin_43940314/article/details/129792090
PBD相关的软件
Houdini中的Vellum https://www.sidefx.com/docs/houdini/vellum/index.html
Mile Macklin的Flex(闭源仅N卡) https://github.com/NVIDIAGameWorks/FleX
NVidia的 physX (开源N卡) https://github.com/NVIDIA-Omniverse/PhysX
它是一种数值求解的算法，主要用在固体求解上。可以用于布料、软体、流体、刚体、塑性体等。最早是Matthias Muller 等人在2007年提出的。Miles Macklin在2016年提出XPBD是对PBD的一大改进。PBD优点是速度快，稳定。
PBD初探 https://blog.csdn.net/weixin_43940314/article/details/126065813
算法如下
（物体受到的力分为内力和外力。外力就是重力等不会受到物体形变量影响的力。内力则是弹性力等让物体恢复变形的力。内力都可以写成constraint的形式。PBD不像FEM那样去计算受力，然后积分得到下时刻位置。而是直接通过约束的满足来得到位置的修正量。约束的施加就相当于是内力了。）
核心在于标黄的第9到11步。也就是project constraints
我们可以设计出多种多样的constraints，比如体积约束，距离约束，应变约束，流体密度约束，刚体约束，塑性约束等等…
PBD的魅力正在于此。你可以结合多种多样的约束设计出自己的solver。
另外，碰撞和摩擦也可以设计成constraint。我们后面会讲。
project Constraints的公式为
算出lambda，然后带入第一行的公式计算出距离修正量，最后施加距离修正量即可。
XPBD 文献笔记 https://blog.csdn.net/weixin_43940314/article/details/126686064
我的讲解 https://www.bilibili.com/video/BV17D4y1M76v/
XPBD的算法如下
它与PBD的唯一区别就在于lambda
也就是多了第7步和第9步。
公式为
这里的dlambda和原始PBD中的s相比，只是在分母和分子项分别多了alpha项。假如alpha设置为0，那么就会还原为PBD。alpha是材料参数，是柔度，也就是刚度的倒数。这里的波浪线代表alpha除以dt平方：
XPBD的好处在于：通过引入alpha, 让系统的条件数降低了, 因此更容易收敛。并且相比于原始PBD中的k，alpha是更容易调节的材料参数。原始PBD可以认为是alpha为0的系统，因此是刚度无穷大的系统，条件数比较大。
根据刘天添的fast simulation of mass-spring system 中的论证，XPBD相当于是隐式的欧拉时间积分。这就是其稳定的来源。
主要参考文献
Jan Bender, Matthias Müller and Miles Macklin， A Survey on Position Based Dynamics, 2017 主要是第五章
约束公式我们需要知道两点：
对于弹簧连接的两个点x1和x2，最小化下面的约束相当于让弹簧回到原长d。
n就是法向量
二面角约束就是让夹角回到原角度。
n1和n2是两个面的法向量，phi0是原角度
把n1和n2写成用x表示的公式，得到
二面角约束的梯度稍微复杂一点
由于只涉及相对位置，所以不妨我们这里假设p1位置为零点。
最终的距离修正量为
(Bender2014)
除了二面角，还有一种Bending约束，就是在对角线上新增加一个弹簧。然后利用弹簧的距离约束来实现弯曲约束。
其中Q是local Hessian bending energy
Q是一个4x4的矩阵，是常数，可以被precompute。
A0 and A1 是 adjacent triangles 的面积
K是
对于任一点q与三角面之间的碰撞，让两者保证一个最小间距h即可。h就是布料的厚度。
n是面法向量。把n写开，得到公式
另外，值得注意的是三角面有两个方向（q在内部和q在外部，面法向量指向外部）
假设一个粒子的位置为x。先找到所有与其有可能碰撞的平面法向量n，然后施加如下约束。drest是初始情况下的间距。
两个粒子之间保持最小距离。这可以用来模拟沙子。
（Macklin2014）
摩擦比较特殊，没有给出约束，但是给出了距离修正量。我们假设点i与点j之间发生了摩擦。
第一行是静摩擦，第二行是动摩擦。mu_s是静摩擦系数，mu_k是动摩擦系数。
其中delta x 垂直 代表的是x与其摩擦的平面之间的切向滑移的距离（也就是去掉了垂直分量）
其中x*是当前的当前的候选位置，包含了之前所有约束的距离修正量。x是时间步开始的时候的位置。
n是摩擦平面的法向。
与其摩擦的另一个点j的距离修正量为。
很简单，就是让四面体回复到原始体积。
如果是二维，退化为面积约束。
气球就是一块封闭的布料。内部有压力。
（Diziol2011）
一般体积约束用四面体。假如能用三角面也能做体积约束，就会节约很多计算量。
核心思想仍然是让体积回复原始体积
但是，我们可以用散度定理将体积分转换为面积分。其中A是面积。i是三角形编号。n是面法向。
其中上面加一横线表示对所有包含点i的三角面的area weighted normal 。
值得注意的是，这里wi不再是质量倒数，而是
g与l分别代表global 和 local
alpha是可控的参数
还有一个单独控制local的公式。使用参数beta
碰撞物体外面的空气也做成网格，保证其不翻转，就能做碰撞。
上面的是2D，下面的是3D.
这种方法的最大好处是解决entagle state问题。也就是即使压扁成一个面，也能恢复原状。
（Bender14）
这套做法基本上和FEM差不多。都是推导出应变-应力-应变能。几乎可以看做只是把FEM那一套换了个皮。
约束为应变能（也就是弹性势能）
其中V是原体积。psi是能量密度。Es是strain energy。
F是deformation gradient。
约束梯度为
其中P是PK1.
其中C是胡克定律中的刚度矩阵
epsilon是应变。如果采用非线性应变，也就是格林应变
那么这种本构叫做St VK 本构模型
在FEM中，本构模型是可以换的。比如可以换为NeoHooken模型。
这里deformation gradient可以用当前位置计算出来。
其中D是shape matrix。Ds是deformed shape matrix, Dm是undeformed shape matrix。
其中X代表的rest postion。也就是初始位置。所以Dm可以预计算。
(Muller14)
另外一种方案不是计算出能量密度来，而是直接去限制应变。
其中S是
对角项对应拉伸变形
非对角项对应剪切变形
于是就有了以上两个constraint
为了让stretch constraint 的 gradient 为线性的，对第一个公式进行修改：
为了让shear与stretch分开，对第二行公式进行修改
f是F的某一列的列向量。
绳子是一个在三维中的一维物体。但是它是能够产生三维的变形的。绳子能产生拉伸、剪切、扭转、弯曲变形。（stretch/shear/torsion/bending）
首先是stretch+shear
其中p1和p2是一小段绳子的两个顶点。q是四元数，用于描述绳子的方位。d3是沿着这段绳子中心线的切向量。或者说是沿着绳子的截面的法向量，它是四元数的函数。
对应的距离修正量为
其次是Bending+torsion
其中第一行那个花体符号其实是Im(对应latex \Im), 就是Imaginary，也就是虚部。即取四元数的虚部。
其中Omega是Darboux vector。它的物理意义类似于角速度。
其中alpha是保证Darboux vector unique的，因为四元数q和-q都能产生相同的Darboux Vector。
位置修正量为
另外记得每次更新后都要对q归一化。
其中wq本来应该是转动惯量的逆，但是这里简单地用inverse mass 来代替了。
(Deul14)
（注意这个不是shape matching）
描述刚体与粒子系统的不同之处在于，我们要求解的变量不再是粒子的位置，而是刚体的平移和旋转。平移有3个自由度，旋转有3个自由度。
我们以铰链系统为例。它是定义在相对坐标系上的。
某个铰链点的位置根据平移和旋转来决定。左边的P是世界坐标。右边的是物体坐标。
与PBD相同，从对constraint的taylor展开开始。
其中那个花体符号是theta，是描述转轴的向量。它也叫exponential map。
应该注意的是，与之前粒子的PBD不同，增加了theta。所以这里的约束的梯度grad C写成J。
这里举一个例子：用铰链来描述刚体运动。
对于球形铰链连接的两个点p1和p2
其中J为
对于球形铰链来说，第一项就是单位阵
第二项中的第一部分也是单位阵
第三部分则复杂一些，可以参考 Grassia1998
接下来就是PBD的常规流程
根据以下公式得到lambda
然后得到更新量（只不过不仅需要平移，也需要旋转）
最后更新世界坐标系下的速度和角速度。
最简单的流体其实是设定一个粒子之间的最小间距约束。这会产生类似沙子的效果。granular material，颗粒物质。
如果想模拟水，则可以使用PBF。PBF是基于SPH的密度约束的。它几乎是PCISPH套了个壳子。
定义密度约束
密度通过采样得到
gradC只需要利用预计算的核函数梯度
位置修正量
请看我的笔记
https://blog.csdn.net/weixin_43940314/article/details/125961503
https://blog.csdn.net/weixin_43940314/article/details/125952851
https://blog.csdn.net/weixin_43940314/article/details/125913088
shape matching的核心思想就是让变形了的物体回复到原来的shape。
通过一个矩阵A来记忆这个shape。而通过还原的多少来分别模拟刚体（立即完全还原）、弹性体（完全还原但过程缓慢）
和塑性体（部分还原且过程缓慢）
下图中R代表rotation, c代表center of mass translation
最小化的约束为
记忆原始形状的A为
As是对称的，无旋转。
Ar则需要极分解来提取出旋转
要恢复的goal position为
最终的位置修正量为
alpha是0到1的控制参数。可以让物体不回复到原状，制造塑性效果。
Region-Based Shape Matching
Fast Summation
FastLSM
Adaptive Lattices
oriented particles
Plastic Deformation
Large Elasto-Plastic Deformation
方案1：
在Gauss-Seidel的并行中，由于一个顶点会被多个四面体共享，所以在写入位置的时候，会导致race condition。一种方案是利用atomic_add。但是这样会导致很慢。
方案2：
另一种方案是做graph coloring。也就是给约束分组，让互不关联的部分并行运行。比如1-2 3-4 5-6分为一组。2-3 4-5 6-7分为一组。先处理一组，再处理另一组。组内可以并行。也就是一组相当于一个阶段。颜色数决定了需要分多少阶段来运行。
方案3：
还有一种方案是使用Jacobian solver。也就是求解约束以后都先存储到相应的约束对应的position delta之中。最后在施加位置修正量。这种方案比起graph coloring的好处是负载平衡。
方案4：
Jacobian的缺点是它的收敛速度比Gauss-Seidel慢得多，尤其是系统非正定的时候。为了解决这个问题，可以使用欠松弛方法。比如约束平均(constraint-averaging)或者是质量分裂(mass-spliting)方法。基于约束平均的方法如下：
其中omega是松弛因子，建议取1到2之间。而ni代表的是这个粒子i被几个constrait所同时影响。delta xi是一次迭代完全结束之后对于粒子i的位置修正量。由于平均本身就已经做了欠松弛，所以这里omega建议不要取小于1的值，因为已经没有必要进一步欠松弛。
方案5：
还有一种方案是混合GS和Jacobian方法。
一种简单的混合方案是先对k-1个颜色进行GS迭代，然后对剩余的constraint进行Jacobian迭代。
在该框架中，所有一切都用粒子表示。
在shape matching时，首先voxelize一切物体，然后在每个voxel中放置粒子。最后用shape matching可以模拟刚体。物体之间可以增加额外的约束来进行链接。比如可以用刚体和布料连接出来一个降落伞。
每个粒子额外存储一个phase属性。同一个phase内部的粒子之间无碰撞。
约束之间可以组合。比如可以用流体+刚体组成一个融化的模型。
在dev分支
https://github.com/chunleili/tiPBD/tree/dev
这是taichi pbd库。
需要taichi>= 1.4.1和meshtaichi
使用方法：
在data/scene中添加场景文件，在model中添加几何文件（现在仅支持tetgen格式），然后从根目录运行main.py。
python main.py
