第五講：彈性體模擬-數(shù)值方法

每個(gè)人都身懷天賦掰曾，但如果用會(huì)不會(huì)爬樹來衡量一條魚璃哟，那它到死都只會(huì)認(rèn)為自己愚蠢。 —— 愛因斯坦

這一講是緊接著上一講的內(nèi)容肄方，開始討論如何通過數(shù)值方法來完成對(duì)彈性體形變的模擬，另外蹬癌，順帶提一句权她，這部分內(nèi)容跟游戲開發(fā)中的物理引擎具有較高關(guān)聯(lián)。

上一講降到逝薪，我們的彈性勢(shì)能可以表示為：
$P(\epsilon) = \sum_ {i, j, k, l = 1} ^ 3 C_{ijkl} \epsilon_{ij} \epsilon_{kl}$

值得澄清的是隅要，這里是省略了高階項(xiàng)之后的彈性勢(shì)能，即這個(gè)公式是在線性形變的情況下成立的董济，如果要考慮更為復(fù)雜的形變步清，后面還要加上更多的高階項(xiàng)。

下面來介紹，我們?cè)谟?jì)算機(jī)中要如何計(jì)算出彈性體形變之后產(chǎn)生的internal force廓啊，大致思路是先計(jì)算出形變梯度矩陣（Deformation Gradient）欢搜，之后根據(jù)這個(gè)矩陣算出Green Strain Tensor，之后再算出勢(shì)能表達(dá)式谴轮，最后對(duì)勢(shì)能求導(dǎo)炒瘟，就得到了對(duì)應(yīng)的conservative force。

首先來看下第步，我們要如何計(jì)算形變梯度矩陣疮装，這里先考慮3D軟體的情況，2D軟體的情況如布料等在后面再考慮粘都。

首先斩个，要想在計(jì)算機(jī)中進(jìn)行模擬，就需要將3D軟體離散化驯杜，2D表面離散化是通過將之拆解成三角形來完成受啥，而3D Volume的離散則是將之拆分成一個(gè)個(gè)的四面體。

形變發(fā)生后鸽心，單個(gè)四面體就會(huì)從material space映射到deformed space滚局，即四面體本身發(fā)生了形變，如下圖所示：

緊接著給出兩個(gè)假設(shè)：

四面體足夠小顽频，那么假設(shè) $\phi(X_i) = x_i$ 藤肢，我們可以進(jìn)行如下的推導(dǎo)：

$x_1 = \phi(X_0 + X_1 - X_0) = \phi(X_0 + X_{10})$

由于四面體足夠小，因此 $X_{10}$ 滿足微分中的一個(gè)極小的局部變化的要求糯景，可以采用泰勒展開：

$x_1 = \phi(X_0) + \frac{\partial \phi}{\partial X} X_{10} +... = x_0 + F X_{10} + ...$

在線性形變的情況下嘁圈，我們可以省略高階項(xiàng)，我們可以得到：

$x_1 - x_0 = F X_{10}$

上述公式中的 $x_1, x_0$ 我們是知道的（最開始的時(shí)候蟀淮， $x_i = X_i$ 最住，所以自然是知道的，之后怠惶，經(jīng)過一個(gè)timestep迭代涨缚，最新的位置經(jīng)過后續(xù)的計(jì)算，自然也是知道的）策治，而形變前的數(shù)據(jù) $X_{10}$ 我們也知道脓魏，因此就可以計(jì)算出F，另外通惫，我們前面也說過茂翔，F(xiàn)是一個(gè)與位置有關(guān)的變量，要精確考慮的話履腋，四面體上的每一點(diǎn)對(duì)應(yīng)的F都應(yīng)該是不同的珊燎，不過這里我們可以假設(shè)：

