Question: The program provided is a simulation of particles with the SPH method (hydrodynamics of smoothed particles). The code is derived from the MIT licensed implementation by Lucas V. Schuermann, who explains in detail the idea behind an SPH fluid simulation on his blog. The simulation presented on this blog is kept simple for educational purposes. . . This greatly limits the number of particles that can be added to it! Your task will be to optimize this code to be able to run it with as many particles as possible.
Files & Code:
ParticleManager.cpp : (that's the file where you need to code)
#include #include #include #include #include
#include "Constants.h" #include "ParticleManager.h"
using namespace std;
Particle::Particle(double _x, double _y) : x(_x), y(_y), vx(0.f), vy(0.f), fx(0.f), fy(0.f), rho(0.f), p(0.f) { }
ParticleManager::ParticleManager() { ax = 0; ay = GRAVITY;
renderMode = "particles"; }
void ParticleManager::init(unsigned long n) { cout
particles.clear(); particles.reserve(n);
while(particles.size()
double x = rand() / (double) (RAND_MAX) * SCREEN_WIDTH; double y = rand() / (double) (RAND_MAX) * SCREEN_HEIGHT;
double centerDist = sqrt( pow(x - SCREEN_WIDTH / 2.0, 2) + pow(y - SCREEN_HEIGHT / 2.0, 2));
if(centerDist
void ParticleManager::addBlock(double center_x, double center_y) { for(int i=0; i
if(x >= 0 && x = 0 && y
cout
void ParticleManager::addOne(double x, double y) { particles.push_back(Particle(x, y)); cout
void ParticleManager::setGravity(int direction) { switch(direction) { case DOWN: ax = 0; ay = +GRAVITY; break; case UP: ax = 0; ay = -GRAVITY; break; case RIGHT: ax = +GRAVITY; ay = 0; break; default: ax = -GRAVITY; ay = 0; } }
void ParticleManager::explode() { for(auto &p : particles) { p.vx = rand()/(double)RAND_MAX * 10000 - 5000; p.vy = rand()/(double)RAND_MAX * 10000 - 5000; } }
void ParticleManager::integrate(double dt) { for(auto &p : particles) { // forward Euler integration if(p.rho != 0 && p.fx == p.fx && p.fy == p.fy) { p.vx += dt*p.fx/p.rho; p.vy += dt*p.fy/p.rho; }
p.x += dt*p.vx; p.y += dt*p.vy;
// enforce boundary conditions if(p.x - PARTICLE_RADIUS SCREEN_WIDTH) { p.vx *= BOUND_DAMPING; p.x = SCREEN_WIDTH - PARTICLE_RADIUS; } if(p.y - PARTICLE_RADIUS SCREEN_HEIGHT) { p.vy *= BOUND_DAMPING; p.y = SCREEN_HEIGHT - PARTICLE_RADIUS; } } }
void ParticleManager::computeDensityPressure() { // Pour chaque particule for(auto &pi : particles) { pi.rho = 0.f;
// Search for all particules that contribute to the pressure and density // for(auto &pj : particles) { double distance = sqrt( pow(pj.x - pi.x, 2) + pow(pj.y - pi.y, 2));
if(distance
void ParticleManager::computeForces() { // Pour chaque particule for(auto &pi : particles) { double pressure_x = 0.f; double pressure_y = 0.f;
double viscosity_x = 0.f; double viscosity_y = 0.f;
// Calculate the sum of the pressure and viscosity // applied by the other particules for(auto &pj : particles) { if(&pi == &pj) continue;
double r = sqrt( pow(pj.x - pi.x, 2) + pow(pj.y - pi.y, 2));
if(r
// compute viscosity force contribution viscosity_x += VISC * MASS * (pj.vx - pi.vx) / pj.rho * VISC_LAP * (H-r); viscosity_y += VISC * MASS * (pj.vy - pi.vy) / pj.rho * VISC_LAP * (H-r); } }
pi.fx = pressure_x + viscosity_x + ax * pi.rho; pi.fy = pressure_y + viscosity_y + ay * pi.rho; } }
void ParticleManager::update(unsigned long dt) {
// TODO : calculate the grid here
computeDensityPressure(); computeForces(); integrate(dt / 10000.0f); }
void ParticleManager::renderParticles(SDL_Renderer* renderer) { SDL_SetRenderDrawColor(renderer, 230, 120, 0, 100); SDL_Rect r;
// Draw particles for(long unsigned int i=0; i void ParticleManager::renderGrid(SDL_Renderer* renderer) {
// TODO : Draw the lines to make the grids cout
}
void ParticleManager::renderCells(SDL_Renderer* renderer) {
// TODO : Draw the boxes in different colors // Represent the number of particle in each box // // Use le calcul : // int alpha = nb_particles * 255 / 5; // if(alpha > 255) alpha = 255; // SDL_SetRenderDrawColor(renderer, 0, 0, 255, alpha); // // To effect the color in the box cout
}
void ParticleManager::render(SDL_Renderer* renderer) { if(renderMode.find("particle") != string::npos) renderParticles(renderer);
if(renderMode.find("grid") != string::npos) { renderCells(renderer); renderGrid(renderer); } }
Eile Edit View Project Build Debug Test Analyze Tools Extensions Window Help Search (Ctrl+Q) Attach... - a E 17.11 Constants.h 6 X ^ Miscellaneous Files (Global Scope) 1 #pragma once 2 3 #define WINDOW_TITLE "Smoothed-particle hydrodynamics simulation" 4 #define SCREEN_WIDTH 720 5 #define SCREEN_HEIGHT 480 6 7 #define DOWN O 8 #define UP 1 9 #define LEFT 2 10 #define RIGHT 3 11 Game.cpp -X Miscellaneous Files 1 E include 2 Ninclude 3 4 Ninclude 5 6 Ninclude "Game.h" 7 using namespace std; 9 10 EiGame::Game(): 11 n_keepPlaying (true) 12 13 14 ET/ 16 Main game Loop 17 18 See 19 http://gameprogrammingpatterns.com/game-loop.html https://gafferongames.com/post/fix_your timestep/ 21 22 Eivoid Gane::loop(SDL_Renderer renderer) { SDL_Event event; 26 27 unsigned int lastFrameTime = SDL_GetTicks() - AsPerStep: unsigned int timeAccumulator - @ while (@_keepPlaying) { unsigned int frameTime - SDL_GetTicks(); unsigned int dt = frameTime - lastFrameTime; // clamp dt to avoid the spiral of death when the game lags if (dt > 188) dt = 188; 30 31 32 33 34 35 36 37 38 39 40 41. 42 43 timeAccumulator + dt; // -- Events -- while (SDL_PollEvent&event)) { handleInput(event); // -- Update -- while (timeAccumulator >= msPerStep) { update(msPerStep); timeAccumulator - asPerStep: updatetimeAccumulator); timeAccumulator = 0; 45 46 47 48 49 se 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 render(renderer); int wait = frameTime + msPerStep - SDL_GetTicks(); 1NL if (wait e) { cout Ninclude 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 Ninclude *Constants.h" Ninclude "GameSPH.h" using namespace std; Fint main(int argc, char* args[]) { ET if (SDL_Init(SDL_INIT_VIDEO) Ninclude Ninclude 28 21 22 // Particle data structure Eclass Particle { public: Particle( double, double); double x,y; // Position double vx, vy; // Velocity double fx, fy; // Total forces double rho; // Density double p; // Pressure 26 27 E: Eclass ParticleManager { public: const double REST_DENS = 288.0; // rest density const double GAS_CONST = 280.8; // const for equation of state const double H = 16.0; // kernel radius // const double HSQ = pow(H, 2.0); // radius2 for optimization const double MASS = 65.8; // assume all particles have the same mass const double VISC = 25.0; // viscosity constant const double GRAVITY = 28888; 37 // size of a particle const double PARTICLE_RADIUS = H / 4.8; // smoothing kernels defined in Mller and their gradients const double POLY6 = 315.8/(65.8*M_PI*pow(H, 9.@)); const double SPIKY_GRAD = -45.8/(M_PI*pow(H, 6.0)); const double VISC_LAP - 45.8/(M_PI*pow(H, 6.0)); // simulation paraneters const double BOUND DAMPING = -0.9; ParticleManager(); void init(unsigned long); void addone(double, double); void addBlock (double, double); void setGravity(int); void explode(); void update(unsigned long); void render (SDL_Renderer*); std::string render Mode; private: double ax, ay; // Gravity void integrate(double dt); void computeDensityPressure(); void computeForces(); std::vector particles; 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 void render Particles (SDL_Renderer*); void renderGrid(SDL_Renderer*); void renderCells(SDL_Renderer*); // TODO: utiliser une grille pour acclrer la simulation std::vector<:vector>> grid; Eile Edit View Project Build Debug Test Analyze Tools Extensions Window Help Search (Ctrl+Q) Attach... - a E 17.11 Constants.h 6 X ^ Miscellaneous Files (Global Scope) 1 #pragma once 2 3 #define WINDOW_TITLE "Smoothed-particle hydrodynamics simulation" 4 #define SCREEN_WIDTH 720 5 #define SCREEN_HEIGHT 480 6 7 #define DOWN O 8 #define UP 1 9 #define LEFT 2 10 #define RIGHT 3 11 Game.cpp -X Miscellaneous Files 1 E include 2 Ninclude 3 4 Ninclude 5 6 Ninclude "Game.h" 7 using namespace std; 9 10 EiGame::Game(): 11 n_keepPlaying (true) 12 13 14 ET/ 16 Main game Loop 17 18 See 19 http://gameprogrammingpatterns.com/game-loop.html https://gafferongames.com/post/fix_your timestep/ 21 22 Eivoid Gane::loop(SDL_Renderer renderer) { SDL_Event event; 26 27 unsigned int lastFrameTime = SDL_GetTicks() - AsPerStep: unsigned int timeAccumulator - @ while (@_keepPlaying) { unsigned int frameTime - SDL_GetTicks(); unsigned int dt = frameTime - lastFrameTime; // clamp dt to avoid the spiral of death when the game lags if (dt > 188) dt = 188; 30 31 32 33 34 35 36 37 38 39 40 41. 42 43 timeAccumulator + dt; // -- Events -- while (SDL_PollEvent&event)) { handleInput(event); // -- Update -- while (timeAccumulator >= msPerStep) { update(msPerStep); timeAccumulator - asPerStep: updatetimeAccumulator); timeAccumulator = 0; 45 46 47 48 49 se 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 render(renderer); int wait = frameTime + msPerStep - SDL_GetTicks(); 1NL if (wait e) { cout Ninclude 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 Ninclude *Constants.h" Ninclude "GameSPH.h" using namespace std; Fint main(int argc, char* args[]) { ET if (SDL_Init(SDL_INIT_VIDEO) Ninclude Ninclude 28 21 22 // Particle data structure Eclass Particle { public: Particle( double, double); double x,y; // Position double vx, vy; // Velocity double fx, fy; // Total forces double rho; // Density double p; // Pressure 26 27 E: Eclass ParticleManager { public: const double REST_DENS = 288.0; // rest density const double GAS_CONST = 280.8; // const for equation of state const double H = 16.0; // kernel radius // const double HSQ = pow(H, 2.0); // radius2 for optimization const double MASS = 65.8; // assume all particles have the same mass const double VISC = 25.0; // viscosity constant const double GRAVITY = 28888; 37 // size of a particle const double PARTICLE_RADIUS = H / 4.8; // smoothing kernels defined in Mller and their gradients const double POLY6 = 315.8/(65.8*M_PI*pow(H, 9.@)); const double SPIKY_GRAD = -45.8/(M_PI*pow(H, 6.0)); const double VISC_LAP - 45.8/(M_PI*pow(H, 6.0)); // simulation paraneters const double BOUND DAMPING = -0.9; ParticleManager(); void init(unsigned long); void addone(double, double); void addBlock (double, double); void setGravity(int); void explode(); void update(unsigned long); void render (SDL_Renderer*); std::string render Mode; private: double ax, ay; // Gravity void integrate(double dt); void computeDensityPressure(); void computeForces(); std::vector particles; 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 void render Particles (SDL_Renderer*); void renderGrid(SDL_Renderer*); void renderCells(SDL_Renderer*); // TODO: utiliser une grille pour acclrer la simulation std::vector<:vector>> grid
