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

简介

      一个模拟变形物体最简单额方法就是将其表示为弹簧质点系统（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为帧号，带入之前的公式得：

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

整个过程的伪代码如下：

// initialization forall particles i         initialize xi , vi and mi endfor// simulation loop loop        forall particles i               fi ← fg + fcoll + ∑ f(xi , vi , x j , v j )        endfor        forall particles i             vi ← vi + ∆t fi /mi             xi ← xi + ∆t vi         endfor         display the system every nth time endloop

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

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

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

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

模拟一根绳子

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

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

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

       代码清单：

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

mass.h

#ifndef MASS_H#define MASS_H#include <glm/glm.hpp>#include <iostream>class Mass{    public:        Mass(float m);        ~Mass();        void applyForce(glm::vec3 force);        void init();        void simulate(float dt);        void setPos(glm::vec3 p);        void setVel(glm::vec3 v);        glm::vec3 getPos();        glm::vec3 getVel();        float getM();    protected:    private:        // The mass value        float m;        // Position in space        glm::vec3 pos;        // Velocity        glm::vec3 vel;        // Force applied on this mass at an instance        glm::vec3 force;};#endif // MASS_H

mass.cpp

#include "mass.h"Mass::Mass(float m){    //ctor    this->m = m;}Mass::~Mass(){    //dtor}void Mass::applyForce(glm::vec3 f){    this->force += f;}void Mass::init(){    force.x = 0;    force.y = 0;    force.z = 0;}void Mass::simulate(float dt){    vel += (force/m) * dt;    pos += vel * dt;}glm::vec3 Mass::getPos(){    return pos;}void Mass::setPos(glm::vec3 p){    pos = p;}void Mass::setVel(glm::vec3 v){    this->vel = v;}glm::vec3 Mass::getVel(){   return this->vel;}float Mass::getM(){    return this->m;}

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

spring.h

#ifndef SPRING_H#define SPRING_H#include "mass.h"#include <iostream>class Spring{    public:        Spring();        Spring(Mass* mass1, Mass* mass2,               float springConstant, float springLength,                float frictionConstant);        ~Spring();        void solve();    protected:    private:        Mass* mass1;        Mass* mass2;        float springConstant;        float restLength;        float frictionConstant;};#endif // SPRING_H

spring.cpp

#include "spring.h"Spring::Spring(){    //ctor}Spring::~Spring(){    //dtor}Spring::Spring(Mass* mass1, Mass* mass2,               float springConstant, float springLength,                float frictionConstant){    this->springConstant = springConstant;//set the springConstant    this->restLength = springLength;//set the springLength    this->frictionConstant = frictionConstant;//set the frictionConstant    this->mass1 = mass1;//set mass1    this->mass2 = mass2;}void Spring::solve(){    glm::vec3 springVector = mass1->getPos() - mass2->getPos();    float r = glm::length(springVector);    glm::vec3 force(0);    if(0 != r)    {        force += (springVector / r) * (r - restLength) * (-springConstant);    }    force += -(mass1->getVel() - mass2->getVel())*frictionConstant;    mass1->applyForce(force);//force is applied to mass1    mass2->applyForce(-force);}

最重要的Simulator

ropesimulator.h

