MLS-MPM 实时物理仿真：从数学简化到源码实现

简介本资源为基于移动最小二乘物质点法MLS-MPM的高性能物理仿真求解器源码对应 ACM 图形交易SIGGRAPH 2018论文工作支持切割、位移不连续性及双向刚体耦合CPIC采用 C 编写适合图形学、计算力学与物理引擎方向的研究者及中高级开发者研读与二次开发。压缩包共 175 个文件约 600KB以 57 个 h 头文件与 24 个 cpp 源文件为核心辅以 23 个 py 脚本、20 个 c 文件及 Makefile、configure.ac 等构建配置另含许可证与说明文档目录结构清晰。目前已有 963 人学习下载。读者可获取完整求解器实现理解 MLS-MPM 算法在切割与刚体耦合场景下的工程落地方式并借助构建脚本快速编译验证为图形学仿真研究与项目复现提供可靠参考。1. 物质点法为什么能在实时里跑起来从 MLS-MPM 的数学简化说起如果你写过基于网格的流体或弹性体仿真大概体会过那种“网格畸变到没法算”的绝望。拉格朗日网格在超大形变下会缠成一团欧拉网格又很难追踪材料的历史状态。物质点法MPM把两者拼在一起粒子携带质量、速度、形变梯度背景网格只负责算力再插值回去。听起来很美但传统 MPM 每一步都要在网格上做大量积分实时应用里根本跑不动。MLS-MPM 的核心贡献是把移动最小二乘MLS近似直接嵌进物质点法的动量方程里让原本需要复杂应力发散计算的步骤退化成一次仿射速度场的构造。这个简化不是近似掉精度而是把数学结构重新组织使得每一步的粒子到网格、网格到粒子的传递可以用更少的浮点运算完成。对做实时物理的开发者来说这意味着同样的硬件预算下粒子数可以翻几倍或者帧率从个位数拉到可交互区间。搞懂它为什么能 work比会调一个求解器参数重要得多。2. 从连续介质到粒子网格MLS-MPM 的底层机制拆解2.1 物质点法的三件套粒子、网格、形变梯度物质点法里每个粒子至少存三样东西质量、速度、形变梯度。质量决定动量交换速度是插值对象形变梯度记录材料从初始构型到当前构型的局部变形。背景网格不存永久数据每一步开始前清空粒子把质量、动量、应力相关的量散射到网格节点上网格解出新的速度场再插值回粒子。这个循环里最贵的一步是应力计算因为形变梯度要先做极分解或奇异值分解再根据本构模型算柯西应力最后求应力发散。传统 MPM 在这一步的浮点操作数很高而且粒子数一多散射和收集的随机访存也会成为瓶颈。MLS-MPM 的切入点就在这里。它不再显式地算应力发散而是利用移动最小二乘近似把动量方程里的应力项改写成仿射速度场的贡献。具体来说粒子对网格节点的影响不再只是一个标量权重乘以速度而是一个仿射矩阵乘以相对位置。这个仿射矩阵直接由形变梯度和本构模型决定省掉了中间应力张量的构造和发散计算。常见的实现思路是把本构模型给出的第一皮奥拉-基尔霍夫应力直接映射成一个 3x3 的仿射矩阵然后粒子到网格的散射就变成一次矩阵向量乘。这一步的数学推导在原始论文里有完整过程但作为实现者你只需要知道原本的应力发散项被吸收进了仿射速度场的构造里计算量降了一个档次。2.2 仿射速度场为什么一次矩阵乘就能替代应力发散仿射速度场的意思是粒子附近的速度分布可以用一个线性变换来近似。在 MLS 近似下粒子对周围空间的影响不是各向同性的而是沿着形变梯度的方向被拉伸或压缩。这个各向异性正好对应材料内部的应力方向。传统方法要先把形变梯度转成应力再求应力在空间上的变化率也就是发散。MLS-MPM 发现如果直接把形变梯度构造成一个仿射矩阵那么动量方程里的应力项可以写成这个矩阵乘以粒子到节点的相对位置。这样一来应力发散的计算就被一个矩阵向量乘替代了。这个替换之所以成立是因为 MLS 近似本身就是在最小化重构误差的意义下用局部仿射函数去拟合速度场。当这个仿射函数被选为形变梯度的函数时动量方程里的应力项自然就变成了仿射矩阵的形式。换句话说不是先算应力再算发散而是把两步合并成一步。对代码来说这意味着粒子到网格的散射循环里每个粒子对每个邻近节点的贡献从“标量权重 × 速度 应力发散项”变成了“标量权重 × 速度 仿射矩阵 × 相对位置”。后者少了一次张量发散计算而且仿射矩阵可以在粒子循环外预先算好。2.3 粒子到网格与网格到粒子两次插值的对称性MLS-MPM 的每一步有两个插值阶段粒子到网格P2G和网格到粒子G2P。P2G 把粒子的质量、动量、仿射矩阵贡献散射到网格节点上网格节点累加后得到节点质量和动量再除以质量得到节点速度。G2P 则用同一套权重函数把网格节点的速度插值回粒子同时更新粒子的形变梯度。这两个阶段用的权重函数必须一致否则动量不守恒。常见的权重函数是二次 B 样条或三次 B 样条它们决定了粒子影响的范围。二次 B 样条每个维度影响 3 个节点三维就是 27 个节点三次 B 样条影响 4 个节点三维 64 个。节点越多精度越高但计算量也越大。MLS-MPM 的巧妙之处在于P2G 阶段散射的仿射矩阵贡献在 G2P 阶段会以对称的方式被重新收集。具体来说P2G 时粒子把仿射矩阵乘以权重再乘以相对位置累加到节点动量上G2P 时节点速度插值回粒子同时粒子的形变梯度更新里也包含节点速度的梯度项。这个梯度项正好对应 P2G 时散射的仿射矩阵。这种对称性保证了每一步的动量守恒和角动量守恒不会因为插值而引入人为的能量漂移。实现时P2G 和 G2P 的循环结构几乎一样只是累加的方向相反。如果你在调试时发现仿真爆炸八成是这两个阶段的权重函数不一致或者相对位置的计算符号搞反了。3. 源码级实现从形变梯度到网格速度的完整链路3.1 粒子数据结构与内存布局在 C 实现里粒子通常用结构体数组AoS或数组结构体SoA来存。AoS 写起来直观但散射到网格时随机访存多缓存命中率低。SoA 把质量、速度、形变梯度分开存P2G 循环里可以按需加载缓存友好。一个典型的粒子结构体长这样struct Particle { float mass; // 质量决定动量交换 float volume; // 初始体积用于计算密度 Vec3 velocity; // 当前速度 Mat3 deformation_gradient; // 形变梯度 F Mat3 affine; // MLS-MPM 的仿射矩阵由 F 和本构模型算出 Vec3 position; // 当前位置 };deformation_gradient是核心状态每一步 G2P 后更新。affine是派生量在 P2G 前根据本构模型算好。常见的本构模型有固定共旋、neo-Hookean、以及各种弹塑性模型。固定共旋最简单适合刚体或低变形弹性体neo-Hookean 能处理大变形但需要做极分解。MLS-MPM 论文里给出的仿射矩阵构造对 neo-Hookean 材料来说是把第一皮奥拉-基尔霍夫应力直接映射成 3x3 矩阵省掉了柯西应力的计算。内存布局上如果粒子数在百万级SoA 的收益很明显。P2G 循环里每个粒子要访问质量、速度、仿射矩阵和位置SoA 可以让这些数组分别连续预取器能更好地工作。但 SoA 的缺点是代码可读性下降而且本构模型计算需要同时访问多个数组。一个折中方案是把频繁一起访问的字段打包成小结构体比如把mass和volume放一起把velocity和position放一起。这样既保留了一定的局部性又不至于让代码太散。3.2 P2G 散射仿射矩阵如何累加到网格P2G 的核心循环是遍历所有粒子对每个粒子找到它影响的网格节点把质量、动量和仿射矩阵贡献累加上去。下面是一个简化版的 P2G 代码假设用二次 B 样条权重三维网格// 网格节点数组每一步开始前清零 std::vectorfloat grid_mass(num_nodes, 0.0f); std::vectorVec3 grid_momentum(num_nodes, Vec3(0,0,0)); for (const auto p : particles) { // 计算粒子在网格上的基准索引和相对位置 Vec3 base p.position / dx - Vec3(0.5f); IVec3 base_idx floor(base); Vec3 fx base - base_idx; // 二次 B 样条权重每个维度 3 个节点 float w[3]; for (int i 0; i 3; i) { w[i] quadratic_bspline(fx[i] - (i - 1)); } // 遍历 3x3x3 个邻近节点 for (int i 0; i 3; i) { for (int j 0; j 3; j) { for (int k 0; k 3; k) { float weight w[i] * w[j] * w[k]; if (weight 0.0f) continue; // 节点索引 int ni base_idx.x i - 1; int nj base_idx.y j - 1; int nk base_idx.z k - 1; int idx ni nj * grid_size.x nk * grid_size.x * grid_size.y; // 节点到粒子的相对位置 Vec3 dpos Vec3(ni, nj, nk) * dx - p.position; // 质量散射 grid_mass[idx] weight * p.mass; // 动量散射速度项 仿射矩阵项 grid_momentum[idx] weight * p.mass * p.velocity; grid_momentum[idx] weight * p.affine * dpos; } } } }这段代码里p.affine * dpos就是 MLS-MPM 替代应力发散的关键项。dpos是节点位置减去粒子位置方向从粒子指向节点。p.affine在粒子循环外根据形变梯度和本构模型算好通常是一个 3x3 矩阵。注意动量散射里速度项和仿射项是分开累加的因为速度项对应平流仿射项对应内部应力。权重函数quadratic_bspline必须满足归一性和对称性否则动量不守恒。二次 B 样条的公式是当 |x| 0.5 时0.75 - x^2当 0.5 |x| 1.5 时0.5 * (1.5 - |x|)^2否则为 0。这个函数在实现时要小心边界尤其是粒子靠近网格边界时节点索引可能越界需要做 clamp 或边界处理。3.3 网格速度求解与 G2P 收集P2G 结束后网格节点上有了质量和动量。节点速度就是动量除以质量但要注意质量为零的节点要跳过否则会除零。然后就是 G2P 阶段把网格速度插值回粒子同时更新形变梯度。G2P 的循环结构和 P2G 几乎一样只是累加方向反过来for (auto p : particles) { Vec3 base p.position / dx - Vec3(0.5f); IVec3 base_idx floor(base); Vec3 fx base - base_idx; float w[3]; for (int i 0; i 3; i) { w[i] quadratic_bspline(fx[i] - (i - 1)); } Vec3 new_velocity(0,0,0); Mat3 velocity_gradient(0.0f); // 速度梯度用于更新形变梯度 for (int i 0; i 3; i) { for (int j 0; j 3; j) { for (int k 0; k 3; k) { float weight w[i] * w[j] * w[k]; if (weight 0.0f) continue; int ni base_idx.x i - 1; int nj base_idx.y j - 1; int nk base_idx.z k - 1; int idx ni nj * grid_size.x nk * grid_size.x * grid_size.y; if (grid_mass[idx] 0.0f) continue; Vec3 node_velocity grid_momentum[idx] / grid_mass[idx]; new_velocity weight * node_velocity; // 速度梯度节点速度乘以权重梯度 Vec3 dpos Vec3(ni, nj, nk) * dx - p.position; velocity_gradient outer_product(node_velocity, dpos) * weight; } } } p.velocity new_velocity; // 更新形变梯度F (I dt * velocity_gradient) * F p.deformation_gradient (Mat3::identity() velocity_gradient * dt) * p.deformation_gradient; }这里velocity_gradient是速度对空间的梯度用来更新形变梯度。outer_product是两个向量的外积得到一个 3x3 矩阵。更新公式F (I dt * velocity_gradient) * F是形变梯度演化的标准形式注意矩阵乘法顺序velocity_gradient在左F在右。如果顺序反了仿真会表现出奇怪的旋转而不是正确的变形。G2P 结束后粒子的位置也要更新p.position p.velocity * dt。然后进入下一步网格清零重新 P2G。3.4 本构模型与仿射矩阵的构造仿射矩阵p.affine的构造取决于本构模型。以固定共旋模型为例它假设材料只发生小变形应力与形变梯度的对称部分成正比。构造过程是先算形变梯度的极分解F R * S其中R是旋转S是对称拉伸。然后仿射矩阵可以写成affine 2 * mu * (F - R) * F^T lambda * log(J) * I其中mu和lambda是拉梅参数J是F的行列式。这个公式看起来复杂但实现时就是几次矩阵乘法和一次极分解。极分解可以用迭代法或特征值分解对于 3x3 矩阵特征值分解的数值稳定性更好。对于 neo-Hookean 材料仿射矩阵的构造类似但应力项不同。MLS-MPM 论文里给出的统一形式是affine stress * F^T其中stress是第一皮奥拉-基尔霍夫应力。不同本构模型的区别只在于stress怎么算。固定共旋的stress是2 * mu * (F - R) lambda * log(J) * F^{-T}neo-Hookean 的stress是mu * (F - F^{-T}) lambda * log(J) * F^{-T}。实现时你可以把本构模型做成一个函数输入形变梯度输出仿射矩阵。这样切换材料只需要换函数不用改 P2G 和 G2P 的循环。4. 避坑与排查MLS-MPM 实现里最容易翻车的几个点4.1 仿真爆炸权重函数不归一或符号错误现象仿真跑几步后粒子速度急剧增大最终变成 NaN。原因通常是权重函数不满足归一性或者 P2G 和 G2P 用了不同的权重函数。二次 B 样条在整数点上的权重和必须为 1如果实现时边界处理不当比如粒子靠近网格边界时节点索引越界被截断权重和就会小于 1导致动量不守恒。解决方法是在 P2G 和 G2P 里用同一个权重函数并且在计算节点索引时做边界检查越界的节点直接跳过但权重和仍然要归一。如果粒子在边界附近可以临时扩大网格范围或者把粒子位置 clamp 到安全区域。4.2 形变梯度更新顺序错误导致旋转异常现象物体在旋转时表现出奇怪的拉伸或剪切而不是刚体旋转。原因通常是形变梯度更新公式里矩阵乘法顺序反了。正确的顺序是F_new (I dt * velocity_gradient) * F_oldvelocity_gradient在左。如果写成F_old * (I dt * velocity_gradient)旋转会被错误地应用到材料坐标系导致非物理变形。解决方法是检查更新公式确保velocity_gradient左乘。另外velocity_gradient的计算里outer_product(node_velocity, dpos)的顺序也不能反node_velocity在前dpos在后。4.3 网格质量为零导致除零现象仿真在某个区域突然出现 NaN或者速度异常大。原因通常是网格节点质量为零但代码里仍然做了momentum / mass。质量为零的节点可能出现在粒子影响范围之外或者粒子散射时权重为零。解决方法是在 G2P 循环里加判断if (grid_mass[idx] 0.0f) continue;。另外P2G 结束后可以遍历网格把质量小于某个阈值的节点标记为无效G2P 时跳过。阈值一般取1e-8左右太小了没意义太大了会丢失细节。4.4 粒子数过多导致缓存失效现象粒子数增加到一定程度后帧率不升反降或者性能提升远低于预期。原因是 P2G 和 G2P 的随机访存模式导致缓存命中率低。解决方法把粒子数据从 AoS 改成 SoA让质量、速度、仿射矩阵分别连续存储。P2G 循环里按需加载预取器能更好地工作。另外网格节点的访问也可以优化比如把网格分成块每个块内的节点连续存储粒子散射时尽量访问同一个块。如果粒子数在百万级以上还可以考虑用空间哈希把粒子分组每组只影响局部网格区域减少随机访存。4.5 时间步长过大导致数值不稳定现象仿真在高速运动或大变形时爆炸但降低速度后正常。原因是显式积分的时间步长超过了 CFL 条件。MLS-MPM 是显式求解器时间步长受网格尺寸和材料波速限制。CFL 条件要求dt dx / c其中c是材料波速通常取sqrt((lambda 2*mu) / rho)。如果dt太大应力波在一个时间步内传播超过一个网格单元数值解会发散。解决方法是根据材料参数和网格尺寸计算最大允许步长或者用子步进每个渲染帧里跑多个小步每步满足 CFL 条件。子步进的代价是计算量增加但稳定性有保证。5. 设计取舍为什么 MLS-MPM 选择牺牲一点精度换实时性MLS-MPM 的数学简化不是没有代价的。传统 MPM 在应力发散的计算上更精确尤其是对于高度非线性的本构模型比如塑性流动和断裂。MLS-MPM 把应力发散吸收进仿射速度场相当于在动量方程里做了一个局部线性化。这个线性化在大多数弹性体和低塑性材料上表现很好但在极端变形下比如材料被拉伸到原长的数倍仿射近似的误差会累积导致仿真结果偏软或偏硬。原始论文里做了大量对比实验结论是在实时应用的精度要求下MLS-MPM 的误差可以接受而性能提升是数量级的。另一个取舍是权重函数的选择。二次 B 样条计算量小但精度低于三次 B 样条。三次 B 样条每个维度影响 4 个节点三维 64 个节点计算量是二次的 2.37 倍但精度提升在视觉上不一定明显。大多数实时实现选二次 B 样条因为它在精度和性能之间平衡得最好。如果你做的是离线渲染或科学计算可以换三次 B 样条但实时应用里二次足够。还有一个取舍是网格分辨率。网格越细粒子影响范围越小计算量越低但粒子到网格的散射会更分散缓存命中率下降。网格越粗每个粒子影响的节点越多计算量上升但散射更集中。常见的做法是让网格尺寸和粒子间距大致相等这样每个粒子影响 3x3x3 个节点计算量和缓存局部性都比较好。如果粒子间距远小于网格尺寸很多粒子会散射到同一组节点造成写冲突如果粒子间距远大于网格尺寸每个粒子影响的节点太多计算量爆炸。我自己的习惯是在实现 MLS-MPM 时先把本构模型和权重函数固定下来用一个小规模场景比如几千个粒子验证动量守恒和能量漂移。如果能量在几百步内漂移超过 5%说明权重函数或更新公式有问题。验证通过后再逐步增加粒子数观察帧率和缓存命中率的变化。从那以后我每次调参都会先跑一遍守恒性检查再去看视觉效果。希望帮到你。本文还有配套的精品资源点击获取