Stress hessian computation in FEM

介绍

在进行软体模拟时，如果使用牛顿法计算最优的下降方向，需要计算能量密度函数\(\Psi\)关于位置\(\vec {\mathbf{x}}\) 的Hessian矩阵，即\(\frac{\partial^2 \Psi}{\partial \vec {\mathbf{x}}^2}\)。其中\(\vec {\mathbf{x}}\)是一个四面体的四个顶点的位置。

\[ \begin{equation} \frac{\partial^2 \Psi}{\partial {\mathbf{x}^2}} = \text{vec}(\frac{\partial F}{\partial {\mathbf{x}}})^T \text{vec}(\frac{\partial P}{\partial F}) \text{vec}(\frac{\partial F}{\partial {\mathbf{x}}}) \end{equation} \]

涉及到的几个张量的维度如下：

• \({\mathbf{x}} = \begin{bmatrix} \vec x_1 & \vec x_2 & \vec x_3 & \vec x_4 \end{bmatrix} \in \mathbb{R}^{12}\)
• \(D_s = \begin{bmatrix} \vec x_2 - \vec x_1 & \vec x_3 - \vec x_1 & \vec x_4 - \vec x_1 \end{bmatrix} \in \mathbb{R}^{3\times 3}\)
• \(D_m^{-1} = \begin{bmatrix} \vec x_{20} - \vec x_{10} & \vec x_{30} - \vec x_{10} & \vec x_{40} - \vec x_{10} \end{bmatrix}^{-1} \in \mathbb{R}^{3\times 3}\)
• \(F = D_s D_m^{-1} \in \mathbb{R}^{3\times 3}\)
• \(P = \frac{\partial \Psi}{\partial F} \in \mathbb{R}^{3\times 3}\)
• \(\frac{\partial P}{\partial F} \in \mathbb{R}^{3\times 3\times 3\times 3}\)
• \(\frac{\partial F}{\partial \vec {\mathbf{x}}} \in \mathbb{R}^{12\times 3\times 3}\)是一个常数张量

预先计算

\[ \text{newIdm} = \begin{bmatrix} -1 & -1 & -1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} Dm^{-1}. \]

方便直接用其他三个顶点的信息求出\(\vec{x_1}\)的信息。

一些不变量：

• \(I_C = ||F||_F^2 = \text{tr}(F^TF)\)
• \(II_C = ||F^TF||_F^2\)
• \(III_C = \det(F^TF) = J\)

本文参考若干论文，具体说明如何进行\(\frac{\partial \Psi ^2}{\partial \vec {\mathbf{x}}^2}\)的计算。

Robust Quasistatic Finite Elements and Flesh Simulation

本节参考论文Robust Quasistatic Finite Elements的7,8节。

在这篇blog里，已经介绍了如何计算\(\frac{\partial \Psi}{\partial x}\)，令其为\(G\)。

\[ \frac{\partial^2 \Psi}{\partial x^2} = \frac{\partial G}{\partial x} = \frac{\partial F}{\partial x} : \frac{\partial G}{\partial F}. \]

计算步骤

整个计算过程分为四步： 1. 计算\(\frac{\partial P}{\partial F}\)。 2. 计算\(\frac{\partial G}{\partial F}\)。 3. 计算\(\frac{\partial F}{\partial x}\)。 4. 计算\(\frac{\partial G}{\partial x}\)。

计算\(\frac{\partial P}{\partial x}\)

论文提出，可以先用SVD分解\(F\)，得到\(F = U\Sigma V^T\)，然后：

\[ \begin{equation} \delta P = U \left\{\frac{\partial P}{\partial F} \Bigg|_{\hat F} : U^T \delta F V\right\} V^T \end{equation} \]

\[ \text{vec}(\frac{\partial P}{\partial F} \Bigg|_{\hat F}) = \begin{bmatrix} A[1][1] & 0 & 0 & 0 & A[1][2] & 0 & 0 & 0 & A[1][3] \\ 0 & B_{12}[1][1] & 0 & B_{12}[1][2] & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & B_{13}[1][1] & 0 & 0 & 0 & B_{13}[1][2] & 0 & 0 \\ 0 & B_{12}[2][1] & 0 & B_{12}[2][2] & 0 & 0 & 0 & 0 & 0 \\ A[2][1] & 0 & 0 & 0 & A[2][2] & 0 & 0 & 0 & A[2][3] \\ 0 & 0 & 0 & 0 & 0 & B_{23}[1][1] & 0 & B_{23}[1][2] & 0 \\ 0 & 0 & B_{13}[2][1] & 0 & 0 & 0 & B_{13}[2][2] & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & B_{23}[2][1] & 0 & B_{23}[2][2] & 0 \\ A[3][1] & 0 & 0 & 0 & A[3][2] & 0 & 0 & 0 & A[3][3] \end{bmatrix}. \] A,B矩阵的具体形式在论文第8节有详细说明，使用了一些不变量\(I_C, II_C, III_C\)以及\(\frac{\partial \Psi ^2}{\partial (I_C, II_C, III_C)^2}\)。

\[ \frac{\partial \Psi ^2}{\partial (I_C, II_C, III_C)^2} = \begin{bmatrix} \frac{\partial \Psi ^2}{\partial I_C^2} & \frac{\partial \Psi ^2}{\partial II_C \partial I_C} & \frac{\partial \Psi ^2}{\partial III_C \partial I_C} \\ \frac{\partial \Psi ^2}{\partial I_C \partial II_C} & \frac{\partial \Psi ^2}{\partial II_C^2} & \frac{\partial \Psi ^2}{\partial III_C \partial II_C} \\ \frac{\partial \Psi ^2}{\partial I_C \partial III_C} & \frac{\partial \Psi ^2}{\partial II_C \partial III_C} & \frac{\partial \Psi ^2}{\partial III_C^2} \end{bmatrix} \]

