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

原创 2016年05月31日 16:25:27
邻域搜索的效率化

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

而空间分割法的出现,使邻域搜索实现了效率化。 空间分割法是一种,把希望检索到的物体存在的空间以格子等方式分隔开,只计算自身及相邻分割领域内包含物体的距离,可以使计算时间大幅减少的方法。空间分割法的处理顺序是,
  1. 物体的记录
    1. 分割空间
    2. 计算出各物体所属的分割领域
  2. 邻域搜索
    1. 计算出 搜索中心坐标所属区域
    2. 列出自身和相邻区域的物体
    3. 计算列出物体间的距离 

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

空间分割法的难点在于空间以怎样的方式分割。其代表性方法为,
  • 等间距格子(下图左)
  • 八叉树,四叉树构造(下图右)
  • kD树构造
等等。下图的四叉树构造,与等间距网格相比,检索效率更高,但是物体记录却更费事。使用粒子法的时候,为了使每一帧粒子都能实现动态移动,每一帧都要进行区域记录。为此,分割·记录较为容易的等间距网格法则会被经常使用。等间距网格,与有分层结构的其他分割方法相比,更容易进行记录同步化处理。 
particle_grid.gif particle_octree.gif
図1 等间距格子(左)和四叉树(左)

GPU的实现


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

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

particle_grid_2.gif
図2 网格和粒子

第一步进行粒子记录处理。顺序如下。
  1. 计算各粒子的网格哈希值。哈希值是各网格单元特有的,可以用来处理网格单元的位置信息。 这里,从网格单元的位置(i,j)计算hash值为
    hash = i+j・nx
    
    (nx,ny)是网格单元的数量,根据图2,各区域的左下角写的数字就是hash值(nx=2)。 哈希值储存在GridParticleHash序列中。同时,在SortedIndex中保存粒子的索引值。(GridParticleHash,SortedIndex的size = 粒子数)
    GridParticleHash 2 0 0 1 1 0 1 2
    SortedIndex 0 1 2 3 4 5 6 7
  2. 以网格哈希值作为key排序。这里使用(radix sort) 排序。
    GridParticleHash 0 0 0 1 1 1 2 2
    SortedIndex 1 2 5 3 4 6 0 7
  3. 储存GridParticleHash中值相同的连续区域(相同分割区域)的起始与终止索引(CellStart好CellEnd序列)。 首先,以0xffffffff 初始化序列(size=index数)。
    CellStart 0xff 0xff 0xff 0xff
    CellEnd 0xff 0xff 0xff 0xff
    搜索排序后GridParticleHash各元素i的前一元素(i-1)的hash值。
    GridParticleHash[i] != GridParticleHash[i-1]
    
    这样的话,以i作为起始点CellStart[GridParticleHash[i]], i-1作为前一个网格单元的终点CellEnd[GridParticleHash[i-1]]储存。
    CellStart 0 3 6 0xff
    CellEnd 2 5 7 0xff
    这时,使用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是运行结果。
sph_result1.jpg sph_result2.jpg sph_result3.jpg
图3 左到右:粒子,表面mesh,折射渲染
sph_result4.jpg

源码链接

使用FLTK工具包,在Visual Studio 2010 + CUDA5.0环境下,
所需库
freeglut, GLEW, CUDA, boost, libpng, zlib, libjpeg, ftgl, FLTK

版权声明:本文为博主原创文章,未经博主允许不得转载。

相关文章推荐

SPH算法简介(四):Hello,SPH

SPH算法简介(四):Hello,SPH 2011年04月2日 |本网站遵守CC版权协议 转载请注明出自www.thecodeway.com     上几节,我们推导出一大推复杂无比的公式,似...

Delphi7高级应用开发随书源码

  • 2003年04月30日 00:00
  • 676KB
  • 下载

教你如何学习自动化测试(QTP)转 自席飞剑

软件测试行业经过这么多年的发展,如今测试行业对从业者的要求是越来越高,不再仅仅局限于要求会写测试用例、会细致的执行测试、能有效的发现系统缺陷等;越来越多的企业对应聘者本身的技能要求也越来越高,招聘信息...

计算机视觉、机器学习相关领域论文和源代码大集合--持续更新……

注:下面有project网站的大部分都有paper和相应的code。Code一般是C/C++或者Matlab代码。 最近一次更新 注:下面有project网站的大部分都有paper和相应的code...

【CG物理模拟系列】流体模拟--粒子法之SPH(代码讲解)

WCSPH,PCISPH,IISPH等研究方法,其本质都是以非压缩性为目标,求解Navier-Stokes方程。 本文以WCSPH为例,讲解下SPH方法代码的实现。 代码讲解 sph_type...

【CG物理模拟系列】流体模拟--粒子法之SPH法的加权函数计算

SPH法的加权函数 Poly6 kernel[1] 用于密度计算等。由于r只存在于2次项中,可以省去平方根的计算。 梯度为, 拉普拉斯算子为, 如图所示, ...

【CG物理模拟系列】流体模拟--粒子法之SPH(理论)

SPH法介绍 SPH(Smoothed Particle Hydrodynamics)是早年银河系碰撞,天体形成等宇宙物理学模拟所使用的方法[Lucy1977],近年来被应用到流体,热等其他现...

【CG物理模拟系列】流体模拟--粒子法之MPS法(理论)

MPS法  前面的文章里我们讲过SPH曾是为了处理压缩性流体问题而提出的方法,与之相对,这一篇来说说用粒子法处理非压缩性流体的研究方法--Moving Particle Semi-implicit...

基于SPH的流体模拟实践和一些技巧总结

SPH的流体模拟是目前大多数游戏所采用的模拟流体方法,特点是简单,十分容易实现,相比与基于Grid的Eulerian方法更加简单和高速,本文主要介绍一下使用SPH的流体模拟中一些常用的技巧和数据结构。...

【CG物理模拟系列】弹性体模拟--Position-based法之Shape Matching

Position-based法 Position-based法与传统的力学基础方法不同,根据构成物体的顶点等元素间的约束条件(Constraint),直接变更顶点的位置坐标的方法。        ...
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:【CG物理模拟系列】流体模拟--粒子法之SPH(实现)
举报原因:
原因补充:

(最多只允许输入30个字)