四面體內(nèi)部的F是一個(gè)常量，那么F就變成了一個(gè)piece-wise constant（可以聯(lián)想分段函數(shù)），當(dāng)然俐末，這里我們其實(shí)還隱藏了一個(gè)假設(shè)料按，那就是從 $X_1$ 到 $X_0$ 的變化是一個(gè)線性變化（想象一下兩個(gè)頂點(diǎn)之間有一個(gè)彈簧），在這種假設(shè)下卓箫，F(xiàn)才會(huì)是一個(gè)常量载矿，否則可能會(huì)是一個(gè)多項(xiàng)式，目前市面上的物理引擎基本上都是按照線性變化來模擬的烹卒，雖然精度上不能做到完全保證闷盔，但是考慮到性能跟效果的平衡，這是一種最好的選擇了旅急。

在上述假設(shè)下逢勾，我們可以有更多的等式：

$x_2 - x_0 = F X_{20}$

$x_3 - x_0 = F X_{30}$

將前面的三組公式放在一起，我們就得到了下面的矩陣乘法形式藐吮，從而推導(dǎo)出F：

$[x_1 - x_0, x_2 - x_0, x_3 - x_0] = F [X_{10}, X_{20}, X_{30}]$

簡(jiǎn)化一下：

$[x] = F [X]$

$[x_1 - x_0, x_2 - x_0, x_3 - x_0]$ 跟 $[X_{10}, X_{20}, X_{30}]$ 都是一個(gè)3x3的矩陣溺拱，根據(jù)上面公式可以推導(dǎo)出：

$F(x) = [x][X]^{-1}$

這里是其中一個(gè)四面體的F，其他四面體也會(huì)有對(duì)應(yīng)的F谣辞，根據(jù)F我們就可以計(jì)算Green Strain Tensor迫摔，這也是每個(gè)四面體一份的：

$\epsilon = F^TF -I$

以及勢(shì)能 $P(\epsilon(x))$ ，對(duì)勢(shì)能進(jìn)行求導(dǎo)（具體來說泥从，consititutive model有公式的話句占，可以手動(dòng)推導(dǎo)其彈性力的表達(dá)式），我們就能得到internal force（彈性力）：

$f_{int} = \frac{\partial P}{\partial x}$

這里躯嫉，勢(shì)能是物體共有的一個(gè)參數(shù)纱烘，是一個(gè)標(biāo)量，而 $x$ 則是物體上所有四面體的頂點(diǎn)組成的一個(gè)向量 $[x_1, x_2, ... x_n]$ 祈餐，那么經(jīng)過求導(dǎo)后得到的彈性力也自然是一個(gè)n維的向量擂啥，每一維代表的是每個(gè)頂點(diǎn)的受力。

上面是內(nèi)部的力昼弟，外面來自于重力或者碰撞的力啤它，我們用外力 $f_{ext}$ 表示。

$f = f_{int} + f_{ext} = m a$

$a = [x_1, x_2, ... x_n]''$

其中的質(zhì)量m是一個(gè)很大的vector舱痘，在上面的公式下，我們就可以用之前說過的implicit euler或者explicit euler完成時(shí)間步長(zhǎng)積分离赫，從而算出更新的位置向量芭逝。

上面這種用離散的四面體來模擬軟體的彈性形變的方法叫做Finite Element Method（FEM），這種方法可以通過不斷提高分辨率（縮小四面體的尺寸）來保證準(zhǔn)確性渊胸，不過同時(shí)也會(huì)導(dǎo)致計(jì)算復(fù)雜度增加旬盯，所以游戲或者很多電影特效模擬不是用的這種方法，而是基于優(yōu)化的方法，Optimization-Based Methods胖翰，這種方法可以實(shí)現(xiàn)快速模擬接剩。

先來回顧一下隱式歐拉方法，速度跟位置向量的演變可以用下面公式來表達(dá)：

$q_{n+1} = q_n +h v_{n+1}$

$v_{n+1} = v_n +h M^{-1}(f_{int}(q_{n+1}) + f_{ext})$

這里的兩個(gè)向量都是包含了多個(gè)四面體的相關(guān)數(shù)據(jù)的萨咳，是一個(gè)3xn的矩陣。

