[游戏物理探秘]三维隐式有限元模拟弹性物体及unity实现

量子计算9 · 发表于 2022-1-12 11:44

花了几个月终于搞定了三维隐式积分有限元，特此纪念。
三维演示视频
分别用python和unity的computeShader写成，unity场景为Implicit3D.unity。源码在
和gitee上的
接下来会涉及到非常多的公式。如果你有任何不懂的地方，欢迎问我，不过线上交流效率很低，如果你在长沙，也欢迎当面来互相交流。我一个人研究这些东西有点孤单，而且这玩意我花了几个月才弄明白，实在有太多经验想要分享了。然后推荐两个非常棒的参考资料“Dynamic Deformables: Implementation and Production Practicalities  ”和“FEM Simulation of 3D Deformable Solids: A practitioner’s guide to theory, discretization and model reduction. ”
简介

弹性物体模拟的话，简单来说就是物体不希望自己发生变形。如果外力让自己发生变形的话，那么外力就会转变成物体自身的能量，存储在物体内部。更专业的名词就是内能或者应变能。如果力没有持续作用的话，弹性物体会利用将这种能量将自己变回初始的状态。
因此模拟的主要步骤就是，先计算物体的变形程度，再由物体的变形程度计算出总能量，再对总能量求导算出每个顶点的力，然后再用这个力来更新速度，再用速度来更新位置。
我发现学习这种数值计算最有效的方法就是手算例子。比如对于一个二维的三角形，三个顶点的最开始坐标分别是

$\bold X_0 = \begin{bmatrix} 0 \\ 0\end{bmatrix} \qquad \bold X_1 = \begin{bmatrix} 1 \\ 0\end{bmatrix} \qquad \bold X_2 = \begin{bmatrix} 0 \\ 1\end{bmatrix} \tag{1}$
如果我拉伸这个三角形，让它的三个顶点变成了如下的位置

$\bold x_0 = \begin{bmatrix} 0 \\ 0\end{bmatrix} \qquad \bold x_1 = \begin{bmatrix} 2 \\ 0\end{bmatrix} \qquad \bold x_2 = \begin{bmatrix} 0 \\ 1\end{bmatrix} \tag{2}$
这时候三角形的面积改变了。根据我们现实中的经验，当我们用力拉伸一个面包时，然后不再用力，那么面包就会变回原来的形状。虽然这是段废话，但是这确实是我们先现实生活的中的经验。
现在假设这个三角形发生了平移，坐标变为如下，

$\bold x_0 = \begin{bmatrix} 0 \\ 1\end{bmatrix} \qquad \bold x_1 = \begin{bmatrix} 1 \\ 1\end{bmatrix} \qquad \bold x_2 = \begin{bmatrix} 0 \\ 2\end{bmatrix} \tag{3}$
那么这个三角形算不算变形呢？需不需要变回原来的形状呢。想想现实中的例子，当你把一个面包从手上放到桌子上，它是否会跳回你的手上呢？
第三个的例子是旋转。假如现在三个顶点的坐标如下，也就是旋转了90度

$\bold x_0 = \begin{bmatrix} 0 \\ 0\end{bmatrix} \qquad \bold x_1 = \begin{bmatrix} 0 \\ 1\end{bmatrix} \qquad \bold x_2 = \begin{bmatrix} 1 \\ 0\end{bmatrix} \tag{4}$
现在又算不算变形呢？譬如将手上的面包旋转90度，它会不会自己旋转回去呢？
根据我们在地球上观察到的物理法则，上面两个问题的答案都是“不会”。因此，我们推断，在物体平移和旋转的时候，没有变形，物体内部并没有能量产生。但是如果物体的长度，面积或体积改变，那么就会产生能量。
那么如何衡量物体的变形程度呢？我们会用到有限元中常用的变形梯度(deformation gradient)，这是一个矩阵。

$d\bold x = \bold F d\bold X \qquad \bold F = \frac{d \bold x}{d\bold X} \tag{5}$
小x是现在顶点的位置，大X是顶点初始位置，F是变形梯度。其中dX是对顶点的微分

$d\bold X = \begin{bmatrix} \bold X_1 - \bold X_0 & \bold X_2 - \bold X_0\end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix} \qquad d \bold X^{-1} = \begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix} \tag{6}$
dx也就是相似的计算，纯位移的情况下，三角形三个顶点的相对位置不发生变化，以刚才说到的平移为例，