#include "ropesimulator.h"RopeSimulator::RopeSimulator(){    //ctor}RopeSimulator::~RopeSimulator(){    //dtor}RopeSimulator::RopeSimulator(        int numOfMasses,//1. the number of massesfloat m,//2. weight of each massfloat springConstant,//3. how stiff the springs arefloat springLength,//4. the length that a spring does not exert any forcefloat springFrictionConstant,//5. inner friction constant of springglm::vec3 g,//6. gravitational accelerationfloat airFrictionConstant,//7. air friction constantfloat groundRepulsionConstant,//8. ground repulsion constantfloat groundFrictionConstant,//9. ground friction constantfloat groundAbsorptionConstant,//10. ground absorption constantfloat groundHeight//11. height of the ground (y position)){    this->numOfMasses = numOfMasses;    this->gravitation = g;    this->airFrictionConstant = airFrictionConstant;    this->groundFrictionConstant = groundFrictionConstant;    this->groundRepulsionConstant = groundRepulsionConstant;    this->groundAbsorptionConstant = groundAbsorptionConstant;    this->groundHeight = groundHeight;    this->masses = new Mass*[numOfMasses];    for (int count = 0; count < numOfMasses; ++count)// We will step to every pointer in the array    {        masses[count] = new Mass(m);        masses[count]->init();     }    for (int index = 0; index < numOfMasses; ++index)//To set the initial positions of masses loop with for(;;)    {        masses[index]->setPos(glm::vec3(index * springLength,0,0));    }    springs = new Spring*[numOfMasses - 1];    for (int index = 0; index < numOfMasses - 1; ++index)//to create each spring, start a loop    {        //Create the spring with index "a" by the mass with index "a" and another mass with index "a + 1".        springs[index] = new Spring(masses[index], masses[index + 1],        springConstant, springLength, springFrictionConstant);    }}void RopeSimulator::release(){    for (int count = 0; count < numOfMasses; ++count)// we will delete all of them{delete(masses[count]);masses[count] = NULL;}delete(masses);masses = NULL;    for (int index = 0; index < numOfMasses - 1; ++index)//to delete all springs, start a loop{delete(springs[index]);springs[index] = NULL;}delete(springs);springs = NULL;}void RopeSimulator::simulate(float dt){    for (int count = 0; count < numOfMasses; ++count)// We will iterate every mass    { masses[count]->simulate(dt);           // std::cout<<"mass "<<count<<" pos: "<<masses[count]->getVel().x<<masses[count]->getVel().y<<masses[count]->getVel().z<<std::endl;    }    ropeConnectionPos += ropeConnectionVel * dt;//iterate the positon of ropeConnectionPos    if (ropeConnectionPos.y < groundHeight)//ropeConnectionPos shall not go under the ground    {        ropeConnectionPos.y = groundHeight;        ropeConnectionVel.y = 0;    }    masses[0]->setPos(ropeConnectionPos);//mass with index "0" shall position at ropeConnectionPos    masses[0]->setVel(ropeConnectionVel);//the mass's velocity is set to be equal to ropeConnectionVel}void RopeSimulator::setRopeConnectionPos(glm::vec3 p){    this->ropeConnectionPos = p;}void RopeSimulator::setRopeConnectionVel(glm::vec3 v){    this->ropeConnectionVel = v;}float RopeSimulator::getGroundHeight(){    return this->groundHeight;}int RopeSimulator::getNumOfMasses(){    return this->numOfMasses;}Mass* RopeSimulator::getMass(int index){   if (index < 0 || index >= numOfMasses)// if the index is not in the arrayreturn NULL;// then return NULLreturn masses[index];}void RopeSimulator::operate(float dt){    this->resetMassesForce();    this->solve();    this->simulate(dt);}void RopeSimulator::solve(){    for (int index = 0; index < numOfMasses - 1; ++index)//apply force of all springs    {        springs[index]->solve();//Spring with index "a" should apply its force    }    for (int index = 0; index < numOfMasses; ++index)//Start a loop to apply forces which are common for all masses    {        masses[index]->applyForce(gravitation * masses[index]->getM());//The gravitational force        masses[index]->applyForce(-masses[index]->getVel() * airFrictionConstant);//The air friction        if (masses[index]->getPos().y < groundHeight)//Forces from the ground are applied if a mass collides with the ground        {            glm::vec3 v;//A temporary Vector3Dv = masses[index]->getVel();//get the velocityv.y = 0;//omit the velocity component in y direction//The velocity in y direction is omited because we will apply a friction force to create//a sliding effect. Sliding is parallel to the ground. Velocity in y direction will be used//in the absorption effect.masses[index]->applyForce(-v * groundFrictionConstant);//ground friction force is appliedv = masses[index]->getVel();//get the velocityv.x = 0;//omit the x and z components of the velocityv.z = 0;//we will use v in the absorption effect//above, we obtained a velocity which is vertical to the ground and it will be used in//the absorption forceif (v.y < 0)//let's absorb energy only when a mass collides towards the groundmasses[index]->applyForce(-v * groundAbsorptionConstant);//the absorption force is applied//The ground shall repel a mass like a spring.//By "Vector3D(0, groundRepulsionConstant, 0)" we create a vector in the plane normal direction//with a magnitude of groundRepulsionConstant.//By (groundHeight - masses[a]->pos.y) we repel a mass as much as it crashes into the ground.glm::vec3 force = glm::vec3(0, groundRepulsionConstant, 0) *(groundHeight - masses[index]->getPos().y);masses[index]->applyForce(force);//The ground repulsion force is applied}}}void RopeSimulator::resetMassesForce()// this method will call the init() method of every mass{    for (int count = 0; count < numOfMasses; ++count)// We will init() every mass        masses[count]->init();// call init() method of the mass}

ropesimulator.cpp

