[游戏物理探秘]三维有限元的隐式积分方法

pc8888888 · 发表于 2021-12-20 09:20

我们慢慢从一维推导起吧。继续手写算例，找开源库。
一维显式方法

假设有两个点，它们都是一维的，起始位置为

$x_0 = 1 \qquad x_1 =3 \tag{1}$
那么对初始位置微分得

$d\bold X = \bold D_m = x_1 - x_0 = 2\qquad d\bold X^{-1} = \bold D_m^{-1} = 0.5 \tag{2}$
假设现在两个顶点位置运动到了

$x_0 = 1 \qquad x_1 = 5 \tag{3}$
那么继续对位置微分，并算出变形梯度

$d\bold x = \bold D_s = x_1 - x_0 = 4 \qquad \bold F = \frac{d\bold x}{d\bold X} = 2 \tag{4}$
一维的变形梯度，就是长度的缩放倍数。初始两个点的距离是2，现在变成了4，那么相当于距离变成了两倍，那么变形梯度就是2。那么Green应变和第一pk应力计算如下，假设mu = lambda = 2

$\bold E = (\bold F^T \bold F - \bold I)/2 = \frac{3}{2} \qquad \bold P = \bold F(2\mu \bold E + \lambda tr(\bold E)\bold I) = 18 \tag{5}$
如果你没搞清楚Green应变和第一pk应力到底是什么，也没关系。大意就是只要它们的数值越大，表明物体变形的越厉害。当它们的数值是零，代表物体没变形。如果你真的想弄清楚这两个玩意是什么也很简单，把你在网上或本地图书馆找到的所有有关弹性力学和连续介质力学的资料和书都看一遍就行了。
接下计算节点力矩阵H

$H = -W \bold P \bold D_m^{-T} = - 2 * 18 * 0.5 = -18 \tag{6}$
其推导方式如下

$\bold f_1 = \bold H_{00} = -18 \qquad \bold f_0 = \bold f_1 = 18 \tag{7}$
显式方法的话，只需要将力乘上时间，除以质量就可以放到速度上，然后再让速度乘上时间就可以加到位置上。但是显式方法说白了就是纯靠猜，一点点靠近结果，不仅慢还不稳定，时间步长一大就容易猜过头。因此接下来直接上隐式。
一维隐式方法

一维隐式代码完整代码在github
或gitee上

首先把显式动力学公式弄上来。下一时刻的位置等于这一时刻的位置，加上下一时刻的速度乘上时间差

$\bold x ^{t+1} = \bold x ^t + \Delta t \bold x^{t+1} \tag{8}$
然后，下一时刻的速度，等于这一时刻的速度加上这一时刻的力乘上时间差除以质量

$\bold M \bold v^{t+1} = \bold M \bold v^{t} + \Delta t F(\bold x^{t}) \tag{9}$
M 就是质量矩阵。那么显式的积分就是只要知道这一时刻的位置和速度，就能知道下一时刻的速度，进而知道下一时刻的位置。但是这种方法精度很差，所以我们要将第二个公式换为

$\bold M \bold v^{t+1} = \bold M \bold v^{t} + \Delta t F(\bold x^{t+1}) \tag{10}$
也就是下一时刻的速度，等于这一时刻的速度加上下一时刻时刻的力乘上时间差除以质量
但是下一时刻的力和下一时刻的速度我们都不知道，这样一来就成了先有鸡还是现有蛋的问题。或许那个工作经验的表情包更能表达这种关系

隐式一阶欧拉积分，就是将第二个动力学公式近似为

$\bold M \bold v^{t+1} = \bold M \bold v^t + \Delta t F(\bold x^t + \Delta t \bold v^{t+1}) \approx \bold M \bold v^t + \Delta t F(\bold x^T) + \Delta t^2 \bold K|_{\bold x^t}\cdot \bold v^{t+1} \tag{11}$
这里的矩阵K就是力对速度的位置求导，或者叫tangent stiffness matrix。重新排列上式

$(\bold M - \Delta t^2 \bold K|_{\bold x^t})\bold v^{t+1} = \bold M \bold v^t + \Delta t F(\bold x^t) \tag{12}$
那么这样就成了一个线性方程组，其中

$\bold A \bold x = \bold b \tag{13}$
其中

$\bold A = \bold M - \Delta t^2 \bold K|_{\bold x^t} \qquad \bold x = \bold v^{t+1} \qquad \bold b = \bold M \bold v^t + \Delta t F(\bold x^t) \tag{14}$
上面公式来自&#34;Point Based Animation of Elastic, Plastic and Melting Objects&#34;。另一种展示方式是vegafem中的ImplicitBackwaradEulerSparese.cpp中对此解释非常清楚
// build effective stiffness:
// Keff = M + h D + h^2 * K
// compute force residual, store it into aux variable qresidual
// this is semi-implicit Euler
// qresidual = h * (-D qdot - fint + fext - h * K * qdot))
// for fully implicit Euler
// qresidual = M (qvel_1 - qvel) + h * (-D qdot - fint + fext - K * (q_1 - q + h qdot) ))
不过每个矩阵分别长什么样子呢？如果是n 个 1维点，那么A 矩阵就是n x 1 = n 行 n 列，
如果是n 个二维点，那么A 矩阵就是n x 2 = 2n 行 2n 列，如果是n 个三维点，那么A 矩阵就是n x 2 = 3n 行 3n 列。M 是质量矩阵，也是对角矩阵，代表每个顶点的质量。F 是顶点力，我们在显式步骤里计算出来了。最难搞定的是那个K 矩阵，也就是力对顶点位置求导。
并且以后这种向量对向量求导，矩阵对向量求导的公式多了去了。
还是回到一开始显式方法的那个例子，显然弹簧两个点离得越远，它们之间的力越大。但是，当第二个顶点每朝第一个顶点远离一个单位时，力究竟应该增加多大呢？这其实就是求导，结果就是

$\delta \bold f = \frac{\partial \bold f}{\partial \bold x}|_{\bold x^*} \cdot \delta \bold x \tag{15}$
我们虽然不知道力对顶点求导是多大，但我们可以反推回去，力和hessian 矩阵有关，也就是

$\delta \bold H = \begin{bmatrix} \delta \bold f_1\end{bmatrix} \tag{16}$
而hessian矩阵又和piola 力有关