$d\bold x = \begin{bmatrix} \bold x_1 - \bold x_0 & \bold x_2 - \bold x_0\end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix} \tag{7}$
不过有的书上对变形梯度是这么写的

$\bold F = \bold D_s \bold D_m^{-1} \qquad \bold D_s = d \bold x \qquad \bold D_m = d\bold X \tag{8}$
Ds 就是现在顶点位置就是dx，Dm就是初始顶点位置就是dX。最终变形梯度计算结果如下

$\bold F = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \tag{9}$
这就是二维情况下，纯位移的形变梯度，是一个单位矩阵。现在我们要通过一些操作，利用变形梯度计算能量，也就是把这个东西变成数字零。纯平移的变形梯度就是单位矩阵。
接下来我们要设计一个能量公式，在纯平移和旋转的时候，得出的结果为零。最容易想到的就是直接让矩阵各个数字的平方相加，然后减去维度数

$\psi = ||\bold F||_F^2 = \sum_{i=0}^2 \sum_{j=0}^2 f_{ij}^2 -2 = 1^2 + 0^2 + 0^2 + 1 ^2 - 2= 0 \tag{10}$
看起来似乎没问题，但如果三角形面积变小的话，产生的能量就是负数了，也就是外力并没有转换为内能，而是凭空消失了。说明能量公式设计起来并不容易。
线弹性

二维显式线弹性的python完整代码在
和gitee上的
那么既然之前的能量公式有问题，那么我们就换个公式呗。也就是大名鼎鼎的线弹性模型(linear elasticity)

$\psi = \mu \varepsilon :\varepsilon + \frac{\lambda}{2}tr^2(\varepsilon) \tag{11}$
其中varespilon 是应变张量(strain tensor, green strain)

$\varepsilon = (\bold F + \bold F^T)/2 - \bold I \tag{12}$
当变形梯度为单位矩阵时，应变张量就是全零矩阵。相当于没有应变。如果你对这些不熟悉，应该看看弹性力学和连续介质力学相关的书籍。在我的代码中，是这么写的
F = np.dot(Ds_new,minv)
strain = (F + F.T) * 0.5 - np.identity(2)
doubleInner = strain[0,0]*strain[0,0] + strain[1,0]*strain[1,0]
            + strain[0,1]*strain[0,1] + strain[1,1]*strain[1,1]
energy = doubleInner * mu + la * 0.5 * np.trace(strain) ** 2在bullet3库中的https://github.com/bulletphysics/bullet3/blob/master/src/BulletSoftBody/btDeformableLinearElasticityForce.h#L208，能量是这么写的
   double elasticEnergyDensity(const btSoftBody::TetraScratch& s)
   {
      double density = 0;
      btMatrix3x3 epsilon = (s.m_F + s.m_F.transpose()) * 0.5 - btMatrix3x3::getIdentity();
      btScalar trace = epsilon[0][0] + epsilon[1][1] + epsilon[2][2];
      density += m_mu * (epsilon[0].length2() + epsilon[1].length2() + epsilon[2].length2());
      density += m_lambda * trace * trace * 0.5;
      return density;
   }这种模型的优点是计算简单友好，适用于处理小变形的情况。比方说，这个三角形的x轴被拉长，顶点变成了下面这种方式

$\bold x_0 = \begin{bmatrix} 0 \\ 0\end{bmatrix} \qquad \bold x_1 = \begin{bmatrix} 2 \\ 0\end{bmatrix} \qquad \bold x_2 = \begin{bmatrix} 0 \\ 1\end{bmatrix} \tag{13}$
此时dx和变形梯度如下

$d\bold x = \begin{bmatrix} 2 & 0 \\ 0 & 1 \end{bmatrix} \qquad \bold F = \begin{bmatrix} 2 & 0 \\ 0 & 1 \end{bmatrix} \tag{14}$
那么应变张量和能量就是

$\varepsilon = \begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix} \qquad \psi = \mu + \frac{\lambda}{2} \tag{15}$
很好能量算出来了，我们就要根据能量的导数来算力了。先把stress的计算公式甩上来。先让能量对变形梯度求导得到应力，然后应力对结点求导得到真正的力。
接下来要介绍一种很特别的应力(stress)。普通的应力计算公式如下

