C++物理引擎构建：从DOP架构到GJK碰撞检测的实战指南

1. 项目概述为什么选择C构建物理引擎如果你正在读这篇文章大概率和我一样对游戏、动画或者机器人仿真背后的“魔法”感到着迷。屏幕上那些布料随风飘动、刚体碰撞翻滚、流体奔腾流淌的画面其核心驱动力就是一个高性能的物理引擎。市面上有成熟的引擎如PhysX、Bullet、Havok为什么我们还要从零开始用C造轮子这就像厨师学做菜最终目标未必是开餐馆但亲手从切菜、调味到颠勺走一遍你对“火候”和“食材”的理解会截然不同。自己动手构建一个物理引擎是深入理解计算机图形学、数值计算和高性能编程最硬核、也最有效的方式。C在这个领域几乎是唯一的选择这不是情怀而是现实需求。物理模拟是计算密集型的每一帧都要处理成千上万个物体的状态更新、碰撞检测和约束求解。这要求对内存布局有绝对的控制权避免GC停顿能进行极致的性能优化SIMD指令、缓存友好设计并且需要与图形API如DirectX、Vulkan进行高效、低延迟的交互。C的零成本抽象、手动内存管理以及成熟的生态如Eigen数学库、Intel TBB并行库使其成为不二之选。这个项目我们将从最基础的牛顿定律出发一步步搭建一个支持刚体动力学和基础碰撞检测的物理引擎过程中你会深刻体会到数据结构、算法和硬件架构是如何协同工作的。2. 核心架构设计与思路拆解构建一个物理引擎远不止是写几个物理公式那么简单。它更像是在设计一个精密的、实时运转的微型世界。你需要决定这个世界如何被描述数据表示时间如何流逝积分器物体如何互动碰撞与约束以及如何高效地管理这一切空间划分与内存。一个糟糕的架构设计会在你试图加入新功能比如软体或流体时让你寸步难行。2.1 面向数据的设计 vs 面向对象的设计这是第一个关键抉择。传统的面向对象OOP思维会让我们自然地想到一个RigidBody类里面封装了位置、旋转、速度、质量等属性以及Update()、ApplyForce()等方法。这在逻辑上很清晰但在性能上可能是灾难。现代CPU的性能瓶颈主要在于内存访问缓存未命中和指令并行化SIMD。OOP方式导致每个刚体的数据分散在堆内存中遍历所有刚体更新位置时CPU缓存里塞满了无关的虚函数表指针和其他成员真正需要连续处理的“位置数据”却无法被高效加载。面向数据的设计DOP是解决方案。它的核心思想是按数据的使用模式来组织内存。对于物理引擎这意味着使用SoAStructure of Arrays而非AoSArray of Structures。我们不为每个刚体创建一个结构体而是创建多个大的、连续的数组或向量分别存储所有刚体的位置、旋转、速度、质量等。// AoS (面向对象/直觉的方式缓存不友好) struct RigidBody { Vec3 position; Quat rotation; Vec3 velocity; float mass; ... }; std::vectorRigidBody bodies; // SoA (面向数据的方式缓存友好) struct RigidBodyWorld { std::vectorVec3 positions; std::vectorQuat rotations; std::vectorVec3 velocities; std::vectorfloat inverseMasses; // 常用倒数避免除法 // ... };当积分器需要更新所有位置时position velocity * dt它只需要在positions和velocities这两个连续的内存块上循环CPU的预取器可以高效工作甚至可以使用SIMD指令一次处理4个或8个数据。这是我们实现高性能的基石。2.2 组件化与系统架构基于DOP思想我们可以采用一种类似ECS实体-组件-系统但更轻量化的架构。我们将物理世界中的实体视为一个ID索引这个索引用于从各个SoA数组中查找对应的数据组件位置、速度等。然后我们定义一系列“系统”System每个系统在一个时间步内对所有相关组件执行特定操作。积分系统遍历所有具有速度和位置的实体更新它们的位置和方向。碰撞检测系统遍历所有碰撞体进行粗测Broad Phase和精测Narrow Phase生成碰撞点信息。约束求解系统处理碰撞约束、关节约束等计算并施加冲量或修正位置。力场系统应用全局力如重力、风力。这种架构的优点是高内聚、低耦合。添加一个新的力类型如磁力或新的碰撞形状如胶囊体你只需要编写对应的系统或组件数据而无需修改核心循环或其他系统。这对于引擎的长期维护和扩展至关重要。3. 数学基础与核心数据结构实现物理引擎本质上是数学在代码中的舞蹈。没有坚实的数学库一切都是空中楼阁。但记住我们追求的是正确性和性能的平衡。3.1 自定义数学库的必要性与优化虽然可以使用glm或Eigen但为了极致性能和深度控制我强烈建议实现一个精简的自定义数学库。你不需要实现所有功能只需核心部分Vec2, Vec3, Vec4用于表示点、向量、颜色。关键是要支持SIMD。对于x86/64架构可以使用__m128(SSE) 类型来同时存储和操作一个Vec4或两个Vec2。确保所有数据是16字节对齐的以发挥SIMD最大效能。Mat3, Mat4用于变换。存储时考虑列优先Column-Major这与OpenGL/DirectX的标准一致也便于SIMD操作列向量。Quaternion用于表示旋转。这是刚体旋转的核心它能避免万向节锁且插值平滑。实现四元数的乘法、共轭、标准化以及从/到旋转矩阵/欧拉角的转换。一个高性能的Vec3实现示例使用SSE intrinsics#include xmmintrin.h // SSE struct alignas(16) Vec3 { // 16字节对齐 union { struct { float x, y, z, _w; }; // _w为填充保证大小为16字节 __m128 simd; }; Vec3 operator(const Vec3 other) const { Vec3 result; result.simd _mm_add_ps(simd, other.simd); return result; } float Dot(const Vec3 other) const { __m128 mul _mm_mul_ps(simd, other.simd); // 水平相加mul [x1*x2, y1*y2, z1*z2, _] mul _mm_hadd_ps(mul, mul); // [xy, z_, xy, z_] mul _mm_hadd_ps(mul, mul); // [xyz_, ...] return _mm_cvtss_f32(mul); } // ... 其他运算 };注意使用内部函数intrinsics虽然能提升性能但会牺牲可移植性。在项目初期可以先使用普通浮点数运算保证正确性在性能分析确定瓶颈后再进行SIMD优化。3.2 刚体状态表示与存储刚体的状态需要完整描述其在世界中的运动。在我们的SoA结构中需要以下核心数组positions:Vec3世界坐标。orientations:Quaternion世界旋转。linearVelocities:Vec3线速度。angularVelocities:Vec3角速度旋转轴方向为旋转轴长度为角速度弧度值。inverseMasses:float质量的倒数。在物理计算中F ma经常被转化为a F * (1/m)使用倒数可以避免大量除法运算。inverseInertiaTensorsLocal:Mat3局部空间下惯性张量的逆。惯性张量描述了物体绕其质心旋转的难易程度使用逆矩阵同样是为了求解效率。此外还需要形状数据如包围盒、球体半径等这些可以单独管理通过刚体ID进行关联。4. 动力学核心积分器实现详解积分器是引擎的“心跳”它负责将力和速度转化为位置和旋转的变化。选择哪种积分器是在速度、稳定性和精度之间的权衡。4.1 显式欧拉法简单但不稳定这是最直观的方法用当前时刻的速度和加速度来更新下一时刻的状态。void ExplicitEuler(RigidBodyWorld world, float dt) { for (size_t i 0; i world.count; i) { // 计算合力产生的加速度 (a F * invMass) Vec3 acceleration world.forces[i] * world.inverseMasses[i]; // 更新速度 world.linearVelocities[i] acceleration * dt; // 更新位置 world.positions[i] world.linearVelocities[i] * dt; // 角速度更新类似... // 清除累积力 world.forces[i] Vec3(0.0f); } }它的优点是计算量小。但缺点非常致命能量会随着模拟而增加导致物体速度越来越快系统爆炸。这在弹簧系统或快速旋转的刚体中尤其明显。因此显式欧拉法在严肃的物理引擎中很少用于核心动力学更新。4.2 半隐式欧拉法Symplectic Euler游戏引擎的宠儿这是游戏物理引擎中最常用的方法。它的巧妙之处在于用更新后的速度去更新位置。void SymplecticEuler(RigidBodyWorld world, float dt) { for (size_t i 0; i world.count; i) { // 1. 用当前力更新速度“隐式”体现在这里其实不是真正的隐式但结构类似 world.linearVelocities[i] world.forces[i] * world.inverseMasses[i] * dt; // 2. 用*新*速度更新位置 world.positions[i] world.linearVelocities[i] * dt; // 角速度部分更新角速度然后用新角速度更新旋转 Vec3 angularAccel world.inverseInertiaTensorsWorld[i] * world.torques[i]; world.angularVelocities[i] angularAccel * dt; // 将角速度增量转换为四元数增量并更新方向 Quat deltaRot Quat::FromAngularVelocity(world.angularVelocities[i], dt); world.orientations[i] (deltaRot * world.orientations[i]).Normalized(); world.forces[i] Vec3(0.0f); world.torques[i] Vec3(0.0f); } }虽然它仍然是一阶精度但它具有辛Symplectic性质这意味着在保守力系统如重力、弹簧中它能很好地保持系统的总能量在正确值附近振荡而不是发散。这对于长期模拟的稳定性至关重要也是它被广泛采用的原因。4.3 龙格-库塔法RK4高精度需求的选择对于需要高精度模拟的场景如航天器轨道、精密机械仿真四阶龙格-库塔法是标准选择。它通过在一个时间步内进行四次采样和加权平均获得了更高的精度。// 伪代码思路定义一个函数 f(state, t) 返回状态的导数速度加速度 // 然后计算k1, k2, k3, k4最终 state state (dt/6)*(k12k22k3k4)RK4的计算量是欧拉法的四倍但精度也高得多。在游戏物理中除非有特殊的高精度需求否则半隐式欧拉法在性能和稳定性上的平衡是更优解。实操心得时间步长的处理固定时间步长Fixed Timestep是物理模拟的黄金法则。不要让物理更新的时间步长dt直接等于每一帧的渲染间隔deltaTime。因为deltaTime会波动导致模拟不一致时快时慢和数值不稳定。正确做法是在每一帧累积真实流逝的时间到一个累加器中然后以固定的physicsDt如1/60秒为单位执行一次或多次完整的物理更新直到累加器被消耗完。可能还会剩下一点时间留到下一帧。这保证了物理世界在任何硬件上都以恒定的速率演化。5. 碰撞检测系统的深度实现碰撞检测是物理引擎中最复杂的部分之一通常分为两阶段Broad Phase粗测和 Narrow Phase精测。其目标是以最快的速度排除不可能碰撞的物体对然后对可能碰撞的对进行精确的几何相交测试。5.1 粗测阶段空间划分算法当场景中有成千上万个物体时两两进行精确碰撞检测O(n²)复杂度是不可行的。粗测阶段通过空间数据结构来快速找到潜在的碰撞对。均匀网格Uniform Grid将空间划分为均匀的立方体单元格。每个物体根据其包围盒被放入一个或多个单元格。只需检查同一单元格或相邻单元格内的物体对。实现简单对于物体均匀分布的场景效率很高。但空单元格会造成内存浪费且物体大小差异大时效果不佳。四叉树/八叉树Quadtree/Octree自适应细分空间。递归地将空间划分为四个2D或八个3D子区域直到每个区域内的物体数量低于阈值。能很好地适应不均匀的物体分布内存利用率高。但树的构建和更新动态物体开销较大。包围盒层次结构BVH一种基于物体而非空间的划分方法。自底向上地将相邻的物体用更大的包围盒包裹起来形成一棵二叉树。BVH的查询效率通常很高特别适合光线追踪和静态场景。对于动态场景需要重构或更新BVH有SAH表面积启发式等优化构建算法。我的选择与实现建议对于通用的、包含动态物体的游戏物理引擎动态AABB树是一个稳健的起点。它专门为动态场景设计通过插入、删除和更新轴对齐包围盒AABB来维护一棵二叉树能高效地支持物体的移动。开源库如Box2D就使用了动态AABB树。实现它需要处理树的旋转保持平衡和节点膨胀为动态物体预留空间。5.2 精测阶段几何相交测试一旦粗测返回了潜在的碰撞对列表精测就需要进行精确的几何计算。我们需要为每种形状组合实现相交测试函数。球体 vs 球体最简单比较圆心距离和半径之和。AABB vs AABB检查在三个轴上的投影是否重叠。球体 vs AABB找到AABB上距离球心最近的点计算该点到球心的距离。胶囊体 vs 胶囊体可以转化为线段到线段的最短距离计算。凸多面体 vs 凸多面体使用分离轴定理SAT或Gilbert–Johnson–Keerthi (GJK) 算法。SAT更直观但只适用于凸体。GJK算法更通用、高效并且可以配合EPAExpanding Polytope Algorithm算法计算出穿透深度和方向即碰撞法向量和穿透深度这是求解碰撞响应所必需的。GJK算法精要GJK的核心思想是使用闵可夫斯基差Minkowski Difference。两个形状A和B的闵可夫斯基差定义为A - B {a - b | a∈A, b∈B}。关键定理是如果两个凸体相交则它们的闵可夫斯基差包含原点。GJK算法使用一个单纯形点、线段、三角形或四面体来迭代逼近闵可夫斯基差并检查原点是否被该单纯形包围。它通常只需要几次迭代就能得出是否相交的结论。注意事项数值精度问题在几何计算中浮点数误差是永恒的敌人。在判断两个物体是否“刚好接触”时使用绝对的相等或与零比较是危险的。必须使用一个容差epsilon例如1e-6。在GJK中当单纯形非常接近原点但未包含时算法可能振荡。一个实用的技巧是引入一个小的“膨胀值”或者在判断原点与单纯形的关系时使用容差。5.3 碰撞信息收集Manifold Generation检测到碰撞后我们需要生成详细的碰撞信息接触流形Contact Manifold供后续的约束求解器使用。一个流形通常包含接触点物体实际接触的一个或多个点对于面接触通常取1-4个代表性点。碰撞法线从物体A指向物体B的单位向量表示分离的方向。穿透深度物体相互嵌入的深度。对于简单形状如球体这些信息很容易计算。对于复杂形状如多边形EPA算法可以在GJK判断相交后计算出最浅的穿透深度和方向并找到一个接触点。为了稳定性通常需要生成多个接触点来防止物体旋转时抖动。6. 约束求解与碰撞响应这是让物理世界看起来“坚实”的关键一步。检测到碰撞后我们不能简单地把物体移开那样会显得很生硬。我们需要计算一个瞬间的冲量Impulse或一个位置修正来模拟真实的碰撞反弹和摩擦。6.1 基于冲量的动力学Impulse-Based Dynamics这是游戏物理中最主流的方法。它的原理是在碰撞发生的瞬间应用一个冲量J来即时改变物体的速度满足碰撞后的速度约束如恢复系数。计算相对速度在接触点处计算两个物体在该点的相对速度。计算法向冲量根据相对速度在碰撞法线方向的分量、恢复系数弹性和质量计算阻止穿透并产生反弹的冲量大小。计算切向冲量摩擦根据相对速度在切向的分量、摩擦系数和质量计算摩擦力冲量大小但大小受限于库仑摩擦定律不能超过法向力乘以摩擦系数。应用冲量将计算出的冲量分别应用到两个物体的线速度和角速度上。这种方法计算量小适合实时模拟。但它是一种“惩罚性”方法因为冲量是在穿透发生后应用的对于堆叠的物体可能会产生轻微的抖动或“过冲”。6.2 基于位置的动力学Position-Based Dynamics, PBDPBD采取了一种不同的哲学它直接操作物体的位置通过迭代求解来满足约束条件如“这两个点之间的距离应该为L”。对于碰撞约束就是“两个物体不能穿透”。在每一帧先进行预测步根据速度更新到一个临时位置。检测碰撞生成“非穿透约束”。通过一个或多个迭代同时求解所有约束包括碰撞、关节等直接修正物体的位置。根据修正后的位置和原始位置反算出新的速度。PBD的优点是非常稳定尤其对于布料、软体等堆叠和拉伸约束多的场景不易爆炸。缺点是由于是直接修正位置物理准确性如动量守恒不如冲量法严格且迭代次数会影响质量和性能。我的建议对于刚体主导的游戏物理顺序冲量求解器Sequential Impulse是工业标准Box2D, Bullet都使用它。它把碰撞和关节约束都统一表示为速度层面的约束然后在一个时间步内多次迭代通常10次左右依次求解每个约束逐渐逼近满足所有约束的解。这种方法在稳定性、性能和物理真实性之间取得了很好的平衡。6.3 约束求解器的实现要点实现一个简单的顺序冲量求解器你需要维护一个约束列表。每个约束如一个接触点提供以下信息雅可比矩阵 J描述了物体速度变化如何影响约束误差。有效质量 M_eff等于(J * M^{-1} * J^T)^{-1}其中M是质量矩阵。它综合了两个物体的质量和惯性张量表示推动这个约束的“难度”。偏差b当前约束的误差如穿透深度除以时间步长转化为速度偏差。在每次迭代中对每个约束计算当前相对速度v_rel在约束方向上的分量J*v。计算需要施加的冲量增量lambdalambda (b - J*v) / M_eff。根据冲量上下限如摩擦约束钳制lambda。将冲量J^T * lambda应用到两个物体的速度上。实操心得恢复系数与摩擦的混合处理恢复系数弹性和摩擦系数不是简单的标量。在堆叠的盒子中如果你把恢复系数设得过高盒子会永远弹跳下去。一个常见的技巧是使用“速度阈值”当碰撞的相对速度很小时将恢复系数设置为0模拟非弹性碰撞。对于摩擦将静摩擦和动摩擦分开处理能获得更真实的效果当切向相对速度低于某个阈值时使用静摩擦系数阻止运动超过阈值时使用动摩擦系数滑动摩擦。7. 性能优化与高级技巧当基础功能跑通后性能就成了下一个战场。一个未经优化的物理引擎在物体数量上百时可能就会卡顿。7.1 内存访问优化这是最立竿见影的优化。我们已经通过SoA迈出了第一步。接下来确保数据对齐使用alignas(16)或alignas(64)缓存行大小来确保数组起始地址和关键数据结构如Vec4对齐这对SIMD和缓存性能至关重要。缓存友好循环避免在紧密循环中进行随机内存访问。例如在碰撞检测的粗测阶段构建一个需要处理的物体对列表然后精测阶段顺序遍历这个列表而不是在精测函数内部再去查找物体数据。热/冷数据分离将频繁访问的数据位置、速度和偶尔访问的数据渲染用的网格句柄、用户数据分开存储。这增加了缓存中“热数据”的密度。7.2 并行化计算物理模拟中的很多任务可以并行。任务并行将不同的系统分配到不同的线程。例如一个线程处理力积分另一个线程处理碰撞检测。但需要注意任务间的依赖关系如碰撞检测必须在积分之后。数据并行这是更常用的模式。例如更新所有刚体位置的操作是彼此独立的可以很容易地用OpenMP或Intel TBB进行并行循环。#include tbb/parallel_for.h tbb::parallel_for(size_t(0), world.count, [](size_t i) { world.linearVelocities[i] world.forces[i] * world.inverseMasses[i] * dt; world.positions[i] world.linearVelocities[i] * dt; // 注意这里需要将forces[i]清零但并行写同一个变量没问题因为每个i唯一。 });SIMD并行在单个核心内使用SSE/AVX指令同时处理4个或8个刚体的数据。这需要你的SoA数据布局非常规整并且算法能够向量化。手动编写SIMD内部函数很繁琐编译器自动向量化有时效果不佳。可以考虑使用像ISPCIntel SPMD Program Compiler这样的语言它能以更直观的方式编写并行计算内核。7.3 惰性计算与脏标记不是所有数据都需要每帧更新。例如物体的世界空间惯性张量矩阵是由局部惯性张量通过旋转矩阵变换得到的。如果物体在本帧没有旋转那么这个变换就不需要重复计算。我们可以为每个物体设置一个“脏标记”Dirty Flag。只有当物体的旋转被更新时才标记其世界惯性张量为“脏”。在需要用到这个数据的地方如计算角加速度时检查标记如果为脏则重新计算并清除标记。这能节省大量不必要的计算。8. 调试、测试与常见问题排查构建物理引擎的过程就是与无数个Bug斗争的过程。没有良好的调试工具你会寸步难行。8.1 可视化调试这是最重要的调试手段。你需要能够实时看到刚体的位置和朝向绘制坐标系或简单的网格模型。速度向量从物体中心画一条箭头。碰撞形状绘制AABB、球体、胶囊体线框。接触点和法线在碰撞点画一个小点并沿着法线方向画一条短线。约束如弹簧画一条线关节画出连接点。我通常会在引擎中集成一个简单的即时模式Immediate Mode渲染器或者使用像Dear ImGui这样的库来绘制调试图形。将物理状态可视化很多问题如穿透、错误的速度方向会一目了然。8.2 单元测试与场景测试为数学库向量、矩阵、四元数运算编写全面的单元测试。使用已知的测试用例例如旋转组合、矩阵求逆等。对于碰撞检测函数创建一些特定的场景两个球体刚好相切、一个盒子完全在另一个里面、边缘碰撞等验证函数返回的结果是否正确。创建一个“测试场”场景里面包含各种典型的物理情况斜坡上的盒子、摆动的链条、堆叠的物体、高速运动的子弹等。运行这个场景观察长期模拟的稳定性。8.3 常见问题与解决实录能量爆炸物体飞走原因积分器不稳定如用了显式欧拉、时间步长太大、约束求解迭代次数不足。排查首先切换到半隐式欧拉法并大幅减小时间步长如从1/60改为1/240。如果问题解决再逐步调大时间步长并增加约束求解迭代次数找到稳定边界。物体抖动特别是堆叠时原因这是约束求解的经典难题。可能是接触点生成不稳定每帧接触点集变化、求解器迭代次数太少、或者接触容差设置不当。解决实现接触持久化Contact Persistence。将上一帧的接触点信息缓存下来如果这一帧在差不多位置还能检测到接触就复用旧的接触点信息这能大大提高稳定性。同时增加位置修正Baumgarte Stabilization或误差减少参数ERP来温和地修正穿透而不是完全依赖速度冲量。旋转异常陀螺失效或过度旋转原因四元数没有及时标准化导致误差累积或者角速度更新到旋转的积分公式有误。排查在每次更新旋转后强制调用Normalize()。确保从角速度积分到四元数增量的公式正确dq/dt 0.5 * ω * q其中ω是角速度四元数形式。性能突然下降原因可能是动态AABB树失衡或者空间划分数据结构需要重构。排查添加性能计数器监控每帧的碰撞对数量、约束求解迭代时间。实现AABB树的平衡性检查当性能下降时触发一次完整的树重构。构建物理引擎是一个庞大的工程但遵循从简到繁、逐步迭代的原则并始终将正确性和稳定性放在首位你最终会获得一个强大且高度可控的核心工具。这个过程带给你的远不止是一个可运行的库更是对实时仿真、数值计算和高性能编程的深刻理解。当你看到自己编写的代码让虚拟世界遵循着物理规律运转时那种成就感是无与伦比的。