基于Position Based Dynamics的布料模拟

之前做过一版简单GPU布料模拟 - GPU布料模拟入门。里面使用的方式是传统的基于力的算法，即:

根据位置约束计算出质点受力
根据受力计算出加速度
利用加速度迭代速度
利用速度迭代位置

如此往复，但这种模拟方式需要比较小的迭代步长才能保持稳定，比较耗费计算量。因此这次准备使用Position Based Dynamics算法来进行模拟，然后使用Unity的Job System版本实现CPU并行计算。 Position Based Dynamics算法的优势在于，其可以拥有很大的迭代步长，且有优秀的稳定性。而JobSystem配合Burst编译，其运算性能也是非常强劲的。

目前实现的功能有:

外力
- 重力
- 风力
内部约束
- 距离约束
- 弯曲约束
碰撞约束
- 与Sphere、Box、Capsule的碰撞
固定约束
- 可以将布料固定到任意位置
对场景物体施加反作用力

暂未实现的功能

布料自碰撞
布料破坏
GPU模拟版本
连续碰撞检测

先放个效果图吧

1. Position Based Dynamics简介

PBD的算法流程图如下

图中x为当前质点位置,p为预测质点位置。可以看到:

算法在每一帧开始时，会根据外力和阻尼先于计算一个位置(p)
根据当前位置(x)和预测位置(p)，去进行碰撞检测，如果发现有碰撞，就生成一个碰撞约束
迭代求解所有约束，最终会得到新预测位置p
利用p - x 来计算速度v，并将p赋值给x
最后再计算一遍速度(处理碰撞反弹之类速度变换)

关于外部力的定义:

像重力、风、可以视为外部力，但是像质点之间的相互作用我们通常视为内部约束，而非力。同样外部物体对质点的碰撞，也视为碰撞约束而不是外部力。

2. 布料模型构建

在物理上，我们通常可以将布料视作由一个个质点组成的结构。质点之间满足一定的约束条件，来模拟布料内部的作用力，如下图所示:

在传统的基于力的模拟中，质点之间的约束通常是用胡克定律以弹簧振子的形式去模拟。但是在基于位置的方案中，我们将会采用另外一种约束公式。

那么具体到实际渲染中，质点代表的是什么呢？

一种方案是，让质点代表骨骼，然后布料的Mesh顶点以权重的方式绑定到骨骼上。这种方案的优势是质点的数量和布料顶点的数量相互独立，也即物理模拟的精度和渲染的精度相互独立。但劣势在于需要手动对布料做骨骼绑定。

另一种方案是直接让质点与布料模型顶点一一对应，然后让三角形的边作为顶点之间的约束。好处是无需额外的骨骼绑定工作，我们可以针对任意Mesh进行布料模拟(也没那么任意，还是需要满足一定的拓扑条件的)，坏处是如果模型顶点较多，将比较耗费性能。

这里作为Demo项目，将直接采用第二种方案。

想要进行布料模拟的Mesh，必须满足以下条件:

一条边只被两个三角形共享
不封闭

3. 数据结构

因为要使用Unity JobSystem来进行并行计算，所以我们将采用平坦的结构来保存位置、速度之类的信息。

这里我们将使用NativeList来保存所有质点的位置、速度、预测位置、质量等信息。位置的初始信息可以直接从Mesh中读取，但质点的质量需要通过一定算法生成。

private NativeList<float3> _vertices;
private NativeList<float3> _velocities;
private NativeList<float3> _predictPositions;
private NativeList<float> _masses;

另外这里还引入了Unity.Mathematics这个库，里面定义了很多与shader类似的数据结构，例如float4这种。这个数学库是Unity专门针对SIMD进行优化的。

3.1 质量生成算法

通常来说，我们可以给一块布指定其密度(kg/m^2)。这样的话，我们就可以根据Mesh的面积来计算质量。具体来说就是:

遍历Mesh的所有三角面
计算每个三角面的面积，乘以密度得到质量
将质量三等份，分别累加给每个顶点

这样我们就得到了每个顶点的质量。

代码如下:

//var indices;
//var vertices;
//...
_masses = new NativeList<float>(this.vertexCount,Allocator.Persistent);
_masses.Resize(this.vertexCount,NativeArrayOptions.ClearMemory);
for(var i = 0; i < indices.Length / 3; i++){
    var offset = i * 3;
    var i0 = indices[offset];
    var i1 = indices[offset + 1];
    var i2 = indices[offset + 2];
    var v0 = vertices[i0];
    var v1 = vertices[i1];
    var v2 = vertices[i2];
    var area = IntersectUtil.GetArea(v0,v1,v2);
    var m = area * _setting.density;
    var m3 = m / 3;
    _masses[i0] += m3;
    _masses[i1] += m3;
    _masses[i2] += m3;
}

4. 预测位置计算

为了使用JobSystem进行并行计算，我们可以实现一个如下的结构:

struct PositionEstimateJob : IJobParallelFor
{
    [ReadOnly] 
    public NativeArray<float3> positions;
    [ReadOnly] 
    public NativeArray<float3> velocities;
    [ReadOnly]
    public NativeArray<float3> normals;
    [ReadOnly]
    public NativeArray<float> masses;
    [WriteOnly]
    public NativeArray<float3> predictPositions;
    public float3 fieldForce;
    public float damper;
    public float dt;

    public void Execute(int index)
    {
      var p = positions[index];
      var v = velocities[index];
      var m = masses[index];
      //....
    }
}

JobSystem将会在多线程上并发执行Execute函数，在这里将会针对每个质点执行一次Execute函数，index为当前执行的任务索引(也即质点索引)。由于在这一步每个质点只更新自己的预测位置，不存在Race Condition，因此可以进行无锁并发。

在这里，我们输入的数据有:

positions - 位置
velocities - 速度
normals - 法线
masses - 质量
fieldForce - 外部场力,例如风力
damper - 阻尼系数
dt - 迭代时间

输出为:

predictPositions

位置预测公式如下:

$$ \begin{aligned} \vec{v_1} = \vec{v_0} + (\vec{g} + \frac{\vec{F_{e}}}{m}) \Delta t \newline v_1 = v_1 * \max(1 - \frac{k_d}{m} \Delta t,0) \newline p_1 = p_0 + v_1 * \Delta t \end{aligned} $$

这三条公式，第一步为预测速度，第二步对速度施加阻尼，其中kd为阻尼系数，第三步利用速度预测位置。

g为重力加速度，$F_e$为作用于质点的外力，在本例里我们只用到了风力。风对质点的作用力，需要投影到质点的法线方向，这也是为什么这个Job的输入需要有法线信息。

随风舞动的效果：

5. 碰撞计算

这里我们将实现针对球体(Sphere)、方体(Box)、胶囊(Capsule)三种基础类型的碰撞检测。为简单起见，目前将采用离散的检测方式，即判断质点是否在相关碰撞体内部。

算法流程为:

for position, predictPosition,index in vertices:
  for collider in colliders:
    concat = CheckCollision(position,predictPosition,collider)
    if concat:
      AddCollisionConstraint(index,concat.position,concat.normal)

可以看到，以上伪代码的具体思路为:

遍历所有质点，取其位置和预测位置，与所有的碰撞体做检测
如果发现进入碰撞体内部，则计算出碰撞体表面距离当前质点最近的一个点
对当前质点施加一个碰撞约束，让其在后面的约束迭代中强行离开碰撞体，回到表面去。

5.1 几种基础碰撞体的数据结构

首先是球体，这是最简单的一种碰撞体

public struct SphereDesc{
    public float3 center;
    public float radius;
}

然后是平行六面体:

public struct BoxDesc{
    public float3 min;
    //3个边轴,xyz为normalzied的朝向，w为长度
    public float4 ax; 
    public float4 ay;
    public float4 az;
}

可以看到我们使用了xyz-min的角、和3个边向量来定义这个Box，严格来说，它能描述三维空间中的任意平行六面体。其中三个边向量，我们用float4保存，并且做了一些预计算。其xyz分量为单位向量，w分量向量长度。相关的预计算可以加速后面的碰撞检测。