$\delta \bold H = - W \delta \bold P (\bold F;\delta \bold F)\bold D_m^{-T} \tag{17}$
也就是我们要先把piola力对位置求导算出来，但piola力又和应变E有关，应变又和变形梯度有关，所以我们首先应当计算变形梯度对位置求导。变形梯度及其导数是这么写得

$\bold F = \bold D_s \bold D_m^{-1} \qquad \delta \bold F = (\delta \bold D_s)\bold D_m^{-1} \tag{18}$
上面这个式子在不同的讲义上样子还不同。比如这里我参考的是“Dynamic Deformables: Implementation and Production Practicalities  ”。但在太极图形课上https://www.bilibili.com/video/BV1eY411x7mK?p=5&t=390.0，它长这样

接下来就要手算了。还记得一维Ds 长什么样吗，非常简单，就是

$\bold D_s = [x_1 - x_0] \tag{19}$
那么我们现在干脆把它算出来吧。很多教程都不怎么举例子，导致我有时候想了个例子想验证一下都不知道结果是否正确。但是本系列教程肯定会多多举例子吧。

$\frac{\partial \bold D_s}{\partial x_0} = - 1 \qquad \frac{\partial \bold D_s}{\partial x_1} = 1 \tag{20}$
请想清楚这里求导的意义，这代表x0 每前进1，那么Ds 就会减少1。而如果x1每前进1，那么Ds就增加1。而Ds正是x1与x0之间的距离。并且求导之后得到的东西是一个有2x1x1x1 = 2个元素的四阶张量。第一个2是因为一个线段有两个顶点，第一个1是因为这两个顶点的位置都是一维的，之后两个1也是因为一维的关系。
至于四阶张量，把它想象成有四个坐标轴的矩阵就行了。或者别把它看成张量，就想象每个顶点每个坐标轴都对应一个1x1的矩阵。
Ds 对位置求导后，不是1就是-1。哪怕三维也是如此。因此可以直接把它弄成全局变量，或者甚至不要Ds求导写出来也可。然后是变形梯度求导

$\delta \bold F = \delta \bold D_s \bold D_m^{-1} \qquad \frac{\partial \bold F}{\partial x_0} = \frac{\partial \bold D_s}{\partial x_0} \bold D_m^{-1} = -\frac{1}{2} \qquad \frac{\partial \bold F}{\partial x_1} = \frac{1}{2} \tag{21}$
这里变形梯度与位置仍然是线性的，如果x0能够前进2，那么变形梯度计算如下

$\bold F^{t+1} = \bold F^{t} + \Delta x_0 \frac{\partial \bold F}{\partial x_0} = 2 + (2 * (-\frac{1}{2})) = 1 \tag{22}$
物理意义很明确，就是线段的长度就是初始长度，变形梯度也变成单位矩阵了。变形梯度求导，在这里也是2x1x1x1的四阶张量。
然后应变求导

$\delta \bold E = \frac{1}{2}(\delta \bold F^T \bold F + \bold F^T \delta \bold F) \qquad \frac{\partial \bold E}{\partial x_0} = - 1 \qquad \frac{\partial \bold E}{\partial x_1} = 1 \tag{23}$
这里应变与位置的关系不是线性的，但是当x0前进1，可以用一阶导粗略得算出应变减少3/4，虽然有误差，但是误差是导数不够精确造成的。物理意义仍然没错。
然后算第一pk应力的导，使用的是stvk模型

$\delta \bold P = \delta \bold F [2\mu \bold E + \lambda tr(\bold E)\bold I] + \bold F(2\mu \delta \bold E + \lambda tr(\bold \delta \bold E)\bold I) \\ \frac{\partial \bold P}{\partial x_0} = -\frac{1}{2}(2*2*\frac{3}{2} + 2*\frac{3}{2}) + 2(2*2* (-1) + 2 * (-1)) = -33/2 \\ \frac{\partial \bold P}{\partial x_1} = 33/2 \tag{24}$
这里第一pk力与位置的关系仍然不是线性的。但是按照方法可以大致看出来，每当x0前进，第一pk应力就会减少。反之亦然。
然后算节点力矩阵的hessian 矩阵

$\delta \bold H = - W (\delta \bold P)\bold D_m^{-T} \qquad \frac{\partial \bold H}{\partial x_0} = 33/2 \qquad \frac{\partial \bold H}{\partial x_1} = -33/2 \tag{25}$
这里的dD, dF, dE, dP, dH 都是2x1x1x1 的四阶张量。节点力对位置求导我们实际上就算出来了

$\delta \vec f_1 = \delta \vec h_0 = \frac{\partial \bold H}{\partial x_0} = 33/2\qquad \delta \vec f_0 = -\delta \vec h_0 = -33/2 \tag{26}$
务必记住节点力对位置求导的物理意义。在我们的例子中，每当顶点0前进1，那么弹簧长度变短，力变小，具体幅度是-33/2。每当顶点1前进1，那么弹簧长度变长，力变大，具体增幅是33/2。
接下来就要用这个组装矩阵了。不过我思考了很久，组装矩阵这部分应该怎么写。如果你完全没接触过有限元，对K矩阵也不熟悉的话，最好的办法是一步一步看所有你能找到的有关有限元和数值仿真的书和资料，用数值计算从杆单元，梁单元，壳单元，薄板单元之类的慢慢算起。并且这个过程需要花费大量时间和精力，并且极其无聊，一眼望去全是公式和代码，每天都在算这个杆弯了多少度，那个杆又伸长了多少，毫无视觉效果。我相信很多人看不懂games103和太极图形课，也是因为卡在了有限元和数值仿真上。不过跨过这道坎就可以模拟出很多好玩的东西了。
当然如果你有其他什么好想法，也欢迎和我交流。线上交流效率不高，也欢迎来长沙或别的地方和我当面交流。
不过现在我还是尽量尝试解释清楚，怎么把节点力放到K矩阵里面去。
首先，现在我们的K矩阵的行数，等于节点的数量乘以节点的自由度，或者节点的数量乘以维度。线性方程的每一行都是一个含有未知量的等式。我们每次填写矩阵的一列。首先填写K矩阵的第一列，需要用到之前说到的2x1个1x1矩阵中的第一个1x1矩阵。
把它放在对应位置

$\qquad \bold K = \begin{bmatrix}* & * \\ \delta \vec h_0 & *\end{bmatrix}= \begin{bmatrix} * & * \\ 33/2 & *\end{bmatrix} \tag{27}$
但是我们这里一列之中还剩下对角线上的元素没有填写，此时需要动量守恒公式。每一列的所有节点力相加得零，也就是

