update

run in debug mode can somewhat avoid a segfault.
Don't know why.

src/Makefile

@@ -1,5 +1,5 @@
 CXX=llvm-g++
-COMPILE=-c -std=c++17 \
+COMPILE=-g -c -std=c++17 \
 -I/usr/local/Cellar/glew/2.1.0/include \
 -I/usr/local/Cellar/glfw/3.3/include \
 -I/usr/local/Cellar/eigen/3.3.7/include/eigen3 \

src/constants.h

@@ -40,6 +40,7 @@ const static glm::vec2 worldOrigin(0.f, 0.f);
 const static float sdfCellSize = H_INV;
 const static int sdfWidth = X_GRID;
 const static int sdfHeight = Y_GRID;
+const static float NARROW_BAND = 1.f;
 
 /* ----- TRANSFER ----- */
 #if INTERPOLATION == 1

src/main.cpp

@@ -101,7 +101,7 @@ int main(int argc, char **argv) {
     // AddParticles();
 
     // P2G, compute grid forces, etc.
-    // Update();
+    Update();
 
     // Display frame at desired rate
     if (t_count % (int)(DT_render / DT) == 0) {
@@ -112,43 +112,42 @@ int main(int argc, char **argv) {
 
       glfwSwapBuffers(window);
 
-      // #if RECORD_VIDEO
-      //       std::string dir = "../result/sim";
-      //
-      //       // zero padding
-      //       std::string num = std::to_string(frameNumber);
-      //       num = std::string(4 - num.length(), '0') + num;
-      //       std::string output = dir + num + ".png";
-      //
-      //       // use opencv to save a frame
-      //       // the frame size becomes twice from frame 2,
-      //       // thus multiplying 2 is needed.
-      //       // don't know why
-      //       cv::Mat outImg(Y_WINDOW * 2, X_WINDOW * 2, CV_8UC3);
-      //       glPixelStorei(GL_PACK_ALIGNMENT, 4);
-      //       glReadPixels(0, 0, X_WINDOW * 2, Y_WINDOW * 2, GL_BGR,
-      //       GL_UNSIGNED_BYTE,
-      //                    (GLvoid *)outImg.data);
-      //       cv::flip(outImg, outImg, 0);
-      //       cv::imwrite(output, outImg);
-      //
-      //       std::cout << "Frame " << frameNumber << " saved." << '\n';
-      //
-      //       frameNumber++;
-      // #endif
+#if RECORD_VIDEO
+      std::string dir = "../result/sim";
+
+      // zero padding
+      std::string num = std::to_string(frameNumber);
+      num = std::string(4 - num.length(), '0') + num;
+      std::string output = dir + num + ".png";
+
+      // use opencv to save a frame
+      // the frame size becomes twice from frame 2,
+      // thus multiplying 2 is needed.
+      // don't know why
+      cv::Mat outImg(Y_WINDOW * 2, X_WINDOW * 2, CV_8UC3);
+      glPixelStorei(GL_PACK_ALIGNMENT, 4);
+      glReadPixels(0, 0, X_WINDOW * 2, Y_WINDOW * 2, GL_BGR, GL_UNSIGNED_BYTE,
+                   (GLvoid *)outImg.data);
+      cv::flip(outImg, outImg, 0);
+      cv::imwrite(output, outImg);
+
+      std::cout << "Frame " << frameNumber << " saved." << '\n';
+
+      frameNumber++;
+#endif
 
       glfwPollEvents();
     } // end outer if
 
     // Don't forget to reset grid every frame !!
-    // Simulation->ResetGrid();
+    Simulation->ResetGrid();
 
     t_count++;
   } // end while
 
-  // #if RECORD_VIDEO
-  //   images2video();
-  // #endif
+#if RECORD_VIDEO
+  images2video();
+#endif
 
   glfwTerminate();
 #endif

src/node.h

@@ -11,7 +11,7 @@ class Node {
   /* Data */
   float Mi; // Node mass
 
-  Vector2f Xi;     // Node position
+  Vector2f Xi;     // Node position (index, also world space position)
   Vector2f Vi;     // Node velocity, before update and after force
   Vector2f Vi_col; // Node velocity, after collision
   Vector2f Vi_fri; // Node velocity, after friction

src/particle.h

@@ -191,7 +191,7 @@ class Snow : public Particle {
     float VOL = 2 * PI * R_BALL * R_BALL / static_cast<float>(NP);
     float MASS = VOL * RHO_snow / 100.0f;
 
-    Vector2f v = Vector2f(0.f, -9.8f); // Initial velocity
+    Vector2f v = Vector2f(0.f, -9.8f * 10.f); // Initial velocity
     Matrix2f a = Matrix2f(0);
 
     // get position from image

src/sdf.h

@@ -59,7 +59,7 @@ class Polygon {
 
     polys.push_back(square);
     polys.push_back(object2);
-    polys.push_back(object3);
+    // polys.push_back(object3);
 
     return polys;
   }

src/solver.cpp

@@ -125,6 +125,7 @@ float Solver::getDistance(glm::vec2 p) {
     return 9999.f;
   } else {
     int idx = idx_x + idx_y * (sdfWidth - 1);
+    // std::cout << idx << '\n';
     return nodes[idx].sdfDist;
   }
 }
@@ -159,6 +160,84 @@ glm::vec2 Solver::getGradient(glm::vec2 p) {
   return gd;
 }
 
+void Solver::applySdfCollision(Node &node) {
+  // node world space position
+  glm::vec2 pos(node.Xi[0], node.Xi[1]);
+  glm::vec2 v(node.Vi[0], node.Vi[1]);
+
+  // std::cout << glm::to_string(pos) << '\n';
+
+  // distance
+  float dist = getDistance(pos);
+  glm::vec2 n;
+  // Polygon *co; // collision object
+  bool isCollisionOn = false;
+
+  glm::vec2 vco(0.f, 0.f);
+
+  // 基于 sdf 的碰撞检测
+  // 和边界处的碰撞检测
+  // 唯一不同的是法向量 n 的计算方式
+  // 而碰撞后的速度，采用同样��处理
+  // for sdf collision detection
+  // narrow band threshold
+  if (glm::abs(dist) < NARROW_BAND) {
+    n = -getGradient(pos);
+    // co = getPolygon(pos); // get collision object
+    isCollisionOn = true;
+
+    // object velocity
+    // linear
+    // vec2 vlin = co->v;
+    glm::vec2 vlin(0.f, 0.f);
+
+    // rotational
+    // vec2 center = (co->lb + co->rt) * 0.5f;
+    // vec2 r = pos - center;
+    //
+    // vec2 vrot;
+    // vrot.x = -r.y / length(r);
+    // vrot.y = r.x / length(r);
+    glm::vec2 vrot(0.f, 0.f);
+
+    // for visualization convenience
+    float scale = 0.1f;
+
+    // vrot *= co->omega * length(r) * scale;
+
+    vco = vlin + vrot;
+
+    // relative velocity
+    glm::vec2 vrel = v - vco;
+
+    float vnLength = glm::dot(n, vrel);
+    // std::cout << vnLength << ", " << isCollisionOn << '\n';
+
+    // if not separating
+    if (vnLength < 0.f && isCollisionOn) {
+      glm::vec2 vt = vrel - vnLength * n;
+      float mu = 0.55f;
+
+      if (length(vt) <= -mu * vnLength) {
+        vrel = glm::vec2(0.f, 0.f);
+      } else {
+        vrel = vt + mu * vnLength * normalize(vt);
+      }
+
+      // back to particle velocity
+      v = vrel + vco;
+
+      // std::cout << v.x << ", " << v.y << '\n';
+
+      node.Vi[0] = v.x;
+      node.Vi[1] = v.y;
+
+      // std::cout << "node velocity: " << node.Vi[0] << ", " << node.Vi[1]
+      //           << '\n';
+    }
+  }
+}
+
 /* -----------------------------------------------------------------------
 |					MATERIAL POINT METHOD ALGORITHM
 |
@@ -233,7 +312,9 @@ void Solver::UpdateNodes() {
       nodes[i].Vi += nodes[i].Fi;
 
       // Apply collisions and frictions
-      nodes[i].NodeCollisions();
+      applySdfCollision(nodes[i]); // sdf-based
+
+      nodes[i].NodeCollisions(); // borders
 #if FRICTION
       nodes[i].NodeFrictions();
 #else

src/solver.h

@@ -43,6 +43,7 @@ class Solver {
   float nearest_distance(glm::vec2, Polygon &);
   float getDistance(glm::vec2);
   glm::vec2 getGradient(glm::vec2);
+  void applySdfCollision(Node &);
 
   void Draw(); // Draw particles, border and nodes (if selected)
   void WriteToFile(