26. 物理シミュレーション

第14章では、物理ベースアニメーションの道具箱を揃えました。
陽的オイラー法と半陰的オイラー法、ベルレ法とルンゲ＝クッタ法、爆発するバネ、バネ質点系の布、トルクと慣性テンソルを持つ剛体、ブロードフェーズとナローフェーズ、反発係数と力積、そして流体の格子法と粒子法です。
実際のエンジンが使う手法の殆どに名前を付けて、其処で止まりました。
此の章では、其の中身を開けます。

硬いバネを陰的オイラー法で解いて何故決して爆発しないのかを確かめ、剛体を正しく回転させ、分離軸判定で箱同士の接触を見つけ、ゲームエンジンと同じやり方で力積と摩擦で解決します。
其の後、リアルタイムの布の殆どを支える位置ベース動力学とXPBDに移ります。
流体では、ナビエ・ストークス方程式を項ごとに読み、Stable Fluidsを2×2の格子で手計算し、SPHの密度を求め、PIC、FLIP、APIC、MPMを比べます。
最後に髪、破壊、GPUでの物理、決定性を扱います。

此の章は第2部で検定の範囲を超えますが、第14章のシラバス項目物理ベースアニメーションから直接育った内容で、26.5節の終わりの「検定のポイント」で何処に繋がるかを示します。
キャラクタ、スキニング、ラグドールは第25章、煙や炎を関与媒質としてレンダリングする方法は第28章で扱います。
数学は線型代数シリーズの内容です。内積と外積（第4章と第5章）、行列の積と転置（第7章）、逆行列（第10章）、クォータニオン（第13章）を使います。

26.1 時間積分を再び

硬さ

14.2節は、爆発するバネからの3つの逃げ道で終わりました。
ステップを小さくする、減衰を加える、陰的積分を使う、の3つです。
時間の尺度が大きく異なる動きを含むシミュレーションを硬い（スティッフな）系と呼びます。
殆ど伸びない（伸びの振動は数ミリ秒）のにゆっくり揺れる（数秒）布、車のサスペンション、重い箱の山などです。
陽的な手法は、誰にも見えなくても一番速い動きを解像しなければならないので、シーンで一番硬いバネが全体の時間刻みを決めてしまいます。

14.2節のバネを$k=100$ N/m、$m=1$ kgにすると、$\omega =\sqrt{k/m}=10$ rad/s、周期は$2\pi /10\approx 0.63$秒です。
半陰的オイラー法には、つまり秒が必要です。
$h=0.25$秒では爆発し、$x$は1、−5.25、21.3、−85.3、341.3と進みます。
同じバネに付いた10 gの布の粒子なら$\omega =100$ rad/sで、秒が必要です。
本物の布は更に硬いです。

陰的オイラー法の計算例

陰的オイラー法（後退オイラー法）は、全てをステップの終わりで評価します：

未知数が両辺に現れるので、各ステップで方程式を解く事になります。
線型なバネ1つなら手で解けます。
1つ目の式を2つ目に代入すると

計算例。
$x_0=1$、$v_0=0$、$h=0.25$、${\omega }^2=100$なので、$h{\omega }^2=25$、$1+h^2{\omega }^2=1+6.25=7.25$です。

• ステップ1：$v_1=\frac{0-25\cdot 1}{7.25}=-3.44828$、$x_1=\frac{1+0.25\cdot 0}{7.25}=0.13793$。
• ステップ2：$v_2=\frac{-3.44828-25\cdot 0.13793}{7.25}=-0.95125$、$x_2=\frac{0.13793+0.25\cdot (-3.44828)}{7.25}=-0.09988$。

1つ目の式でステップ1を確かめると、$x_1=x_0+h{v}_1=1-0.86207=0.13793$です。
エネルギー$E=\frac{1}{2}m{v}^2+\frac{1}{2}k{x}^2$は50 Jから始まり、1ステップ後は$0.5\cdot {3.44828}^2+50\cdot {0.13793}^2=5.9453+0.9512=6.8966$ J、丁度$\frac{50}{7.25}$です。
毎ステップ、エネルギーは$1+h^2{\omega }^2$で割られます。
同じ数を掛ける陽的オイラー法（14.2節）の鏡像です。
其の為、陰的オイラー法はステップがどれだけ大きくても爆発出来ません。
代償は図に見えています。
本物のバネが保つエネルギーを失うのです。

$h=0.1$秒（$h\omega =1$）ではエネルギーが毎ステップ半分になり、1秒（10ステップ）後には${0.5}^{10}\approx 0.001$しか残りません。
永遠に揺れる筈のバネが止まってしまいました。
60 Hzのステップでは1ステップの損失は$\frac{1}{1+100/3600}=0.973$倍ですが、1秒後には矢張り${0.973}^{60}\approx 0.19$です。
布なら此の人工的な減衰はむしろ歓迎される事が多い（布は元々減衰します）ですが、振り子や軌道では間違いです。
此れが取引です。
安定だが、数値的に減衰する（数値減衰）。

安定領域

何故一方の手法は必ず増え、もう一方は必ず減るのでしょうか。
数値解析では、全ての積分法をテスト方程式

で調べます。
$\lambda$は複素数です。
解$y_0{e}^{\lambda t}$は、実部が負なら減衰し、$\lambda$が純虚数なら振動します。
手法を適用すると$y_{n+1}=R(z)y_n$（$z=h\lambda$）となり、$|R(z)|\le 1$なら安定です。
其の様な$z$の集合が其の手法の安定領域です。

• 陽的オイラー法：$y_{n+1}=y_n+h\lambda y_n$なので$R(z)=1+z$。$-1$を中心とする半径1の円板の内側で安定です。
• 陰的オイラー法：$y_{n+1}=y_n+h\lambda y_{n+1}$なので$R(z)=\frac{1}{1-z}$。$+1$を中心とする半径1の円板の外側で安定で、此れは左半平面全体を含みます。此の性質を持つ手法をA安定と呼びます。

減衰の無いバネは、$\lambda =\pm i\omega$の2つのテスト方程式の様に振る舞います（$(x,v)$を$(v,-{\omega }^{2}x)$に移す2×2行列の固有値です）。
すると$z=\pm ih\omega$は虚軸上にあり、

此れらの2乗が、既に見たエネルギーの倍率です。
虚軸は陽的オイラー法の円板と0でしか接しないので、陽的オイラー法は減衰の無いバネに対してどんな時間刻みでも不安定です。

減衰は負の実数の$\lambda$です。
1 kgに付いた$c=200$ N·s/mのダンパーは$\lambda =-c/m=-200$なので、60 Hzでは$z=-3.33$です。
陽的オイラー法は1ステップごとに速度を$1+z=-2.33$倍し、エネルギーを加えるダンパーになってしまいます（$h\le 2m/c=0.01$秒が必要です）。
陰的オイラー法は$\frac{1}{1+3.33}=0.23$倍で、穏やかに減衰します。
（半陰的オイラー法のバネに対するステップ行列は、の間、長さが丁度1の固有値を持ちます。此れが14.2節の有界なエネルギーです。）

陰的中点法、$R(z)=\frac{1+z/2}{1-z/2}$は、丁度左半平面で安定で、虚軸上では$|R|=1$です。
同じ解く手間で、線型なバネのエネルギーを正確に保ちます。

実際の陰的積分

本物のシーンには何千もの粒子と非線型な力があります（バネの力は両端を結ぶ方向に依存します）。
すると${𝐯}_{n+1}={𝐯}_n+h{𝐌}^{-1}𝐟({𝐱}_{n+1},{𝐯}_{n+1})$は大きな非線型の連立方程式になり、ニュートン法で解きます。
現在の状態の周りで力を線型化し、線型の連立方程式を1回か数回解くのです。
BaraffとWitkin（1998年、14.2節で引用）は、1ステップに1回のニュートン反復を行います：

$𝐌$は（対角の）質量行列、$\frac{\mathrm{\partial }𝐟}{\mathrm{\partial }𝐱}$は剛性行列で、全ての位置が少し動いた時に全ての力がどう変わるかを表します。
1次元のバネ（$\frac{\mathrm{\partial }f}{\mathrm{\partial }x}=-100$、$\frac{\mathrm{\partial }f}{\mathrm{\partial }v}=0$）では$7.25\mathrm{\Delta }v=-25$となり、上と同じ$v_1$が得られます。
布では、行列はバネの端ごとに3×3のブロックを持つので巨大ですが疎で、共役勾配法で反復的に解きます。
映画の布や髪のソルバーは此の様に動き、ゲームでは主に26.3節のより安価な位置ベースの手法を使います。

26.2 剛体の詳細

状態とワールド座標の慣性テンソル

14.5節では、剛体に位置$𝐱$、速度$𝐯$、姿勢、角速度$𝝎$を与え、慣性テンソル$𝐈$でトルクと角加速度を結び付けました。
剛体シミュレータの出来を左右する細部が2つあります。

1つ目に、$𝐈$が一定なのは物体自身の座標系の中だけです。
質量$m$、局所座標のx、y、z軸に沿った辺の長さ$w$、$h$、$d$の箱では、対角行列になります：

物体が行列$𝐑$で回転すると、ワールドから見た慣性は変わります：

