calculus and analysis - Can Mathematica help me evaluate an integral over disjoint disks $I = int_{D_1} int_{D_2} log|x-y| dy dx$?

I want to evaluate an integral that involves two disjoint unit disks $D_1$ and $D_2$. $D_1$ is centered at $(-2,0)$ and $D_2$ is centered at $(0,2)$. The integral I want to compute is

$$I = \int_{D_1} \int_{D_2} \log|x-y| dy dx.$$

I looked at the in-built Python integration methods and also the quadpy library but although they have lots of options for integration over a single disk, I couldn't find anything that can help me with integrating over two disjoint disks.

Is it possible to evaluate this integral in Mathematica? I don't need an optimum method, I just need to obtain the value of this integral.

Answer

NIntegrate[
 Log[Norm[{x1, x2} - {y1, y2}]], 
 {x1, x2} ∈ Disk[{-2, 0}, 1], 
 {y1, y2} ∈ Disk[{2, 0}, 1]
 ]

13.6822

(the actual result being 13.682176919165677),

In order to enter ∈, just type Esc e l Esc.

In order to increase the precision, use the option PrecisionGoal.

Edit

Another possibility that relieves Mathematica from the need to discretize the disks and that allows her to use higher-order quadrature formulas is to employ polar coordinates on each of the disks:

NIntegrate[
 Log[Norm[{r1 Cos[θ1] - 2,r1 Sin[θ1]} - {r2 Cos[θ2] + 2,r2 Sin[θ2]}]] r1 r2,

 {r1, 0, 1}, {θ1, -Pi, Pi},
 {r2, 0, 1}, {θ2, -Pi, Pi},
 Method -> "LocalAdaptive"
 ]