$\lim _{\Delta S -> P} \frac{\Delta F}{\Delta A} \tag{16}$
大意就是应力等于力除以面积。当面积越小越准确。然而我们之前说了弹性物体的两种状态，也就是初始未变形的状态，以及现在已经变形的状态。我们可以很方便给初始未变形的状态上的某个区域某个大小的力，但是对于已经变形的弹性物体，我们就很难计算到底是要给哪个区域，哪些顶点施加力。
因此，我们还是直接说，我们施加的力，是基于初始的未变形的状态的，那么这种应力就是First Piola Kirchhoff stress。
这种解释来自&#34;Continuum Mechanics and Linear Elasticity - An Applied Mathematics Introduction by Ciprian D. Coman&#34;第3.8节

并且“FEM Simulation of 3D Deformable Solids: A practitioner’s guide to theory, discretization and model reduction  ” 第12页也有类似的表述

不过first piola kirchhoff 的缺点就是它是不对称的，所以不便于计算。具体来说，如果矩阵是对称的话，我们不仅能减少一半的计算量，还能玩各种操作来加快计算。first piola kirchhoff 力的另一种计算方式是能量密度对变形梯度求导

$\frac{\partial \Psi}{\partial \bold F} = \bold P(\bold F) \tag{17}$
Psi 是刚才算出的能量密度，F是变形梯度，P 是first piola krichhoff stress 。对于线弹性力来说，求导如下

$\bold P(\bold F) = \mu(\bold F + \bold F^T - 2\bold I) + \lambda tr(\bold F - \bold I)\bold I \tag{18}$
并且计算物理中关于这种矩阵求导的非常非常多。在我的代码是这么写的
piola = mu * (F + F.T - 2 * np.identity(2)) + la * (F[0,0] - 1 + F[1,1] - 1) * np.identity(2)又比如Bullet3库是这么写的，地址在https://github.com/bulletphysics/bullet3/blob/master/src/BulletSoftBody/btDeformableLinearElasticityForce.h#L390
void firstPiola(const btSoftBody::TetraScratch& s, btMatrix3x3& P)
   {
      btMatrix3x3 corotated_F = s.m_corotation.transpose() * s.m_F;

      btMatrix3x3 epsilon = (corotated_F + corotated_F.transpose()) * 0.5 - btMatrix3x3::getIdentity();
      btScalar trace = epsilon[0][0] + epsilon[1][1] + epsilon[2][2];
      P = epsilon * btScalar(2) * m_mu + btMatrix3x3::getIdentity() * m_lambda * trace;
   }再手算一遍，假设mu和lambda都是2，那么线弹性能量就是3，那么第一piola krichhoff  应力

$\bold P(\bold F) = 2\mu \bold E + \lambda tr(\bold E)\bold I \\ = (2\mu \begin{bmatrix} 1 & 0 \\ 0 & 0\end{bmatrix} + \lambda\begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix}) = \begin{bmatrix} 6 & 0 \\ 0 & 2\end{bmatrix} \tag{19}$
节点力的计算

一个三角形有三个顶点，为了动量守恒，它们的节点力相加必须为零，也就是

$\vec f_0 + \vec f_1 + \vec f_2 = 0 \tag{20}$
节点力的来源于能量对节点位置求导，我们先将能量变成能量密度乘以面积或体积，然后将偏微分分成能量密度对变形梯度求导，以及变形梯度对节点位置求导。前者就是之前的first piola kirchhoff力，后者根据变形梯度的定义更可以直接得到。

$- \vec f = \frac{\partial E}{\partial \bold x} = W\frac{\partial \Psi}{\partial \bold x}= \frac{\partial \Psi}{\partial \bold F}\frac{\partial \bold F}{\partial \bold x} = W \bold P \bold D_m^{-T} \tag{21}$
但是我们通常这么写，先算出第一个和第二个节点的力，然后根据动量守恒算出第一个节点的力

$\bold H = \begin{bmatrix} \vec f_1 & \vec f_2 \end{bmatrix} = -W\bold P(\bold F)\bold D_m^{-T} \qquad \vec f_0 = -\vec f_1 - \vec f_2 \tag{22}$
推导如下

这种节点力还要不同的写法，比如在“Robust Quasistatic Finite Elements and Flesh Simulation  ” 是这么写的

$\bold g_i = -\bold P (A_1 \bold N_1 + A_2 \bold N_2 + A_3 \bold N_3)/3 \qquad \bold g_0 = -(\bold g_1 + \bold g_2 + \bold g_3) \tag{23}$
其中AN是三角形的法向量根据面积的权重，P是之前求出的piola应力。也可以写成另一种方法