计算\(\frac{\partial G}{\partial F}\)

\[ G = PB_m = P(-\text{Vol}D_m^{-T})\\ \frac{\partial G}{\partial F} = -\text{Vol}\frac{\partial P}{\partial F}D_m^{-T} \]

这里需要代入上面的公式1逐个计算\(\frac{\partial G}{\partial F_{ij}}\)。因为G是四面体的力，因此只需要求出\(\vec x_2,\vec x_3,\vec x_4\)上的力，\(\vec x_1\)通过取其他三个顶点的和的负值得到。\(\frac{\partial G}{\partial F_{ij}}\)也同理。

由矩阵微分可知，

\[ \begin{equation} \frac{\partial G}{\partial F_{ij}} = -\text{Vol}\cdot U \left\{\frac{\partial P}{\partial F} \Bigg|_{\hat F} : U^T (e_i \otimes e_j) V\right\} V^T D_m^{-T} \end{equation} \]

双点积可以转换成向量化后矩阵相乘的形式，即：

\[ \frac{\partial G}{\partial x} = \frac{\partial F}{\partial x} : \frac{\partial G}{\partial F} = \frac{\partial F_{jk}}{\partial x_i} e_i \otimes e_j \otimes e_k : \frac{\partial G_n}{\partial F_{lm}} e_l \otimes e_m \otimes e_n = \frac{\partial F_{jk}}{\partial x_i} \frac{\partial G_n}{\partial F_{jk}} e_i \otimes e_n \]

因为\(F\in \mathbb{R}^{3\times 3}\)，所以\(\frac{\partial G}{\partial F_{ij}}\in \mathbb{R}^{12}\)。下面约定按照列主序的方式将\(\frac{\partial G_i}{\partial F}\)展开成一个12维向量：

\[ \frac{\partial G_i}{\partial F} = \begin{bmatrix} \frac{\partial G_i}{\partial F_{11}} & \frac{\partial G_i}{\partial F_{12}} & \frac{\partial G_i}{\partial F_{13}} \\ \frac{\partial G_i}{\partial F_{21}} & \frac{\partial G_i}{\partial F_{22}} & \frac{\partial G_i}{\partial F_{23}} \\ \frac{\partial G_i}{\partial F_{31}} & \frac{\partial G_i}{\partial F_{32}} & \frac{\partial G_i}{\partial F_{33}} \\ \end{bmatrix}\\ \text{vec}(\frac{\partial G_i}{\partial F}) = \begin{bmatrix} \frac{\partial G_{x1x}}{\partial F_{11}}\\ \frac{\partial G_{x1x}}{\partial F_{21}}\\ \frac{\partial G_{x1x}}{\partial F_{31}}\\ \frac{\partial G_{x1x}}{\partial F_{12}}\\ \frac{\partial G_{x1x}}{\partial F_{22}}\\ \frac{\partial G_{x1x}}{\partial F_{32}}\\ \frac{\partial G_{x1x}}{\partial F_{13}}\\ \frac{\partial G_{x1x}}{\partial F_{23}}\\ \frac{\partial G_{x1x}}{\partial F_{33}} \end{bmatrix}\\ \text{vec}(\frac{\partial G}{\partial F}) = \begin{bmatrix} \frac{\partial G_{x1x}}{\partial F_{11}} & \frac{\partial G_{x1y}}{\partial F_{11}} & \frac{\partial G_{x1z}}{\partial F_{11}} & \frac{\partial G_{x2x}}{\partial F_{11}} & \frac{\partial G_{x2y}}{\partial F_{11}} & ... & \frac{\partial G_{x4x}}{\partial F_{11}} & \frac{\partial G_{x4y}}{\partial F_{11}} & \frac{\partial G_{x4z}}{\partial F_{11}} \\ \frac{\partial G_{x1x}}{\partial F_{21}} & \frac{\partial G_{x1y}}{\partial F_{21}} & \frac{\partial G_{x1z}}{\partial F_{21}} & \frac{\partial G_{x2x}}{\partial F_{21}} & \frac{\partial G_{x2y}}{\partial F_{21}} & ... & \frac{\partial G_{x4x}}{\partial F_{21}} & \frac{\partial G_{x4y}}{\partial F_{21}} & \frac{\partial G_{x4z}}{\partial F_{21}} \\ \frac{\partial G_{x1x}}{\partial F_{31}} & \frac{\partial G_{x1y}}{\partial F_{31}} & \frac{\partial G_{x1z}}{\partial F_{31}} & \frac{\partial G_{x2x}}{\partial F_{31}} & \frac{\partial G_{x2y}}{\partial F_{31}} & ... & \frac{\partial G_{x4x}}{\partial F_{31}} & \frac{\partial G_{x4y}}{\partial F_{31}} & \frac{\partial G_{x4z}}{\partial F_{31}} \\ \frac{\partial G_{x1x}}{\partial F_{12}} & \frac{\partial G_{x1y}}{\partial F_{12}} & \frac{\partial G_{x1z}}{\partial F_{12}} & \frac{\partial G_{x2x}}{\partial F_{12}} & \frac{\partial G_{x2y}}{\partial F_{12}} & ... & \frac{\partial G_{x4x}}{\partial F_{12}} & \frac{\partial G_{x4y}}{\partial F_{12}} & \frac{\partial G_{x4z}}{\partial F_{12}} \\ \frac{\partial G_{x1x}}{\partial F_{22}} & \frac{\partial G_{x1y}}{\partial F_{22}} & \frac{\partial G_{x1z}}{\partial F_{22}} & \frac{\partial G_{x2x}}{\partial F_{22}} & \frac{\partial G_{x2y}}{\partial F_{22}} & ... & \frac{\partial G_{x4x}}{\partial F_{22}} & \frac{\partial G_{x4y}}{\partial F_{22}} & \frac{\partial G_{x4z}}{\partial F_{22}} \\ \frac{\partial G_{x1x}}{\partial F_{32}} & \frac{\partial G_{x1y}}{\partial F_{32}} & \frac{\partial G_{x1z}}{\partial F_{32}} & \frac{\partial G_{x2x}}{\partial F_{32}} & \frac{\partial G_{x2y}}{\partial F_{32}} & ... & \frac{\partial G_{x4x}}{\partial F_{32}} & \frac{\partial G_{x4y}}{\partial F_{32}} & \frac{\partial G_{x4z}}{\partial F_{32}} \\ \frac{\partial G_{x1x}}{\partial F_{13}} & \frac{\partial G_{x1y}}{\partial F_{13}} & \frac{\partial G_{x1z}}{\partial F_{13}} & \frac{\partial G_{x2x}}{\partial F_{13}} & \frac{\partial G_{x2y}}{\partial F_{13}} & ... & \frac{\partial G_{x4x}}{\partial F_{13}} & \frac{\partial G_{x4y}}{\partial F_{13}} & \frac{\partial G_{x4z}}{\partial F_{13}} \\ \frac{\partial G_{x1x}}{\partial F_{23}} & \frac{\partial G_{x1y}}{\partial F_{23}} & \frac{\partial G_{x1z}}{\partial F_{23}} & \frac{\partial G_{x2x}}{\partial F_{23}} & \frac{\partial G_{x2y}}{\partial F_{23}} & ... & \frac{\partial G_{x4x}}{\partial F_{23}} & \frac{\partial G_{x4y}}{\partial F_{23}} & \frac{\partial G_{x4z}}{\partial F_{23}} \\ \frac{\partial G_{x1x}}{\partial F_{33}} & \frac{\partial G_{x1y}}{\partial F_{33}} & \frac{\partial G_{x1z}}{\partial F_{33}} & \frac{\partial G_{x2x}}{\partial F_{33}} & \frac{\partial G_{x2y}}{\partial F_{33}} & ... & \frac{\partial G_{x4x}}{\partial F_{33}} & \frac{\partial G_{x4y}}{\partial F_{33}} & \frac{\partial G_{x4z}}{\partial F_{33}} \end{bmatrix}_{9\times 12}\\ \]