前面說過隱式歐拉最終會(huì)變成一個(gè)高維方程，很難求解张弛，這里考慮將之轉(zhuǎn)換成一個(gè)優(yōu)化問題草添。

首先，internal force是勢(shì)能的導(dǎo)數(shù)給出如下：

$f_{int}(q_{n+1}) = - \nabla_q W$

根據(jù)上面的公式舀凛，我們有：

$q_{n+1} = q_n +h v_{n}+h^2 M^{-1}(f_{int}(q_{n+1}) + f_{ext})$

這邊希望構(gòu)造出一個(gè)函數(shù) $J(q_{n+1})$ 俊扳，并且使得這個(gè)函數(shù)的極值：

$min J(q_{n+1}) \rightarrow \nabla J =0$

達(dá)成的條件跟前面的公式相一致，也就是說猛遍，希望這個(gè)函數(shù)的極值對(duì)應(yīng)的解 $q_{n+1}$ 就正好是上面一個(gè)等式的解馋记，而這個(gè)函數(shù)可以給出為（比較容易推導(dǎo)？）：

$J = \frac{1}{2 h^2}||M^{\frac{1}{2}}(q_{n+1} - S_n)||_F^2 + W(q_{n+1})$

其中

$|| A ||_F = \sum_{i,j} A_{ij} ^ 2$

$S_n = q_n + h v_n + h^2 M^{-1} f_{ext}$

從而我們可以知道：
$q_{n+1} - S_n = h^2 M^{-1} f_{int}(q_{n+1})$

前面說到懊烤，這個(gè)函數(shù)提出的目的是為了求得前面implicit timestep等式的等價(jià)解抗果，這里的一個(gè)疑惑就是，同樣是解方程奸晴，前面implicit方程不好解冤馏，這里的方程同樣也不好解，那么這種變化的意義在哪里呢寄啼？

這里的作用在于公式中的后面部分 $W(q_{n+1})$ 逮光，這個(gè)部分表示的是系統(tǒng)的勢(shì)能，前面說過彈性形變的勢(shì)能公式 $P(\epsilon)$ 是通過一系列的物理推導(dǎo)得出的墩划，是精確解涕刚，這種求解起來會(huì)比較麻煩，但是在游戲中乙帮，我們可以使用非精確解杜漠，只要保證模擬的結(jié)果看起來真實(shí)就好，這樣一來察净，我們就能夠選擇一些能夠讓這個(gè)求解過程變得容易的公式代入這里的勢(shì)能項(xiàng)中驾茴，從而達(dá)到簡(jiǎn)化問題的目的。

那么這個(gè)勢(shì)能公式要怎么來設(shè)計(jì)呢氢卡？我們知道锈至，彈性形變導(dǎo)致的勢(shì)能總體來說，可以分成兩種译秦，分別是體積變化導(dǎo)致的形變（如從體積為1變成體積為2）峡捡，或者是扭曲導(dǎo)致的形變击碗，如一個(gè)cube變成平行六面體等。

前者勢(shì)能可以表示為：
$(det(F) -1)^2$

det是矩陣的行列式们拙，這個(gè)公式的意義在于稍途，當(dāng)沒有形變的時(shí)候，矩陣就是單位矩陣砚婆，這個(gè)勢(shì)能的結(jié)果就是0

后者勢(shì)能可以表示為：

$(tr(F^T F) -3 ) ^2$

tr（trace）是矩陣的跡械拍，一個(gè)n×n矩陣A的主對(duì)角線（從左上方至右下方的對(duì)角線）上各個(gè)元素的總和被稱矩陣A的跡（也就等于所有特征值的和）。同樣射沟，當(dāng)沒有形變時(shí)殊者，這個(gè)F就是單位矩陣，所以跡就是3验夯，那么勢(shì)能依然為0猖吴。

