物理引擎开发实战（一）如何把一个盒子停在平台上？
大家可能用过很多的物理引擎，但是可能自己动手写过的人不是很多。虽说文章名字是“物理引擎”，但实际上笔者这篇东西将主要针对如何自己做一个物理demo，不会加入太多框架的东西，防止demo在没有必要的地方过度复杂，影响阅读，只会在系列结束的时候稍微加入规划的引擎框架。
为什么写这个东西呢是因为笔者在年青时候意气风发也不管别人怎么搞，自己就是要往自己的引擎里面写个物理引擎。写得很差留下了残念，这次好好地学下别人怎么搞再写一遍，没有必要的话就不把它搞成引擎了。
每一章节我都会给出模拟的核心代码。同时因为这个demo主要是针对物理的模拟，其中渲染部分的实现我就不在文章里赘述了，如果大家有什么问题可以留言或者联系我.
把一个盒子停在平台上可能比你想的要困难一点。为了防止“步子迈太大，扯着蛋”，我们先列出最基础需要哪些东西，然后我们再给它增加细节：
1. 抽象刚体
2. 建立迭代循环
3. 检测碰撞
4. 碰撞处理
然后我们一块一块地说下
抽象刚体
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | class RigidBody : public PointCloud { public: std::string name; //用于表示这个刚体是不是“静态的”，不能被推动，比如我们的地板// bool isStatic = false; //质量// float mass = 1; //1除以质量// float oneDivMass = 1; //如果等于0那完全非弹性碰撞，1完全弹性碰撞// float impulseCoefficient = 1; //转动惯量，决定一个物体在xyz三个轴上面有多难转动// Matrix4 inertia; //转动惯量矩阵的逆矩阵，计算上有用// Matrix4 inertiaInverse; //摩擦系数// float linearFrictionCoefficient = 0.3; //这里是我们把刚体的一些“易变”参数封装好几份，以后有用// RigidData datas[RDI_count]; } |
这个结构体用于抽象刚体，现在demo里面每个实体Entity持有一个RigidBody
大家可以注意到它有一个RigidData的数组
对于RigidBody值得一提的是inertia(惯性)，或者我们标准点称它为inertia tensor（惯性张量）。可以搜索一下它的意义
https://en.wikipedia.org/wiki/Moment_of_inertia
数值形态的它用于2D平面
矩阵形态的它用于3D空间
简单的说：
动量等于质量乘以速度
冲量等于动量的改变量
角动量等于惯性张量乘以角速度
角冲量等于角动量的改变量
我们需要计算物体在碰撞后角速度的变化量，所以我们需要这个惯性张量。不同的物体的惯性张量不同，它是一个3x3的矩阵，和物体的质量分布有关。Demo里我们就简单点，给出这样一个例子:
1. 立方体
2. 质量均匀
3. 质心在物体中心，demo中也就等于实体的位置，本地坐标系的原点
所以我们有这样的惯性矩阵：
代码上大概是这个样子：
1 2 3 4 5 6 7 8 9 10 | //注意我们的cube认为是单位长度的 1 // float len = entity->scale.x * 1; float hei = entity->scale.y * 1; float wid = entity->scale.z * 1; float mass = entity->rigidBody->isStatic ? 99999999.0f : 1; Matrix4 cubeInertia(mass*0.083333f * (hei*hei + wid*wid),0,0,0, 0,mass*0.083333f * (len*len + wid*wid),0,0, 0,0,mass*0.083333f * (len*len + hei*hei),0, 0,0,0,1); |
可以看出这个东西是Matrix3
这里笔者用了个Matrix4这是因为懒得去再写个Matrix3
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 | class RigidData { public: //这一次模拟帧，物体被施加了多大的力// Vector3 force; //线速度// Vector3 velocity; //加速度// Vector3 acceleration; //角速度// Vector3 angularVel; //位置// Vector3 position; //朝向// Quaternion rotation; //缩放// Vector3 scale; //能量// float energy = 0; //是不是“静置”状态，主要用来省性能和增加稳定性// bool isDormant = false; Matrix4 MakeTransMatrix() const; bool CheckIfCanbeDormant(); }; |
RigidData结构体用于存放RigidBody在运算过程中的数据信息
为什么要多份以后的章节可能会涉及（在笔者没有改变主意的情况下）
然后我们更新速度类信息的方法大概是这样
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | Entity* entity = *iterABegin; RigidBody* rigidBody = entity->rigidBody; RigidData& fromData = rigidBody->datas[from]; RigidData& toData = rigidBody->datas[to]; if(!rigidBody->isStatic) { if(!fromData.isDormant) { Vector3 combinedForce = fromData.force + (sVec3AxisNegY * rigidBody->mass * G_ACC); toData.acceleration = combinedForce * rigidBody->oneDivMass; toData.velocity = fromData.velocity + (toData.acceleration * deltaTime); toData.angularVel = fromData.angularVel; } else { memcpy(&toData, &fromData, sizeof(RigidData)); } } else { memcpy(&toData, &fromData, sizeof(RigidData)); } |
合力等于各种力的叠加
加速度等于力除以质量
速度等于原来的速度加上加速度乘以时间
角速度我们暂不考虑空气阻力所以它就等于它自己
为了验证系统本身的稳定性现阶段我们并没有运用damp
然后我们在更新位置信息的时候大概是这样
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | Entity* entity = *iterABegin; RigidBody* rigidBody = entity->rigidBody; RigidData& fromData = rigidBody->datas[from]; RigidData& toData = rigidBody->datas[to]; if(!rigidBody->isStatic) { if(!fromData.isDormant) { toData.position = fromData.position + (toData.velocity * deltaTime); float angularVelVal = toData.angularVel.length(); float angularDeltaVal = angularVelVal * deltaTime; Vector3 angularVelAxis = angularDeltaVal > EPSILON_FLT ? toData.angularVel / angularVelVal : sVec3AxisY; Quaternion deltaRotQuat(angularDeltaVal, angularVelAxis); toData.rotation = deltaRotQuat * fromData.rotation; toData.scale = fromData.scale; } } |
位移我们并没有用标准的匀变速运动公式而是直接认为每个时间片里面都是匀速运动
朝向改变我们利用了四元数的叠加，不清楚的同学看这个链接
http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/arithmetic/
这里我们就不多描述数学计算的部分了
建立循环
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 | void PhysicsRoutine::Update(double deltaTime, long frame) { if(!SceneMgr::IsValid()) { return; } //初始化物理模拟，把实体的位置同步给刚体// InitRoutine(); //这里我们更新物体的受力、速度信息// UpdateVelocities(deltaTime, RDI_real, RDI_real); //更新当前碰撞信息// UpdateManifolds(RDI_real); //用于增加稳定性// WarmStart(deltaTime, RDI_real); //处理所有的碰撞，迭代多次，求得稳定的结果// for(int i=0; i < 10; i++) { HandleContacts(deltaTime, RDI_real, RDI_real); } //更新物体的位置信息// UpdatePosAndRots(deltaTime, RDI_real, RDI_real); //把物体的位置信息同步给实体，用于绘制// Finalize(deltaTime, RDI_real, RDI_real); } |
更新物体的运动信息
更新物体的碰撞信息
搞了个WarmStart的事情
迭代多次来处理碰撞
更新物体的位置信息
把物体的位置信息刷给实体
这里大家可以先盲从笔者的循环实现，在后期的demo中，这个循环可能有更多的信息加入。大家可以暂时先忽略UpdateManifolds的实现和WarnStart的实现，大致理解为：
1. Manifold用来跨帧记录谁和谁碰成什么样子
2. WarmStart是一种工程手段，用于“帧一致性”高的时候帮助碰撞处理更快的得到正确的结果。“帧一致性”是指在时间间隔较小的时候，两帧之间的数据没有突变，这种理解带来了很多的算法，比如verlet也是其中之一
3. 这两个东西都是用来增加稳定性的
这里我们展示一下Manifold的样子
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 | #define CONTACT_POINT_COUNT 5 #define CONTACT_DRIFTING_THRESHOLD 0.004f class ContactPoint { public: //A物体的碰撞点在其局部坐标系下的位置// Vector3 localWittnessPointA; //B物体的碰撞点在其局部坐标系下的位置// Vector3 localWittnessPointB; //A物体的碰撞点在在世界坐标系下的位置// Vector3 globalWittnessPointA; //B物体的碰撞点在在世界坐标系下的位置// Vector3 globalWittnessPointB; //碰撞法线，B指向A// Vector3 normal; //碰撞法线，A指向B// Vector3 negativeNormal; //碰撞切线1，摩擦力的方向// Vector3 tangent; //碰撞切线2，摩擦力2的方向// Vector3 tangent2; //碰撞深度// float penetrationDistance; //法线方向上的总碰撞冲量// float totalNormalImpulse; //切线1方向上的总碰撞冲量// float totalTangentImpulse; //切线2方向上的总碰撞冲量// float totalTangentImpulse2; }; class ContactManifold { public: //实体A// Entity* entityA; //实体B// Entity* entityB; //碰撞点们// ContactPoint contactPoints[CONTACT_POINT_COUNT]; int contactPointCount = 0; |
可以看出ContactManifold就是用来记录EntityA在contactPoints这些地方和EntityB碰在一起了。
基于冲量的刚体物理模拟
好，在大致介绍完数据结构后是时候说下我们到底怎么搞刚体模拟了。
首先我们是基于冲量的模拟，以防大家忘记了冲量的意思，我这里简单描述下：
是动量的改变量
是速度改变的直接原因
其实这里就解释了为什么我们不去基于力或者其它东西来做模拟，因为“简单”,如果我们使用基于力的实现，那我们可能需要以下公式：
可以看出，如果基于力或者加速度，我们还需要时间参数
而冲量公式我们是直接得出速度，省了些麻烦
那好，基于冲量了我们该怎么模拟呢？这个网页有比较通俗易懂的描述https://en.wikipedia.org/wiki/Collision_response
这里我做个更加白话的总结：
1. 碰撞发生过后无非是物体的线速度和角速度被改变了
2. 线速度改变的方式是这个公式描述的
V代表线速度
J代表冲量
M代表质量
n代表法向量
3. 角速度的改变方式是这个公式描述的
代表角速度
J代表冲量
I代表惯性张量
r代表从质心到碰撞点的向量
n代表法向量
4. 我们还需要相对速度
Vr代表相对速度
V2V1各代表物体2物体1的线速度
21各代表物体2物体1的角速度
5. 才能求出冲量
Vr代表相对速度
m代表物体的质量
n代表法向量
I代表惯性张量
r代表质心到碰撞点的向量
6. 现在问题来了：挖掘机到底哪家…不不不这个r该怎么得到？
这就带来了我们的新问题“碰撞检测”
碰撞检测
用了这么大个标题来说碰撞检测，说明东西还是比较多的，这里笔者简单说下。基于之前的需求，我们需要碰撞检测回答两个问题：
1. 两个物体间有没有碰撞
2. 这两个物体是在什么位置发生碰撞的
这两个问题回答起来其实并不是很容易，我们先说个简单的：
球，大家都知道球的碰撞是计算他们的球心的距离是不是小于两个球半径的合。碰撞的点一定是他们球心连线在半径位置上的点。很简单没有问题。
但是如果它是个2D或者3D的方块呢？或者说是任意形状呢？我们先说方块，难道我们去for一定的步长检测CubeA的一个点是否在CubeB中吗？先不问你怎么判断一个点是不是在CubeB中，这个精度是不是其实就不可行？
这里我想说有算法解决这个问题，而且我非常喜欢这个算法，这个算法的名字叫GJK，它可以回答第一个问题：两个物体间有没有碰撞？
GJK的信息大家可以看下这个网页http://www.dyn4j.org/2010/04/gjk-gilbert-johnson-keerthi/。这里我概括一下它的思想，只支持凸包，GJK算法基于一种叫做明可夫斯基空间差的概念，别被名字吓到，其实就是两个物体包含的所有点的位置相减构成的点云。GJK的美丽之处就在于它不去检测一个物体的点是否在另一个物体中而是去判断这个点云里面有没有包含原点，是的，就是有没有包含原点，因为如果两个物体有交叠，那他们肯定有共有点，一减就等于vec3(0,0,0)即原点。然后它优雅的地方还在于不是去穷举两个物体的点集来相减，而是尝试在这个空间差的点云中去找一个“最简形（2D是三角形，3D是4面体）”，如果我能证明这个最简形能框住原点那我就证明了这个云包含了原点进而证明了两物体相交。
这个算法的实现细节如果大家想讨论的话我们专门做个专题来说。这里就不冲淡主题了。笔者在章节结束附上一个简单实现。
好，我们知道物体有没有碰撞了，然后呢？我们还差一个问题要回答：这两个物体是在什么位置发生碰撞的？这个时候GJK就不够了，就需要它的姐妹算法EPA，大家可以在这个网页了解一下它的信息http://www.dyn4j.org/2010/05/epa-expanding-polytope-algorithm/，这个算法的3D版本笔者年轻时并没能把它实现健壮了，还是比较复杂的。这次怕把做demo的时间拖太长就没有自己来实现了，demo全部完成之后我们再回来挑战。这里我们就说下它的思想，它在GJK发现碰撞的情况下，从GJK最后包含原点的最简形出发，反向找到在点云上最接近的边界，这个地方就被认为是首次碰撞发生的地方。
大家可以观察下RigidBody继承的PointCloud其实就是为了它们的核心方法Support(Vector3 direction);实现的
1 2 3 4 5 | class PointCloud { public: virtual Vector3 getFarthestVectAtDir(const Vector3& dir, const Matrix4* trans) const = 0; }; |
由于这两个算法必须在实现上配套，不然数据结构的差异将带来不必要的难度。这里笔者就很猥琐的偷窃了Bullet引擎的实现，在里面做了一些修改让它能跑在这边的数据结构上。请大家先别纠结这个没有牙齿的行为，我们先回到物理模拟的重点上来。
碰撞处理
现在我们有r了对吧，不知道大家还记得我们的r不？这个时候我们终于凑出了需要的所有变量，可以进行碰撞冲量计算了：
较多麻袋！我们是不是忽略了一个重要的问题：摩擦力呢？好，我们现在顺势给出摩擦力的处理方法：“Coulomb friction model”。这个网页上也给出具体的描述https://en.wikipedia.org/wiki/Collision_response
简单说：用这个模型我们把摩擦力抽象成为了切线方向上的冲量，然后就可以直接套用冲量公式，只是把里面的法线改为切线。
t代表tangent即切线
怎么突然感觉我们就有了所有需要的信息了？是不是感觉自己可以动手写了？是不是好像突然世界就美丽了？
貌似我们可以按照这样的方式模拟：
1. 更新物体速度数据
2. 检测有没有碰撞，有得出碰撞信息
3. 对每个碰撞状态处理法线和切线的冲量变化
经过笔者2011年的实践后确实证明了这个方法可行！
吗？
非也！第一版就否定了, too young too naïve
如果我们严格地按照这个思路来做的话你会发现你的盒子在地面上抖得就和猫王一样，然后要不了多久就沉到地下去了。为什么？
1. 一个碰撞点是不够的
一个稳定的支撑同一时间至少需要3个支撑点，不然EPA算法在每一帧得出一个碰撞点并处理，就会给这个盒子带来不断的旋转量，从而导致抖动
2. 同一帧提取多个碰撞点并不好，这个idea在笔者以前写的时候用过，bullet引擎也用过，那个时候叫做perturb，基本思想是在发现出现碰撞的时候，将这个实体向几个方向倾斜一下尽量求出多个不同的碰撞点来构成支撑，其实这种伪造碰撞点的方式效果并不太好
3. 这个时候我们需要的是ContactManifold，并且以以下规则更新它里面存放的碰撞点信息：
a) 每一帧只提取一个碰撞点，然后尝试将它添加进入这个集合
b) 一个集合里面最多只保留4个碰撞点
c) 当超过4个点时我们根据以下规则决定碰撞点去留
i. 保留撞入深度最深的碰撞点
ii. 根据碰撞点在AB两个物体上的距离是否仍然在对穿或者接触决定是否需要删除
iii. 选取构成四边形面积最大的碰撞点
注意这种筛选碰撞点的思想来源于bullet。
其实最经典的思路是：
a) 保留撞入深度最深的碰撞点
b) 找一个距离这个点最远的点
c) 找一个距离这两个点连线最远的点
d) 找一个距离这个三角形距离最远的点，如果在内部就丢弃
第二种方式其实原理上应该效果更好，但是带来了两个比较耗的求距离的计算，没有bullet的这个简化算法高效。放在demo里的话我还得跟大家讲怎么算距离（好麻烦），所以果断偷窃Bullet的算法：
这样我们在几个物理帧之内就可以得到一个相对稳定的结果了。
4. 对于缓慢下沉，还需要解决浮点数精度带来的问题，因为我们每次都是在发现碰撞的时候才来处理它，那意味着我们每次处理它的时候只要碰撞带来的冲量不稍微大于能让它停止继续穿入的量的话它就会因为0.99999f乘以value的原因不断因误差下沉。那如何做到“刚好稍微大于能让它继续穿入的量”呢？方法有很多，比如工程代码中的方式其实是在穿入距离超过一定量的时候，把穿入的距离按一定比例反馈给冲量，这种方式叫做“Baumgarte”，只运用它的话在非堆叠的时候效果不错，但是堆叠就不够了，当然还有其它在处理碰撞点上效果更好一些的方法，不过那是后话，我们下一章再说
1 2 3 4 5 | float biasPenetrationDepth = 0.0; if ((-contact.penetrationDistance) > PENETRATION_TOLERANCE && deltaTime > 0) { biasPenetrationDepth = 0.8f * (PENETRATION_VELOCITY_BIAS / deltaTime) * fmax(0.0f, float((-contact.penetrationDistance) - PENETRATION_TOLERANCE)); } |
10. 然后我们需要对运用于物体的冲量进行修正。为什么？因为并不是每一次冲量计算都会给我带来把两个物体分离的结果。如果我们不管这些冲量是否大于0都去运用的话，很可能导致模拟让物体反而互相贴近；而如果我们完全忽略这些因为物体正在离开另一物体而带来的“修正性负冲量”的话，无疑我们的模拟就放大了不稳定，这个时候Erin Catto大神在他的sequential impulse中提到，我们应该累计每个点的冲量结果，把它们累加起来，让它不能小于0，然后每次运用相对于这个累加的冲量的该变量的话就兼有了两边的好处。
1 2 3 | float newImpulse = fmax(contact.totalNormalImpulse + collisionImpulseVal, 0.0f); collisionImpulseVal = newImpulse - contact.totalNormalImpulse; contact.totalNormalImpulse = newImpulse; |
14. 然而这样还是不够的，我们还需要热启动“WarmStart”，也是Erin Catto提出的用于增加物理系统稳定性的改进。思想是，我们把这个累计下来的上一帧的总冲量结果在在开始新帧的碰撞处理前提前在物体上运用一遍，这样能显著增加系统的稳定，在多个物体堆叠的时候效果尤为明显。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 | void PhysicsRoutine::WarmContact(RigidDataIndex to, ContactManifold* manifold, ContactPoint& contact, double deltaTime) { Entity* entityA = manifold->entityA; RigidBody* rigidBodyA = entityA->rigidBody; RigidData& dataA = rigidBodyA->datas[to]; Entity* entityB = manifold->entityB; RigidBody* rigidBodyB = entityB->rigidBody; RigidData& dataB = rigidBodyB->datas[to]; //let normal point from a -> b Vector3 normal = contact.negativeNormal; Vector3 relativePosA = contact.globalWittnessPointA - dataA.position; Vector3 relativePosB = contact.globalWittnessPointB - dataB.position; //-- collision impulse begin ----------------------------------------------------------------------------------------------------------------------------------------------------// float collisionImpulseVal = contact.totalNormalImpulse; Vector3 collisionImpulse = normal * collisionImpulseVal; if(!rigidBodyA->isStatic) { dataA.velocity = dataA.velocity - (collisionImpulse * rigidBodyA->oneDivMass); dataA.angularVel = dataA.angularVel - (collisionImpulseVal * rigidBodyA->inertiaInverse.transform3x3(relativePosA * normal)); } if(!rigidBodyB->isStatic) { dataB.velocity = dataB.velocity + (collisionImpulse * rigidBodyB->oneDivMass); dataB.angularVel = dataB.angularVel + (collisionImpulseVal * rigidBodyB->inertiaInverse.transform3x3(relativePosB * normal)); } //-- collision impulse end ----------------------------------------------------------------------------------------------------------------------------------------------------// //-- friction impulse begin ----------------------------------------------------------------------------------------------------------------------------------------------------// Vector3 tangent = contact.tangent; float frictionImpulseVal = contact.totalTangentImpulse; Vector3 frictionImpulse = frictionImpulseVal * tangent; if(!rigidBodyA->isStatic) { dataA.velocity = dataA.velocity - frictionImpulse * rigidBodyA->oneDivMass; dataA.angularVel = dataA.angularVel - frictionImpulseVal * rigidBodyA->inertiaInverse.transform3x3(relativePosA * tangent); } if(!rigidBodyB->isStatic) { dataB.velocity = dataB.velocity + frictionImpulse * rigidBodyB->oneDivMass; dataB.angularVel = dataB.angularVel + frictionImpulseVal * rigidBodyB->inertiaInverse.transform3x3(relativePosB * tangent); } //-- friction impulse end ----------------------------------------------------------------------------------------------------------------------------------------------------// } |
我们为什么可以这么搞？还记得之前提到的“帧一致性”吗?
66. 摩擦力也需要clamp一下来防止为系统带来额外的能量
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 | //-- collision impulse begin --------------------------------------------------------------------------------// float biasPenetrationDepth = 0.0; if ((-contact.penetrationDistance) > PENETRATION_TOLERANCE && deltaTime > 0) { biasPenetrationDepth = 0.8f * (PENETRATION_VELOCITY_BIAS / deltaTime) * fmax(0.0f, float((-contact.penetrationDistance) - PENETRATION_TOLERANCE)); } float collisionImpulseValA = (relativePosA * normal).dot(rigidBodyA->inertiaInverse.transform3x3(relativePosA * normal)); float collisionImpulseValB = (relativePosB * normal).dot(rigidBodyB->inertiaInverse.transform3x3(relativePosB * normal)); float collisionImpulseVal = -(1.000f + combinedImpulseCoefficient + biasPenetrationDepth) * relativeVelocity.dot(normal) / (rigidBodyA->oneDivMass + rigidBodyB->oneDivMass + collisionImpulseValA + collisionImpulseValB); float newImpulse = fmax(contact.totalNormalImpulse + collisionImpulseVal, 0.0f); collisionImpulseVal = newImpulse - contact.totalNormalImpulse; contact.totalNormalImpulse = newImpulse; Vector3 collisionImpulse = normal * collisionImpulseVal; if(!rigidBodyA->isStatic) { dataA.velocity = dataA.velocity - (collisionImpulse * rigidBodyA->oneDivMass); dataA.angularVel = dataA.angularVel - (collisionImpulseVal * rigidBodyA->inertiaInverse.transform3x3(relativePosA * normal)); } if(!rigidBodyB->isStatic) { dataB.velocity = dataB.velocity + (collisionImpulse * rigidBodyB->oneDivMass); dataB.angularVel = dataB.angularVel + (collisionImpulseVal * rigidBodyB->inertiaInverse.transform3x3(relativePosB * normal)); } //-- collision impulse end -----------------------------------------------------------------------------------// |
法线方向上的冲量
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | //-- friction impulse begin ----------------------------------------------------------------------// Vector3 relativeVelocityTangent = relativeVelocity - (relativeVelocity.dot(normal) * normal); float relativeVelocityOnTangent = relativeVelocityTangent.normalizeF(); contact.tangent = relativeVelocityTangent; Vector3 tangent = relativeVelocityTangent; float frictionImpulseValA = (relativePosA * tangent).dot(rigidBodyA->inertiaInverse.transform3x3(relativePosA * tangent)); float frictionImpulseValB = (relativePosB * tangent).dot(rigidBodyB->inertiaInverse.transform3x3(relativePosB * tangent)); float frictionImpulseVal = -1.000f * relativeVelocityOnTangent / (rigidBodyA->oneDivMass + rigidBodyB->oneDivMass + frictionImpulseValA + frictionImpulseValB); frictionImpulseVal = fmax(fabs(combinedFrictionCoefficient * collisionImpulseVal), fabs(frictionImpulseVal)) * (frictionImpulseVal > 0 ? 1.0f : -1.0f); Vector3 frictionImpulse = frictionImpulseVal * tangent; if(!rigidBodyA->isStatic) { dataA.velocity = dataA.velocity - frictionImpulse * rigidBodyA->oneDivMass; dataA.angularVel = dataA.angularVel - frictionImpulseVal * rigidBodyA->inertiaInverse.transform3x3(relativePosA * tangent); } if(!rigidBodyB->isStatic) { dataB.velocity = dataB.velocity + frictionImpulse * rigidBodyB->oneDivMass; dataB.angularVel = dataB.angularVel + frictionImpulseVal * rigidBodyB->inertiaInverse.transform3x3(relativePosB * tangent); } //-- friction impulse end -------------------------------------------------------------------// |
摩擦力带来的冲量
这个时候我们的模拟就一般能看看了：
可以观察到还有一些穿透
Yeah!
大功告成了！
天下太平了！
吗？
非也…万里长征第一步都没走完
不过那是后话了，我们下一章再说
下一章叫
如何把一堆盒子停在平台上