$\qquad \bold K = \begin{bmatrix}-\delta \vec h_0 & * \\ \delta \vec h_0 & *\end{bmatrix}= \begin{bmatrix} -33/2 & * \\ 33/2 & *\end{bmatrix} \tag{28}$
完成第一列。注意我们只用到了公式7中的dhdx0，而没有用到dhdx1。接下来我们要用dhdx1把第二列添上。

$\qquad \bold K = \begin{bmatrix}* & *\\ * & \vec h_1 \end{bmatrix}= \begin{bmatrix} * & *\\ * & -33/2\end{bmatrix} \tag{29}$
然后填写最后一个元素

$\qquad \bold K = \begin{bmatrix}* & -\vec h_1 \\ * & \vec h_1\end{bmatrix}= \begin{bmatrix} * & 33/2 \\ * & -33/2\end{bmatrix} \tag{30}$
最后填充如下

$\bold K = \begin{bmatrix} -33/2 & 33/2\\ 33/2 & -33/2\end{bmatrix} \tag{31}$
再回忆一遍我们要求的线性方程组

$\bold A \Delta \bold x = (\bold M - \bold K dt^2)\Delta \bold x =\Delta t \bold f \tag{32}$
简单起见，假设dt = 1, 质量为1，那么写出来就是

$(\begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix} - \begin{bmatrix} -33/2 & 33/2\\ 33/2 & -33/2\end{bmatrix})\begin{bmatrix} \Delta x_0 \\ \Delta x_1\end{bmatrix} = \begin{bmatrix}18 \\ -18\end{bmatrix} \tag{33}$
最后解的下一时刻的速度就是

$\begin{bmatrix} \Delta x_0 \\ \Delta x_1\end{bmatrix} = \begin{bmatrix} 0.529\\ -0.529\end{bmatrix} \tag{34}$
分析一下，线段上的两个点现在处在伸长状态，所以它们现在为了减小能量，需要相互靠近，很合理。接下来再整理一遍流程

算法解释如下
初始化步骤省略
算法第一行，外层主循环，每个时间步都要执行一遍
算法第三行，内层循环，遍历所有三角形/(四面体)
算法第四行，计算每个元素现在的Ds或者叫dx，用于计算变形梯度
算法第五行，Ds对节点位置求导
算法第六行，计算变形梯度
算法第七行，计算变形梯度对顶点位置求导
算法第八行，计算first piola kirchhoff 力对顶点位置求导
算法第九到第十一行，计算节点力对顶点位置求导
while(time < timeFinal):

Kmat = np.zeros((2,2))
# 算法第四行
Ds = x1 - x0
# 算法第六行
F = Ds * Dminv
E = (F * F - 1)/2
trE = E
P = F*(2*mu*E + la*trE)
H = - W * P * Dminv
f1 = H
f0 = - f1

# 算法第五行
dD0 = -1
dD1 = 1
# 算法第七行
dF0 = dD0 * Dminv
dF1 = dD1 * Dminv

dE0 = (dF0 * F + F * dF0) / 2
dE1 = (dF1 * F + F * dF1) / 2
# 算法第八行
dP0 = dF0 * (2*mu*E + la*trE) + F*(2*mu*dE0 + la*dE0)
dP1 = dF1 * (2*mu*E + la*trE) + F*(2*mu*dE1 + la*dE1)
# 算法第九到十一行
dH0 = - W * dP0 * Dminv
Kmat[1,0] = dH0
Kmat[0,0] =  - dH0
dH1 = - W * dP1 * Dminv
Kmat[0,1] = - dH1
Kmat[1,1] = dH1

A = mass * np.identity(2) - Kmat * dt * dt
b = np.array([v0,v1]) * mass + dt * np.array([f0,f1])
x = np.dot(np.linalg.inv(A),b)

v0 = x[0]
v1 = x[1]
x0 = x0 + v0 * dt
x1 = x1 + v1 * dt

xt[0,time] = x0
xt[1,time] = x1
time += 1二维隐式

先看看二维隐式方法能够有多弹的表现

完整演示视频在下方。你也可以在我的unity代码中的Implicit2D中找到这个场景。
python代码则在
和
仍然假设三角形三个点初始位置，以及Dm如下

$\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} \qquad \tag{25}$
那么Dm如下

$\bold D_m = \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{26}$
现在假设三角形三个顶点移动到如下位置

$\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} \qquad \tag{27}$
那么Ds以及变形梯度为

$\bold D_s = \begin{bmatrix} \bold x_1 - \bold x_0 \,|\, \bold x_2 - \bold x_0\end{bmatrix} = \begin{bmatrix} 2 & 0 \\ 0 & 1\end{bmatrix} \qquad \bold F = \bold D_s \bold D_m^{-1} = \begin{bmatrix} 2 & 0 \\ 0 & 1\end{bmatrix} \tag{28}$
代码如下
p0 = node_pos[element_idx[ie,0]]
p1 = node_pos[element_idx[ie,1]]
p2 = node_pos[element_idx[ie,2]]
dx = np.array([[p1[0] - p0[0],p2[0] - p0[0]],
                     [ p1[1] - p0[1],p2[1] - p0[1]]])

# 形变梯度
F = np.dot(dx,element_minv[ie])格林应变和第一pk应力如下，仍然假设mu = 2, lambda = 2

$\bold E = (\bold F^T \bold F - \bold I)/2 = \begin{bmatrix} 3/2 & 0 \\ 0 & 0\end{bmatrix} \qquad \bold P = \bold F(2\mu \bold E + \lambda tr(\bold E)\bold I) = \begin{bmatrix} 18 & 0 \\ 0 & 3\end{bmatrix} \tag{29}$
并且hessian 矩阵

$\bold H = - area \bold P(\bold F)\bold D_m^{-T} = \begin{bmatrix} -9 & 0 \\ 0 & -3/2\end{bmatrix} \tag{30}$
代码如下
E = (np.dot(F.T,F)- np.identity(2)) * 0.5
piola = np.dot(F, 2 * mu * E + la * (E[0,0] + E[1,1]) * np.identity(2))
doubleInner = E[0,0]*E[0,0] + E[1,0]*E[1,0] + E[0,1]*E[0,1] + E[1,1]*E[1,1]
energy = doubleInner * mu + la / 2 * (E[0,0] + E[1,1])**2
H = - area * np.dot(piola,element_minv[ie].transpose())那么可以计算三个节点力

