Skip to content

Commit

Permalink
SPH working fine
Browse files Browse the repository at this point in the history
  • Loading branch information
@AsulconS
AsulconS committed Aug 11, 2020
1 parent 8112f7a commit 754c5ba
Show file tree
Hide file tree
Showing 8 changed files with 430 additions and 89 deletions.
24 changes: 24 additions & 0 deletions include/HSGIL/math/mUtils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,30 @@ HSGIL_API int absolute(int val);
*/
HSGIL_API float absolute(float val);

/**
* @brief Function that returns the power of an integer value with an exponent
*
* @param val
* @return float
*/
HSGIL_API int npow(int val, unsigned int exp);

/**
* @brief Function that returns the power of an unsigned integer value with an exponent
*
* @param val
* @return float
*/
HSGIL_API unsigned int npow(unsigned int val, unsigned int exp);

/**
* @brief Function that returns the power of a floating point value with an exponent
*
* @param val
* @return float
*/
HSGIL_API float npow(float val, unsigned int exp);

/**
* @brief Clamping function to limit an integer value between 2 ranges
*
Expand Down
2 changes: 1 addition & 1 deletion include/HSGIL/window/formWindow.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,7 +173,7 @@ class HSGIL_API FormWindow : public IWindow
*
* @return float
*/
virtual float getAspectRatio() override;
virtual float getAspectRatio() const override;

protected:
/**
Expand Down
2 changes: 1 addition & 1 deletion include/HSGIL/window/iWindow.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ class HSGIL_API IWindow
*
* @return float
*/
virtual float getAspectRatio() = 0;
virtual float getAspectRatio() const = 0;

protected:
/**
Expand Down
2 changes: 1 addition & 1 deletion include/HSGIL/window/renderingWindow.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ class HSGIL_API RenderingWindow : public IWindow
*
* @return float
*/
virtual float getAspectRatio() override;
virtual float getAspectRatio() const override;