總結(jié)可以看出，我們需要的勢(shì)能公式其實(shí)就是一個(gè)約束條件挥转，這個(gè)約束條件或者說約束函數(shù)在無形變的時(shí)候輸出為0海蔽，在有形變的時(shí)候輸出的是一個(gè)正值（或者說在無形變的時(shí)候取得極小值0），比如其形式為：

$C(q) = 0$

在這個(gè)條件下绑谣，我們就可以定義勢(shì)能為：

$W = \frac{1}{2} C^T(q) A C(q)$

上式中的A是一個(gè)對(duì)角陣党窜，C(q)是一個(gè)向量，對(duì)于其中的每個(gè)元素：

$C_i = det(F_i) - 1$

這里希望有：

$C_i = 0$

這里的i用來指示四面體的序號(hào)借宵，這里翻譯一下幌衣，就是說我們希望每個(gè)四面體的形變梯度矩陣的行列式為1，也就是說每個(gè)四面體都是體積不變的壤玫。

下面來介紹一下Position-based Dynamics方法豁护，回到前面的公式：

$J = \frac{1}{2 h^2}||M^{\frac{1}{2}}(q_{n+1} - S_n)||_F^2 + W(q_{n+1})$

這里，我們將勢(shì)能用脈沖函數(shù)（實(shí)際上是脈沖函數(shù)取反欲间，不過為了描述方便楚里，這里用脈沖函數(shù)來表示）來表示：

$J = \frac{1}{2 h^2}||M^{\frac{1}{2}}(q_{n+1} - S_n)||_F^2 + \sum_i \delta_{C_i}(q_{n+1})$

其中脈沖函數(shù)：

$\delta_{C_i}(q_{n+1})= \begin{cases} +\infty & {C_i(q_{n+1}) \neq 0}\\ 0 & {C_i(q_{n+1}) = 0} \end{cases}$

之所以這樣設(shè)計(jì)，就是為了滿足前面說過的猎贴，我們希望函數(shù)J取得最小值班缎，那么就需要保證 $C_i(q_{n+1}) = 0$ ，這里的i依然是四面體的序號(hào)她渴。

繼續(xù)下去达址，Position-based Dynamics方法會(huì)采用Gauss-Siedel Iteration（假設(shè)1）進(jìn)行求解，即對(duì)于后面的眾多四面體而言惹骂，我們會(huì)先只考慮一個(gè)四面體作用下的求解苏携，并將結(jié)果代入原公式之后考慮第二個(gè)四面體作用下的求解，循環(huán)往復(fù)对粪。如果只考慮單個(gè)四面體作用右冻，原公式就變成了：

$J = \frac{1}{2 h^2}||M^{\frac{1}{2}}(q_{n+1} - S_n)||_F^2 + \delta_{C_i}(q_{n+1})$

再定義：
$\Delta q = q_{n+1} - S_n$

公式就變成：

$J = \frac{1}{2 h^2}||M^{\frac{1}{2}}\Delta q||_F^2 + \delta_{C_i}(q_{n+1})$

將公式轉(zhuǎn)換一下，我們想要求解J的最小值著拭，就等價(jià)于在條件：

$C_i(S_n + \Delta q) = C_i(q_{n+1}) = 0$

下求解

$\frac{1}{2 h^2}||M^{\frac{1}{2}}\Delta q||_F^2$

的最小值纱扭。

之后，再假設(shè) $\Delta q$ 足夠欣苷凇（假設(shè)2）乳蛾，那么前面的條件就可以做泰勒展開為：

$C_i(S_n + \Delta q) = C_i(S_n) + \nabla^T C_i \Delta q = 0$

而這個(gè)就可以用上一講的拉格朗日乘子來求解，也就是說鄙币，在上面條件下的最小值求解就變成了：

$\Delta q = - M^{-1} \nabla C(q) \frac{C(q)}{||M^{-\frac{1}{2}} \nabla C||_F^2}$

