I try to reformulate the commentaries of fgp as an answer.
First consider a probability space $(\Omega, \mathcal{A}, \mathbb{P}')$ and a random variable uniformly distributed on $[0,2\pi]$:
$$X: (\Omega, \mathcal{A}, \mathbb{P}') \rightarrow ([0,2\pi],\mathcal{B}([0,2\pi])).$$
Furthermore consider the parametrization of the unit circle
$$p: [0,2\pi] \rightarrow S^1; \quad x \mapsto e^{i\cdot x}$$
which is continuous.
Now consider the space $(S^1,\mathcal{B}(S^1))$ where $\mathcal{B}(S^1)$ is the $\sigma$-algebra generated by the open sets of $S^1$.
We define $\mathbb{P}$ to be the probability measure induced by the map
$$p\circ X: (\Omega, \mathcal{A}, \mathbb{P}') \rightarrow (S^1,\mathcal{B}(S^1)); \quad \omega \mapsto e^{i X(\omega)}.$$
Then we have
\begin{equation}
P(Z\in A) = \frac{1}{2\pi i}\int_{A}{\frac{1}{z}}dz \qquad \forall A \in \mathcal{B}(S^1)
\end{equation}
or more generally
$$E[f(Z)]=\frac{1}{2\pi i}\int_{S^1}{\frac{f(z)}{z}}dz$$
for any measurable, bounded function $f: S^1 \rightarrow \mathbb{R}$
So we must be careful what the measurable space really is. In the case of the "density" in the question we are considering the probability space $(S^1,\mathcal{B}(S^1), \mathbb{P})$ and on this we get can give the value of $\mathbb{P}$ via the above formula.
Hence the above is not a density with respect to the Lebesgue-measure on $\mathbb{C}$, but a way of formulating the value with the help of the parametrization of (or - to be more precise - a contour integral along) the unit circle $S^1$.
We are computing the probability that two random points in the unit disk (here random means with respect to a uniform probability distribution) lie within a distance $\leq\frac{1}{2}$. We may assume without loss of generality that the first point lies on the segment joining the origin with the point $(1,0)$, with a probability density function supported on $[0,1]$ and given by $f(z)=2z$. If we assume to know the position of the first point, the second point lies within a distance $\leq\frac{1}{2}$ iff it belongs to the original disk and to a smaller disk. If $z\in[0,1]$ and $g(z)$ is the Lebesgue measure of the set:
$$ E(z) = \left\{(x,y): x^2+y^2\leq 1, (x-z)^2+y^2\leq\frac{1}{4}\right\}, $$
the wanted probability is just:
$$ P=\frac{1}{\pi}\int_{0}^{1} f(z)g(z)\,dz = \frac{2}{\pi}\int_{0}^{1} z\cdot g(z)\,dz.$$
So the problem boils down to computing $g(z)$, that has a rather ugly closed form in terms of the $\arcsin$ function. Numerically,
$$ P\approx 19.728\%.$$
Best Answer
Answer: On the second line of your derivation.
Let's do it correctly: $$ \begin{eqnarray} F_U\left(u\right) &=& \Pr\left(\sqrt{X^2+Y^2} \leqslant u \right) \\ &=& \int_S f_{X,Y}\left(x,y\right) I\left(\sqrt{X^2+Y^2} \leqslant u\right) \mathrm{d}x\mathrm{d}y \\ &=& \left.\int_{-1}^{1} \left(\int_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}} \frac{1}{\pi} I\left(\sqrt{X^2+Y^2} \leqslant u\right) \mathrm{d} y \right)\mathrm{d} x \right|_{\text{assuming } 0<u \leqslant u} \\ &=& \int_{-u}^{u} \left(\int_{-\sqrt{u^2-x^2}}^{\sqrt{u^2-x^2}} \frac{1}{\pi} \mathrm{d} y \right)\mathrm{d} x \\ &=& \left.\int_{-u}^{u} \frac{2}{\pi} \sqrt{u^2-x^2} \mathrm{d}x \right|_{x = z \cdot u} \\ &=& u^2 \int_{-1}^{1} \frac{2}{\pi} \sqrt{1-z^2} \mathrm{d}z = u^2 \end{eqnarray} $$
As a side note, it would have been easier to work in polar coordinates here.