流体模拟（二）

SPH算法实现1：

由于我们计算每个粒子的状态时，都需要获得在它光滑核半径内（邻域内）的所有粒子信息。我们如果遍历每个粒子，计算欧式距离的话，那开销就过于庞大了。因此我们可以将我们的空间划分成多个网格。

如上图所示，我们可以把空间划分成均匀网格，从而使我们在获取每个粒子的邻域粒子时，只需要访问该粒子的周围网格获取邻域粒子。这时候我们可以使用哈希表的方式，给每个空间网格赋予编号，根据网格号在哈希表中获取粒子索引。当然，一个网格中可能存在多个粒子。我们则可以建立一个冲突链表，在哈希冲突的时候，将粒子索引通过链表的方式连接。也可以

我们可以将每个网格设置为2倍的光滑核半径大小，这样可以只遍历4个网格就可以获取邻域粒子。方法如下图：

我们在访问当前粒子p的时候，我们根据其位置找到距离其最近的网格点m，由于网格长度为2r（光滑核半径）。所以我们只需要在m周围的四个网格便可以搜索到所有的领域粒子，而不是去遍历粒子p点周围的9个网格。即我们只需要遍历s点所在的网格以及x，y轴方向各加1的网格。如果在3维情况，则只需要遍历2*2*2，8个网格。

根据上述所说。我们的流体系统便需要一个粒子缓存类，来保存所有的粒子信息。一个空间网格划分类，包含空间网格的编号以及哈希表信息。由于在每帧迭代所有的粒子信息都会改变，每个粒子的邻域信息都会改变，所以我们还需要建立一个邻接表类，来保存每个粒子的邻域粒子的信息。

