Instability (slipping) of stacked boxes

Please don't post Bullet support questions here, use the above forums instead.
The_O_King
Posts: 1
Joined: Mon Mar 25, 2019 8:08 pm

Instability (slipping) of stacked boxes

Post by The_O_King »

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 slipping 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 as well as the baumgarte stabilization factors that determine the bias to apply to the impulse however that has not resolved the issue.

Collision Resolution code:

Code: Select all

                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:

Code: Select all

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/blo ... System.cpp and https://github.com/The-O-King/L3beh/blo ... mUtils.hpp