而Gauss-siedel Iteration的問題是肃叶，需要進(jìn)行多次迭代才能收斂，從而導(dǎo)致計(jì)算比較慢十嘿，性能較差因惭，這也是Position-based Dynamics（PBD）的缺點(diǎn)，如果Iteration數(shù)目不夠就會(huì)導(dǎo)致一些模擬效果質(zhì)量較低绩衷，比如布料會(huì)看起來像彈簧一樣生硬蹦魔。

我們來回顧下之前的內(nèi)容：

最開始我們是以隱式歐拉開頭，因?yàn)檎f到咳燕，想要求得下一步的位置勿决，需要求解一個(gè)十分復(fù)雜的非線性方程，不但計(jì)算量巨大招盲，而且由計(jì)算機(jī)來計(jì)算的話無法給出準(zhǔn)確解低缩；
于是，我們將問題轉(zhuǎn)換成了一個(gè)最優(yōu)化問題曹货，通過PBD+迭代的方式來求得近似解咆繁，但是這個(gè)方法的問題在于迭代時(shí)間較長(zhǎng)，很難快速得到結(jié)果控乾。

那么很自然的問題就是么介，是否有一種方法是居于兩者之間的，一方面我們不需要求解一個(gè)十分復(fù)雜的方程蜕衡，而是只需要求解一個(gè)簡(jiǎn)單的方程壤短，且這個(gè)方程可以通過少數(shù)的迭代就能輸出結(jié)果的。

這個(gè)方法就是物理界常用的Projective Dynamics慨仿，實(shí)際上這種方法跟PBD是處于同一框架下得久脯，不同的是PBD這里的internel energy函數(shù)是一個(gè)脈沖函數(shù)之和，而Projective Dynamics方法的internal energy函數(shù)則是：
$W(q_{n+1}) = d(p, q_{n+1}) + \delta_c(p)$

這個(gè)式子中一共引入了兩個(gè)新的符號(hào)镰吆，一個(gè)是變量p帘撰，這個(gè)變量代表的是某個(gè)粒子（如前面的四面體）的一個(gè)目標(biāo)狀態(tài)（比如，我們可以選擇四面體不發(fā)生形變時(shí)的狀態(tài)為一個(gè)目標(biāo)狀態(tài)万皿，當(dāng)然我們也可以選擇其他狀態(tài)作為目標(biāo)狀態(tài)摧找，具體如何選核行，后面再說），而d則是表達(dá)p跟q差異性的函數(shù)蹬耘，可以是一個(gè)非常簡(jiǎn)單的函數(shù)芝雪，比如就取兩者的距離作為函數(shù)。

整個(gè)函數(shù)是什么意思呢综苔，我們這里來舉個(gè)例子說明一下惩系，假如我們選擇某個(gè)函數(shù)（simple distance measure）作為上面的函數(shù)d：
$d(p, q_{n+1}) = || F(q_{n+1}) - R_0 ||_F^2 = || A \cdot q_{n+1} - R_0 ||_F^2$

這里的F跟前面一樣，還是deformation gradient如筛， $R_0$ 是一個(gè)沒有形變堡牡，但是經(jīng)過一個(gè)旋轉(zhuǎn)的狀態(tài)，兩者都是3x3的矩陣杨刨，其中 $R_0$ 可以通過一個(gè)叫做Shape Matching的算法求得晤柄。

另外一項(xiàng) $\delta_c(p)$ 是一個(gè)脈沖函數(shù)，這個(gè)脈沖函數(shù)的取值受到約束條件c的影響拭嫁，當(dāng) $c(p) = 0$ 的時(shí)候可免，取值為0，否則取值為無窮大做粤，而由于這個(gè)函數(shù)的存在浇借，要想使得整體的internal energy最小，就必須得滿足這里的約束條件怕品，從而使得脈沖函數(shù)取值為0妇垢。

