论文阅读笔记：Descent methods for elastic body simulation on the GPU (源代码及实现细节)
材料来源于 Descent methods for elastic body simulation on the GPU, ACMTransactions on Graphics (TOG), 2016.
0. 概述
在本论文中，提出了一种***。下面将详细介绍该方法的源代码及实现细节，并对照论文中的理论方法进行逐一说明。
1. 仿真算法总体框架及流程
本方法的仿真流程如 * 所示，采用迭代法求解每一帧的模型网格形变位移。在 \(t+1\) 时刻，其主要包括以下几个步骤：1）初始化，即给定迭代的初始预测值；2）计算搜索方向；3，根据搜索方向、步长等更新网格节点位置；4、调整步长，并判断是否收敛；5、某些加速方法。后续将一一详细介绍。
（算法伪代码，待更新）
2. 源代码及实现细节
该论文的作者公开了源代码。对照论文中的理论公式和源代码，进一步分析源代码的实现流程及细节。核心代码位于源文件 CUDA_HYPER_TET_MESH.h 中，主要计算流程为函数 float Update(TYPE t, int iterations, TYPE dir[]) 。
2.1 初始化迭代
理论方法
作者在论文中提到，our system choose the constant acceleration approach to initialize \(\boldsymbol{q}^{(0)}\) by default. 即采用的是恒定加速度的方式，设 \(\boldsymbol{q}^{(0)} = \boldsymbol{q}_{t} + h\boldsymbol{v}_{t} + \eta h(\boldsymbol{v}_{t} - \boldsymbol{v}_{t-1})\)，其中 \(\eta = 0.2\) 。
此外，在计算中需要使用中间变量 \(\boldsymbol{s}_{t} = \boldsymbol{q}_{t} + h\boldsymbol{v}_{t}\) ，在初始化阶段一并做了计算。
代码实现
迭代流程开始时，需要先给定预测值。在源代码中，是通过如下代码实现的：
__global__ void Update_Kernel(float* X, float* V, const float* prev_V, float* S, const float *fixed, const float *more_fixed, const float *offset_X, float* fixed_X, const float t, const int number, const float dir_x, const float dir_y, const float dir_z)
{
int i = blockDim.x * blockIdx.x + threadIdx.x;
if(i>=number) return;
if(more_fixed[i])
{
//fixed_X[i*3+0]=X[i*3+0]+dir_x;
//fixed_X[i*3+1]=X[i*3+1]+dir_y;
//fixed_X[i*3+2]=X[i*3+2]+dir_z;
fixed_X[i*3+0]=offset_X[i*3+0]+dir_x;
fixed_X[i*3+1]=offset_X[i*3+1]+dir_y;
fixed_X[i*3+2]=offset_X[i*3+2]+dir_z;
}
//V[i*3+1]-=2*t;
//Calculate S
S[i*3+0]=X[i*3+0]+V[i*3+0]*t;
S[i*3+1]=X[i*3+1]+V[i*3+1]*t;
S[i*3+2]=X[i*3+2]+V[i*3+2]*t;
//Initialize position
X[i * 3 + 0] += (V[i * 3 + 0] + (V[i * 3 + 0] - prev_V[i * 3 + 0])*0.2)*t;
X[i * 3 + 1] += (V[i * 3 + 1] + (V[i * 3 + 1] - prev_V[i * 3 + 1])*0.2)*t;
X[i * 3 + 2] += (V[i * 3 + 2] + (V[i * 3 + 2] - prev_V[i * 3 + 2])*0.2)*t;
}
可以看到，在迭代开始时，网格节点位置的初始值设定为 X += (V + (V - prev_V)*0.2)*t; ，即 X = X + t*V + 0.2*t*(V - prev_V) 。
此外，对于 \(\boldsymbol{s}_{t}\) 同样做了计算，即 S = X + V*t 。
2.2 计算搜索方向
初始化之后，便开始迭代过程。在每一个迭代步中，首先要计算搜索方向，然后根据搜索方向和步长更新网格节点位置，最后，再重新调整步长。下面，主要介绍搜索方向的计算实现细节。
理论方法
（待补充）
代码实现
计算搜索方向的代码位移 void Evaluate(TYPE stepping, TYPE t, bool update_C) 中，主要为其中的 Compute_FM_Kernel(...) 函数。具体为：
__global__ void Compute_FM_Kernel(const float* X, const int* Tet, const float* inv_Dm, const float* Vol, float* lambda, float* Force, float* C, float* ext_C, float *E,
const int model, float stiffness_0, float stiffness_1, float stiffness_2, float stiffness_3, float stiffness_p, const int tet_number, const float lower_bound, const float upper_bound, const bool update_C=true)
{
(省略部分代码)
atomicAdd(&Force[p0 + 0], force[ 0]);
atomicAdd(&Force[p0 + 1], force[ 1]);
atomicAdd(&Force[p0 + 2], force[ 2]);
atomicAdd(&Force[p1 + 0], force[ 3]);
atomicAdd(&Force[p1 + 1], force[ 4]);
atomicAdd(&Force[p1 + 2], force[ 5]);
atomicAdd(&Force[p2 + 0], force[ 6]);
atomicAdd(&Force[p2 + 1], force[ 7]);
atomicAdd(&Force[p2 + 2], force[ 8]);
atomicAdd(&Force[p3 + 0], force[ 9]);
atomicAdd(&Force[p3 + 1], force[10]);
atomicAdd(&Force[p3 + 2], force[11]);
（省略部分代码）
if(update_C==false) return;
//Now compute the stiffness matrix
float alpha[3][3], beta[3][3], gamma[3][3];
float vi[3], vj[3], r[3];
for(int i=0; i<3; i++)
for(int j=i; j<3; j++)
{
alpha[i][j]=2*dEdI+4*(Sigma[i]*Sigma[i]+Sigma[j]*Sigma[j])*dEdII;
beta[i][j]=4*Sigma[i]*Sigma[j]*dEdII-2*III*dEdIII/(Sigma[i]*Sigma[j]);
vi[0]=2*Sigma[i];
vi[1]=4*Sigma[i]*Sigma[i]*Sigma[i];
vi[2]=2*III/Sigma[i];
vj[0]=2*Sigma[j];
vj[1]=4*Sigma[j]*Sigma[j]*Sigma[j];
vj[2]=2*III/Sigma[j];
dev_Matrix_Vector_Product_3(&H[0][0], vj, r);
gamma[i][j]=DOT(vi, r)+4*III*dEdIII/(Sigma[i]*Sigma[j]);
}
float dGdF[12][9];
//G is related to force, according to [TSIF05], (g0, g3, g6), (g1, g4, g7), (g2, g5, g8), (g9, g10, g11)
float dF[9], temp0[9], temp1[9];
for(int i=0; i<9; i++)
{
memset(&dF, 0, sizeof(float)*9);
dF[i]=1;
dev_Matrix_Product_3(dF, V, temp0);
dev_Matrix_T_Product_3(U, temp0, temp1);
//diagonal
temp0[0]=(alpha[0][0]+beta[0][0]+gamma[0][0])*temp1[0]+gamma[0][1]*temp1[4]+gamma[0][2]*temp1[8];
temp0[4]=gamma[0][1]*temp1[0]+(alpha[1][1]+beta[1][1]+gamma[1][1])*temp1[4]+gamma[1][2]*temp1[8];
temp0[8]=gamma[0][2]*temp1[0]+gamma[1][2]*temp1[4]+(alpha[2][2]+beta[2][2]+gamma[2][2])*temp1[8];
//off-diagonal
temp0[1]=alpha[0][1]*temp1[1]+beta[0][1]*temp1[3];
temp0[3]=alpha[0][1]*temp1[3]+beta[0][1]*temp1[1];
temp0[2]=alpha[0][2]*temp1[2]+beta[0][2]*temp1[6];
temp0[6]=alpha[0][2]*temp1[6]+beta[0][2]*temp1[2];
temp0[5]=alpha[1][2]*temp1[5]+beta[1][2]*temp1[7];
temp0[7]=alpha[1][2]*temp1[7]+beta[1][2]*temp1[5];
dev_Matrix_Product_T_3(temp0, V, temp1);
dev_Matrix_Product_3(U, temp1, temp0);
dev_Matrix_Product_T_3(temp0, &inv_Dm[t*9], &dGdF[i][0]);
}
//Transpose dGdF
float temp;
for(int i=0; i<9; i++) for(int j=i+1; j<9; j++)
SWAP(dGdF[i][j], dGdF[j][i]);
for(int j=0; j< 9; j++)
{
dGdF[ 9][j]=-dGdF[0][j]-dGdF[1][j]-dGdF[2][j];
dGdF[10][j]=-dGdF[3][j]-dGdF[4][j]-dGdF[5][j];
dGdF[11][j]=-dGdF[6][j]-dGdF[7][j]-dGdF[8][j];
}
float new_idm[4][3];
new_idm[0][0]=-inv_Dm[t*9+0]-inv_Dm[t*9+3]-inv_Dm[t*9+6];
new_idm[0][1]=-inv_Dm[t*9+1]-inv_Dm[t*9+4]-inv_Dm[t*9+7];
new_idm[0][2]=-inv_Dm[t*9+2]-inv_Dm[t*9+5]-inv_Dm[t*9+8];
new_idm[1][0]= inv_Dm[t*9+0];
new_idm[1][1]= inv_Dm[t*9+1];
new_idm[1][2]= inv_Dm[t*9+2];
new_idm[2][0]= inv_Dm[t*9+3];
new_idm[2][1]= inv_Dm[t*9+4];
new_idm[2][2]= inv_Dm[t*9+5];
new_idm[3][0]= inv_Dm[t*9+6];
new_idm[3][1]= inv_Dm[t*9+7];
new_idm[3][2]= inv_Dm[t*9+8];
atomicAdd(&C[p0*3+0], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[ 9][0], 4, 3, 3, 0, 0));
atomicAdd(&C[p0*3+4], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[10][0], 4, 3, 3, 0, 1));
atomicAdd(&C[p0*3+8], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[11][0], 4, 3, 3, 0, 2));
atomicAdd(&C[p1*3+0], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[ 0][0], 4, 3, 3, 1, 0));
atomicAdd(&C[p1*3+4], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[ 3][0], 4, 3, 3, 1, 1));
atomicAdd(&C[p1*3+8], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[ 6][0], 4, 3, 3, 1, 2));
atomicAdd(&C[p2*3+0], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[ 1][0], 4, 3, 3, 2, 0));
atomicAdd(&C[p2*3+4], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[ 4][0], 4, 3, 3, 2, 1));
atomicAdd(&C[p2*3+8], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[ 7][0], 4, 3, 3, 2, 2));
atomicAdd(&C[p3*3+0], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[ 2][0], 4, 3, 3, 3, 0));
atomicAdd(&C[p3*3+4], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[ 5][0], 4, 3, 3, 3, 1));
atomicAdd(&C[p3*3+8], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[ 8][0], 4, 3, 3, 3, 2));
}
其他细节
（待补充）
2.3 更新网格节点位置
接下来，将根据搜索方向、迭代步长等，更新网格节点的位置。
理论方法
在第 \(k+1\) 个迭代步中，网格节点的位置更新为
\[\begin{aligned}
\boldsymbol{q}^{(k+1)} &= \boldsymbol{q}^{(k)} + \alpha^{(k)}\Delta\boldsymbol{q}^{(k)} \\
&= \boldsymbol{q}^{(k)} - \alpha^{(k)}\text{diag}^{-1}(\boldsymbol{H}^{(k)})\boldsymbol{g}^{(k)}
\end{aligned}\]
代码实现
这部分代码主要位于 Constraint_1_Kernel(...) 中，具体如下：
__global__ void Constraint_1_Kernel(const float* M, const float* X, const float* prev_X, const float* V, float* E, float* G, float* P, const float* S, float* next_X,
const float* fixed, const float* more_fixed, const float* fixed_X, float *F, const float* C, const float* ext_C, const float stepping, const int number, const float t, const float inv_t, const float gravity)
{
int i = blockDim.x * blockIdx.x + threadIdx.x;
if(i>=number) return;
float oc = M[i]*inv_t*inv_t;
float c = M[i]*inv_t*inv_t+fixed[i]+more_fixed[i];
// Get Force
float error[3];
error[0]=oc*(S[i*3+0]-X[i*3+0])+(c-oc)*(fixed_X[i*3+0]-X[i*3+0])+F[i*3+0];
error[1]=oc*(S[i*3+1]-X[i*3+1])+(c-oc)*(fixed_X[i*3+1]-X[i*3+1])+F[i*3+1]+gravity*M[i];
error[2]=oc*(S[i*3+2]-X[i*3+2])+(c-oc)*(fixed_X[i*3+2]-X[i*3+2])+F[i*3+2];
// Update Energy
float energy=0;
energy+=oc*(S[i*3+0]-X[i*3+0])*(S[i*3+0]-X[i*3+0]); // kinetic energy
energy+=oc*(S[i*3+1]-X[i*3+1])*(S[i*3+1]-X[i*3+1]); // kinetic energy
energy+=oc*(S[i*3+2]-X[i*3+2])*(S[i*3+2]-X[i*3+2]); // kinetic energy
energy+=(c-oc)*(fixed_X[i*3+0]-X[i*3+0])*(fixed_X[i*3+0]-X[i*3+0]); // fixed energy
energy+=(c-oc)*(fixed_X[i*3+1]-X[i*3+1])*(fixed_X[i*3+1]-X[i*3+1]); // fixed energy
energy+=(c-oc)*(fixed_X[i*3+2]-X[i*3+2])*(fixed_X[i*3+2]-X[i*3+2]); // fixed energy
energy*=0.5f;
energy+=-gravity*M[i]*X[i*3+1]; // gravity energy
E[i]+=energy;
float cx=C[i*9+0]+c+ext_C[i];
float cy=C[i*9+4]+c+ext_C[i];
float cz=C[i*9+8]+c+ext_C[i];
// Update Gradient
G[i]=error[0]*error[0]+error[1]*error[1]+error[2]*error[2];
// Update Force
F[i*3+0]=error[0];
F[i*3+1]=error[1];
F[i*3+2]=error[2];
// Update position
next_X[i*3+0]=X[i*3+0]+error[0]*stepping/cx;
next_X[i*3+1]=X[i*3+1]+error[1]*stepping/cy;
next_X[i*3+2]=X[i*3+2]+error[2]*stepping/cz;
}
由上可知，在每个迭代步中，网格节点的位置更新为 next_X = X + error*stepping/cx ，其中，X 为上一个迭代步中网格节点位置，stepping 是步长，即 \(\alpha^{(k)}\) ，error 是优化函数的梯度，即 \(\boldsymbol{g}^{(k)}\) ，cx 应该就是 \(\text{diag}(\boldsymbol{H}^{(k)})\) 。下面将分别详细介绍：
由上述代码可知，X 和 stepping 是位置和步长，容易理解。
error 是优化函数的梯度，即 \(\boldsymbol{g}^{(k)}\) ，具体计算为：
\[\begin{aligned}
- \boldsymbol{g}^{(k)} &= - \nabla\epsilon(\boldsymbol{q}^{(k)}) \\
&= \frac{1}{h^2}(\boldsymbol{q}_{t} + h\boldsymbol{v}_{t} - \boldsymbol{q}^{(k)})\boldsymbol{M} - \nabla E(\boldsymbol{q}^{(k)}) \\
&= \frac{1}{h^2}(\boldsymbol{s}_{t} - \boldsymbol{q}^{(k)})\boldsymbol{M} - \nabla E(\boldsymbol{q}^{(k)})
\end{aligned}\]
其中， \(\boldsymbol{s}_{t} = \boldsymbol{q}_{t} + h\boldsymbol{v}_{t}\) 。具体的计算代码为：
S[i*3+0]=X[i*3+0]+V[i*3+0]*t;
//
float oc = M[i]*inv_t*inv_t;
error[0]=oc*(S[i*3+0]-X[i*3+0])+(c-oc)*(fixed_X[i*3+0]-X[i*3+0])+F[i*3+0];
其中，fixed_X 部分可暂时忽略；F 即为之前计算得到的弹性力。
cx 是 Hessian 矩阵的对角元 \(\text{diag}(\boldsymbol{H}^{(k)})\) ，具体计算为：
\[\text{diag}(\boldsymbol{H}^{(k)}) = \text{diag}(E(\boldsymbol{q}^{(k)})) + \frac{1}{h^2}\boldsymbol{M}
\]
其计算代码为：
float c = M[i]*inv_t*inv_t+fixed[i]+more_fixed[i];
float cx=C[i*9+0]+c+ext_C[i];
同样，此处的 fixed ext_C 等先暂时忽略；C 为之前计算得到的刚度矩阵的对角元。
另外，energy 是另外用于评估收敛的量，这里暂且不讨论。
