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

流体模拟（二）

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);
        p