论文阅读笔记：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\) 个迭代步中，网格节点的位置更新为

代码实现

这部分代码主要位于 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)}\) ，具体计算为：

其中， \(\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)})\) ，具体计算为：

其计算代码为：

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 是另外用于评估收敛的量，这里暂且不讨论。

本栏目推荐文章
• 代码随想录 day18 找树左下角的值 路径总和 从中序与后序遍历序列构造二叉树
• 用jacoco统计JAVA项目测试代码覆盖率
• （坚持每天写算法）基础算法复习与学习part1基础算法1-7——高精度减法（处理t=1和t>1代码的写法，t为操作次数）
• GDB调试之源代码查看与管理(六)
• 写代码时如何合理的画图表达逻辑结构
• 十行python代码实现文件去重，去除重复文件的脚本
• Java小细节之缩进中if 和else的对应关系
• [代码随想录] 第四天
• idea 中java代码修改后运行代码不生效
• 一行代码创建网页（html）

源代码 simulation 细节 Descent methods源代码simulation细节descent approximation gradient descent method 细节 代码 可读性 细节 代码 细节 过程 代码 问题 接口 铺垫 细节 代码 地雷 细节 代码 方法 ecapa_tdnn细节 代码 论文 descent 爬虫 细节 代码 文件