I am looking for explanations of algorithms to adaptively sample a function of two variables $f(x,y)$, in a given domain $x_0\le x \le x_1$, $y_0\le y \le y_1$. Intuitively, I want to sample more densely regions were $f(x,y)$ has more variation. I know almost nothing a priori about the function, except that inside the domain $x_0\le x \le x_1$, $y_0\le y \le y_1$, there maybe regions where it is infinite (that I obviously don't want to sample, so the algorithm should detect that too), and that everywhere else, the function is continuous. These singular regions are "smooth sets".
Looks like you probably would like to steal some ideas from this discussion about implementing something similar in Julia. It's hard to cleanly summarize since there's a lot of different details, but essentially it's an algorithm which does a simple grid and then adaptively adds points recursively in the regions where the linear approximation is worse. It's like a course grid adaptivity in finite element methods where you take more grid points and see how far your interpolation was off at the grid points, and if it was a lot, you accept the new grid points and recursively re-apply the coarsening operation. Then details ensue, but it's all open source with a related discussion so you can check it out all the gory details. Taking that from 1D to 2D likely involves new code but the same principles.
the answer maybe this https://scicomp.stackexchange.com/a/13048/25606
It sounds like you want the ** multi-precision plot function from SymPy**, which is capable of plotting arbitrary black-box functions over a given range.