弹性有限元与超弹性物体模拟方法
在引入有限元之前，先简单介绍相关的物理理论。在后面的部分，使用粗体符号如 $\boldsymbol{x}$ 表示向量或矩阵(张量), 未加粗的符号为标量, 如 $x$。
1、形变（Deformation）
当弹性材料发生形变，其上点$x$会移动到新的位置，我们形变映射(Deformation map) $\boldsymbol \phi$ 表达此关系,它是一个向量到向量的映射。有：
$$\boldsymbol x_{deformed} = \boldsymbol \phi( \boldsymbol x_{rest})$$
为了更好地描述形变，通常使用形变梯度(deformation gradient)来表示这个形变。形变梯度定义为：
$$\boldsymbol F = \frac{\partial \boldsymbol \phi( \boldsymbol x_{rest})}{\partial \boldsymbol x_{rest}} = \frac{\boldsymbol x_{deformed}}{\boldsymbol x_{rest}} \tag{1}$$
这里 $\boldsymbol F$是一个n阶张量，2D问题就是一个二阶张量(2 x 2 矩阵)，3D就是一个三阶张量(3 x 3矩阵)。
$\boldsymbol F$的行列式（通常用$\boldsymbol J$表示）也是非常有用的，因为它的表征无限小体积的变化。它通常表示为:
$$\boldsymbol J = det(\boldsymbol F) \tag{2}$$
$\boldsymbol J$表示在无限小体积（面积）相对原始体积（面积）的比率。例如，在刚性运动（旋转和平移）中这很好理解， $\boldsymbol F$是旋转矩阵且$\boldsymbol J = 1$。请注意，单位矩阵也是旋转矩阵。 $\boldsymbol J > 1$表示体积(面积)增加，$\boldsymbol J < 1$表示体积(面积)减小。
$\boldsymbol J = 0$意味着体积变成0。在真实世界中，这将不会发生。然而，数值上获得这样的$\boldsymbol F$是可能的。在3D中，这表明物质被压缩使其成为一个面或一条线或是一个无体积的点。 $\boldsymbol J < 0$意味着物体被反转。考虑2D中的一个三角形， $\boldsymbol J < 0$意味着某个顶点穿过了它相对的边，然后面积变成了负值。
2、弹性（Elasticity）
2.1 应变势能（Strain energy）与应力（stress）
对于超弹性材料,其应力与应变关系由一个应变能量密度函数定义:
$$\Psi = \Psi(\boldsymbol{F}) \tag{3}$$
直观地理解：$\Psi $ 是惩罚形变的隐函数。这个势函数与物体抵抗形变所产生的力有关，这个力叫应力(stress)，它是材料内部的一种弹性力。而之前 (1) 定义的形变梯度$\boldsymbol F$也叫做应变（strain）。
知道了势能密度之后，对于一个物体，其总的弹性势能就是：
$$E = \int_{\Omega} \boldsymbol{\Psi(F)} d \boldsymbol{X} \tag{4}$$
因此，物体在任意一点的受力可以通过势能的负梯度来计算：
$$ \boldsymbol {f(x)} = - \frac{\partial E}{\partial \boldsymbol{x}} \tag{5}$$
应力代表无限小体积（面积）的材料成分对其附近施加的内力。为了表示应力，定义了不同的度量来表示：
First Piola-Kirchhof(PK1)应力 $\boldsymbol {P(F)} = \frac{\partial \boldsymbol{\Psi}( \boldsymbol{F})}{\partial \boldsymbol{F}} $
Kirchhoff 应力： $\boldsymbol \tau$
柯西（cauchy）应力: $\boldsymbol \sigma$
三种应力都是一个张量，维度与应变 $\boldsymbol F$ 一样。三种应变可以相互转化：
$$\boldsymbol \tau = J\boldsymbol \sigma = \boldsymbol{PF}^{T}, \boldsymbol P = J \boldsymbol {\sigma F}^{-T} \tag{6}$$
图形学模拟里一般常用PK1应力与柯西应力。
2.2 弹性模量（Elastic moduli）
物理模拟中，常常使用各向同性模型，指的沿物体任意一个方向施加形变，沿该方向对抗该形变的力都是一样的。为了描述各向同性物体的弹性，引入弹性模量来度量弹性，有以下几种：
Young’s modulus $E = \frac{\sigma}{\varepsilon}$
Bulk modulus $K = - V \frac{dP}{dV}$
Poisson’s ratio $\nu \in [0 , 0.5]$
Lam$\acute{e}$’s first parameter $\mu$
Lam$\acute{e}$’s second parameter $\lambda$
他们之间也可以相互转化,
$$K = \frac{E}{3(1 - 2\nu)}, \lambda = \frac{E\nu}{(1 + \nu)(1 - 2\nu)}, \mu = \frac{E}{2(1 + \nu)} \tag{7}$$
图形学模拟里一般常常设定 $E$ 与 $\nu$ 的值, 其他值通过这两个算出来。
2.3 超弹性材料模型
图形学里常常使用一个确定的函数来描述弹性势能与应变间的关系，进而将应变与应力联系起来。一般常用的线性弹性模型有：
（1）Neo-Hookean模型
$$\boldsymbol{\Psi}( \boldsymbol{F}) = \frac{\mu}{2}\sum_{i}[(\boldsymbol{FF^{T}})_{ii}-1]-\mu log(J) + \frac{\lambda}{2}log^{2}(J) \tag{8}$$
$$\boldsymbol{P}( \boldsymbol{F}) = \frac{\partial \Psi}{\partial \boldsymbol F} = \mu(\boldsymbol {F-F^{T}}) + \lambda log(J) \boldsymbol{F^{-T}} \tag{9}$$
（2） (Fixed) Corotated模型
$$\boldsymbol{\Psi}( \boldsymbol{F}) = \mu \sum_{i}(\boldsymbol{\sigma_{i}-1})^{2}- \frac{\lambda}{2}(J - 1)^{2} \tag{10}$$
$$\boldsymbol{P}( \boldsymbol{F}) = \frac{\partial \Psi}{\partial \boldsymbol F} = 2\mu(\boldsymbol {F-R}) + \lambda (J - 1)J\boldsymbol{F^{-T}} \tag{11}$$
公式（10）里的 $\sigma_{i}$ 是应变 $\boldsymbol F$的奇异值。公式（11)里的 $R$ 是 $F$ 进行极分解 $\boldsymbol {F = RS}$ 后得到的。
3、线性有限元方法（Linear finite element method)
有限元方法是一种 Galerkin 离散化方法，使用连续偏微分方程的弱形式来构建离散方程。直观地讲，就是讲一个连续的物体或表面划分为一个个微元，也就是element,在一个个微元上根据物理性质来构建方程逼近连续的求解结果。
在物理模拟中，常常使用三角形（2D）或四面体（3D）来作为相应的有限元element。
引入有限元之后，公式（4）可以重写成：
$$E = \sum_{e}E^{e}(x) = \sum_{e} \int_{\Omega_{e}} \boldsymbol{\Psi(F)} d \boldsymbol{X} = \sum_{e} V_{e} \boldsymbol{\Psi(F_{e})} \tag{12}$$
相应地，每个微元上的点受力可以写成：
$$ \boldsymbol {f(x)} = - \frac{\partial E(x)}{\partial \boldsymbol{x}}
= - \sum_{e} \frac{\partial E^{e}(x)}{\partial \boldsymbol{x}}
= - \sum_{e} V_{e}\frac{\partial \boldsymbol\Psi(\boldsymbol F_{e})}{\partial {\boldsymbol F_{e}}} \frac{\partial \boldsymbol F_{e}}{\partial {\boldsymbol x}}
= - \sum_{e} V_{e} \boldsymbol P(\boldsymbol F_{e}) \frac{\partial \boldsymbol F_{e}}{\partial {\boldsymbol x}} \tag{13}$$
3.1 四面体有限元与三角形有限元
线性有限元（用于弹性）假设形变映射是仿射的，因此变形梯度 $\boldsymbol F$ 在单个element内是常量。
对于四面体微元，有
表示成
$$\boldsymbol D_{S}= \boldsymbol F \boldsymbol D_{m} \tag{14}$$
其中 $\boldsymbol D_{s} = \left[
\begin{matrix}
\boldsymbol x_{1} - \boldsymbol x_{4} & \boldsymbol x_{2} - \boldsymbol x_{4} & \boldsymbol x_{3} - \boldsymbol x_{4}
\end{matrix}
\right]$ 被称为deformed shape矩阵, $\boldsymbol D_{m} = \left[
\begin{matrix}
\boldsymbol X_{1} - \boldsymbol X_{4} & \boldsymbol X_{2} - \boldsymbol X_{4} & \boldsymbol X_{3} - \boldsymbol X_{4}
\end{matrix}
\right]$ 被称为reference shape矩阵。它是一个常数矩阵，只与初始状态有关。
四面体的体积可以计算为：
$$V_{e} = \frac{1}{6}|det(\boldsymbol D_{m})| \tag{15}$$
由$\boldsymbol D_{m}$ 可逆， 可以求出形变梯度 $\boldsymbol F$ :
$$\boldsymbol F = \boldsymbol D_{s} \boldsymbol D_{m}^{-1} \tag{16}$$
为了计算四面体每个顶点的受力，需要计算 $\frac{\partial \boldsymbol F_{e}}{\partial {\boldsymbol x}}$， 这里直接给出受力的计算结果。具体推导看这里。
$$\boldsymbol H = [\boldsymbol f_{i}, \boldsymbol f_{2}, \boldsymbol f_{3}] = -V_{e} \boldsymbol P(\boldsymbol F) \boldsymbol D_{m}^{-T} \tag{17}$$
这里的 $\boldsymbol P(\boldsymbol F)$ 可以直接使用公式（9）或（11）的结果。$\boldsymbol H$分别是1，2，3三个顶点的受力，另一个顶点4的受力为 $\boldsymbol f_{4} = - \boldsymbol f_{1} - \boldsymbol f_{2} - \boldsymbol f_{3}$ 。
基于之上的讨论，一个可以用于计算四面体有限元弹性力的伪代码流程如下：
对于三角形网格，其推导是类似的，只是少了一个维度。
3.2 显式时间积分
基于上一节的讨论，已经得到了计算每个每个顶点上的受力。使用半隐式欧拉迭代方法更新点的位置和速度。
$$\boldsymbol v_{t+1,i} = \boldsymbol v_{t,i} + \Delta t \frac{\boldsymbol f_{t,i} }{m_{i}} \tag{18}$$
$$\boldsymbol x_{t+1,i} = \boldsymbol x_{t,i} + \Delta t \boldsymbol v_{t+1,i}\tag{19}$$
3.3 隐式时间积分
隐式时间积分的更新迭代如下：
$$[\boldsymbol I - \Delta t^{2}\boldsymbol M ^{-1}\frac{\partial \boldsymbol f}{\partial \boldsymbol x}(\boldsymbol x_{t})]\boldsymbol v_{t+1} = \boldsymbol v_{t} + \Delta t \boldsymbol M^{-1}\boldsymbol f(\boldsymbol x_{t}) \tag{20}$$
这还需要计算力的导数 $\frac{\partial \boldsymbol f}{\partial \boldsymbol x}$, $-\frac{\partial \boldsymbol f}{\partial \boldsymbol x}$被称为stiffness matrix，在实际迭代求解中我们并不用真的构建这个矩阵，我们只需要知道它和某个向量 $ \boldsymbol w$ 的乘积结果即可进行迭代。下面的算法是力的微分$\delta \boldsymbol f = \frac{\partial \boldsymbol f}{\partial \boldsymbol x} \delta \boldsymbol x$的计算方式：
具体推导。
4、引用
[1] E. Sifakis and J. Barbic (2012). “FEM simulation of 3D deformable solids: a practitioner’s guide to theory, discretization and model reduction”. In: Acm siggraph 2012 courses, pp. 1–50.
[2] C. Jiang et al. (2016). “The material point method for simulating continuum materials”. In: ACM SIGGRAPH 2016 Courses, pp. 1–52.