胶囊体:

public struct CapsuleDesc{
    public float3 c0;
    public float4 axis; //xyz为单位化的方向，w为长度
    public float radius;
}

Capsule我们可以将其看成两个半球 + 一个圆柱体的组合，如下图所示:

以上的结构数据均定义在世界坐标空间，关于如何判定一个点是否在这些碰撞体的内部，这里将不详细展开，不然又是很长篇。总之，在IntersectUtil这个类里，我们实现了相关的函数，以球体为例

public static bool GetClosestSurfacePoint(float3 p, SphereDesc sphere, out ConcatInfo concatInfo)

GetClosestSurfacePoint它接收一个点p坐标，和一个碰撞体描述(这里是SphereDesc)，如果发现点在碰撞体内部，就返回true,并且输出一个ConcatInfo结构。

ConcatInfo里保存了碰撞体表面距离p最近一个点的position和normal信息。

我们的目的不仅仅是让这些碰撞体来约束布料，同时布料也要反作用于这些碰撞体(假如是RigidBody的话)。

因此这里额外定义一个刚体描述结构:

public struct RigidbodyDesc{
    public float mass;
    public float bounciness;
    public float3 velocity;
}

然后将RigidBody和Collider组合在一起:

public struct RigidbodyColliderDesc<T> where T:IColliderDesc{
    public int entityId;
    public T collider;
    public RigidbodyDesc rigidbody;
}

其中entityId后面会用来查找Unity对象。

然后定义个结构把所有需要与布料进行碰撞检测的数据整合在一起:

public struct CollidersGroup{
    private NativeList<RigidbodyColliderDesc<SphereDesc>> _spheres;
    private NativeList<RigidbodyColliderDesc<BoxDesc>> _boxes;
    private NativeList<RigidbodyColliderDesc<CapsuleDesc>> _capsules;
}

5.2 CollisionJob

接下来我们要用多线程任务并行执行碰撞检测和反作用力计算。

看一下这个Job的输入输出数据结构:


public struct CollisionJob : IJobParallelFor
{
    [ReadOnly]
    public float dt;
    [ReadOnly]
    public NativeArray<float3> positions;
    [ReadOnly]
    public NativeArray<float3> predictPositions;
    [ReadOnly]
    public NativeArray<float> masses;
    [ReadOnly]
    public CollidersGroup collidersGroup;
    [WriteOnly]
    public NativeArray<CollisionConstraintInfo> collisionConstraints;
    public NativeArray<ConstraintType> constraintTypes;
    public NativeList<RigidBodyForceApply>.ParallelWriter rigidBodyForceApplies;
}

它会根据当前质点的predictPosition和position信息，与CollidersGroup中的所有碰撞体进行碰撞检测，如果检测到有碰撞发生，则会做如下事情:

修正predictPosition，避免碰撞发生
将双方速度投影到接触点法线方向，利用动量守恒、动能公式重新计算双方在法线方向的速度
将对Collider的速度修正贡献写入到rigidBodyForceApplies中，供后续实际作用到unity rigidboy对象。
将对质点的约束信息写入到collisionConstraints中，供后续约束迭代阶段使用.

那么在法线方向上，碰撞引起的速度改变要怎么计算呢？这其实是高中的物理知识。两个物体碰撞，要符合动量守恒定律，即:

$m_0v_0 + m_1v_1 = m_0v_0' + m_1v_1'$。

当两个物体是完全弹性碰撞时，他们在动能上要符合以下公式(即动能守恒):

$m_0v_0^2+m_1v_1^2=m_0v_0'^2+m_1v_1'^2$

而如果是完全非弹性碰撞，则动能不守恒，但最终两个物体的速度一致，即:

$v_0' == v_1'$

那么介于这两者之间的情况，我们通常用一个弹性系数C来描述。C = 1表示完全弹性碰撞，C = 0表示完全非弹性碰撞，完全的碰撞速度计算公式如下:

