I don't have a paper, but rather a

link you can look over. This gives a decent explaination, but does assume some familiarity with the subject.

The block solver is an attempt to make the iterative solver more stable. By definition the iterative solver iterates over all the contacts and solves them one by one. Naturally, if a body has two contacts (in 2D), these will be solved separately. This isn't really an issue since the solver should converge to the global solution anyway, but solving both of these at the same time can increase numeric stability and the rate of convergence.

So the block solver, rather than solving one body's contact constraints (max of 2 in 2D) one by one, it solves them together using a mini-MLCP. This is possible with 2D and 2 constraints since there's only 4 possible variations of the complimentary conditions. These are the "ifs" you see in the code.

Take a look at

Condition Numbers. In this context, the maxCondition is an number selected to determine whether we can solve (or invert) the matrix. What was noticed is that once the two contacts of a body were close enough, there tends to be a large amount of error (manifested as oscillation or large impulses for example). Using a condition number we can quickly avoid these situations and solve just one of the constraints.

William