本文整理了布料模拟实践过程中学习到的各种理论方法，记录了实践过程中遇到的问题和解决办法。主要基于Unity完成了基于显示欧拉法和PBD方法的布料模拟，并使用ComputeShader对模拟过程的加速。

实现思路

思路概述

计算机图形学中的物理模拟方法可以笼统概括为三个阶段：初始化阶段，物理模拟阶段，更新数据阶段。

初始化阶段中需要构建物理模型。在开始模拟前需要先理解一块布的结构，一块布中存在结构力、剪切力以及弯曲力（图），根据这一结构特点可构建质点弹簧模型，据此模型计算质点间的相互作用完成模拟。

物理模拟阶段的具体实现方法很多，分为两个主要流派，分别是基于物理的模拟和非物理的模拟。

物理模拟方法以欧拉法为主。在质点弹簧系统中，根据胡克定律计算质点间的力，更新质点加速度、速度、位置。欧拉法具体有显示欧拉法、半隐式欧拉法和隐式欧拉法，区别在于显示欧拉法使用t时刻受到的合力更新t+1时刻的位置，隐式欧拉法使用t+1时刻受到的合力更新t+1时刻的位置，半隐式欧拉法结合显示欧拉法和隐式欧拉法，通过t时刻受到的合力更新更新t+1时刻的速度，使用该速度更新t+1时刻的位置。显示欧拉的方程是线性的，隐式欧拉的方程是非线性的。通常来说求解非线性方程，可以使用数值分析，比如将非线性方程转换后可以用牛顿迭代的方法求解最小值，或者使用龙格库塔方法进行求解。从稳定性角度来讲，显示欧拉法是条件稳定的，在时间步长较小时可以保持稳定，但当时间步长增大误差就会累积从而失去稳定性。此外随着模拟时间增加，误差也会累积使得出现模拟错误。隐式欧拉法理论上是无条件稳定的，但由于通常使用线性方式近似求解所以也是条件稳定的，仍需控制时间步长进行才能保证模拟稳定运行。

非物理模拟方法中常见的有基于位置的动力学模拟方法（Position-Based Dynamics，PBD）。主要思路是通过欧拉法预测下一时刻位置，再施加约束进行修正，PBD方法不仅可以用到布料模拟中，而且可以用到软体模拟等模拟中。作者给出了PBD方法的主要思路和框架，只需要确定约束方程，并将约束方程加入到约束列表中解约束即可。从稳定性方面看，与欧拉法相比PBD方法更为稳定，但也不是无条件稳定的，其稳定性也受迭代次数、时间步长的影响。

数据更新阶段需要将更新好的位置、速度重新赋给模型顶点，然后绘制图形。

质点弹簧系统的构建

结构化网格

结构化的网格手动计算每个顶点的初始位置，构建三角形，以及三角形索引。然后为每个三角形每一条边边构建一个弹簧。在此之前，为了防止重复构建弹簧，需要剔除重复边。

非结构化网格

非结构化的网格直接读取Mesh中的顶点、三角形索引等数据。再为每个三角形的每一条边构建一条弹簧。同样地为了防止重复构建弹簧，需要剔除重复边。非结构化网格可以用于任意网格的模拟。

物理模拟过程

实践中实现了显示欧拉和PBD方法的物理模拟。下面给出伪代码。

显示欧拉

foreall vertices i
	velocityt+1 = velocityt + dt*at;
 	positiont+1 = positiont + dt*velocityt;
endfor

PBD方法的物理模拟

PBD方法的模拟流程如下：

高斯-赛德尔迭代法（Gauss-Seidel Iteration）

高斯-赛德尔迭代法每次迭代中立即使用求解得到新值进行计算。在模拟中体现为遍历每一条边，通过调整每条边两顶点位置，以及边所属的正方格实现对预测位置的修正。

问题：如果一个正方形网格中只构建一条剪切弹簧，那么会出现布料向一侧变形的错误结果。这是由于Gauss-Seidel迭代以弹簧顺序进行，调整下一根弹簧后上一根弹簧调整过的顶点也会发生移动。

