Efficient GPU Cloth Simulation with Non-distance Barriers and Subspace Reuse

一句话总结

本文提出一套面向高分辨率布料的 GPU 仿真框架，用”非距离障碍(NDB)+ 部分 CCD”替换掉绝大部分昂贵的连续碰撞检测，用”子空间复用”消除迭代求解器难以处理的低频误差，再配合”残差前推”缓解时间预算受限时的阻尼伪影，整体比现有 GPU 布料仿真快约一个数量级，同时保证全程无穿透。

研究背景

高分辨率布料能刻画出丰富的褶皱与折痕，但自由度(DOF)激增会严重拖累仿真性能，而交互式应用(游戏、虚拟试衣、走秀设计)对每步的时间预算有硬性上限。

当前主流的无穿透方案受 IPC(Incremental Potential Contact)启发：把不等式约束转成障碍势能，在一对相互靠近的图元之间引入随距离减小而无限增强的斥力。这类方法稳健，但存在两个核心瓶颈：

• 碰撞处理昂贵：基于距离的障碍(本文称之为 DBB, distance-based barrier)依赖图元间的真实距离 \(d_i\)，需要在每次非线性迭代里反复做连续碰撞检测(CCD)求解首次碰撞时刻(TOI)，对稠密网格几乎所有三角形都可能接触，CCD 成为整条流水线的主导开销。
• 系统求解低频误差难消：GPU 上常用的 Jacobi 类迭代求解器擅长平滑高频误差，却对低频残差束手无策；而模型降阶/多重网格虽能处理低频，但系统矩阵随布料姿态非线性变化，逐姿态重建子空间不现实。

本文围绕这两点做深度优化。

方法

整体框架

方法建立在 Projective Dynamics(PD)之上。隐式积分下每步求解变分问题：

\[\arg\min_{x} E = I(x,\dot{x}) + \Psi(x) + B(x)\]

其中惯性项 \(I = \tfrac{1}{2h^2}\lVert M^{1/2}(x-z)\rVert^2\)，\(\Psi\) 为弹性势能，\(B\) 为碰撞障碍势。PD 把它拆成可高度并行的局部步(逐约束投影目标位置 \(y_i\))与全局步线性求解：

\[\left(\frac{M}{h^2} + \sum_i w_i S_i^\top A_i^\top A_i S_i\right) x = \frac{M}{h^2} z + \sum_i w_i S_i^\top A_i^\top B_i y_i\]

其中权重 \(w_i\) 表示约束的优先级。本文的三项核心贡献分别改造了碰撞障碍(NDB + 部分 CCD)、全局求解(子空间复用)、以及时间预算受限时的收尾(残差前推)。

flowchart TD
    A[预计算: 静止姿态全局矩阵 H 的特征分解<br/>大子空间 U(r̄维) / 复用子空间 V(r=30维)] --> B[Warm start: 在 U 子空间求解 I+Ψ, 3~5 次 LG 迭代]
    B --> C[Full CCD 建立初始碰撞约束 B]
    C --> D[外循环]
    D --> E[内循环: 局部投影 Ψ 与 B]
    E --> F[全局步: V 子空间复用求解 + Chebyshev A-Jacobi]
    F --> G[Partial CCD 更新 NDB 权重 w_i = k·K^{a_i}]
    G --> E
    E --> H[Full CCD + line search filtering 保证无穿透]
    H --> D
    D --> I{退出 TOI 过小?}
    I -->|是| J[Residual Forwarding: 估计虚拟力 δf 传给下一步]
    I -->|否| K[下一时间步]
    J --> K

关键设计 1：非距离障碍 NDB

放弃对真实距离 \(d_i\) 的依赖，改用碰撞事件的”生命周期”来自适应加权。若某接触约束 \(C_i\) 连续 \(a_i\) 次 LG 迭代都处于激活状态，则其权重按指数增长：

\[w_i = B_i^{NDB} = k\,K^{a_i}\]

