Fast Algorithm for Stiffness Matrix Computing in Deformable Object Simulation

Constructing Stiffness matrix is quite a time-consuming part in FEM-based deformable object simulation. This blog provides a fast method to compute stiffness matrix for neo-hookean model, with minimal redundancy in computing. This algorithm has been implemented in IGsim. Beyond neo-hookean model, this algorithm can also be generailized to other hyperelastic models.

1. Background

1.1 Continuum Mechanics Theory

$X \in \Omega \subseteq R^3$ is the location of material in undeformed position, and the deformation function $\phi: R^3 \to R^3$ transfer material from its original position to current position. Define deformation gradient $\mathbf{F}$ as:

$$\mathbf{F} = \frac{\partial(\phi_1, \phi_2, \phi_3)}{\partial(X_1, X_2, X_3)}, ~ \mathbf{F}_{ij} = \frac{\partial \phi_i}{\partial X_j}$$

thus $\mathbf{F}$ is a tensor field of size $(3, 3)$ on $\Omega$.

Deformation energy is solely depend on the extent of deformation, thus we define deformation energy density $\Psi$ as a mapping from $\mathbf{F}$ to a non-negative real number. And the whole deformation energy is $E = \int_{\Omega} \Psi(\mathbf F)$. A typical energy density function is neo-hookean model:

$$\Psi(\mathbf{F}; \mu, \lambda) = \frac{\mu}{2}[ tr(\mathbf{F}) - \log( \det \mathbf{F}) - 3] + \frac{\lambda}{8} \log^2(\det \mathbf{F}).$$

Piola tensor $\mathbf P$ is defined as $\mathbf P = \frac{\partial \Psi(\mathbf F)}{\partial \mathbf F}$, for the convience in calculating inner force. Please refer to Finite Element Method for Soft Deformable Object Simulation for more information on deviation of Piola tensor and its property.

1.2 Simulation with Finite Elements

To enable FEM-based simulation, tetrahedronlization is performed. After tetrahedronlization and under linear assumption, the deformation gradient $\mathbf{F}$ in each tetrahedron is constant and can be derived by $\mathbf{F} = \mathbf{D}_s \mathbf{D}_{m}^{-1},$ where

$$\begin{split} &\mathbf{D}_m = [\mathbf{X}_0 - \mathbf{X}_3, \mathbf{X}_1 - \mathbf{X}_3, \mathbf{X}_2 - \mathbf{X}_3], \\ &\mathbf{D}_s ~ = [\mathbf{x}_0 - \mathbf{x}_3, \mathbf{x}_1 - \mathbf{x}_3, \mathbf{x}_2 - \mathbf{x}_3], \end{split}$$

and $\mathbf{X}_i$ and $\mathbf{x}_i$ are original and deformed positions of 4 corners of a tetrahedron respectively.

The inner force on each node of tetrahedron can be drived through

$$\mathbf H = [\mathbf{f}_0, \mathbf{f}_1, \mathbf{f}_2] = -W \mathbf{P}(\mathbf{F}) \mathbf{D}_{m}^{-T}.$$

So far it seems quite clean and efficient to compute elastic forces, which is actually the first deriative of deformation energy $\frac{\partial E}{\partial \mathbf{D}_s}$. However, Neotown’s method requires second derivative of optimization objective, namely the derivative of elastic forces. So now the problem comes to how to calculate $\frac{\partial \mathbf H}{\partial \mathbf{D}_s}$ with high efficiency.

2. Fast Stiffness Matrix Computing for Neohookean Model

2.1 Mathematical Derivation

The Piola tensor of neo-hookean model is as following:

$$\mathbf{P}(\mathbf{F}) = \mu (\mathbf{F} - \mathbf{F}^{-T}) + \lambda \log(J) \mathbf{F}^{-T},$$

where $J = \det \mathbf F$. So we can get $\delta\mathbf{P}$ with respect to $\delta\mathbf{F}$

$$\delta\mathbf{P}(\mathbf{F}) = \mu \delta \mathbf{F} + [\mu - \lambda \log(J)] \mathbf{F}^{-T} \delta \mathbf{F}^{T} \mathbf{F}^{-T} + \lambda tr(\mathbf{F}^{-1}\delta\mathbf{F})\mathbf{F}^{-T},$$

and also $\delta \mathbf H$

