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

OpenGL进阶(十五) - 弹簧质点系统（Mass Spring Systems）简介 一个模拟变形物体最简单额方法就是将其表示为弹簧质点系统（Mass Spring Systems），因为方法简单，这个框架对于学习物理模拟和时间积分的各种技术都比较理想。一个弹簧质点包含了一系列由多个弹簧连接起来的质点，这样的系统的物理属性非常直接，模拟程序也很容易编写。 但是简单是有代价了：1.物体的行为依赖于弹簧系统的设置方法；2.很难通过调整弹簧系数来得到想要的结果；3.弹簧质点网格不能直接获取体效果。 在很多的应用中这些缺点可以忽略，在这种场合下，弹簧质点网格是最好的选择，因为够快够简单。 今天首先说明一下Mass Spring Systems的理论基础，然后来实现绳子的real-time模拟。物理模型 假设有N个质点，质量为mi，位置为xi，速度为vi , 1iN. 这些质点由一组弹簧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为帧号，带入之前的公式得：这个就是显式的欧拉解法，下一时刻的状态完全由当前状态决定。整个过程的伪代码如下：cppview plaincopy1. /initialization2. forallparticlesi3. initializexi,viandmi4. endfor5. /simulationloop6. loop7. forallparticlesi8. fifg+fcoll+f(xi,vi,xj,vj)9. endfor10. forallparticlesi11. vivi+tfi/mi12. xixi+tvi13. endfor14. displaythesystemeverynthtime15. endloopfg表示重力，fcoll表示碰撞外力。 显式积分是最简单的积分方法，不过有一个严重的问题 它只能在步长很小的时候稳定，这就意味着在两个可视帧中间需要进行多个模拟步骤，不稳定的主要原因是欧拉方法仅仅依赖于当前状态，它假设在整个步长中，外力是个常量。 假设一个场景：两个质点和一根弹簧，轻轻拉伸后松手，两个小球向中间靠拢，这时候如果时间步长过大，由于力 在一个步长中大小不变，最终结果可能导致弹簧越弹越开，最后被拉断。 更多关于显式与隐式积分，参考显式方法与隐式方法。模拟一根绳子 在布料的模拟中，大部分都使用的弹簧质点系统。 整个模拟的的过程就和伪代码中的描述的一样。环境：Ubuntu12.04 codeblocks10.04 glm9.5 代码清单：Mass类：质点类，表示质点，主要属性是质量，速度，位置mass.hcppview plaincopy1. #ifndefMASS_H2. #defineMASS_H3. #include4. #include5. classMass6. 7. public:8. Mass(floatm);9. Mass();10. voidapplyForce(glm:vec3force);11. voidinit();12. voidsimulate(floatdt);13. voidsetPos(glm:vec3p);14. voidsetVel(glm:vec3v);15. 16. glm:vec3getPos();17. glm:vec3getVel();18. floatgetM();19. protected:20. private:21. /Themassvalue22. floatm;23. /Positioninspace24. glm:vec3pos;25. /Velocity26. glm:vec3vel;27. /Forceappliedonthismassataninstance28. glm:vec3force;29. 30. ;31. 32. #endif/MASS_Hmass.cppcppview plaincopy1. #includemass.h2. 3. Mass:Mass(floatm)4. 5. /ctor6. this-m=m;7. 8. 9. Mass:Mass()10. 11. /dtor12. 13. 14. voidMass:applyForce(glm:vec3f)15. 16. this-force+=f;17. 18. 19. voidMass:init()20. 21. force.x=0;22. force.y=0;23. force.z=0;24. 25. 26. voidMass:simulate(floatdt)27. 28. vel+=(force/m)*dt;29. pos+=vel*dt;30. 31. 32. glm:vec3Mass:getPos()33. 34. returnpos;35. 36. 37. voidMass:setPos(glm:vec3p)38. 39. pos=p;40. 41. 42. voidMass:setVel(glm:vec3v)43. 44. this-vel=v;45. 46. 47. glm:vec3Mass:getVel()48. 49. returnthis-vel;50. 51. 52. 53. floatMass:getM()54. 55. returnthis-m;56. Spring类:弹簧类，一个弹簧，连着两个质点，还有弹簧的一些属性。spring.hcppview plaincopy1. #ifndefSPRING_H2. #defineSPRING_H3. #includemass.h4. #include5. classSpring6. 7. public:8. Spring();9. Spring(Mass*mass1,Mass*mass2,10. floatspringConstant,floatspringLength,11. floatfrictionConstant);12. Spring();13. 14. voidsolve();15. 16. protected:17. private:18. Mass*mass1;19. Mass*mass2;20. 21. floatspringConstant;22. floatrestLength;23. floatfrictionConstant;24. 25. ;26. 27. #endif/SPRING_Hspring.cppcppview plaincopy1. #includespring.h2. 3. Spring:Spring()4. 5. /ctor6. 7. 8. Spring:Spring()9. 10. /dtor11. 12. 13. Spring:Spring(Mass*mass1,Mass*mass2,14. floatspringConstant,floatspringLength,15. floatfrictionConstant)16. 17. this-springConstant=springConstant;/setthespringConstant18. this-restLength=springLength;/setthespringLength19. this-frictionConstant=frictionConstant;/setthefrictionConstant20. 21. this-mass1=mass1;/setmass122. this-mass2=mass2;23. 24. 25. voidSpring:solve()26. 27. glm:vec3springVector=mass1-getPos()-mass2-getPos();28. floatr=glm:length(springVector);29. 30. glm:vec3force(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);/forceisappliedtomass138. mass2-applyForce(-force);39. 最重要的Simulatorropesimulator.hcppview plaincopy1. #includeropesimulator.h2. 3. RopeSimulator:RopeSimulator()4. 5. /ctor6. 7. 8. RopeSimulator:RopeSimulator()9. 10. /dtor11. 12. 13. RopeSimulator:RopeSimulator(14. intnumOfMasses,/1.thenumberofmasses15. floatm,/2.weightofeachmass16. floatspringConstant,/3.howstiffthespringsare17. floatspringLength,/4.thelengththataspringdoesnotexertanyforce18. floatspringFrictionConstant,/5.innerfrictionconstantofspring19. glm:vec3g,/6.gravitationalacceleration20. floatairFrictionConstant,/7.airfrictionconstant21. floatgroundRepulsionConstant,/8.groundrepulsionconstant22. floatgroundFrictionConstant,/9.groundfrictionconstant23. floatgroundAbsorptionConstant,/10.groundabsorptionconstant24. floatgroundHeight/11.heightoftheground(yposition)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=newMass*numOfMasses;36. 37. for(intcount=0;countinit();41. 42. 43. for(intindex=0;indexsetPos(glm:vec3(index*springLength,0,0);46. 47. 48. 49. springs=newSpring*numOfMasses-1;50. for(intindex=0;indexnumOfMasses-1;+index)/tocreateeachspring,startaloop51. 52. /Createthespringwithindexabythemasswithindexaandanothermasswithindexa+1.53. springsindex=newSpring(massesindex,massesindex+1,54. springConstant,springLength,springFrictionConstant);55. 56. 57. 58. voidRopeSimulator:release()59. 60. for(intcount=0;countnumOfMasses;+count)/wewilldeleteallofthem61. 62. delete(massescount);63. massescount=NULL;64. 65. 66. delete(masses);67. masses=NULL;68. 69. for(intindex=0;indexnumOfMasses-1;+index)/todeleteallsprings,startaloop70. 71. delete(springsindex);72. springsindex=NULL;73. 74. 75. delete(springs);76. springs=NULL;77. 78. 79. voidRopeSimulator:simulate(floatdt)80. 81. for(intcount=0;countsimulate(dt);84. /std:coutmasscountpos:getVel().xgetVel().ygetVel().zstd:endl;85. 86. 87. ropeConnectionPos+=ropeConnectionVel*dt;/iteratethepositonofropeConnectionPos88. 89. if(ropeConnectionPos.ysetPos(ropeConnectionPos);/masswithindex0shallpositionatropeConnectionPos95. masses0-setVel(ropeConnectionVel);/themasssvelocityissettobeequaltoropeConnectionVel96. 97. 98. 99. voidRopeSimulator:setRopeConnectionPos(glm:vec3p)100. 101. this-ropeConnectionPos=p;102. 103. 104. voidRopeSimulator:setRopeConnectionVel(glm:vec3v)105. 106. this-ropeConnectionVel=v;107. 108. 109. floatRopeSimulator:getGroundHeight()110. 111. returnthis-groundHeight;112. 113. 114. intRopeSimulator:getNumOfMasses()115. 116. returnthis-numOfMasses;117. 118. 119. Mass*RopeSimulator:getMass(intindex)120. 121. if(index=numOfMasses)/iftheindexisnotinthearray122. returnNULL;/thenreturnNULL123. returnmassesindex;124. 125. 126. voidRopeSimulator:operate(floatdt)127. 128. this-resetMassesForce();129. this-solve();130. this-simulate(dt);131. 132. 133. voidRopeSimulator:solve()134. 135. for(intindex=0;indexsolve();/Springwithindexashouldapplyitsforce138. 139. 140. 141. for(intindex=0;indexapplyForce(gravitation*massesindex-getM();/Thegravitationalforce144. massesindex-applyForce(-massesindex-getVel()*airFrictionConstant);/Theairfriction145. 146. if(massesindex-getPos().ygetVel();/getthevelocity151. v.y=0;/omitthevelocitycomponentinydirection152. 153. /Thevelocityinydirectionisomitedbecausewewillapplyafrictionforcetocreate154. /aslidingeffect.Slidingisparalleltotheground.Velocityinydirectionwillbeused155. /intheabsorptioneffect.156. massesindex-applyForce(-v*groundFrictionConstant);/groundfrictionforceisapplied157. 158. v=massesindex-getVel();/getthevelocity159. v.x=0;/omitthexandzcomponentsofthevelocity160. v.z=0;/wewillusevintheabsorptioneffect161. 162. /above,weobtainedavelocitywhichisverticaltothegroundanditwillbeusedin163. /theabsorptionforce164. 165. if(v.yapplyForce(-v*groundAbsorptionConstant);/theabsorptionforceisapplied167. 168. /Thegroundshallrepelamasslikeaspring.169. /ByVector3D(0,groundRepulsionConstant,0)wecreateavectorintheplanenormaldirection170. /withamagnitudeofgroundRepulsionConstant.171. /By(groundHeight-massesa-pos.y)werepelamassasmuchasitcrashesintotheground.172. glm:vec3force=glm:vec3(0,groundRepulsionConstant,0)*173. (groundHeight-massesindex-getPos().y);174. 175. massesindex-applyForce(force);/Thegroundrepulsionforceisapplied176. 177. 178. 179. 180. 181. voidRopeSimulator:resetMassesForce()/thismethodwillcalltheinit()methodofeverymass182. 183. for(intcount=0;countinit();/callinit()methodofthemass185. ropesimulator.cppcppview plaincopy1. #includeropesimulator.h2. 3. RopeSimulator:RopeSimulator()4. 5. /ctor6. 7. 8. RopeSimulator:RopeSimulator()9. 10. /dtor11. 12. 13. RopeSimulator:RopeSimulator(14. intnumOfMasses,/1.thenumberofmasses15. floatm,/2.weightofeachmass16. floatspringConstant,/3.howstiffthespringsare17. floatspringLength,/4.thelengththataspringdoesnotexertanyforce18. floatspringFrictionConstant,/5.innerfrictionconstantofspring19. glm:vec3g,/6.gravitationalacceleration20. floatairFrictionConstant,/7.airfrictionconstant21. floatgroundRepulsionConstant,/8.groundrepulsionconstant22. floatgroundFrictionConstant,/9.groundfrictionconstant23. floatgroundAbsorptionConstant,/10.groundabsorptionconstant24. floatgroundHeight/11.heightoftheground(yposition)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=newMass*numOfMasses;36. 37. for(intcount=0;countinit();41. 42. 43. for(intindex=0;indexsetPos(glm:vec3(index*springLength,0,0);46. 47. 48. 49. springs=newSpring*numOfMasses-1;50. for(intindex=0;indexnumOfMasses-1;+index)/tocreateeachspring,startaloop51. 52. /Createthespringwithindexabythemasswithindexaandanothermasswithindexa+1.53. springsindex=newSpring(massesindex,massesindex+1,54. springConstant,springLength,springFrictionConstant);55. 56. 57. 58. voidRopeSimulator:release()59. 60. for(intcount=0;countnumOfMasses;+count)/wewilldeleteallofthem61. 62. delete(massescount);63. massescount=NULL;64. 65. 66. delete(masses);67. masses=NULL;68. 69. for(intindex=0;indexnumOfMasses-1;+index)/todeleteallsprings,startaloop70. 71. delete(springsinde