$\bold G = \bold P \bold B \qquad \bold G = (\bold g_1,\bold g_2, \bold g_3) \qquad \bold B_m = (\bold b_1,\bold b_2,\bold b_3) = -V\bold D_m^{-T} \tag{24}$
这些公式在我的代码中是这么写的
gradC1 = np.array([H[0,0],H[1,0]])
gradC2 = np.array([H[0,1],H[1,1]])
gradC0 = - gradC1 - gradC2

node_force[element_idx[ie,0],:] += gradC0
node_force[element_idx[ie,1],:] += gradC1
node_force[element_idx[ie,2],:] += gradC2其中bullet3库是这么写的，代码在https://github.com/bulletphysics/bullet3/blob/master/src/BulletSoftBody/btDeformableLinearElasticityForce.h#L263
btMatrix3x3 force_on_node123 = psb->m_tetraScratches[j].m_corotation * P * tetra.m_Dm_inverse.transpose();
btVector3 force_on_node0 = force_on_node123 * grad_N_hat_1st_col;

btSoftBody::Node* node0 = tetra.m_n[0];
btSoftBody::Node* node1 = tetra.m_n[1];
btSoftBody::Node* node2 = tetra.m_n[2];
btSoftBody::Node* node3 = tetra.m_n[3];
size_t id0 = node0->index;
size_t id1 = node1->index;
size_t id2 = node2->index;
size_t id3 = node3->index;

// elastic force
btScalar scale1 = scale * tetra.m_element_measure;
force[id0] -= scale1 * force_on_node0;
force[id1] -= scale1 * force_on_node123.getColumn(0);
force[id2] -= scale1 * force_on_node123.getColumn(1);
force[id3] -= scale1 * force_on_node123.getColumn(2);又比如positionBasedDynamics库是这么算的，地址在https://github.com/InteractiveComputerGraphics/PositionBasedDynamics/blob/master/PositionBasedDynamics/PositionBasedDynamics.cpp#L1046
void PositonBasedDynamics::computeGradCGreen(Real restVolume,
const Matrix3r &invRestMat, const Matrix3r &sigma, Vector3r *J)
{
   Matrix3r H;
   Matrix3r T;
   T = invRestMat.transpose();
   H = sigma * T * restVolume;
   J[0][0] = H(0, 0);
   J[1][0] = H(0, 1);
   J[2][0] = H(0, 2);

   J[0][1] = H(1, 0);
   J[1][1] = H(1, 1);
   J[2][1] = H(1, 2);

   J[0][2] = H(2, 0);
   J[1][2] = H(2, 1);
   J[2][2] = H(2, 2);

   J[3] = -J[0] - J[1] - J[2];
}继续手推，

$\bold H = -W\bold P(\bold F)d\bold X^{-T} = -\frac{1}{2}\begin{bmatrix} 6 & 0 \\ 0 & 2\end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix} = \begin{bmatrix} -3 & 0 \\ 0 & -1\end{bmatrix} \tag{25}$
也就是

$\vec f_0 = \begin{bmatrix} 3 \\ 1\end{bmatrix} \qquad \vec f_1 = \begin{bmatrix} -3 \\ 0\end{bmatrix} \qquad \vec f_2 = \begin{bmatrix} 0 \\ -1\end{bmatrix} \tag{26}$
最后，再整理一遍流程

