*位置ベース流体(Position Based Fluids) [#k29c9bfc]
粒子法を用いた流体シミュレーションにおいて重要な要素となるのが
非圧縮性の強制である．
非圧縮性とは流体の密度が流れによって空間的にも時間的にも変化しないことを表し，
水のような液体はもちろん，
流速が音速より十分小さいならば空気も非圧縮性流体と考えることができる
(ちなみに音速に近くなると空気は圧縮する，
音速に達した飛行機が発する衝撃波いわゆるソニックブームは空気が圧縮することで引き起こされる)．
非圧縮性のための様々な方法(MPS,WCSPH,PCISPH,IISPHなど)が提案されているが，
それらに比べて非常に高速かつ安定した計算が可能なのが位置ベース流体法(Position Based Fluids, 以下PBF)¬e{Macklin2013:M. Macklin and M. Muller, "Position based fluids" ACM Trans. Graph., 32, pp.104:1-104:12, 2013.};
である．
PBFは位置ベース法(Position-based method)¬e{Muller2007:M. Muller, B. Heidelberger, M. Hennix and J. Ratcliff, "Position based dynamics", J. Vis. Comun. Image Represent., 18, pp.109-118, 2007.};
の一種である．
位置ベース法では従来の力学ベースの方法と異なり，位置(物体の空間座標値)を満たすべき制約に従って直接変更する方法である
力学ベースでは支配方程式に従い加速度->速度->位置と変更していくため，
それぞれで積分に伴う誤差の蓄積が発生するが位置に直接制約をかける位置ベース法はそれがない．
ただし，頂点や粒子の位置関係に基づく制約をかけるためどちらかというと幾何学的な方法で，
あまり物理的ではないのではないかという議論もある．
**位置ベース法の基本 [#v82f65a4]
PBFの前に基本となる位置ベース法の考え方を説明する．
まず，移動する計算点&ref(pbf.eq1.gif,nolink,70%);に対して，その座標値そのものに制約条件&ref(pbf.eq2.gif,nolink,70%);を掛けることを考える．
&ref(pbf.eq3.gif,nolink,70%);の条件が常に満たされているとすると，
粒子が微少量&ref(pbf.eq4.gif,nolink,70%);移動した後もこの条件は満たされているということになるので，
&ref(pbf.eq5.gif,nolink,70%);となる．
これをテイラー展開してやれば，
#ref(pbf.eq6.gif,nolink,70%)
となる．2次以降の項を無視すれば&ref(pbf.eq4.gif,nolink,70%);に関する方程式が得られる．
#ref(pbf.eq7.gif,nolink,70%)
&ref(pbf.eq8.gif,nolink,70%);は制約として定義しているのでそこから計算できる．
いま，&ref(pbf.eq1.gif,nolink,70%);を粒子の中心座標とすると求めるべき未知数の数は
粒子数を&ref(pbf.eq9.gif,nolink,70%);としたら3次元空間の場合&ref(pbf.eq10.gif,nolink,70%);となる．
一方で&ref(pbf.eq11.gif,nolink,70%);が各粒子ごとに計算される(つまり&ref(pbf.eq12.gif,nolink,70%);)ならば，
上記の方程式の数は&ref(pbf.eq9.gif,nolink,70%);なので未知数&ref(pbf.eq10.gif,nolink,70%);に対して少なすぎて解が定まらない
(疑似逆行列などを使ってできないことはない)．
そこで&ref(pbf.eq4.gif,nolink,70%);に対して別の制約をかける．
計算点&ref(pbf.eq1.gif,nolink,70%);が運動方程式に従って移動していると考えると，
ここでの&ref(pbf.eq4.gif,nolink,70%);はその運動とは別の制約&ref(pbf.eq11.gif,nolink,70%);を満たすための移動となる．
そうすると元の運動(平行移動と回転)を阻害しないような移動方向，
言い換えると平行移動と回転運動量を保存するような移動方向を定める必要がある．
制約&ref(pbf.eq11.gif,nolink,70%);がこれらの運動の状態に全く作用されない(依存しない)と考えるならば，
その傾き&ref(pbf.eq13.gif,nolink,70%);はこれらの運動方向に対して垂直な方向となる．
つまり，&ref(pbf.eq4.gif,nolink,70%);の方向を&ref(pbf.eq13.gif,nolink,70%);となるような制約を掛けてやればよいことになる．
これはちょうど制約付きの最適化問題になるのでラグランジュ乗数&ref(pbf.eq14.gif,nolink,70%);を導入して，
#ref(pbf.eq15.gif,nolink,70%)
とする．
これを上の&ref(pbf.eq4.gif,nolink,70%);に関する方程式に代入して整理すると，
#ref(pbf.eq16.gif,nolink,70%)
となる．求めた&ref(pbf.eq14.gif,nolink,70%);を&ref(pbf.eq4.gif,nolink,70%);の式に代入すれば制約を満たすための移動ベクトルが分かるので，
それで計算点位置を修正してやればよい．
複数の計算点が絡む場合は1回の修正では制約を満たすことは難しいので，
上記の式による制約を条件を満たすまで反復処理する．
**位置ベース法の粒子法への適用 [#b2340bb0]
さて，話を粒子法の場合に戻そう．
まず，非圧縮性条件が満たすべき制約(Constraint)を考える．
流体が圧縮しないということはその密度&ref(pbf.eq17.gif,nolink,70%);が初期密度&ref(pbf.eq18.gif,nolink,70%);から変化しないということであり，
つまり，条件としては&ref(pbf.eq19.gif,nolink,70%);である．これを変形して位置&ref(pbf.eq1.gif,nolink,70%);に関する制約&ref(pbf.eq11.gif,nolink,70%);として定式化すると，&br;
&ref(pbf.eq20.gif,nolink,70%); ・・・(1)&br;
となる．
ここでは粒子法の一種であるSPH法を使ってある粒子&ref(pbf.eq21.gif,nolink,70%);の密度&ref(pbf.eq22.gif,nolink,70%);を求める．
#ref(pbf.eq23.gif,nolink,70%)
ここで&ref(pbf.eq24.gif,nolink,70%);は粒子&ref(pbf.eq21.gif,nolink,70%);の近傍粒子の集合，&ref(pbf.eq25.gif,nolink,70%);は粒子質量，&ref(pbf.eq26.gif,nolink,70%);は有効半径である．
近傍粒子数を&ref(pbf.eq9.gif,nolink,70%);とすると粒子&ref(pbf.eq21.gif,nolink,70%);の密度に関する制約は近傍粒子座標&ref(pbf.eq27.gif,nolink,70%);の
関数となるので，&ref(pbf.eq14.gif,nolink,70%);の計算式は以下のようになる．&br;
&ref(pbf.eq28.gif,nolink,70%); ・・・(2)&br;
&ref(pbf.eq29.gif,nolink,70%);は密度&ref(pbf.eq22.gif,nolink,70%);の式を式(1)に代入すれば計算できる．
&ref(pbf.eq30.gif,nolink,70%);も同様で，SPH法における勾配値は単に重み関数&ref(pbf.eq31.gif,nolink,70%);の勾配をとればよいので，
#ref(pbf.eq32.gif,nolink,70%)
となる．
この式で注意しないといけないのは&ref(pbf.eq31.gif,nolink,70%);は&ref(pbf.eq33.gif,nolink,70%);と&ref(pbf.eq34.gif,nolink,70%);(と&ref(pbf.eq26.gif,nolink,70%);)に関する関数なので，
&ref(pbf.eq35.gif,nolink,70%);で微分したときに0にならないのは&ref(pbf.eq36.gif,nolink,70%);と&ref(pbf.eq37.gif,nolink,70%);の場合だけである．
そうすると上記の式はさらに書き下すことができて，&br;
&ref(pbf.eq38.gif,nolink,70%); ・・・(3)&br;
となる．&ref(pbf.eq36.gif,nolink,70%);の時は&ref(pbf.eq39.gif,nolink,70%);に関わらず&ref(pbf.eq40.gif,nolink,70%);は&ref(pbf.eq35.gif,nolink,70%);を含むのでそのまま計算すればよい．
&ref(pbf.eq41.gif,nolink,70%);の場合は&ref(pbf.eq42.gif,nolink,70%);となったときだけ&ref(pbf.eq40.gif,nolink,70%);を含むのでそれだけ計算する
(マイナスの符号がついているのは&ref(pbf.eq43.gif,nolink,70%);なので)．
また，粒子質量&ref(pbf.eq25.gif,nolink,70%);がすべての粒子で一定ならば&ref(pbf.eq44.gif,nolink,70%);の外に出せて，
最終的に&ref(pbf.eq4.gif,nolink,70%);の計算ではすべて打ち消されるので計算上は省略してもよい．
式(2)と式(3)により各粒子のラグランジュ乗数&ref(pbf.eq45.gif,nolink,70%);が計算されるので，
これらから各粒子の位置修正量&ref(pbf.eq46.gif,nolink,70%);を求める．
このとき近傍粒子への影響と近傍粒子からの影響が非対称にならないようする(SPH法の圧力項等と同じ)．&br;
&ref(pbf.eq47.gif,nolink,70%); ・・・(4)&br;
先ほども説明したように質量一定ならばこちらの&ref(pbf.eq48.gif,nolink,70%);も省略できる．
&ref(pbf.eq45.gif,nolink,70%);と&ref(pbf.eq49.gif,nolink,70%);の平均をとるならば&ref(pbf.eq50.gif,nolink,70%);をつけるところだが，
SOR法の加速緩和係数と同じ役割を果たすようで，ない方が収束は速くなる．
**張力の安定化 [#za0cd3e4]
SPH法では近傍粒子数が少ない場合に負の圧力が発生して粒子同士が過度に集まってしまう
(particle clusteringやparticle stackingと呼ばれている)現象が発生する．
PBFではこれを解決するためにMonaghan
¬e{Monaghan2000:J. Monaghan, "SPH without a Tensile Instability", Journal of Computational Physics, 159, pp.290-311, 2000.};
が提案した以下のような項を式(4)に追加している．
#ref(pbf.eq51.gif,nolink,70%)
式中のパラメータの値としては&ref(pbf.eq52.gif,nolink,70%);,&ref(pbf.eq53.gif,nolink,70%);,&ref(pbf.eq54.gif,nolink,70%);あたりがよいようである．
式(4)の&ref(pbf.eq55.gif,nolink,70%);の部分を&ref(pbf.eq56.gif,nolink,70%);とするのだが，
Monaghanの元の論文では&ref(pbf.eq57.gif,nolink,70%);は加速度なので位置修正に直接使うのはおかしくなるので，
下記の私の実装では&ref(pbf.eq58.gif,nolink,70%);を掛けてある．
**境界粒子 [#fb3d5633]
前述のparticle stackingの問題は固体境界付近でも発生する．
粒子と固体境界の単純な衝突処理だけだと固体内に粒子が存在しないため，
境界付近に粒子が集まってしまう．
それを解決する一つの方法が固体境界内に位置固定の粒子(境界粒子)を散布することである．
通常，境界粒子は有効半径に合わせて複数層配置するが，
Akinciらの方法¬e{Akinci2012:N. Akinci, M. Ihmsen, G. Akinci, B. Solenthaler and M. Teschner, "Versatile rigid-fluid coupling for incompressible SPH", ACM Trans. Graph., ACM, 31, pp.62:1-62:8, 2012.};
を使えば，仮想的な体積を事前に計算しておくことで1層だけで済ますことができる．
境界に沿って粒子を配置した後，境界粒子だけを使って以下の仮想体積を計算して格納しておく．
#ref(pbf.eq59.gif,nolink,70%)
そして，PBFでも計算時に境界粒子の質量を&ref(pbf.eq60.gif,nolink,70%);としてやればよい
(詳しい式などは[[こちらの論文:http://slis.tsukuba.ac.jp/~fujis/pdf/iieej2015_fujis.pdf]]を参照)．
*実装例 [#fc7280bb]
Visual Studio 2012 + CUDA7.5の環境で作成したコードは以下．
#ref(rx_pbf.zip)
-ビルドするのに必要なライブラリ
freeglut, GLEW, CUDA, (OpenVDB)
各ライブラリについては[[ライブラリのインストール]]を参照．
-簡単な説明＆注意事項
--GUIツールキットとしてFLTKを用いている．
--rx_fltk_glcanvas.hの55行目
#define RX_USE_GPU
をコメントアウトするとCPUで計算する．
--OpenVDBはポリゴンモデルを読み込んでその内部に境界粒子を配置する場合に必要．使用する場合はプロプロセッサ定義にRX_USE_OPENVDBを追加する．
--モデルの読み込み自体は[[3Dモデルファイルの入出力]]のrx_model.libを使っている．
--Visual Studioから実行する場合は，作業ディレクトリを"../../bin"にしておく．
--"s"キー or Start/Stopボタンでシミュレーションスタート/ストップ．
--F1〜F9キー or Sceneメニューからシーン切替可能．binフォルダに"sph_scene_*.cfg"という名前のファイルがあり，
これがシーンを記述するファイル．デフォルトでは4つのシーンが入っている．
--CUDAはv7.5で動作を確認している．
--ポリゴンモデルを固体境界とした場合に境界粒子配置のためにOpenVDBを用いている(正確には境界粒子を配置するための陰関数場生成のため)．デフォルトではOpenVDB使用をOFFにしているので，ポリゴン境界に粒子を配置したい場合はプリプロセッサでRX_USE_OPENVDBを定義するか，rx_sph_solid_poly.cppの先頭でdefineすること．
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