这样每一行即上面提到的逐个计算出的\(\frac{\partial G}{\partial F_{ij}}\)。因为这个计算过程只会计算出\(\vec x_2,\vec x_3,\vec x_4\)上的导数，所以\(\vec x_1\)上的导数可以通过取其他三个顶点的和的负值得到：

\[ \frac{\partial G_{x1x}}{\partial F_{ij}} + \frac{\partial G_{x2x}}{\partial F_{ij}} + \frac{\partial G_{x3x}}{\partial F_{ij}} + \frac{\partial G_{x4x}}{\partial F_{ij}} = 0\\ \frac{\partial G_{x1y}}{\partial F_{ij}} + \frac{\partial G_{x2y}}{\partial F_{ij}} + \frac{\partial G_{x3y}}{\partial F_{ij}} + \frac{\partial G_{x4y}}{\partial F_{ij}} = 0\\ \frac{\partial G_{x1z}}{\partial F_{ij}} + \frac{\partial G_{x2z}}{\partial F_{ij}} + \frac{\partial G_{x3z}}{\partial F_{ij}} + \frac{\partial G_{x4z}}{\partial F_{ij}} = 0\\ \]

计算\(\frac{\partial F}{\partial x}\)

这个在计算梯度的时候已经计算过了，它是：

\[ \text{vec}(\frac{\partial F}{\partial x})^T = \text{newIdm} \otimes I_{3\times 3} \]

计算\(\frac{\partial G}{\partial x}\)

\[ \frac{\partial G}{\partial x} = \text{vec}(\frac{\partial F}{\partial x})^T \text{vec}(\frac{\partial G}{\partial F}) \]

取第i列：

\[ \frac{\partial G_i}{\partial x} = \text{vec}(\frac{\partial F}{\partial x})^T \text{vec}(\frac{\partial G_i}{\partial F}) \\ =\text{newIdm} \otimes I_{3\times 3} \text{vec}(\frac{\partial G_i}{\partial F})\\ \]

由Kronecker积的性质:

给定矩阵\(A,B,X\)：

\[ B^T \otimes A \cdot \text{vec}(X) = \text{vec}(AXB) \]

令\(B^T = \begin{bmatrix} -1 & -1 & -1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}D_m^{-1}, X = \frac{\partial G_i}{\partial F}, A = I_{3\times 3}\)，则：

\[ \begin{equation} \frac{\partial G_i}{\partial x} = \text{vec}(\frac{\partial G_i}{\partial F}(\text{newIdm}\cdot D_m^{-1})^T) \end{equation} \]

代码实现

C++

这一节以论文Descent Methods for Elastic Body Simulation on the GPU提供的代码中的Compute_FM_Kernel函数为例说明上面的计算过程，里面后半段代码(if(update_C==false) return;之后)计算了hessian矩阵。

