Hi everyone,

I've implemented a simple 3D rigid body physics engine using a substepping-based XPBD solver as described by Muller et al: https://matthias-research.github.io/pages/publications/PBDBodies.pdf​​. For narrowphase collision detection, I'm using the separating axis test to generate 4 point contact manifolds following Dirk Gregorius' GDC 2015 and 2013 talks: http://media.steampowered.com/apps/valve/2015/DirkGregorius_Contacts.pdf​.​

The problem I'm running into is that I don't know how to properly leverage the 4 contact points when performing the positional update to resolve a contact (Muller et al assumes one point per contact). Currently, I just iterate through all 4 contact points sequentially for each contact, but this causes simple cases like a cube falling on a flat plane to exhibit an improper torque on the bounce (the resulting rotation depends on which cube corner contact is processed first):

Ideally, I'm looking for a way to derive a single positional update from the 4 points in my contact manifold. For the video above, obviously the ideal positional update would be to the middle of the bottom face of the cube, which would resolve all 4 contacts with no rotation.

Intuitively, it seems like in general it should be possible to carefully pick a single point in between the 4 contact points, where the proper correction magnitude applied to that point ensures all 4 contacts have been resolved (with no extra work done). Unfortunately, I'm not sure how to formulate these positional constraints. There seems to be quite a bit of discussion on this forum from a few years ago for doing this for impulses, in particular the direct enumeration method used by Box2D. However it seems like this would be quite expensive in 3D (2^4 cases), and the rigid bodies XPBD paper doesn't go into quite enough detail for me to fully understand how to build the necessary LCP in the first place for positional updates. Is there any way to simplify the problem leveraging the fact that all 4 contacts have the same normal or some property of positional updates?

I'd appreciate any suggestions / pointers!

I wouldn't recommend trying to optimize the 4 contacts to a single point. That seems unlikely to be stable in practice, any small perturbation and it will become unstable. Stability relies on having 4 contacts that apply impulses which cancel out except to counteract gravity.

The most practical method for computing the contact resolution is the sequential impulse solver, which is relatively simple and used by most real-time physics engines. It involves repeatedly calculating pairwise impulses for all of the contact points (there may be multiple contacts per pair of bodies), with some clamping and warm-starting so that it converges more quickly. With this sort of method, it helps convergence and stability to apply the impulses in a random order in each iteration, so that no contact will always dominate the others (i.e. if always applied later in the iteration).

I'd recommend looking into these for more information: Bullet Physics library, Erin Catto, Erwin Coumans.

Thanks for the response. However, my solver doesn't use impulses at all, it's based on (Extended) Position Based Dynamics (first paper link in OP). My primary question here is how I can leverage the benefits of 4 point contact manifolds within the PBD framework.

You solve it exactly like in an impulse solver. You can even accumulate the pushes (lambda) and clamp the incremental pushes against the accumulated push. Pretty much the same like in an impulse solver. But there lies also the problem with this approach. While in an impulse solver the accumulated impulses converge against the actual contact forces and ultimately find an equilibrium (balance) with the forces acting on the body, the accumulated pushes go to zero (e.g. the penetration is resolved). This makes it difficult to warm-start the solver on the position level. Finally, friction is not really a position constraint and good friction is important for good contact and stacking.

The question you need to ask yourself is what problem does (or do you hope does) XPBD with rigid bodies solve for you. I use a two pass solver where I first solve on the velocity level and then resolve penetration on the position level. We shipped this in Half-Life - Alyx and it worked great. This quite a common and well understood approach for rigid body dynamics which is both fast and gets better simulation results then velocity level solving alone . Box2D does this for over a decade now. So this is really not new.

In my opinion there is not the one solver approach. Ultimately there is a pottpourri of mathematical and physical tools and you need to choose based on your use case. Say you are asked to write a cloth simulator. The choice of solver is based on the product and team your have. E.g. if you are making a League of Legends times game where each character has only a small pixel coverage you probably only need a simple ‘wiggle’ solver. For the hero selection in the same game you might want a more animation coupled approach since this is more about presentation and you want some artist control here. Finally, for virtual clothing you might take a more sophisticated simulation approach since you don't know what cloth you need to simulate and there are no animators involved in the cloth simulation.

