IT WOOOOOORKS <3 <3

src/sph/sph_system.cpp

@@ -11,8 +11,16 @@ namespace sph {
                 for (float x = 0; x < 2; x += 0.2f)
                     if (particles.size() < N)
                     {
+
                         float x_off = glm::linearRand(.0f, 1.f);
-                        particles.emplace_back(new Particle(x + x_off, y, z));
+//                        particles.emplace_back(new Particle(x + x_off, y, z));
+                        particles.emplace_back(
+                                new Particle(
+                                        glm::linearRand(1.f, 1.1f),
+                                        glm::linearRand(1.f, 1.1f),
+                                        glm::linearRand(1.f, 1.1f)
+                                        )
+                                );
                     }
 
         std::cout << "There are " << particles.size() << "particles." <<
@@ -32,7 +40,8 @@ namespace sph {
             for(auto &p_j : particles) {
                 glm::vec3 r_ij = p_j->pos - p_i->pos;
                 // euclidean distance between p_i and p_j
-                float r2 = glm::l2Norm(p_j->pos, p_i->pos);
+                float r = glm::l2Norm(p_j->pos, p_i->pos);
+                float r2 = glm::pow(r, 2.f);
                 
                 // 0 <= r <= d
                 if (r2 <= H2)
@@ -53,19 +62,23 @@ namespace sph {
             for(auto &p_j : particles) {
                 if(&p_i == &p_j) continue;
 
-                float r_ij = glm::length(p_j->pos - p_i->pos);
+                glm::vec3 r_ijv = p_i->pos - p_j->pos;
+                float r_ij = glm::length(r_ijv);
                 
                 if(r_ij <= H) {
                     // compute pressure force contribution
-                    f_pressure += -M * r_ij / (2*p_j->rho) *
-                        POLY6_GRAD_PRESS * glm::pow(H - r_ij, 2.f);
+
+                    f_pressure += -M * (p_i->p + p_j->p) / (2*p_j->rho) *
+                        POLY6_GRAD_PRESS * r_ijv/r_ij * glm::pow(H-r_ij, 2.0f);
+
                     f_viscosity += MU * M * (p_j->v - p_i->v) / p_j->rho *
                         POLY6_GRAD_VISC * (H - r_ij);
                 }
 
                 glm::vec3 f_gravity = G * p_i->rho;
 
-                p_i->f = f_pressure + f_viscosity + f_gravity;
+//                p_i->f = f_pressure + f_viscosity + f_gravity;
+                p_i->f = f_gravity + f_pressure;// * 0.001f;
             }
         }
     }
@@ -76,18 +89,22 @@ namespace sph {
             p->v += DT * p->f / p->rho;
             p->pos += DT * p->v;
 
+            const float max_b = 1.2f;
+            const float min_b = 0.9f;
+
             //boundary conditions
-            if (p->pos.x - EPS <= 0) {
+            if (p->pos.x - EPS <= min_b) {
                 p->v.x *= BOUND_DAMPING;
-                p->pos.x = EPS;
+                p->pos.x = min_b + EPS;
             }
 
-            if(p->pos.x + EPS > 2) {
+            if(p->pos.x + EPS > max_b) {
                 p->v.x *= BOUND_DAMPING;
-                p->pos.x = 2 - EPS;
+                p->pos.x = max_b - EPS;
             }
 
-            if(p->pos.y - EPS <= 0) {
+            // ground
+            if(p->pos.y - EPS <= 0.0f) {
                 p->v.y *= BOUND_DAMPING;
                 p->pos.y = EPS;
             }
@@ -98,15 +115,17 @@ namespace sph {
             }
 
             //TODO DEPTH CONSTANT VALUE
-            if(p->pos.z - EPS <= 0) {
+            if(p->pos.z - EPS <= min_b) {
                 p->v.z *= BOUND_DAMPING;
-                p->pos.z = EPS;
+                p->pos.z = min_b + EPS;
             }
 
-            if(p->pos.z + EPS >= 2) {
+            if(p->pos.z + EPS >= max_b) {
                 p->v.z *= BOUND_DAMPING;
-                p->pos.z = 2 - EPS;
+                p->pos.z = max_b - EPS;
             }
+
+            p->v *= 0.98;
         }
     }
 
@@ -114,17 +133,17 @@ namespace sph {
         compute_density_and_pressure();
         compute_forces();
 
-        if(t > DT)
-            for(float dt = 0.0f; dt < t; dt+=DT)
-                integrate();
-        else
+//        if(t > DT)
+//            for(float dt = 0.0f; dt < t; dt+=DT)
+//                integrate();
+//        else
             integrate();
     }
 
     void System::render() {
         // TODO MODERN OPENGL!!!!
         glEnable(GL_POINT_SMOOTH);
-        glPointSize(H / 2.f);
+        glPointSize( 1.f);
 
         glColor4f(0.5f, 0.5f, 0.8f, 1.f);
 

src/sph/sph_system.h

@@ -18,18 +18,18 @@ namespace sph {
     // rest density
     const static float REST_DENS = 1000.0f;
     // const for equation of state
-    const static float GAS_CONST = 2000.0f;
+    const static float GAS_CONST = 3.0f;
     // kernel radius
-    const static float H = 4.f;
+    const static float H = 0.0457f;
     // (kernel radius)^2 for optimization
     const static float H2 = H * H;
     // mass of particles
-    const static float M = 60.f;
+    const static float M = 0.02f;
     // viscosity constant
     const static float MU = 42.f;
     // TODO decide the place for this
     // integration timestep
-    const static float DT = 0.0005f; 
+    const static float DT = 0.01f;
 
     // epsilon (e.g. at boundary check)
     const static float EPS = 0.0001f;