$\vec f_1 = \begin{bmatrix} -9 \\ 0 \end{bmatrix} \qquad \vec f_2 = \begin{bmatrix} 0 \\ -3/2 \end{bmatrix} \qquad \vec f_0 = - \vec f_1 - \vec f_2 = \begin{bmatrix} 9 \\ 3/2\end{bmatrix} \tag{31}$
稍微判断一下就知道算出是正确的，因为现在三角形面积变得比初始面积大，所以内能就要想着缩小面积。那么第零个节点，也就是左下角那个节点，开始向右上角移动。其它同理。
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然后开始求导。如果我们重新划分一下现在的顶点位置

$\bold x = \begin{bmatrix} \bold x_0 \\ \bold x_1 \\ \bold x_2 \end{bmatrix} = \begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ .. \\ x_5 \end{bmatrix} \quad \bold x_0 = \begin{bmatrix} x_0 \\ x_1 \end{bmatrix}\quad \bold x_1 = \begin{bmatrix} x_2 \\ x_3 \end{bmatrix}\quad \bold x_2 = \begin{bmatrix} x_4 \\ x_5 \end{bmatrix} \tag{32}$
然后想明白Ds 由什么组成

$\bold D_s = [\bold x_1 - \bold x_0 \quad |\quad \bold x_2 - \bold x_0] = \begin{bmatrix} x_2 - x_0 & x_4 - x_0\\ x_3 - x_1 & x_5 - x_1\end{bmatrix} \tag{33}$
此时的dD是3x2x2x2的四阶张量。因此我们可以很方便得求导这个3x2x2x2的张量

$\frac{\partial \bold D_s}{\partial x _0} = \begin{bmatrix} -1 & - 1 \\ 0 & 0\end{bmatrix} \qquad \frac{\partial \bold D_s}{\partial x _1} = \begin{bmatrix} 0 & 0 \\ -1 & -1\end{bmatrix} \\ \frac{\partial \bold D_s}{\partial x _2} = \begin{bmatrix} 1 & 0 \\ 0 & 0\end{bmatrix} \qquad \frac{\partial \bold D_s}{\partial x _3} = \begin{bmatrix} 0 & 0 \\ 1 & 0\end{bmatrix} \\ \frac{\partial \bold D_s}{\partial x _4} = \begin{bmatrix} 0 & 1 \\ 0 & 0\end{bmatrix} \qquad \frac{\partial \bold D_s}{\partial x _5} = \begin{bmatrix} 0 & 0 \\ 0 & 1\end{bmatrix} \tag{34}$
代码如下
dD = np.zeros((6,2,2))
dD[0,:,:] = np.array([[-1,-1],[0,0]])
dD[1,:,:] = np.array([[0,0],[-1,-1]])
dD[2,:,:] = np.array([[1,0],[0,0]])
dD[3,:,:] = np.array([[0,0],[1,0]])
dD[4,:,:] = np.array([[0,1],[0,0]])
dD[5,:,:] = np.array([[0,0],[0,1]])既然dD 是 4阶张量有24个数字，那么dF 也是一样是四阶张量有24个数字。此时的dF 与dD的数值是一样的

$\frac{\partial \bold F}{\partial x _0} = \begin{bmatrix} -1 & - 1 \\ 0 & 0\end{bmatrix} \qquad \frac{\partial \bold F}{\partial x _1} = \begin{bmatrix} 0 & 0 \\ -1 & -1\end{bmatrix} \\ \frac{\partial \bold F}{\partial x _2} = \begin{bmatrix} 1 & 0 \\ 0 & 0\end{bmatrix} \qquad \frac{\partial \bold F}{\partial x _3} = \begin{bmatrix} 0 & 0 \\ 1 & 0\end{bmatrix} \\ \frac{\partial \bold F}{\partial x _4} = \begin{bmatrix} 0 & 1 \\ 0 & 0\end{bmatrix} \qquad \frac{\partial \bold F}{\partial x _5} = \begin{bmatrix} 0 & 0 \\ 0 & 1\end{bmatrix} \tag{35}$
不过我觉得写出这样不好理解，因为变形梯度本来只是一个2x2矩阵，那么对x012345求导，本应该是矩阵中的每个元素对x012345求导，这样的物理就很好理解，如果x2 减去1，也就是最右边那个顶点的x坐标减去一，那么变形梯度就正好变成了单位矩阵，如下

$\hat {\bold F} = \bold F - (\Delta x_2 = 1)\frac{\partial \bold F}{\partial x_2} =\begin{bmatrix} 2 & 0 \\ 0 & 1\end{bmatrix} - \begin{bmatrix} 1 & 0 \\ 0 & 0\end{bmatrix} \tag{36}$
但现在这样写，求导完的结果在上面的式子却被分散在6个矩阵了。所以在“Dynamic Deformables: Implementation and Production Practicalities  ”这篇文章所附代码中，是这样写的

$\bold d = \bold D_m^{-1} = \begin{bmatrix}d_{00} & d_{01} \\ d_{10} & d_{11} \end{bmatrix}= \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \tag{37}$

$\frac{\partial \bold F}{\partial \bold x} = \begin{bmatrix} -d_{00}-d_{01} & 0 & d_{00} & 0 & d_{01} & 0 \\ 0 & -d_{00}-d_{01} & 0 & d_{00} & 0 & d_{01} \\ -d_{10}-d_{11} & 0 & d_{10} & 0 & d_{11} & 0 \\ 0 & -d_{10}-d_{11} & 0 & d_{10} & 0 & d_{11}\end{bmatrix} \tag{38}$
代码请见https://graphics.pixar.com/library/DynamicDeformablesSiggraph2020/ 的 ComputePFPx函数中。其结果和上面的是一样的。第一列对应dFdx0，第二列对应dFdx1，以此类推。

