今天看了Corotational模型的代码。
在Vega中，获得模型内力的方法是先构造一个ForceModel对象，再调用其对应方法。
对于Corotational模型，构造的流程为:
构造CorotationalLinearFEM -> 构造CorotationalLinearFEMForceModel
获取内力时，调用顺序为:
CorotationalLinearFEMForceModel::GetInternalForce ->
CorotationalLinearFEM::ComputeForceAndStiffnessMatrix ->
CorotationalLinearFEM::ComputeForceAndStiffnessMatrixOfSubmesh
其中ComputeForceAndStiffnessMatrixOfSubmesh 计算指定的体网格元素范围（给出下标上下限）的内力。
其中包含以下几种计算的方法:1. 线性（wrap=0）2. 把旋转矩阵抵消后的线性模型(wrap = 1) 3. 准确模型（wrap = 2）
以下是源代码注释:
// computes the internal forces and (warped) stiffness matrix for the entire mesh
// vertex displacements (input) and internal forces (output) must be (pre-allocated) vectors of length 3 * numVertices
// the internal forces are returned with the sign corresponding to f_int(x) on the left side of the equation M * x'' + f_int(x) = f_ext
// i.e., the computed internal forces are *negatives* of the actual physical internal forces acting on the material
// warp:
// 0: no warping (linear FEM)
// 1: stiffness warping (corotational linear FEM with approximate stiffness matrix) [Mueller 2004]
// 2: corotational linear FEM with exact tangent stiffness matrix (see the technical report [Barbic 2012])
说明计算的内力是实际的负值。
先说明后两种的计算方法。
1.得到变形后的点P、 变形前的点M
2.F = P M^-1 得到变形梯度F
3.对F做极分解， F = RS， F正交， S对称正定
4.KElement = R * K * R^TKElement = R * K * R^T 其中K为线性模型的单元刚度矩阵
5.f = R * K * （R^T x - x0）得到每个元素的f
6.若wrap=2 计算tangent stiffness matrix，十分复杂，暂时忽略。
对于线性模型，方法很简单，就是用单元刚度矩阵乘以位移。
下一步证明，wrap=1/2的方法与 Siggraph12的course note FEM Simulation of 3D Deformable Solids给出的构造是等效的。这篇course先说明各个本构模型的能量，再推出每个本构模型力的表示法。
对于线性模型，弹性能量 E = 1/2 * d^T * K * d， K为刚度矩阵，d为每个顶点的位移向量，这一点可以通过把小变形应变张量代入线性模型的能量密度，再离散化得出。整个过程也是刚度矩阵的推导过程。
对于Corotate模型，弹性能量E = 1/2 * d‘^T * K * d’，K仍为刚度矩阵，但d‘为顶点位移抵消了旋转（乘以R^T）后的差值。这一点可以严格证明。证明的方法与刚度矩阵的推导方法类似，是根据F = RS = d(X +u ) / dX, 以及对corotate模型 ε = S - I ，推出 ε 与 du / dX 的表达式，推出 E = 1/2 * d‘^T * K * d’ 。
由此可见，Corotate模型就是把变形映射回旋转前的状态后的线性模型。