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


posted @ 2014-12-18 10:36  dydx  阅读(361)  评论(0编辑  收藏  举报