Thanks for the reply. Maybe I am mistaken, but since this is a double integral, and with Gauss-Legendre quadrature we have evaluation with the same values of each variable, this causes the denominator to be 0 at each evaluation. So even if the singularity is skipped over, then we still have a divergent result. Is it valid, as an approximation, to simply choose slightly different values, eg. I = I + w(i)*w(j)*f[x(i), 0.999...999*x(j)]? Also, are you familiar with any 2D quadrature rules that use different evaluation points intrinsically? Thanks.