（右から左へ読みます。ワールドのベクトルを物体の空間に回し、物体の慣性を掛け、元に戻します。回転では${𝐑}^{-1}={𝐑}^{\mathrm{\top }}$なので、逆行列は安価です。線型代数シリーズ第10章。）
エンジンは${𝐈}_{\text{body}}^{-1}$を一度だけ保存し、毎ステップ${𝐈}_{\text{world}}^{-1}$を作り直します。

計算例。
$m=12$ kg、$2\times 1\times 1$ mの箱では、${𝐈}_{\text{body}}=\mathrm{diag}(1\cdot (1+1),1\cdot (4+1),1\cdot (4+1))=\mathrm{diag}(2,5,5)$ kg·m²です。
z軸周りに45°回転させます。
$c=s=\frac{\sqrt{2}}{2}$とすると、左上の2×2ブロックは

になり、$I_{zz}=5$は変わりません。
此処で、ワールドのx軸周りの回転として角運動量$𝐋=(1,0,0)$を与えます。
角速度は$𝝎={𝐈}_{\text{world}}^{-1}𝐋$です。
2×2ブロックの逆行列はなので

箱は押した軸の周りには回りません。
$𝝎$は$𝐋$から23°ずれています。
$𝐋$と$𝝎$が平行になるのは、3本の慣性主軸（此処では物体自身のx、y、z）の周りだけです。
投げた本がぐらぐら揺れるのは其の為です。

2つ目に、$𝐈$が変わるので、$𝐅=m𝐚$の回転版には余分な項が付きます（オイラーの運動方程式）：

ジャイロ項$𝝎\times (𝐈𝝎)$は、回転するラケットをひっくり返らせます。
此の項は$𝝎$の2次なので、半陰的オイラー法では速く回る細長い物体にエネルギーが注ぎ込まれて爆発する事があります。
其の為エンジンは、此の項を捨てる（安定だが「安定過ぎる」）か、1回のニュートン反復で陰的に積分します。
Erin CattoがGDC 2015の数値解法の講演で示した方法です。

姿勢の積分

14.5節では姿勢のクォータニオンを$𝐪\leftarrow 𝐪+\frac{h}{2}(0,𝝎)𝐪$で更新し、正規化し直しました。
此れはクォータニオンの微分$\dot{𝐪}=\frac{1}{2}(0,𝝎)𝐪$に陽的オイラー法を使った物です。
$𝝎$はワールド座標、積は線型代数シリーズ第13章のハミルトン積です。

計算例。
単位元$𝐪=(1,0,0,0)$（$(w,x,y,z)$の順）から始め、$𝝎=(0,0,2)$ rad/sでz軸周りに$h=0.1$秒回します。
$(0,0,0,2)(1,0,0,0)=(0,0,0,2)$なので、$𝐪+0.05\cdot (0,0,0,2)=(1,0,0,0.1)$です。
長さは$\sqrt{1.01}=1.004988$で、正規化すると$(0.995037,0,0,0.099504)$です。
単位クォータニオン$(\cos \frac{\theta }{2},\sin \frac{\theta }{2}\widehat{𝐚})$は$\theta$の回転なので、此れは$\theta =2\arctan \frac{0.099504}{0.995037}=0.19934$ radです。
正確な答えは$\omega h=0.2$ radです。
物体は0.3%遅く回り、正規化しなければ長さは1ステップごとに1.004988倍に伸びます。
正確な代わりの方法は指数写像です。
$\widehat{𝝎}$周りの$\Vert 𝝎\Vert h$の回転、此処では$(\cos 0.1,0,0,\sin 0.1)$を、$𝐪$に左から掛けます。

接触の生成：分離軸判定

第14章でナローフェーズの仕事を示しました。
接触点、接触法線、めり込み量を見つける事です。
凸多角形や凸多面体では、分離軸定理（SAT）が定番の道具です。
2つの凸形状が重ならないのは、其れらの射影（「影」）が重ならない軸が存在する時、そして其の時に限ります。
2次元の多角形では両方の辺の法線を調べれば十分で、箱なら各箱の2本の軸です。

中心$𝐜$、単位軸$𝐮$、$𝐯$、半分の大きさ$e_u$、$e_v$の箱を単位軸$𝐧$に射影すると、$𝐜\cdot 𝐧$を中心とする半径

の区間になり、2つの箱は其の軸上で$r_A+r_B-|({𝐜}_B-{𝐜}_A)\cdot 𝐧|$だけ重なります。
負の値は隙間です。

計算例。
箱Aは原点にあり、軸に平行で、半分の大きさは$(2,1)$です。
箱Bは半分の大きさ$(1,0.5)$で、軸が${𝐮}_B=(0.8,0.6)$、${𝐯}_B=(-0.6,0.8)$になる様に回転しています（約36.9°）。
4本の軸で、Aの半径は2、1、$2\cdot 0.8+1\cdot 0.6=2.2$、$2\cdot 0.6+1\cdot 0.8=2.0$、Bの半径は$1\cdot 0.8+0.5\cdot 0.6=1.1$、$1\cdot 0.6+0.5\cdot 0.8=1.0$、1、0.5です。

Bの中心が$(3,1.8)$の時、${𝐮}_B$に沿った距離は$3\cdot 0.8+1.8\cdot 0.6=3.48$で、$2.2+1=3.2$より大きいです。
隙間が0.28あるので、2つの箱は離れています。
Aの2本の軸（AABBの判定が見るのは其れだけです）では両方とも重なっているにも関わらずです。
Bが$(2.5,1.5)$の時は4本全ての軸で重なるので、箱は衝突しています。
重なりが最小の軸が接触法線で、$𝐧={𝐮}_B=(0.8,0.6)$です（$({𝐜}_B-{𝐜}_A)\cdot {𝐮}_B>0$なので、AからBへ向きます）。
重なり0.3がめり込み量で、Bを$0.3𝐧$動かせば離れます。

箱の上に箱が載っている時、法線と深さだけでは足りません。
接触点が1つだけだと、箱がぐらつくからです。
エンジンはクリッピングで、2次元なら最大2点（3次元なら4点）の接触点の集合（コンタクトマニフォールド）を作ります。
軸はBから来たので、Aに向いたBの面が参照面で、外向きの法線は$-{𝐮}_B$です。
其れと最も反対向きのAの面、つまり$(2,-1)$から$(2,1)$までのAの右辺（其の法線$(1,0)$と${𝐮}_B$の内積0.8が最大です）が入射面です。
入射辺を参照面の側面で切り取ります。
第6章の6.3節のサザーランド・ホッジマンのクリッピングと全く同じです。
Bの座標では点$(2,y)$は$v=(𝐩-{𝐜}_B)\cdot {𝐯}_B=0.8y-0.9$を持ち、側面$v\ge -0.5$は辺を$y=0.5$で切ります。
切り取った点のうち、参照面の裏側にある物を残します。
$(2,1)$は$u=(𝐩-{𝐜}_B)\cdot {𝐮}_B=-0.7$なので深さ$1-0.7=0.3$、$(2,0.5)$は$u=-1.0$で深さ0、丁度触れています。
深い点と触れている点の2つで、BがAの中へ回り込むのを防ぐには十分です。

3次元では、2つの箱の候補の軸は15本あります。
各箱の面の法線3本ずつと、Aの辺とBの辺の外積9本です（面が関わらずに辺同士がぶつかる事があるからです）。

GJKとEPA

SATには候補の軸のリストが必要です。
箱や凸多角形なら良いのですが、球、カプセル、任意の凸包では困ります。
GJK法（Gilbert・Johnson・Keerthi、1988年）は、1つの質問に答えられるどんな凸形状でも扱えます。
其の質問がサポート関数です。
与えられた方向$𝐝$に一番遠い形状の点は何処か？
球なら$𝐜+r\widehat{𝐝}$、箱なら$𝐝$の側の角です。

鍵はミンコフスキー差$A-B=\{𝐚-𝐛\}$、Aの点とBの点の全ての差の集合です。
AとBが重なるのは、ある$𝐚=𝐛$がある時、つまり原点が$A-B$の内側にある時です。
其のサポート関数は$s_{A-B}(𝐝)=s_A(𝐝)-s_B(-𝐝)$なので、$A-B$を作る必要はありません。
GJKは$A-B$のサポート点で作る単体（点、線分、三角形、四面体）を原点に向けて育てます。
原点を囲めば形状は重なっていて、サポート点が原点を越えられなければ離れていて、一番近い単体の要素が距離を与えます。
重なっている時は、EPA（expanding polytope algorithm、van den Bergen、2001年）が最後の単体を多面体に膨らませ、原点に一番近い$A-B$の面を見つけます。
其の法線が接触法線、距離がめり込み量です。

衝突している箱と$𝐝={𝐮}_B$では、$s_A\cdot 𝐝=(2,1)\cdot {𝐮}_B=2.2$、$s_B(-𝐝)\cdot 𝐝=({𝐜}_B-{𝐮}_B)\cdot {𝐮}_B=1.9$なので、$A-B$は其の方向に原点を$0.3$越えています。
SATが見つけた深さで、EPAは其れを全方向の最小値として見つけます。

