简介
一个模拟变形物体最简单额方法就是将其表示为弹簧质点系统(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/