$\frac{\partial \bold E}{\partial x _0} = \begin{bmatrix} -2 & - 1 \\ -1 & 0\end{bmatrix} \qquad \frac{\partial \bold E}{\partial x _1} = \begin{bmatrix} 0 & -0.5 \\ -0.5 & -1\end{bmatrix} \\ \frac{\partial \bold E}{\partial x _2} = \begin{bmatrix} 2 & 0 \\ 0 & 0\end{bmatrix} \qquad \frac{\partial \bold E}{\partial x _3} = \begin{bmatrix} 0 & 0.5 \\ 0.5 & 0\end{bmatrix} \\ \frac{\partial \bold E}{\partial x _4} = \begin{bmatrix} 0 & 1 \\ 1 & 0\end{bmatrix} \qquad \frac{\partial \bold E}{\partial x _5} = \begin{bmatrix} 0 & 0 \\ 0 & 1\end{bmatrix}$
现在应变求导，但应变与位移并非线性关系，所以如果像刚才那样直接怼一阶导就会有误差。不过仔细观察大体上还能观察出来，那么就是这些矩阵相加的结果是零矩阵。然后dP没什么好写的，用调试器自己去看吧。不过得说说dH

$\frac{\partial \bold H}{\partial x _0} = \begin{bmatrix} 16.5 & 5.5 \\ 2 & 2\end{bmatrix} \qquad \frac{\partial \bold H}{\partial x _1} = \begin{bmatrix} 2 & 2 \\ 5.5 & 4.5\end{bmatrix} \\ \frac{\partial \bold H}{\partial x _2} = \begin{bmatrix} -16.5 & 0 \\ 0 & -2\end{bmatrix} \qquad \frac{\partial \bold H}{\partial x _3} = \begin{bmatrix} 0 & -2 \\ -5.5 & 0\end{bmatrix} \\ \frac{\partial \bold H}{\partial x _4} = \begin{bmatrix} 0 & -5.5 \\ -2 & 0\end{bmatrix} \qquad \frac{\partial \bold H}{\partial x _5} = \begin{bmatrix} -2 & 0 \\ 0 & -4.5\end{bmatrix}$
继续仔细观察这些矩阵。你能发现什么？因为如果自己不仔细观察的话，也没有那本书会教我怎么观察。嗯，那些教科书都是堆砌公式罢了。首先，hessian矩阵对顶点求导，所有元素相加的和仍然是零。其次，第一个矩阵的第一个元素是16.5，这和一维隐式的结果是一样的。
上面的写成代码如下
for i in range(6):
dF[i,:,:] = np.dot(dD[i,:,:],minv)
d_F = dF[i,:,:]
dE[i,:,:] = (np.dot(d_F.T,F) + np.dot(F.T,d_F))*0.5
d_E = dE[i,:,:]
dP[i,:,:] = np.dot(d_F,2 * mu * E + la * (E[0,0] + E[1,1]) * np.identity(2))
dP[i,:,:] += np.dot(F,2 * mu * d_E + la * (d_E[0,0] + d_E[1,1]) * np.identity(2))
dH[i,:,:] = - area * np.dot(dP[i,:,:],minv.T)继续回想我们的三角形，最右边那个顶点，如果在x轴的坐标减一，那么Hessian矩阵中所有数值将减去dHdx2，也就是将加上一个正数。而Hessian矩阵长这个样子

$\bold H = - area \bold P(\bold F)\bold D_m^{-T} = \begin{bmatrix} -9 & 0 \\ 0 & -3/2\end{bmatrix} \tag{40}$
所以Hessian矩阵将继续接近零矩阵。说明求导的方向是对的。而如果最上面那个顶点在x轴移动，也就是增加减少dHdx4，那么并无助于把hessian矩阵变成零矩阵，这么一想也是合理的，因为上面的顶点在x移动时并不改变三角形的面积。而如果在y轴减去一个正数，也就是让原hessian矩阵减去dHdx5，那么原hessian矩阵也会更加接近零矩阵或单位矩阵，这样也是合理的，因为这样缩小了面积，让面积更接近初始的面积。
现在需要将这些3x2x2x2的四阶张量填到矩阵K中，仍然每次一列，每次使用一个2x2矩阵即可。
不过你得想清楚究竟是把什么东西填到什么东西上去。因为二维三角形三个顶点，每个顶点相互作用就是九种关系，也就是每次处理一个三角形元素的时候，我们需要添加把3x3个2x2个上矩阵放到K矩阵上去。因此我们需要9个2x2矩阵。但是我们的dH 却只有6个2x2矩阵，因此我们需要用动量守恒的关系式，请看下面的代码
第一列相当于这么填的

$\bold K = \begin{bmatrix} -(\partial \bold H/\partial x_0)_{00}-(\partial \bold H/\partial x_0)_{01} & * & * & * & * & * \\ -(\partial \bold H/\partial x_0)_{10}-(\partial \bold H/\partial x_0)_{11} & * & * & * & * & * \\ (\partial \bold H/\partial x_0)_{00} & * & * & * & * & * \\ (\partial \bold H/\partial x_0)_{10} & * & * & * & * & * \\ (\partial \bold H/\partial x_0)_{01} & * & * & * & * & * \\ (\partial \bold H/\partial x_0)_{11} & * & * & * & * & *\end{bmatrix} \tag{41}$
写成代码如下
# 3 个顶点
for n in range(3):
   # 2 个维度
   for d in range(2):
      # 第 idx 列，每列3 x 2 个元素
      idx = n * 2 + d
      # 先填写第一第二个顶点，第零个顶点之后填
      K[2,idx] += dH[idx,0,0]
      K[3,idx] += dH[idx,1,0]
      K[4,idx] += dH[idx,0,1]
      K[5,idx] += dH[idx,1,1]

      K[0,idx] += - dH[idx,0,0] - dH[idx,0,1]
      K[1,idx] += - dH[idx,1,0] - dH[idx,1,1]最后隐式一阶欧拉积分，再把公式写一遍

$(\bold M - \Delta t^2 \bold K|_{\bold x^t})\bold v^{t+1} = \bold M \bold v^t + \Delta t F(\bold x^t) \tag{42}$
积分代码如下
A = mass * np.identity(6) -  K  * dt * dt
b = np.zeros((6))
for n in range(3):
      for d in range(2):
         b[n*2+d] = mass * node_vel[n,d] + dt * node_force[n,d]

x = np.dot(np.linalg.inv(A), b)
for n in range(3):
      for d in range(2):
         node_vel[n,d] = x[n*2+d]
         node_pos[n,d] += node_vel[n,d]*dt接下来是三维隐式，但是和二维差不多。但在此之前，先让我们看看别的开源库是怎么写的
开源代码大赏

