背景：在计算机图形学领域，主流动力学仿真方法是基于力的(force-based) 。累积内力和外力，根据牛顿第二运动定律计算加速度。 并运用时间积分的方法获得物体的速度和位置。
方法：本文提出了一种忽略速度层的方法，直接作用于位置。
优势：具有良好的可控性。可以避免基于力的系统中显式积分方案的过冲问题；此外，可以轻松处理碰撞约束，并且可以通过将点投影到有效位置来完全解决穿透问题。
主要特点：
1） 基于位置的仿真可以控制显式积分并消除典型的不稳定性问题。
2） 在模拟过程中可以直接操纵顶点和对象部分的位置。
3） 我们提出的公式允许在基于位置的设置中处理一般约束 。
4） 基于显式位置的求解器易于理解和实现 。
讲解：* 为核心步骤
第（1）-（3）行：初始化系统
第（5）行：
f
e
x
t
\textbf{f}_{ext}
fext表示无法转换成位置约束的力，我们仅用他来表征重力，因此改写为
v
i
←
v
i
+
Δ
t
g
\mathbf{v}_{i} \leftarrow \mathbf{v}_{i}+\Delta t \mathbf{g}
vi←vi+Δtg
第（6）行：可以对速度进行阻尼处理
第（7）行*：新位置的估计坐标
p
i
\bm{p}_i
pi由显式欧拉积分步骤获得。
第（8）行：除了在仿真过程中从始至终不变的给定约束
C
1
,
…
,
C
M
C_1, \dots , C_M
C1,…,CM，还生成了
M
c
o
l
l
M_{coll}
Mcoll个时变的碰撞约束
第（9）-（11）行*：同时考虑了固定的和碰撞的约束，并对
p
i
\bm{p}_i
pi进行修正以满足这些约束
第（12）-（15）行*：顶点的位置被调整到最优估计位置，并对速度进行更新。
第（16）行：对于碰撞的点会根据摩擦力和恢复系数修正。
求解器的输入是 M + M c o l l M+M_{coll} M+Mcoll个约束和 N N N个新顶点的估计位置 p 1 , … , p N \bm{p}_1, \dots, \bm{p}_N p1,…,pN。求解器的任务修正这些估计值，使其满足这些约束。这些约束往往是非线性的。比如最简单的位置约束： C ( p 1 , p 2 ) = ∣ p 1 − p 2 ∣ − d C(\bm{p}_1,\bm{p}_2) = |\bm{p}_1 - \bm{p}_2| - d C(p1,p2)=∣p1−p2∣−d，以及还包含了诸多不等式。
为了解决这些等式与不等式的集合，我们采用了高斯-赛德尔迭代(Gauss-Seidel-type Iteration, GS)。高斯-赛德尔迭代是数值代数中的一种迭代法，用于求解线性方程组（linear system of equations）。
对于由工程技术中产生的大型稀疏矩阵方程组（阶数很高，但零元素较多，例如求某些偏微分方程数值解所产生的线性方程组），利用迭代法求解此方程组比较合适，在计算机内存和运算两方面，迭代法通常都可利用矩阵中有大量零元素的特点。[参考介绍]：(https://blog.csdn.net/wwanrong/article/details/78927020)。
高斯-赛德尔迭代：
x
i
(
k
+
1
)
=
1
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
+
1
)
−
∑
j
=
i
+
1
n
a
i
j
x
j
(
k
)
)
,
i
=
1
,
2
,
…
,
n
x_{i}^{(k+1)}=\frac{1}{a_{i i}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i+1}^{n} a_{i j} x_{j}^{(k)}\right), \quad i=1,2, \ldots, n
xi(k+1)=aii1(bi−j=1∑i−1aijxj(k+1)−j=i+1∑naijxj(k)),i=1,2,…,n雅克比迭代：
x
i
(
k
+
1
)
=
1
a
i
i
(
b
i
−
∑
j
≠
i
a
i
j
x
j
(
k
)
)
,
i
=
1
,
2
,
…
,
n
x_{i}^{(k+1)}=\frac{1}{a_{i i}}\left(b_{i}-\sum_{j \neq i} a_{i j} x_{j}^{(k)}\right), \quad i=1,2, \ldots, n
xi(k+1)=aii1⎝⎛bi−j=i∑aijxj(k)⎠⎞,i=1,2,…,n这里我们使用了GS的思想独立地一个接一个地求解每一个约束。
投影的目的是为了使移动后的点仍然满足约束。最重要的是保护动量和角动量守恒。 ∑ i m i Δ p i = 0 \sum_{i} m_{i} \Delta \mathbf{p}_{i}=\mathbf{0} ∑imiΔpi=0 和 ∑ i r i × m i Δ p i = 0 \sum_{i} \mathbf{r}_{i} \times m_{i} \Delta \mathbf{p}_{i}=\mathbf{0} ∑iri×miΔpi=0。当投影后约束不成立，则需要添加虚拟力(ghost force), 然而,只有内部约束需要这样保存动量。 碰撞或附着约束对对象产生影响是允许的。（被碰了以后显然不需要再加一个力保持它的稳定）。
假设有 n n n个约束 C C C以及刚度系数 k k k, p = [ p 1 T , … , p n T ] T \bm{p} =[\bm{p}_1^T, \dots, \bm{p}_n^T]^T p=[p1T,…,pnT]T,内部约束 C C C与刚体的状态无关。 C ( p + Δ p ) ≈ C ( p ) + ∇ p C ( p ) ⋅ Δ p = 0 C(\mathbf{p}+\Delta \mathbf{p}) \approx C(\mathbf{p})+\nabla_{\mathbf{p}} C(\mathbf{p}) \cdot \Delta \mathbf{p}=0 C(p+Δp)≈C(p)+∇pC(p)⋅Δp=0意味着方向一致。 Δ p = λ ∇ p C ( p ) \Delta \mathbf{p}=\lambda \nabla_{\mathbf{p}} C(\mathbf{p}) Δp=λ∇pC(p)
求解
λ
\lambda
λ
Δ
p
=
−
C
(
p
)
∣
∇
p
C
(
p
)
∣
2
∇
p
C
(
p
)
\Delta \mathbf{p}=-\frac{C(\mathbf{p})}{\left|\nabla_{\mathbf{p}} C(\mathbf{p})\right|^{2}} \nabla_{\mathbf{p}} C(\mathbf{p})
Δp=−∣∇pC(p)∣2C(p)∇pC(p)这是由单个约束给出的非线性方程的迭代解的常规 Newton-Raphson 步骤。
Δ
p
i
=
−
s
∇
p
i
C
(
p
1
,
…
,
p
n
)
\Delta \mathbf{p}_{i}=-s \nabla_{\mathbf{p}_{i}} C\left(\mathbf{p}_{1}, \ldots, \mathbf{p}_{n}\right)
Δpi=−s∇piC(p1,…,pn)其中
s
=
C
(
p
1
,
…
,
p
n
)
∑
j
∣
∇
p
j
C
(
p
1
,
…
,
p
n
)
∣
2
s=\frac{C\left(\mathbf{p}_{1}, \ldots, \mathbf{p}_{n}\right)}{\sum_{j}\left|\nabla_{\mathbf{p}_{j}} C\left(\mathbf{p}_{1}, \ldots, \mathbf{p}_{n}\right)\right|^{2}}
s=∑j∣∣∇pjC(p1,…,pn)∣∣2C(p1,…,pn)
如果质点有质量则乘以系数
w
i
=
1
/
m
i
w_{i}=1 / m_{i}
wi=1/mi
Δ
p
=
λ
w
i
∇
p
C
(
p
)
\Delta \mathbf{p}=\lambda w_{i} \nabla_{\mathbf{p}} C(\mathbf{p})
Δp=λwi∇pC(p)其中
s
=
C
(
p
1
,
…
,
p
n
)
∑
j
w
j
∣
∇
p
j
C
(
p
1
,
…
,
p
n
)
∣
2
s=\frac{C\left(\mathbf{p}_{1}, \ldots, \mathbf{p}_{n}\right)}{\sum_{j}w_{j}\left|\nabla_{\mathbf{p}_{j}} C\left(\mathbf{p}_{1}, \ldots, \mathbf{p}_{n}\right)\right|^{2}}
s=∑jwj∣∣∇pjC(p1,…,pn)∣∣2C(p1,…,pn)
Δ
p
i
=
−
s
w
i
∇
p
i
C
(
p
1
,
…
,
p
n
)
\Delta \mathbf{p}_{i}=-s w_i \nabla_{\mathbf{p}_{i}} C\left(\mathbf{p}_{1}, \ldots, \mathbf{p}_{n}\right)
Δpi=−swi∇piC(p1,…,pn)
在算法的第（8）行， M M M个约束在物体建模时就已经固定， M c o l l M_{coll} Mcoll在每一个时刻更新，当布料的质点进入物体内部，则添加约束来防止该点进入物体内部。
但是本文讨论的碰撞仅适用于与静态对象的碰撞，因为没有冲量传递。（当动态碰撞的时候需要重新设计）
在算法的第（6）行，可以使用已有的一些阻尼算法进行“减速”，如下算法，
r
i
=
x
i
−
x
c
m
\mathbf{r}_{i}=\mathbf{x}_{i}-\mathbf{x}_{c m}
ri=xi−xcm，
k
d
a
m
p
i
n
g
∈
[
0
…
1
]
k_{damping} \in [0 \dots 1]
kdamping∈[0…1]是阻尼系数，当阻尼系数是1的时候，所有的点就像刚体一样运动。
使用基于位置的方法，将顶点附加到静态或运动对象非常简单。 顶点的位置简单地设置为静态目标位置或在每个时间步更新以与运动对象的位置重合。 为了确保包含该顶点的其他约束不会移动它，它的反质量 wi 设置为零。
我们的布料模拟器接受任意三角形网格作为输入。 我们对输入网格施加的唯一限制是它代表一个流形，即每条边最多由两个三角形共享。网格的每个节点都成为仿真中的一个顶点。 使用者需要提供布料的密度
ρ
[
k
g
/
m
2
]
\rho [kg/m^2]
ρ[kg/m2]。每个顶点的质量是相邻三角形质量之和的1/3。
m
i
=
∑
m
Δ
3
m_i = \frac{\sum m_{\Delta}} {3}
mi=3∑mΔ
拉伸约束：对于每条边，我们为拉伸约束设置了拉伸函数（等式）：
C
stretch
(
p
1
,
p
2
)
=
∣
p
1
−
p
2
∣
−
l
0
,
C_{\text {stretch }}\left(\mathbf{p}_{1}, \mathbf{p}_{2}\right)=\left|\mathbf{p}_{1}-\mathbf{p}_{2}\right|-l_{0},
Cstretch (p1,p2)=∣p1−p2∣−l0,
l
0
l_0
l0是边的初始长度，刚度系数
k
s
t
r
e
t
c
h
k_{stretch}
kstretch是全局变量由读者给定。
弯曲约束：
ϕ
0
\phi_0
ϕ0是初始二面角 ，刚度系数
k
b
e
n
d
k_{bend}
kbend是全局变量由读者给定。
C
bend
(
p
1
,
p
2
,
p
3
,
p
4
)
=
acos
(
(
p
2
−
p
1
)
×
(
p
3
−
p
1
)
∣
(
p
2
−
p
1
)
×
(
p
3
−
p
1
)
∣
⋅
(
p
2
−
p
1
)
×
(
p
4
−
p
1
)
∣
(
p
2
−
p
1
)
×
(
p
4
−
p
1
)
∣
)
−
φ
0
,
\begin{gathered} C_{\text {bend }}\left(\mathbf{p}_{1}, \mathbf{p}_{2}, \mathbf{p}_{3}, \mathbf{p}_{4}\right)= \\ \operatorname{acos}\left(\frac{\left(\mathbf{p}_{2}-\mathbf{p}_{1}\right) \times\left(\mathbf{p}_{3}-\mathbf{p}_{1}\right)}{\left|\left(\mathbf{p}_{2}-\mathbf{p}_{1}\right) \times\left(\mathbf{p}_{3}-\mathbf{p}_{1}\right)\right|} \cdot \frac{\left(\mathbf{p}_{2}-\mathbf{p}_{1}\right) \times\left(\mathbf{p}_{4}-\mathbf{p}_{1}\right)}{\left|\left(\mathbf{p}_{2}-\mathbf{p}_{1}\right) \times\left(\mathbf{p}_{4}-\mathbf{p}_{1}\right)\right|}\right)-\varphi_{0}, \end{gathered}
Cbend (p1,p2,p3,p4)=acos(∣(p2−p1)×(p3−p1)∣(p2−p1)×(p3−p1)⋅∣(p2−p1)×(p4−p1)∣(p2−p1)×(p4−p1))−φ0,
见3.4节
假设每个三角形大小一致，我们使用空间散列来查找顶点三角形碰撞，如果一个点
q
\bm{q}
q经过三角形。
C
(
q
,
p
1
,
p
2
,
p
3
)
=
(
q
−
p
1
)
⋅
(
p
2
−
p
1
)
×
(
p
3
−
p
1
)
∣
(
p
2
−
p
1
)
×
(
p
3
−
p
1
)
∣
−
h
C\left(\mathbf{q}, \mathbf{p}_{1}, \mathbf{p}_{2}, \mathbf{p}_{3}\right)=\left(\mathbf{q}-\mathbf{p}_{1}\right) \cdot \frac{\left(\mathbf{p}_{2}-\mathbf{p}_{1}\right) \times\left(\mathbf{p}_{3}-\mathbf{p}_{1}\right)}{\left|\left(\mathbf{p}_{2}-\mathbf{p}_{1}\right) \times\left(\mathbf{p}_{3}-\mathbf{p}_{1}\right)\right|}-h
C(q,p1,p2,p3)=(q−p1)⋅∣(p2−p1)×(p3−p1)∣(p2−p1)×(p3−p1)−h
h
h
h是布的厚度，如果顶点相对于三角形法线从下方进入，则约束函数必须为
C
(
q
,
p
1
,
p
2
,
p
3
)
=
(
q
−
p
1
)
⋅
(
p
3
−
p
1
)
×
(
p
2
−
p
1
)
∣
(
p
3
−
p
1
)
×
(
p
2
−
p
1
)
∣
−
h
C\left(\mathbf{q}, \mathbf{p}_{1}, \mathbf{p}_{2}, \mathbf{p}_{3}\right)=\left(\mathbf{q}-\mathbf{p}_{1}\right) \cdot \frac{\left(\mathbf{p}_{3}-\mathbf{p}_{1}\right) \times\left(\mathbf{p}_{2}-\mathbf{p}_{1}\right)}{\left|\left(\mathbf{p}_{3}-\mathbf{p}_{1}\right) \times\left(\mathbf{p}_{2}-\mathbf{p}_{1}\right)\right|}-h
C(q,p1,p2,p3)=(q−p1)⋅∣(p3−p1)×(p2−p1)∣(p3−p1)×(p2−p1)−h
C ( p 1 , … , p N ) = ( ∑ i = 1 n triangles ( p t 1 i × p t 2 i ) ⋅ p t 3 i ) − k pressure V 0 C\left(\mathbf{p}_{1}, \ldots, \mathbf{p}_{N}\right)=\left(\sum_{i=1}^{n_{\text {triangles }}}\left(\mathbf{p}_{t_{1}^{i}} \times \mathbf{p}_{t_{2}^{i}}\right) \cdot \mathbf{p}_{t_{3}^{i}}\right)-k_{\text {pressure }} V_{0} C(p1,…,pN)=(i=1∑ntriangles (pt1i×pt2i)⋅pt3i)−kpressure V0
∇ p i C = ∑ j : t 1 j = i ( p t 2 j × p t 3 j ) + ∑ j : t 2 j = i ( p t 3 j × p t 1 j ) + ∑ j : t 3 j = i ( p t 1 j × p t 2 j ) \nabla_{\mathbf{p}_{i}} C=\sum_{j: t_{1}^{j}=i}\left(\mathbf{p}_{t_{2}^{j}} \times \mathbf{p}_{t_{3}^{j}}\right)+\sum_{j: t_{2}^{j}=i}\left(\mathbf{p}_{t_{3}^{j}} \times \mathbf{p}_{t_{1}^{j}}\right)+\sum_{j: t_{3}^{j}=i}\left(\mathbf{p}_{t_{1}^{j}} \times \mathbf{p}_{t_{2}^{j}}\right) ∇piC=j:t1j=i∑(pt2j×pt3j)+j:t2j=i∑(pt3j×pt1j)+j:t3j=i∑(pt1j×pt2j)