下面的代码计算了\(\text{vec}(\frac{\partial P}{\partial F} \Bigg|_{\hat F})\)里的alpha,beta,gamma。

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]);
   }

下面的代码用公式2计算了\(\frac{\partial G}{\partial F}\)。遍历\(F\)的9个位置，比如i=1时取\(\delta F = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\)。

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]);
}

因为这样求出的每一个dGdF[i]是\(\frac{\partial G}{\partial F_{ij}}\)，但是后面要取出\(\frac{\partial G_i}{\partial F}\)，所以需要转置一下，顺便把\(\vec x_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];
}

最后根据公式3计算hessian。先计算\(\text{newIdm}\)。 刚刚通过转置已经让dGdF[i]=\(\frac{\partial G_i}{\partial F}\)，把他当成3x3矩阵计算乘法，再flatten成向量即得到hessian的第i列。这段代码因为只需要hessian的对角部分，所以直接取了对角元来组装对角向量。需要注意的是，因为这个代码项目使用行主序存储矩阵，因此-dev_Matrix_Product_T(&new_idm[0][0], &dGdF[9][0], 4, 3, 3, 0, 0)计算了\(\frac{\partial G_i}{\partial F}(\text{newIdm})^T\)的转置。这样取出的对应元素才是对角元。

   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+0], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[9][0], 4, 3, 3, 0, 0));
atomicAdd(&C[p0+1], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[10][0], 4, 3, 3, 0, 1));
atomicAdd(&C[p0+2], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[11][0], 4, 3, 3, 0, 2));
atomicAdd(&C[p1+0], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[0][0], 4, 3, 3, 1, 0));
atomicAdd(&C[p1+1], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[3][0], 4, 3, 3, 1, 1));
atomicAdd(&C[p1+2], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[6][0], 4, 3, 3, 1, 2));
atomicAdd(&C[p2+0], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[1][0], 4, 3, 3, 2, 0));
atomicAdd(&C[p2+1], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[4][0], 4, 3, 3, 2, 1));
atomicAdd(&C[p2+2], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[7][0], 4, 3, 3, 2, 2));
atomicAdd(&C[p3+0], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[2][0], 4, 3, 3, 3, 0));
atomicAdd(&C[p3+1], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[5][0], 4, 3, 3, 3, 1));
atomicAdd(&C[p3+2], -dev_Matrix_Product_T(&new_idm[0][0], &dGdF[8][0], 4, 3, 3, 3, 2));

Matlab

下面的代码略去了\(\frac{\partial G}{\partial F}\)的计算，用三种方法计算了hessian矩阵。一种是原始的张量双点积(dGdx)，一种是向量化后的大矩阵乘法(dGdx2)，一种是类似上面代码的逐行求出(dGdx3)。

pgpf = sym('pgpf', [3, 3, 12]);
x = sym('x', [12, 1]);

Ds = [x(4)-x(1) x(7)-x(1) x(10)-x(1); x(5)-x(2) x(8)-x(2) x(11)-x(2); x(6)-x(3) x(9)-x(3) x(12)-x(3)];
PDsPx = sym(zeros(12, 3, 3));
for i = 1:3
    for j = 1:3
        for k = 1:12
            PDsPx(k, i, j) = diff(Ds(i, j), x(k));
        end
    end
end

dGdF129 = reshape(pgpf(:),9,12);
dGdF129 = reshape(PDsPx(:),12,9) * dGdF129(1:9,1:9);

dGdF = sym(zeros(3, 3, 12));
for i = 1:3
    for j = 1:3
        dGdF(j, i, :) = dGdF129(:,(i - 1) * 3 + j);
    end
end

DmInv = sym('DmInv', [3, 3]);
dFdx = tensorDot(PDsPx, DmInv);
dGdx = doubleDotProd(dFdx, dGdF);
dGdx2 = sym(zeros(12, 12));
new_idm = [-1 -1 -1; 1 0 0; 0 1 0; 0 0 1] * DmInv;
for j = 1:12
    dGdx2(:,j)= reshape(reshape(dGdF129(j,:),3,3)* transpose(new_idm),12,1);
end
dGdx3 = kron(new_idm, eye(3)) * transpose(dGdF129);
disp(dGdx-dGdx3);
disp(dGdx-dGdx2);

tensorDot和doubleDotProd代码在这里。

Nonlinear Material Design Using Principal Stretches

本节参考论文Nonlinear Material Design Using Principal Stretches的3.2节。

Neo-Hookean模型

以Bonet and Wood (2008)-style Neo-Hookean为例：

\[ \Psi(F)_{\text{BW08}} = \frac{\mu}{2}(I_C-3) - \mu \ln J + \frac{\lambda}{2}(\ln J)^2 \]

可以对\(F\)进行Polar SVD分解(翻转\(\Sigma\)符号使得\(U,V\)行列式为正)：

\[ F = U\Sigma V^T \]

则\(\Psi\)可以有以下形式：

\[ \begin{aligned} \Psi(F)_{\text{BW08}} &=\hat\Psi(\Sigma)_{\text{BW08}} \\ &=\frac{\mu}{2}(\sum_{i=1}^3 \sigma_i^2 - 3) - \mu \ln(\sigma_1\sigma_2\sigma_3) + \frac{\lambda}{2}(\ln(\sigma_1\sigma_2\sigma_3))^2 \end{aligned} \]

那么：

\[ P = U \hat P V^T \]

其中：