算法解释如下
算法第二行是初始化，只需要对元素，也就是每个三角形(或四面体)执行一遍
算法第三行，计算每个元素的Dm或者叫dX，用于之后计算变形梯度
算法第四行，对Dm求逆得到Bm，用于之后计算变形梯度
算法第五行，计算三角形的面积(或四面体的体积)，用于之后计算节点力
算法第八行，外层主循环，每个时间步都要执行一遍
算法第九行，将节点力设置为零
算法第十行，内层循环，遍历所有三角形/(四面体)
算法第十一行，计算每个元素的Ds或者叫dx，用于计算变形梯度
算法第十二行，计算变形梯度
算法第十三行，计算first piola kirchhoff 应力
算法第十四到第十六行，计算真正的节点力
那么最简单的二维线弹性显式求解如下
import numpy as np
# 初始化三角形初始位置
node_pos = np.array([[0,0],[1,0],[0,1]],dtype = float)
# 顶点的位置梯度，对应上面算法第三行
Ds = np.array([[node_pos[1,0] - node_pos[0,0],node_pos[2,0] -
node_pos[0,0]],
[node_pos[1,1] - node_pos[0,1],node_pos[2,1] -
node_pos[0,1]]])
# 求逆，用于准备计算形变梯度，对应上面算法第四行
minv = np.linalg.inv(Ds)
# 假设某一时刻，三角形变化到了这样的位置
node_pos = np.array([[0,0],[2,0],[0,1]],dtype = float)
time = 0
timeFinal = 100
areat = np.zeros((timeFinal))
while(time < timeFinal):
   time += 1
   # 形变梯度中的分子，对应上面算法第11行
   Ds_new = np.array([[node_pos[1,0] - node_pos[0,0],node_pos[2,0]
   - node_pos[0,0]],
   [node_pos[1,1] - node_pos[0,1],node_pos[2,1] -
   node_pos[0,1]]])
   # 形变梯度，对应上面算法第12行
   F = np.dot(Ds_new,minv)
   # 应力，也就是varepsilon
   strain = (F + F.T) * 0.5 - np.identity(2)
   # lame常数
   mu = 2
   # lame常数
   la = 2
   #
   doubleInner = strain[0,0]*strain[0,0] + strain[1,0]*strain[1,0]
   doubleInner += strain[0,1]*strain[0,1] + strain[1,1]*strain[1,1]
   # 线弹性的能量计算公式
   energy = doubleInner * mu + la * 0.5 * np.trace(strain) ** 2
   #first piola kirchhoff stress，对应上面算法第13行
   piola = mu * (F + F.T - 2 * np.identity(2)) + la * (F[0,0] - 1 + F[1,1] - 1) * np.identity(2)
   # 三角形面积
   area = 0.5
   # 计算hessian矩阵对应上面算法第14行
   H = - area * np.dot(piola,minv.transpose())
   # 计算节点力，对应上面算法第15和第16行
   gradC1 = np.array([H[0,0],H[1,0]])
   gradC2 = np.array([H[0,1],H[1,1]])
   gradC0 = - gradC1 - gradC2
   invMass = 1
   dt = 0.1
   # 判断是否收敛
   sumGradC = invMass * (gradC0[0]**2 + gradC0[1]**2)
   sumGradC += invMass * (gradC1[0]**2 + gradC1[1]**2)
   sumGradC += invMass * (gradC2[0]**2 + gradC2[1]**2)
   if sumGradC < 1e-10:
      break
   # 校正位置，方法来源于PositionBasedDynamics
   node_pos[0,:] += dt * energy / sumGradC * invMass * gradC0
   node_pos[1,:] += dt * energy / sumGradC * invMass * gradC1
   node_pos[2,:] += dt * energy / sumGradC * invMass * gradC2
   areat[time - 1] = 0.5 * (node_pos[0,0] * (node_pos[1,1] - node_pos[2,1])
                  + node_pos[1,0] * (node_pos[2,1] - node_pos[0,1])
                  + node_pos[2,0] * (node_pos[0,1] - node_pos[1,1]))
   但是不要高兴得太早了，如果使用上面的代码，你会发现甚至算不出正确的结果，无法收敛。因为线弹性忽略了一个非常重要的东西，旋转。因旋转并未造成面积改变。举个例子，仍然是之前那个三角形，不过现在逆时针旋转90度

$\bold x_0 = \begin{bmatrix} 0 \\ 0\end{bmatrix} \qquad \bold x_1 = \begin{bmatrix} 0 \\ 1\end{bmatrix} \qquad \bold x_2 = \begin{bmatrix} -1 \\ 0\end{bmatrix} \tag{27}$
此时dx和变形梯度如下

$d\bold x = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix} \qquad \bold F = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix} \tag{28}$
那么应变张量和能量就是

$\varepsilon = \begin{bmatrix} -1 & -1 \\ 1 & -1 \end{bmatrix} \qquad \psi = 4\mu + 4 \frac{\lambda}{2} \tag{29}$
因为mu和lambda不为零，那么能量也必不为零，但是，你所做的，仅仅是旋转了物体而已，却也凭空产生能量。所以用线弹性模型很少在各种物理模拟库中见到。
STVK模型

