EN中文
开发者社区 > 程序

OpenGL进阶(十五)：弹簧质点系统（Mass Spring Systems）

发表于2017-09-26
评论1 9k浏览

简介


      一个模拟变形物体最简单额方法就是将其表示为弹簧质点系统（Mass Spring Systems），因为方法简单，这个框架对于学习物理模拟和时间积分的各种技术都比较理想。一个弹簧质点包含了一系列由多个弹簧连接起来的质点，这样的系统的物理属性非常直接，模拟程序也很容易编写。

      但是简单是有代价了：

1.物体的行为依赖于弹簧系统的设置方法；

2.很难通过调整弹簧系数来得到想要的结果；

3.弹簧质点网格不能直接获取体效果。

      在很多的应用中这些缺点可以忽略，在这种场合下，弹簧质点网格是最好的选择，因为够快够简单。

      今天首先说明一下Mass Spring Systems的理论基础，然后来实现绳子的real-time模拟。


物理模型

      假设有N个质点，质量为mi，位置为xi，速度为vi , 1<i<N.

      这些质点由一组弹簧S连接，弹簧参数为（i,  j, lo, ks, kd）. i,j 为连接的弹簧质点，lo为弹簧完全放松时的长度，ks为弹簧弹性系数，kd为阻尼系数，由胡科定律知弹簧施加在两顶点上的力可以表示为：


         由受力守恒知 fi fj = 0. fi和fj的大小和弹簧伸长成正比关系。

        对于阻尼的计算，


          大小和速度成正比，并且符合力守恒，则对于一个质点，其受力方程为：



模拟

在模拟算法中，牛顿第二定律是关键，f = ma。在已知质量和外力的情况下，通过 a=f/m 可以得到加速度，将二阶常微分方程写成两个一阶方程：


可以得到解析解：



初始状态为 v(to) = vo, x(to)=xo.

积分将时间t内所有的变化加和，模拟的过程就是从 to 开始不断地计算 x(t) 和 v(t)，然后更新质点的位置。

最简单的数值解法就是利用有限的差分近似微分：


其中delta t 是为时间步长，t为帧号，带入之前的公式得：



这个就是显式的欧拉解法，下一时刻的状态完全由当前状态决定。

整个过程的伪代码如下：

[cpp] view plain copy
  1. // initialization  
  2.  forall particles i  
  3.          initialize xi , vi and mi  
  4.  endfor  
  5. // simulation loop  
  6.  loop  
  7.         forall particles i  
  8.                fi ← fg   fcoll   ∑ f(xi , vi , x j , v j )  
  9.         endfor  
  10.         forall particles i  
  11.              vi ← vi   ∆t fi /mi  
  12.              xi ← xi   ∆t vi  
  13.          endfor  
  14.          display the system every nth time  
  15.  endloop  

fg表示重力，fcoll表示碰撞外力。

        显式积分是最简单的积分方法，不过有一个严重的问题 ——它只能在步长很小的时候稳定，这就意味着在两个可视帧中间需要进行多个模拟步骤，不稳定的主要原因是欧拉方法仅仅依赖于当前状态，它假设在整个步长中，外力是个常量。

        假设一个场景：两个质点和一根弹簧，轻轻拉伸后松手，两个小球向中间靠拢，这时候如果时间步长过大，由于力 在一个步长中大小不变，最终结果可能导致弹簧越弹越开，最后被拉断。

        更多关于显式与隐式积分，参考显式方法与隐式方法。


模拟一根绳子

       在布料的模拟中，大部分都使用的弹簧质点系统。

       整个模拟的的过程就和伪代码中的描述的一样。

环境：Ubuntu12.04 codeblocks10.04 glm9.5

       代码清单：

Mass类：质点类，表示质点，主要属性是质量，速度，位置

