inverseRadialSmoothedCap::usage = "It is not possible to do a Fourier series periodic extension of the inverse radial function V(r) = -β/r, because the integral blows up at the origin.   It blows up away from the origin too, but for a periodic extension where the function is confined to a finite range around the origin that doesn't matter.

Here we introduce a modification to V(r) = -β/r, using a cubic polynomial in r to cap it.  We fit to match conditions:

V(ϵ) = -β/ϵ
V'(ϵ) = β/(ϵ^2)
V(0) = V(ϵ) - α/ϵ = -(α + β)/ϵ
V'(0) = 0." ;

Clear[epsilon, alpha, beta, vPlusCubic]
vPlusCubic[r_, epsilon_, alpha_, beta_] =
Module[{vv, ww, s, v, aa, bb},
vv[tr_, te_, talpha_, tbeta_] = aa ( tr - te)^3/te^4 + bb(tr - te)^2/ te^3 + tbeta(tr - te)/te^2 - tbeta/te ;

ww[tr_] = D[vv[tr, epsilon, alpha, beta], tr] ;

s= Solve[ {-(beta + alpha)/epsilon == vv[0, epsilon, alpha, beta], 0 == ww[0]}, {aa, bb}]  // Flatten ;

vv[r, epsilon, alpha, beta] /. s 
] ;


Clear[ vr ]
vr = OverVector[r] // TraditionalForm ;
inverseRadialSmoothedCap[r_,epsilon_,alpha_,beta_]:=
If [ Norm[r]<epsilon, vPlusCubic[Norm[r],epsilon,alpha,beta], -beta/Norm[r] ] ;

Column[ {
Row[{"vPlusCubic(r) = v(|",vr ,"|<ϵ) = ", vPlusCubic[r, ϵ, α, β] // TraditionalForm}],
Row[{"v(", vr," = (x, y)) = ", inverseRadialSmoothedCap[{x, y}, ϵ, α, β] // TraditionalForm}],
Row[{"v(0) = ",
Simplify[inverseRadialSmoothedCap[{0, 0}, ϵ, α, β] , ϵ > 0] // TraditionalForm}]
}]

Now, let’s look at a specific instance of such a function, setting α = β = 1, and ϵ = 0.5.  In two dimensions this function looks like

ClearAll[alpha, beta, epsilon, v]
alpha = 1 ;
beta = 1 ;
epsilon = 1 ;

v[{x_, y_}] = inverseRadialSmoothedCap[ {x, y}, epsilon, alpha, beta ] ;
(*v[{x,y}] // TraditionalForm*)
(*v[{0, 0}]*)

Timing[
With[{range = 5},
Plot3D[ 
v[{x, y}],
{x, -range, range}, 
{y, -range, range}, 
PlotRange -> {Full, Full, {v[{0, 0}]- 0.5, 0}},

(* http://mathematica.stackexchange.com/questions/13657/plot-showing-discontinuity-where-it-shouldnt *)
Exclusions->None
]
]
]

Clear[ft, a, b, am, gm, g1, g2, uOfx, vOfx]
a = { 5, 13 } ;
b = { 15, 3 } ;
(*a = { 20, 0 } ;
b = { 0, 20 } ;*)
am = { a, b } ;
gm =2 Pi Inverse[ am ] ;
g1 = gm[[All, 1]] ;
g2 = gm[[All, 2]] ;

uOfx[r_] = r . g1/(2 Pi) ;
vOfx[r_] = r . g2 /(2 Pi) ;
ft[{r_, s_}] = FourierSeries[v[a r + b s],{r,s},{3, 3},FourierParameters->{1, 2 Pi}]  ;
ft[{r, s}] // TraditionalForm

(* no good ... Mathematica didn't evaluate the Fourier integral *)

ft[{r, s}] /. { r -> uOfx[{x,y}], s -> vOfx[{x,y}] }