Hi everyone,

I have a sequential impulse solver. When I have a box collide to the ground, I notice that there is a small velocity occurring when two boxes collide. One of the boxes being ground with inverse-mass = 0, and the other being a regular box. Gravity was set to 9.81 m/s^2.

I was expecting the box's linear and angular velocity to be zero, but it just ends up being a really small value, which bothers me since I want it to be stationary. I have a feeling there is a mistake being made here, but not sure where I'm going wrong. Would appreciate a second pair of eyes to help me, or some tips on how to debug.

	f32 solveConstraint(RigidBody* body1, RigidBody* body2, vec3 i1, vec3 i2, vec3 j1, vec3 j2) {
		return 
			dot(i1, body1->linearVel) +
			dot(i2, body1->angularVel) +
			dot(j1, body2->linearVel) +
			dot(j2, body2->angularVel);
	}

	void computeAndUpdateByFriction(Contact* contact, RigidBody* body1, RigidBody* body2, vec3 tangent, f32& Pt) {
		vec3 raXt1 = cross(contact->relContactPoint1, tangent);
		vec3 rbXt1 = cross(contact->relContactPoint2, tangent);
		f32 tangentConstraint = solveConstraint(body1, body2, -tangent, -raXt1, tangent, rbXt1);

		f32 effectiveMassTangent1 = body1->inverseMass + body2->inverseMass
			+ dot(-raXt1, body1->worldInverseInertia * (-raXt1))
			+ dot(rbXt1, body2->worldInverseInertia * (rbXt1));

		f32 frictionLambdaDelta = (-tangentConstraint / effectiveMassTangent1);
		f32 Pt0 = Pt;
		Pt = glm::clamp(Pt0 + frictionLambdaDelta, -contact->contactPointMass, contact->contactPointMass);
		frictionLambdaDelta = Pt - Pt0;

		vec3 frictionVec = tangent * frictionLambdaDelta;
		body1->linearVel -= body1->inverseMass * frictionVec;
		body1->angularVel -= body1->worldInverseInertia * (raXt1 * frictionLambdaDelta);

		body2->linearVel += body2->inverseMass * frictionVec;
		body2->angularVel += body2->worldInverseInertia * (rbXt1 * frictionLambdaDelta);
	}

	void ApplyImpulse(Contact* contact, RigidBody* body1, RigidBody* body2, f32 deltaTime, int iteration) {
		// constrant contact, need to constraing bodies so the velocity and angular velocity no longer penetrate
		// Use some extra energy to push the colliders apart
		vec3 raXn = cross(contact->relContactPoint1, contact->normal);
		vec3 rbXn = cross(contact->relContactPoint2, contact->normal);

		f32 separatingSpeed = dot((-body1->linearVel - cross(body1->angularVel, contact->relContactPoint1) + body2->linearVel + cross(body2->angularVel, contact->relContactPoint2)), contact->normal);
		f32 penetrationDepth = contact->penetration;

		f32 const allowedPenetration = 0.001f;
		f32 const restitutionSlop = 0.05f;
		f32 const baumgarteTerm = 0.1f;
		f32 b = (glm::max(0.0f, penetrationDepth - allowedPenetration) * baumgarteTerm  * ( 1.0f / deltaTime));

		b += contact->restitution * glm::max(separatingSpeed - 0.05f, 0.0f);
		
		/*
		f32 constraintVal = dot(-1.0f * contact->normal, body1->linearVel) +
			dot(-1.0f * raXn, body1->angularVel) +
			dot(contact->normal, body2->linearVel) +
			dot(rbXn, body2->angularVel) + b;
			*/
		f32 constraintVal = -separatingSpeed + b;

		// calculate effective mass JM^-1J^T
		f32 effectiveMass = (body1->inverseMass) + (body2->inverseMass);
		effectiveMass += dot(-raXn, (body1->worldInverseInertia * (-raXn)));
		effectiveMass += dot(rbXn, (body2->worldInverseInertia * rbXn));
		 
		f32 deltaLambda = constraintVal / effectiveMass;
		f32 Pn0 = contact->Pn;
		contact->Pn = glm::max(0.0f, Pn0 + deltaLambda);
		deltaLambda = contact->Pn - Pn0;

		// calculate deltaV = M^-1 * J * Lambda
		vec3 contactNormalLambda = contact->normal * deltaLambda;

	//	DrawArrowC(contact->relContactPoint1, -contact->normal, vec3(deltaLambda));

		body1->linearVel -= body1->inverseMass * contactNormalLambda;
		body1->angularVel -= body1->worldInverseInertia * (raXn * deltaLambda);

		body2->linearVel += body2->inverseMass * contactNormalLambda;
		body2->angularVel += body2->worldInverseInertia * (rbXn * deltaLambda);

		// friction calculation
		computeAndUpdateByFriction(contact, body1, body2, contact->tangent1, contact->Pt1);
		computeAndUpdateByFriction(contact, body1, body2, contact->tangent2, contact->Pt2);
	}

The velocities (samples):

Body 0 :: Velocity: (-0.000166105, -1.88632e-06, 8.73014e-05), Avel: (3.46468e-06, 7.21794e-05, -4.87758e-07)

Body 0 :: Velocity: (-0.000166063, -2.20508e-06, 8.72632e-05), Avel: (3.62967e-06, 7.21843e-05, -3.53976e-07)

Body 0 :: Velocity: (-0.000166246, -2.49605e-06, 8.70451e-05), Avel: (4.30431e-06, 7.21454e-05, -1.0112e-06)

Body 0 :: Velocity: (-0.000166064, -2.05297e-06, 8.73081e-05), Avel: (3.53673e-06, 7.21999e-05, -4.3973e-07)

Little experience with equation solving, but I know the theory.

Welcome to the world of digital floating point !!

The behavior you describe looks like quite normal to me, ie it's expected. There are two reasons for this.

1. Unlike true math that you learn from books etc, floating point numbers in a computer has a finite precision. Computations with floating point numbers also have finite precision. As a result, most values in a floating point numbers are approximations. They are close to the value they should represent, but not exactly the same.

2. Solvers don't magically find the solution. Instead they start at some point and then iteratively improve the result.

Obviously, if you work with approximated values, iterating to the real answer becomes more difficult when you're more near the solution. At some point the deviations caused by having approximated answers due to point 1 could become dominant, and your iteration routine would then get steered by random noise because the computer cannot express all values exactly.

To avoid that, solvers have a precision as input. If you get an answer from the solver, the difference between the real answer and the answer that you get should be within the precision.

And that's about it. You can find much more information at the internet, one search term I'd recommend is “what every programmer should know about floating point”. The documentation of the solver may also be useful to check.

But in short:

Assume that most results in floating point are (close) approximations. They are hardly ever “the real” answer that true match derives.

This has as consequence that checking equality (a == b) should be mostly assumed to never hold, even if a solver or even true math says they should be. [[In your case, x == 0 thus.]]

Thanks Alberth, yea I figured it might be floating point issues, but I just don't know how to determine if it really is just something that's normal for an iterative solver, or something I'm doing wrong. The uncertainty is what gets me confused.

Thanks for the advice, appreciate it!