\[ \hat P = \begin{bmatrix} \frac{\partial \hat\Psi}{\partial \sigma_{1}} & 0 & 0 \\ 0 & \frac{\partial \hat\Psi}{\partial \sigma_{2}} & 0 \\ 0 & 0 & \frac{\partial \hat\Psi}{\partial \sigma_{3}} \end{bmatrix} \]

对于Neo-Hookean模型：

\[ \hat P_{ii} = \mu(\sigma_i-\frac{1}{\sigma_i}) + \lambda \ln(\sigma_1\sigma_2\sigma_3)\frac{1}{\sigma_i} \]

接下来计算\(\frac{\partial P}{\partial F}\)。 \[ \begin{aligned} \frac{\partial P}{\partial F_{ij}} &= \frac{\partial U}{\partial F_{ij}} \hat P V^T + U \frac{\partial \hat P}{\partial \sigma_{k}} \frac{\partial \sigma_{k}}{\partial F_{ij}} V^T + U \hat P \frac{\partial V^T}{\partial F_{ij}}\\ \end{aligned} \]

其中有三个未知量，分别是\(\frac{\partial U}{\partial F_{ij}}\)，\(\frac{\partial V^T}{\partial F_{ij}}\)和\(\frac{\partial \sigma_{k}}{\partial F_{ij}}\)。根据论文Nonlinear Material Design Using Principal Stretches的3.2节，可以使用\(\frac{\partial F}{\partial F_{ij}}\)来计算这三个未知量。

\[ \begin{aligned} \frac{\partial F}{\partial F_{ij}} &= \frac{\partial U}{\partial F_{ij}} \Sigma V^T + U \frac{\partial \Sigma}{\partial F_{ij}} V^T + U \Sigma \frac{\partial V^T}{\partial F_{ij}}\\ U^T \frac{\partial F}{\partial F_{ij}} V &= U^T \frac{\partial U}{\partial F_{ij}} \Sigma + \frac{\partial \Sigma}{\partial F_{ij}} + \Sigma \frac{\partial V^T}{\partial F_{ij}} V\\ \end{aligned} \]

已知\(U^T U = I\)，\(V^T V = I\)，有：

\[ \begin{aligned} \frac{\partial U^T U}{\partial F_{ij}} & = \frac{\partial U^T}{\partial F_{ij}} U + U^T \frac{\partial U}{\partial F_{ij}} \\ & = (U^T \frac{\partial U}{\partial F_{ij}} )^T + U^T \frac{\partial U}{\partial F_{ij}} = 0\\ \frac{\partial V^T V}{\partial F_{ij}} & = \frac{\partial V^T}{\partial F_{ij}} V + V^T \frac{\partial V}{\partial F_{ij}}\\ & = \frac{\partial V^T}{\partial F_{ij}} V + (V^T \frac{\partial V}{\partial F_{ij}})^T = 0\\ \end{aligned} \]

所以，\(\frac{\partial U}{\partial F_{ij}}\)，\(\frac{\partial V^T}{\partial F_{ij}}\)都是反对称矩阵。对角线元素都为0。假设 \[ \begin{aligned} \Sigma &= \begin{bmatrix} \sigma_1 & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 & \sigma_3 \end{bmatrix}\\ U^T\frac{\partial F}{\partial F_{ij}} V&= \begin{bmatrix} f_{11} & f_{12} & f_{13} \\ f_{21} & f_{22} & f_{23} \\ f_{31} & f_{32} & f_{33} \end{bmatrix},\\ U^T\frac{\partial U}{\partial F_{ij}} &= \begin{bmatrix} 0 & u_{12} & u_{13} \\ -u_{12} & 0 & u_{23} \\ -u_{13} & -u_{23} & 0 \end{bmatrix},\\ \frac{\partial V^T}{\partial F_{ij}}V &= \begin{bmatrix} 0 & v_{12} & v_{13} \\ -v_{12} & 0 & v_{23} \\ -v_{13} & -v_{23} & 0 \end{bmatrix} \end{aligned} \]

则：

\[ \begin{aligned} U^T \frac{\partial U}{\partial F_{ij}} \Sigma + \Sigma \frac{\partial V^T}{\partial F_{ij}} V &= \begin{bmatrix} 0 & \sigma_2 u_{12} + \sigma_1 v_{12} & \sigma_3 u_{13} + \sigma_1 v_{13} \\ -\sigma_1 u_{12} - \sigma_2 v_{12} & 0 & \sigma_3 u_{23} + \sigma_2 v_{23} \\ -\sigma_1 u_{13} - \sigma_3 v_{13} & -\sigma_2 u_{23} - \sigma_3 v_{23} & 0 \end{bmatrix} \end{aligned} \]

可以看出，\(U^T \frac{\partial U}{\partial F_{ij}} \Sigma + \Sigma \frac{\partial V^T}{\partial F_{ij}} V\)对\(\frac{\partial F}{\partial F_{ij}}\)的对角元没有贡献，所以： \[ \frac{\partial \Sigma}{\partial F_{ij}} = \begin{bmatrix} f_{11} & 0 & 0 \\ 0 & f_{22} & 0 \\ 0 & 0 & f_{33} \end{bmatrix} \]

即\(\frac{\partial \sigma_{k}}{\partial F_{ij}} = f_{kk}\)，\(k=1,2,3\)。以及

