You aren’t setting any tolerance in these integrals, so it is using the default relative tolerance of \approx 10^{-8}. This is a problem here because the integral fooint is actually zero, so a relative tolerance |err| ≤ rtol * |int| can never be satisfied (unless the error is exactly zero), so it will keep refining the integrals until it hits the limits of machine precision (which will take a long time in 2d).
Passing an absolute tolerance atol=1e-8 to both quadgk calls causes the calculation to run very quickly.
(You might also look into IteratedIntegration.jl for performing nested quadgk calls more efficiently. Or HCubature.jl via a change of variables to map to a rectangular domain. Or, if you are really just integrating low-degree polynomials and similar very smooth functions, use a tensor-product of fixed-order Gaussian quadrature rules, or SimplexQuad.jl, or similar.)