Found an error: My code occasionally subdivided quadrilaterals incorrectly.

Blue: interior
Green: boundary triangles
Red: part of boundary triangles that is actually used

After fixing above problem, I get following convergence plot.

I have also made the grid symmetric wrt axes and removed the code to zero-mean the solution.