\[ \begin{bmatrix} 0 & \sigma_2 u_{12} + \sigma_1 v_{12} & \sigma_3 u_{13} + \sigma_1 v_{13} \\ -\sigma_1 u_{12} - \sigma_2 v_{12} & 0 & \sigma_3 u_{23} + \sigma_2 v_{23} \\ -\sigma_1 u_{13} - \sigma_3 v_{13} & -\sigma_2 u_{23} - \sigma_3 v_{23} & 0 \end{bmatrix} = \begin{bmatrix} 0 & f_{12} & f_{13} \\ f_{21} & 0 & f_{23} \\ f_{31} & f_{32} & 0 \end{bmatrix} \] 三组对称位置的元素分别可以构成一个二元一次方程组。

\[ \begin{aligned} \begin{bmatrix} \sigma_2 & \sigma_1 \\ \sigma_1 & \sigma_2 \end{bmatrix} \begin{bmatrix} u_{12} \\ v_{12} \end{bmatrix} &= \begin{bmatrix} f_{12} \\ -f_{21} \end{bmatrix},\\ \begin{bmatrix} \sigma_3 & \sigma_1 \\ \sigma_1 & \sigma_3 \end{bmatrix} \begin{bmatrix} u_{13} \\ v_{13} \end{bmatrix} &= \begin{bmatrix} f_{13} \\ -f_{31} \end{bmatrix},\\ \begin{bmatrix} \sigma_3 & \sigma_2 \\ \sigma_2 & \sigma_3 \end{bmatrix} \begin{bmatrix} u_{23} \\ v_{23} \end{bmatrix} &= \begin{bmatrix} f_{23} \\ -f_{32} \end{bmatrix} \end{aligned} \]

总结一下计算\(\frac{\partial^2 \Psi}{\partial \vec{x}^2}\)的步骤：

1. 计算\(F = U\Sigma V^T\)，得到\(\sigma_1,\sigma_2,\sigma_3,U,V\).
2. 计算\(\hat P\).
3. 遍历\(i,j\)，
   1. 计算\(\frac{\partial F}{\partial F_{ij}}\)，取\(\frac{\partial F}{\partial F_{ij}}\)的对角元素作为\(\frac{\partial \sigma_{k}}{\partial F_{ij}}\)，\(k=1,2,3\).
   2. 用其他元素代入上面的三个方程组，解出\(u_{12},v_{12},u_{13},v_{13},u_{23},v_{23}\)，即得到\(\frac{\partial U}{\partial F_{ij}}\)和\(\frac{\partial V^T}{\partial F_{ij}}\).
   3. 使用公式\(\frac{\partial P}{\partial F_{ij}} = \frac{\partial U}{\partial F_{ij}} \hat P V^T + U \frac{\partial \hat P}{\partial \sigma_{k}} \frac{\partial \sigma_{k}}{\partial F_{ij}} V^T + U \hat P \frac{\partial V^T}{\partial F_{ij}}, k=1,2,3\).
4. \(\frac{\partial^2 \Psi}{\partial \vec{x}^2} = \text{vec}^T(\frac{\partial F}{\partial \vec{x}}) \text{vec}(\frac{\partial P}{\partial F}) \text{vec}(\frac{\partial F}{\partial \vec{x}})\)

Corotated linear elasticity

\[ \begin{aligned} F &= U\Sigma V^T\\ R &= UV^T, S = V\Sigma V^T\\ \Psi(F) &= \mu ||F-R||_F^2 + \frac{\lambda}{2}(\text{tr}^2(R^TF-I))\\ P &= \frac{\partial \Psi}{\partial F} = 2\mu(F-R) + \lambda \text{tr}(R^TF-I)R\\ \frac{\partial P}{\partial F} &= 2\mu (I - \frac{\partial R}{\partial F}) + \lambda (R \otimes R + \text{tr}(R^TF-I)\frac{\partial R}{\partial F})\\ &=2\mu (I - \frac{\partial R}{\partial F}) + \lambda [F \otimes R + \text{tr}(R^TF-I)\frac{\partial R}{\partial R}] : \frac{\partial R}{\partial F}\\ \end{aligned} \]

其中用到了\(\frac{\partial R}{\partial F}\)，因为

\[ \frac{\partial R}{\partial F_{ij}} = \frac{\partial U}{\partial F_{ij}} V^T + U \frac{\partial V^T}{\partial F_{ij}} \]

上一节的推导中可以通过\(\frac{\partial F}{\partial F_{ij}}\)计算\(\frac{\partial U}{\partial F_{ij}}\)和\(\frac{\partial V^T}{\partial F_{ij}}\)，所以可以直接计算\(\frac{\partial R}{\partial F}\)。

使用上面公式计算\(\frac{\partial P}{\partial F}\)之后，计算Hessian矩阵的方法就和上一节一样了。

A Better Way For Isotropic Solids

这一节参考Siggraph courseDynamic deformables: implementation and production practicalities (now with code!)的5.5节。 一些符号定义：

\[ F = \begin{bmatrix} f_0 \Bigg| f_1 \Bigg| f_2 \end{bmatrix}=U\Sigma V^T=RS\\ S = V\Sigma V^T\\ R = UV^T\\ J=\text{det}(F)\\ \hat x = \begin{bmatrix} 0 & -x_2 & x_1 \\ x_2 & 0 & -x_0 \\ -x_1 & x_0 & 0 \end{bmatrix}\\ \]

定义三个新的不变量：

\[ \begin{aligned} I_1 &= \text{tr}(S)=\text{tr}(\Sigma)\\ I_2 &= \text{tr}(F^TF)=\text{tr}(S^T S)\\ I_3 &= J = \text{det}(F),\\ \end{aligned} \]

他们的梯度和Hessian矩阵分别是：