在這個(gè)例子中，我們選取了p為經(jīng)過某個(gè)旋轉(zhuǎn) $R_0$ 后的狀態(tài)肉康，而在這個(gè)狀態(tài)下闯估，可以滿足上面的脈沖函數(shù)結(jié)果為0，那么最終的internal energy的計(jì)算就變成了一個(gè)非常簡(jiǎn)單的矩陣的線性計(jì)算吼和。

這里涨薪，我們還可以選擇其他的p，比如：
$d(p, q_{n+1}) = || F(q_{n+1}) - F_0 ||_F^2$

這里的 $F_0$ 也是一個(gè)deformation gradient炫乓，表示的是最接近于當(dāng)前形變狀態(tài)的deformation gradient刚夺，但是在這個(gè)gradient的作用下，形變后的四面體的體積跟原始四面體保持一致（這種對(duì)于動(dòng)物形體的模擬比較有用末捣，比如肚子在受擊之后體積保持不變）侠姑， $F_0$ 的計(jì)算同樣可以通過脈沖函數(shù)求得，比如只需要找到體積不變的形變的約束條件箩做，就能得到對(duì)應(yīng)的deformation gradient莽红。

在這種internal energy的作用下，最優(yōu)解的求解就變成了如下的函數(shù)形式：

$J = \frac{1}{2 h^2}||M^{\frac{1}{2}}(q_{n+1} - S_n)||_F^2 + d(p, q_{n+1}) + \delta_c(p)$

可以看到邦邦，p只出現(xiàn)在后兩項(xiàng)安吁，q則只出現(xiàn)在前兩項(xiàng)醉蚁，從而我們可以通過有限的迭代來求解這個(gè)最優(yōu)化問題。

而由于脈沖函數(shù)的存在柳畔，p的求解只需要保證脈沖函數(shù)為0即可馍管，因此整個(gè)問題可以快速變成q變量下的最優(yōu)解郭赐，這個(gè)過程叫做local projection薪韩，即找到限制條件c（constraint）來保證脈沖函數(shù)為0，從而將問題從兩個(gè)變量降維成一個(gè)變量捌锭。

在求得p之后俘陷，就可以來求解q以獲得整個(gè)函數(shù)的最優(yōu)解：

$J = \frac{1}{2 h^2}||M^{\frac{1}{2}}(q_{n+1} - S_n)||_F^2 + ||p - q_{n+1})||_2^2$

而由于上面的函數(shù)實(shí)際上都是q的一個(gè)二次函數(shù)，而其導(dǎo)數(shù)則是q的一次函數(shù)观谦，從而只需要將導(dǎo)數(shù)等于0拉盾，就能得到最終的解，這一步叫做global solve豁状。

前面說過捉偏，projective Dynamics是一個(gè)迭代的過程，上面展示的local projection跟global solve就是這個(gè)過程中的僅有兩個(gè)步驟泻红，這兩個(gè)步驟需要循環(huán)往復(fù)的去執(zhí)行（因?yàn)閜是local projection的夭禽，也就是一個(gè)局部解，這個(gè)解并不一定就對(duì)應(yīng)著最優(yōu)解谊路，而是需要不斷調(diào)整q來逼近最優(yōu)解）讹躯。

這個(gè)方法的迭代數(shù)目比PBD會(huì)少很多，目前UE缠劝、Houndini等軟件使用的是X-PBD方案潮梯，是在PBD方案上的一種更新方案，但是其迭代次數(shù)依然比Projective Dynamics方法要多惨恭，之所以不用Projective Dynamics方法秉馏，是因?yàn)榫€性系統(tǒng)的求解在維數(shù)很大（比如q是一個(gè)很大的向量）的時(shí)候會(huì)存在比較大的困難。

最后編輯于：2022.05.23 23:55:19

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者

