I have been working on implementing a 3D physics engine for a game engine I am developing for a school project. I have been following the work done by both Erin Catto and Box2D as well as Randy Gaul and qu3e, a sequential impulses solver.

The problem arises when two boxes are stacked on top of each other. In a demo, I have a ground plane that is totally flat and I have the two boxes fall directly on top of each other, in which they bounce a little resulting in them being slightly offset from each other when coming to "rest" on top of each other. Over time however, the top box begins to slide off the bottom box slowly, even though they should be at rest.

There isn't really any jitter, I believe this may have something to do with the normal of the contact manifold I am using for the collision with the two boxes being just slightly offset from being directly up in the y-direction when the two boxes are on top of one another (and as a result the basis vectors that I calculate for the tangent direction calculations that I use for friction are slightly offset). I am running my simulation at the moment with 4 iterations on the sequential impulse solver, but have run it with up to 15 iterations and I still see this sliding artifact.

One thing that I have noticed that is implemented in qu3e is the idea of reprojecting the friction when updating a contact point that existed in a previous frame, however I am unsure if that will make a difference with the problem that I am dealing with right now.

I have also tried changing the number of iterations in the solver, increased/decreased the friction of boxes in the scene, as well as the baumgarte stabilization factors that determine the bias to apply to the impulse however that has not resolved the issue. I have also attached a video that shows what is happening in my engine right now.

Collision Resolution code:

                float invMass1 = pc1.invMass, invMass2 = pc2.invMass;
                glm::mat3 invInertia1 = pc1.invInertia, invInertia2 = pc2.invInertia;

                // Pre Solving Computations
                float e = min(pc1.restitutionCoefficient, pc2.restitutionCoefficient);
                float mu = sqrt(pc1.friction * pc2.friction);
                for (int w = 0; w < res->numContacts; w++){
                    contactPoint cp = res->points[w];
                    res->points[w].r1 = cp.point - tc1.worldPosition;
                    res->points[w].r2 = cp.point - tc2.worldPosition;
                    glm::vec3 crossPoint1 = glm::cross(res->points[w].r1, res->normal);
                    glm::vec3 crossPoint2 = glm::cross(res->points[w].r2, res->normal);
                    float kNormal = invMass1 + invMass2;
                    kNormal += glm::dot(crossPoint1, invInertia1 * crossPoint1) + glm::dot(crossPoint2, invInertia2 * crossPoint2);
                    res->points[w].massNormal = 1.0f / kNormal;
                    res->points[w].bias = -PERCENT / deltaTime * max(0.0f, res->points[w].penetrationDist - THRESH);

                    for (int z = 0; z < 2; z++){
                        glm::vec3 crossPoint1T = glm::cross(res->points[w].r1, res->tangent[z]);
                        glm::vec3 crossPoint2T = glm::cross(res->points[w].r2, res->tangent[z]);
                        float kTangent = invMass1 + invMass2;
                        kTangent += glm::dot(crossPoint1T, invInertia1 * crossPoint1T) + glm::dot(crossPoint2T, invInertia2 * crossPoint2T);
                        res->points[w].massTangent[z] = 1.0f / kTangent;
                    }

                    //Accumulated Impulses Warm Starting
                    glm::vec3 P = res->points[w].Pn * res->normal + res->points[w].Pt[0] * res->tangent[0] + res->points[w].Pt[1] * res->tangent[1];
                    pc1.velocity -= invMass1 * P;
                    pc1.angularVelocity -= invInertia1 * glm::cross(res->points[w].r1, P);
                    pc2.velocity += invMass2 * P;
                    pc2.angularVelocity += invInertia2 * glm::cross(res->points[w].r2, P);
                }

                for (int q = 0; q < 4; q++){
                    for (int w = 0; w < res->numContacts; w++){
                        // Normal Impulse
                        glm::vec3 relativeVelocity = pc2.velocity - pc1.velocity + glm::cross(pc2.angularVelocity, res->points[w].r2) - glm::cross(pc1.angularVelocity, res->points[w].r1);
                        float velocityAlongNormal = glm::dot(relativeVelocity, res->normal);
                        float j = (-velocityAlongNormal + res->points[w].bias) * res->points[w].massNormal;
                        
                        float temp = res->points[w].Pn;
                        res->points[w].Pn = min(temp + j, 0.0f);
                        j = res->points[w].Pn - temp;

                        glm::vec3 P = j * res->normal;
                        pc1.velocity -= invMass1 * P;
                        pc1.angularVelocity -= invInertia1 * glm::cross(res->points[w].r1, P);
                        pc2.velocity += invMass2 * P;
                        pc2.angularVelocity += invInertia2 * glm::cross(res->points[w].r2, P);

                        // Friction Impulse
                        relativeVelocity = pc2.velocity - pc1.velocity + glm::cross(pc2.angularVelocity, res->points[w].r2) - glm::cross(pc1.angularVelocity, res->points[w].r1);
                        for (int z = 0; z < 2; z++){
                            float velocityAlongTangent = glm::dot(relativeVelocity, res->tangent[z]); 
                            float j = -velocityAlongTangent * res->points[w].massTangent[z];
                            float maxJ = mu * res->points[w].Pn;

                            float temp = res->points[w].Pt[z];
                            res->points[w].Pt[z] = glm::clamp(temp + j, maxJ, -maxJ); // Negate the larger one because Pn is negative
                            j = res->points[w].Pt[z] - temp;

                            glm::vec3 P = j * res->tangent[z];
                            pc1.velocity -= invMass1 * P;
                            pc1.angularVelocity -= invInertia1 * glm::cross(res->points[w].r1, P);
                            pc2.velocity += invMass2 * P;
                            pc2.angularVelocity += invInertia2 * glm::cross(res->points[w].r2, P);
                        }
                    }
                }
            }

