I'm currently writing a physics simulation in 3D, that uses sequential impulses. I've managed to implement normal impulses, but I can't get friction impulses to work. As I understand I first get two tangent vectors, that are perpendicular to each other and the normal. Then I calculate the friction impulse and clamp it by total friction accumulator in the range of -FrictionCoefficient * NormalImpulse to FrictionCoefficient * NormalImpulse. I tried applying friction at the center of manifold, it seemed to work, but even with high friction coefficient values it slides a lot, so I don't consider it a solution.

Here's my code for tangents calculation:

glm::vec3 first_tangent{data.contactVelocity - collision_data.normal * glm::dot(data.contactVelocity, collision_data.normal)};
	float lat_rel_vel{ glm::length2(first_tangent) };
	glm::vec3 second_tangent{};
	if(lat_rel_vel > 0.0f) {
		first_tangent *= 1.0f / glm::sqrt(lat_rel_vel);
		second_tangent = glm::cross(first_tangent, collision_data.normal);
		second_tangent = glm::normalize(second_tangent);
	}
	else {
		GetTwoTangential(collision_data.normal, first_tangent, second_tangent);
	}

(GetTwoTangential function is the exact copy of btPlaneSpace1)

And here's my impulse calculation:

float frictional_mass1{ totalInverseMass +
			glm::dot(first_tangent, glm::cross(first->GetInverseWorldTensor() * glm::cross(data.relativeFirst, first_tangent), data.relativeFirst) +
							  glm::cross(second->GetInverseWorldTensor() * glm::cross(data.relativeSecond, first_tangent), data.relativeSecond)) };
	float j1{ -glm::dot(data.contactVelocity, first_tangent) };
	j1 /= frictional_mass1;

	const float tempFriction1{ accumulatedFrictions[currentManifold].first_tangent };
	accumulatedFrictions[currentManifold].first_tangent = glm::clamp(accumulatedFrictions[currentManifold].first_tangent + j1,
		-0.8f * accumulatedImpulses[currentManifold],
		0.8f * accumulatedImpulses[currentManifold]);
	const float deltaFriction1{ accumulatedFrictions[currentManifold].first_tangent - tempFriction1 };
	glm::vec3 fullFriction_1{deltaFriction1 * first_tangent};

	first->AppyLinearImpulse(-fullFriction_1);
	second->AppyLinearImpulse(fullFriction_1);

	first->AppyAngularImpulse(glm::cross(data.relativeFirst, -fullFriction_1));
	second->AppyAngularImpulse(glm::cross(data.relativeSecond, fullFriction_1));

	float frictional_mass2{ totalInverseMass +
		glm::dot(second_tangent, glm::cross(first->GetInverseWorldTensor() * glm::cross(data.relativeFirst, second_tangent), data.relativeFirst) +
						  glm::cross(second->GetInverseWorldTensor() * glm::cross(data.relativeSecond, second_tangent), data.relativeSecond)) };
	float j2{ -glm::dot(data.contactVelocity, second_tangent) };
	j2 /= frictional_mass1;

	const float tempFriction2{ accumulatedFrictions[currentManifold].second_tangent };
	accumulatedFrictions[currentManifold].second_tangent = glm::clamp(accumulatedFrictions[currentManifold].second_tangent + j2,
		-0.8f * accumulatedImpulses[currentManifold],
		0.8f * accumulatedImpulses[currentManifold]);
	const float deltaFriction2{ accumulatedFrictions[currentManifold].second_tangent - tempFriction2 };
	glm::vec3 fullFriction_2{deltaFriction2 * second_tangent};

	first->AppyLinearImpulse(-fullFriction_2);
	second->AppyLinearImpulse(fullFriction_2);

	first->AppyAngularImpulse(glm::cross(data.relativeFirst, -fullFriction_2));
	second->AppyAngularImpulse(glm::cross(data.relativeSecond, fullFriction_2));

(0.8f is my friction coefficient, I use it as a placeholder for now)

So here's how my collision response works:

1. Get true normal impulse( using sequential impulse solver to get normal impulses at each contact point)
2. Get friction impulse( also using sequential impulses at each contact point)
3. (Normal and friction impulses are calculated separately, first normal then friction)
4. Integrate velocities

I don't see any obvious issues with your code after a quick read. Have you verified that the tangent frame is orthonormal? Have you tried negating the impulses in case of sign error?

Here is an excerpt from my working code, which may be useful to you.

// Solver setup
{
    // Precompute linear and angular impulse vectors and compute the jacobian along the solver directions.
    c.jacobian = c.frame.computeImpulseVectors<direction>( c.r1, c.r2, velocity1, velocity2, mass1, mass2 );
    c.invJacobian = math::select( c.jacobian > 0.0f, math::reciprocal( c.jacobian ), SIMDFloat4(0.0f) );
    
    // Apply the accumulated impulse from the last frame.
    c.frame.applyImpulses<direction>( c.accumulatedImpulseX, c.accumulatedImpulseY, c.accumulatedImpulseZ, velocity1, velocity2 );
    
    // Compute the bounce speed along the normal direction due to the coefficient of restitution.
    if ( c.material.elasticity == 0.0f )
        c.bounceSpeed = 0.0f;
    else
    {
        const Float speedZ = c.frame.z.getRelativeSpeed<direction>( c.r1, c.r2, velocity1, velocity2 );
        c.bounceSpeed = -speedZ * PhysicsMaterial::getElasticityForSpeed( c.material.elasticity, speedZ,
                                                        c.material.captureVelocity, c.material.yieldVelocity );
    }
}

// Solver iteration
{
    // Compute the relative speed along the contact frame directions.
    const SIMDFloat4 r1 = c.r1;
    const SIMDFloat4 r2 = c.r2;
    const Float speedX = c.frame.x.getRelativeSpeed<direction>( r1, r2, velocity1, velocity2 );
    const Float speedY = c.frame.y.getRelativeSpeed<direction>( r1, r2, velocity1, velocity2 );
    const Float speedZ = c.frame.z.getRelativeSpeed<direction>( r1, r2, velocity1, velocity2 );
    
    // Determine an impulse that makes the relative speed the bounce speed.
    const Float jZNew = (speedZ - c.bounceSpeed)*c.invJacobian[2];
    const Float jZ = ImpulseFrame::accumulateImpulsePositive( c.accumulatedImpulseZ, jZNew );
    
    // If the magnitude of the relative speed along both tangent directions is less than
    // the static friction threshold, use static friction. Otherwise, use kinetic friction.
    const Bool staticFriction = (math::square(speedX) + math::square(speedY) < c.material.frictionThreshold2);
    const Float jtMax = c.accumulatedImpulseZ*(staticFriction ? c.material.frictionS : c.material.frictionK);
    
    // Friction clamped to pyramid
    const Float jXNew = speedX*c.invJacobian[0];
    const Float jX = ImpulseFrame::accumulateImpulseClamped( c.accumulatedImpulseX, jXNew, -jtMax, jtMax );
    const Float jYNew = speedY*c.invJacobian[1];
    const Float jY = ImpulseFrame::accumulateImpulseClamped( c.accumulatedImpulseY, jYNew, -jtMax, jtMax );
    
    // Apply the impulse along each solver frame direction.
    c.frame.applyImpulses<direction>( jX, jY, jZ, velocity1, velocity2 );
}

// Clamp impulse
T accumulateImpulseClamped( T& jAccumulated, const T& jNew, const T& jMin, const T& jMax )
{
    const T jOld = jAccumulated;
    const T jAcc = math::clamp( jOld - jNew, jMin, jMax );
    const T j = jAcc - jOld;
    jAccumulated = jAcc;
    return j;
}

@Aressera Thanks for quick answer. I fixed it by implementing solver function from qu3e and now it works perfectly. As it turned out all my calculations were right, just the order of calculation and so on were a bit wrong.