人面猴
序言：七十年代末脱羡，一起剝皮案震驚了整個(gè)濱河市萝究，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌轻黑，老刑警劉巖糊肤，帶你破解...
沈念sama閱讀 212,816評(píng)論 6贊 492
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場(chǎng)離奇詭異氓鄙，居然都是意外死亡馆揉，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 90,729評(píng)論 3贊 385
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門抖拦，熙熙樓的掌柜王于貴愁眉苦臉地迎上來升酣，“玉大人舷暮，你說我怎么就攤上這事∝眩” “怎么了下面？”我有些...
開封第一講書人閱讀 158,300評(píng)論 0贊 348
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵，是天一觀的道長(zhǎng)绩聘。經(jīng)常有香客問我沥割，道長(zhǎng)，這世上最難降的妖魔是什么凿菩？我笑而不...
開封第一講書人閱讀 56,780評(píng)論 1贊 285
?港島之戀（遺憾婚禮）
正文為了忘掉前任机杜，我火速辦了婚禮，結(jié)果婚禮上衅谷，老公的妹妹穿的比我還像新娘椒拗。我一直安慰自己，他們只是感情好获黔，可當(dāng)我...
茶點(diǎn)故事閱讀 65,890評(píng)論 6贊 385
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布蚀苛。她就那樣靜靜地躺著，像睡著了一般玷氏。火紅的嫁衣襯著肌膚如雪堵未。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 50,084評(píng)論 1贊 291
城市分裂傳說
那天预茄，我揣著相機(jī)與錄音兴溜，去河邊找鬼。笑死耻陕，一個(gè)胖子當(dāng)著我的面吹牛拙徽，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播诗宣，決...
沈念sama閱讀 39,151評(píng)論 3贊 410
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼膘怕，長(zhǎng)吁一口氣：“原來是場(chǎng)噩夢(mèng)啊……” “哼！你這毒婦竟也來了召庞？” 一聲冷哼從身側(cè)響起岛心，我...
開封第一講書人閱讀 37,912評(píng)論 0贊 268
萬榮殺人案實(shí)錄
序言：老撾萬榮一對(duì)情侶失蹤，失蹤者是張志新（化名）和其女友劉穎篮灼，沒想到半個(gè)月后忘古，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
沈念sama閱讀 44,355評(píng)論 1贊 303
?護(hù)林員之死
正文獨(dú)居荒郊野嶺守林人離奇死亡诅诱，尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點(diǎn)故事閱讀 36,666評(píng)論 2贊 327
?白月光啟示錄
正文我和宋清朗相戀三年髓堪，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
茶點(diǎn)故事閱讀 38,809評(píng)論 1贊 341
活死人
序言：一個(gè)原本活蹦亂跳的男人離奇死亡干旁，死狀恐怖驶沼，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情争群，我是刑警寧澤回怜，帶...
沈念sama閱讀 34,504評(píng)論 4贊 334
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布，位于F島的核電站换薄，受9級(jí)特大地震影響玉雾，放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜专控，卻給世界環(huán)境...
茶點(diǎn)故事閱讀 40,150評(píng)論 3贊 317
男人毒藥：我在死后第九天來索命
文/蒙蒙一抹凳、第九天我趴在偏房一處隱蔽的房頂上張望。院中可真熱鬧伦腐，春花似錦、人聲如沸失都。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,882評(píng)論 0贊 21
一樁弒父案，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽粹庞。三九已至咳焚，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間庞溜，已是汗流浹背革半。一陣腳步聲響...
開封第一講書人閱讀 32,121評(píng)論 1贊 267
情欲美人皮
我被黑心中介騙來泰國(guó)打工，沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留流码，地道東北人又官。一個(gè)月前我還...
沈念sama閱讀 46,628評(píng)論 2贊 362
代替公主和親
正文我出身青樓，卻偏偏與公主長(zhǎng)得像漫试，于是被迫代替她去往敵國(guó)和親六敬。傳聞我的和親對(duì)象是個(gè)殘疾皇子，可洞房花燭夜當(dāng)晚...
茶點(diǎn)故事閱讀 43,724評(píng)論 2贊 351

第五講：彈性體模擬-數(shù)值方法

推薦閱讀更多精彩內(nèi)容