# Physics engine and squishing of stacked objects

I have a general question that seems to apply to most physics Engines (Box2D, Unity, Matter.js, …). I'm trying to make a tower of rectangular objects that are stacked on the [static/mass=infinity] ground, and the more blocks I have, the more they "squish".

I put together a test case in CodePen, which is reproduced here:

var canvas = document.getElementById("canvas");
var context = canvas.getContext("2d");
var width = canvas.width;
var height = canvas.height;

//physics vars:
var phys = Matter;
var bod = Matter.Body;
var bodies = Matter.Bodies;
var engine;
var render;
var world;
var optionsForBodies = { // Florian
restitution: 1,
};

var boxNumber = 0;

// --- TEMP ---
const keyState = {};
window.onkeydown = function (event) {
keyState[event.code] = true;
};
window.onkeyup = function (event) {
delete keyState[event.code];
};
// --- /TEMP ---

var floor = phys.Bodies.rectangle(100, 225, 200, 50, {isStatic: true});
var blist = [floor];
for (var i = 0; i < 13; i++) {
// for (var i = 0; i < 2; i++) {
var box = bodies.rectangle(100, 200 - i * 20, 56.8, 20, optionsForBodies);
blist.push(box);
boxNumber += 1;
}
}

function drawMatter(name){
var bodies = phys.Composite.allBodies(world);
for(var i = 0; i < bodies.length; i++)
{
if(!bodies[i].render.visible)
continue;
var verts = bodies[i].vertices;
if(verts.length <= 0)
continue;
context.fillStyle = "#aaa";
context.strokeStyle = "#000";
context.lineWidth = 2;

var ofsx = 0, ofsy = 0;
context.beginPath();
context.moveTo(verts[0].x + ofsx, verts[0].y + ofsy);
for(var j = 1; j < verts.length; j ++)
context.lineTo(verts[j].x + ofsx, verts[j].y + ofsy);
context.lineTo(verts[0].x + ofsx, verts[0].y + ofsy);
context.fill();
context.stroke();
}
}

function updateMatter(ms){
phys.Engine.update(engine, ms);
}

function clrCanvas(){
context.fillStyle = "#fff";
context.fillRect(0, 0, width, height);
}

var inter;

function init(){
initMatter();
clrCanvas();
requestAnimationFrame(step);
}

function initMatter(){
engine = Matter.Engine.create();

world = engine.world;
world.gravity.y = 3.0;
render = Matter.Render.create({
canvas: canvas,
engine: engine,
options: {
width: 400,
height: 300
}
});
}

var frameNo = 0;
function step() {
var SCALER = keyState.KeyA ? 3 : 1;
var SPEED = 16 * SCALER;
if (frameNo++ % SCALER === 0) {
updateMatter(SPEED);
clrCanvas();
drawMatter();
}
requestAnimationFrame(step);
}

init();
<script src="https://cdnjs.cloudflare.com/ajax/libs/matter-js/0.11.0/matter.min.js"></script>
<body>
<canvas id="canvas" width="400px" height="400px"></canvas>
</body>

To me what's happening is not obvious and I'd love an explanation in human words, since in real life if you stacked crates they would not squish inside each other, but it's a constant with EVERY physics engine, despite playing with restitution, slop, damping, density and many parameters.

Ideally I'd like those crates to be dense enough to absorb a big part of the shock (=not bounce much) of each other and never penetrate each other.

What concepts am I missing?

In the codepen above, try to select the canvas and then press the 'A' key, which reduces the framerate (increasing the dt from 16ms to 50ms) which increases the forces computed for the objects at each iteration and therefore penetrating even more each other.

Thanks!!

• Bennett Foddy gave a great talk on game physics at GDC 2015 where he describes why this happens (see the section "Collisions are resolved one pair of bodies at a time"). He talks about some mitigations like playing with the solver iteration count or fixed timestep, but ultimately this will always be present to some extent in any conventional realtime physics system. Depending on your needs, we might be able to combine colliders or add helper colliders to tame the stack - try telling us more about the intended gameplay effect. – DMGregory May 17 '17 at 14:27

Sorry, I didn't know how to mark your comment as helpful. Turned out to be, I watched the GDC talk with great interest and along with my research, I could grasp why what I'm looking to do is just gonna happen with collisions done in pairs.

Basically, think about crates who are squished one in each other, each having the same weight and mass; the crate on top has as much the right to say to the crate below it that it should move down (to not collide with it) as the crate below it will move the crate above it to not collide anymore. This "dance" could continue almost forever, and it's the core reason of why I experienced the issue.

We can apparently not solve it unless we have an infinite number of iterations (mathematically). In practice, it is solved partially by raising the number of iterations of the solver significantly. Or by avoiding this design. For iterations, something around double the number of crates is giving a much better result, but you will never have anything perfect and you've gotta make a compromise.

I played with the friction and restitution too, but expect nothing miraculous. It was enough for my game, even if that crate design unfortunately drains on the CPU way more than one'd initially think.

I’m surprised that this question has such few views and votes. Stacking is a fundamental pathological use case of physics engines that implement a simple contact solver model. It showcases how the principle behind simulating the dynamics on an individual per-contact basis breaks down when long chains or complex interactions of contact are at play.

In the typical situation of a general purpose engine and a stack of boxes, I’d analogize the situation with a chain of nodes that are connected. In both situations the forces computed between them are integrated one at a time only at their contact boundaries.