接触拘束と力積

14.6節では、正面衝突を法線に沿った1つの力積で解決しました。
回転があると、問題になるのは接触点の速度$𝐯+𝝎\times 𝐫$です。
$𝐫$は重心から接触点へのオフセットです。
接触拘束は、此の点が相手の物体の中へ動いてはならない事、$v_n=(𝐯+𝝎\times 𝐫)\cdot 𝐧\ge 0$、そして接触は押せても決して引けない事、法線方向の力積${\lambda }_n\ge 0$を表します。
$𝐫$で方向$𝐧$に沿った力積$\lambda$は、両方の速度を変えます：

従って$v_n$は$\lambda /m_{\text{eff}}$だけ変わります。
$m_{\text{eff}}$は有効質量です：

（2つ目の物体が静止物でなければ、其の分の同じ項も加えます。）
近づく動きを完全に止める（$e=0$）には、${\lambda }_n=-m_{\text{eff}}{v}_n$を加えます。

摩擦はクーロン摩擦の法則に従います。
接線方向の力積は、法線方向の力積の高々$\mu$倍、$|{\lambda }_t|\le \mu {\lambda }_n$です。
3次元では、許される接触の力積は法線の周りの半頂角$\arctan \mu$の摩擦円錐を作ります。
エンジンは普通、2本の接線方向で近似し、四角錐にします。
円錐の内側では接触は止まっていて（静止摩擦）、表面上では滑ります。

計算例。
木箱（2次元、$m=1$ kg、$I=\frac{m(1^2+1^2)}{12}=\frac{1}{6}$ kg·m²の単位正方形）が$𝐯=(0,-2)$ m/sで落ち、右下の角の下にある細い支柱に受け止められます。
$𝐫=(0.5,-0.5)$、$𝐧=(0,1)$、$𝐭=(1,0)$、$\mu =0.5$、$e=0$です。
2次元では外積$𝐫\times 𝐧=r_x{n}_y-r_y{n}_x$は数で、此処では$0.5$です。
又、$𝝎\times 𝐫=\omega (-r_y,r_x)=\omega (0.5,0.5)$です。

• 有効質量：法線は$\frac{1}{m_{\text{eff}}}=1+{0.5}^2\cdot 6=2.5$、接線も同じです（$𝐫\times 𝐭=0-(-0.5)(1)=0.5$）。
• 法線：$v_n=-2$なので${\lambda }_n=2/2.5=0.8$。すると$𝐯=(0,-1.2)$、$\omega =6\cdot 0.5\cdot 0.8=2.4$ rad/s。角は$v_n=-1.2+2.4\cdot 0.5=0$で止まります。
• 摩擦：角は$v_t=0+2.4\cdot 0.5=1.2$ m/sで滑っています。止めるには${\lambda }_t=-1.2/2.5=-0.48$が必要ですが、クーロンの法則は$\mu {\lambda }_n=0.4$までしか許さないので${\lambda }_t=-0.4$で、角は滑り続けます。すると$𝐯=(-0.4,-1.2)$、$\omega =2.4+6\cdot 0.5\cdot (-0.4)=1.2$。

ところが今度は、角が再び支柱の中へ動きます。
$v_n=-1.2+1.2\cdot 0.5=-0.6$です。
摩擦の力積が、法線の力積が頼りにしていた回転を奪ったのです。
1つの拘束を解くと別の拘束が壊れるので、反復します。

シーケンシャルインパルスと射影ガウス＝ザイデル法

ゲームエンジンは拘束を1つずつ訪れ、今の拘束を直す力積を加え、全体を何周か繰り返します。
Erin Cattoは此れをシーケンシャルインパルスとしてGDCの講演（「Iterative Dynamics with Temporal Coherence」2005年、「Fast and Simple Physics using Sequential Impulses」2006年）で広めました。
Box2Dのソルバーです。
数学的には射影ガウス＝ザイデル法（PGS）で、全ての接触の連立方程式にガウス＝ザイデル反復を行い、各未知数を更新の後に許される範囲へ射影し戻します。

鍵となる規則は、ステップ全体で接触ごとに累積インパルスを持ち、各増分ではなく累積の値をクランプする事です。
増分は負になっても構いません（大き過ぎた力積を一部取り戻すのです）。
法線の合計が$\ge 0$、摩擦の合計が$\pm \mu {\lambda }_n$の内側に留まりさえすれば良いのです。
木箱の続きです：

誤差は1反復ごとに0.3倍になります。
厳密な答えは、両方の条件を一緒に解いて得られます。
接触は滑るので${\lambda }_t=-\mu {\lambda }_n=-0.5{\lambda }_n$、すると$\omega =6\cdot 0.5({\lambda }_n+{\lambda }_t)=1.5{\lambda }_n$、$v_n=-2+{\lambda }_n+0.5\cdot 1.5{\lambda }_n=-2+1.75{\lambda }_n=0$から${\lambda }_n=\frac{8}{7}$です。
角は$v_t=-\frac{4}{7}+0.5\cdot \frac{12}{7}=\frac{2}{7}$ m/sで右に滑り、摩擦は正しく左に押します。

接触が多いと収束は遅くなります（高い山では、一番上の箱の重さが全ての接触を通って下まで伝わらなければなりません）。
其の為エンジンは決まった回数、典型的には4〜20回だけ反復し、残りの誤差を受け入れます。

#include <algorithm>
#include <cstdio>
#include "linalg.hh"

// 2D cross products: r x n is a scalar, w x r is a vector.
float cross2(Vec2 a, Vec2 b) { return a.x * b.y - a.y * b.x; }
Vec2 crossW(float w, Vec2 r) { return {-w * r.y, w * r.x}; }

struct Body {
  float invMass, invInertia;
  Vec2 v;
  float w;  // angular velocity (rad/s, counter-clockwise)
};

// One contact point at offset r from the center, normal n, tangent t, against a static object.
struct Contact {
  Vec2 r, n, t;
  float mu;
  float Pn = 0.0f, Pt = 0.0f;  // accumulated impulses (keep them for warm starting)
};

void applyImpulse(Body &b, Vec2 r, Vec2 dir, float lambda) {
  b.v = b.v + (lambda * b.invMass) * dir;
  b.w += b.invInertia * cross2(r, dir) * lambda;
}

void solveContact(Body &b, Contact &c) {
  // normal: effective mass 1 / (1/m + (r x n)^2 / I)
  float rn = cross2(c.r, c.n);
  float kn = b.invMass + rn * rn * b.invInertia;
  float vn = dot(b.v + crossW(b.w, c.r), c.n);
  float old = c.Pn;
  c.Pn = std::max(old - vn / kn, 0.0f);  // clamp the accumulated impulse: never pull
  applyImpulse(b, c.r, c.n, c.Pn - old);
  // friction: Coulomb, |Pt| <= mu * Pn
  float rt = cross2(c.r, c.t);
  float kt = b.invMass + rt * rt * b.invInertia;
  float vt = dot(b.v + crossW(b.w, c.r), c.t);
  old = c.Pt;
  c.Pt = std::clamp(old - vt / kt, -c.mu * c.Pn, c.mu * c.Pn);
  applyImpulse(b, c.r, c.t, c.Pt - old);
}

int main() {
  Body crate{1.0f, 6.0f, {0.0f, -2.0f}, 0.0f};  // m = 1 kg, I = 1/6 kg m^2
  Contact c{{0.5f, -0.5f}, {0.0f, 1.0f}, {1.0f, 0.0f}, 0.5f};
  for (int it = 1; it <= 5; ++it) {
    solveContact(crate, c);
    float vn = dot(crate.v + crossW(crate.w, c.r), c.n);
    printf("iteration %d: Pn %.4f  Pt %.4f  v (%.4f, %.4f)  w %.4f  vn %.4f\n", it, c.Pn, c.Pt,
           crate.v.x, crate.v.y, crate.w, vn);
  }
  return 0;
}

上の表の反復1〜4と、5回目の${\lambda }_n=1.1401$、$v_n=-0.0049$を出力します。

ジョイントとモーター

第25章の25.6節のラグドールは、ジョイントで繋がれた剛体です。ジョイントは、接触と同じ周回の中で解かれる拘束です。
ボールジョイントは各物体のアンカー点を一緒に保ちます。アンカー点の速度${𝐯}_A+{𝝎}_A\times {𝐫}_A$と${𝐯}_B+{𝝎}_B\times {𝐫}_B$が等しくなければならず、x、y、zに沿った3つの式を其々上の有効質量で解きますが、クランプはしません。ジョイントは押すだけでなく引く事も出来るからです。
ヒンジは2つの回転の拘束を加え、回転軸を1本だけ自由にします。可動域の制限は接触の様な不等式で、限界に達した時だけ働きます。
モーターは目標を持つ速度の拘束で、例えば「ヒンジ軸周りの相対角速度が${\omega }_{\text{target}}$」です。其の累積インパルスを$\pm {\tau }_{max}h$にクランプして、強さを有限にします。PD制御は代わりに、26.3節のコンプライアンスの様に、目標の角度へ向かう柔らかいバネとして働きます。

めり込み、ウォームスタート、積み重ね、スリープ