其中 \(k\) 为初始权重，\(K\) 为增长率底数(实验取 \(K=2\))。当碰撞长时间不被解决，\(w_i \to \infty\)，全局步就会以最高优先级满足该约束的无穿透目标位置。NDB 不再是严格意义上的障碍函数(不显式依赖 \(x\))，但同样能把变量约束在可行域内。代价是：迭代次数是离散计数，难以精细区分邻近顶点的斥力，需要更多迭代。

关键设计 2：部分 CCD（Partial CCD）

由于 NDB 只需判断碰撞是否”仍然激活”(布尔值)，不需要精确 TOI，于是把 CCD 的三次方程求根问题降级为一系列点积计算。基于碰撞的必要条件，构造查询函数：

\[Q(\lambda) = \big(p_2^1(\lambda) - p_1^1(\lambda)\big)\cdot\big(p_2^0(\lambda) - p_1^0(\lambda)\big)\]

只要在参数域 \(\Omega_\lambda\) 上采样点足够密、且所有采样点的 \(Q(\lambda)>0\)，即可判定无碰撞。论文用拉格朗日余项定理推出充分采样间隔的理论上界：

\[\rho^* = \frac{\left(\tfrac{1}{\alpha}-1\right)(H_0)^2}{2\sqrt{2}\,L(H_1+2L)}\]

实测中在 RTX 3090 上，部分 CCD 比原始 IPC 求解三次方程快 6000 倍以上，比先进多项式求解器快约 2000 倍；对 5 万 DOF 以上的高分辨率网格，往往只需 3 个采样点就工作良好。部分 CCD 作为 warm start，仍配合低频次的完整 CCD 兜底保证无穿透。

关键设计 3：子空间复用（Subspace Reuse）

核心观察：低频应变传播与碰撞诱导的高频形变近乎正交，且静止姿态构建的低频模态描述的是全局应变分布，对高频碰撞形变不敏感。因此可以复用静止姿态全局矩阵 \(H\) 的最小 \(r\) 个特征向量组成子空间 \(U\)，即便碰撞把矩阵改成 \(H+\Delta H\) 也照样用：

\[U^\top(H+\Delta H)U\,q = U^\top(b - Hx - \Delta H x)\]

先在子空间内解掉低频误差，得到理想初值 \(\tilde{x}=X+Uq\)，再用聚合 Jacobi(A-Jacobi)迭代平滑剩余的高频残差。关键工程技巧是把局部投影中碰撞约束的 \(A_i=B_i=\mathrm{Id}\)，使 \(\Delta H\) 成为对角矩阵，于是子空间更新可写成秩一矩阵的加权和：

\[U^\top \Delta H\, U = \sum_{j=1}^{k} w_j\, U_{c_j}\otimes U_{c_j} \in \mathbb{R}^{r\times r}\]

其中每个 \(U_{c_j}\otimes U_{c_j}\) 可预计算。这使得子空间矩阵更新比用 cuBLAS 直接算 \(U^\top\Delta H U\) 快 50~150 倍，整个子空间求解(取 \(r=30\))耗时不到 0.1 ms。实测子空间复用平均减少后续 A-Jacobi 迭代约 70%。实现上用两个子空间：大维度 \(\bar{r}\) 做无碰撞 warm start，小维度 \(r=30\) 做碰撞在环的复用。

关键设计 4：残差前推（Residual Forwarding）

交互应用有硬性时间预算，未收敛就被迫结束会因小步长 line search 过滤导致过阻尼/锁死。RF 把当前步未解决的残差力估计为一个”幽灵外力” \(\delta f\) 传给下一步：

\[z = x^* + h\dot{x}^* + h^2 M^{-1}(f_{ext}+\delta f)\]

最优 \(\delta f\) 通过最小化 \(\lVert \tfrac{h^2}{2}M^{-1}\delta f - (\nabla^2 E(x^*))^{-1}f_r\rVert\) 得到，其中 \((\nabla^2 E)^{-1}f_r\) 等价于再多做一次牛顿求解，可用子空间复用高效近似。RF 跳过 line search，允许非碰撞顶点在动量和弹性下继续运动，从而缓解阻尼锁死。

