Corotational 模型代码
今天看了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模型就是把变形映射回旋转前的状态后的线性模型。