速度の力積は物体が更にめり込むのを止めますが、既にあるめり込みは残り、溜まって行きます。
定番の対策、バウムガルテ安定化は、$v_n=0$の代わりに、深さに比例した小さな離れる速度を求めます：

許容めり込み（スロップ）$d_{\text{slop}}$は数ミリの重なりを残し、静止した接触が触れたり離れたりしてちらつかない様にします。
$\beta$は0.1〜0.3程度で、残りを数ステップかけて穏やかに直します。
$\beta =0.2$、$h=\frac{1}{60}$秒、深さ2 cm、スロップ0.5 cmなら、目標は$0.2\cdot 60\cdot 0.015=0.18$ m/sです。
此のバイアスはエネルギーを加える（箱が「ポン」と弾ける）ので、多くのエンジンは位置を別のパスで直すか、26.3節のコンプライアンスの様なソフトな拘束を使います。

ウォームスタート：接触はフレームを跨いで存続する（「Aの頂点2とBの面1」の様な特徴IDで対応付けます）ので、各ステップは前のステップの累積インパルスを加える所から始めます。
静止した山では力積は殆ど変わらないので、反復は殆ど答えの所から始まります。
剛体エンジンの定番のベンチマークである安定した積み重ねには、更に良い接触点の集合、十分な反復回数、地面から上へ順に解く事が必要です。

スリープ：並進と回転の速さが小さな閾値より下にある時間が暫く（例えば約0.5秒）続いた物体は凍結され、何かが触れるまで処理を飛ばされます。
触れ合う物体は一緒に眠り一緒に起きるアイランドを作ります。
典型的なゲームのレベルでは殆ど全てが殆ど常に眠っていて、其のおかげで何千もの物体を扱えます。

すり抜け（14.6節）に対しては、連続衝突検出が衝突の時刻を求めます。
例えば保守的前進法（conservative advancement）は、GJKの距離が安全だと保証する分だけ進む事を繰り返します。

26.3 位置ベース動力学

考え方

力と力積は速度を変え、位置が其れに付いて行きます。
位置ベース動力学（PBD、Müller・Heidelberger・Hennix・Ratcliff、VRIPHYS 2006、論文誌版は2007年）は此れを逆にします。
拘束が成り立つ様に位置を直接動かし、位置がどれだけ動いたかから速度を求めるのです。
Thomas Jakobsenの「Advanced Character Physics」（2001年、IO InteractiveのHitman: Codename 47の為の物）のベルレ法の布やラグドールを体系化した物です。
各ステップは次の通りです：

1. 予測：${𝐯}_i\leftarrow {𝐯}_i+h{w}_i{𝐟}_{\text{ext},i}$、次に${𝐩}_i={𝐱}_i+h{𝐯}_i$。$w_i=1/m_i$は逆質量で、固定された粒子は$w_i=0$です。
2. 拘束の射影：決まった回数だけ、全ての拘束（距離、曲げ、衝突）を訪れ、其れが成り立つ様に其の粒子の${𝐩}_i$をガウス＝ザイデル式に動かします。
3. 更新：${𝐯}_i=({𝐩}_i-{𝐱}_i)/h$、${𝐱}_i={𝐩}_i$。其の後、衝突した粒子に摩擦と反発を適用します。

爆発しようがありません。
射影は粒子を拘束の誤差より遠くへは決して動かさず、速度は実際に起きた動きそのものだからです。
此の頑丈さの為に、ゲームの布、ロープ、ソフトボディ、パーティクルはPBDが主流になりました。

距離拘束の射影

拘束は、0であるべき関数$C({𝐩}_1,\dots ,{𝐩}_n)$です。
距離$d$を保つべき2つの粒子なら$C=\Vert {𝐩}_1-{𝐩}_2\Vert -d$です。
PBDは粒子を$C$の勾配（$C$が一番速く変わる方向）に沿って動かします。
重い粒子が少なく動き、運動量が保存される様に、逆質量で重み付けします：

距離拘束では${\mathrm{\nabla }}_{{𝐩}_1}C=\widehat{𝐧}=\frac{{𝐩}_1-{𝐩}_2}{\Vert {𝐩}_1-{𝐩}_2\Vert }$、${\mathrm{\nabla }}_{{𝐩}_2}C=-\widehat{𝐧}$で、どちらも長さ1なので

計算例（14.4節のバネを硬い棒として）。
${𝐩}_1=(0,0)$、${𝐩}_2=(3,4)$、自然長$d=4$です。
距離は5、$C=1$、$\widehat{𝐧}=(-0.6,-0.8)$です。

• 等しい質量、$w_1=w_2=1$：$\mathrm{\Delta }{𝐩}_1=-0.5\cdot 1\cdot (-0.6,-0.8)=(0.3,0.4)$、$\mathrm{\Delta }{𝐩}_2=(-0.3,-0.4)$。新しい位置は$(0.3,0.4)$と$(2.7,3.6)$で、距離は$\sqrt{{2.4}^2+{3.2}^2}=4$。
• ${𝐩}_1$を固定、$w_1=0$：${𝐩}_2$だけが$(-0.6,-0.8)$丸ごと動き、$(2.4,3.2)$へ。
• $w_1=1$、$w_2=3$（2つ目の粒子が3倍軽い）：${𝐩}_1$は補正の4分の1動いて$(0.15,0.2)$へ、${𝐩}_2$は4分の3動いて$(2.55,3.4)$へ。重心は動きません。

1つの拘束は1回の射影で正確に解けます。
しかし布には粒子を共有する何千もの拘束があり、1つを直すと隣が壊れるので、PBDは上のシーケンシャルインパルスと全く同じ様に全体を何周もします（非線型のガウス＝ザイデル反復です）。

硬さが反復回数に依存する

拘束を柔らかくするには、PBDは各補正に硬さ$k\in [0,1]$を掛けます。
しかし補正は反復ごとに掛かるので、$n$回の反復の後に残る誤差は$(1-k{)}^n$です。
$k=0.5$なら、1回で0.5、4回で0.0625、16回で0.000015です。
同じ材質が、反復回数を増やすと硬くなり、時間刻みでも変わってしまいます。
PBDの論文は$k\prime =1-(1-k{)}^{1/n}$で補います（$k=0.5$、$n=4$なら$k\prime =0.159$で、$(1-0.159{)}^4=0.5$）。
此れで反復回数への依存は直りますが、時間刻みへの依存は直りません。
下の布が其の効果を示しています。
同じシミュレーション（2つの角から吊るした13×13粒子、1.2 mの布、$k=1$の構造とせん断の距離拘束、60 Hzで3秒）を、1ステップあたり1回、4回、16回の反復で行った物です。

「硬い」拘束でも、途中で止まる反復ソルバーは布を伸びたままにします。
特に、下の全ての重さを支える固定点の近くです。
反復が少なければゴム、多ければ木綿になり、結果は解像度にも依存します。

XPBD：コンプライアンス

XPBD（extended position-based dynamics、Macklin・Müller・Chentanez、Motion in Games 2016）は、各拘束に物理的なコンプライアンス$\alpha$を与えて此れを直します。
コンプライアンスは硬さの逆数（バネなら$\alpha =1/k$ m/N、硬い棒なら0）です。
更に拘束ごとにラグランジュ乗数$\lambda$を累積します。
此れは拘束の力積です。
$\tilde{\alpha }=\alpha /h^2$として、各射影は次を計算します：

$\lambda$は各ステップの始めに0に戻します。
$\alpha =0$なら、此れは丁度PBDの射影です。

計算例。
同じ2つの粒子を、今度は$k=100$ N/mのバネで繋ぎます。
$\alpha =0.01$ m/Nで、$h=0.1$秒なら$\tilde{\alpha }=0.01/0.01=1$です。
1回目の反復：$\mathrm{\Delta }\lambda =\frac{-1-0}{2+1}=-\frac{1}{3}$なので、$\mathrm{\Delta }{𝐩}_1=1\cdot (-0.6,-0.8)\cdot (-\frac{1}{3})=(0.2,0.2667)$、$\mathrm{\Delta }{𝐩}_2=(-0.2,-0.2667)$。
新しい長さは$4.3333$なので$C=\frac{1}{3}$です。
2回目の反復：$\mathrm{\Delta }\lambda =\frac{-\frac{1}{3}-1\cdot (-\frac{1}{3})}{3}=0$で、何回反復しても何も変わりません。
答えは射影ではなく釣り合いです。
各粒子は慣性が行きたかった所から$\frac{1}{3}$ m引き離されていて、其れには$m\cdot \frac{1/3}{h^2}=33.3$ Nの力が要ります。
一方$\frac{1}{3}$ m伸びたバネは$kC=100\cdot \frac{1}{3}=33.3$ Nで引きます。
此れはバネの陰的オイラー法の1ステップを射影で解いた物で、結果はもう反復回数に依存しません。
反復を増やすと、大きな系が同じ答えに近づくだけです。

Macklinらは後に（「Small Steps in Physics Simulation」、SCA 2019）、同じ計算量なら、1回の反復の小さなサブステップを沢山行う方が、少ないステップで沢山反復するより良い事を示しました。

#include <cstdio>
#include "linalg.hh"

