企业官网建设流程全解析

简介：本资源为基于移动最小二乘物质点法（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::vector<float> grid_mass(num_nodes, 0.0f); std::vector<Vec3> 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_old，velocity_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%，说明权重函数或更新公式有问题。验证通过后再逐步增加粒子数，观察帧率和缓存命中率的变化。从那以后，我每次调参都会先跑一遍守恒性检查，再去看视觉效果。希望帮到你。

本文还有配套的精品资源，点击获取