$$ \begin{aligned} v_0' = \frac{(c + 1)m_1v_1 + v_0(m_0-cm_1)}{m_0+m_1} \newline v_1' = \frac{(c + 1)m_0v_0 + v_0(m_1-cm_0)}{m_0+m_1} \end{aligned} $$

总之，利用以上的公式，我们就成功的完成了法线方向的碰撞反馈计算。

实际上，如果两个物体表面存在摩擦力，那么在垂直于法线的另一个方向上也会产生力的相互作用。但是在当前的这个版本里，我们暂且忽略一下吧。

6. 内部约束

在布料模拟中，存在如下几种约束:

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

下面依次来说明.

6.1 距离约束

距离约束就是两个质点之间的距离，必须保持被约束在一定范围之内。

通常在布料初始化的时候，会记录下存在距离约束关系的质点之间的距离，称为RestLength(也即图中的d)。在后续模拟中，当两个质点间距离小于RestLength时，我们就要将其推开，反之则要将他们拉拢。也即他们之间存在一根无形的弹簧。

那么两个质点之间的距离约束，要符合怎样的迭代公式呢？我们要如何为p1和p2分别计算出它们的修正量delta_p1和delta_p2呢？

这里直接给出结论公式:

$$ \begin{aligned} \Delta \vec{p_1} = - \frac{m_2}{m_1 + m_2}(|p_1 - p_2| - d) \frac{p_1 - p_2}{|p_1-p_2|} \newline \Delta \vec{p_2} = \frac{m_1}{m_1 + m_2}(|p_1 - p_2| - d) \frac{p_1 - p_2}{|p_1-p_2|} \newline \end{aligned} $$

原论文[1] Matthias Müller, Position Based Dynamics, 2006 中给出了更一般的约束关系的公式推导，有兴趣可以看看，主要是用了牛顿-拉夫逊迭代法。

在本例中，我们会根据三角面信息，为所有存在边连接的点建立距离约束。即:

private void BuildDistConstraints(){
    var edges = _meshModifier.edges;
    foreach(var e in edges){
        this.AddDistanceConstraint(e.vIndex0,e.vIndex1);
    }
}

距离约束结构定义:

public struct DistanceConstraintInfo{
    public float restLength;
    public int vIndex0;
    public int vIndex1;
}

6.2 弯曲约束

弯曲约束可以用来控制相邻两个面片的对折程度，这符合现实中布料的特性。例如当我们将一块布料对折放置在平面上，其对折处会有一段弧形隆起。这就是布料内部结构存在的抵抗弯曲的力。下图展示了一个弯曲约束的效果:

在布料建模中，弯曲约束可以由如下一幅图定义，这里我直接拿论文中的图:

两个相邻的三角面，共4个顶点，组成一个约束公式。

两个三角面的法线点乘，给出了三角面的夹角信息，在后续的模拟迭代中，我们只要约束这个夹角就行了。论文中同样给出了完整的计算公式推导过程，这里直接截图贴上:

细细品味吧。

弯曲约束的结构定义:

public struct BendConstraintInfo{
    public int vIndex0;
    public int vIndex1;
    public int vIndex2;
    public int vIndex3;
    public float rest; // angle
}

6.3 固定约束

顾名思义，这种约束是将某个质点固定在某个坐标上，令其无法移动。所以这个约束的实现很简单，我们只要让这个质点位置不变，速度为0即可。

6.4 碰撞约束

在Part5中，我们提到了碰撞计算。但是在碰撞计算阶段并不会去直接修改质点的位置，而是生成一个碰撞约束，交由后续约束迭代阶段进行质点位置计算。

目前碰撞约束的结构很简单:

public struct CollisionConstraintInfo{
    public float3 concatPosition;
    public float3 normal;
    public float3 velocity;
}

在约束迭代阶段，质点会将自己的位置尽量往concatPosition去靠近。碰撞约束作为最高级约束，在其生效时，其他的约束应该不起作用。 velocity保存了碰撞计算得到的质点速度。

6.5 约束迭代步骤

考虑到Unity JobSystem的特性，目前实现的约束迭代流程如下:

Distance Constraint J

PBDClothLearn

Install / Use

README