Games 103 作业三
作业三的内容主要就是实现一下FVM。我们按照文档中的步骤,第一步就是去独立地更新mesh的速度和位置,在初始化每个顶点的受力时,需要考虑到重力的影响。
for(int i=0 ;i<number; i++)
{
//TODO: Add gravity to Force.
Force[i] = new Vector3(0, -9.8f * mass, 0);
}
for(int i=0; i<number; i++)
{
//TODO: Update X and V here.
V[i] += Force[i] * dt / mass;
V[i] *= damp;
X[i] += V[i]*dt;
}
接下来要考虑mesh的碰撞处理。在作业给的场景中,就只有一个简单的floor。
floor的世界坐标是(0, -3, 0),那么这里就直接简化碰撞处理了:
//TODO: (Particle) collision with floor.
if (X[i].y < -3f)
{
V[i].y += (-3f - X[i].y) / dt;
X[i].y = -3f;
}
下一步,我们要在Start方法中预先计算好常量矩阵
D
m
\textbf{D}_m
Dm以及它的逆。根据PPT的第25页,有
D
m
=
[
X
10
X
20
X
30
]
\textbf{D}_m = [\textbf{X}_{10} \textbf{X}_{20} \textbf{X}_{30}]
Dm=[X10X20X30]
每个四面体都有这样的一个矩阵:
//TODO: Need to allocate and assign inv_Dm
inv_Dm = new Matrix4x4[tet_number];
for(int i = 0; i < tet_number; i++)
{
Matrix4x4 Dm = Build_Edge_Matrix(i);
inv_Dm[i] = Dm.inverse;
}
Matrix4x4 Build_Edge_Matrix(int tet)
{
Matrix4x4 ret=Matrix4x4.zero;
//TODO: Need to build edge matrix here.
Vector4 e1 = X[Tet[tet * 4 + 1]] - X[Tet[tet * 4 + 0]];
Vector4 e2 = X[Tet[tet * 4 + 2]] - X[Tet[tet * 4 + 0]];
Vector4 e3 = X[Tet[tet * 4 + 3]] - X[Tet[tet * 4 + 0]];
ret.SetColumn(0, e1);
ret.SetColumn(1, e2);
ret.SetColumn(2, e3);
ret.SetColumn(3, new Vector4(0, 0, 0, 1));
return ret;
}
接下来,计算Deformation gradient F和Green strain G,计算的方法也在PPT第25页:
F
=
[
x
10
x
20
x
30
]
D
m
−
1
\textbf{F} = [\textbf{x}_{10} \textbf{x}_{20} \textbf{x}_{30}] \textbf{D}_m^{-1}
F=[x10x20x30]Dm−1
G = 1 2 ( F T F − I ) \textbf{G} = \dfrac{1}{2}(\textbf{F}^T\textbf{F} - \textbf{I}) G=21(FTF−I)
//TODO: Deformation Gradient
Matrix4x4 m = Build_Edge_Matrix(tet);
Matrix4x4 F = m * inv_Dm[tet];
//TODO: Green Strain
Matrix4x4 G = Matrix_Mul(0.5f, Matrix_Sub(F.transpose * F, Matrix4x4.identity));
作业里使用的是the second Piola-Kirchho stress,也给出了相应的公式:
S
=
2
s
1
G
+
s
0
t
r
(
G
)
I
\textbf{S} = 2s_1\textbf{G} + s_0tr(\textbf{G})\textbf{I}
S=2s1G+s0tr(G)I
//TODO: Second PK Stress
Matrix4x4 S = Matrix_Add(Matrix_Mul(2.0f * stiffness_1, G),
Matrix_Mul(stiffness_0 * Matrix_Trace(G), Matrix4x4.identity));
那么根据ppt,就可以计算出四面体上每个顶点的力了:
P
=
FS
\textbf{P} = \textbf{F}\textbf{S}
P=FS
[ f 1 f 2 f 3 ] = − 1 6 d e t ( D m − 1 ) P D m − T [\textbf{f}_1 \textbf{f}_2 \textbf{f}_3] = - \dfrac{1}{6det(\textbf{D}_m^{-1})}\textbf{P}\textbf{D}_m^{-\textbf{T}} [f1f2f3]=−6det(Dm−1)1PDm−T
f 0 = − f 1 − f 2 − f 3 \textbf{f}_0 = -\textbf{f}_1 - \textbf{f}_2 - \textbf{f}_3 f0=−f1−f2−f3
//TODO: Elastic Force
Matrix4x4 fm = Matrix_Mul(-1.0f / (6.0f * inv_Dm[tet].determinant),
F * S * inv_Dm[tet].transpose);
f1 = fm.GetColumn(0);
f2 = fm.GetColumn(1);
f3 = fm.GetColumn(2);
Force[Tet[tet * 4 + 0]] -= f1 + f2 + f3;
Force[Tet[tet * 4 + 1]] += f1;
Force[Tet[tet * 4 + 2]] += f2;
Force[Tet[tet * 4 + 3]] += f3;
其实这个时候房子已经可以跳动起来了,只是。。。
模拟很不稳定,房子直接飞了。作业中为了避免这一情况,使用了拉普拉斯平滑。我们需要统计每个顶点所邻接的所有其他顶点,计算出一个邻接的平均速度,然后进行加权:
void Laplacian_Smoothing(float blend)
{
for (int i = 0; i < number; i++)
{
V_sum[i] = Vector3.zero;
V_num[i] = 0;
}
for(int tet=0; tet<tet_number; tet++)
{
Vector3 sum = V[Tet[tet * 4 + 0]] + V[Tet[tet * 4 + 1]] + V[Tet[tet * 4 + 2]] + V[Tet[tet * 4 + 3]];
V_sum[Tet[tet * 4 + 0]] += sum;
V_sum[Tet[tet * 4 + 1]] += sum;
V_sum[Tet[tet * 4 + 2]] += sum;
V_sum[Tet[tet * 4 + 3]] += sum;
V_num[Tet[tet * 4 + 0]] += 4;
V_num[Tet[tet * 4 + 1]] += 4;
V_num[Tet[tet * 4 + 2]] += 4;
V_num[Tet[tet * 4 + 3]] += 4;
}
for(int i = 0; i < number; i++)
{
V[i] = (V[i] + blend * V_sum[i] / V_num[i]) / (1 + blend);
}
}
我们取blend为0.1,再看看现在的效果:
基础部分就完成了。下面再看一下附加题部分,其实就是根据PPT的第35页,换一种方法计算P。
首先是利用SVD得到三个矩阵:
Matrix4x4 U = Matrix4x4.zero;
Matrix4x4 S = Matrix4x4.zero;
Matrix4x4 V = Matrix4x4.zero;
svd.svd(F, ref U, ref S, ref V);
根据公式:
P
=
U
d
i
a
g
(
∂
W
∂
λ
0
,
∂
W
∂
λ
1
,
∂
W
∂
λ
2
)
V
T
P = \textbf{U} diag(\dfrac{\partial W}{\partial \lambda_0}, \dfrac{\partial W}{\partial \lambda_1}, \dfrac{\partial W}{\partial \lambda_2}) \textbf{V}^\textbf{T}
P=Udiag(∂λ0∂W,∂λ1∂W,∂λ2∂W)VT
那么问题就只剩计算这三个偏导数。而W是关于
I
C
I_C
IC和
I
I
C
II_C
IIC的函数(注意,这里PPT给的公式有误)
W
=
s
0
8
(
I
C
−
3
)
2
+
s
1
4
(
I
I
C
−
2
I
C
+
3
)
W = \dfrac{s_0}{8}(I_C - 3)^2 + \dfrac{s_1}{4}(II_C - 2I_C + 3)
W=8s0(IC−3)2+4s1(IIC−2IC+3)
而
I
C
=
λ
0
2
+
λ
1
2
+
λ
2
2
I_C = \lambda_0^2 + \lambda_1^2 + \lambda_2^2
IC=λ02+λ12+λ22
I I C = λ 0 4 + λ 1 4 + λ 2 4 II_C = \lambda_0^4 + \lambda_1^4 + \lambda_2^4 IIC=λ04+λ14+λ24
又有
∂
W
∂
λ
=
∂
W
∂
I
C
∂
I
C
∂
λ
+
∂
W
∂
I
I
C
∂
I
I
C
∂
λ
\dfrac{\partial W}{\partial \lambda} = \dfrac{\partial W}{\partial I_C} \dfrac{\partial I_C}{\partial \lambda} + \dfrac{\partial W}{\partial II_C} \dfrac{\partial II_C}{\partial \lambda}
∂λ∂W=∂IC∂W∂λ∂IC+∂IIC∂W∂λ∂IIC
所以只要分别计算上述式子中右侧的偏导数,就能得到最终的P了:
float lambda0 = S[0, 0];
float lambda1 = S[1, 1];
float lambda2 = S[2, 2];
float Ic = lambda0 * lambda0 + lambda1 * lambda1 + lambda2 * lambda2;
float dWdIc = 0.25f * stiffness_0 * (Ic - 3f) - 0.5f * stiffness_1;
float dWdIIc = 0.25f * stiffness_1;
float dIcdlambda0 = 2f * lambda0;
float dIcdlambda1 = 2f * lambda1;
float dIcdlambda2 = 2f * lambda2;
float dIIcdlambda0 = 4f * lambda0 * lambda0 * lambda0;
float dIIcdlambda1 = 4f * lambda1 * lambda1 * lambda1;
float dIIcdlambda2 = 4f * lambda2 * lambda2 * lambda2;
float dWd0 = dWdIc * dIcdlambda0 + dWdIIc * dIIcdlambda0;
float dWd1 = dWdIc * dIcdlambda1 + dWdIIc * dIIcdlambda1;
float dWd2 = dWdIc * dIcdlambda2 + dWdIIc * dIIcdlambda2;
Matrix4x4 diag = Matrix4x4.zero;
diag[0, 0] = dWd0;
diag[1, 1] = dWd1;
diag[2, 2] = dWd2;
diag[3, 3] = 1;
Matrix4x4 P = U * diag * V.transpose;
最终效果如下,感觉两种计算方式的效果差不多: