I'm having a strange problem with the friction impulses part of my impulse solver using split impulses as described by Erin Catto.

I'm using biased velocities (no energy added), and when applying all velocities at the end of the frame I sum the current velocities and the biased velocities and integrate with the sum (for linear and angular), then I clear the biased velocities.

What I do it is exactly as described by the slides, and the results are pretty good, but I'm facing with strange friction behaviours.

I'm not using contact caching yet, but I don't think it affects the sliding behaviour of the objects (only jittering and warm starting), so, each frame the contacts are cleared.

That's what I'm doing:

• Update Contact Positions (since I'm not using contact caching this doesn't need to be considered now right?);
• Test Collisions (compute local contact points, relative contact points, normal, penetration, tangent basis as described here), etc.;
• Apply Forces (integrate velocities);
• Do Warm Start (pre-step integrated);
• For each iteration... Apply Impulses (5 iterations);
• Apply Velocities (integrate positions);

Am I missing something? I'm only worried about velocites constraints for the moment.

Since code helps sometimes that's how I'm solving for lambdas (starting from 4):

4. Do Warm Start:

inline void WarmStart(float _fInvDt) {
		for (unsigned int I = 0U; I < m_ui32TotalContacts; ++I) {
			CM_CONTACT& cContact = m_cContacts[I];

			const CVector3& vPoint = cContact.vPointOnB;
			const CVector3& vNormal = cContact.vNormalB;

			//Inverse normal mass:
			float fNormalMassSum = cContact.prbRigidBody0->InvMass() + cContact.prbRigidBody1->InvMass();
			if (!cContact.prbRigidBody0->IsStatic()) {
				fNormalMassSum += (cContact.prbRigidBody0->m_mWorldInvInertia * cContact.vR0.Cross(vNormal)).Cross(cContact.vR0).Dot(vNormal);
			}
			if (!cContact.prbRigidBody1->IsStatic()) {
				fNormalMassSum += (cContact.prbRigidBody1->m_mWorldInvInertia * cContact.vR1.Cross(vNormal)).Cross(cContact.vR1).Dot(vNormal);
			}
			cContact.fNormalMass = ML_ONE / fNormalMassSum;

			for (unsigned int I = 0U; I < 2U; ++I) {
				const CVector3& vTangent = cContact.vTangent[I];
				float fTanMassSum = cContact.prbRigidBody0->InvMass() + cContact.prbRigidBody1->InvMass();
				if (!cContact.prbRigidBody0->IsStatic()) {
					fTanMassSum += (cContact.prbRigidBody0->m_mWorldInvInertia * cContact.vR0.Cross(vTangent)).Cross(cContact.vR0).Dot(vTangent);
				}
				if (!cContact.prbRigidBody1->IsStatic()) {
					fTanMassSum += (cContact.prbRigidBody1->m_mWorldInvInertia * cContact.vR1.Cross(vTangent)).Cross(cContact.vR1).Dot(vTangent);
				}
				cContact.fTangentMass[I] = ML_ONE / fTanMassSum;
			}

#define BAUMGART 0.2f
#define ALLOWED_PENETRATION 0.01f

			cContact.fBias = -BAUMGART * _fInvDt * CMathLib::Min(ML_ZERO, -cContact.fDistance + ALLOWED_PENETRATION);

			CVector3 vP = vNormal * cContact.fPN;

			for (unsigned int I = 0U; I < 2U; ++I) {
				vP += cContact.vTangent[I] * cContact.fPT[I];
			}

			if (!cContact.prbRigidBody0->IsStatic()) {
				cContact.prbRigidBody0->m_vLinVel += vP * cContact.prbRigidBody0->InvMass();
				cContact.prbRigidBody0->m_vAngVel += cContact.prbRigidBody0->m_mWorldInvInertia * cContact.vR0.Cross(vP);
			}

			if (!cContact.prbRigidBody1->IsStatic()) {
				cContact.prbRigidBody1->m_vLinVel -= vP * cContact.prbRigidBody1->InvMass();
				cContact.prbRigidBody1->m_vAngVel -= cContact.prbRigidBody1->m_mWorldInvInertia * cContact.vR1.Cross(vP);
			}
		}
	}

For each iteration... Apply Impulses (5 iterations):

CM_CONTACT& cContact = m_cContacts[I];

			const CVector3& vPoint = cContact.vPointOnB;
			const CVector3& vNormal = cContact.vNormalB;

			//Tangent:
			{
				CVector3 vVR0 = cContact.prbRigidBody0->m_vLinVel + cContact.prbRigidBody0->m_vAngVel.Cross(cContact.vR0);
				CVector3 vVR1 = cContact.prbRigidBody1->m_vLinVel + cContact.prbRigidBody1->m_vAngVel.Cross(cContact.vR1);
				CVector3 vDV = vVR0 - vVR1;

				for (unsigned int I = 0U; I < 2U; ++I) {
					const CVector3& vTangent = cContact.vTangent[I];
					float fTangentMass = cContact.fTangentMass[I];

					float fLambda = -vDV.Dot(vTangent) * fTangentMass;
					float fMaxLambda = 0.5f * cContact.fPN;

					float fPT0 = cContact.fPT[I];
					cContact.fPT[I] = CMathLib::Clamp(fPT0 + fLambda, -fMaxLambda, fMaxLambda);
					fLambda = cContact.fPT[I] - fPT0;
					CVector3 vPT = vTangent * fLambda;

					if (!cContact.prbRigidBody0->IsStatic()) {
						cContact.prbRigidBody0->m_vLinVel += vPT * cContact.prbRigidBody0->InvMass();
						cContact.prbRigidBody0->m_vAngVel += cContact.prbRigidBody0->m_mWorldInvInertia * cContact.vR0.Cross(vPT);
					}

					if (!cContact.prbRigidBody1->IsStatic()) {
						cContact.prbRigidBody1->m_vLinVel -= vPT * cContact.prbRigidBody1->InvMass();
						cContact.prbRigidBody1->m_vAngVel -= cContact.prbRigidBody1->m_mWorldInvInertia * cContact.vR1.Cross(vPT);
					}
				}
			}

			//Normal:
			{
				CVector3 vVR0 = cContact.prbRigidBody0->m_vLinVel + cContact.prbRigidBody0->m_vAngVel.Cross(cContact.vR0);
				CVector3 vVR1 = cContact.prbRigidBody1->m_vLinVel + cContact.prbRigidBody1->m_vAngVel.Cross(cContact.vR1);
				CVector3 vDV = vVR0 - vVR1;
				float fLambda = -vDV.Dot(vNormal) * cContact.fNormalMass;

				float fPN0 = cContact.fPN;
				cContact.fPN = CMathLib::Max(fPN0 + fLambda, ML_ZERO);
				fLambda = cContact.fPN - fPN0;
				CVector3 vPN = vNormal * fLambda;

				if (!cContact.prbRigidBody0->IsStatic()) {
					cContact.prbRigidBody0->m_vLinVel += vPN * cContact.prbRigidBody0->InvMass();
					cContact.prbRigidBody0->m_vAngVel += cContact.prbRigidBody0->m_mWorldInvInertia * cContact.vR0.Cross(vPN);
				}

				if (!cContact.prbRigidBody1->IsStatic()) {
					cContact.prbRigidBody1->m_vLinVel -= vPN * cContact.prbRigidBody1->InvMass();
					cContact.prbRigidBody1->m_vAngVel -= cContact.prbRigidBody1->m_mWorldInvInertia * cContact.vR1.Cross(vPN);
				}
			}

			//Biased Normal:
			{
				CVector3 vVR0B = cContact.prbRigidBody0->m_vBiasLinVel + cContact.prbRigidBody0->m_vBiasAngVel.Cross(cContact.vR0);
				CVector3 vVR1B = cContact.prbRigidBody1->m_vBiasLinVel + cContact.prbRigidBody1->m_vBiasAngVel.Cross(cContact.vR1);
				CVector3 vVNB = vVR0B - vVR1B;
				float fLambda = (-vVNB.Dot(vNormal) + cContact.fBias) * cContact.fNormalMass;
				float fPNB0 = cContact.fPNB;
				cContact.fPNB = CMathLib::Max(fPNB0 + fLambda, ML_ZERO);
				fLambda = cContact.fPNB - fPNB0;

				CVector3 vPNB = vNormal * fLambda;

				if (!cContact.prbRigidBody0->IsStatic()) {
					cContact.prbRigidBody0->m_vBiasLinVel += vPNB * cContact.prbRigidBody0->InvMass();
					cContact.prbRigidBody0->m_vBiasAngVel += cContact.prbRigidBody0->m_mWorldInvInertia * cContact.vR0.Cross(vPNB);
				}

				if (!cContact.prbRigidBody1->IsStatic()) {
					cContact.prbRigidBody1->m_vBiasLinVel -= vPNB * cContact.prbRigidBody1->InvMass();
					cContact.prbRigidBody1->m_vBiasAngVel -= cContact.prbRigidBody1->m_mWorldInvInertia * cContact.vR1.Cross(vPNB);
				}
			}

(SOLVED)

The problem was related with the two points of contacts relative to the rigid body positions (R0 and R1). I was not calculating the values with respect to the same contact point, instead I was calculating the R0 as ContactPosition - ClosestPointOnA, and R1 as ContactPosition - ClosestPointOnB.