Finally, you cannot solve the 4 contact points on block with e.g. simple direct enumeration. A manifold of four contact points is over-constrained and you would need to use e.g. SVD. This would be very expensive and probably would not add much to the overall simulation quality in a correctly implemented impulse solver. Impulse solvers work very well for contacts and stacking.

Thank you very much for the response Dirk!

The question you need to ask yourself is what problem does (or do you hope does) XPBD with rigid bodies solve for you. I use a two pass solver where I first solve on the velocity level and then resolve penetration on the position level. We shipped this in Half-Life - Alyx and it worked great. This quite a common and well understood approach for rigid body dynamics which is both fast and gets better simulation results then velocity level solving alone . Box2D does this for over a decade now. So this is really not new.

I'm working on an engine for training game-playing agents with algorithms like reinforcement learning (RL) in physics-based environments. Basically environments like OpenAI's Hide and Seek (https://openai.com/research/emergent-tool-use) but ultimately more complex.​ The primary goals for the rigid body engine are performance (for high-throughput training) and stability (because RL loves to learn to break the physics engine if it can). I chose XPBD because Muller et al's arguments around stability were convincing, and the paper does a good job of laying out all the components of the solver (contacts, joints, etc) concisely. The literature around impulse-based solvers seemed a bit more scattered when I was starting out, although if I had found Erin Catto's presentations when I started I may have made different choices. I think XPBD can provide what I need, although I'm starting to question the performance somewhat - all the narrowphase invocations for substepping are quite expensive alongside needing to continually recompute changing contact positions inside the solver loop.

Dirk Gregorius said:

You solve it exactly like in an impulse solver. You can even accumulate the pushes (lambda) and clamp the incremental pushes against the accumulated push….

Finally, you cannot solve the 4 contact points on block with e.g. simple direct enumeration. A manifold of four contact points is over-constrained and you would need to use e.g. SVD. This would be very expensive and probably would not add much to the overall simulation quality in a correctly implemented impulse solver. Impulse solvers work very well for contacts and stacking.

This makes a lot of sense, thanks. Thinking about your response more, I've realized I was incorrect about the root cause of the problem I'm seeing. While there are some minor positional / rotational errors as a result of processing the contact points independently, at the first substep where the cube hits the ground there is (correctly) barely any rotation. The actual problem is that while computing restitution an incorrect angular velocity is applied to the cube.

My solver (following Muller et al) does a second pass after the positional updates where I iterate over all the contacts and correct the velocities computed by PBD to account for restitution and dynamic friction (3.6 in the paper). This is where I'm going wrong while handling multiple points per contact: when computing restitution, the solver first cancels out the current velocity (the default velocity computed by PBD from positional deltas) in the normal direction for each contact point and then replaces it with the velocity accounting for restitution. I think the problem is that since I update the velocity after each contact point, later points in the contact effectively are trying to cancel out the normal component of the velocity from prior contact points before computing their own restitution. Perhaps I can replace this with a 2 step process where I first cancel out the PBD velocities and then apply the restitution impulses…

Solving on the position level is a non-linear problem. Conceptually you use an outer Newton loop and an inner Gauss-Seidel loop. These loops are then combined into one. When you iterate on the position level you need to update all constraints and Jacobians plus the inertia tensor when you project since they are functions of the positions. This makes things quite expensive. In practice most implementations I have seen keep Jacobians constant. This does not work very well for joints from my experience. Mathematically, I think (!), this is more like a damped Newton and I am not sure what the convergence guarantees are here for the rigid body problem.

For the problem you describe I would use a two phase solver as I described above and install block solvers where I can. Block solvers can be just matrix solves for pairwise joints or even more complex chains or hierarchies of joints. Most of this can be found in Box2D. My recommendation for you would be to start from Box2D as well and port it to 3D as needed. This will give you a clear path forward. As a nice side effect you can easily test every variation you want in this framework. Velocity level with Baumgarte stabilization, two phase solvers and position projection like XPBD. For references I suggest everything E. Catto published over the years. You can find all of on Box2D.org. You also will get an implementation of CCD and Continuous Physics (CP) with an implementation of continuous sweeps with angular motion.

HTH,

-Dirk