比如上篇就介绍了的比如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#L781，也计算了用于隐式计算的dF, dE, dP之类的。
case MATERIAL_TYPE_StVK:
   {
      EigenMatrix3 E = 0.5 * (F.transpose()*F - EigenMatrix3::Identity());
      EigenMatrix3 Finv = F.inverse();
      EigenMatrix3 FinvT = Finv.transpose();

      for (unsigned int i = 0; i < 3; i++)
      {
         for (unsigned int j = 0; j < 3; j++)
         {
               deltaF = dFdF(i, j);
               deltaE = 0.5 * (deltaF.transpose() * F + F.transpose() * deltaF);
               dPdF(i, j) = deltaF * (2 * m_mu * E + m_lambda * E.trace() * EigenMatrix3::Identity()) + F * (2 * m_mu * deltaE + m_lambda*deltaE.trace()*EigenMatrix3::Identity());
               ...同样也是上篇介绍的耶鲁大学图形组助理教授Theodore Kim的开源库Cubica库，项目地址http://www.tkim.graphics/cubica/。在stvk.cpp中，隐式计算的结点对结点求导是这样的，
MATRIX STVK::stiffness(TET& tet, bool diagonal)
{
computeStresses(tet);
computePFPu(tet, _pF_pu);
MATRIX product = _pf_pF * _pF_pu;
return product;
}f 就是力，p 是 partial ，为偏导数的意思。F是变形梯度。u是变形梯度。所以上面的代码实际上就是

$\frac{\partial \bold f}{\partial \bold x} = \frac{\partial \bold f}{\partial \bold F}\frac{\partia \bold F}{\partial \bold x} \tag{43}$
同样cubica库利用maple自动生成c++代码，展开了全部的矩阵运算，导致stvk.cpp有7000多行代码。虽然大部分代码都是自动生成的，但也足可见隐式积分的计算量有多大了。你没有看错

另一个全部展开矩阵计算的库是物理模拟库是bartels，地址在https://github.com/dilevin/Bartels/tree/master/src。
又比如超级无敌的vegafem库中的isotropicHyperelasticFEM.cpp，隐式积分中的矩阵K是这么计算。地址在https://github.com/starseeker/VegaFEM/blob/master/libraries/isotropicHyperelasticFEM/isotropicHyperelasticFEM.cpp#L734。
void IsotropicHyperelasticFEM::ComputeTetK(int el, double K[144], int clamped)
{
/*
   dP/dF is a column major matrix, but is stored as a 1D vector

   | dP_11/dF_11  dP_11/dF_12  dP_11/dF_13  dP_11/dF_21 ... dP_11/dF_33 |
   | dP_12/dF_11  dP_12/dF_12  dP_12/dF_13  dP_12/dF_21 ... dP_12/dF_33 |
   |                            ...                                  |
   | dP_33/dF_11  dP_33/dF_12  dP_33/dF_13  dP_33/dF_21 ... dP_33/dF_33 |
*/
double dPdF[81]; //in 9x9 matrix format
double dGdF[81]; //in 9x9 matrix format

Compute_dPdF(el, dPdF, clamped);
Compute_dGdF(&(areaWeightedVertexNormals[4 * el + 0]), &(areaWeightedVertexNormals[4 * el + 1]),
            &(areaWeightedVertexNormals[4 * el + 2]), dPdF, dGdF);
//dF_dU was already computed by the constructor before calling this function
double * dFdU = &dFdUs[108 * el];

// K is stored column-major (however, it doesn&#39;t matter because K is symmetric)
for (int row=0; row<9; row++)
{
   for (int column=0; column<12; column++)
   {
   double result = 0;
   for (int inner=0; inner<9; inner++)
   {
   //dGdF is 9x9, and dFdU is 9x12
   result += dGdF[9 * row + inner]*dFdU[12 * inner + column];
   }
   K[12 * column + row] = result;
   }
}

//The last three columns are combinations of the first nine columns.
//The reason is that the nodal force of the 4th vertex equals to
//the minus of the sum of the 1st, 2nd, and 3rd vertices (see p3
//section 4 of [Irving 04]
for (int row = 0; row < 12; row++)
{
   //10th column
   K[12 * row +  9] = -K[12 * row + 0] - K[12 * row + 3] - K[12 * row + 6];
   //11th column
   K[12 * row + 10] = -K[12 * row + 1] - K[12 * row + 4] - K[12 * row + 7];
   //12th column
   K[12 * row + 11] = -K[12 * row + 2] - K[12 * row + 5] - K[12 * row + 8];
}
}
vegafem在注释给出的公式是

$\bold K = \frac{\partial \bold G}{\partial \bold u} = \frac{\partial \bold G}{\partial \bold u}\frac{\partial \bold G}{\partial \bold u} = (\frac{\partial \bold P}{\partial \bold F}\bold B_m) \frac{\partial \bold F}{\partial \bold u} \tag{44}$
其中G是节点力。然而这公式在不同的地方也长得不一样。在上篇介绍过的论文“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{45}$
其中AN是当前元素，即三角形的法向量根据面积的权重。并且节点0的力是节点123的力相加并乘以-1。另一种方法是

$\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{46}$
在“FEM Simulation of 3D Deformable Solids: A practitioner’s guide to theory, discretization and model reduction.” 是这么写的

$\bold H = \begin{bmatrix} \vec f_1 & \vec f_2 & \vec f_3\end{bmatrix} = -W\bold P(\bold F)\bold D_m^{-T} \qquad \vec f_0 = -\vec f_1 - \vec f_2 - \vec f_3 \tag{47}$
所以G和H都是一样的东西，唯一的区别就是符号不同。【我一直觉得英文只有26个字母真是太少了，表达数学公式时完全不够用，加上希腊字母也不多】
伯克利大学的开源库arcsim库，隐式积分是这么写的，地址在https://github.com/zhou13/arcsim/blob/master/src/physics.cpp#L542
   // Expand Mx&#39;&#39; + Cx&#39; + Kx = f using backward Euler we have
   // 1/Dt [M + Dt C + Dt^2 (K - J)] Dv = F + (Dt J - C - Dt K) x&#39; - K x
   // where x&#39;&#39;(t + Dt) = (x&#39;(t + Dt) - x&#39;(t)) / Dt = Dv / Dt
   // For first step we have
   // M Dv/Dt = F (x + Dt (v + Dv))
   // A Dv = b
   // A = M - Dt^2 DF/Dx
   // b = Dt (F(x) + Dt DF/Dx v))
   SpMat<Mat3x3> A(nn, nn);
   vector<Vec3> b(nn, Vec3(0));

   for (size_t n = 0; n < nn; n++) {
      A(n, n) += Mat3x3(nodes[n]->m) - dt * dt * Jext[n];
      b[n] += dt * (fext[n] + dt * Jext[n] * nodes[n]->v);
   }
   ...
   vector<Vec3> dv = mat3x3_linear_solve(A, b);然后是泽森科工的ZENO开源库，隐式一阶欧拉积分为在https://github.com/zenustech/zeno/blob/master/projects/FEM/src/integrator/backward_euler_integrator.cpp，代码如下
int BckEulerIntegrator::EvalElmObjDerivJacobi(const TetAttributes attrs,
      const std::shared_ptr<MuscleForceModel>& force_model,
      const std::shared_ptr<DiricletDampingModel>& damping_model,
      const std::vector<Vec12d>& elm_states,
      FEM_Scaler* elm_obj,Vec12d& elm_deriv,Mat12x12d& elm_H,bool filtering) const{
   Vec12d u2 = elm_states[2];
   Vec12d u1 = elm_states[1];
   Vec12d u0 = elm_states[0];
   FEM_Scaler h = _dt;
   Vec12d v2 = (u2 - u1)/h;

   FEM_Scaler PhiI,PsiE,PsiD;
   Vec9d dPsiE,dPsiD;
   Mat9x9d ddPsiE,ddPsiD;

   FEM_Scaler vol = attrs._volume;
   FEM_Scaler m = vol * attrs._density / 4;  // nodal mass

   Mat3x3d F,L;
   ComputeDeformationGradient(attrs._Minv,u2,F);
   ComputeDeformationGradient(attrs._Minv,v2,L);
   force_model->ComputePhiDerivHessian(attrs._activation,
      attrs._fiberWeight,
      attrs._fiberOrient,
      attrs._E,
      attrs._nu,
      F,
      PsiE,dPsiE,ddPsiE,filtering);
   damping_model->ComputePhiDerivHessian(attrs._d,L,PsiD,dPsiD,ddPsiD);

   Vec12d y = u2 - 2*u1 + u0 - h*h*_gravity.replicate(4,1);
   PhiI = y.squaredNorm() * m / 2 / h;
   *elm_obj = PhiI + h * PsiE * vol + h * h * PsiD * vol;

   const Mat9x12d& dFdX = attrs._dFdX;
   elm_deriv = m * y / h + dFdX.transpose() * (h*dPsiE + h*dPsiD) * vol;

   elm_H = Mat12x12d::Identity() * m / h + vol * dFdX.transpose() * (h*ddPsiE + ddPsiD) * dFdX;

   return 0;
}来猜猜看，上面所给出代码，每行代码分别对应了哪个公式呢？
三维隐式

不管怎么说，先来看看我们是如何推导Ds对顶点求导，与二维一样。

$\bold x_0 = \begin{bmatrix} x_0 \\ x_1 \\x_2\end{bmatrix} \qquad \bold x_1 = \begin{bmatrix} x_3 \\ x_4 \\x_5\end{bmatrix} \qquad \bold x_2 = \begin{bmatrix} x_6\\ x_7 \\x_8\end{bmatrix} \qquad \bold x_3 = \begin{bmatrix} x_9 \\ x_{10} \\x_{11}\end{bmatrix} \tag{48}$
那么Ds计算方式如下

$\bold Ds = \begin{bmatrix} x_3 - x_0 & x_6 - x_0 & x_9 - x_0 \\ x_4 - x_1 & x_7 - x_1 & x_{10} - x_1 \\ x_5 - x_2 & x_8 - x_2 & x_{11} - x_2\end{bmatrix} \tag{48}$
对每个顶点求导，得到的dD是4x3x3x3的四阶向量。注意我的顶点顺序与“FEM Simulation of 3D Deformable Solids: A practitioner’s guide to theory, discretization and model reduction.  ”这个讲义力的并不一样

$\frac{\partial \bold D_s}{\partial x_0} = \begin{bmatrix} -1 & - 1 & - 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix} \quad \frac{\partial \bold D_s}{\partial x_1} = \begin{bmatrix} 0 & 0 & 0 \\-1 & - 1 & - 1 \\ 0 & 0 & 0\end{bmatrix} \quad \frac{\partial \bold D_s}{\partial x_2} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\-1 & - 1 & - 1 \end{bmatrix} \\ \frac{\partial \bold D_s}{\partial x_3} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix} \quad \frac{\partial \bold D_s}{\partial x_4} = \begin{bmatrix} 0 & 0 & 0 \\1 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix} \quad \frac{\partial \bold D_s}{\partial x_5} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\1 &0 & 0 \end{bmatrix} \\ \frac{\partial \bold D_s}{\partial x_6} = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix} \quad \frac{\partial \bold D_s}{\partial x_7} = \begin{bmatrix} 0 & 0 & 0 \\0 & 1 & 0 \\ 0 & 0 & 0\end{bmatrix} \quad \frac{\partial \bold D_s}{\partial x_8} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\0 & 1 & 0 \end{bmatrix} \\ \frac{\partial \bold D_s}{\partial x_9} = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix} \quad \frac{\partial \bold D_s}{\partial x_{10}} = \begin{bmatrix} 0 & 0 & 0 \\0 & 0 & 1 \\ 0 & 0 & 0\end{bmatrix} \quad \frac{\partial \bold D_s}{\partial x_{11}} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\0 & 0 & 1 \end{bmatrix} \\$
接下来和二维几乎是一样的。同样我们也可把dFdx换成一个9x12个矩阵，而不是四阶张量。就像“Dynamic Deformables: Implementation and Production Practicalities  ”所附代码那样，把矩阵全部展开。cubica库也将矩阵运算全部展开了。首先是Dminv求逆