雅可比迭代法（Jacobi Iteration）

雅可比迭代法又称同时迭代法，通过将线性方程的对角矩阵变换为对角矩阵及其他部分，可在每次迭代中直接求得解。体现在模拟过程中为不再需要预测位置而是直接通过解方程得到模拟位置。雅可比迭代法相比高斯赛德尔方法对并行计算更加友好，且不容易发生布料形变的问题。

实现方法

数据结构

显示欧拉

Gauss-Seidel

Jacobi

布料构建

布料模拟

欧拉法

计算位置
计算速度
碰撞响应

PBD Gauss-Seidel

距离约束
弯曲约束
碰撞约束

PBD Jacobi

约束方法
碰撞响应

GPU模拟

计算着色器（Compute Shader）

DirectX 11、OpenGL 4.3后新增了计算着色器。GPU除了可以进行传统光栅化管线还可以利用GPU的架构进行并行运算，完成后处理效果、粒子系统等功能。Compute Shader是独立于Vertex Shader和Fragment Shader的管线。

写入冲突与原子操作

在调用计算着色器中的函数需要在CPU端进行，使用Dispatch()调用计算着色器函数，用户可以指定运行的线程组数量。一个线程组包含多个线程，运行时GPU上的多个线程组就会同时无序运行。如果两个弹簧结构中有同一个顶点的索引，在向该顶点写入数据时可能会有多个线程向同一地址的数据写入的情况，使得写入失败。

对于这种情况可以考虑原子操作，在CS中对于int和uint类型的数据可以使用一系列InterLock()函数避免该情况，但因为位置、速度等都是三维浮点数就需要做出调整。这里参考了别人的解决方法，大概是将要写入的数据转为uint类型，计入临时Buffer中，最后更新数据阶段再从Buffer中读取。

实现框架

基于Unity实现，渲染部分由Unity绘制。C#脚本负责初始化数据及设置计算着色器的Buffer，计算着色器负责处理并行

效果显示

测试设备

CPU：12th Gen Intel® i9-12900HX
GPU: NVIDIA GeForce RTX 3080 Ti

效果演示

请移步B站： https://www.bilibili.com/video/BV1qabLeME4p/

未完成的工作

未完成布料的自相交

虽然PBD文献中给出了布料自相交的约束公式，但每个顶点需要和每个三角形进行碰撞检测放在GPU上就不好实现。看到了一些方法比如可以用BVH树加快这一过程，但在GPU中构建树又是难度较高的事。因此这部分还没来得及完成。

未完成在GPU中直接绘制布料

用户可以调用DrawProceduralIndirect()在GPU程序化绘制几何图形，以此节约从GPU读取Buffer数据到CPU这一过程的时间。个人一直掌握不好这个函数的用法，常常绘制出不合预期的图形或根本绘制不出图形。故还没有应用到项目中。要不然在GPU端的模拟帧率还可以继续提升。

参考文献及资料

[1] Matthias Müller, Heidelberger B, Hennix M, etal. Position Based Dynamics[J]. VRIPHYS, 2006(3).
[2] WU L, WU B, YANG Y, etal. A Safe and Fast Repulsion Method for GPU-based Cloth Self Collisions[J]. ACM, 2020, 1(1).
[3]GAMES103
[4] 迭代法求解线性方程（简单迭代法雅可比(Jacobi)迭代法高斯-塞德尔(Gauss-Seidel)迭代法逐次超松弛(SOR)迭代法）数值计算方法_高斯-雅可比-CSDN博客
[5] PBD（Position Based Dynamics）学习笔记_pbd算法-CSDN博客
[6] wlgys8/GPUClothSimulationLearn: Unity GPU布料物理模拟入门 (github.com)
[7] wlgys8/PBDClothLearn: Cloth Simulation by Position Based Dynamics + Unity Job System (github.com)
[8] Solicey/UnityPBDClothSimulationOnGPU: Cloth simulation based on Position Based Dynamics (PBD) in Unity, using compute shaders (github.com)