学习教程来自:GAMES201：高级物理引擎实战指南2020
以下大部分图片来自教程PPT，仅作为笔记用于学习和分享，侵删


笔记内容大多为课程内容的翻译和转述，外加一些自己的理解，若有不正确的地方恳请大家交流和指正

笔记

MPM和FEM都属于Galerkin methods，MPM不同于FEM方法，其中没有元素。因此MPM ∈ Element-free Galerkin (EFG)

1. Moving Least Squares MPM (MLS-MPM)

移动最小二乘MPM
Y. Hu et al. (2018). “A moving least squares material point method with displacement discontinuity and two-way rigid body coupling”. In: ACM Transactions on Graphics (TOG) 37.4, pp. 1–14

1.1 APIC

3C. Jiang, C. Schroeder, and J. Teran (2017). “An angular momentum conserving affine-particle-in-cell method”. In: Journal of Computational Physics 338, pp. 137–164.

对比PIC方法额外维护了变量C（2x2或3x3的矩阵），记录了粒子周围的affine速度场（ax+b中的a项）


</br>

  1. P2G：求网格上的动量和质量，其中增加affine的部分
  2. Grid operation：从动量求出速度，求解泊松方程（pressure projection）得到散度为0的速度场
  3. G2P：gather速度到粒子上，更新粒子位置，求新的矩阵C

1.2 MLS-MPM


</br>

对比APIC改进的地方：

  1. P2G：将形变梯度引入动量的计算，增加了一个弹力项
  2. Grid operation：计算弹性物体时，pressure projection替换为边界条件的约束计算得到速度
  3. G2P：无改动

1.2.1 形变梯度的计算


</br>

使用矩阵C近似了速度的梯度得到


</br>

1.2.2 弹力项

一个累加到动量上的冲量(impulse)：


</br>

其中f为对弹性势能求导


</br>

求导过程：

  1. （5）带入势能公式
  2. （6）对x hat的求导=对速度的求导乘以τ
  3. （7）链式求导法则带入对F和C的求导
  4. （8）求导，第一项为PK1 stress定义，后2项中右乘的矩阵求导后为矩阵转置
  5. （9）化简（8）


</br>

1.2.3 边界条件 BC

在网格上边做分步积分和散度定理来enforce
3种边界类型：

  1. sticky：粘性的。速度直接归0
  2. slip：滑移边界。速度-边界法线方向的分量，靠近和离开边界时均有阻力
  3. separate：分离边界。只保留靠近边界时的阻力


</br>

其他情况：

  1. 重力对v的作用，施加在边界计算之前
  2. 包含其他碰撞的物体（移动的边界）时，v取相对速度
  3. 边界的摩擦

1.3 MLS-MPM对比MPM的优点

  1. 重用APIC的C矩阵作为形变梯度中∇v的近似，减少了浮点运算
  2. 将形变的更新从G2P阶段提前到P2G阶段，减少了带宽消耗
  3. 在P2G阶段计算动量时，其中APIC和MLS-MPM的动量项可以合并，节省了浮点运算

2. Constitutive Models 本构模型

2.1 一些在MPM框架下实现的材料模型

  1. 弹性物体：NeoHookean & Corotated
  2. 流体：EOS(Equation-of-States) 状态方程法
  3. 弹塑性物体（如雪、沙子）：Yield criteria: ad-hoc boxing, Cam-clay, Drucker-prager, NACC，…

本构模型实现时需要在MPM中注意的点：形变更新的计算、压强的计算

2.2 弹性物体

形变更新部分Deformation update：


</br>

PK1 stress部分（见Games201学习笔记2：拉格朗日视角2）：


</br>

2.3 流体

2.3.1 常规解法

用体积变化比（Volume ratio）求压强p，K为bulk modulus


</br>


</br>

得到Cauchy stress


</br>

数值不稳定性：在求行列式时，可能出现catastrophic cancellation（灾难性消除），指当2个巨大的浮点数相减时出现的精度问题，在流体模型中可能出现，在弹性模型中一般不会出现

如何避免：在粒子中增加对J的存储，而非使用形变梯度F求出新的J


</br>

2.3.2 Lazy solution

使用弹性模型corotated model，把参数µ设置为0

2.4 弹塑性物体

奇异值分解singular value decompositions(SVD)：将实数矩阵分解，U和V为正交矩阵，Σ为对角矩阵，其中Σ对焦线上的值为奇异值（singular values）


</br>

SVD的几何意义：矩阵变换分解为3步，旋转、缩放、旋转


</br>

SVD在形变中的意义：分解形变梯度为3个矩阵，其中塑性形变中可剔除旋转的部分

求SVD中3个矩阵的方法：
A. McAdams et al. (2011). Computing the singular value decomposition of 3x3 matrices with minimal branching and elementary floating point operations. Tech. rep. University of Wisconsin-Madison Department of Computer Sciences.

在计算时，对SVD中的3个矩阵约定如下：

  1. U和V的行列式为1
  2. Σ中元素绝对值降序
  3. 奇异值中的最小值可以为复数

2.4.1 形变梯度

将相变梯度拆分为弹性和塑性2部分，其中只保留弹性部分的势能密度


</br>

一个例子：
1A. Stomakhin et al. (2013). “A material point method for snow simulation”. In: ACM Transactions on Graphics (TOG) 32.4, pp. 1–10.


</br>

过程分析：

  1. 先更新弹性形变梯度作为$\hat{F}^{n+1}_{p}$
  2. 对$\hat{F}^{n+1}_{p}$求SVD矩阵
  3. 裁剪掉Σ中拉伸和压缩超出范围的值
  4. 用裁剪后的Σ作为弹性形变梯度，裁剪掉的值作为塑性形变梯度

3. MPM中使用拉格朗日力

Lagrangian forces in MPM

将MPM中的粒子作为有限元方法的顶点，计算采用FEM的势能模型，不再需要像上边一样计算形变梯度
好处：

  1. 对比FEM：能处理自碰撞
  2. 对比MPM：由于使用了Mesh而不是Grid，避免MPM中由于采样密度不够高导致的 Numerical fracture（数值断裂）
  3. 更容易耦合MPM和FEM
赏