// XPBD distance constraint C = |p1 - p2| - d between two particles with inverse masses w1, w2.
// alpha is the compliance (1 / stiffness), h the time step, lambda the accumulated multiplier.
// With alpha = 0 this is exactly the PBD projection.
void solveDistance(Vec3 &p1, Vec3 &p2, float w1, float w2, float d, float alpha, float h,
                   float &lambda) {
  Vec3 delta = p1 - p2;
  float len = length(delta);
  if (len < 1e-9f || w1 + w2 == 0.0f) return;
  Vec3 n = delta / len;             // gradient of C with respect to p1
  float C = len - d;
  float alphaTilde = alpha / (h * h);
  float dLambda = (-C - alphaTilde * lambda) / (w1 + w2 + alphaTilde);
  lambda += dLambda;
  p1 = p1 + (w1 * dLambda) * n;
  p2 = p2 - (w2 * dLambda) * n;
}

int main() {
  const float alphas[] = {0.0f, 0.01f};  // rigid (PBD) and a spring of k = 100 N/m
  for (float alpha : alphas) {
    Vec3 p1{0, 0, 0}, p2{3, 4, 0};
    float lambda = 0.0f;
    for (int it = 1; it <= 2; ++it) {
      solveDistance(p1, p2, 1.0f, 1.0f, 4.0f, alpha, 0.1f, lambda);
      printf("alpha %.2f, iteration %d: p1 (%.4f, %.4f)  p2 (%.4f, %.4f)  "
             "length %.4f  lambda %.4f\n",
             alpha, it, p1.x, p1.y, p2.x, p2.y, length(p2 - p1), lambda);
    }
  }
  return 0;
}

$\alpha =0$ではPBDの結果（長さ4.0000、$\lambda =-0.5$）、$\alpha =0.01$ではXPBDの結果（長さ4.3333、$\lambda =-0.3333$）を出力し、どちらでも2回目の反復は何も変えません。

PBDの布

PBDの布は14.4節のバネ質点系の布と同じ粒子の格子を使い、バネの代わりに拘束を置きます：

• 伸び：構造の辺に沿った距離拘束（せん断に対しては、せん断の辺か三角形の面積の拘束）。
• 曲げ：曲げバネの様に粒子を1つ飛ばす距離拘束か、隣り合う三角形の間の二面角の拘束。後者は伸びと独立です。
• 衝突：破られている間だけ射影する不等式の拘束$C\ge 0$。例えば体の表面や符号付き距離場（第21章）に対する$C=(𝐩-𝐪)\cdot 𝐧-\text{thickness}$や、自己衝突の為の粒子と三角形の拘束です。

ベルレ法の布は同じ物をベルレ法で書いた物です。
速度は${𝐱}_n-{𝐱}_{n-1}$に隠れているので、位置を動かして拘束を緩和すれば、速度は勝手に更新されます。

PBDよりずっと前に、Provot（Graphics Interface 1995）は伸び易いバネ質点系の布をひずみ制限で直しました。
各ステップの後、臨界の伸び率（彼の吊るした布では10%）より伸びたバネは、両端を丁度其の限界まで戻します。
両端が自由なら両方を、片方が固定なら自由な方だけを動かします。
此れは不等式$\Vert {𝐩}_1-{𝐩}_2\Vert \le 1.1d$のPBDの射影で、布のソルバーは今でも柔らかいバネと此の様な硬い限界を組み合わせます。

シェイプマッチング

静止形状に戻る柔らかい物体（ゼリー、アニメの潰しと伸ばし）には、シェイプマッチング（Müller・Heidelberger・Teschner・Gross、「Meshless Deformations Based on Shape Matching」、SIGGRAPH 2005）が使え、メッシュは全く要りません。
各ステップで、静止位置${𝐪}_i$を現在の位置${𝐱}_i$に一番良く合わせる剛体の動き（回転$𝐑$と平行移動）を求め、各粒子を目標位置${𝐠}_i$へ割合$\alpha$だけ引き寄せます：

${𝐜}_0$は静止時の重心、$𝐑$は$𝐀$の回転部分（極分解$𝐀=𝐑𝐒$）です。
2次元では回転角は単純に$\theta =\mathrm{atan2}(A_{21}-A_{12},A_{11}+A_{22})$です。

計算例。
4つの単位質量の正方形が$𝐪=(-1,-1),(1,-1),(1,1),(-1,1)$で静止していて、${𝐜}_0=(0,0)$です。
其れが90°回転し、へこみました：$𝐱=(1,-1),(1,1),(-0.4,0.4),(-1,-1)$。
重心は$𝐜=\frac{(0.6,-0.6)}{4}=(0.15,-0.15)$で、オフセット${𝐱}_i-𝐜$は$(0.85,-0.85)$、$(0.85,1.15)$、$(-0.55,0.55)$、$(-1.15,-0.85)$です。
${𝐪}_i$との外積（直積）を足すと：

なので$\theta =\mathrm{atan2}(6.8,0)=90\text{°}$です。
へこみに騙されませんでした。
目標は、静止形状を90°回して$𝐜$へ動かした物で、$(1.15,-1.15)$、$(1.15,0.85)$、$(-0.85,0.85)$、$(-0.85,-1.15)$です。
$\alpha =0.5$なら、へこんだ角は$(-0.4,0.4)$から目標へ半分動き、$(-0.625,0.625)$になります。

$𝐑$の代わりに行列$𝐀$全体で合わせると伸びの変形が許され、重なり合うクラスタに分けると局所的な揺れが出ます。

有限要素法（概要）

バネ質点系とPBDは速いですが、硬さがメッシュに依存し、「此れはゴムだ」と言う事が難しいです。
有限要素法（FEM）は代わりに連続体力学から出発します。
材料を四面体（2次元では三角形）に分け、各要素の中で変形は線型だと仮定し、変形勾配

で表します。
${𝐃}_m$の列は静止時の要素の辺ベクトル、${𝐃}_s$の列は同じ辺の今のベクトルです。
ひずみは変形を、応力は面積あたりの内力を測り、線型弾性は行列の形のフックの法則で両者を結びます。
其処で使うのがヤング率（ゴムで約0.01〜0.1 GPa、鋼で200 GPa）とポアソン比（伸ばした時にどれだけ細くなるか、ゴムでは0.5近く）です。

計算例。
静止時の角が$(0,0)$、$(1,0)$、$(0,1)$の2次元の三角形は${𝐃}_m=𝐈$です。
x方向に$(1.1,0)$まで伸ばすと、$𝐅=\mathrm{diag}(1.1,1)$、線型ひずみは$𝜺=\frac{1}{2}(𝐅+{𝐅}^{\mathrm{\top }})-𝐈=\mathrm{diag}(0.1,0)$です。
今度は90°回すだけにして、$(0,0)$、$(0,1)$、$(-1,0)$にします。
、$𝜺=\mathrm{diag}(-1,-1)$で、全く変形していない三角形が100%圧縮された事になります。
其の為、線型のFEMは回転する物体に存在しない力を生みます。
共回転FEMは、シェイプマッチングと同じ極分解で先に回転を取り除き（$𝐅=𝐑𝐒$、ひずみは$𝐒-𝐈$から）、非線型のモデルはどんな回転でも0になるグリーンのひずみ$\frac{1}{2}({𝐅}^{\mathrm{\top }}𝐅-𝐈)$を使います。
要素の力は26.1節の疎な剛性行列に組み立てられ、陰的に積分されます。
映画では筋肉や肉にFEMを使い、執筆時点ではPhysXがGPUのFEMソフトボディを提供しています。

26.4 流体の詳細

ナビエ・ストークス方程式を項ごとに

速度場$𝐮(𝐱,t)$、圧力$p$、一定の密度$\rho$、動粘性係数$\nu$の流体について、14.7節の非圧縮のナビエ・ストークス方程式は

です。
$\mathrm{\nabla }$（ナブラ）は偏微分を並べたベクトル$\left(\frac{\mathrm{\partial }}{\mathrm{\partial }x},\frac{\mathrm{\partial }}{\mathrm{\partial }y},\frac{\mathrm{\partial }}{\mathrm{\partial }z}\right)$で、各項は1つの物理的な効果です：

• $-(𝐮\cdot \mathrm{\nabla })𝐮$、移流：流体が自分の速度を運びます。渦の原因で、非線型の難しい部分です。
• $-\frac{1}{\rho }\mathrm{\nabla }p$、圧力：流体は圧力の高い所から低い所へ押されます。圧力は、流体を非圧縮に保つ為に必要な値になります。
• $\nu {\mathrm{\nabla }}^2𝐮$、粘性：速度が隣へ拡散します。${\mathrm{\nabla }}^2$（ラプラシアン）は「隣の平均引く自分」です。蜂蜜は$\nu$が大きく、空気はとても小さいです。
• $𝐟$：重力、浮力、かき混ぜる手などの外力。
• $\mathrm{\nabla }\cdot 𝐮=0$：発散、つまり体積あたりの正味の流出が0です。セルに流れ込んだ物は流れ出なければなりません。此れが非圧縮性です。

