Fluid Particle Grid

Question

So recently I have been looking at some fluid simulation techniques, specifically Smoothed Particle Hydrodynamics. One of the key components is a grid to avoid O(n^2) searches. I have been having some trouble implement a grid of my own.

I have been able to kind of divide the particles into a grid, and then query a particle ID to get a list of neighbors. Great! However upon further analysis I was receiving list much bigger than I expected, and sometimes exceeding the total number of particles in the simulation, sometimes there would be 140 in a single cell

I have been able to determine that the problem maybe in the UpdateGrid() method, as it seems to be allocating 24,000 particles instead 400. I am really quite lost. Any help would be appreciated.

#include "FluidScene.h"

CFluidScene::CFluidScene(void)
{
    GRAVITY = Vector2D(0, -2);
    EPSILON = pow(10, -6)/ 2;

    CELL_SPACE = 0.15;
}


CFluidScene::~CFluidScene(void)
{
}

void CFluidScene::Start()
{
    std::cout << "fluid scene loading...\n";

    Width = 9;
    Height = 12;

    //make the grid... here we go
    int gridx = Width / CELL_SPACE;
    int gridy = Height / CELL_SPACE;
    Grid = new std::vector<int>[gridx, gridy];

    //spawn some particles
    for (int y = 1; y < 21; y++)
    {
        for (int x = 1; x < 21; x++)
        {
            Positions.push_back(Vector2D((float)x * CELL_SPACE, ((float)y * CELL_SPACE + Height / 1.5)));
            Velocities.push_back(Vector2D::Zero());
            Forces.push_back(Vector2D::Zero());
        }
    }

    CFluidRenderer::Start(&Positions);
    CScene::Start();
}

void CFluidScene::Update(float dTime)
{
    Solve(dTime);
// test getting the neighbours
    std::vector<int> r = GetNeighbours(19);
    if (r.size() > 0)
        std::cout << Positions[GetNeighbours(19)[0]].y << "\n";

    std::cout << r.size() << std::endl;;
    //system("cls");
    CScene::Update(dTime);
}

void CFluidScene::Solve(float dTime)
{
    UpdateGrid();
    CalculateForces();
    UpdateParticles(dTime);
}

void CFluidScene::CalculateForces()
{
    for (int i = 0; i < Positions.size(); i++)
    {
        Forces[i] += GRAVITY;
    }
}

void CFluidScene::UpdateParticles(float dTime)
{
    for (int i = 0; i < Positions.size(); i++)
    {
        Velocities[i] += Forces[i];
        Forces[i] = Vector2D::Zero();
        Positions[i] += Velocities[i] * (dTime * dTime) * 0.99;
        BoundaryCollisions(Positions[i]);
    }
}

void CFluidScene::BoundaryCollisions(Vector2D &pos)
{
    if (pos.x < 0)
    {
        pos.x = 0 + EPSILON;
    }
    if (pos.x > Width)
    {
        pos.x = Width - EPSILON;
    }
    if (pos.y < 0)
    {
        pos.y = 0 + EPSILON;
    }
    if (pos.y > Height)
    {
        pos.y = Height - EPSILON;
    }
}

void CFluidScene::UpdateGrid()
{
    int w = Width / CELL_SPACE;
    int h = Height / CELL_SPACE;
    Grid = new std::vector<int>[w, h];
    int i;
    for (i = 0; i < Positions.size(); i++)
    {
        int gridx = (int)Positions[i].x / CELL_SPACE;
        int gridy = (int)Positions[i].y / CELL_SPACE;
        Grid[gridx, gridy].push_back(i);
    }
}


std::vector<int> CFluidScene::GetNeighbours(int i)
{
    std::vector<int> ret;

    int gridx = (int)Positions[i].x / CELL_SPACE;
    int gridy = (int)Positions[i].y / CELL_SPACE;

    int ax, ay;

    for (int y = -1; y < 2; y++)
    {
        for (int x = -1; x < 2; x++)
        {
            ax = gridx + x;
            if (ax < 0)
                ax = 0;
            ay = gridy + y;
            if (ay < 0)
                ay = 0;
            int j;
            for (j = 0; j < Grid[ax, ay].size(); j++)
            {
                ret.push_back(Grid[ax, ay][j]);
            }
        }
    }

    return ret;
}

Answer 1

You allocate Grid like this: Grid = new std::vector<int>[gridx, gridy];

C++ is not C#, this does not quite what you expect.

gridx, gridy is evaluated (thanks to the comma operator) to just gridy, so you are allocating an array of gridy std::vector<int>s.

The same thing happens when you try to access Grid[x,y], you are actually accessing Grid[y].

Also, you are newing a bunch of vectors on each update, without deleting the old ones anywhere, so you are leaking memory, too. It's not really necessary to allocate the grid again and again, or at all.

You could just use something like std::vector<std::vector<std::vector<int>>> for your grid, and access grid cells like this: Grid[x][y]

You can clear the stored ints like this:

for (int x = 0; x < gridx; ++x)
{
    for (int y = 0; y < gridy; ++y)
    {
        Grid[x][y].clear();
    }
}