SPH法の実装(GPU実装含む)
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
SPH法をCPU,GPUで実装した例を示す．
----
#contents
----
*近傍探索の効率化 [#ude709f6]
SPHなどの粒子法では近傍粒子からの影響を考慮して，支配方程...
近傍探索を効率化する手法として空間分割法がある．空間分割...
+物体の登録
++空間を分割
++各物体が所属する分割領域を計算
+近傍探索
++探索中心座標が所属する領域を計算
++自身と隣接する領域の物体をリストアップ
++リストアップされた物体との距離を算出
となり，大きく登録と探索の2つに分かれる．登録はシーン中の...
空間分割法で問題となるのはどのように空間を分割するかであ...
-等間隔格子(下図左)
-八分木,四分木構造(下図右)
-kD木構造
などがあげられる．下の二つは分割領域を階層構造にして，領...
一方，粒子法では各ステップごとに粒子がダイナミックに動く...
そのため，分割・登録が容易な等間隔格子がよく用いられる．
等間隔格子の場合は，階層構造を持つほかの領域分割法に比べ...
CENTER:&ref(particle_grid.gif); &ref(particle_octree.gif);
CENTER:図1 等間隔格子(左)と四分木(左)
*GPUでの実装 [#e2bd734c]
[[Particle Simulation using CUDA (PDF):http://developer.d...
を参考にしてSPH法へ拡張する．
最初にParticle Simulation using CUDAで用いられている近傍...
以下の図を例として説明する．
CENTER:&ref(particle_grid_2.gif);
CENTER:図2 グリッドセルとパーティクル
ステップの最初にパーティクルの登録作業を行う．登録手順は...
+各パーティクルのグリッドハッシュ値を計算．ハッシュ値は各...
hash = i+j・nx
とする．(nx,ny)はグリッドセルの数で，図2において各領域の...
ハッシュ値はGridParticleHash配列に格納する．また，後でソ...
|GridParticleHash|2|0|0|1|1|0|1|2|
| SortedIndex|0|1|2|3|4|5|6|7|
+グリッドハッシュ値をキーとしてソートする．ソートには[[基...
|GridParticleHash|0|0|0|1|1|1|2|2|
| SortedIndex|1|2|5|3|4|6|0|7|
+GridParticleHashで同じ値(同じ分割領域)が続く部分の始まり...
まず，配列(サイズ=グリッドセル数)を0xffffffffで初期化する．
|CellStart|0xff|0xff|0xff|0xff|
| CellEnd|0xff|0xff|0xff|0xff|
ソート済みのGridParticleHashの各要素iに対して一つ前の要素...
GridParticleHash[i] != GridParticleHash[i-1]
ならば，iを始点としてCellStart[GridParticleHash[i]]に，
i-1を一つ前のグリッドセルの終点としてCellEnd[GridParticle...
|CellStart|0|3|6|0xff|
| CellEnd|2|5|7|0xff|
このときパーティクル位置などを格納した配列をSortedIndexを...
SortedPos[i] = Pos[SortedIndex[i]];
近傍探索フェーズの手順は以下．
+パーティクル位置 (or 任意の座標値) x のグリッドセルを算出
+周囲のグリッドを含めて以下の処理を実行
++グリッドハッシュ値hash計算
++CellStart,CellEndを参照して(CellStart[hash],CellEnd[has...
++for(int j = start; j <= end; ++j) でループ
+++近傍候補パーティクル位置pjを参照
pj=SortedPos[j]
+++|pj-x|<hならそのパーティクルは近傍
SPHで用いる場合は，近傍探索フェースの最後で見つかった近傍...
また，グリッド幅は有効半径hに基づいて決定する．
各パーティクルの密度値を計算する例を以下に示す．
#code(C){{
/*!
* 与えられたセル内のパーティクルとの距離から密度を計算
* @param[in] gridPos グリッド位置
* @param[in] index パーティクルインデックス
* @param[in] pos 計算座標
* @param[in] cell パーティクルグリッドデータ
* @return セル内のパーティクルから計算した密度値
*/
__device__
float calDensityCell(int3 gridPos, uint i, float3 pos0, r...
{
uint gridHash = calcGridHash(gridPos);
// セル内のパーティクルのスタートインデックス
uint startIndex = cell.dCellStart[gridHash];
float h = params.EffectiveRadius;
float dens = 0.0f;
if(startIndex != 0xffffffff){ // セルが空でないかのチェ...
// セル内のパーティクルで反復
uint endIndex = cell.dCellEnd[gridHash];
for(uint j = startIndex; j < endIndex; ++j){
//if(j == i) continue;
float3 pos1 = make_float3(cell.dSortedPos[j]);
float3 rij = pos0-pos1;
float r = length(rij);
if(r <= h){
float q = h*h-r*r;
dens += params.Mass*params.Wpoly6*q*q*q;
}
}
}
return dens;
}
/*!
* パーティクル密度計算(カーネル関数)
* @param[out] newDens パーティクル密度
* @param[out] newPres パーティクル圧力
* @param[in] cell パーティクルグリッドデータ
*/
__global__
void sphCalDensity(float* newDens, float* newPres, rxPart...
{
// パーティクルインデックス
uint index = __mul24(blockIdx.x,blockDim.x)+threadIdx.x;
if(index >= cell.uNumParticles) return;
float3 pos = make_float3(cell.dSortedPos[index]); // パ...
float h = params.EffectiveRadius;
// パーティクル周囲のグリッド
int3 grid_pos0, grid_pos1;
grid_pos0 = calcGridPos(pos-make_float3(h));
grid_pos1 = calcGridPos(pos+make_float3(h));
// 周囲のグリッドも含めて近傍探索，密度計算
float dens = 0.0f;
for(int z = grid_pos0.z; z <= grid_pos1.z; ++z){
for(int y = grid_pos0.y; y <= grid_pos1.y; ++y){
for(int x = grid_pos0.x; x <= grid_pos1.x; ++x){
int3 n_grid_pos = make_int3(x, y, z);
dens += calDensityCell(n_grid_pos, index, pos, cell);
}
}
}
// 密度値を結果に書き込み
uint oIdx = cell.dSortedIndex[index];
newDens[oIdx] = dens;
}
}}
デバイス関数calDensityCellはセル内のパーティクルを参照し...
*表面メッシュ生成 [#ca575a8b]
パーティクルをそのままレンダリングしたのでは流体には見え...
特に透明な液体をレンダリングする際には，液体表面の位置と...
ここでは3DCGにおいて一般的に用いられている三角形メッシュ...
三角形メッシュによる液体表面抽出にはMC(Marching Cubes)法...
MC法では，陰関数で表現されたフィールドから等値面を抽出す...
手順としては，
+立方体グリッド(Cube)を空間に配置
+各Cubeの辺上で陰関数値が0となる点を線型補間や二分法，ニ...
+Cubeに含まれる点の配置と数から三角形メッシュを生成
グリッドに含まれるメッシュ数やメッシュの構成はテーブルを...
より詳しい手順やソースコード，テーブルは以下のサイトを参...
-http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/
GPUでの実装は[[CUDAでMC法(改造版)]]参照．
*実装結果 [#k4a9f2ef]
GeForceGTX580の環境でテストした．
パーティクル数が約25000で，SPH計算だけなら200fps，メッシ...
ちなみにCPUでも同様の方法で実装したところ，SPH計算だけで1...
(メッシュ化には上記のGPU実装MC法を用いたが，GPUへのデータ...
図3,4に実行結果を示す．
CENTER:&ref(sph_result1.jpg,,50%); &ref(sph_result2.jpg,...
CENTER:図3 左からパーティクル，表面メッシュ，屈折レンダリ...
CENTER:&ref(sph_result4.jpg,,50%);
CENTER:図4 ポリゴンで表された谷形状内での流れ
*ソースコード [#cf51a873]
***CUDA9.1用ソースコード [#rebfb11d]
Visual Studio 2015 + CUDA9.1の環境で動作確認したコードは...
#ref(rx_sph_v2.1.zip)
-必要ライブラリや説明は下記のCUDA5用のものとほぼ同じだが...
-動作確認済み環境 : Visual Studio 2015 (Version 14.0.2543...
***CUDA5用ソースコード [#y2c41c83]
Visual Studio 2010 + CUDA5.0の環境で作成したコードは以下．
#ref(rx_sph_v2.zip)
-ビルドするのに必要なライブラリ
freeglut, GLEW, CUDA, boost, libpng, zlib, libjpeg, ftgl...
各ライブラリについては[[ライブラリのインストール]]を参照．
-簡単な説明＆注意事項
--GUIツールキットとしてFLTKを用いている．
--rx_fltk_glcanvas.hの59行目
#define RX_USE_GPU
をコメントアウトするとCPUで計算する．
--[[Anisotropic Kernel]]によるメッシュ生成にも対応してい...
--Visual Studioから実行する場合は，作業ディレクトリを"../...
--CUDA5のSDKにあるヘッダファイル(helper_*)を用いているの...
C:\ProgramData\NVIDIA Corporation\CUDA Samples\v5.0\comm...
にインストールされる．
--"s"キー or Start/Stopボタンでシミュレーションスタート/...
--F1〜F9キー or Sceneメニューからシーン切替可能．binフォ...
これがシーンを記述するファイル．デフォルトでは7つのシーン...
--bin/cubeフォルダには何も入っていないが，屈折レンダリン...
サンプルとしてterragenで作成したものを以下に置く．
#ref(cubemap_texture_example.zip)
--[[OBJ形式]]のファイルを読み込んで，固体オブジェクトにす...
--研究用に作ったコードを削って軽くした物なので，いろいろ...
終了行:
SPH法をCPU,GPUで実装した例を示す．
----
#contents
----
*近傍探索の効率化 [#ude709f6]
SPHなどの粒子法では近傍粒子からの影響を考慮して，支配方程...
近傍探索を効率化する手法として空間分割法がある．空間分割...
+物体の登録
++空間を分割
++各物体が所属する分割領域を計算
+近傍探索
++探索中心座標が所属する領域を計算
++自身と隣接する領域の物体をリストアップ
++リストアップされた物体との距離を算出
となり，大きく登録と探索の2つに分かれる．登録はシーン中の...
空間分割法で問題となるのはどのように空間を分割するかであ...
-等間隔格子(下図左)
-八分木,四分木構造(下図右)
-kD木構造
などがあげられる．下の二つは分割領域を階層構造にして，領...
一方，粒子法では各ステップごとに粒子がダイナミックに動く...
そのため，分割・登録が容易な等間隔格子がよく用いられる．
等間隔格子の場合は，階層構造を持つほかの領域分割法に比べ...
CENTER:&ref(particle_grid.gif); &ref(particle_octree.gif);
CENTER:図1 等間隔格子(左)と四分木(左)
*GPUでの実装 [#e2bd734c]
[[Particle Simulation using CUDA (PDF):http://developer.d...
を参考にしてSPH法へ拡張する．
最初にParticle Simulation using CUDAで用いられている近傍...
以下の図を例として説明する．
CENTER:&ref(particle_grid_2.gif);
CENTER:図2 グリッドセルとパーティクル
ステップの最初にパーティクルの登録作業を行う．登録手順は...
+各パーティクルのグリッドハッシュ値を計算．ハッシュ値は各...
hash = i+j・nx
とする．(nx,ny)はグリッドセルの数で，図2において各領域の...
ハッシュ値はGridParticleHash配列に格納する．また，後でソ...
|GridParticleHash|2|0|0|1|1|0|1|2|
| SortedIndex|0|1|2|3|4|5|6|7|
+グリッドハッシュ値をキーとしてソートする．ソートには[[基...
|GridParticleHash|0|0|0|1|1|1|2|2|
| SortedIndex|1|2|5|3|4|6|0|7|
+GridParticleHashで同じ値(同じ分割領域)が続く部分の始まり...
まず，配列(サイズ=グリッドセル数)を0xffffffffで初期化する．
|CellStart|0xff|0xff|0xff|0xff|
| CellEnd|0xff|0xff|0xff|0xff|
ソート済みのGridParticleHashの各要素iに対して一つ前の要素...
GridParticleHash[i] != GridParticleHash[i-1]
ならば，iを始点としてCellStart[GridParticleHash[i]]に，
i-1を一つ前のグリッドセルの終点としてCellEnd[GridParticle...
|CellStart|0|3|6|0xff|
| CellEnd|2|5|7|0xff|
このときパーティクル位置などを格納した配列をSortedIndexを...
SortedPos[i] = Pos[SortedIndex[i]];
近傍探索フェーズの手順は以下．
+パーティクル位置 (or 任意の座標値) x のグリッドセルを算出
+周囲のグリッドを含めて以下の処理を実行
++グリッドハッシュ値hash計算
++CellStart,CellEndを参照して(CellStart[hash],CellEnd[has...
++for(int j = start; j <= end; ++j) でループ
+++近傍候補パーティクル位置pjを参照
pj=SortedPos[j]
+++|pj-x|<hならそのパーティクルは近傍
SPHで用いる場合は，近傍探索フェースの最後で見つかった近傍...
また，グリッド幅は有効半径hに基づいて決定する．
各パーティクルの密度値を計算する例を以下に示す．
#code(C){{
/*!
* 与えられたセル内のパーティクルとの距離から密度を計算
* @param[in] gridPos グリッド位置
* @param[in] index パーティクルインデックス
* @param[in] pos 計算座標
* @param[in] cell パーティクルグリッドデータ
* @return セル内のパーティクルから計算した密度値
*/
__device__
float calDensityCell(int3 gridPos, uint i, float3 pos0, r...
{
uint gridHash = calcGridHash(gridPos);
// セル内のパーティクルのスタートインデックス
uint startIndex = cell.dCellStart[gridHash];
float h = params.EffectiveRadius;
float dens = 0.0f;
if(startIndex != 0xffffffff){ // セルが空でないかのチェ...
// セル内のパーティクルで反復
uint endIndex = cell.dCellEnd[gridHash];
for(uint j = startIndex; j < endIndex; ++j){
//if(j == i) continue;
float3 pos1 = make_float3(cell.dSortedPos[j]);
float3 rij = pos0-pos1;
float r = length(rij);
if(r <= h){
float q = h*h-r*r;
dens += params.Mass*params.Wpoly6*q*q*q;
}
}
}
return dens;
}
/*!
* パーティクル密度計算(カーネル関数)
* @param[out] newDens パーティクル密度
* @param[out] newPres パーティクル圧力
* @param[in] cell パーティクルグリッドデータ
*/
__global__
void sphCalDensity(float* newDens, float* newPres, rxPart...
{
// パーティクルインデックス
uint index = __mul24(blockIdx.x,blockDim.x)+threadIdx.x;
if(index >= cell.uNumParticles) return;
float3 pos = make_float3(cell.dSortedPos[index]); // パ...
float h = params.EffectiveRadius;
// パーティクル周囲のグリッド
int3 grid_pos0, grid_pos1;
grid_pos0 = calcGridPos(pos-make_float3(h));
grid_pos1 = calcGridPos(pos+make_float3(h));
// 周囲のグリッドも含めて近傍探索，密度計算
float dens = 0.0f;
for(int z = grid_pos0.z; z <= grid_pos1.z; ++z){
for(int y = grid_pos0.y; y <= grid_pos1.y; ++y){
for(int x = grid_pos0.x; x <= grid_pos1.x; ++x){
int3 n_grid_pos = make_int3(x, y, z);
dens += calDensityCell(n_grid_pos, index, pos, cell);
}
}
}
// 密度値を結果に書き込み
uint oIdx = cell.dSortedIndex[index];
newDens[oIdx] = dens;
}
}}
デバイス関数calDensityCellはセル内のパーティクルを参照し...
*表面メッシュ生成 [#ca575a8b]
パーティクルをそのままレンダリングしたのでは流体には見え...
特に透明な液体をレンダリングする際には，液体表面の位置と...
ここでは3DCGにおいて一般的に用いられている三角形メッシュ...
三角形メッシュによる液体表面抽出にはMC(Marching Cubes)法...
MC法では，陰関数で表現されたフィールドから等値面を抽出す...
手順としては，
+立方体グリッド(Cube)を空間に配置
+各Cubeの辺上で陰関数値が0となる点を線型補間や二分法，ニ...
+Cubeに含まれる点の配置と数から三角形メッシュを生成
グリッドに含まれるメッシュ数やメッシュの構成はテーブルを...
より詳しい手順やソースコード，テーブルは以下のサイトを参...
-http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/
GPUでの実装は[[CUDAでMC法(改造版)]]参照．
*実装結果 [#k4a9f2ef]
GeForceGTX580の環境でテストした．
パーティクル数が約25000で，SPH計算だけなら200fps，メッシ...
ちなみにCPUでも同様の方法で実装したところ，SPH計算だけで1...
(メッシュ化には上記のGPU実装MC法を用いたが，GPUへのデータ...
図3,4に実行結果を示す．
CENTER:&ref(sph_result1.jpg,,50%); &ref(sph_result2.jpg,...
CENTER:図3 左からパーティクル，表面メッシュ，屈折レンダリ...
CENTER:&ref(sph_result4.jpg,,50%);
CENTER:図4 ポリゴンで表された谷形状内での流れ
*ソースコード [#cf51a873]
***CUDA9.1用ソースコード [#rebfb11d]
Visual Studio 2015 + CUDA9.1の環境で動作確認したコードは...
#ref(rx_sph_v2.1.zip)
-必要ライブラリや説明は下記のCUDA5用のものとほぼ同じだが...
-動作確認済み環境 : Visual Studio 2015 (Version 14.0.2543...
***CUDA5用ソースコード [#y2c41c83]
Visual Studio 2010 + CUDA5.0の環境で作成したコードは以下．
#ref(rx_sph_v2.zip)
-ビルドするのに必要なライブラリ
freeglut, GLEW, CUDA, boost, libpng, zlib, libjpeg, ftgl...
各ライブラリについては[[ライブラリのインストール]]を参照．
-簡単な説明＆注意事項
--GUIツールキットとしてFLTKを用いている．
--rx_fltk_glcanvas.hの59行目
#define RX_USE_GPU
をコメントアウトするとCPUで計算する．
--[[Anisotropic Kernel]]によるメッシュ生成にも対応してい...
--Visual Studioから実行する場合は，作業ディレクトリを"../...
--CUDA5のSDKにあるヘッダファイル(helper_*)を用いているの...
C:\ProgramData\NVIDIA Corporation\CUDA Samples\v5.0\comm...
にインストールされる．
--"s"キー or Start/Stopボタンでシミュレーションスタート/...
--F1〜F9キー or Sceneメニューからシーン切替可能．binフォ...
これがシーンを記述するファイル．デフォルトでは7つのシーン...
--bin/cubeフォルダには何も入っていないが，屈折レンダリン...
サンプルとしてterragenで作成したものを以下に置く．
#ref(cubemap_texture_example.zip)
--[[OBJ形式]]のファイルを読み込んで，固体オブジェクトにす...
--研究用に作ったコードを削って軽くした物なので，いろいろ...
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