protected:
/**
Expand Down
167 changes: 82 additions & 85 deletions main.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,26 +1,25 @@
#include <HSGIL/hsgil.hpp>
#include <particle.hpp>

#include <random>
#include <vector>

#include <particle.hpp>

// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
struct SIM_State
{
GLuint VAO;
GLuint VBO;
float* vertexData;

Particle* particles;
unsigned int nParticles;

float density;
float gasConstant;
GLuint stride;
std::vector<float> vertexData;

float h;
float h2;
std::vector<Particle> particles;

float timeStep;
float restDensity;
float mass;
float viscosity;
float gasStiffness;
float supportRadius;

float damping;
float margin;
Expand All @@ -33,33 +32,33 @@ struct SIM_State
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
// Kernels
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
float W(const float q, const float h)
float poly6DefaultKernel(const gil::Vec3f r, const float h)
{
return q * 315.0f / (65.0f * gil::constants::PI * pow(h, 9.0f));
return (315.0f / (64.0f * gil::constants::PI * pow(h, 9.0f))) * pow(h * h - pow(gil::module(r), 2.0f), 3.0f);
}

float W1(const float q, const float h)
gil::Vec3f spikyGradientKernel(const gil::Vec3f r, const float h)
{
return q * -45.0f / (gil::constants::PI * pow(h, 6.0f));
return (-45.0f / (gil::constants::PI * pow(h, 6.0f))) * gil::normalize(r) * pow(h - gil::module(r), 2.0f);
}

float W2(const float q, const float h)
float viscosityLaplacianKernel(const gil::Vec3f r, const float h)
{
return q * 45.0f / (gil::constants::PI * pow(h, 6.0f));
return (45.0f / (gil::constants::PI * pow(h, 6.0f))) * (h - gil::module(r));
}
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------

// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
// Euler Solver
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
void eulerIntegrate(const SIM_State& sim, const float step)
void eulerIntegrate(SIM_State& sim)
{
for(unsigned int i = 0; i < sim.nParticles; ++i)
for(unsigned int i = 0; i < sim.particles.size(); ++i)
{
Particle& pi = sim.particles[i];

pi.v += step * pi.f / sim.mass;
pi.r += step * pi.v;
pi.v += sim.timeStep * pi.f / pi.density;
pi.r += sim.timeStep * pi.v;

if(pi.r.x - sim.margin < 0.0f)
{
Expand All @@ -76,13 +75,11 @@ void eulerIntegrate(const SIM_State& sim, const float step)
pi.v.y *= sim.damping;
pi.r.y = sim.margin;
}
/*
if(pi.r.y + sim.margin > sim.boundaryHeight)
{
pi.v.y *= sim.damping;
pi.r.y = sim.boundaryHeight - sim.margin;
}
*/
if(pi.r.z - sim.margin < 0.0f)
{
pi.v.z *= sim.damping;
Expand All @@ -99,7 +96,7 @@ void eulerIntegrate(const SIM_State& sim, const float step)

// Init Functions
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
void initGLParams(const SIM_State& sim, gil::RenderingWindow& window, gil::Shader& shader) // Fix this Internally (const Window)
void initGLParams(const SIM_State& sim, const gil::RenderingWindow& window, gil::Shader& shader)
{
shader.use();

Expand All @@ -109,57 +106,62 @@ void initGLParams(const SIM_State& sim, gil::RenderingWindow& window, gil::Shade
glClearColor(0.8f, 0.8f, 0.8f, 1.0f);
shader.setFloat("pointSize", 32.0f);

glm::vec3 viewPos {128.0f, 128.0f, 256.0f};
//glm::vec3 viewPos {1.0f, 1.0f, 2.0f};
glm::vec3 viewPos {0.4f, 0.8f, 1.6f};
glm::mat4 view = glm::lookAt(viewPos, glm::vec3{0.0f, 0.0f, 0.0f}, glm::vec3{0.0f, 1.0f, 0.0f});
glm::mat4 projection = glm::perspective(45.0f, window.getAspectRatio(), 0.1f, 1000.0f);
shader.setMat4("view", view);
shader.setMat4("projection", projection);
shader.setVec3("particleColor", {1.0f, 0.13f, 0.0f});
}

void initSPH(SIM_State& sim)
{
sim.vertexData = new float[6 * sim.nParticles];
sim.particles = new Particle[sim.nParticles];

Particle p;
gil::Vec3f pos;
unsigned int p {0};
for(pos.x = sim.boundaryWidth * 0.0f; pos.x < sim.boundaryWidth * 0.5f; pos.x += sim.h * 0.6f)
for(pos.x = sim.margin; pos.x < sim.boundaryWidth * 0.5f; pos.x += sim.supportRadius * 0.6f)
{
for (pos.y = sim.boundaryHeight * 0.0f; pos.y < sim.boundaryHeight * 0.5f; pos.y += sim.h * 0.6f)
for(pos.y = sim.margin; pos.y < sim.boundaryHeight * 0.5f; pos.y += sim.supportRadius * 0.6f)
{
for (pos.z = sim.boundaryDepth * 0.0f; pos.z < sim.boundaryDepth * 0.5f; pos.z += sim.h * 0.6f)
for(pos.z = sim.margin; pos.z < sim.boundaryDepth * 0.5f; pos.z += sim.supportRadius * 0.6f)
{
sim.particles[p].r = pos;
sim.particles[p].v = {8.0f, 32.0f, 8.0f};
sim.particles[p].color = {1.0f, 0.13f, 0.0f};
sim.vertexData[p * 3 + 3] = 1.0f;
sim.vertexData[p * 3 + 4] = 0.13f;
sim.vertexData[p * 3 + 5] = 0.0f;
++p;
p.r = pos;
p.v = {0.0f, 0.0f, 0.0f};
p.color = {1.0f, 0.13f, 0.0f};
sim.particles.push_back(std::move(p));

// Positions
sim.vertexData.push_back(0.0f);
sim.vertexData.push_back(0.0f);
sim.vertexData.push_back(0.0f);
// Colors
sim.vertexData.push_back(1.0f);
sim.vertexData.push_back(0.13f);
sim.vertexData.push_back(0.0f);
}
}
}
sim.nParticles = p;
sim.stride = 6;

glGenVertexArrays(1, &sim.VAO);
glGenBuffers(1, &sim.VBO);

glBindVertexArray(sim.VAO);

glBindBuffer(GL_ARRAY_BUFFER, sim.VBO);
glBufferData(GL_ARRAY_BUFFER, 6 * sim.nParticles * sizeof(float), sim.vertexData, GL_STATIC_DRAW);
glBufferData(GL_ARRAY_BUFFER, sim.stride * sim.particles.size() * sizeof(float), sim.vertexData.data(), GL_STATIC_DRAW);

// Position Attrib
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 6 * sizeof(float), (void*)0);
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, sim.stride * sizeof(float), (void*)0);
glEnableVertexAttribArray(0);

// Particle Color
glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 6 * sizeof(float), (void*)(3 * sizeof(float)));
glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, sim.stride * sizeof(float), (void*)(3 * sizeof(float)));
glEnableVertexAttribArray(1);

glBindVertexArray(0);

std::cout << "Initialized with " << sim.nParticles << " particles" << std::endl;
std::cout << "Initialized with " << sim.particles.size() << " particles" << std::endl;
}
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------

Expand All @@ -172,22 +174,19 @@ int main()
}

SIM_State sim;
sim.nParticles = 100000;

sim.density = 1000.0f;
sim.gasConstant = 2000.0f;

sim.h = 16.0f;
sim.h2 = sim.h * sim.h;
sim.timeStep = 0.01f;
sim.restDensity = 998.29f;
sim.mass = 0.02f;
sim.viscosity = 3.5f;
sim.gasStiffness = 3.0f;
sim.supportRadius = 0.0457f;

sim.mass = 64.0f;
sim.viscosity = 250.0f;

sim.margin = sim.h;
sim.margin = sim.supportRadius;
sim.damping = -0.5f;
sim.boundaryWidth = 160.0f;
sim.boundaryHeight = 160.0f;
sim.boundaryDepth = 160.0f;
sim.boundaryWidth = 0.6f;
sim.boundaryHeight = 0.6f;
sim.boundaryDepth = 0.6f;

