2009年5月30日土曜日

粒子法のプログラム第5回(力の計算)


今回は、粒子法のプログラムでどのように力を計算しているかを説明します。

力の計算は次のコードになります。


void calculate_force( Particles* ps )
{
double pterm, vterm, r, c;
vec dr(3), force(3), fcurr(3);

FOR_EACH_PARTICLES( ps, p )
{
force(0) = 0.0; force(1) = 0.0; force(2) = 0.0;
FOR_EACH_PARTICLES( ps, pj )
{
if ( p == pj ) continue;
dr = ( p->r - pj->r ) * SPH_SIMSCALE;
r = norm( dr );
if ( H > r )
{
c = H - r;
pterm = -0.5 * c * SpikyKern * ( p->p + pj->p ) / r;
vterm = LapKern * SPH_VISC;
fcurr = pterm * dr + vterm * ( pj->v - p->v );
fcurr *= c * p->rho * pj->rho;
force += fcurr;
}
}
p->f = force;
}
}
それぞれの粒子について、ナビエ=ストークス方程式の圧力項、粘性項にあたる力を計算しています。

影響半径内にある自分以外の粒子に対して、距離に応じて重みを足し合わせています。SpikyKernとLapKernがカーネルの係数です。

力の計算も、密度と圧力の計算と同様に、近傍粒子探索を工夫したり対称性を利用することで、高速化できます。

次回は、どのように粒子の位置を更新するかという説明になります。

粒子法のプログラム
粒子法のプログラム第1回(概要)
粒子法のプログラム第2回(プログラムの大枠)
粒子法のプログラム第3回(データ構造)
粒子法のプログラム第4回(密度と圧力の計算)
粒子法のプログラム第5回(力の計算)
粒子法のプログラム第6回(境界条件と粒子位置の更新)
粒子法のプログラム最終回(粒子の出力)

その他の解説エントリ
SPHによる巻き波のシミュレーション1
SPHによる巻き波のシミュレーション2
SPHによる巻き波のシミュレーション3
このあとやりたいこと

固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション

流体シミュレーションに関するエントリ
【粒子法】粒子を流体としてレンダリング
3次元の粒子法シミュレーション
粒子法のシーンを2倍のサイズにしてみたが…
粒子法のシーンを2倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation

剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション(プレビュー)
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies

動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧


--

2009年5月24日日曜日

目覚めを良くするための最良の方法


最近、枕を変えてみた。

これまではふつうのやわらかい枕。それを、四角くて固めのい草枕に変えた。

すると、目覚めがすごくすっきりするようになった。

枕でここまで変わるのかー

びっくり。これから蒸し暑くなるので、その点でもグッド!


2009年5月23日土曜日

粒子法のプログラム第4回(密度と圧力の計算)


今回は、粒子法のプログラムでどのように密度と圧力を計算しているかを説明します。

粒子の密度と圧力の計算は次のコードになります。
void calculate_density_and_pressure( Particles* ps )
{
double r2, c, sum, H2;
vec dr;

H2 = H*H;

FOR_EACH_PARTICLES( ps, p )
{
sum = 0.0;
FOR_EACH_PARTICLES( ps, pj )
{
if ( p == pj ) continue;
dr = (p->r - pj->r) * SPH_SIMSCALE;
r2 = norm2( dr );
if ( H2 > r2 )
{
c = H2 - r2;
sum += c * c * c;
}
}
p->rho = sum * SPH_PMASS * Poly6Kern;
p->p   = ( p->rho - SPH_RESTDENSITY ) * SPH_INTSTIFF;
p->rho = 1.0 / p->rho;  // take inverse for later calculation
}
}

それぞれの粒子について、影響半径内にある自分以外の粒子に対して、距離に応じて重みを足し合わせています。そのあと、足し合わせた重みに、粒子質量についての係数とカーネル係数をかけています。密度が計算された後、圧力を密度から計算します。

上のやり方だと、計算量のオーダーは、各粒子について他の粒子それぞれとの計算になるので、粒子の数をnとするとO(n^2)となります。しかし、影響半径外の粒子との計算は不要なので、他粒子すべてと計算する必要は本来ありません。なので、近傍粒子をあらかじめマップに入れておき、それを参照するようにすれば、計算量をO(n)にすることができます。

また、粒子pと粒子pjについて対称なので、その対称性を利用すれば、計算量をさらに半分にすることが可能です。

密度と圧力は、このように計算されています。いったん素直な実装になっていますが、上の2つの点で高速化することができます。

次回は、粒子にかかる力をどのように計算しているかという説明です。

粒子法のプログラム
粒子法のプログラム第1回(概要)
粒子法のプログラム第2回(プログラムの大枠)
粒子法のプログラム第3回(データ構造)
粒子法のプログラム第4回(密度と圧力の計算)
粒子法のプログラム第5回(力の計算)
粒子法のプログラム第6回(境界条件と粒子位置の更新)
粒子法のプログラム最終回(粒子の出力)

その他の解説エントリ
SPHによる巻き波のシミュレーション1
SPHによる巻き波のシミュレーション2
SPHによる巻き波のシミュレーション3
このあとやりたいこと

固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション

流体シミュレーションに関するエントリ
【粒子法】粒子を流体としてレンダリング
3次元の粒子法シミュレーション
粒子法のシーンを2倍のサイズにしてみたが…
粒子法のシーンを2倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation

剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション(プレビュー)
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies

動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧


--

2009年5月19日火曜日

【粒子法】粒子を流体としてレンダリング


粒子法では流体を粒子として表現します。

そのため、そのまま描画すると、この投稿のように球のあつまりとして描画することになります。

これを、流体として見えるようにレンダリングしてみました。

マーチングキューブで陰曲面をポリゴン化し、レイトレーシングでレンダリングしています。

フォトンマップも使っているのでコースティクスも表現されています。

屈折率を高めにしたのでキラキラ。




■固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション

■関連する記事
粒子法のプログラム第1回(概要)
粒子法のプログラム第2回(プログラムの大枠)
粒子法のプログラム第3回(データ構造)
粒子法のプログラム第4回(密度と圧力の計算)
粒子法のプログラム第5回(力の計算)
粒子法のプログラム第6回(境界条件と粒子位置の更新)
粒子法のプログラム最終回(粒子の出力)
【粒子法】粒子を流体としてレンダリング
3次元の粒子法シミュレーション
粒子法のシーンを2倍のサイズにしてみたが…
粒子法のシーンを2倍のサイズにしてみた

粒子法(SPH)のプログラム
粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。

粒子法(SPH)のプログラム一覧



動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧


--