\[ \begin{aligned} g_1 &= \text{vec}(R)\\ g_2 &= \text{vec}(2F)\\ g_{J} &= \text{vec}(\begin{bmatrix} f1 \times f2 \Bigg| f2 \times f0 \Bigg| f0 \times f1 \end{bmatrix})\\ H_1 &= \sum_{i=1}^3 \lambda_i \text{vec}(Q_i) \text{vec}(Q_i)^T\\ H_2 &= 2 I_{9\times 9}\\ H_{J} &= \begin{bmatrix} 0 & -\hat{f}_2 & \hat{f}_1 \\ \hat{f}_2 & 0 & -\hat{f}_0 \\ -\hat{f}_1 & \hat{f}_0 & 0 \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 & 0 & F_{22} & -F_{12} & 0 & -F_{21} & F_{11} \\ 0 & 0 & 0 & -F_{22} & 0 & F_{02} & F_{21} & 0 & -F_{01} \\ 0 & 0 & 0 & F_{12} & -F_{02} & 0 & -F_{11} & F_{01} & 0 \\ 0 & -F_{22} & F_{12} & 0 & 0 & 0 & 0 & F_{20} & -F_{10} \\ F_{22} & 0 & -F_{02} & 0 & 0 & 0 & -F_{20} & 0 & F_{00} \\ -F_{12} & F_{02} & 0 & 0 & 0 & 0 & F_{10} & -F_{00} & 0 \\ 0 & F_{21} & -F_{11} & 0 & -F_{20} & F_{10} & 0 & 0 & 0 \\ -F_{21} & 0 & F_{01} & F_{20} & 0 & -F_{00} & 0 & 0 & 0 \\ F_{11} & -F_{01} & 0 & -F_{10} & F_{00} & 0 & 0 & 0 & 0 \end{bmatrix} \end{aligned} \] 其中： \[ \lambda_0 = \frac{2}{\sigma_x + \sigma_y}\\ \lambda_1 = \frac{2}{\sigma_y + \sigma_z}\\ \lambda_2 = \frac{2}{\sigma_z + \sigma_x}\\ Q_0 = \frac{1}{\sqrt{2}}U\begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}V^T\\ Q_1 = \frac{1}{\sqrt{2}}U\begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & -1 & 0 \end{bmatrix}V^T\\ Q_2 = \frac{1}{\sqrt{2}}U\begin{bmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \end{bmatrix}V^T\\ \] 这样可以得到\(\text{vec}(\frac{\partial P}{\partial F})\)：

\[ \begin{aligned} \text{vec}(\frac{\partial^2 \Psi}{\partial F^2}) &= \frac{\partial^2 \Psi}{\partial I_1^2} g_1 g_1^T + \frac{\partial \Psi}{\partial I_1} H_1 + \frac{\partial^2 \Psi}{\partial I_2^2} g_2 g_2^T + \frac{\partial \Psi}{\partial I_2} H_2 + \frac{\partial^2 \Psi}{\partial I_3^2} g_3 g_3^T + \frac{\partial \Psi}{\partial I_3} H_3\\ &= \sum_{i=1}^3 \frac{\partial^2 \Psi}{\partial I_i^2} g_i g_i^T + \frac{\partial \Psi}{\partial I_i} H_i \end{aligned} \]

最后计算hessian：

\[ \begin{aligned} \frac{\partial^2 \Psi}{\partial \vec{x}^2} &= \text{vec}(\frac{\partial F}{\partial \vec{x}})^T \text{vec}(\frac{\partial P}{\partial F}) \text{vec}(\frac{\partial F}{\partial \vec{x}})\\ &=\text{newIdm} \otimes I_{3\times 3} \cdot \text{vec}(\frac{\partial P}{\partial F}) \cdot (\text{newIdm})^T \otimes I_{3\times 3} \end{aligned} \]

依然可以用Kronecker积的性质简化计算。将\(\text{vec}(\frac{\partial P}{\partial F})\)写成9个列向量： \[ \text{vec}(\frac{\partial P}{\partial F}) = \begin{bmatrix} \text{vec}(P_0) & \text{vec}(P_1) & \text{vec}(P_2) & ... &\text{vec}(P_7) & \text{vec}(P_8) \end{bmatrix} \]

创建一个12x9的矩阵tmp，每一列等于\(\text{vec}(P_i \cdot\text{newIdm})\)。

再创建一个9x9的矩阵，每一行等于tmp的一行reshape成3x3矩阵再乘以\(\text{newIdm}^T\)。因为hessian矩阵是对称的，所以这个9x9矩阵就是hessian。这样只需要\((9+12)\times (3\times 3 \times 4)\)次乘法。直接计算hessian需要\(12\times 9 \times 9 + 12 \times 12 \times 9\)次乘法。

计算hessian的matlab代码：

dPdF = sym('dpdf', [3, 3, 3, 3]);
vec_dPdF = reshape(dPdF,9,9);
d2Psidx2 = sym(zeros(12,12));
for i = 1:9
    d2Psidx22(:,i)=reshape(reshape(vec_dPdF(:,i),3,3)*transpose(new_idm),12,1);
end
for i = 1:12
    d2Psidx22(i,:)=reshape(reshape(d2Psidx22(i,1:9),3,3)*transpose(new_idm),12,1);
end

计算hessian的C++代码，里面可以直接用\(D_m^{-1}\)计算然后手动算出其余项：

