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)}MDelta 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，
加法后面惩罚图元体积的变化，以接近现实情况
再根据重力、接触力等即可模拟物体运动