流体が運ぶ量（煙の密度、温度、染料）は全て、移流の部分だけに従います：$\frac{\mathrm{\partial }s}{\mathrm{\partial }t}=-(𝐮\cdot \mathrm{\nabla })s$。

Stable Fluids

14.7節で名前を出したJos StamのStable Fluids（SIGGRAPH 1999）は、此の方程式を格子の上で演算子分割で解きます。
1ステップで項を1つずつ順に適用し、其々にどんな時間刻みでも安定な手法を使うのです。

1. 力を加える：$𝐮\leftarrow 𝐮+h𝐟$。
2. セミラグランジュ法で移流する：格子点の新しい値を求めるには、速度に沿って流体が来た所$𝐱-h𝐮(𝐱)$まで逆に辿り、其処で古い場を補間します。
3. 陰的に拡散する：$(𝐈-\nu h{\mathrm{\nabla }}^2){𝐮}_{\text{new}}=𝐮$を解きます（又陰的オイラー法なので、どんな$\nu$と$h$でも大丈夫です）。煙や水では、段階2の数値拡散だけで既に十分過ぎるので、此の段階を省く事が多いです。
4. 射影する：圧力のポアソン方程式を解いて其の勾配を引き、発散を取り除きます。此れはヘルムホルツ・ホッジ分解を使います。どんなベクトル場も発散の無い部分と勾配の和で、射影は前者だけを残します。

計算例：セミラグランジュ法の移流。
間隔1のセルの1次元の列が、$x=0,\dots ,5$に煙の密度$s=(0,0,1,1,0,0)$を持ち、速度は何処でも$u=0.6$、$h=1$です。
$x=4$のセルでは、流体は$4-0.6=3.4$、$s_3=1$と$s_4=0$の間から来ました。
$0.6\cdot 1+0.4\cdot 0=0.6$です。
$x=3$は2.4から来て$0.6\cdot 1+0.4\cdot 1=1$、$x=2$は1.4から来て$0.6\cdot 0+0.4\cdot 1=0.4$です。
新しい列は$(0,0,0.4,1,0.6,0)$です。
塊は右へ0.6セル動き、合計は2のままで、縁がぼけました。
ずっと大きなステップ$uh=2.6$で、もっと長い列$(0,0,1,1,0,0,0,0,0)$にすると、結果は$(0,0,0,0,0.4,1,0.6,0,0)$で、矢張りまともです。
新しい値は全て古い値の補間なので、古い値の範囲から決して出られません。
此の段階が爆発出来ないのは其の為です。
素朴な陽的差分法（「風上差分」、$s_i\leftarrow s_i-uh(s_i-s_{i-1})$）は、$uh=0.6$では同じ$(0,0,0.4,1,0.6,0)$を与えますが、$uh=2.6$では$(0,0,-1.6,1,2.6,0,\dots )$になります。
負の煙と、増える尖りです。
陽的な手法には$uh/\mathrm{\Delta }x\le 1$、CFL条件（Courant、Friedrichs、Lewyにちなむ）が必要です。
セミラグランジュ法の代償は数値散逸です。
ぼかしが毎ステップ繰り返され、渦は消え、煙はどろっと重たく見えます。

GPUでは、移流はセルごとに1回のテクスチャの読み出しで、双線型補間はサンプラーが行います（コンピュートシェーダーはGLSLシリーズの第15章で出て来ます）：

Texture2D<float2> velocity;
Texture2D<float4> source;
SamplerState linearClamp;

// Semi-Lagrangian advection: fetch the value from where the fluid came from.
// uv: this cell's texture coordinate, h: time step, u is stored in cells per second.
float4 advect(float2 uv, float h, float2 invGridSize) {
  float2 u = velocity.SampleLevel(linearClamp, uv, 0);
  float2 back = uv - h * u * invGridSize;
  return source.SampleLevel(linearClamp, back, 0);
}

layout(set = 0, binding = 0) uniform sampler2D velocity;
layout(set = 0, binding = 1) uniform sampler2D source;

// Semi-Lagrangian advection: fetch the value from where the fluid came from.
// uv: this cell's texture coordinate, h: time step, u is stored in cells per second.
vec4 advect(vec2 uv, float h, vec2 invGridSize) {
  vec2 u = textureLod(velocity, uv, 0.0).xy;
  vec2 back = uv - h * u * invGridSize;
  return textureLod(source, back, 0.0);
}

Texture2D<float2> velocity;
Texture2D<float4> source;
SamplerState linearClamp;

// Semi-Lagrangian advection: fetch the value from where the fluid came from.
// uv: this cell's texture coordinate, h: time step, u is stored in cells per second.
float4 advect(float2 uv, float h, float2 invGridSize) {
  float2 u = velocity.SampleLevel(linearClamp, uv, 0);
  float2 back = uv - h * u * invGridSize;
  return source.SampleLevel(linearClamp, back, 0);
}

2×2の格子での圧力投影

格子のソルバーは普通、速度をセルの面に、圧力をセルの中心に置きます。
此れがMAC格子（スタッガード格子）で、セルの発散と面を挟む圧力の差が自然に出て来ます。
格子の間隔1、$\rho =1$、$h=1$とすると：

• セルの発散は（右の面のu − 左の面のu）＋（上の面のv − 下の面のv）です。
• 射影は各面を挟む圧力の差を引きます：$u_{\text{face}}\leftarrow u_{\text{face}}-(p_{\text{right}}-p_{\text{left}})$、vも上と下で同様です。
• 結果の発散が0になるには、各セルで圧力が${\sum }_{\text{fluid neighbors}}(p_{\text{nb}}-p)=\text{divergence}$を満たさなければなりません。壁の面は速度0で、項を加えません。

計算例。
固い壁の箱の中に4つのセルがあり、下の2つのセルの間の面だけが速度$u=1$を持ちます（流体が右へ押されています）。
左下のセルは流体を失い（発散+1）、右下は得て（−1）、上のセルは0です。
各セルには流体の隣が2つあるので、1つのセルの更新は$p=\frac{\text{sum of neighbor pressures}-\text{divergence}}{2}$です。
全て0から始め、左下、右下、左上、右上の順でガウス＝ザイデル法を1周します：

• $p_{00}=\frac{0+0-1}{2}=-0.5$
• $p_{10}=\frac{p_{00}+p_{11}-(-1)}{2}=\frac{-0.5+0+1}{2}=0.25$
• $p_{01}=\frac{p_{00}+p_{11}-0}{2}=-0.25$
• $p_{11}=\frac{p_{10}+p_{01}-0}{2}=0$

2周目は何も変えないので、此の小さな系は解けました（圧力は定数の分だけ不定で、差だけが意味を持ちます）。
新しい面の速度は、下が$1-(0.25-(-0.5))=0.25$、上が$0-(0-(-0.25))=-0.25$、左の列が$0-(-0.25-(-0.5))=-0.25$、右の列が$0-(0-0.25)=0.25$です。
左下のセルを確かめます。
右の面から0.25出て、上の面から−0.25出る（つまり流入）ので発散は0で、他の3つのセルも同じ様に釣り合います。

右への一押しが、反時計回りの小さな渦になりました。
閉じた箱の中では、流体は回るしかないのです。
元の速さは4分の1しか残りません。
射影は発散の無い部分を残し、残りを捨てるからです。
本物の格子は${10}^5$〜${10}^8$個のセルを持ち、ポアソン方程式を解く部分（前処理付き共役勾配法、マルチグリッド、GPUでは決まった回数のヤコビ反復）が計算時間の大半を占めます。
此の格子では、素朴なヤコビ法（全てのセルを古い値から一斉に更新）は実は2つの状態の間を永遠に行き来します。
各回で3分の2だけ動く重み付きヤコビ法は収束します（下の面は0.333、0.333、0.259、0.259、0.251、…と0.25に近づきます）。

煙と炎

煙は流れと一緒に移流される2つの場、煙の密度$s$と温度$T$、そして浮力を加えます（Fedkiw・Stam・Jensen、「Visual Simulation of Smoke」、SIGGRAPH 2001）：

重い煙は沈み、熱いガスは昇ります。
同じ論文は渦度閉じ込めも導入しました。
渦（渦度、$𝐮$の回転）を測り、小さな渦を回し直す力を加えて、数値散逸が消した細部を取り戻すのです。
炎は燃料又は反応の場、燃料が燃える所の膨張、そして温度から決まる色（第3章の黒体の色）を加えます。

上の結果は玩具版（72×112セル、セミラグランジュ法の移流、浮力、少しの渦度閉じ込め、1ステップに60回のヤコビ反復）で、密度を明るさとして表示しています。
煙の中で光が散乱する様子をレンダリングするのはボリュームレンダリングで、第28章で扱います。

SPH：カーネルと密度

SPH（14.7節）では、ある点の全ての量は近くの粒子の重み付き和で、平滑化カーネル$W(r,h)$を使います。
カーネルは距離0で最大で、影響半径$h$で0になります（此処の$h$は半径で、時間刻みではありません）。
Müller・Charypar・Gross（2003年）は、仕事ごとに別のカーネルを選びました：

