Projective Dynamics 的local step实现推导

介绍

Projective Dynamics是一种用于软体模拟的方法，算法分为local step和global step两个部分。其中local step可以对于每个四面体约束并行计算，global step只需要求解一个线性方程组，而他的矩阵非常特殊，是一个Gram矩阵，因此可以预先用Cholesky分解。本文主要介绍local step的实现的推导过程。用Corotated strain model的CUDA的实现作为例子。

背景

这部分内容参考SIGGRAPH 2012 Course FEM Simulation of 3D Deformable Solids。

能量密度函数(Energy Density Function)

\(\Psi[\Phi;\vec X]\)是一个能量密度函数，用来衡量在\(\vec X\)附近的一个无穷小体积\(dV\)的应变能。

Corotated linear elasticity

Corotated linear elasticity is a constitutive model that attempts to combine the simplicity of the stress-deformation relationship in a linear material with just enough nonlinear characteristics to secure rotational invariance.

对于corotated linear elasticity，\(\Psi\)可以写成： \[ \Psi[F]=\mu||F-R||^2_F+\frac{\lambda}{2}(\text{tr}^2(R^TF)-I) \] 对于四面体的四个顶点\(\vec q_1,\vec q_2,\vec q_3,\vec q_4\)，在任意时刻， \[ F=\begin{bmatrix}\vec q_1-\vec q_4 & \vec q_2-\vec q_4 & \vec q_3-\vec q_4\end{bmatrix}D_m^{-1} \] 其中 \[ D_m = \begin{bmatrix}\vec q_{10}-\vec q_{40} & \vec q_{20}-\vec q_{40} & \vec q_{30}-\vec q_{40}\end{bmatrix} \] 是rest shape时计算得到的。一般这个会预先算好。

Projective Dynamics

Local step

\[ \min_{p_i} \frac{w_i}{2}||A_i S_i q -B_ip_i||^2_F+\delta_{C_i}(p_i) \]

Global step

\[ (\frac{M}{h^2}+\sum_{i=1}^n w_i S_i^T A_i^T A_i S_i)q=\frac{M}{h^2}s_n+\sum_{i=1}^n w_i S_i^T A_i^T B_ip_i \]

推导

根据论文，local step只需要选择想要最小化的约束，将投影结果算出后累加到每个顶点上即可。这里选择最小化能量密度函数，即： \[ \min_{p_i} \frac{w_i}{2}||F_i^T-R_i^T||^2_F+\delta_{C_i}(p_i) \] 从上面的公式可以看出，\(B_i\)是一个单位矩阵。

转置是为了方便凑出\(S_i^TA_i^Tp_i\)，这样就能将\(A_iS_i q\)对应到\(F^T\)，\(p_i\)对应到\(R^T\)。

接下来只需要凑出\(S_i^TA_i^Tp_i\)即可。在我们的实现里，\(q\)是\(\mathbb{R}^{n\times 3}\)的矩阵（但是存储方式是flatten成\(\mathbb{R}^{3n}\)的）。现在只考虑局部的一个四面体的四个顶点，因此我们把\(q\)当做一个\(\mathbb{R}^{4\times 3}\)的矩阵： \[ q=\begin{bmatrix}\vec q_1^T\\\vec q_2^T\\\vec q_3^T\\\vec q_4^T\end{bmatrix} \]

首先凑出\(F\)，根据背景部分的公式： \[ F= (\begin{bmatrix} 1 & 0 & 0 & -1 \\ 0 & 1 & 0 & -1 \\ 0 & 0 & 1 & -1 \end{bmatrix} q)^T D_m^{-1} \]

\[ F^T= D_m^{-T}\begin{bmatrix} 1 & 0 & 0 & -1 \\ 0 & 1 & 0 & -1 \\ 0 & 0 & 1 & -1 \end{bmatrix} q \] 这样就得到： \[ A_i= D_m^{-T},\\ S_i= \begin{bmatrix} 1 & 0 & 0 & -1 \\ 0 & 1 & 0 & -1 \\ 0 & 0 & 1 & -1 \end{bmatrix} \]

