FEM 把物体分为多个基本图元，如四面体。
我们可以用四面体的顶点集合来表示这个物体。
1. FEM物理引擎工作流
Loop:
通过受力更新顶点速度
碰撞检测<先导顶点位置更新>
顶点位置更新
2. 顶点位置更新
为了达到更好的效果使用后向欧拉法，
x表示位置，t表示时间，可以导出两个公式：
下时刻位置=这时刻位置+这时刻到下时刻的速度*时间差；
这时刻到下时刻的速度=这时刻速度+物体下时刻位置对应受力/质量*时间差
（比较准确的速度表示方法是向后差分）
对应数学表达
\[x_{t+1} = x_t + v_{t+1}\Delta t \v_{t+1} = v_t + \frac{f(x_{t+1})}{M}\Delta t
\]
消去x_{t+1}
\[v_{t+1} = v_t + \frac{f(x_t + v_{t+1}\Delta t)}{M}\Delta t
\]
我们在x_t点泰勒展开一阶来近似
\[v_{t+1} = v_t + (f(x_t) + \frac{\partial f}{\partial x}(x_t) + v_{t+1}\Delta t) \frac{\Delta t}M
\]
化简得
\[(1 - \frac{\frac{\partial f}{\partial x}(x_t) \Delta t^2}M)v_{t+1} = v_t + \frac{f(x_t)}M\Delta t
\]
这里只要求一个自变量v_{t+1}，使用迭代法求解
3. 处理物体图元间的内力
一个四面体有四个顶点，可以看做一个 4 * 3 的矩阵 X
现在我们可以通过一个 3 * 3 的变形矩阵 R 以及平移向量 b 映射，得到变形后的四面体
\[X_{old} R + b = X_{new}
\]
变形的内力指向形变减少的方向，可以通过偏导计算变形位移，并用胡克定律计算
但是这里并不直接计算内力，而是用一个势能场来计算势能 E，以更好地泛化问题。
\[E = \psi(\frac{\partial X_{new}}{\partial X_{old}}) = \psi(GR^T)
\]
上面的 G 表示一个形状与 X 一样的全为 1 的矩阵
这里可以获得整个物体的内力势能，即所有四面体的势能之和:
\[U = \sum_i E_i
\]
由于内力会作用在整个物体势能减小的方向，
因此，物体全局势能 U 对每一个顶点的负梯度就等效于该点的受力
\[f_i = - \frac{\partial U}{\partial x_i}
\]
有限元模型就是在探讨如何选择势能场的问题，常用的比如:
- 超弹性材料 Neo-Hookean 应用于橡胶、生物组织:
\[F = GR^T \J = \frac{V_{new}}{V_{old}} = |R| \\psi(F) = \frac{\mu}{2}\sum_i[(F^TF)_{ii}-1] - \mu log(J) + \frac{\lambda}2 log^2(J)
\]
- 小形变 Corotated 应用于金属、类刚体等:
\[\psi(F) = \mu \sum_i[(\sigma_i-1)^2] + \frac{\lambda}2(J-1)^2
\]
这里需要一个自动微分求解库来完成操作，可以使用Torch，
加法后面惩罚图元体积的变化，以接近现实情况
再根据重力、接触力等即可模拟物体运动