initSPH(sim);

Expand All @@ -202,66 +201,66 @@ int main()
window.pollEvents();

// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
// Compute Density-Pressure
// Compute Mass-Density and Pressure
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
for(unsigned int i = 0; i < sim.nParticles; ++i)
for(unsigned int i = 0; i < sim.particles.size(); ++i)
{
Particle& pi = sim.particles[i];

pi.density = 0;
for(unsigned int j = 0; j < sim.nParticles; ++j)
for(unsigned int j = 0; j < sim.particles.size(); ++j)
{
Particle& pj = sim.particles[j];
float r2 {gil::module(pj.r - pi.r) * gil::module(pj.r - pi.r)};
gil::Vec3f r {pi.r - pj.r};

if(r2 < sim.h2)
if(gil::module(r) < sim.supportRadius)
{
pi.density += sim.mass * W((sim.h2 - r2) * (sim.h2 - r2) * (sim.h2 - r2), sim.h);
pi.density += sim.mass * poly6DefaultKernel(r, sim.supportRadius);
}
// ^ This makes it better (I don't know why profe :'v) ...sim.mass * W(r, sim.h)... antigua version
}
//pi.density += 64.0f;
pi.pressure = sim.gasConstant * (sim.particles[i].density - sim.density);
pi.pressure = sim.gasStiffness * (sim.particles[i].density - sim.restDensity);
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
}

// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
// Compute Forces
// Compute Internal and External Forces
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
for(unsigned int i = 0; i < sim.nParticles; ++i)
for(unsigned int i = 0; i < sim.particles.size(); ++i)
{
Particle& pi = sim.particles[i];

gil::Vec3f fp {0.0f, 0.0f, 0.0f};
gil::Vec3f fv {0.0f, 0.0f, 0.0f};
gil::Vec3f pressureForce {0.0f, 0.0f, 0.0f};
gil::Vec3f viscosityForce {0.0f, 0.0f, 0.0f};
gil::Vec3f gravityForce {0.0f, -gil::constants::GAL, 0.0f};

for(unsigned int j = 0; j < sim.nParticles; ++j)
for(unsigned int j = 0; j < sim.particles.size(); ++j)
{
if(&pi == &sim.particles[j])
{
continue;
}

Particle& pj = sim.particles[j];
float r {gil::module(pj.r - pi.r)};
gil::Vec3f r {pi.r - pj.r};

if(r < sim.h)
if(gil::module(r) < sim.supportRadius)
{
fp += -1.0f * gil::normalize(pj.r - pi.r) * sim.mass * (pi.pressure + pj.pressure) / (2.0f * pj.density) * W1((sim.h - r) * (sim.h - r), sim.h); // <- Lo mismo aqui
fv += sim.viscosity * sim.mass * ((pj.v - pi.v) / pj.density) * W2(sim.h - r, sim.h);
pressureForce += ((pi.pressure / pow(pi.density, 2)) + (pj.pressure / pow(pj.density, 2))) * sim.mass * spikyGradientKernel(r, sim.supportRadius);
viscosityForce += (pj.v - pi.v) * (sim.mass / pj.density) * viscosityLaplacianKernel(r, sim.supportRadius);
}
}

gil::Vec3f g {0.0f, -gil::constants::GAL * 12000, 0.0f};
gil::Vec3f fg {g * pi.density};
pi.f = fp + fv + fg;
pressureForce *= -pi.density;
viscosityForce *= sim.viscosity;
gravityForce *= sim.restDensity;
pi.f = pressureForce + viscosityForce + gravityForce;
}
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------

// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
// Integrate by Euler 1th Order
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
eulerIntegrate(sim, timer.getDeltaTime());
timer.getDeltaTime();
eulerIntegrate(sim);
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------

// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
Expand All @@ -271,19 +270,19 @@ int main()
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

glm::mat4 model;
for(unsigned int i = 0; i < sim.nParticles; ++i)
for(unsigned int i = 0; i < sim.particles.size(); ++i)
{
Particle& pi = sim.particles[i];

model = glm::mat4(1.0f);
model = glm::translate(model, glm::vec3{pi.r.x, pi.r.y, pi.r.z});
shader.setMat4("model", model);

float zColorDepth {1.0f - (pi.r.z / sim.boundaryDepth)};
float zColorDepth {pi.r.z / sim.boundaryDepth};
shader.setVec3("colorDepth", {zColorDepth, zColorDepth, zColorDepth});

glBindVertexArray(sim.VAO);
glDrawArrays(GL_POINTS, 0, sim.nParticles);
glDrawArrays(GL_POINTS, 0, sim.particles.size());
glBindVertexArray(0);
}

Expand All @@ -296,7 +295,5 @@ int main()
glDeleteVertexArrays(1, &sim.VAO);
glDeleteBuffers(1, &sim.VBO);

delete[] sim.vertexData;
delete[] sim.particles;
return 0;
}
Loading

0 comments on commit 754c5ba

Please sign in to comment.