二维显式stvk的代码在
和gitee上的
线弹性模型失败的原因来源于它认为旋转也是变形，所以会产生能量。因此我们要重新设计一个能量公式，让旋转不产生能量。最简单的非线性模型就是St. Venant-Kirchhoff  能量模型，它长下面这样

$\Psi(\bold F) = \mu \bold E :\bold E + \frac{\lambda}{2}tr^2(\bold E) \tag{30}$
同时把测量应变的矩阵换成了green lagrange strain，如下

$\bold E = \frac{1}{2}(\bold F^T \bold F - \bold I) \tag{31}$
这种应变测量方法也不会把旋转当作变形。这种模型的缺点之后再说。现在只要知道使用这种模型，那么物体旋转时就不会产生能量，甩线弹性几条街就行了。
此时stvk模型的first piola kirchhoff 应力计算如下

$\bold P(\bold F) = \bold F(2\mu \bold E + \lambda tr(\bold E)\bold I) \tag{32}$
仍然使用之前的代码，然后仅仅需要改变三行代码，就能将线弹性模型换成stvk模型
# Green Strain，也就是E
E = (np.dot(F,F.T)- np.identity(2)) * 0.5
# Stvk的能量计算公式
energy = doubleInner * mu + la / 2 * (E[0,0] + E[1,1])**2
#first piola kirchhoff stress
piola = np.dot(F, 2 * mu * E + la * (E[0,0] + E[1,1]) * np.identity(2))最后计算的各顶点位置和形变梯度如下

$\bold x_0 =\begin{bmatrix} 0.427 \\ -0.072\end{bmatrix} \qquad \bold x_1 =\begin{bmatrix} 1.395 \\ 0.177\end{bmatrix} \qquad \bold x_2 =\begin{bmatrix} 0.177 \\ 0.895\end{bmatrix} \qquad \bold F = \begin{bmatrix} 0.968 & -0.250 \\ 0.250 & 0.968\end{bmatrix} \tag{33}$
接下来看看有哪些开源库也写了stvk模型吧。
比如Quasi-Newton-Methods-For-Real-Time-Simulation库，也是是太极图形课主讲刘天添老师的仓库。在其中的constraint_tet.cpp中，地址为https://github.com/ltt1598/Quasi-Newton-Methods-for-Real-time-Simulation-of-Hyperelastic-Materials/blob/db9a9e3ebd6eec08a5682c0bc5f7c1cc3f6f024d/GenPD/GenPD/source/constraint_tet.cpp#L266
有这么一段代码
   case MATERIAL_TYPE_StVK:
   {
      EigenMatrix3 I = EigenMatrix3::Identity();
      EigenMatrix3 E = 0.5*(F.transpose()*F - I);
      P = F * (2 * m_mu*E + m_lambda*E.trace() * I);
      e_this = m_mu*E.squaredNorm() + 0.5*m_lambda*std::pow(E.trace(), 2);
      ScalarType J = F.determinant();
      if (J < 1)
      {
         P += -m_kappa / 24 * std::pow((1 - J) / 6, 2) * J * F.inverse().transpose();
         e_this += m_kappa / 12 * std::pow((1 - J) / 6, 3);
      }
   }著名的positionBasedDynamics库https://github.com/InteractiveComputerGraphics/PositionBasedDynamics/blob/master/PositionBasedDynamics/PositionBasedDynamics.cpp#L1110中
const Real trace = epsilon(0, 0) + epsilon(1, 1) + epsilon(2, 2);
   const Real ltrace = lambda*trace;
   sigma = epsilon * 2.0*mu;
   sigma(0, 0) += ltrace;
   sigma(1, 1) += ltrace;
   sigma(2, 2) += ltrace;
   sigma = F * sigma;

   Real psi = 0.0;
   for (unsigned char j = 0; j < 3; j++)
      for (unsigned char k = 0; k < 3; k++)
         psi += epsilon(j, k) * epsilon(j, k);
   psi = mu*psi + static_cast<Real>(0.5)*lambda * trace*trace;
   energy = restVolume * psi;另一份不错的代码来自耶鲁大学图形组助理教授Theodore Kim的开源库Cubica。他的个人主页http://www.tkim.graphics/，Cubica主页http://www.tkim.graphics/cubica/。嗯，这些教授有很多代码虽然开源了，但并不放到github上，必须通过教授主页才能找到。所以为了找到这些代码真的非常花时间。