粒子缓存类

    struct Point{glm::vec3 pos;//点位置float density;//密度float pressure;//压力glm::vec3 accel;//加速度glm::vec3 velocity;//速度glm::vec3 velocity_eval;//最后速度int next;//指向下一个点的索引};class PointBuffer{public:PointBuffer();void reset(unsigned int capacity);//重置粒子点缓冲unsigned int size()const{return m_fluidCounts;}//返回粒子个数//获取索引为index的点Point* get(unsigned int index){return m_FluidBuf+index;}const Point* get(unsigned int index) const {return m_FluidBuf+index;}Point* addPointReuse();//缓冲中加入新的点并返回virtual ~PointBuffer();private:Point* m_FluidBuf;//粒子点缓存unsigned int m_fluidCounts;//粒子个数unsigned int m_BufCapcity;//粒子容量};

首先创建单个粒子的信息结构体。保存位置，密度，压力，加速度，速度以及最终速度。而next便是我们防止哈希冲突的一个简单链表。若空则赋值-1。这里我使用名为GLM的数学库来处理矩阵以及向量的操作。

接下来是粒子缓存类，定义了三个私有变量：所有的粒子缓存，粒子个数以及粒子容量。设定了重置粒子缓存，以及加入粒子函数和获取相应索引粒子的公有函数。

具体实现如下：

    PointBuffer::PointBuffer():m_FluidBuf(0),m_fluidCounts(0),m_BufCapcity(0){}PointBuffer::~PointBuffer(){delete[] m_FluidBuf;m_FluidBuf=nullptr;}void PointBuffer::reset(unsigned int capacity){m_BufCapcity=capacity;if(m_FluidBuf!=0){//当点缓存不为空，则清空点缓存delete [] m_FluidBuf;m_FluidBuf=nullptr;}if(m_BufCapcity>0)//给点缓存分配空间m_FluidBuf=new Point[m_BufCapcity]();m_fluidCounts=0;//设置点数量为0}Point* PointBuffer::addPointReuse(){if(m_fluidCounts>=m_BufCapcity){if(m_BufCapcity*2>ELEM_MAX){//当超过上限时，返回最后值return m_FluidBuf+m_fluidCounts-1;}m_BufCapcity*=2;Point* new_data=new Point[m_BufCapcity]();memcpy(new_data, m_FluidBuf, m_fluidCounts*sizeof(Point));delete [] m_FluidBuf;m_FluidBuf=new_data;}//新的点Point* point=m_FluidBuf+m_fluidCounts++;point->pos=glm::vec3(0);point->next=0;point->velocity=glm::vec3(0);point->velocity_eval=glm::vec3(0);point->pressure=0;point->density=0;point->accel=glm::vec3(0);return point;}

这里在每次添加粒子的时候，若粒子数量超过分配缓存的容量（m_BufCapcity）而还没超过粒子最大容量（ELEM_MAX）时，则扩充至双倍容量。而且在每次加入新的点的时候，都给其赋初值0。这里我们的粒子根据每次调用addPointReuse（）函数都会拥有自己的编号，编号为加入缓存的顺序。

空间网格类

class fBox3{public:fBox3():min(glm::vec3(0)),max(glm::vec3(0)){}fBox3(glm::vec3 aMin,glm::vec3 aMax):min(aMin),max(aMax){}glm::vec3 min,max;};class GridContainer{public:GridContainer(){}~GridContainer(){}//空间细分void init(const fBox3& box,float sim_scale,float cell_size,float border,int* rex);void insertParticles(PointBuffer*pointBuffer);void findCells(const glm::vec3& p,float radius,int *gridCell);void findTwoCells(const glm::vec3& p,float radius,int *gridCell);int findCell(const glm::vec3& p);int getGridData(int gridIndex);const glm::ivec3* getGridRes()const{return &m_GridRes;}const glm::vec3* getGridMin(void) const { return &m_GridMin; }const glm::vec3* getGridMax(void) const { return &m_GridMax; }const glm::vec3* getGridSize(void) const { return &m_GridSize; }int getGridCellIndex(float px,float py,float pz);//获取对应xyz下的网格索引float getDelta(){return m_GridDelta.x;}private://空间网格std::vector<int> m_gridData;//表格信息（储存网格内的当前的粒子）glm::vec3 m_GridMin;//表格左下角glm::vec3 m_GridMax;//表格右上角glm::ivec3 m_GridRes;//表格规格（n * m * l）glm::vec3 m_GridSize;//表格c大小glm::vec3 m_GridDelta;//表格偏移量float m_GridCellSize;//一个格子大小（通常为2倍的光滑核半径）};

我们首先创建一个BOX3类，让我们的流体在一个立方体空间内模拟。我们空间网格类包含私有成员：m_gridData表示网格内的粒子，若该网格有多个粒子，则该粒子是哈希链表的头节点。m_Gridmin,m_GridMax为网格最左下角和右上角。m_GridRes为网格规格（如10*10*10）。m_GridSize表示每个网格实际尺寸。m_GridDelta为网格偏移量。

函数包括插入粒子，寻找对应pos的单元格（findCell），获取领域单元格（findcells），获取两倍邻域单元格（findTwocells用于各项异性网格重建），获取相应网格索引内的粒子索引（getGridData有多个粒子返回next表中的第一个）。

实现代码如下：

int GridContainer::getGridData(int gridIndex){if(gridIndex<0||gridIndex>=m_gridData.size())return -1;return m_gridData[gridIndex];}int GridContainer::getGridCellIndex(float px, float py, float pz){int gx=(int)((px-m_GridMin.x)*m_GridDelta.x);int gy=(int)((py-m_GridMin.y)*m_GridDelta.y);int gz=(int)((pz-m_GridMin.z)*m_GridDelta.z);return (gz*m_GridRes.y+gy)*m_GridRes.x+gx;}void GridContainer::init(const fBox3 &box, float sim_scale, float cell_size, float border,int* rex){float world_cellsize=cell_size/sim_scale;m_GridMin=box.min;m_GridMin-=border;m_GridMax=box.max;m_GridMax+=border;m_GridSize = m_GridMax;m_GridSize -= m_GridMin;m_GridCellSize=world_cellsize;//网格规格m_GridRes.x=(int)ceil(m_GridSize.x/world_cellsize);m_GridRes.y=(int)ceil(m_GridSize.y/world_cellsize);m_GridRes.z=(int)ceil(m_GridSize.z/world_cellsize);//将网格大小调整为单元格大小的倍数m_GridSize.x=m_GridRes.x*cell_size/sim_scale;m_GridSize.y=m_GridRes.y*cell_size/sim_scale;m_GridSize.z=m_GridRes.z*cell_size/sim_scale;//计算偏移量m_GridDelta=m_GridRes;m_GridDelta/=m_GridSize;int gridTotal=(int)(m_GridRes.x*m_GridRes.y*m_GridRes.z);rex[0]=m_GridRes.x*8;rex[1]=m_GridRes.y*8;rex[2]=m_GridRes.z*8;m_gridData.resize(gridTotal);}void GridContainer::insertParticles(PointBuffer *pointBuffer){std::fill(m_gridData.begin(), m_gridData.end(), -1);Point* p=pointBuffer->get(0);for(unsigned int n=0;n<pointBuffer->size();n++,p++){int gs=getGridCellIndex(p->pos.x, p->pos.y, p->pos.z);//每个网格内的点划分为一个链表(m_gridData[gs]是该网格中链表的头节点)if(gs>=0&&gs<m_gridData.size()){p->next=m_gridData[gs];m_gridData[gs]=n;}else p->next=-1;}}int GridContainer::findCell(const glm::vec3 &p){int gc=getGridCellIndex(p.x, p.y, p.z);if(gc<0||gc>=m_gridData.size())return -1;return gc;}void GridContainer::findCells(const glm::vec3 &p, float radius, int *gridCell){for(int i=0;i<8;i++)gridCell[i]=-1;//计算当前粒子点光滑核所在网格范围int sph_min_x=((-radius+p.x-m_GridMin.x)*m_GridDelta.x);int sph_min_y=((-radius+p.y-m_GridMin.y)*m_GridDelta.y);int sph_min_z=((-radius+p.z-m_GridMin.z)*m_GridDelta.z);if ( sph_min_x < 0 ) sph_min_x = 0;if ( sph_min_y < 0 ) sph_min_y = 0;if ( sph_min_z < 0 ) sph_min_z = 0;//获取8个网格gridCell[0]=(sph_min_z*m_GridRes.y+sph_min_y)*m_GridRes.x+sph_min_x;gridCell[1]=gridCell[0]+1;gridCell[2]=gridCell[0]+m_GridRes.x;gridCell[3]=gridCell[2]+1;if(sph_min_z+1<m_GridRes.z){gridCell[4]=gridCell[0]+m_GridRes.y*m_GridRes.x;gridCell[5]=gridCell[4]+1;gridCell[6]=gridCell[4]+m_GridRes.x;gridCell[7]=gridCell[6]+1;}if(sph_min_x+1>=m_GridRes.x){gridCell[1]=-1;gridCell[3]=-1;gridCell[5]=-1;gridCell[7]=-1;}if(sph_min_y>=m_GridRes.y){gridCell[2]=-1;gridCell[4]=-1;gridCell[6]=-1;gridCell[8]=-1;}}void GridContainer::findTwoCells(const glm::vec3& p,float radius,int *gridCell){for(int i=0;i<64;i++)gridCell[i]=-1;//计算当前粒子点光滑核所在网格范围int sph_min_x=((-radius+p.x-m_GridMin.x)*m_GridDelta.x);int sph_min_y=((-radius+p.y-m_GridMin.y)*m_GridDelta.y);int sph_min_z=((-radius+p.z-m_GridMin.z)*m_GridDelta.z);if ( sph_min_x < 0 ) sph_min_x = 0;if ( sph_min_y < 0 ) sph_min_y = 0;if ( sph_min_z < 0 ) sph_min_z = 0;int base=(sph_min_z*m_GridRes.y+sph_min_y)*m_GridRes.x+sph_min_x;for(int z=0;z<4;z++){for(int y=0;y<4;y++){for(int x=0;x<4;x++){if((sph_min_x+x>=m_GridRes.x)||(sph_min_y+y>=m_GridRes.y||(sph_min_z+z>=m_GridRes.z)))gridCell[16*z+4*y+x]=-1;elsegridCell[16*z+4*y+x]=base+(z*m_GridRes.y+y)*m_GridRes.x+x;}}}}

在初始化中，设定整体空间的大小，网格的规格。由于处于3D空间，我们的findcells函数利用指针gridCell返回当前粒子的周围8个网格的索引。

这里主要提一下插入粒子函数：

void GridContainer::insertParticles(PointBuffer *pointBuffer){std::fill(m_gridData.begin(), m_gridData.end(), -1);Point* p=pointBuffer->get(0);for(unsigned int n=0;n<pointBuffer->size();n++,p++){int gs=getGridCellIndex(p->pos.x, p->pos.y, p->pos.z);//每个网格内的点划分为一个链表(m_gridData[gs]是该网格中链表的头节点)if(gs>=0&&gs<m_gridData.size()){p->next=m_gridData[gs];m_gridData[gs]=n;}else p->next=-1;}}

该函数是将所有的粒子放入网格，生成对应的信息。首先，我们给m_gridData的所有索引赋初值-1。然后根据粒子编号，从第一个粒子pointBuffer->get(0)开始加入网格哈希表。在下面的循环中，我们首先获得当前粒子所在网格的索引。然后根据该索引查找该网格里的粒子。由于在刚开始，网格类都没有粒子，因此第一个加入的粒子，其next则为-1，然后该网格索引的粒子保存为该粒子。若下次在有粒子的位置属于该网格，则这次粒子获取的next则是上次网格保存的粒子索引，然后替换网格保存的索引为新的粒子。因此在多次循环之后，最后进入网格的粒子则为哈希表内冲突链表的头指针了。

邻接表类

#define MAX_NEIGHBOR_COUNTS 80class NeighborTable{public:NeighborTable();void reset(unsigned short pointCounts);//重置邻接表void point_prepare(unsigned short ptIndex);//预备点数据bool point_add_neighbor(unsigned short ptIndex,float distance);//给当前点添加邻接表void point_commit();//递交点给邻接表int getNeighborCounts(unsigned short ptIndex){return m_pointExtraData[ptIndex].neighborCounts;}//获取邻接表中的点个数void getNeighborInfo(unsigned short ptIndex,int index, unsigned short& neighborIndex, float& neighborDistance);//获取索引ptIndex的邻接表中第index个点的数据~NeighborTable();private:union PointExtraData{struct{unsigned neighborDataOffset:24;//偏移unsigned neighborCounts:8;//个数};};PointExtraData* m_pointExtraData;//邻接表信息unsigned int m_pointCounts;//粒子数unsigned int m_pointCapacity;//粒子容量unsigned char* m_neighborDataBuf;// 邻接表的数据缓存unsigned int m_dataBufSize;//bytes 数据缓存尺寸unsigned int m_dataBufOffset;//bytes 数据缓存中的偏移//当前点点数据unsigned short m_currPoint;//索引int m_currNeighborCounts;//邻居点数量unsigned short m_currNeightborIndex[MAX_NEIGHBOR_COUNTS];//邻居中点的索引float m_currNrighborDistance[MAX_NEIGHBOR_COUNTS];//邻居中点的距离void _growDataBuf(unsigned int need_size);//扩容};

这里的m_neighborDataBuf采用了一个数据块去储存邻接表内邻域粒子的索引和距离信息。而m_pointExtraData保存了额外的邻接表信息，包括偏移和个数。由于我们的所有邻接表信息在一个unsigned char的一个数据块中，所以我们在访问每个邻接粒子的时候需要知道它在数据块中的偏移位置，以及它邻接粒子的个数。

然后我们看一下实现的代码：

NeighborTable::NeighborTable():m_pointExtraData(0),m_pointCounts(0),m_pointCapacity(0),m_neighborDataBuf(0),m_dataBufSize(0),m_dataBufOffset(0),m_currPoint(0),m_currNeighborCounts(0){}NeighborTable::~NeighborTable(){if(m_pointExtraData){delete [] m_pointExtraData;m_pointExtraData=nullptr;}if(m_neighborDataBuf){delete [] m_neighborDataBuf;m_neighborDataBuf=nullptr;}}void NeighborTable::reset(unsigned short pointCounts){int a=sizeof(PointExtraData);if(pointCounts>m_pointCapacity){if(m_pointExtraData){delete [] m_pointExtraData;m_pointExtraData=nullptr;}m_pointExtraData=new PointExtraData[a*pointCounts]();m_pointCapacity=m_pointCounts;}m_pointCounts=pointCounts;memset(m_pointExtraData, 0, a*m_pointCapacity);m_dataBufOffset=0;}//准备点void NeighborTable::point_prepare(unsigned short ptIndex){m_currPoint=ptIndex;m_currNeighborCounts=0;}bool NeighborTable::point_add_neighbor(unsigned short ptIndex, float distance){if(m_currNeighborCounts>=MAX_NEIGHBOR_COUNTS)return false;m_currNeightborIndex[m_currNeighborCounts]=ptIndex;m_currNrighborDistance[m_currNeighborCounts]=distance;m_currNeighborCounts++;return true;}void NeighborTable::point_commit(){if(m_currNeighborCounts==0)return;unsigned int index_size=m_currNeighborCounts*sizeof(unsigned short);unsigned int distance_size=m_currNeighborCounts*sizeof(float);//扩大空间if(m_dataBufOffset+index_size+distance_size>m_dataBufSize)_growDataBuf(m_dataBufOffset+index_size+distance_size);//设置邻居数据m_pointExtraData[m_currPoint].neighborCounts=m_currNeighborCounts;m_pointExtraData[m_currPoint].neighborDataOffset=m_dataBufOffset;//复制索引点的信息到数据缓存中memcpy(m_neighborDataBuf+m_dataBufOffset, m_currNeightborIndex, index_size);m_dataBufOffset+=index_size;memcpy(m_neighborDataBuf+m_dataBufOffset, m_currNrighborDistance, distance_size);m_dataBufOffset+=distance_size;}void NeighborTable::_growDataBuf(unsigned int need_size){unsigned int newSize=m_dataBufSize>0?m_dataBufSize:1;while (newSize<need_size) {newSize*=2;}if (newSize<2024) {newSize=1024;}unsigned char* newBuf=new unsigned char[newSize]();if(m_neighborDataBuf){memcpy(newBuf, m_neighborDataBuf, m_dataBufSize);delete [] m_neighborDataBuf;}m_neighborDataBuf=newBuf;m_dataBufSize=newSize;}//ptIndex是点的索引，index是该点邻居表内的邻居索引,获取的neighborIndex为该邻居的默认索引void NeighborTable::getNeighborInfo(unsigned short ptIndex, int index, unsigned short &neighborIndex, float &neighborDistance){PointExtraData neighData=m_pointExtraData[ptIndex];unsigned short* indexBuf=(unsigned short*)(m_neighborDataBuf+neighData.neighborDataOffset);float* distanceBuf=(float*)(m_neighborDataBuf+neighData.neighborDataOffset+sizeof(unsigned short)*neighData.neighborCounts);neighborIndex=indexBuf[index];neighborDistance=distanceBuf[index];}

point_prepare函数用于每次在处理当前粒子之前，先给其赋初值。point_add_neighbor函数则是将每个邻域粒子与当前粒子的欧式距离在光滑核半径内时，我们保存索引和距离信息。point_commit函数将我们保存的邻接粒子信息储存至邻接表内。

到这里我们的三个流体的结构类已经完成了，之后只要根据SPH的数学概念便可以完成流体系统。

SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现1相关推荐

SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现2
流体模拟(二) SPH算法实现2: 在前面一节我们完成了粒子缓存类,网格类和邻接表类.我们现在可以正式的整合在我们的流体系统类中了. 流体系统类 class FluidSystem{public:Fl ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（3）-光滑核函数
流体模拟(二) 光滑核函数: sph中涉及的光滑核可以理解为:在一定的光滑核半径内,所受的力受距离权重的影响,距离越近所受影响越大.其表现形式如图所示. 这里我们便可以将流体看成一个个粒子的集合,每一 ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（2）-粒子受力分析
流体模拟(二) SPH算法的粒子受力分析: SPH算法的基本设想,就是将连续的流体想象成一个个相互作用的微粒,这些粒子相互影响,共同形成了复杂的流体运动.其实现的原理则是我们在初始空间里创建多个粒子, ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（1）-数学原理
流体模拟(二) SPH算法的数学原理: 标量场和矢量场如果空间区域内任意一点P,都有一个确定的数f(P),则称这个空间区域内确定了一个标量场,如果空间区域内任意一点P,都有一个确定的向量vF(P), ...
SPH（光滑粒子流体动力学）流体模拟实现六：Position Based Fluid（PBF）
SPH(光滑粒子流体动力学)流体模拟实现六:Position Based Fluid(PBF) PBF方法和前篇提到的PCISPH方法类似,都属于迭代矫正法.PCISPH是通过迭代预测压力,通过压力变 ...
SPH（光滑粒子流体动力学）流体模拟实现七：屏幕空间流体渲染（SSF）
SPH(光滑粒子流体动力学)流体模拟实现七:屏幕空间流体渲染(SSF) 之前都是用的Marching Cube重建流体表面的方法.近来为了做对比实验,尝试了屏幕空间流体渲染的方法.发现屏幕空间的方法不 ...
SPH（光滑粒子流体动力学）流体模拟实现：算法总览
流体模拟(一) 流体模拟算法总体流程: 流体现象广泛存在于自然界.日常生活以及工业生产中,对流体的模拟即流体动画, 一直是基于物理的动画以及计算机图形学的重要研究内容.目前, 基于物理模拟的流体动画 ...
SPH（光滑粒子流体动力学）流体模拟实现三：Marching Cube算法（1）
流体模拟(三) Marching Cube算法(1) 我们在实现流体表面重建时,需要事先在空间中划分网格,我们的流体系统正好已经完成了此项工作.其次利用Marching Cube算法计算出构成表面的 ...
SPH（光滑粒子流体动力学）流体模拟实现三：Marching Cube算法（2）
流体模拟(三) Marching Cube算法(2) 我们在之前的流体系统类里新加入一些函数和成员,用来引用我们的MC类,便可以获得生成有表面的流体模型了,效果如图: 添加简单的天空盒以及Phong氏 ...

SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现1