どちらも半径$h$の球の上で積分すると1です。
poly6は$r^2$しか要らない（平方根が要らない）ので密度に使います。
しかし其の傾きは中心で0なので、とても近い2つの粒子は圧力の力を感じず、固まってしまいます。
其の為、圧力の力には中心で傾きが最大のspikyカーネルを使います。
（粘性には、ラプラシアンが正の3つ目のカーネルを使います。）

計算例：密度。
$h=1$（カーネル半径を単位として）では$W_{\text{poly6}}(r)=\frac{315}{64\pi }(1-r^2{)}^3$、$\frac{315}{64\pi }=1.5667$です。
ある粒子が距離0.5、0.5、0.8に隣を持ち、質量は全て1です。
自分自身（$r=0$）も数えます：

状態方程式$p=k(\rho -{\rho }_0)$が密度を圧力に変えます。
$k=2$、静止密度${\rho }_0=2$なら$p=1.923$で、此の混み合った粒子は隣を押し退けます。
圧力の力は${𝐟}_i=-{\sum }_j{m}_j\frac{p_i+p_j}{2{\rho }_j}\mathrm{\nabla }{W}_{\text{spiky}}({𝐱}_i-{𝐱}_j)$で、2つの粒子が互いに等しく押し合う様に対称にしてあります。

此の「弱圧縮性」のSPHは、水が縮まない様に大きな$k$が必要で、其れは系を硬くします。
後の手法はソルバーで非圧縮性を課し（PCISPH、IISPH、DFSPH）、Position Based Fluids（Macklin・Müller、SIGGRAPH 2013）は${\rho }_i/{\rho }_0-1=0$を粒子ごとのPBDの拘束にします。
執筆時点のPhysXのGPU流体は此の方式です。
隣は、セルの大きさ$h$の均一格子で探します。

PIC、FLIP、APIC

格子は圧力が得意で移流が苦手、粒子は其の逆なので、ハイブリッド法は毎ステップ両方を使います：

1. P2G：粒子の速度を格子へ移します（近くの面へのカーネルで重み付けした平均）。
2. 格子の上で力と圧力投影を解きます。
3. G2P：速度を粒子へ戻し、粒子を格子の速度場の中で動かします。

問題は、どう戻すかです。
PIC（particle-in-cell）は各粒子に補間した新しい格子の速度を与えます。
格子より細かい物は全て平均されて消えるので、滑らかで安定ですが、散逸がとても大きいです。
FLIP（fluid-implicit-particle、Brackbill・Ruppel、1986年）は、格子の速度の変化だけを粒子自身の速度に足します。
粒子の細部が残りますが、ノイズも残ります。
ZhuとBridson（「Animating Sand as a Fluid」、SIGGRAPH 2005）は両方を非圧縮の流れに合わせてグラフィックスに持ち込み、数値粘性を調整する為の重み付き平均を提案しました（砂にはPIC、水にはFLIP）。

計算例。
粒子が2.4 m/sを持っています。
其の位置での格子の速度は、P2Gの後が2.0 m/s、圧力を解いた後が1.5 m/sです。
PIC：粒子は1.5になります。
FLIP：$2.4+(1.5-2.0)=1.9$。
例えばFLIP 95%、PIC 5%の混合：$0.95\cdot 1.9+0.05\cdot 1.5=1.88$。

APIC（the affine particle-in-cell method、Jiang・Schroeder・Selle・Teran・Stomakhin、SIGGRAPH 2015）は、各粒子の周りで速度がどう変わるかを表す小さな行列（局所的にアフィンな速度）を粒子に持たせ、両方向に移します。
回転が保たれ、移す時に角運動量が正確に保存され、FLIPのノイズ無しにPICの散逸が取り除かれます。

MPM（物質点法）は1990年代の計算力学から来た手法で、ハイブリッドの考えを固体に広げます。
粒子は質量、速度、そして26.3節のFEMの様な変形勾配$𝐅$を持ち、格子が毎ステップ応力の力を計算します。
メッシュを作る事も作り直す事も無いので、材料は裂け、流れ、合わさり、積み重なれます。
Stomakhin・Schroeder・Chai・Teran・Selle（「A Material Point Method for Snow Simulation」、SIGGRAPH 2013）はディズニーの「アナと雪の女王」（Frozen）の雪に此れを使い、今では砂、泥、泡など、半分固体で半分流体の材料に使われています。

水のレンダリング

シミュレーションが与えるのは粒子か格子だけで、絵には表面と材質が必要です：

• 表面：粒子から作った陰関数をマーチングキューブ法（第9章の9.3節）でポリゴン化するか、格子の上のレベルセット（符号付き距離場、第21章）を使います。ゲームでは粒子を深度バッファに描いて、画面空間でシェーディングする事も多いです。
• シェーディング：フレネル反射（第22章）、屈折、深い水に色を付ける吸収（第28章）、コースティクス、そして追加の粒子としてのしぶき、泡、気泡。
• 海は普通シミュレーションしません。測定されたスペクトルを持つ多数の波の和をFFTで計算すると、タイル状に繰り返せる高さ場が得られます。

26.5 髪、破壊、実際のシミュレーション

髪

14.4節では、髪の1本を伸びと曲げのバネを持つ粒子の鎖とし、ガイドヘアと補間で残りを作るとしました。
髪は殆ど伸びず（硬いバネには最悪の場合です）、ねじれます。
2つの方向が主流です：

• 位置ベースの鎖。 「Follow the leader」（Müller・Kim・Chentanez、VRIPHYS 2012）は、根元から先端への1回のパスで、各粒子を1つ前の粒子から正しい距離へ射影します。髪は正確に伸びなくなり、生じる偏りは速度の補正で隠します。
• ロッド。 映画の髪は、ねじれを含む弾性ロッドのモデルを陰的に解きます。

ボロノイによる破壊

14.7節では、破壊を、壊れる拘束で接着した事前に割った破片として説明しました。
破片は普通ボロノイ図から作ります。
種の点をばらまき（衝撃が予想される点の近くに多く撒くと、其処では細かい破片、遠くでは大きな塊になります）、物体の各点を一番近い種に割り当てます。
各ボロノイ領域は、他の種ごとに1つの半空間の共通部分で、其の境界は2つの種の垂直二等分線です。
其の為、全ての破片は凸で、凸の衝突形状やGJKにぴったりです。

計算例。
種${𝐬}_1=(0,0)$と${𝐬}_2=(4,2)$。
両方から等距離の点は$x^2+y^2=(x-4{)}^2+(y-2{)}^2$を満たし、整理すると$8x+4y=20$、つまり$2x+y=5$です。
点$(1,1)$はなので${𝐬}_1$に属します（確認：距離の2乗は2と$9+1=10$）。
領域を切り出すには、物体の多角形（又はメッシュ）を各二等分線で、第6章の半平面のクリッピングと全く同じ様に切ります。

GPUでの物理

GPUには何千もの独立した仕事が必要です（第29章）。
パーティクル、格子の流体（移流はセルごとに1回の読み出し、ヤコビ法はセルごとに1回のステンシル）、SPH、PBDは上手く対応します。
分岐の多いナローフェーズと逐次的な拘束の解法を持つ剛体は、もっと難しいです。
此の章のガウス＝ザイデル法のソルバーは本質的に逐次的なので（各拘束は隣の最新の値を使います）、GPUは2つの答えを使います：

• ヤコビ法：各拘束は古い位置から補正を計算し、補正を粒子ごとに平均します。完全に並列ですが、収束は遅くなります。
• グラフ彩色：拘束を、同じ粒子を共有する拘束が無いグループ（色）に分け、各グループを並列に、グループ同士を順に実行します。此れは正確なガウス＝ザイデル法のままです。26.3節の布は8色を使いました（水平、垂直、対角、逆対角の拘束を、其々列の偶奇（垂直の物は行の偶奇）で2つに分けた物）。

執筆時点では、NVIDIA PhysXの様なエンジンが剛体、PBDの粒子（流体、粒状体）、FEMのソフトボディをCUDAのGPUで動かし、多くのゲームが布、髪、パーティクルの効果をコンピュートシェーダーで動かしています。

決定性とネットワーク

同じ入力がビット単位で同じ結果を与える時、シミュレーションは決定的です。
当たり前に聞こえますが、浮動小数点の演算は結合的ではありません。
単精度では$({10}^8-{10}^8)+1=1$ですが、${10}^8+(-{10}^8+1)=0$です。
$-{10}^8+1$が丸められて$-{10}^8$に戻るからです。
足す順番の違い（スレッドの順序、SIMDの幅）、掛け算と足し算を1つのFMA命令にまとめるコンパイラ、異なるsinは最後のビットを変える事があり、カオス的な系（崩れる山、壊れる壁）は其れを数秒で全く違う結果に増幅します。

1台のマシンでの決定性には、固定時間刻み（14.8節）、物体と接触の決まった順番（ソートし、ハッシュ表の順番に頼らない）、初期化されていないデータが無い事が必要です。
プラットフォームを跨ぐ決定性には、更にコンパイラの浮動小数点の設定の管理が必要です。
執筆時点では、例えばJolt Physicsにはプラットフォームを跨いで同じ結果を出す事を目指すビルドオプションが有ります。

マルチプレイヤーのゲームは2つの戦略から選びます：