在stvk.cpp中，有这么一段代码
MATRIX3 STVK::secondPiolaKirchhoff(MATRIX3& F, bool diagonal)
{
   MATRIX3 C = F.transpose() * F;

   // Green-Lagrange strain tensor
   MATRIX3 E = 0.5 * (C - MATRIX3::I());

// 2nd Piola-Kirchoff stress tensor
   MATRIX3 S = (_lambda * trace(E) * MATRIX3::I()) + (2 * _mu * E);

   return S;
}另外，为了提升代码效率，Cubica库非常丧心病狂地用maple和matlab生产了将矩阵运算全部展开的c++代码。【在大佬那里也许是常规操作，但我真的是第一次见到】，原因是从节点位置矩阵，到计算出节点力的过程中，要经过相当多的矩阵乘法，转置，相加相减。其中有些矩阵是对称的，有的矩阵有很多零，导致很多乘法和加法是重复的，还有很多乘零和加零的操作，尤其是使用隐式方法的时候。所以自动生成的代码就是将这些重复的没必要的加法和乘法去掉了。部分代码如下
const double a = nodes[0][0]; // 节点位置
const double b = nodes[0][1];
const double c = nodes[0][2];
const double d = nodes[1][0];
...
const double m = matInv(0,0); // 计算变形梯度用到的Dm或者dX
const double n = matInv(0,1);
const double o = matInv(0,2);
...
const double t1 = d-a;
const double t3 = g-a;
const double t5 = j-a;
const double t7 = t1*m+t3*p+t5*s;
...
const double t55 = t10*o+t12*r+t14*u;
const double t56 = t55*t55;
const double t61 = t19*o+t21*r+t23*u;
const double t62 = t61*t61;
const double t65 = lam*(t8/2.0+t17/2.0+t26/2.0-3.0/2.0+t32/2.0+t38/2.0+t44/2.0+t50/2.0+t56/2.0+t62/2.0);
...
forces[0][0] = -t65*(t7*t66+t31*t68+t49*t70)-mu*(2.0*t75*t66+t81*(t66*t31+t7*t68)+t90*(t66*t49+t7*t70)+2.0*t97*t68+t103*(t68*t49+t31*t70)+2.0*t110*t70);
forces[0][1] = -t65*(t16*t66+t37*t68+t55*t70)-mu*(2.0*t121*t66+t81*(t66*t37+t16*t68)+t90*(t66*t55+t16*t70)+2.0*t134*t68+t103*(t68*t55+t37*t70)+2.0*t142*t70);stvk模型虽然是不错的模型，但是到现在为止我们还用的是显式积分方法，顶点数多，力大，时间步长大的话，就很容易算过头，然后穿模还拉不回来。所以下一篇将使用隐式积分方法，并在unity上实现cpu二维隐式和gpu三维隐式方法。
下一篇：

fwalker · 发表于 2022-1-12 11:52

好棒，支持一下

LiteralliJeff · 发表于 2022-1-12 11:55

很专业当时偷懒了推了一个简单拟合函数实现但简单只需要2天实现3d 胶体 https://zhuanlan.zhihu.com/p/361126215

mypro334 · 发表于 2022-1-12 11:56

不错的工作

Doris232 · 发表于 2022-1-12 12:02

矩阵相乘运算展开到元素的加减乘除能提高效率？

BlaXuan · 发表于 2022-1-12 12:11

从节点位置矩阵，到计算出节点力的过程中，要经过相当多的矩阵乘法，转置，相加相减。其中有些矩阵是对称的，有的矩阵有很多零，导致很多乘法和加法是重复的，还有很多乘零和加零的操作，尤其是使用隐式方法的时候。所以自动生成的代码就是将这些重复的没必要的加法和乘法去掉了。

Zephus · 发表于 2022-1-12 12:11

这是非线性连续介质力学里面的东西，游戏现在要做柔性体了么，你们游戏方向有要招聘力学方向博士么

ainatipen · 发表于 2022-1-12 12:17

也许你可以试试泽森科工这家公司，他们可能需要力学方向的，网址在zenus

XGundam05 · 发表于 2022-1-12 12:18

需要的

TheLudGamer · 发表于 2022-1-12 12:24

这家公司据我了解是需要程序员多一些，核心力学算法都是现成的，好像不是那么需要力学方向的

		自动登录	找回密码
密码			立即注册

[笔记] [游戏物理探秘]三维隐式有限元模拟弹性物体及unity实现

本帖子中包含更多资源