mass.h

  1. #ifndef MASS_H  
  2. #define MASS_H  
  3. #include <glm/glm.hpp>  
  4. #include <iostream>  
  5. class Mass  
  6. {  
  7.     public:  
  8.         Mass(float m);  
  9.         ~Mass();  
  10.         void applyForce(glm::vec3 force);  
  11.         void init();  
  12.         void simulate(float dt);  
  13.         void setPos(glm::vec3 p);  
  14.         void setVel(glm::vec3 v);  
  15.   
  16.         glm::vec3 getPos();  
  17.         glm::vec3 getVel();  
  18.         float getM();  
  19.     protected:  
  20.     private:  
  21.         // The mass value  
  22.         float m;  
  23.         // Position in space  
  24.         glm::vec3 pos;  
  25.         // Velocity  
  26.         glm::vec3 vel;  
  27.         // Force applied on this mass at an instance  
  28.         glm::vec3 force;  
  29.   
  30. };  
  31.   
  32. #endif // MASS_H  

mass.cpp

  1. #include "mass.h"  
  2.   
  3. Mass::Mass(float m)  
  4. {  
  5.     //ctor  
  6.     this->m = m;  
  7. }  
  8.   
  9. Mass::~Mass()  
  10. {  
  11.     //dtor  
  12. }  
  13.   
  14. void Mass::applyForce(glm::vec3 f)  
  15. {  
  16.     this->force  = f;  
  17. }  
  18.   
  19. void Mass::init()  
  20. {  
  21.     force.x = 0;  
  22.     force.y = 0;  
  23.     force.z = 0;  
  24. }  
  25.   
  26. void Mass::simulate(float dt)  
  27. {  
  28.     vel  = (force/m) * dt;  
  29.     pos  = vel * dt;  
  30. }  
  31.   
  32. glm::vec3 Mass::getPos()  
  33. {  
  34.     return pos;  
  35. }  
  36.   
  37. void Mass::setPos(glm::vec3 p)  
  38. {  
  39.     pos = p;  
  40. }  
  41.   
  42. void Mass::setVel(glm::vec3 v)  
  43. {  
  44.     this->vel = v;  
  45. }  
  46.   
  47. glm::vec3 Mass::getVel()  
  48. {  
  49.    return this->vel;  
  50. }  
  51.   
  52.   
  53. float Mass::getM()  
  54. {  
  55.     return this->m;  
  56. }  


Spring类:弹簧类，一个弹簧，连着两个质点，还有弹簧的一些属性。

spring.h

  1. #ifndef SPRING_H  
  2. #define SPRING_H  
  3. #include "mass.h"  
  4. #include <iostream>  
  5. class Spring  
  6. {  
  7.     public:  
  8.         Spring();  
  9.         Spring(Mass* mass1, Mass* mass2,  
  10.                float springConstant, float springLength,  
  11.                 float frictionConstant);  
  12.         ~Spring();  
  13.   
  14.         void solve();  
  15.   
  16.     protected:  
  17.     private:  
  18.         Mass* mass1;  
  19.         Mass* mass2;  
  20.   
  21.         float springConstant;  
  22.         float restLength;  
  23.         float frictionConstant;  
  24.   
  25. };  
  26.   
  27. #endif // SPRING_H  

spring.cpp

[cpp] view plain copy
  1. #include "spring.h"  
  2.   
  3. Spring::Spring()  
  4. {  
  5.     //ctor  
  6. }  
  7.   
  8. Spring::~Spring()  
  9. {  
  10.     //dtor  
  11. }  
  12.   
  13. Spring::Spring(Mass* mass1, Mass* mass2,  
  14.                float springConstant, float springLength,  
  15.                 float frictionConstant)  
  16. {  
  17.     this->springConstant = springConstant;                                   //set the springConstant  
  18.     this->restLength = springLength;                                     //set the springLength  
  19.     this->frictionConstant = frictionConstant;                               //set the frictionConstant  
  20.   
  21.     this->mass1 = mass1;                                                 //set mass1  
  22.     this->mass2 = mass2;  
  23. }  
  24.   
  25. void Spring::solve()  
  26. {  
  27.     glm::vec3 springVector = mass1->getPos() - mass2->getPos();  
  28.     float r = glm::length(springVector);  
  29.   
  30.     glm::vec3 force(0);  
  31.     if(0 != r)  
  32.     {  
  33.         force  = (springVector / r) * (r - restLength) * (-springConstant);  
  34.     }  
  35.     force  = -(mass1->getVel() - mass2->getVel())*frictionConstant;  
  36.   
  37.     mass1->applyForce(force);                                                    //force is applied to mass1  
  38.     mass2->applyForce(-force);  
  39. }  

最重要的Simulator

ropesimulator.h