template <typename Scalar>
__host__ __device__ void ComputeHessian(const Scalar* DmInv, const Matrix9<Scalar>& d2PsidF2, Matrix12<Scalar>& H) {
    for (size_t i = 0; i < 9; i++)
    {
        H[3][i] = d2PsidF2[0][i] * DmInv[0] + d2PsidF2[3][i] * DmInv[3] + d2PsidF2[6][i] * DmInv[6];
        H[4][i] = d2PsidF2[1][i] * DmInv[0] + d2PsidF2[4][i] * DmInv[3] + d2PsidF2[7][i] * DmInv[6];
        H[5][i] = d2PsidF2[2][i] * DmInv[0] + d2PsidF2[5][i] * DmInv[3] + d2PsidF2[8][i] * DmInv[6];;
        H[6][i] = d2PsidF2[0][i] * DmInv[1] + d2PsidF2[3][i] * DmInv[4] + d2PsidF2[6][i] * DmInv[7];
        H[7][i] = d2PsidF2[1][i] * DmInv[1] + d2PsidF2[4][i] * DmInv[4] + d2PsidF2[7][i] * DmInv[7];
        H[8][i] = d2PsidF2[2][i] * DmInv[1] + d2PsidF2[5][i] * DmInv[4] + d2PsidF2[8][i] * DmInv[7];
        H[9][i] = d2PsidF2[0][i] * DmInv[2] + d2PsidF2[3][i] * DmInv[5] + d2PsidF2[6][i] * DmInv[8];
        H[10][i] = d2PsidF2[1][i] * DmInv[2] + d2PsidF2[4][i] * DmInv[5] + d2PsidF2[7][i] * DmInv[8];
        H[11][i] = d2PsidF2[2][i] * DmInv[2] + d2PsidF2[5][i] * DmInv[5] + d2PsidF2[8][i] * DmInv[8];
        H[0][i] = -H[3][i] - H[6][i] - H[9][i];
        H[1][i] = -H[4][i] - H[7][i] - H[10][i];
        H[2][i] = -H[5][i] - H[8][i] - H[11][i];
    }
    Scalar temp[9];
    for (size_t i = 0; i < 12; i++)
    {
        memcpy(temp, H[i].data(), 9 * sizeof(Scalar));
        H[i][3] = temp[0] * DmInv[0] + temp[3] * DmInv[3] + temp[6] * DmInv[6];
        H[i][4] = temp[1] * DmInv[0] + temp[4] * DmInv[3] + temp[7] * DmInv[6];
        H[i][5] = temp[2] * DmInv[0] + temp[5] * DmInv[3] + temp[8] * DmInv[6];
        H[i][6] = temp[0] * DmInv[1] + temp[3] * DmInv[4] + temp[6] * DmInv[7];
        H[i][7] = temp[1] * DmInv[1] + temp[4] * DmInv[4] + temp[7] * DmInv[7];
        H[i][8] = temp[2] * DmInv[1] + temp[5] * DmInv[4] + temp[8] * DmInv[7];
        H[i][9] = temp[0] * DmInv[2] + temp[3] * DmInv[5] + temp[6] * DmInv[8];
        H[i][10] = temp[1] * DmInv[2] + temp[4] * DmInv[5] + temp[7] * DmInv[8];
        H[i][11] = temp[2] * DmInv[2] + temp[5] * DmInv[5] + temp[8] * DmInv[8];
        H[i][0] = -H[i][3] - H[i][6] - H[i][9];
        H[i][1] = -H[i][4] - H[i][7] - H[i][10];
        H[i][2] = -H[i][5] - H[i][8] - H[i][11];
    }
}

下面是课程里提到的例子：

Neo-Hookean模型

Bonet and Wood (2008)-style Neo-Hookean：

\[ \begin{aligned} \Psi(F)_{\text{BW08}} &= \frac{\mu}{2}(||F||_F^2-3) - \mu \ln J + \frac{\lambda}{2}(\ln J)^2\\ &= \frac{\mu}{2}(I_2-3) - \mu \ln I_3 + \frac{\lambda}{2}(\ln I_3)^2\\ \frac{\partial \Psi}{\partial I_1} &= 0, \frac{\partial^2 \Psi}{\partial I_1^2} = 0\\ \frac{\partial \Psi}{\partial I_2} &= \frac{\mu}{2}, \frac{\partial^2 \Psi}{\partial I_2^2} = 0\\ \frac{\partial \Psi}{\partial I_3} &= \frac{\lambda \text{ln} I_3 - \mu}{I_3}, \frac{\partial^2 \Psi}{\partial I_3^2} = \frac{\lambda(1-\text{ln}I_3)+\mu}{I_3^2}\\ \text{vec}(\frac{\partial P}{\partial F}) &= \mu I_{9\times 9} +\frac{\lambda(1-\text{ln}I_3)+\mu}{I_3^2} g_3 g_3^T + \frac{\lambda \text{ln} I_3 - \mu}{I_3} H_3\\ \end{aligned}\\ \]

ARAP模型

\[ \begin{aligned} \Psi(F)_{\text{ARAP}} &= ||F-R||_F^2 \\ &= I_2 - 2I_1 + 3\\ \frac{\partial \Psi}{\partial I_1} &= -2, \frac{\partial^2 \Psi}{\partial I_1^2} = 0\\ \frac{\partial \Psi}{\partial I_2} &= 1, \frac{\partial^2 \Psi}{\partial I_2^2} = 0\\ \frac{\partial \Psi}{\partial I_3} &= 0, \frac{\partial^2 \Psi}{\partial I_3^2} = 0\\ \text{vec}(\frac{\partial P}{\partial F}) &= 2 I_{9\times 9} - 2 H_1 \end{aligned} \]