简介
一个模拟变形物体最简单额方法就是将其表示为弹簧质点系统（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/
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
