Matrix reduction to improve the numerical stability

Hi, I am having some numerical stability issue when solving with some less perfect geometries. I tried the same example on Sfepy and I found the problem can actually by solved without convergence issue. I found the problem is when applying the Dirichlet BCs, jaxfem is just zeroing the correspoding rows and columns in dres/du or A, while Sfepy would do the matrix reduction.

I wonder if there is a particular reason not to implement the matrix reduction? I plan to give it a try, since I do need to improve the numerical stability where it is affecting my performance heavily.