记录流体模拟的学习历程：（一）基于opengl的弹簧质点模拟系统

记录流体模拟的学习历程：（一）基于opengl的质点弹簧系统

最近开始学习的流体仿真，发现参考资料真的比较少，所以在此记录本人学习流体仿真的历程，希望以后可以坚持更新到我能完美的模拟出各种流体效果，也希望能给来看的人提供一些帮助。

最近在读《fluid engine development》这本书，真的很好用，读完了第一章，想基于作者给出的算法做出一个质点弹簧系统，但是由于个人原因，不想配置作者书上源码的各种配置，所以没有用作者的矩阵库，而是用的OpenGL常用的glm库。
同时推荐一个B站的OpenGL教学视频，这也是我学习OpenGL所看的视频。最好的OpenGL教程之一。
虽然是英文版，如果听不懂只需要照着他的代码跟着写就可以了，基本不会报错的。

下面先放一个我模拟的效果视频。（顺便问下怎么做gif好用）
这个系统我只分析了重力和弹力，并没有添加阻力等，因为将系统实现出来，别的力再加进去就变得十分简单。
重力F = ma;
弹力根据胡克定律 F = -k*（L-L0）；
在这里L是当前的弹簧长度，L0是弹簧未发生形变时的长度。K为胡克系数。
下面看具体的代码实现。

private:
 //假定的每个质点的质量
 float mass = 1.0f;
 //弹簧不被拉伸时的长度
 float restLength = 1.0f;
 //胡克定律系数
 float stiffness = 10.0f;
 //小球的数量
 size_t numberOfPoints;
public:
 //各个小球的位置信息
 std::vector<glm::vec3> positions;
 //各个小球的速度信息
 std::vector<glm::vec3> velocities;
 //各个小球的受力信息（所受到的全部力）
 std::vector<glm::vec3> forces;
 //各个边的信息
 std::vector<Edge> edges;

 MassSpringSystem(size_t numberOfPoints = 10);
 void calculateForce();
 //更新状态，每刻的加速度，速度 ，位置等。
 void updateStates(float timeIntervalInSeconds);

具体就看注释吧。。。

MassSpringSystem::MassSpringSystem(size_t numberOfPoints)
 :numberOfPoints(10)
{
 size_t numberOfEdges = numberOfPoints - 1;
 //初始化各个数组的长度
 positions.resize(numberOfPoints);
 velocities.resize(numberOfPoints);
 forces.resize(numberOfPoints);
 edges.resize(numberOfEdges);
 for (size_t i = 0; i < numberOfPoints; i++) {
  //static_cast:将i的数据类型转化为float
  //这个就是弄了一下各个点的初始位置
  positions[i].x = static_cast<float>(i);
  positions[i].y = static_cast<float>(i);
 }
 //初始化每条边
 for (size_t i = 0; i < numberOfEdges; i++)
 {
  edges[i] = Edge{i,i+1};
 }
};
//计算受力情况
void MassSpringSystem::calculateForce() {
 size_t numberOfPoints = positions.size();
 size_t numberOfEdges = edges.size();
 glm::vec3 gravity(0.0f, -9.8f, 0.0f);
 for (size_t i = 0; i < numberOfPoints; i++)
 {
  //重力
  forces[i] = mass * gravity;
 }
 for (size_t i = 0; i < numberOfEdges; i++) {
  size_t pointIndex0 = edges[i].first_point;
  size_t pointIndex1 = edges[i].end_point;
  //记录两个质点的位置信息
  glm::vec3 pos0(positions[pointIndex0]);
  glm::vec3 pos1(positions[pointIndex1]);
  glm::vec3 r(pos0 - pos1);
  float distance = glm::distance(pos0,pos1);
  if (distance > 0) {
    glm::vec3 elasticForce
   (-stiffness * (distance - restLength) * glm::normalize(r));
   //两个相连接的点受力方向相反。
   forces[pointIndex0] += elasticForce;
   forces[pointIndex1] -= elasticForce;
  }
 }
}
void MassSpringSystem::updateStates(float timeIntervalInSeconds) {
 for (size_t i = 0; i < numberOfPoints; i++)
 {
 //将第一个物质点固定不动
  if (i == 0)
   continue;
   //计算加速度
  glm::vec3 newAcceleration(forces[i] / mass);
  //计算新的速度 V = V0 + at
  glm::vec3 newVelocity(velocities[i] + 
   (timeIntervalInSeconds * newAcceleration));
  //计算新的位置X = X0 + vt
  glm::vec3 newPosition(positions[i] +
   (timeIntervalInSeconds * newVelocity));
  velocities[i] = newVelocity;
  positions[i] = newPosition;
 }

注释写的比较粗糙，但是我相信你们都可以看懂。。

最后就是绘制了，首先计算出每帧的时间间隔

GLfloat currentFrame = glfwGetTime();
        deltaTime = currentFrame - lastFrame;
        lastFrame = currentFrame;

然后进行绘制

void drawPoint(std::vector<glm::vec3> positions) 
{     /* Draw a point */
     glClearColor(255, 250, 250, 1); 
     glClear(GL_COLOR_BUFFER_BIT);
     glPointSize(5.0f);    
     glBegin(GL_POINTS);    
     glColor3f(0, 0, 0.0);     
     glVertex2f(0, -1);     
     for (int i = 0; i < 10; i++) {        
        glVertex2f(positions[i].x, positions[i].y);     
     }    
     glEnd();     
     glBegin(GL_LINES);
     glLineWidth(2);//设置线段宽度
     glColor3f(1, 0, 0);
     for (int i = 0; i < 9; i++) {
        glVertex2f(positions[i].x, positions[i].y);
        glVertex2f(positions[i + 1].x, positions[i + 1].y);
     }
     glEnd();
}

因为我也是刚开始学，刚开始接触这方面，所以代码肯定会有很多地方不合理，需要优化，同时也希望自己不断进步。

接下来我会在这个基础上做一个布料模拟，希望一个月时间能够完成吧。