[cpp] view plain copy
  1. #include "ropesimulator.h"  
  2.   
  3. RopeSimulator::RopeSimulator()  
  4. {  
  5.     //ctor  
  6. }  
  7.   
  8. RopeSimulator::~RopeSimulator()  
  9. {  
  10.     //dtor  
  11. }  
  12.   
  13. RopeSimulator::RopeSimulator(  
  14.         int numOfMasses,                                //1. the number of masses  
  15.         float m,                                        //2. weight of each mass  
  16.         float springConstant,                           //3. how stiff the springs are  
  17.         float springLength,                             //4. the length that a spring does not exert any force  
  18.         float springFrictionConstant,                   //5. inner friction constant of spring  
  19.         glm::vec3 g,                            //6. gravitational acceleration  
  20.         float airFrictionConstant,                      //7. air friction constant  
  21.         float groundRepulsionConstant,                  //8. ground repulsion constant  
  22.         float groundFrictionConstant,                   //9. ground friction constant  
  23.         float groundAbsorptionConstant,                 //10. ground absorption constant  
  24.         float groundHeight                              //11. height of the ground (y position)  
  25.         )  
  26. {  
  27.     this->numOfMasses = numOfMasses;  
  28.     this->gravitation = g;  
  29.     this->airFrictionConstant = airFrictionConstant;  
  30.     this->groundFrictionConstant = groundFrictionConstant;  
  31.     this->groundRepulsionConstant = groundRepulsionConstant;  
  32.     this->groundAbsorptionConstant = groundAbsorptionConstant;  
  33.     this->groundHeight = groundHeight;  
  34.   
  35.     this->masses = new Mass*[numOfMasses];  
  36.   
  37.     for (int count = 0; count < numOfMasses;  count)        // We will step to every pointer in the array  
  38.     {  
  39.         masses[count] = new Mass(m);  
  40.         masses[count]->init();  
  41.      }  
  42.   
  43.     for (int index = 0; index < numOfMasses;  index)            //To set the initial positions of masses loop with for(;;)  
  44.     {  
  45.         masses[index]->setPos(glm::vec3(index * springLength,0,0));  
  46.   
  47.     }  
  48.   
  49.     springs = new Spring*[numOfMasses - 1];  
  50.     for (int index = 0; index < numOfMasses - 1;  index)            //to create each spring, start a loop  
  51.     {  
  52.         //Create the spring with index "a" by the mass with index "a" and another mass with index "a   1".  
  53.         springs[index] = new Spring(masses[index], masses[index   1],  
  54.         springConstant, springLength, springFrictionConstant);  
  55.     }  
  56. }  
  57.   
  58. void RopeSimulator::release()  
  59. {  
  60.     for (int count = 0; count < numOfMasses;  count)        // we will delete all of them  
  61.         {  
  62.             delete(masses[count]);  
  63.             masses[count] = NULL;  
  64.         }  
  65.   
  66.         delete(masses);  
  67.         masses = NULL;  
  68.   
  69.     for (int index = 0; index < numOfMasses - 1;  index)        //to delete all springs, start a loop  
  70.         {  
  71.             delete(springs[index]);  
  72.             springs[index] = NULL;  
  73.         }  
  74.   
  75.         delete(springs);  
  76.         springs = NULL;  
  77. }  
  78.   
  79. void RopeSimulator::simulate(float dt)  
  80. {  
  81.     for (int count = 0; count < numOfMasses;  count)        // We will iterate every mass  
  82.     {  
  83.          masses[count]->simulate(dt);  
  84.            // std::cout<<"mass "<<count<<" pos: "<<masses[count]->getVel().x<<masses[count]->getVel().y<<masses[count]->getVel().z<<std::endl;  
  85.     }  
  86.   
  87.     ropeConnectionPos  = ropeConnectionVel * dt;    //iterate the positon of ropeConnectionPos  
  88.   
  89.     if (ropeConnectionPos.y < groundHeight)          //ropeConnectionPos shall not go under the ground  
  90.     {  
  91.         ropeConnectionPos.y = groundHeight;  
  92.         ropeConnectionVel.y = 0;  
  93.     }  
  94.     masses[0]->setPos(ropeConnectionPos);                //mass with index "0" shall position at ropeConnectionPos  
  95.     masses[0]->setVel(ropeConnectionVel);                //the mass's velocity is set to be equal to ropeConnectionVel  
  96.   
  97. }  
  98.   
  99. void RopeSimulator::setRopeConnectionPos(glm::vec3 p)  
  100. {  
  101.     this->ropeConnectionPos = p;  
  102. }  
  103.   
  104. void RopeSimulator::setRopeConnectionVel(glm::vec3 v)  
  105. {  
  106.     this->ropeConnectionVel = v;  
  107. }  
  108.   
  109. float RopeSimulator::getGroundHeight()  
  110. {  
  111.     return this->groundHeight;  
  112. }  
  113.   
  114. int RopeSimulator::getNumOfMasses()  
  115. {  
  116.     return this->numOfMasses;  
  117. }  
  118.   
  119. Mass* RopeSimulator::getMass(int index)  
  120. {  
  121.    if (index < 0 || index >= numOfMasses)     // if the index is not in the array  
  122.             return NULL;                            // then return NULL  
  123.         return masses[index];  
  124. }  
  125.   
  126. void RopeSimulator::operate(float dt)  
  127. {  
  128.     this->resetMassesForce();  
  129.     this->solve();  
  130.     this->simulate(dt);  
  131. }  
  132.   
  133. void RopeSimulator::solve()  
  134. {  
  135.     for (int index = 0; index < numOfMasses - 1;  index)        //apply force of all springs  
  136.     {  
  137.         springs[index]->solve();                     //Spring with index "a" should apply its force  
  138.   
  139.     }  
  140.   
  141.     for (int index = 0; index < numOfMasses;  index)                //Start a loop to apply forces which are common for all masses  
  142.     {  
  143.         masses[index]->applyForce(gravitation * masses[index]->getM());               //The gravitational force  
  144.         masses[index]->applyForce(-masses[index]->getVel() * airFrictionConstant);    //The air friction  
  145.   
  146.         if (masses[index]->getPos().y < groundHeight)     //Forces from the ground are applied if a mass collides with the ground  
  147.         {  
  148.             glm::vec3 v;                                //A temporary Vector3D  
  149.   
  150.                 v = masses[index]->getVel();                     //get the velocity  
  151.                 v.y = 0;                                //omit the velocity component in y direction  
  152.   
  153.                 //The velocity in y direction is omited because we will apply a friction force to create  
  154.                 //a sliding effect. Sliding is parallel to the ground. Velocity in y direction will be used  
  155.                 //in the absorption effect.  
  156.                 masses[index]->applyForce(-v * groundFrictionConstant);      //ground friction force is applied  
  157.   
  158.                 v = masses[index]->getVel();                     //get the velocity  
  159.                 v.x = 0;                                //omit the x and z components of the velocity  
  160.                 v.z = 0;                                //we will use v in the absorption effect  
  161.   
  162.                 //above, we obtained a velocity which is vertical to the ground and it will be used in  
  163.                 //the absorption force  
  164.   
  165.                 if (v.y < 0)                         //let's absorb energy only when a mass collides towards the ground  
  166.                     masses[index]->applyForce(-v * groundAbsorptionConstant);        //the absorption force is applied  
  167.   
  168.                 //The ground shall repel a mass like a spring.  
  169.                 //By "Vector3D(0, groundRepulsionConstant, 0)" we create a vector in the plane normal direction  
  170.                 //with a magnitude of groundRepulsionConstant.  
  171.                 //By (groundHeight - masses[a]->pos.y) we repel a mass as much as it crashes into the ground.  
  172.                 glm::vec3 force = glm::vec3(0, groundRepulsionConstant, 0) *  
  173.                     (groundHeight - masses[index]->getPos().y);  
  174.   
  175.                 masses[index]->applyForce(force);            //The ground repulsion force is applied  
  176.             }  
  177.   
  178.         }  
  179. }  
  180.   
  181. void RopeSimulator::resetMassesForce()                              // this method will call the init() method of every mass  
  182. {  
  183.     for (int count = 0; count < numOfMasses;  count)        // We will init() every mass  
  184.         masses[count]->init();                       // call init() method of the mass  
  185. }  

