符号约定
-
离散时间节点 t 1 , t 2 , . . . , t n , . . . t_1,t_2,...,t_n,... t1,t2,...,tn,... 和固定时间间隔 h = t n − t n − 1 h=t_n-t_{n-1} h=tn−tn−1
-
系统节点数目 m m m ,系统弹簧数目 s s s
-
系统状态 q n ∈ R 3 m \boldsymbol{q}_n\in\mathbb{R}^{3m} qn∈R3m,描述在 t n t_n tn 时刻每一个节点的空间位置
-
质量矩阵 M ∈ R 3 m × 3 m \boldsymbol{M}\in\mathbb{R}^{3m\times 3m} M∈R3m×3m,是一个对角矩阵
-
速度 v n ∈ R 3 m \boldsymbol{v}_n\in\mathbb{R}^{3m} vn∈R3m,描述在 t n t_n tn 时刻每一个节点的速度向量
-
外力 f ( q n ) : R 3 m → R 3 m \boldsymbol{f}(\boldsymbol{q}_n):\mathbb{R}^{3m}\rightarrow\mathbb{R}^{3m} f(qn):R3m→R3m
-
系统能量 E : R 3 m → R E:\mathbb{R}^{3m}\rightarrow\mathbb{R} E:R3m→R,有 f = − ∇ E \boldsymbol{f}=-\nabla E f=−∇E
公式推导
根据牛顿第二定律得到关于系统状态的隐式积分
q
n
+
1
=
q
n
+
h
v
n
+
1
v
n
+
1
=
v
n
+
h
M
−
1
f
(
q
n
+
1
)
(1)
\begin{aligned} \boldsymbol{q}_{n+1}&=\boldsymbol{q}_n+h\boldsymbol{v}_{n+1} \\ \boldsymbol{v}_{n+1}&=\boldsymbol{v}_n+h\boldsymbol{M}^{-1}\boldsymbol{f}(\boldsymbol{q}_{n+1}) \end{aligned} \tag{1}
qn+1vn+1=qn+hvn+1=vn+hM−1f(qn+1)(1)
我们消去
v
\boldsymbol{v}
v 可以得到
q
n
+
1
−
2
q
n
+
q
n
−
1
=
h
2
M
−
1
f
(
q
n
+
1
)
(3)
\boldsymbol{q}_{n+1}-2\boldsymbol{q}_n+\boldsymbol{q}_{n-1}=h^2\boldsymbol{M}^{-1}\boldsymbol{f}(\boldsymbol{q}_{n+1}) \tag{3}
qn+1−2qn+qn−1=h2M−1f(qn+1)(3)
物理仿真的目的是根据已知条件求解
q
n
+
1
\boldsymbol{q}_{n+1}
qn+1,因此我们可以简记
x
:
=
q
n
+
1
\boldsymbol{x}:=\boldsymbol{q}_{n+1}
x:=qn+1,
y
:
=
2
q
n
−
q
n
−
1
\boldsymbol{y}:=2\boldsymbol{q}_n-\boldsymbol{q}_{n-1}
y:=2qn−qn−1,那么方程就变成了
M
(
x
−
y
)
=
h
2
f
(
x
)
(4)
\boldsymbol{M}(\boldsymbol{x}-\boldsymbol{y})=h^2\boldsymbol{f}(\boldsymbol{x}) \tag{4}
M(x−y)=h2f(x)(4)
式中
y
\boldsymbol{y}
y 可以看成是惯性,如果没有外力,那么下一刻系统状态将会变成
q
n
+
(
q
n
−
q
n
−
1
)
\boldsymbol{q}_n+(\boldsymbol{q}_n-\boldsymbol{q}_{n-1})
qn+(qn−qn−1)。要注意的是我们平时容易把
y
\boldsymbol{y}
y 当成函数来看,但这里方程(4)中只有
x
\boldsymbol{x}
x 和
f
(
x
)
\boldsymbol{f}(\boldsymbol{x})
f(x) 是未知的,其它都是已知的常量。
构造辅助函数
g
(
x
)
=
1
2
(
x
−
y
)
T
M
(
x
−
y
)
+
h
2
E
(
x
)
(5)
g(\boldsymbol{x})=\frac{1}{2}(\boldsymbol{x}-\boldsymbol{y})^{T}\boldsymbol{M}(\boldsymbol{x}-\boldsymbol{y})+h^2E(\boldsymbol{x}) \tag{5}
g(x)=21(x−y)TM(x−y)+h2E(x)(5)
求解(4)等价于求解(5)的极小值点。
典型的求解非线性函数极值点的方法是牛顿法:
x n + 1 = x n − H e s s g ( x n ) − 1 ∇ g ( x n ) \boldsymbol{x}_{n+1}=\boldsymbol{x}_n-Hess_g(\boldsymbol{x}_n)^{-1}\nabla g(\boldsymbol{x}_n) xn+1=xn−Hessg(xn)−1∇g(xn)
但是牛顿法每一次迭代中,海森矩阵都会变化,因此每次迭代都要求一次海森矩阵的逆,非常费时间。
接下来就是该算法的精华所在。
根据胡克定律,系统能量(弹性势能)表示为:
E
(
x
)
=
1
2
∑
i
=
1
s
k
i
(
∣
∣
p
i
1
−
p
i
2
∣
∣
−
r
)
2
(6)
E(\boldsymbol{x})=\frac{1}{2}\sum\limits_{i=1}^{s}k_i(||\boldsymbol{p}_{i_1}-\boldsymbol{p}_{i_2}||-r)^2 \tag{6}
E(x)=21i=1∑ski(∣∣pi1−pi2∣∣−r)2(6)
其中
k
i
k_i
ki 是弹性系数,
p
i
1
\boldsymbol{p}_{i_1}
pi1 和
p
i
2
\boldsymbol{p}_{i_2}
pi2 是弹簧的端点坐标,
i
1
,
i
2
∈
{
1
,
2
,
.
.
.
,
m
}
i_1,i_2\in\{1,2,...,m\}
i1,i2∈{1,2,...,m},
r
r
r 是弹簧原长。但是这个式子是不好处理的,因此论文中引入辅助变量
d
i
\boldsymbol{d}_i
di,将弹性势能等价表示为:
E
(
x
)
=
min
d
∈
U
1
2
∑
i
=
1
s
k
i
∣
∣
p
i
1
−
p
i
2
−
d
i
∣
∣
2
(7)
E(\boldsymbol{x})=\min\limits_{\boldsymbol{d}\in U}\frac{1}{2}\sum\limits_{i=1}^sk_i||\boldsymbol{p}_{i_1}-\boldsymbol{p}_{i_2}-\boldsymbol{d}_i||^2 \tag{7}
E(x)=d∈Umin21i=1∑ski∣∣pi1−pi2−di∣∣2(7)
其中
U
=
{
(
d
1
,
.
.
.
,
d
s
)
∈
R
3
s
:
∣
∣
d
i
∣
∣
=
r
i
}
U=\{(\boldsymbol{d}_1,...,\boldsymbol{d}_s)\in\mathbb{R}^{3s}:||\boldsymbol{d}_i|| = r_i \}
U={(d1,...,ds)∈R3s:∣∣di∣∣=ri} (论文里这里有一点纰漏)。
(6)和(7)等价的原因是
(
∣
∣
p
1
−
p
2
∣
∣
−
r
)
2
=
min
∣
∣
d
∣
∣
=
r
∣
∣
p
1
−
p
2
−
d
∣
∣
2
(||\boldsymbol{p}_1-\boldsymbol{p}_2||-r)^2=\min\limits_{||\boldsymbol{d}||=r}||\boldsymbol{p}_1-\boldsymbol{p}_2-\boldsymbol{d}||^2
(∣∣p1−p2∣∣−r)2=∣∣d∣∣=rmin∣∣p1−p2−d∣∣2 ,这个很好理解,一个点距离球面最近的点,就是这个点与球心连线(的延长线)和球面的交点。变换成 (7) 式之后,论文中直接华丽地把公式改写成了矩阵形式,我这里做一点点简单推导:
∣
∣
p
1
−
p
2
−
d
∣
∣
2
=
∣
∣
p
1
∣
∣
2
+
∣
∣
p
2
∣
∣
2
+
∣
∣
d
∣
∣
2
−
2
p
1
⋅
p
2
−
2
p
1
⋅
d
+
2
p
2
⋅
d
=
∣
∣
p
1
−
p
2
∣
∣
2
−
2
(
p
1
−
p
2
)
⋅
d
+
r
2
\begin{aligned} ||\boldsymbol{p}_1-\boldsymbol{p}_2-\boldsymbol{d}||^2 &= ||\boldsymbol{p}_1||^2 + ||\boldsymbol{p}_2||^2 + ||d||^2 - 2\boldsymbol{p}_1\cdot\boldsymbol{p}_2 - 2\boldsymbol{p}_1\cdot\boldsymbol{d} + 2\boldsymbol{p}_2\cdot\boldsymbol{d} \\ &= ||\boldsymbol{p}_1-\boldsymbol{p}_2||^2 - 2(\boldsymbol{p}_1-\boldsymbol{p}_2)\cdot\boldsymbol{d} + r^2 \end{aligned}
∣∣p1−p2−d∣∣2=∣∣p1∣∣2+∣∣p2∣∣2+∣∣d∣∣2−2p1⋅p2−2p1⋅d+2p2⋅d=∣∣p1−p2∣∣2−2(p1−p2)⋅d+r2
由于是优化问题,常数项
r
2
r^2
r2 可以丢弃,因此式中剩下一个很像二次函数的东西,于是就可以得到论文中的华丽变换(论文中这个公式其实不是严格等号,差了一个
r
2
r^2
r2 常数项,加上 min 会好一些):
min
d
∈
U
1
2
∑
i
=
1
s
k
i
∣
∣
p
i
1
−
p
i
2
−
d
i
∣
∣
2
=
min
d
∈
U
1
2
x
T
L
x
−
x
T
J
d
(8)
\min\limits_{\boldsymbol{d}\in U}\frac{1}{2}\sum\limits_{i=1}^sk_i||\boldsymbol{p}_{i_1}-\boldsymbol{p}_{i_2}-\boldsymbol{d}_i||^2 = \min\limits_{\boldsymbol{d}\in U} \frac{1}{2}\boldsymbol{x}^T\boldsymbol{L}\boldsymbol{x}-\boldsymbol{x}^T\boldsymbol{J}\boldsymbol{d} \tag{8}
d∈Umin21i=1∑ski∣∣pi1−pi2−di∣∣2=d∈Umin21xTLx−xTJd(8)
其中
x
=
(
p
1
,
.
.
.
,
p
m
)
\boldsymbol{x}=(\boldsymbol{p}_1,...,\boldsymbol{p}_m)
x=(p1,...,pm),矩阵
L
∈
R
3
m
×
3
m
,
J
∈
R
3
m
×
3
s
\boldsymbol{L}\in\mathbb{R}^{3m\times 3m},\boldsymbol{J}\in\mathbb{R}^{3m\times3s}
L∈R3m×3m,J∈R3m×3s 定义如下:
L
=
(
∑
i
=
1
s
k
i
A
i
A
i
T
)
⨂
I
3
J
=
(
∑
i
=
1
s
k
i
A
i
S
i
T
)
⨂
I
3
(9)
\boldsymbol{L}=(\sum\limits_{i=1}^s k_i \boldsymbol{A}_i\boldsymbol{A}_i^T) \bigotimes \boldsymbol{I}_3 \\ \boldsymbol{J}=(\sum\limits_{i=1}^s k_i \boldsymbol{A}_i\boldsymbol{S}_i^T) \bigotimes \boldsymbol{I}_3 \tag{9}
L=(i=1∑skiAiAiT)⨂I3J=(i=1∑skiAiSiT)⨂I3(9)
其中
A
i
∈
R
m
\boldsymbol{A}_i\in\mathbb{R}^m
Ai∈Rm 是关于弹簧端点索引
i
1
,
i
2
i_1,i_2
i1,i2 的向量,有
A
i
,
i
1
=
1
,
A
i
,
i
2
=
−
1
A_{i,i_1}=1, \ A_{i,i2}=-1
Ai,i1=1, Ai,i2=−1,其它向量元素都是 0 ;而
S
i
∈
R
s
\boldsymbol{S}_i\in\mathbb{R}^s
Si∈Rs 是关于弹簧索引
i
i
i 的向量,有
S
i
,
i
=
1
S_{i,i}=1
Si,i=1,其它向量元素都是 0 ;
⨂
\bigotimes
⨂ 是克罗内克积。
如果算上外力
f
e
x
t
∈
R
3
m
\boldsymbol{f}_{ext}\in\mathbb{R}^{3m}
fext∈R3m,系统能量可以表示为:
E
(
x
)
=
min
d
∈
U
1
2
x
T
L
x
−
x
T
J
d
−
x
T
f
e
x
t
(10)
E(\boldsymbol{x})=\min\limits_{\boldsymbol{d}\in U} \frac{1}{2}\boldsymbol{x}^T\boldsymbol{L}\boldsymbol{x}-\boldsymbol{x}^T\boldsymbol{J}\boldsymbol{d} - \boldsymbol{x}^T\boldsymbol{f}_{ext} \tag{10}
E(x)=d∈Umin21xTLx−xTJd−xTfext(10)
将(10)代入(5)得到最终的优化问题
min
x
∈
R
3
m
,
d
∈
U
1
2
x
T
(
M
+
h
2
L
)
x
−
h
2
x
T
J
d
−
h
2
x
T
f
e
x
t
−
1
2
x
T
M
y
−
1
2
y
T
M
x
(11)
\min\limits_{\boldsymbol{x}\in\mathbb{R}^{3m},\boldsymbol{d}\in U} \frac{1}{2}\boldsymbol{x}^T(\boldsymbol{M}+h^2\boldsymbol{L})\boldsymbol{x} - h^2\boldsymbol{x}^T\boldsymbol{J}\boldsymbol{d} - h^2\boldsymbol{x}^T\boldsymbol{f}_{ext} - \frac{1}{2}\boldsymbol{x}^T\boldsymbol{M}\boldsymbol{y} - \frac{1}{2}\boldsymbol{y}^T\boldsymbol{M}\boldsymbol{x} \tag{11}
x∈R3m,d∈Umin21xT(M+h2L)x−h2xTJd−h2xTfext−21xTMy−21yTMx(11)
式中舍弃了一个常数项
1
2
y
T
M
y
\frac{1}{2}\boldsymbol{y}^T\boldsymbol{M}\boldsymbol{y}
21yTMy 。
数值计算
论文中描述的求解(11)的方法是块坐标下降法,需要计算若干个迭代,每次迭代有两个步骤:local step 和 global step ,分别说明如下:
-
local step 是固定 x \boldsymbol{x} x 计算 d \boldsymbol{d} d,其实就是把 d \boldsymbol{d} d 投影到弹簧原长的位置,具体操作为遍历每个弹簧,找到其两个端点 p i 1 \boldsymbol{p}_{i_1} pi1 和 p i 2 \boldsymbol{p}_{i_2} pi2,记 p i 12 = p i 1 − p i 2 \boldsymbol{p}_{i_{12}}=\boldsymbol{p}_{i_1}-\boldsymbol{p}_{i_2} pi12=pi1−pi2,那么就更新 d i = r i p i 12 ∣ ∣ p i 12 ∣ ∣ \boldsymbol{d}_i = r_i \frac{\boldsymbol{p}_{i_{12}}}{||\boldsymbol{p}_{i_{12}}||} di=ri∣∣pi12∣∣pi12 。
-
global step 是固定 d \boldsymbol{d} d 计算 x \boldsymbol{x} x,也就是求一个凸二次优化问题。 M + h 2 L \boldsymbol{M}+h^2\boldsymbol{L} M+h2L 是一个对称正定矩阵,求解(11)的最小值可以令其梯度等于 0 ,也就是:
( M + h 2 L ) x = h 2 J d + M y + h 2 f e x t (12) (\boldsymbol{M}+h^2\boldsymbol{L})\boldsymbol{x} = h^2\boldsymbol{J}\boldsymbol{d}+\boldsymbol{My}+h^2\boldsymbol{f}_{ext} \tag{12} (M+h2L)x=h2Jd+My+h2fext(12)
因此该步骤就等价于求解一个线性方程组。因为 M + h 2 L \boldsymbol{M}+h^2\boldsymbol{L} M+h2L 一直是不变的,所以可以在初始化的时候先做稀疏矩阵 Cholesky 分解( A = L L T A=LL^T A=LLT),可节省大量时间。
如果考虑速度阻尼 α \alpha α,只需令 y = q n + α ( q n − q n − 1 ) \boldsymbol{y}=\boldsymbol{q}_n+\alpha(\boldsymbol{q}_n-\boldsymbol{q}_{n-1}) y=qn+α(qn−qn−1) 即可。
demo
我基于 fast mass spring 用OpenGL(GLFW)做了布料仿真。布料的配置过程是这样的:
- 正方形布料,边节点数目为 n n n ,总节点数目为 n 2 n^2 n2 ( n > 1 n>1 n>1);
- 布料宽度为 w w w,相邻节点间隔为 w / ( n − 1 ) w/(n-1) w/(n−1);
- 布料弹簧数目为 6 n 2 − 10 n + 2 6n^2-10n+2 6n2−10n+2,其中 structural 类型的弹簧 2 n ( n − 1 ) 2n(n-1) 2n(n−1) 个,shearing类型的弹簧 2 ( n − 1 ) 2 2(n-1)^2 2(n−1)2 个,bending类型的弹簧 2 n ( n − 2 ) 2n(n-2) 2n(n−2) 个;
设置悬挂节点质量为无穷大,外力为0,得到如下效果:
![]() |
![]() |
![]() |
![]() |
该算法在 Intel i7-11700K CPU 单核上,3600个节点,2万多个弹簧,每一帧迭代 10次,运行帧率可以超过60帧,非常丝滑~