$\bold d = \begin{bmatrix} d_{00} & d_{01} & d_{02} \\ d_{10} & d_{11} & d_{12} \\ d_{20} & d_{21} & d_{22}\end{bmatrix} = \bold D_m^{-1} \tag{50}$
那么pFpx就是下面的矩阵。也有相当多的开源库是这么算的。

$\frac{\partial \bold F}{\partial \bold x} = \begin{bmatrix} -d_{00} - d_{01} - d_{02} & 0 & 0 & d_{00} & 0 & 0 &d_{01} & 0 & 0 & d_{02} & 0 & 0 \\ 0 & -d_{00} - d_{01} - d_{02} & 0 & 0 & d_{00} & 0 & 0 &d_{01} & 0 & 0 & d_{02} & 0 %C\ 0 & 0 & -d_{00} - d_{01} - d_{02} & 0 & 0 & d_{00} & 0 & 0 &d_{01} & 0 & 0 & d_{02} \\ d_{10} - d_{11} - d_{12} & 0 & 0 & d_{10} & 0 & 0 &d_{11} & 0 & 0 & d_{12} & 0 & 0 \\ 0 & -d_{10} - d_{11} - d_{12} & 0 & 0 & d_{10} & 0 & 0 &d_{11} & 0 & 0 & d_{12} & 0 \\ 0 & 0 & -d_{10} - d_{11} - d_{12} & 0 & 0 & d_{10} & 0 & 0 &d_{11} & 0 & 0 & d_{12} \\ d_{20} - d_{21} - d_{22} & 0 & 0 & d_{20} & 0 & 0 &d_{21} & 0 & 0 & d_{22} & 0 & 0 \\ 0 & -d_{20} - d_{21} - d_{22} & 0 & 0 & d_{20} & 0 & 0 &d_{21} & 0 & 0 & d_{22} & 0 \\ 0 & 0 & -d_{20} - d_{21} - d_{22} & 0 & 0 & d_{20} & 0 & 0 &d_{21} & 0 & 0 & d_{22} \end{bmatrix}$
接下来就是二维的完全一样了，至少对于stvk模型来说是这样的。python完整代码在
https://github.com/clatterrr/PhysicsSimulationPlayground/blob/master/ImplictElasticFem/StvkImplicitMultiple3D.py
https://gitee.com/clatterrr/PhysicsSimulationPlayground/blob/master/ImplictElasticFem/StvkImplicit3D.py
Unity 在computerShader实现的坑