ropesimulator.cpp

[cpp] view plain copy
  1. #include "ropesimulator.h"  
  2.   
  3. RopeSimulator::RopeSimulator()  
  4. {  
  5.     //ctor  
  6. }  
  7.   
  8. RopeSimulator::~RopeSimulator()  
  9. {  
  10.     //dtor  
  11. }  
  12.   
  13. RopeSimulator::RopeSimulator(  
  14.         int numOfMasses,                                //1. the number of masses  
  15.         float m,                                        //2. weight of each mass  
  16.         float springConstant,                           //3. how stiff the springs are  
  17.         float springLength,                             //4. the length that a spring does not exert any force  
  18.         float springFrictionConstant,                   //5. inner friction constant of spring  
  19.         glm::vec3 g,                            //6. gravitational acceleration  
  20.         float airFrictionConstant,                      //7. air friction constant  
  21.         float groundRepulsionConstant,                  //8. ground repulsion constant  
  22.         float groundFrictionConstant,                   //9. ground friction constant  
  23.         float groundAbsorptionConstant,                 //10. ground absorption constant  
  24.         float groundHeight                              //11. height of the ground (y position)  
  25.         )  
  26. {  
  27.     this->numOfMasses = numOfMasses;  
  28.     this->gravitation = g;  
  29.     this->airFrictionConstant = airFrictionConstant;  
  30.     this->groundFrictionConstant = groundFrictionConstant;  
  31.     this->groundRepulsionConstant = groundRepulsionConstant;  
  32.     this->groundAbsorptionConstant = groundAbsorptionConstant;  
  33.     this->groundHeight = groundHeight;  
  34.   
  35.     this->masses = new Mass*[numOfMasses];  
  36.   
  37.     for (int count = 0; count < numOfMasses;  count)        // We will step to every pointer in the array  
  38.     {  
  39.         masses[count] = new Mass(m);  
  40.         masses[count]->init();  
  41.      }  
  42.   
  43.     for (int index = 0; index < numOfMasses;  index)            //To set the initial positions of masses loop with for(;;)  
  44.     {  
  45.         masses[index]->setPos(glm::vec3(index * springLength,0,0));  
  46.   
  47.     }  
  48.   
  49.     springs = new Spring*[numOfMasses - 1];  
  50.     for (int index = 0; index < numOfMasses - 1;  index)            //to create each spring, start a loop  
  51.     {  
  52.         //Create the spring with index "a" by the mass with index "a" and another mass with index "a   1".  
  53.         springs[index] = new Spring(masses[index], masses[index   1],  
  54.         springConstant, springLength, springFrictionConstant);  
  55.     }  
  56. }  
  57.   
  58. void RopeSimulator::release()  
  59. {  
  60.     for (int count = 0; count < numOfMasses;  count)        // we will delete all of them  
  61.         {  
  62.             delete(masses[count]);  
  63.             masses[count] = NULL;  
  64.         }  
  65.   
  66.         delete(masses);  
  67.         masses = NULL;  
  68.   
  69.     for (int index = 0; index < numOfMasses - 1;  index)        //to delete all springs, start a loop  
  70.         {  
  71.             delete(springs[index]);  
  72.             springs[index] = NULL;  
  73.         }  
  74.   
  75.         delete(springs);  
  76.         springs = NULL;  
  77. }  
  78.   
  79. void RopeSimulator::simulate(float dt)  
  80. {  
  81.     for (int count = 0; count < numOfMasses;  count)        // We will iterate every mass  
  82.     {  
  83.          masses[count]->simulate(dt);  
  84.            // std::cout<<"mass "<<count<<" pos: "<<masses[count]->getVel().x<<masses[count]->getVel().y<<masses[count]->getVel().z<<std::endl;  
  85.     }  
  86.   
  87.     ropeConnectionPos  = ropeConnectionVel * dt;    //iterate the positon of ropeConnectionPos  
  88.   
  89.     if (ropeConnectionPos.y < groundHeight)          //ropeConnectionPos shall not go under the ground  
  90.     {  
  91.         ropeConnectionPos.y = groundHeight;  
  92.         ropeConnectionVel.y = 0;  
  93.     }  
  94.     masses[0]->setPos(ropeConnectionPos);                //mass with index "0" shall position at ropeConnectionPos  
  95.     masses[0]->setVel(ropeConnectionVel);                //the mass's velocity is set to be equal to ropeConnectionVel  
  96.   
  97. }  
  98.   
  99. void RopeSimulator::setRopeConnectionPos(glm::vec3 p)  
  100. {  
  101.     this->ropeConnectionPos = p;  
  102. }  
  103.   
  104. void RopeSimulator::setRopeConnectionVel(glm::vec3 v)  
  105. {  
  106.     this->ropeConnectionVel = v;  
  107. }  
  108.   
  109. float RopeSimulator::getGroundHeight()  
  110. {  
  111.     return this->groundHeight;  
  112. }  
  113.   
  114. int RopeSimulator::getNumOfMasses()  
  115. {  
  116.     return this->numOfMasses;  
  117. }  
  118.   
  119. Mass* RopeSimulator::getMass(int index)  
  120. {  
  121.    if (index < 0 || index >= numOfMasses)     // if the index is not in the array  
  122.             return NULL;                            // then return NULL  
  123.         return masses[index];  
  124. }  
  125.   
  126. void RopeSimulator::operate(float dt)  
  127. {  
  128.     this->resetMassesForce();  
  129.     this->solve();  
  130.     this->simulate(dt);  
  131. }  
  132.   
  133. void RopeSimulator::solve()  
  134. {  
  135.     for (int index = 0; index < numOfMasses - 1;  index)        //apply force of all springs  
  136.     {  
  137.         springs[index]->solve();                     //Spring with index "a" should apply its force  
  138.   
  139.     }  
  140.   
  141.     for (int index = 0; index < numOfMasses;  index)                //Start a loop to apply forces which are common for all masses  
  142.     {  
  143.         masses[index]->applyForce(gravitation * masses[index]->getM());               //The gravitational force  
  144.         masses[index]->applyForce(-masses[index]->getVel() * airFrictionConstant);    //The air friction  
  145.   
  146.         if (masses[index]->getPos().y < groundHeight)     //Forces from the ground are applied if a mass collides with the ground  
  147.         {  
  148.             glm::vec3 v;                                //A temporary Vector3D  
  149.   
  150.                 v = masses[index]->getVel();                     //get the velocity  
  151.                 v.y = 0;                                //omit the velocity component in y direction  
  152.   
  153.                 //The velocity in y direction is omited because we will apply a friction force to create  
  154.                 //a sliding effect. Sliding is parallel to the ground. Velocity in y direction will be used  
  155.                 //in the absorption effect.  
  156.                 masses[index]->applyForce(-v * groundFrictionConstant);      //ground friction force is applied  
  157.   
  158.                 v = masses[index]->getVel();                     //get the velocity  
  159.                 v.x = 0;                                //omit the x and z components of the velocity  
  160.                 v.z = 0;                                //we will use v in the absorption effect  
  161.   
  162.                 //above, we obtained a velocity which is vertical to the ground and it will be used in  
  163.                 //the absorption force  
  164.   
  165.                 if (v.y < 0)                         //let's absorb energy only when a mass collides towards the ground  
  166.                     masses[index]->applyForce(-v * groundAbsorptionConstant);        //the absorption force is applied  
  167.   
  168.                 //The ground shall repel a mass like a spring.  
  169.                 //By "Vector3D(0, groundRepulsionConstant, 0)" we create a vector in the plane normal direction  
  170.                 //with a magnitude of groundRepulsionConstant.  
  171.                 //By (groundHeight - masses[a]->pos.y) we repel a mass as much as it crashes into the ground.  
  172.                 glm::vec3 force = glm::vec3(0, groundRepulsionConstant, 0) *  
  173.                     (groundHeight - masses[index]->getPos().y);  
  174.   
  175.                 masses[index]->applyForce(force);            //The ground repulsion force is applied  
  176.             }  
  177.   
  178.         }  
  179. }  
  180.   
  181. void RopeSimulator::resetMassesForce()                              // this method will call the init() method of every mass  
  182. {  
  183.     for (int count = 0; count < numOfMasses;  count)        // We will init() every mass  
  184.         masses[count]->init();                       // call init() method of the mass  
  185. }  

主要就是这三个类了，注释比较清楚。最终效果：


如社区发表内容存在侵权行为，您可以点击这里查看侵权投诉指引

标签：

进阶OpenGLListVIAIRArray