Collision Detection Code:

void BoxVsBox(int c1, TransformComponent &tc1, int c2, TransformComponent &tc2, collisionInfo &res, World* mWorld){
    BoxColliderComponent cc1 = mWorld->getComponent<BoxColliderComponent>(c1);
    BoxColliderComponent cc2 = mWorld->getComponent<BoxColliderComponent>(c2);
    glm::mat4 orientation1 = tc1.getOrientation();
    glm::mat4 orientation2 = tc2.getOrientation();
    OBB obb1, obb2;
    obb1.AxisX = glm::vec3(orientation1[0]);
    obb1.AxisY = glm::vec3(orientation1[1]);
    obb1.AxisZ = glm::vec3(orientation1[2]);
    obb1.Half_size = cc1.halfSize;
    obb2.AxisX = glm::vec3(orientation2[0]);
    obb2.AxisY = glm::vec3(orientation2[1]);
    obb2.AxisZ = glm::vec3(orientation2[2]);
    obb2.Half_size = cc2.halfSize;

    glm::vec3 RPos = (tc2.worldPosition + cc2.offset) - (tc1.worldPosition + cc1.offset);
    glm::vec3 facePlanes[6] = { obb1.AxisX, obb1.AxisY, obb1.AxisZ, obb2.AxisX, obb2.AxisY, obb2.AxisZ };
    glm::vec3 edgePlanes[9] = { 
        glm::normalize(glm::cross(obb1.AxisX,obb2.AxisX)), 
        glm::normalize(glm::cross(obb1.AxisX,obb2.AxisY)),
        glm::normalize(glm::cross(obb1.AxisX,obb2.AxisZ)),
        glm::normalize(glm::cross(obb1.AxisY,obb2.AxisX)), 
        glm::normalize(glm::cross(obb1.AxisY,obb2.AxisY)), 
        glm::normalize(glm::cross(obb1.AxisY,obb2.AxisZ)),
        glm::normalize(glm::cross(obb1.AxisZ,obb2.AxisX)), 
        glm::normalize(glm::cross(obb1.AxisZ,obb2.AxisY)), 
        glm::normalize(glm::cross(obb1.AxisZ,obb2.AxisZ))
    };
    float faceDistances[6] = { 
        separatingDistance(RPos, facePlanes[0], obb1, obb2), 
        separatingDistance(RPos, facePlanes[1], obb1, obb2), 
        separatingDistance(RPos, facePlanes[2], obb1, obb2),
        separatingDistance(RPos, facePlanes[3], obb1, obb2),
        separatingDistance(RPos, facePlanes[4], obb1, obb2),
        separatingDistance(RPos, facePlanes[5], obb1, obb2) 
    };
    float edgeDistances[9] = {
        separatingDistance(RPos, edgePlanes[0], obb1, obb2),
        separatingDistance(RPos, edgePlanes[1], obb1, obb2),
        separatingDistance(RPos, edgePlanes[2], obb1, obb2),
        separatingDistance(RPos, edgePlanes[3], obb1, obb2),
        separatingDistance(RPos, edgePlanes[4], obb1, obb2),
        separatingDistance(RPos, edgePlanes[5], obb1, obb2),
        separatingDistance(RPos, edgePlanes[6], obb1, obb2),
        separatingDistance(RPos, edgePlanes[7], obb1, obb2),
        separatingDistance(RPos, edgePlanes[8], obb1, obb2)
    };

    float minPenetrationFace = -FLT_MAX;
    int faceMin = 7;
    bool collided = true;
    for (int x = 0; x < 6; x++){
        if (faceDistances[x] > 0){
            collided = false;
        }
        if (faceDistances[x] > minPenetrationFace){
            minPenetrationFace = faceDistances[x];
            faceMin = x;
        }
    }

    float minPenetrationEdge = -FLT_MAX;
    int edgeMin = 10;
    if (collided){
        for (int x = 0; x < 9; x++){
            if (edgeDistances[x] > 0){
                collided = false;
            }
            else if (edgeDistances[x] > minPenetrationEdge){
                minPenetrationEdge = edgeDistances[x];
                edgeMin = x;
            }
        }
    }
    res.distBetweenObjects = minPenetrationEdge > minPenetrationFace ? minPenetrationEdge : minPenetrationFace;

    if (collided){
        if (minPenetrationEdge * .95 > minPenetrationFace + .01){
            //Edge Collision
            res.normal = edgePlanes[edgeMin];
            if (glm::dot(RPos, res.normal) < 0)
                res.normal *= -1;
            
            int axisEdge1 = edgeMin / 3;
            int axisEdge2 = edgeMin % 3;
            glm::vec3 edge1[2];
            glm::vec3 edge2[2];
            getEdge(orientation1, obb1, res.normal, edge1, axisEdge1, res.points[0].contactID.inI);
            getEdge(orientation2, obb2, -res.normal, edge2, axisEdge2, res.points[0].contactID.inR);
            res.numContacts = 1;
            res.points[0].point = computeContactPointEdges(edge1, edge2);
            res.points[0].penetrationDist = abs(minPenetrationEdge);
            res.points[0].contactID.key = 2;
            res.normal *= -1;
            computeBasis(res.normal, res.tangent, res.tangent + 1);
        }
        else{
            //Face Collision
            glm::vec3 referenceFaceNormal = facePlanes[faceMin];
            int referenceFaceOpposite = false;
            if (glm::dot(RPos, referenceFaceNormal) < 0){
                referenceFaceNormal *= -1;
                referenceFaceOpposite = true;
            }

            glm::vec3 incidentFaceNormal;
            clipVertex incidentFaceVertices[4];
            clipVertex referenceFaceVertices[4];

            if (faceMin < 3){
                float min = FLT_MAX;
                bool opposite = false;
                int axis = 0;
                for (int x = 0; x < 3; x++){
                    float dotRes = glm::dot(facePlanes[3+x], referenceFaceNormal);
                    if (dotRes < min){
                        min = dotRes;
                        incidentFaceNormal = facePlanes[3+x];
                        opposite = false;
                        axis = x;
                    }
                    if (-dotRes < min){
                        min = -dotRes;
                        incidentFaceNormal = -facePlanes[3+x];
                        opposite = true;
                        axis = x;
                    }
                }
                getVerticesFromFaceNormal(incidentFaceVertices, obb2, orientation2, axis, opposite);
                getVerticesFromFaceNormal(referenceFaceVertices, obb1, orientation1, faceMin, referenceFaceOpposite);
                clipVertices(referenceFaceVertices, referenceFaceNormal, incidentFaceVertices, incidentFaceNormal, res.points, res.numContacts);
                res.normal = referenceFaceNormal;
            }
            else{
                referenceFaceNormal *= -1;
                referenceFaceOpposite = !referenceFaceOpposite;
                float min = FLT_MAX;
                bool opposite = false;
                int axis = 0;
                for (int x = 0; x < 3; x++){
                    float dotRes = glm::dot(facePlanes[x], referenceFaceNormal);
                    if (dotRes < min){
                        min = dotRes;
                        incidentFaceNormal = facePlanes[x];
                        opposite = false;
                        axis = x;
                    }
                    if (-dotRes < min){
                        min = -dotRes;
                        incidentFaceNormal = -facePlanes[x];
                        opposite = true;
                        axis = x;
                    }
                }
                getVerticesFromFaceNormal(incidentFaceVertices, obb1, orientation1, axis, opposite);
                getVerticesFromFaceNormal(referenceFaceVertices, obb2, orientation2, faceMin % 3, referenceFaceOpposite);
                clipVertices(referenceFaceVertices, referenceFaceNormal, incidentFaceVertices, incidentFaceNormal, res.points, res.numContacts);
                res.normal = referenceFaceNormal * -1.0f;
            }
            res.normal *= -1.0f;
            computeBasis(res.normal, res.tangent, res.tangent + 1);
        }
    }
}