E.g. as far as the solver is concerned, the top box isn’t interacting with any other box than the second-to-top box. Therefore, in the course of one timestep, only the penetration resolution impulse (or positional correction, as necessary) from the second-to-top box can affect the top box.

So let’s imagine 5 nodes connected with constraints

T=0: * - * - * - * - *


They must be 1 unit of space between them (because for example the blocks are one unit thick). Represented here by the distance of 3 monospace characters.

A force (such as gravity, but let’s assume for the sake of easy analysis this stack is static in a weightless environment and we experience a single impulse) is experienced on the top of the stack pushing it down.

T=1: * * - * - * - *

T=2: * - * * - * - *

T=3: *  * - *  * - *

T=4: * - * - *  *  *


...Etc. Interactions become more diffused.

The point here is that the way the “speed of sound” has to travel through the material shows how it is limited by your timestep and number of iterations. This is where the realism breaks down. Indeed if you double the number of blocks stacked together (and halve their height), you have something that is more or less twice as unstable in the simulation,

In the real world, there is no such thing as a timestep or iterations, and yielding to a squishing force should not depend in any significant amount on the number of discrete structures involved (we do not see that books with a higher count of pages are proportionately more squishy, for example. Nor that books are the squishiest of all things.)

We understand that the correct interaction is that no squishing of the stack happens (provided we are modeling an incompressible block material, of course):

T=1: * - * - * - * - *


In order for it to be possible to have this “correct” computation, the engine has to realize that all of the blocks are connected to each other based on their geometric contact state, and perform a global constraint resolution on that connectivity. In our simplified node representation, node 0 and node 4 are connected through a chain of 3-character spacing constraints, so upon node 0’s receipt of an impulse, the solver must somehow determine that since node 4 will not “give” in that direction, the normal force impulse shall be delivered immediately, canceling out the squish motion.

I think that this is all quite possible to resolve in theory, and the TOI solver for the computation of the time of impact by sweeping in the box2d engine is an example of an engine taking some steps forward (though that one does not do much in the way of directly addressing massive contact interactions such as large stacks: it simply allows for detecting and adding “in-between” iterations at the time of impact, so that bullet-like fast moving objects will not freely warp through objects without interacting with them).

So this has been my approach to explaining the situation around poor stacking performance in typical rigid body simulations and how to start thinking about improving the situation.

The trouble with putting it into practice... Well, simple solvers that only resolve local contact constraints are very beneficial as they give you straightforward implementations that have nice computational scaling. You resolve your contact constraints for all of the points of contact that your collision detector computed, and you get on with the next frame or timestep.

I have not done much exhaustive brainstorming around this, though I have pondered this on and off over the years.

Consider that if you want to actually compute the correct result of resolving forces imposed on stacks or jumbles or piles of objects, it stands to reason that the entire connectivity graph based on contact for that pile of objects must be resolved each frame. This is not in and of itself prohibitive, given that the collision detection would already be required anyway for the basic contact solver, so you just additionally compute the connectivity graph based on the contacting objects.

Then, this graph (one for each group of contacting objects) must be used each step of the simulation to implement and solve a global constraint system (one for each group of contacting objects) to resolve the interacting forces.

It’s clear that the computational complexity scales rather dramatically with the dimension (given that the cardinality of this graph will likely scale with the size of the pile squared or cubed depending on if you’re doing 2D or 3D physics) and further modulated by whatever solver you are using to solve the contact constraint model that you have to have to keep things behaving as they should. Even though the runtime of such models might be subquadratic, remember that the model is based on the contact connectivity, which would be subject to almost arbitrary change each timestep, meaning that the time complexity of setting up the model in order to run it is also very much implicated here. Likely quadratic or worse, with worse yet constants.

In almost all game-type simulations, it is difficult to make a case for the importance of having a high fidelity model of just how a jumble or pile or stack of objects (which would themselves exhibit idealized material properties, e.g. that objects are totally rigid to begin with, reducing the resultant realism anyway) reacts to arbitrary disruption. The result is always highly chaotic, and the computation involved to evaluate the most correct response is always going to be rather prohibitive: A massive constraint system must be modeled out of the connectivity graph for that particular physics state, and then every force impulse that is applied to any object on that group of objects is accumulated for input to that model and evaluated. The good news is that the model can be evaluated one time per timestep, but the overhead of having to initialize the model for every single timestep (as the scene evolves) is most likely to be the performance killer.

So in conclusion, although I referred to the pairwise contact solver model as "basic", it is usually the best fit for the problem. You're really opening up a can of worms trying to move beyond it. But I do believe that under special circumstances given a particular type of scene, a more sophisticated and flexible solver can make great strides on something like stack stability.

But you'd be about as well served pursuing heuristical solutions (such as putting simulated objects to sleep etc) or design changes to address the issue.

• I can make another useful analogy, if you create a chain using links of revolving constraints in these physics simulations, you often end up with a similar situation. The “basic” iterative solver resolves the forces on each chain one at a time, so that when you put any reasonably useful tension on that chain, you begin to have to crank the timestep down a lot, leading to a computationally inefficient simulation, and a still much stretchier chain than you originally wanted. The beautiful hack around this is to add one single maximum distance constraint on the first and last links of your chain. – Steven Lu Jan 16 '19 at 3:57