First you choose two perpendicular vectors (or more depending on how exact you want to approximate the Coloumb friction cone) in your contact plane. Basically something like this:
// Let n be the normal of the contact patch
v32 t1, t2;
MATHvec3Perp( t1, n );
MATHvec3Cross( t2, t1, n );
Then you define the two related velocity constraints;
dC1/dt = ( v1 + w1 x r1 - v2 - w2 x r2 ) * t1
dC2/dt = ( v1 + w1 x r1 - v2 - w2 x r2 ) * t2
Note that friction constraints are non-holonome since there exists no position constraint for them. Practically this means there is no stabilization term (Baumgarte).
The tricky part now is that the friction constraints are related to the corrsponding normals through the Coloumb friction law, that means for the block LCP you get the following conditions:
-mu * lambda_n < lambda_t1 < mu * lambda_n
-mu * lambda_n < lambda_t2 < mu * lambda_n
I don't know how this gets incorporated into your Danzig solver, but Iassume you must make sure to solve the corresponding normal force (impulse) before the related friction constraints. Then you can clamp your friction forces. AFAIK the ODE has an direct solver based on Danzig's algorithm so you might look how they do it there. I personally don't like the way how it is done in the ODE in quickstep using "findex", but you can quickly find easier solutions.
HTH,
-Dirk
PS: I suggest to look up the following papers for more information:
Thesis of Kenny Erleben:
http://www.diku.dk/~kenny/thesis.pdf
Papers of Erin Catto:
http://www.gphysics.com/gdc-2005/
http://www.gphysics.com/gdc-2006/