Looking at my messy result, it is actually a sum of several integrals which Mathematica returned unevaluated.
Doing the same thing numerically (using NIntegrate) gives the result 19.9868 in 0.34 s. (Mathematica 10.0, Windows 8.1). Actually the result has a tiny imaginary part 5.5 *10^-24 as an artifact of the numerical method ('numerical noise').
Where do you get the z from?
1722 seconds on my laptop, for a very messy result. I wonder if we're solving the same integral.
Timing[ Integrate[(x^2 + y^2)*Boole[(x - z)^2 + (y - z)^2 <= 1], {x, -1, 3}, {y, -1, 3}, {z, 0, 2}] ]
took 11.5 seconds on my computer.
Me too, but it's an interesting 3D problem.
Wow! Very good. Thanks!