#include "ropesimulator.h"RopeSimulator::RopeSimulator(){    //ctor}RopeSimulator::~RopeSimulator(){    //dtor}RopeSimulator::RopeSimulator(        int numOfMasses,//1. the number of massesfloat m,//2. weight of each massfloat springConstant,//3. how stiff the springs arefloat springLength,//4. the length that a spring does not exert any forcefloat springFrictionConstant,//5. inner friction constant of springglm::vec3 g,//6. gravitational accelerationfloat airFrictionConstant,//7. air friction constantfloat groundRepulsionConstant,//8. ground repulsion constantfloat groundFrictionConstant,//9. ground friction constantfloat groundAbsorptionConstant,//10. ground absorption constantfloat groundHeight//11. height of the ground (y position)){    this->numOfMasses = numOfMasses;    this->gravitation = g;    this->airFrictionConstant = airFrictionConstant;    this->groundFrictionConstant = groundFrictionConstant;    this->groundRepulsionConstant = groundRepulsionConstant;    this->groundAbsorptionConstant = groundAbsorptionConstant;    this->groundHeight = groundHeight;    this->masses = new Mass*[numOfMasses];    for (int count = 0; count < numOfMasses; ++count)// We will step to every pointer in the array    {        masses[count] = new Mass(m);        masses[count]->init();     }    for (int index = 0; index < numOfMasses; ++index)//To set the initial positions of masses loop with for(;;)    {        masses[index]->setPos(glm::vec3(index * springLength,0,0));    }    springs = new Spring*[numOfMasses - 1];    for (int index = 0; index < numOfMasses - 1; ++index)//to create each spring, start a loop    {        //Create the spring with index "a" by the mass with index "a" and another mass with index "a + 1".        springs[index] = new Spring(masses[index], masses[index + 1],        springConstant, springLength, springFrictionConstant);    }}void RopeSimulator::release(){    for (int count = 0; count < numOfMasses; ++count)// we will delete all of them{delete(masses[count]);masses[count] = NULL;}delete(masses);masses = NULL;    for (int index = 0; index < numOfMasses - 1; ++index)//to delete all springs, start a loop{delete(springs[index]);springs[index] = NULL;}delete(springs);springs = NULL;}void RopeSimulator::simulate(float dt){    for (int count = 0; count < numOfMasses; ++count)// We will iterate every mass    { masses[count]->simulate(dt);           // std::cout<<"mass "<<count<<" pos: "<<masses[count]->getVel().x<<masses[count]->getVel().y<<masses[count]->getVel().z<<std::endl;    }    ropeConnectionPos += ropeConnectionVel * dt;//iterate the positon of ropeConnectionPos    if (ropeConnectionPos.y < groundHeight)//ropeConnectionPos shall not go under the ground    {        ropeConnectionPos.y = groundHeight;        ropeConnectionVel.y = 0;    }    masses[0]->setPos(ropeConnectionPos);//mass with index "0" shall position at ropeConnectionPos    masses[0]->setVel(ropeConnectionVel);//the mass's velocity is set to be equal to ropeConnectionVel}void RopeSimulator::setRopeConnectionPos(glm::vec3 p){    this->ropeConnectionPos = p;}void RopeSimulator::setRopeConnectionVel(glm::vec3 v){    this->ropeConnectionVel = v;}float RopeSimulator::getGroundHeight(){    return this->groundHeight;}int RopeSimulator::getNumOfMasses(){    return this->numOfMasses;}Mass* RopeSimulator::getMass(int index){   if (index < 0 || index >= numOfMasses)// if the index is not in the arrayreturn NULL;// then return NULLreturn masses[index];}void RopeSimulator::operate(float dt){    this->resetMassesForce();    this->solve();    this->simulate(dt);}void RopeSimulator::solve(){    for (int index = 0; index < numOfMasses - 1; ++index)//apply force of all springs    {        springs[index]->solve();//Spring with index "a" should apply its force    }    for (int index = 0; index < numOfMasses; ++index)//Start a loop to apply forces which are common for all masses    {        masses[index]->applyForce(gravitation * masses[index]->getM());//The gravitational force        masses[index]->applyForce(-masses[index]->getVel() * airFrictionConstant);//The air friction        if (masses[index]->getPos().y < groundHeight)//Forces from the ground are applied if a mass collides with the ground        {            glm::vec3 v;//A temporary Vector3Dv = masses[index]->getVel();//get the velocityv.y = 0;//omit the velocity component in y direction//The velocity in y direction is omited because we will apply a friction force to create//a sliding effect. Sliding is parallel to the ground. Velocity in y direction will be used//in the absorption effect.masses[index]->applyForce(-v * groundFrictionConstant);//ground friction force is appliedv = masses[index]->getVel();//get the velocityv.x = 0;//omit the x and z components of the velocityv.z = 0;//we will use v in the absorption effect//above, we obtained a velocity which is vertical to the ground and it will be used in//the absorption forceif (v.y < 0)//let's absorb energy only when a mass collides towards the groundmasses[index]->applyForce(-v * groundAbsorptionConstant);//the absorption force is applied//The ground shall repel a mass like a spring.//By "Vector3D(0, groundRepulsionConstant, 0)" we create a vector in the plane normal direction//with a magnitude of groundRepulsionConstant.//By (groundHeight - masses[a]->pos.y) we repel a mass as much as it crashes into the ground.glm::vec3 force = glm::vec3(0, groundRepulsionConstant, 0) *(groundHeight - masses[index]->getPos().y);masses[index]->applyForce(force);//The ground repulsion force is applied}}}void RopeSimulator::resetMassesForce()// this method will call the init() method of every mass{    for (int count = 0; count < numOfMasses; ++count)// We will init() every mass        masses[count]->init();// call init() method of the mass}

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

参考

Nehe toturial Rope Physics - http://nehe.gamedev.net/tutorial/rope_physics/17006/

Siggragh Real Time Physic Class - http://www.matthiasmueller.info/realtimephysics/