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 masses
		float m,										//2. weight of each mass
		float springConstant,							//3. how stiff the springs are
		float springLength,								//4. the length that a spring does not exert any force
		float springFrictionConstant,					//5. inner friction constant of spring
		glm::vec3 g,							//6. gravitational acceleration
		float airFrictionConstant,						//7. air friction constant
		float groundRepulsionConstant,					//8. ground repulsion constant
		float groundFrictionConstant,					//9. ground friction constant
		float groundAbsorptionConstant,					//10. ground absorption constant
		float 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 array
			return NULL;							// then return NULL
		return 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 Vector3D

				v = masses[index]->getVel();						//get the velocity
				v.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 applied

				v = masses[index]->getVel();						//get the velocity
				v.x = 0;								//omit the x and z components of the velocity
				v.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 force

				if (v.y < 0)							//let's absorb energy only when a mass collides towards the ground
					masses[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 masses
		float m,										//2. weight of each mass
		float springConstant,							//3. how stiff the springs are
		float springLength,								//4. the length that a spring does not exert any force
		float springFrictionConstant,					//5. inner friction constant of spring
		glm::vec3 g,							//6. gravitational acceleration
		float airFrictionConstant,						//7. air friction constant
		float groundRepulsionConstant,					//8. ground repulsion constant
		float groundFrictionConstant,					//9. ground friction constant
		float groundAbsorptionConstant,					//10. ground absorption constant
		float 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 array
			return NULL;							// then return NULL
		return 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 Vector3D

				v = masses[index]->getVel();						//get the velocity
				v.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 applied

				v = masses[index]->getVel();						//get the velocity
				v.x = 0;								//omit the x and z components of the velocity
				v.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 force

				if (v.y < 0)							//let's absorb energy only when a mass collides towards the ground
					masses[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/

  • 17
    点赞
  • 65
    收藏
    觉得还不错? 一键收藏
  • 14
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 14
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值