实验结果

平台为 i7-12700 CPU + RTX 3090 GPU。下表摘录 Table 1 中若干代表性场景(时间单位 ms，FPS 为整体帧率)：

场景	物体数	DOF	单元数	子空间 \(\bar{r}\)/\(r\)	LG 时间(次数)	FPS
扭转布料	1	121K	80K	120/30	47.4 (157)	12.4
动物过河	44	678K	1.1M	120/30	103.2 (47)	7.0
打结	1	310K	203K	120/30	63.6 (74)	11.8
茶壶堆叠	11	500K	338K	120/30	35.9 (31)	16.6
布料搅拌机	8	820K	540K	120/30	160.7 (73)	5.2
覆盖船体	5	930K	941K	120/30	117.8 (51)	7.2
走秀(封面)	2	1.1M	656K	120/30	167.3 (67)	4.8

对比关键结论(数字忠于原文)：

• vs PD-IPC[Lan et al. 2022b]：在异质材料的”动物过河”场景中快 12×；NDB 用部分 CCD 更新障碍权重，减少 60% 以上计算时间；相比 DBB 在扭转测试中节省 25%~30% 的 LG 迭代。
• vs PD-BFGS[Li et al. 2023]：”打结”场景快 130×(PD-BFGS 每帧需数秒)，因为其 B-spline 子空间维度高达约 10000，而本文预计算子空间更新只需 \(r=30\)。
• vs PD-Coulomb[Ly et al. 2020]：”球上悬垂”场景快约 30×(25 FPS)，但差距主要来自 PD-Coulomb 是 CPU 实现；两者方法正交。
• 残差前推：茶壶落布场景每步 50 次迭代时结果几乎与完全收敛的 Newton-IPC 一致，且比 [Li et al. 2021a] 快 500×;迭代数降到 5 次时高频细节丢失、出现橡胶状伪影。

亮点与局限

亮点：

• 把碰撞检测从”解三次方程求 TOI”降级成点积采样(部分 CCD)，配合 NDB 的自适应加权，从根本上砍掉了布料仿真最大的开销。
• 子空间复用巧妙利用”低频模态对高频碰撞不敏感”这一观察，用静止姿态的特征向量应对变化的系统矩阵，并通过秩一分解把更新预计算，让子空间求解 <0.1 ms。
• 兼顾无穿透保证、高分辨率(百万级 DOF)与交互帧率，整体性能领先现有方法约一个数量级。

局限：

• 基于 PD，弹性模型 \(\Psi\) 必须是二次型，难以直接纳入更复杂的同质化或数据驱动织物模型。
• 残差前推是启发式手段，当残差误差过大时会失效，甚至引入不自然的抖动/凸起伪影或导致仿真失败。
• 作为非牛顿方法，对时间步长和终止判据比 projection-Newton 更敏感(见 Table 2，过大步长/过松判据会出现 Jittery 或 Fail)。

延伸思考

• NDB 用”碰撞生命周期”替代”真实距离”，本质是把连续几何量替换成离散迭代计数——这在收敛精度上是妥协(难以精细区分邻近图元)，但换来了巨大的计算收益。这种”用迭代动力学隐式编码约束强度”的思路，是否能推广到摩擦、接触之外的其他不等式约束求解？
• 子空间复用可以理解为全局步的两级多重网格：静止姿态子空间消低频、A-Jacobi 消高频。论文提到该策略对低频运动同样困扰迭代求解器的体积可变形体也高度有效，说明它并非布料专属。
• 局限中提到可与块下降法(block descent)结合以支持通用材料，这暗示”二次型 PD + 子空间复用”与”通用非线性材料”之间存在可衔接的桥梁，是后续通用化的自然方向。
• 部分 CCD 的采样密度理论上界依赖多个动态变量，实践中却用极少采样点(3 个)就够——理论保守性与工程可行性之间的这种落差，值得在其它几何查询问题中借鉴。