So, I am currently implementing a 3D engine with physics and so far so good, but I seem to have a problem stacking boxes. The collision detection and contact information is generated from my MPR algorithm, which does give the correct information. It works fantastic with spheres for example. Boxes work in terms of rotation on impact, etc, but the scenario where I want to stack couple of boxes, they start to jitter, even with disabled friction and restitution and I am not sure how to fix that.

Most of my implementation follows ideas from other engines, so I do have warm starting and sequential impulses, etc. So I am not sure what part of the engine is causing this behavior. Any ideas? Do I have to somehow account for multiple contact points, cause I do see face/face collision, I get 4 different contact points, which makes sense, but they are resolved with the penetration constraint gauss-seidel solver in succession, so is there some logic missing here? Do I have to somehow impulse resolve them all at the same time, like all 4 points?

I never dug too deep into physics algorithms. You might find something useful in the book Real Time Collision Detection. I know some 3rd party physics engines I've used will “put objects to sleep if they are considered not moving”. This concept might help a little bit, but it sounds like you have an actual bug. Post a video of how bad the jitter is.

real-time collision detection is great for the narrow-phase collision primitives, but it doesn't talk much about the actual rigid body simulation.

In general, you have to make sure that your system doesn't gain energy. This can happen through both small numerical errors, and through the specific penetration response/penalty functions.

The big stick to hit this problem with, is damping. Linear damping is somewhat helpful (each frame, set velocity to velocity * 0.99,) but angular damping is crucial! The reason for this is that even a small amount of penetration along the corners of a cube (which are the most likely to penetrate,) can generate significant change in angular momentum. You may want to damp angular movement quite hard (angvel = angvel * 0.97 each step, say!) or use a magnitude-based damping – the larger the angular velocity, the more it gets damped.

Another nice trick is to move penalty impulses towards center of gravity. This will change “rotational bounces” more towards “linear bounces,” so you don't get those crazy-high angular values. When adding a penalty impulse calculate its position as a blend of (penetration-position, center-of-gravity.) Even something drastic like a 50:50 blend can look convincing while being much nicer, numerically.

omg ty so much, i was missing the damping factors ahahah, it's way more stable now and can stack a decent amount before it slowly starts falling, ty so much <3