然后计算\(R^T\)，先对\(F\)做SVD分解，得到\(F=U\Sigma V^T\)，然后\(R=UV^T\)。于是 \[ S_i^TA_i^Tp_i =(D_m^{-T}\begin{bmatrix} 1 & 0 & 0 & -1 \\ 0 & 1 & 0 & -1 \\ 0 & 0 & 1 & -1 \end{bmatrix} )^T(UV^T)^T \] 因为这一项是一个\(\mathbb{R}^{4\times 3}\)的矩阵，所以我们需要将其flatten成\(\mathbb{R}^{12}\)的向量，然后累加到每个顶点上。实际操作中就计算\(S_i^TA_i^Tp_i\)的转置\(\mathbb{R}^{3\times4}\)，然后将每个列向量累加到对应的顶点上即可，也就是说最后得到的是一个\(\mathbb{R}^{3\times4}\)的矩阵： \[ UV^TD_m^{-T}\begin{bmatrix} 1 & 0 & 0 & -1 \\ 0 & 1 & 0 & -1 \\ 0 & 0 & 1 & -1 \end{bmatrix} \]

实现：

__global__ void solveLocal(const float* V0, const float wi, float* xProj, const glm::mat3* DmInvs, const float* qn, const indexType* tets, int numTets)
{
    int index = (blockIdx.x * blockDim.x) + threadIdx.x;
    if (index < numTets)
    {
        const int v0Ind = tets[index * 4 + 0] * 3;
        const int v1Ind = tets[index * 4 + 1] * 3;
        const int v2Ind = tets[index * 4 + 2] * 3;
        const int v3Ind = tets[index * 4 + 3] * 3;
        const glm::vec3 v0 = glm::vec3(qn[v0Ind + 0], qn[v0Ind + 1], qn[v0Ind + 2]);
        const glm::vec3 v1 = glm::vec3(qn[v1Ind + 0], qn[v1Ind + 1], qn[v1Ind + 2]);
        const glm::vec3 v2 = glm::vec3(qn[v2Ind + 0], qn[v2Ind + 1], qn[v2Ind + 2]);
        const glm::vec3 v3 = glm::vec3(qn[v3Ind + 0], qn[v3Ind + 1], qn[v3Ind + 2]);
        const glm::mat3 DmInv = DmInvs[index];

        glm::mat3 R = glm::mat3(v0 - v3, v1 - v3, v2 - v3) * DmInv
        glm::mat3 U;
        glm::mat3 S;
        
        svd(R[0][0], R[1][0], R[2][0], R[0][1], R[1][1], R[2][1], R[0][2], R[1][2], R[2][2],
            U[0][0], U[1][0], U[2][0], U[0][1], U[1][1], U[2][1], U[0][2], U[1][2], U[2][2],
            S[0][0], S[1][0], S[2][0], S[0][1], S[1][1], S[2][1], S[0][2], S[1][2], S[2][2],
            R[0][0], R[1][0], R[2][0], R[0][1], R[1][1], R[2][1], R[0][2], R[1][2], R[2][2]);

        R = U * glm::transpose(R);

        if (glm::determinant(R) < 0)
        {
            R[2] = -R[2];
        }

        const glm::mat4x3 piTAiSi = glm::abs(V0[index]) * wi * R * glm::transpose(DmInv)
            * glm::mat4x3{ 1, 0, 0, 0, 1, 0, 0, 0, 1, -1, -1, -1 };

        atomicAdd(&(xProj[v0Ind + 0]), piTAiSi[0][0]);
        atomicAdd(&(xProj[v0Ind + 1]), piTAiSi[0][1]);
        atomicAdd(&(xProj[v0Ind + 2]), piTAiSi[0][2]);

        atomicAdd(&(xProj[v1Ind + 0]), piTAiSi[1][0]);
        atomicAdd(&(xProj[v1Ind + 1]), piTAiSi[1][1]);
        atomicAdd(&(xProj[v1Ind + 2]), piTAiSi[1][2]);

        atomicAdd(&(xProj[v2Ind + 0]), piTAiSi[2][0]);
        atomicAdd(&(xProj[v2Ind + 1]), piTAiSi[2][1]);
        atomicAdd(&(xProj[v2Ind + 2]), piTAiSi[2][2]);

        atomicAdd(&(xProj[v3Ind + 0]), piTAiSi[3][0]);
        atomicAdd(&(xProj[v3Ind + 1]), piTAiSi[3][1]);
        atomicAdd(&(xProj[v3Ind + 2]), piTAiSi[3][2]);
    }
}