Can't compute the Integral

Use NIntegrate and pattern test the arguments to f1,f2,f3 with NumericQ to ensure NIntegrate does not do any symbolic preprocessing first. You do not need N as it has no effect here. I've given $mu$ and $i$ values of 1, but feel free to change them:

f1[x_?NumericQ, y_?NumericQ, z_?NumericQ, r_?NumericQ] := 
 NIntegrate[
  Cos[a]/(x^2 + y^2 + z^2 + r^2 - 2*r*(x*Cos[a] + y*Sin[a]))^(3/
      2), {a, 0, 2*Pi}];

f2[x_?NumericQ, y_?NumericQ, z_?NumericQ, r_?NumericQ] := 
  NIntegrate[
   Sin[a]/(x^2 + y^2 + z^2 + r^2 - 2*r*(x*Cos[a] + y*Sin[a]))^(3/
       2), {a, 0, 2*Pi}];

f3[x_?NumericQ, y_?NumericQ, z_?NumericQ, r_?NumericQ] := 
  NIntegrate[(r - y*Sin[a] - 
      x*Cos[a])/(x^2 + y^2 + z^2 + r^2 - 
       2*r*(x*Cos[a] + y*Sin[a]))^(3/2), {a, 0, 2*Pi}];

With[{i = 1, μ = 1, r = .1},
 VectorPlot3D[{
   (μ*i*r*z/(4*Pi))*f1[x, y, z, r],
   (μ*i*r*z/(4*Pi))*f2[x, y, z, r],
   (μ*i*r/(4*Pi))*f3[x, y, z, r]
   }, {x, -0.5, 0.5}, {y, -0.5, 0.5}, {z, -0.5, 0.5}]
 ]

With Version 12.1.1, at least, the integrals can be obtained symbolically, and I obtained the first one in less than an hour. Unfortunately, it is an enormous conditional expression. But, there is an easier way. Beginning with the expression,

fint = (x^2 + y^2 + z^2 + r^2 - 2*r*(x*Cos[a] + y*Sin[a]))^(-1/2);

compute

f1int = Simplify[(D[fint, x] - x D[fint, z]/z)/r]
(* Cos[a]/(r^2 + x^2 + y^2 + z^2 - 2 r x Cos[a] - 2 r y Sin[a])^(3/2) *)

f2int = Simplify[(D[fint, y] - y D[fint, z]/z)/r]
(* Sin[a]/(r^2 + x^2 + y^2 + z^2 - 2 r x Cos[a] - 2 r y Sin[a])^(3/2) *)

f3int = Simplify[-D[fint, r]]
(* (r - x Cos[a] - y Sin[a])/(r^2 + x^2 + y^2 + z^2 - 2 r x Cos[a] - 2 r y Sin[a])^(3/2) *)

The terms f1int, f2int, and f3int are recognizable as the integrands of f1, f2, and f3 in the question. Therefore, the three integrals in the question are simply the corresponding combinations of derivatives of the fint integral, which can be obtained in a few minutes. (The second assumption avoids singularities in f.)

f = Integrate[f, {a, 0, 2*Pi}, Assumptions -> (x | y | z | r) ∈ Reals && 
    x^2 + y^2 + z^2 + r^2 + 2 r Sqrt[x^2 + y^2] > 0]
(* (4 EllipticK[(4 r Sqrt[x^2 + y^2])/(r^2 + x^2 + y^2 + 2 r Sqrt[x^2 + y^2] + z^2)])/
   Sqrt[r^2 + x^2 + y^2 + 2 r Sqrt[x^2 + y^2] + z^2] *)

For instance, the f3 integral, obtained in seconds, is

FullSimplify[-D[f, r]]
(* (2 (-(-r^2 + x^2 + y^2 + z^2) (r^2 + x^2 + y^2 + 2 r Sqrt[x^2 + y^2] + z^2) 
   EllipticE[(4 r Sqrt[x^2 + y^2])/(r^2 + x^2 + y^2 + 2 r Sqrt[x^2 + y^2] + z^2)] + 
   ((-r^2 + x^2 + y^2)^2 + 2 (r^2 + x^2 + y^2) z^2 + z^4) EllipticK[(4 r Sqrt[x^2 + y^2])
   /(r^2 + x^2 + y^2 + 2 r Sqrt[x^2 + y^2] + z^2)]))/(r (r^2 + x^2 + y^2 - 2 r 
   Sqrt[x^2 + y^2] + z^2) (r^2 + x^2 + y^2 + 2 r Sqrt[x^2 + y^2] + z^2)^(3/2)) *)

The f1 and f2 integrals are computed similarly from

FullSimplify[(D[f, x] - x D[f, z]/z)/r]

FullSimplify[(D[f, y] - y D[f, z]/z)/r]

also in seconds.