Thank you in advance for your help! If you need any additional information please let me know too!
Additional code for the project can be found here:

https://github.com/The-O-King/L3beh/blob/master/CollisionSystem.cpp

https://github.com/The-O-King/L3beh/blob/master/CollisionSystemUtils.hpp

Bump! I would be really appreciative of some help on this, also if there's anything wrong/missing from my post do let me know! I'd be happy to clarify or fix

Thanks!

You should be able to watch your velocity and trace where the increasing velocity is coming from. It's probably coming from some cross products involved in the friction solver, and something in there might be incorrect. That is my best guess. This may also be a warm starting bug.

Things to try:

• Disable warm starting and jack up the iterations. Does the behavior change?
• Log your rigid body and solver state, so you can view difference between each frame. It is ideal to also log the deltas. This way you can track where the growing velocity comes from.
• Place the boxes very far away from the origin, and observe the behavior. If the behavior is dramatically different, you likely have messed up some cross products in the solver.

You will have to put in some time to do some fairly low level debugging to track the problem down. If you have trouble here, then my next best advice would be to take my code from qu3e and integrate it into your project and make sure it works as intended. Either this, or the reverse, where you integrate your code into qu3e and get a working baseline. When I first started learning about game physics solvers I always started with some of Erin's code as a working baseline, and would incrementally add changes or features. If something broke, I could immediately know what change caused it, which helped to pinpoint my misunderstandings.

