[游戏物理探秘]三维隐式有限元模拟弹性物体及unity实现
花了几个月终于搞定了三维隐式积分有限元，特此纪念。三维演示视频
分别用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. ”
简介
弹性物体模拟的话，简单来说就是物体不希望自己发生变形。如果外力让自己发生变形的话，那么外力就会转变成物体自身的能量，存储在物体内部。更专业的名词就是内能或者应变能。如果力没有持续作用的话，弹性物体会利用将这种能量将自己变回初始的状态。
因此模拟的主要步骤就是，先计算物体的变形程度，再由物体的变形程度计算出总能量，再对总能量求导算出每个顶点的力，然后再用这个力来更新速度，再用速度来更新位置。
我发现学习这种数值计算最有效的方法就是手算例子。比如对于一个二维的三角形，三个顶点的最开始坐标分别是
https://www.zhihu.com/equation?tex=%5Cbold+X_0+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+X_1+%3D+%5Cbegin%7Bbmatrix%7D+1+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+X_2+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+1%5Cend%7Bbmatrix%7D+%5Ctag%7B1%7D
如果我拉伸这个三角形，让它的三个顶点变成了如下的位置
https://www.zhihu.com/equation?tex=%5Cbold+x_0+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_1+%3D+%5Cbegin%7Bbmatrix%7D+2+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_2+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+1%5Cend%7Bbmatrix%7D+%5Ctag%7B2%7D
这时候三角形的面积改变了。根据我们现实中的经验，当我们用力拉伸一个面包时，然后不再用力，那么面包就会变回原来的形状。虽然这是段废话，但是这确实是我们先现实生活的中的经验。
现在假设这个三角形发生了平移，坐标变为如下，
https://www.zhihu.com/equation?tex=%5Cbold+x_0+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+1%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_1+%3D+%5Cbegin%7Bbmatrix%7D+1+%5C%5C+1%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_2+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+2%5Cend%7Bbmatrix%7D+%5Ctag%7B3%7D
那么这个三角形算不算变形呢？需不需要变回原来的形状呢。想想现实中的例子，当你把一个面包从手上放到桌子上，它是否会跳回你的手上呢？
第三个的例子是旋转。假如现在三个顶点的坐标如下，也就是旋转了90度
https://www.zhihu.com/equation?tex=%5Cbold+x_0+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_1+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+1%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_2+%3D+%5Cbegin%7Bbmatrix%7D+1+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Ctag%7B4%7D
现在又算不算变形呢？譬如将手上的面包旋转90度，它会不会自己旋转回去呢？
根据我们在地球上观察到的物理法则，上面两个问题的答案都是“不会”。因此，我们推断，在物体平移和旋转的时候，没有变形，物体内部并没有能量产生。但是如果物体的长度，面积或体积改变，那么就会产生能量。
那么如何衡量物体的变形程度呢？我们会用到有限元中常用的变形梯度(deformation gradient)，这是一个矩阵。
https://www.zhihu.com/equation?tex=d%5Cbold+x+%3D+%5Cbold+F+d%5Cbold+X+%5Cqquad+%5Cbold+F+%3D+%5Cfrac%7Bd+%5Cbold+x%7D%7Bd%5Cbold+X%7D+%5Ctag%7B5%7D
小x是现在顶点的位置，大X是顶点初始位置，F是变形梯度。其中dX是对顶点的微分
https://www.zhihu.com/equation?tex=d%5Cbold+X+%3D+%5Cbegin%7Bbmatrix%7D+%5Cbold+X_1+-+%5Cbold+X_0+%26+%5Cbold+X_2+-+%5Cbold+X_0%5Cend%7Bbmatrix%7D+%3D+%5Cbegin%7Bbmatrix%7D+1+%26+0++%5C%5C+0+%26+1%5Cend%7Bbmatrix%7D+%5Cqquad+d+%5Cbold+X%5E%7B-1%7D+%3D+%5Cbegin%7Bbmatrix%7D+1+%26+0++%5C%5C+0+%26+1%5Cend%7Bbmatrix%7D+%5Ctag%7B6%7D
dx也就是相似的计算，纯位移的情况下，三角形三个顶点的相对位置不发生变化，以刚才说到的平移为例，
https://www.zhihu.com/equation?tex=d%5Cbold+x+%3D+%5Cbegin%7Bbmatrix%7D+%5Cbold+x_1+-+%5Cbold+x_0++%26+%5Cbold+x_2+-+%5Cbold+x_0%5Cend%7Bbmatrix%7D+%3D+%5Cbegin%7Bbmatrix%7D+1+%26+0++%5C%5C+0+%26+1%5Cend%7Bbmatrix%7D+%5Ctag%7B7%7D
不过有的书上对变形梯度是这么写的
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbold+D_s+%5Cbold+D_m%5E%7B-1%7D+%5Cqquad+%5Cbold+D_s+%3D+d+%5Cbold+x+%5Cqquad+%5Cbold+D_m+%3D+d%5Cbold+X+%5Ctag%7B8%7D
Ds 就是现在顶点位置就是dx，Dm就是初始顶点位置就是dX。最终变形梯度计算结果如下
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D+1+%26+0+%5C%5C+0+%26+1+%5Cend%7Bbmatrix%7D+%5Ctag%7B9%7D
这就是二维情况下，纯位移的形变梯度，是一个单位矩阵。现在我们要通过一些操作，利用变形梯度计算能量，也就是把这个东西变成数字零。纯平移的变形梯度就是单位矩阵。
接下来我们要设计一个能量公式，在纯平移和旋转的时候，得出的结果为零。最容易想到的就是直接让矩阵各个数字的平方相加，然后减去维度数
https://www.zhihu.com/equation?tex=%5Cpsi+%3D+%7C%7C%5Cbold+F%7C%7C_F%5E2+%3D+%5Csum_%7Bi%3D0%7D%5E2+%5Csum_%7Bj%3D0%7D%5E2+f_%7Bij%7D%5E2+-2+%3D+1%5E2+%2B+0%5E2+%2B+0%5E2+%2B+1+%5E2+-+2%3D+0+%5Ctag%7B10%7D
看起来似乎没问题，但如果三角形面积变小的话，产生的能量就是负数了，也就是外力并没有转换为内能，而是凭空消失了。说明能量公式设计起来并不容易。
线弹性
二维显式线弹性的python完整代码在
和gitee上的
那么既然之前的能量公式有问题，那么我们就换个公式呗。也就是大名鼎鼎的线弹性模型(linear elasticity)
https://www.zhihu.com/equation?tex=%5Cpsi+%3D+%5Cmu+%5Cvarepsilon+%3A%5Cvarepsilon+%2B+%5Cfrac%7B%5Clambda%7D%7B2%7Dtr%5E2%28%5Cvarepsilon%29+%5Ctag%7B11%7D
其中varespilon 是 应变张量(strain tensor, green strain)
https://www.zhihu.com/equation?tex=%5Cvarepsilon+%3D+%28%5Cbold+F+%2B+%5Cbold+F%5ET%29%2F2+-+%5Cbold+I++%5Ctag%7B12%7D
当变形梯度为单位矩阵时，应变张量就是全零矩阵。相当于没有应变。如果你对这些不熟悉，应该看看弹性力学和连续介质力学相关的书籍。在我的代码中，是这么写的
F = np.dot(Ds_new,minv)
strain = (F + F.T) * 0.5 - np.identity(2)
doubleInner = strain*strain + strain*strain
+ strain*strain + strain*strain
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 + epsilon + epsilon;
density += m_mu * (epsilon.length2() + epsilon.length2() + epsilon.length2());
density += m_lambda * trace * trace * 0.5;
return density;
}这种模型的优点是计算简单友好，适用于处理小变形的情况。比方说，这个三角形的x轴被拉长，顶点变成了下面这种方式
https://www.zhihu.com/equation?tex=%5Cbold+x_0+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_1+%3D+%5Cbegin%7Bbmatrix%7D+2+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_2+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+1%5Cend%7Bbmatrix%7D+%5Ctag%7B13%7D
此时dx和变形梯度如下
https://www.zhihu.com/equation?tex=d%5Cbold+x+%3D+%5Cbegin%7Bbmatrix%7D+2+%26+0+%5C%5C+0+%26+1+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+F+%3D++%5Cbegin%7Bbmatrix%7D+2+%26+0+%5C%5C+0+%26+1+%5Cend%7Bbmatrix%7D+%5Ctag%7B14%7D
那么应变张量和能量就是
https://www.zhihu.com/equation?tex=%5Cvarepsilon+%3D+%5Cbegin%7Bbmatrix%7D+1+%26+0+%5C%5C+0+%26+0+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cpsi+%3D+%5Cmu+%2B+%5Cfrac%7B%5Clambda%7D%7B2%7D+%5Ctag%7B15%7D
很好能量算出来了，我们就要根据能量的导数来算力了。先把stress的计算公式甩上来。先让能量对变形梯度求导得到应力，然后应力对结点求导得到真正的力。
接下来要介绍一种很特别的应力(stress)。普通的应力计算公式如下
https://www.zhihu.com/equation?tex=%5Clim+_%7B%5CDelta+S+-%3E+P%7D+%5Cfrac%7B%5CDelta+F%7D%7B%5CDelta+A%7D+%5Ctag%7B16%7D
大意就是应力等于力除以面积。当面积越小越准确。然而我们之前说了弹性物体的两种状态，也就是初始未变形的状态，以及现在已经变形的状态。我们可以很方便给初始未变形的状态上的某个区域某个大小的力，但是对于已经变形的弹性物体，我们就很难计算到底是要给哪个区域，哪些顶点施加力。
因此，我们还是直接说，我们施加的力，是基于初始的未变形的状态的，那么这种应力就是First Piola Kirchhoff stress。
这种解释来自"Continuum Mechanics and Linear Elasticity - An Applied Mathematics Introduction by Ciprian D. Coman"第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 力的另一种计算方式是能量密度对变形梯度求导
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+%5CPsi%7D%7B%5Cpartial+%5Cbold+F%7D+%3D+%5Cbold+P%28%5Cbold+F%29++%5Ctag%7B17%7D
Psi 是刚才算出的能量密度，F是变形梯度，P 是first piola krichhoff stress 。对于线弹性力来说，求导如下
https://www.zhihu.com/equation?tex=%5Cbold+P%28%5Cbold+F%29+%3D+%5Cmu%28%5Cbold+F+%2B+%5Cbold+F%5ET+-+2%5Cbold+I%29+%2B+%5Clambda+tr%28%5Cbold+F++-+%5Cbold+I%29%5Cbold+I+%5Ctag%7B18%7D
并且计算物理中关于这种矩阵求导的非常非常多。在我的代码是这么写的
piola = mu * (F + F.T - 2 * np.identity(2)) + la * (F - 1 + F - 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 + epsilon + epsilon;
P = epsilon * btScalar(2) * m_mu + btMatrix3x3::getIdentity() * m_lambda * trace;
}再手算一遍，假设mu和lambda都是2，那么线弹性能量就是3，那么第一piola krichhoff应力
https://www.zhihu.com/equation?tex=%5Cbold+P%28%5Cbold+F%29+%3D+2%5Cmu+%5Cbold+E+%2B+%5Clambda+tr%28%5Cbold+E%29%5Cbold+I+%5C%5C+%3D+%282%5Cmu+%5Cbegin%7Bbmatrix%7D+1+%26+0+%5C%5C+0+%26+0%5Cend%7Bbmatrix%7D+%2B+%5Clambda%5Cbegin%7Bbmatrix%7D+1+%26+0+%5C%5C+0+%26+1%5Cend%7Bbmatrix%7D%29+%3D+%5Cbegin%7Bbmatrix%7D+6+%26+0+%5C%5C+0+%26+2%5Cend%7Bbmatrix%7D+%5Ctag%7B19%7D
节点力的计算
一个三角形有三个顶点，为了动量守恒，它们的节点力相加必须为零，也就是
https://www.zhihu.com/equation?tex=%5Cvec+f_0+%2B+%5Cvec+f_1+%2B+%5Cvec+f_2+%3D+0+%5Ctag%7B20%7D
节点力的来源于能量对节点位置求导，我们先将能量变成能量密度乘以面积或体积，然后将偏微分分成能量密度对变形梯度求导，以及变形梯度对节点位置求导。前者就是之前的first piola kirchhoff力，后者根据变形梯度的定义更可以直接得到。
https://www.zhihu.com/equation?tex=-+%5Cvec+f+%3D+%5Cfrac%7B%5Cpartial+E%7D%7B%5Cpartial+%5Cbold+x%7D+%3D+W%5Cfrac%7B%5Cpartial+%5CPsi%7D%7B%5Cpartial+%5Cbold+x%7D%3D+%5Cfrac%7B%5Cpartial+%5CPsi%7D%7B%5Cpartial+%5Cbold+F%7D%5Cfrac%7B%5Cpartial+%5Cbold+F%7D%7B%5Cpartial+%5Cbold+x%7D+%3D+W+%5Cbold+P+%5Cbold+D_m%5E%7B-T%7D+%5Ctag%7B21%7D
但是我们通常这么写，先算出第一个和第二个节点的力，然后根据动量守恒算出第一个节点的力
https://www.zhihu.com/equation?tex=%5Cbold+H+%3D+%5Cbegin%7Bbmatrix%7D+%5Cvec+f_1+%26+%5Cvec+f_2+%5Cend%7Bbmatrix%7D+%3D+-W%5Cbold+P%28%5Cbold+F%29%5Cbold+D_m%5E%7B-T%7D+%5Cqquad+%5Cvec+f_0+%3D+-%5Cvec+f_1+-+%5Cvec+f_2+%5Ctag%7B22%7D
推导如下
这种节点力还要不同的写法，比如在“Robust Quasistatic Finite Elements and Flesh Simulation” 是这么写的
https://www.zhihu.com/equation?tex=%5Cbold+g_i+%3D+-%5Cbold+P+%28A_1+%5Cbold+N_1+%2B+A_2+%5Cbold+N_2+%2B+A_3+%5Cbold+N_3%29%2F3+%5Cqquad+%5Cbold+g_0+%3D+-%28%5Cbold+g_1+%2B+%5Cbold+g_2+%2B+%5Cbold+g_3%29+%5Ctag%7B23%7D
其中AN是三角形的法向量根据面积的权重，P是之前求出的piola应力。也可以写成另一种方法
https://www.zhihu.com/equation?tex=%5Cbold+G+%3D+%5Cbold+P+%5Cbold+B+%5Cqquad+%5Cbold+G+%3D+%28%5Cbold+g_1%2C%5Cbold+g_2%2C+%5Cbold+g_3%29+%5Cqquad+%5Cbold+B_m+%3D+%28%5Cbold+b_1%2C%5Cbold+b_2%2C%5Cbold+b_3%29+%3D+-V%5Cbold+D_m%5E%7B-T%7D+%5Ctag%7B24%7D
这些公式在我的代码中是这么写的
gradC1 = np.array(,H])
gradC2 = np.array(,H])
gradC0 = - gradC1 - gradC2
node_force,:] += gradC0
node_force,:] += gradC1
node_force,:] += gradC2其中bullet3库是这么写的，代码在https://github.com/bulletphysics/bullet3/blob/master/src/BulletSoftBody/btDeformableLinearElasticityForce.h#L263
btMatrix3x3 force_on_node123 = psb->m_tetraScratches.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;
btSoftBody::Node* node1 = tetra.m_n;
btSoftBody::Node* node2 = tetra.m_n;
btSoftBody::Node* node3 = tetra.m_n;
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 -= scale1 * force_on_node0;
force -= scale1 * force_on_node123.getColumn(0);
force -= scale1 * force_on_node123.getColumn(1);
force -= 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 = H(0, 0);
J = H(0, 1);
J = H(0, 2);
J = H(1, 0);
J = H(1, 1);
J = H(1, 2);
J = H(2, 0);
J = H(2, 1);
J = H(2, 2);
J = -J - J - J;
}继续手推，
https://www.zhihu.com/equation?tex=%5Cbold+H+%3D+-W%5Cbold+P%28%5Cbold+F%29d%5Cbold+X%5E%7B-T%7D+%3D+-%5Cfrac%7B1%7D%7B2%7D%5Cbegin%7Bbmatrix%7D+6+%26+0+%5C%5C+0+%26+2%5Cend%7Bbmatrix%7D+%5Cbegin%7Bbmatrix%7D+1+%26+0+%5C%5C+0+%26+1%5Cend%7Bbmatrix%7D+%3D+%5Cbegin%7Bbmatrix%7D+-3+%26+0+%5C%5C+0+%26+-1%5Cend%7Bbmatrix%7D++%5Ctag%7B25%7D+
也就是
https://www.zhihu.com/equation?tex=%5Cvec+f_0+%3D+%5Cbegin%7Bbmatrix%7D+3+%5C%5C+1%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cvec+f_1+%3D+%5Cbegin%7Bbmatrix%7D+-3+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cvec+f_2+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+-1%5Cend%7Bbmatrix%7D+%5Ctag%7B26%7D
最后，再整理一遍流程
算法解释如下
算法第二行是初始化，只需要对元素，也就是每个三角形(或四面体)执行一遍
算法第三行，计算每个元素的Dm或者叫dX，用于之后计算变形梯度
算法第四行，对Dm求逆得到Bm，用于之后计算变形梯度
算法第五行，计算三角形的面积(或四面体的体积)，用于之后计算节点力
算法第八行，外层主循环，每个时间步都要执行一遍
算法第九行，将节点力设置为零
算法第十行，内层循环，遍历所有三角形/(四面体)
算法第十一行，计算每个元素的Ds或者叫dx，用于计算变形梯度
算法第十二行，计算变形梯度
算法第十三行，计算first piola kirchhoff 应力
算法第十四到第十六行，计算真正的节点力
那么最简单的二维线弹性显式求解如下
import numpy as np
# 初始化三角形初始位置
node_pos = np.array([,,],dtype = float)
# 顶点的位置梯度，对应上面算法第三行
Ds = np.array([ - node_pos,node_pos -
node_pos],
- node_pos,node_pos -
node_pos]])
# 求逆，用于准备计算形变梯度，对应上面算法第四行
minv = np.linalg.inv(Ds)
# 假设某一时刻，三角形变化到了这样的位置
node_pos = np.array([,,],dtype = float)
time = 0
timeFinal = 100
areat = np.zeros((timeFinal))
while(time < timeFinal):
time += 1
# 形变梯度中的分子，对应上面算法第11行
Ds_new = np.array([ - node_pos,node_pos
- node_pos],
- node_pos,node_pos -
node_pos]])
# 形变梯度，对应上面算法第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*strain + strain*strain
doubleInner += strain*strain + strain*strain
# 线弹性的能量计算公式
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 - 1 + F - 1) * np.identity(2)
# 三角形面积
area = 0.5
# 计算hessian矩阵对应上面算法第14行
H = - area * np.dot(piola,minv.transpose())
# 计算节点力，对应上面算法第15和第16行
gradC1 = np.array(,H])
gradC2 = np.array(,H])
gradC0 = - gradC1 - gradC2
invMass = 1
dt = 0.1
# 判断是否收敛
sumGradC = invMass * (gradC0**2 + gradC0**2)
sumGradC += invMass * (gradC1**2 + gradC1**2)
sumGradC += invMass * (gradC2**2 + gradC2**2)
if sumGradC < 1e-10:
break
# 校正位置，方法来源于PositionBasedDynamics
node_pos += dt * energy / sumGradC * invMass * gradC0
node_pos += dt * energy / sumGradC * invMass * gradC1
node_pos += dt * energy / sumGradC * invMass * gradC2
areat = 0.5 * (node_pos * (node_pos - node_pos)
+ node_pos * (node_pos - node_pos)
+ node_pos * (node_pos - node_pos))
但是不要高兴得太早了，如果使用上面的代码，你会发现甚至算不出正确的结果，无法收敛。因为线弹性忽略了一个非常重要的东西，旋转。因旋转并未造成面积改变。举个例子，仍然是之前那个三角形，不过现在逆时针旋转90度
https://www.zhihu.com/equation?tex=%5Cbold+x_0+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_1+%3D+%5Cbegin%7Bbmatrix%7D+0+%5C%5C+1%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_2+%3D+%5Cbegin%7Bbmatrix%7D+-1+%5C%5C+0%5Cend%7Bbmatrix%7D+%5Ctag%7B27%7D
此时dx和变形梯度如下
https://www.zhihu.com/equation?tex=d%5Cbold+x+%3D+%5Cbegin%7Bbmatrix%7D+0+%26+-1+%5C%5C+1+%26+0+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D+0+%26+-1+%5C%5C+1+%26+0+%5Cend%7Bbmatrix%7D+%5Ctag%7B28%7D
那么应变张量和能量就是
https://www.zhihu.com/equation?tex=%5Cvarepsilon+%3D+%5Cbegin%7Bbmatrix%7D+-1+%26+-1+%5C%5C+1+%26+-1+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cpsi+%3D+4%5Cmu+%2B+4+%5Cfrac%7B%5Clambda%7D%7B2%7D+%5Ctag%7B29%7D
因为mu和lambda不为零，那么能量也必不为零，但是，你所做的，仅仅是旋转了物体而已，却也凭空产生能量。所以用线弹性模型很少在各种物理模拟库中见到。
STVK模型
二维显式stvk的代码在
和gitee上的
线弹性模型失败的原因来源于它认为旋转也是变形，所以会产生能量。因此我们要重新设计一个能量公式，让旋转不产生能量。最简单的非线性模型就是St. Venant-Kirchhoff能量模型，它长下面这样
https://www.zhihu.com/equation?tex=%5CPsi%28%5Cbold+F%29+%3D+%5Cmu+%5Cbold+E+%3A%5Cbold+E+%2B+%5Cfrac%7B%5Clambda%7D%7B2%7Dtr%5E2%28%5Cbold+E%29+%5Ctag%7B30%7D
同时把测量应变的矩阵换成了green lagrange strain，如下
https://www.zhihu.com/equation?tex=%5Cbold+E+%3D+%5Cfrac%7B1%7D%7B2%7D%28%5Cbold+F%5ET+%5Cbold+F+-+%5Cbold+I%29+%5Ctag%7B31%7D
这种应变测量方法也不会把旋转当作变形。这种模型的缺点之后再说。现在只要知道使用这种模型，那么物体旋转时就不会产生能量，甩线弹性几条街就行了。
此时stvk模型的first piola kirchhoff 应力计算如下
https://www.zhihu.com/equation?tex=%5Cbold+P%28%5Cbold+F%29+%3D+%5Cbold+F%282%5Cmu+%5Cbold+E+%2B+%5Clambda+tr%28%5Cbold+E%29%5Cbold+I%29+%5Ctag%7B32%7D
仍然使用之前的代码，然后仅仅需要改变三行代码，就能将线弹性模型换成stvk模型
# Green Strain，也就是E
E = (np.dot(F,F.T)- np.identity(2)) * 0.5
# Stvk的能量计算公式
energy = doubleInner * mu + la / 2 * (E + E)**2
#first piola kirchhoff stress
piola = np.dot(F, 2 * mu * E + la * (E + E) * np.identity(2))最后计算的各顶点位置和形变梯度如下
https://www.zhihu.com/equation?tex=%5Cbold+x_0+%3D%5Cbegin%7Bbmatrix%7D+0.427+%5C%5C+-0.072%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_1+%3D%5Cbegin%7Bbmatrix%7D+1.395+%5C%5C+0.177%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+x_2+%3D%5Cbegin%7Bbmatrix%7D+0.177+%5C%5C+0.895%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D+0.968+%26+-0.250+%5C%5C+0.250+%26+0.968%5Cend%7Bbmatrix%7D+%5Ctag%7B33%7D
接下来看看有哪些开源库也写了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; // 节点位置
const double b = nodes;
const double c = nodes;
const double d = nodes;
...
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 = -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 = -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三维隐式方法。
下一篇： 好棒，支持一下 很专业 当时偷懒了 推了一个简单拟合函数 实现 但简单只需要2天实现3d 胶体 https://zhuanlan.zhihu.com/p/361126215 不错的工作 矩阵相乘运算 展开到 元素的加减乘除能提高效率？ 从节点位置矩阵，到计算出节点力的过程中，要经过相当多的矩阵乘法，转置，相加相减。其中有些矩阵是对称的，有的矩阵有很多零，导致很多乘法和加法是重复的，还有很多乘零和加零的操作，尤其是使用隐式方法的时候。所以自动生成的代码就是将这些重复的没必要的加法和乘法去掉了。 这是非线性连续介质力学里面的东西，游戏现在要做柔性体了么，你们游戏方向有要招聘力学方向博士么 也许你可以试试泽森科工这家公司，他们可能需要力学方向的，网址在zenus 需要的 这家公司据我了解是需要程序员多一些，核心力学算法都是现成的，好像不是那么需要力学方向的
页:
[1]
2