$$\begin{aligned} \delta \mathbf{H} =& - W \delta \mathbf{P} \mathbf{D}_{m}^{-T} \\ =& - W \{\mu \delta \mathbf{F} + [\mu - \lambda \log(J)] \mathbf{F}^{-T} \delta \mathbf{F}^{T} \mathbf{F}^{-T} + \lambda tr(\mathbf{F}^{-1}\delta\mathbf{F})\mathbf{F}^{-T}\} \mathbf{D}_{m}^{-T} \\ =& - W\mu \delta \mathbf{F} \mathbf{D}_{m}^{-T} \\ &- W [\mu - \lambda \log(J)] \mathbf{F}^{-T} \delta \mathbf{F}^{T} \mathbf{F}^{-T} \mathbf{D}_{m}^{-T} \\ &- W\lambda tr(\mathbf{F}^{-1}\delta\mathbf{F})\mathbf{F}^{-T} \mathbf{D}_{m}^{-T} \\ =& - W\mu \delta \mathbf{D}_s \mathbf{D}_{m}^{-1} \mathbf{D}_{m}^{-T} \\ &- W [\mu - \lambda \log(J)] \mathbf{F}^{-T} \mathbf{D}_{m}^{-T} \delta \mathbf{D}_s^{T} \mathbf{F}^{-T} \mathbf{D}_{m}^{-T} \\ &- W\lambda tr(\mathbf{F}^{-1}\delta \mathbf{D}_s \mathbf{D}_{m}^{-1})\mathbf{F}^{-T} \mathbf{D}_{m}^{-T} \end{aligned}$$

Denote parameters

$$a=-W\mu, ~ b=-W[\mu-\lambda \log(J)], ~ c=-W\lambda,$$

and matrices

$$\mathbf{A} = \mathbf{D}_{m}^{-1} \mathbf{D}_{m}^{-T}, ~ \mathbf{B} = \mathbf{F}^{-T} \mathbf{D}_{m}^{-T} = \mathbf{D}_s^{-T},$$

$\delta \mathbf H$ can be simplified to following expression:

$$\begin{aligned} \delta \mathbf{H} =& a \delta \mathbf{D}_s \mathbf{A} + b \mathbf{B} \delta \mathbf{D}_s^{T} \mathbf{B} + c~tr(\mathbf{F}^{-1}\delta \mathbf{D}_s \mathbf{D}_{m}^{-1}) \mathbf{B} \\ =& a \delta \mathbf{D}_s \mathbf{A} + b \mathbf{B} \delta \mathbf{D}_s^{T} \mathbf{B} + c(\mathbf{B} : \delta \mathbf{D}_s) \mathbf{B} \end{aligned}$$

At last, we can get all elements in $\frac{\partial \mathbf H}{\partial \mathbf F}$ through simple addons and multiplications:

$$\frac{\partial \mathbf{H}_{kl}}{\partial (\mathbf{D}_s)_{ij}} = a \delta^{i}_{k} \mathbf{A}_{jl} + b\mathbf{B}_{kj}\mathbf{B}_{il} + c \mathbf{B}_{ij} \mathbf{B}_{kl}$$

2.2 Code Example

For more information, please refer to elastic_neohookean.cpp from milkpku/IGsim

/* 
Inputs:
	Bm		Mat3, inverse of Dm = [X0-X3, X1-X3, X2-X3]
	Ds 		Mat3, Ds = [x0-x3, x1-x3, x2-x3]
	W 		scalar, volume of tetrahedron
Outputs:
	H[][] 	matrix of size (12, 12), elements of stiffness matrix
 */
using Mat3 = Eigen::Matrix3d

/* compute deformation gradient F */
Mat3 F = Ds * Bm[i]; 

/* accumulate stiffness matrix */
double J = F.determinant();
double a = W(i) * Mu(i);
double b = W(i) * (Mu(i) - Lam(i) * std::log(J));
double c = W(i) * Lam(i);

Mat3 A = Bm[i] * Bm[i].transpose();
Mat3 B = Ds.transpose().inverse();

/* \par H_kl / \par D_sj = a \delta{s, k} A_jl + b B_kj B_sl + c B_sj B_kl */
double H[12][12];

for (int l = 0; l < 3; l++)
for (int k = 0; k < 3; k++)
	for (int j = 0; j < 3; j++)
	for (int s = 0; s < 3; s++)
		H[3 * l + k][3 * j + s] = b * B(k, j) * B(s, l) + c * B(s, j) * B(k, l);

for (int l = 0; l < 3; l++)
	for (int s = 0; s < 3; s++)
	for (int j = 0; j < 3; j++)
		H[3 * l + s][3 * j + s] += a * A(j, l);

for (int s = 0; s < 9; s++) {
	H[9][s] = H[s][9] = -(H[s][0] + H[s][3] + H[s][6]);
	H[10][s] = H[s][10] = -(H[s][1] + H[s][4] + H[s][7]);
	H[11][s] = H[s][11] = -(H[s][2] + H[s][5] + H[s][8]);
}

H[9][9] = -(H[9][0] + H[9][3] + H[9][6]);
H[9][10] = H[10][9] = -(H[9][1] + H[9][4] + H[9][7]);
H[9][11] = H[11][9] = -(H[9][2] + H[9][5] + H[9][8]);
H[10][10] = -(H[10][1] + H[10][4] + H[10][7]);
H[10][11] = H[11][10] = -(H[10][2] + H[10][5] + H[10][8]);
H[11][11] = -(H[11][2] + H[11][5] + H[11][8]);