Hi Randy, before I start I just wanted to thank you for all the tutorials (such as your own blog posts and the Tutsplus articles) you have posted over the years about creating your own engines, they have been an invaluable resource in me getting to where I am today with this project! I also really appreciate the advice  

Also apologies for the delayed response, I didn't have very much time to work on this project this past week but I have done some initial testing and was wondering if there were any thoughts about what I am seeing right now

• I have tried disabling warm starting and increasing iterations and the behavior didn't seem to change very much, still sliding
• I placed the boxes far from the origin and the simulation seemed to be more stable overall which was not what I was expecting

The biggest thing I have observed is the way that the simulation changes when the order of entities iterated over changes. When iterating from bottom to the top of a stacked box stack (floor, box1, box2) the top box is stable with each other while the middle boxes has this strange wobble; when iterating from top to bottom I see the sliding that I describe in my initial post. I have posted a video of this new behavior and I will keep digging to see why this is occurring, but if anyone has some insight I would be happy to hear it!

One thing I am wondering is whether or not I need an island solver (which I assume iterates over the boxes in an island in some optimal order), which I know is implemented in qu3e; but is an island solver necessary for basic simulation stability or should the simulation be stable with or without an island solver?

Islanding code is not necessary. You want to remove as many features as possible to isolate your problem, so adding islanding will only exacerbate your current problem.

I think you are missing to warmstart the friction impulses if I read your code correctly. Note that you need to project the old impulses onto the new tangent directions for the friction warmstarting to work correctly.

Hey there guys, thanks for the help! And Dirk, both your GDC talks on the Separating Axis Theorem extremely helpful when researching and first implementing my code, so I just wanted to thank you as well!

So I have been investigating a little further, and I've tried to simplify everything as much as possible, so I have disabled warm starting and I still observe this behavior. I did find some discussion online about projecting the old tangent impulses onto the new tangent directions for friction warm starting to work, I will definitely look at implementing that once I have a slightly more stable simulation haha, thanks

One thing that I notice that I have changed that is different from other implementations that I have seen is that in my code I do the PreSolve and the Solve steps of the sequential impulse solver immediately after detecting a collision, rather than detecting all collision pairs first and then doing the presolve step on all collision pairs and then the actual solve. Does the order in which I do the Presolve/Solve steps in relation to when I detect all collisions matter? I changed my code to test both these methods however it didn't seem to make a difference (which was my initial hunch), but just wanted a sanity check to make sure.

I also believe I might have a bug in the way that I ID contact points which results in contact points that should be the same between frames not being the same, so I am also investigating this as a potential cause of this issue

Thanks again everyone for your help! Just being able to rubber duck this problem on the forums helps me a ton, I super appreciate it

Using the Sequential Impulses (SI) technique, the ultimate goal of detecting collisions and generating contact manifolds is simply initializing contact constraints. 

I would avoid code indirections as much as possible in an initial implementation in order to keep things down to Earth.

For me PreSolve/PostSolve is a function notifying an user a contact state before/after solving the contact constraint.

Typically one performs operations in the following order.

1. Perform collision detection and contact creation on each pair of shapes. This will generate some temporary constraint solver data such as contact point and normal.
2. Notify the user the contact state before solving the contact constraints by calling PreSolve on a contact callback.
3. Initialize the contact constraints. That is, compute temporary constraint solver data such as constraint masses.
4. Solve contact constraints.
5. Notify the user the contact state after solving the contact constraints by calling PostSolve on a contact callback.

If you have problems with generating contact IDs the you can implement a different strategy. One is check if the old and new contact points are close up to some tolerance. If they are close then they are considered to be similar and you can copy the old impulses into the new contact point. Keep in mind that this is just a quick thing to implement in order to test things out. I would definitely not use it in production but it works if things don't move around much from frame to frame. For example, when simulating a stack of 5 unit boxes.