【CG物理模拟系列】流体模拟--粒子法之SPH（实现）

邻域搜索的效率化

SPH等粒子法，由于需要考虑到邻域粒子带来的影响，通常邻域搜索都会消耗大量时间。如果我们只是单纯的计算所有粒子组合的欧氏距离的话，计算时间只会呈指数增加。

而空间分割法的出现，使邻域搜索实现了效率化。 空间分割法是一种，把希望检索到的物体存在的空间以格子等方式分隔开，只计算自身及相邻分割领域内包含物体的距离，可以使计算时间大幅减少的方法。空间分割法的处理顺序是，

1. 物体的记录
   1. 分割空间
   2. 计算出各物体所属的分割领域
2. 邻域搜索
   1. 计算出 搜索中心坐标所属区域
   2. 列出自身和相邻区域的物体
   3. 计算列出物体间的距离 

如上所述，分成记录和检索两部分。记录只在场景中的物体（粒子）位置变化的时候实行，检索则是在有必要的时候连续进行。

空间分割法的难点在于空间以怎样的方式分割。其代表性方法为，

• 等间距格子(下图左)
• 八叉树，四叉树构造(下图右)
• kD树构造

等等。下图的四叉树构造，与等间距网格相比，检索效率更高，但是物体记录却更费事。使用粒子法的时候，为了使每一帧粒子都能实现动态移动，每一帧都要进行区域记录。为此，分割·记录较为容易的等间距网格法则会被经常使用。等间距网格，与有分层结构的其他分割方法相比，更容易进行记录同步化处理。 

　

図1 等间距格子(左)和四叉树(左)

↑

GPU的实现

参照Particle Simulation using CUDA (PDF)(code源文件在NVIDIA的GPU Computing SDK中) ，来实现SPH法的扩充。

首先，介绍使用Particle Simulation using CUDA进行粒子搜索方法，如图所示。

図2 网格和粒子

第一步进行粒子记录处理。顺序如下。

1. 计算各粒子的网格哈希值。哈希值是各网格单元特有的，可以用来处理网格单元的位置信息。 这里，从网格单元的位置(i,j)计算hash值为

   hash = i+j・nx

   (nx,ny)是网格单元的数量，根据图2，各区域的左下角写的数字就是hash值(nx=2)。 哈希值储存在GridParticleHash序列中。同时，在SortedIndex中保存粒子的索引值。(GridParticleHash,SortedIndex的size = 粒子数)GridParticleHash20011012SortedIndex01234567
2. 以网格哈希值作为key排序。这里使用(radix sort) 排序。GridParticleHash00011122SortedIndex12534607
3. 储存GridParticleHash中值相同的连续区域（相同分割区域）的起始与终止索引(CellStart好CellEnd序列)。 首先，以0xffffffff 初始化序列(size=index数)。CellStart0xff0xff0xff0xffCellEnd0xff0xff0xff0xff搜索排序后GridParticleHash各元素i的前一元素(i-1)的hash值。

   GridParticleHash[i] != GridParticleHash[i-1]

   这样的话，以i作为起始点CellStart[GridParticleHash[i]]， i-1作为前一个网格单元的终点CellEnd[GridParticleHash[i-1]]储存。CellStart0360xffCellEnd2570xff这时，使用SortedIndex对粒子位置等的储存序列进行排序。

   SortedPos[i] = Pos[SortedIndex[i]];

邻域搜索的顺序如下。

1. 算出粒子位置(or 任意坐标) x 网格单元
2. 对包含周围网格在内进行如下处理
   1. 计算网格单元的哈希值hash
   2. 按照CellStart,CellEnd(CellStart[hash],CellEnd[hash])，对网格内粒子的起始与终止索引(start, end)进行处理
   3. for(int j = start; j <= end; ++j)
      1. 邻域候选粒子位置pj如下

         pj=SortedPos[j]

      2. 若|pj-x|<h成立，粒子为邻域粒子。

使用SPH法的时候，通过计算与临近搜索阶段最后找到的临近粒子间的距离来计算kernel。再以有效半径h为依据，判断是否继续执行。 

各粒子的密度计算代码如下所示。

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75

 

/*! * 通过单元内粒子见得距离计算密度 * @param[in] gridPos 网格位置 * @param[in] index 粒子索引 * @param[in] pos 计算坐标 * @param[in] cell 粒子网格数据 * @return 密度值 */__device__float calDensityCell(int3 gridPos, uint i, float3 pos0, rxParticleCell cell){    uint gridHash = calcGridHash(gridPos);     // cell内粒子的开始索引    uint startIndex = cell.dCellStart[gridHash];     float h = params.EffectiveRadius;    float dens = 0.0f;    if(startIndex != 0xffffffff){    // cell非空        // cell内粒子循环        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;} /*! * 粒子密度计算（kernel函数） * @param[out] newDens 粒子密度 * @param[out] newPres 粒子压力 * @param[in]  cell 粒子网格数据 */__global__void sphCalDensity(float* newDens, float* newPres, rxParticleCell cell){    // 粒子检索    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;}

↑

生成表面mesh

此时，如果直接渲染粒子的话，会看不到流体。特别是透明流体渲染的时候，必须要计算出液体表面的位置和法线。这里，我们使用三角形mesh生成液体表面。

通过MC(Marching Cubes)法来实现 液体表面生成三角形mesh 。MC法的处理顺序为，

1. 在空间内部署立方grid(Cube)
2. 用插值法或二分法，牛顿法等 算出各Cube边上隐函数值为0的点
3. 从Cube内包含点的构造和数中生成三角形mesh

实现结果

GeForceGTX580环境下。粒子数约25000个，只计算SPH需200fps，加上面的生成约80-90fps。图3,4是运行结果。

　　

图3 左到右：粒子，表面mesh，折射渲染

源码链接

使用FLTK工具包，在Visual Studio 2010 + CUDA5.0环境下，

所需库

freeglut, GLEW, CUDA, boost, libpng, zlib, libjpeg, ftgl, FLTK

GPU安装源码链接：https://github.com/DonDracula/OpenGL_projects/tree/master/rx_sph(GPU)