首先为了生成四面体模型，你需要一个tetgen，下载到本地并编译运行没什么说的。然后在blender弄好模型，导出为带ascii的ply模型，否则tetgen不认。tetgen处理好后，生成的.node和.ele和.face文件是我们需要的。
其中.node就是顶点的位置信息，而.ele就是四面体索引信息，而.face则是表面三角形的索引信息，不参与数值计算，仅仅用于绘制。但是我们并没有uv信息用于绘制，这就很尴尬了qwq。
在并行计算的时候，在算力的时候，由于是逐四面体计算，肯定有不同的四面体共用同一个顶点，因此需要将顶点力加到同一个顶点上。因此需要原子操作。然而内置原子操作不支持浮点，因此我们可以自己实现一个原子相加浮点的函数。如下
RWByteAddressBuffer  node_force;
...
void atomicAddFloat(RWByteAddressBuffer buf, int idx, float value)
{
   uint i_val = asuint(value);
   uint tmp0 = 0;
   uint tmp1;
   // Raw Bmatyte 4 点浮点数对齐
   uint meidx = idx * 4;
   while (true)
   {
      //如果node_force[idx] == tmp0，
      //那么tmp1 = node_force[idx], node_force[idx] = i_val
      //此时tmp0 == tmp1 == 原始值
      //因为有可能在相加的过程中，别的核也在加，那么自己就被改了
      //所以除非相等，才保证在自己改动的过程别人没改过
      //如果不相等，就要更新别人改过后的值，也就是tmp1
      buf.InterlockedCompareExchange(meidx, tmp0, i_val, tmp1);
      if (tmp0 == tmp1)
         break;
      tmp0 = tmp1;
      i_val = asuint(value + asfloat(tmp1));
   }
}
...
atomicAddFloat(node_force, tidx + 0, gradC0.x);第二个是线性方程组求解，我用的是稠密的共轭梯度。这是一个性能瓶颈，也许改成稀疏的共轭梯度会好很多。并且共轭梯度每次循环都在把两个向量点积成一个数，因此需要用到使用共享内存的parallelReduction。此外为了效率应当尽量展开循环。
groupshared float sharedMem[THREADTOTAL];
void ParallelReduction(uint3 Gid : SV_GroupID, uint3 DTid : SV_DispatchThreadID, uint GI : SV_GroupIndex)
{
   sharedMem[GI] = Source[DTid.x]; // store in shared memory
   GroupMemoryBarrierWithGroupSync(); // wait until everything is transfered from device memory to shared memory

   for (uint s = THREADTOTAL / 2; s > 0; s >>= 1) {
      if (GI < s)
         sharedMem[GI] += sharedMem[GI + s];
      GroupMemoryBarrierWithGroupSync();
   }
   // Have the first thread write out to the output
   if (GI == 0) {
      Results[Gid.x] = sharedMem[0];
   }
}然而有关computeShader的计算技巧的资料并不多，最好的办法是先去学NVIDA的CUDA。这样就能理解很多并行技巧了。另外共轭梯度需要判断什么时候停止迭代，但是如果把数据传回cpu判断，那么帧率就会严重下降。因此我干脆不判断了。也许有更好的方法。
然后是添加碰撞力和重力，或者鼠标拖拽力的时候，一般都是在求解方程前把力加到node_force上。如果是如果求解方程后再手动修改位置和速度，可能容易穿模。
到目前为止，我们实现了三维隐式stvk模型。但除了stvk模型外，还有很多经典的模型，比如常用于肌肉模拟的neohookean模型，AsRigidAsPossible模型，Corotation模型等。除此之外还有很多技巧来优化弹性的表现。这就是接下来的事情的。
不过在此之前，需要首先搞定另一个可怕的东西，三维数值奇异值分解。
上一篇：

NoiseFloor · 发表于 2021-12-20 09:28

就算是隐式积分有时候fem都会爆掉，有时候还要搞点trick fem才能稳定住[飙泪笑] games201上貌似说了隐式积分结合linear search 只能说优化的坑太深了

yukamu · 发表于 2021-12-20 09:30

这图笑死我了[捂脸]

fwalker · 发表于 2021-12-20 09:33

实时的话可以看一下xpbd或者projective dynamics

LiteralliJeff · 发表于 2021-12-20 09:38

讲道理，想求下无水印的图[飙泪笑]

Mecanim · 发表于 2021-12-20 09:47

嗯，我去试试

七彩极 · 发表于 2021-12-20 09:48

看来隐式积分也不是万能的哈哈

RecursiveFrog · 发表于 2021-12-20 09:51

水印取掉了，自取~

Ilingis · 发表于 2021-12-20 09:53

感谢，给你点个赞

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

[笔记] [游戏物理探秘]三维有限元的隐式积分方法

本帖子中包含更多资源