• ロックステップ：全てのマシンがシミュレーション全体を動かし、入力だけを送ります。帯域は僅か（何千ものユニットを持つストラテジーゲームで人気）ですが、完全な決定性が必要です。
• 状態同期：サーバーが正となるシミュレーションを動かし、物体の状態を送ります。クライアントは自分のプレイヤーを予測し（クライアント側予測）、サーバーと食い違ったら巻き戻して入力を再生します（ロールバック）。帯域は多く要りますが、小さな違いは問題になりません。

検定のポイント

以下はシラバス項目物理ベースアニメーション（第14章）に繋がり、検定対策として押さえておきたい区別です：

• 陽的オイラー法は現在の状態を使い、減衰の無いバネにはどんな時間刻みでも不安定です。陰的（後退）オイラー法はステップの終わりを使い、方程式を解く必要があり、どんなステップでも安定ですが人工的に減衰します。半陰的オイラー法はの間安定です。
• 衝突検出はブロードフェーズの後にナローフェーズです。2つの凸形状が離れているのは、射影を分ける軸がある時、そして其の時に限ります（SAT）。応答は接触法線に沿った力積で、摩擦は法線の力積の$\mu$倍までに制限されます。
• 硬いバネのバネ質点系の布には、小さなステップか陰的積分が必要です。位置ベースの手法は、拘束を満たす様に位置を直接動かします。
• 流体：オイラー的な格子法（Stable Fluids：移流、拡散、圧力投影）とラグランジュ的な粒子法（SPH）、そして両方を使うハイブリッド法。

26.6 まとめ

• 陰的オイラー法はバネのエネルギーを1ステップごとに$1+h^2{\omega }^2$で割ります。どんなステップでも安定（A安定）ですが減衰します。陽的オイラー法の安定領域$|1+z|\le 1$は虚軸を含みません。
• 剛体にはワールドの慣性テンソル$𝐑{𝐈}_{\text{body}}{𝐑}^{\mathrm{\top }}$、正規化したクォータニオンの更新、ジャイロ項への注意が必要です。
• SATは分離、法線、深さを与えます（2次元の箱で4軸、3次元で15軸）。クリッピングで接触点の集合を作ります。GJKとEPAはサポート関数を通じてどんな凸形状でも扱えます。
• シーケンシャルインパルス（射影ガウス＝ザイデル法）は累積インパルスをクランプし（${\lambda }_n\ge 0$、$|{\lambda }_t|\le \mu {\lambda }_n$）、ウォームスタート、スロップ付きのバウムガルテのバイアス、スリープに頼ります。ジョイントとモーターも同じ周回で解きます。
• PBDは位置を拘束へ射影します。XPBDのコンプライアンス$\tilde{\alpha }=\alpha /h^2$は硬さを反復回数から独立にします。ひずみ制限、シェイプマッチング、共回転FEMは同じ考えの上に立ちます。
• Stable Fluids：力、セミラグランジュ法の移流、陰的な拡散、圧力投影。SPH：密度にpoly6、圧力にspiky。PIC、FLIP、APICは粒子と格子の間で速度を移し、MPMは雪や砂を扱います。
• 髪は伸びない鎖かロッド、破壊は凸のボロノイの破片、GPUはヤコビ法かグラフ彩色を使います。決定性には固定時間刻み、決まった順番、管理された浮動小数点が必要です。

26.7 用語

26.8 課題

1. ★ 質量$m=1$ kgの物体が$k=400$ N/mのバネに吊るされ、$x_0=0.1$ m、$v_0=0$です。$h=0.05$秒の陰的オイラー法で2ステップ計算し、各ステップ後のエネルギーを求めて下さい。各ステップはエネルギーを何倍にしますか？ 陽的オイラー法なら何倍ですか？
2. ★ 6 kgの箱の辺の長さは、局所座標のx、y、z軸に沿って1、2、3 mです。${𝐈}_{\text{body}}$を求め、次にz軸周りに90°回転した後の${𝐈}_{\text{world}}$と、$𝐋=(1,1,0)$に対する角速度を求めて下さい。其れは$𝐋$と平行ですか？
3. ★★ 箱Aは原点にあって軸に平行で、半分の大きさは$(1.5,1)$です。箱Bは半分の大きさ$(1,0.5)$、軸${𝐮}_B=(0.6,0.8)$、${𝐯}_B=(-0.8,0.6)$で、中心は$(2.3,0.5)$です。4本の軸全てでSATを行い、接触法線と深さを求め、次にBを$(2.8,0.5)$に動かして、どの軸が2つを分けるかを答えて下さい。
4. ★★ 3つの粒子のロープが真っ直ぐ下に垂れています。${𝐩}_1=(0,0)$は固定、${𝐩}_2=(0,-1.5)$、${𝐩}_3=(0,-3)$、他は全て$w=1$、自然長は1です。PBDで2周（拘束1–2、次に2–3）行い、各周の後の長さを求めて下さい。次に両方の拘束を$\tilde{\alpha }=1$のXPBDで1周行って下さい。どちらが自然長に近いですか？ 其れが欠陥ではないのは何故ですか？
5. ★★ 26.4節の2×2の箱で、代わりに上の2つのセルの間の面だけが$u=-1$（左向き）で動く状態から始めます。発散、圧力のガウス＝ザイデル法1周（同じ順番）、4つの新しい面の速度を求めて下さい。全てのセルの発散が0である事を確かめ、渦がどちら向きに回るかを答えて下さい。
6. ★★★ 10個の箱の塔がガタガタ震え、ゆっくり地面に沈んで行きます。ソルバーの反復回数、ウォームスタート、バウムガルテの係数とスロップ、接触の順番、スリープの役割を説明し、最初にどれを変えるか、どんな副作用を予想するかを答えて下さい。

選択問題

問1 ★ 次のうち、硬いバネの積分について正しい物はどれか。

問2 ★★ 次のうち、剛体の接触について正しい物はどれか。

問3 ★★ 次のうち、シミュレーションの手法について正しい物はどれか。

26.9 実践課題

木箱、ロープ、煙の一吹きを手で。
方眼紙、鉛筆、定規、電卓を用意して下さい。
此の章の各種類のソルバーを1ステップずつ行います。

1. 陽的と陰的。 $k=100$ N/m、$m=1$ kgのバネが$x_0=0.2$ m、$v_0=0$から始まります。$h=0.1$秒で、陽的オイラー法と陰的オイラー法（26.1節の閉じた形）で1ステップずつ進め、其々の前後のエネルギーを求めて下さい。
2. SAT。 方眼紙に箱A（中心$(0,0)$、軸に平行、半分の大きさ$(1.5,1)$）と箱B（中心$(1,1.8)$、半分の大きさ$(1,0.5)$、軸${𝐮}_B=(0.6,0.8)$、${𝐯}_B=(-0.8,0.6)$）を描きます。4本の軸全てで重なりを計算し、接触法線と深さを選び、図に法線を描いて下さい。
3. 摩擦付きの力積。 26.2節の木箱（$m=1$、$I=\frac{1}{6}$、$𝐯=(0,-2)$、$𝐫=(0.5,-0.5)$、$𝐧=(0,1)$）を$\mu =0.2$でやり直します。シーケンシャルインパルスの反復を手で2回（法線、次に摩擦、累積の値をクランプ）行い、次に本文と同じ様に滑る場合を厳密に解いて下さい。
4. PBDのロープ。 3つの粒子が$(0,0)$（固定）、$(0,-1)$、$(0,-2)$に静止していて、自然長は1です。$g=10$ m/s²、$h=0.1$秒で位置を予測し、距離の射影を2周（1–2、次に2–3）行い、位置の変化から新しい速度を求めて下さい。
5. 圧力投影。 26.4節の2×2の箱で、左の2つのセルの間の面が$v=1$（上向き）、他は全て0です。発散、圧力のガウス＝ザイデル法1周、新しい面の速度を求め、結果を矢印で描いて下さい。
6. SPHの密度。 $h=1$、単位質量で、$(0,0)$、$(0.5,0)$、$(0,0.5)$、$(0.6,0.8)$の粒子に対する$(0,0)$でのpoly6の密度を求めて下さい。

自分で確認する方法：

• 段階1：エネルギーは2 Jから始まります。陽的オイラー法は其れを2倍に、陰的オイラー法は半分にし（$h\omega =1$）、陰的な結果は$x_1=x_0+h{v}_1$を満たします。
• 段階2：4つの重なりは全て正です。最小はAのy軸上なので、法線は$(0,1)$、深さは0.5未満です。
• 段階3：反復1の後、摩擦の力積は丁度限界$-\mu {\lambda }_n$にあり、残りの$v_n$は反復ごとに縮み、厳密な${\lambda }_n$は0.8と1の間で、接触は滑り続けます。
• 段階4：各周の後、2つ目の長さは丁度1で、1つ目は少し長過ぎ、2周目の後は其の半分になります。動く2つの粒子は同じ小さな下向きの速度を得ます。
• 段階5：全ての面の速さは0.25で、4つのセルの発散は全て0、矢印は輪を作ります（今回は時計回りです）。
• 段階6：$(0.6,0.8)$の粒子は丁度距離$h$にあり、何も寄与しません。密度は2.9より少し小さくなります。