I need to integrate a function in a region around the origin, the problem is that it has a singularity there. Having tried Integrate[] without any sort of response after waiting for a long time, I tried to perform the integration numerically. The code is something like this:

{Lx,Ly}={1/50, 1/50};
Aelem[x_, y_, n_, m_] := Aelem[x, y, n, m] = 
NIntegrate[current[x1, y1, n, m]/Sqrt[(x - x1)^2 + (y - y1)^2], {y1, -Ly, Ly}, {x1, -Lx, Lx}, Exclusions -> {0, 0}, Method -> {"LocalAdaptive"}]

Helem[i_, j_, n_, m_] := Helem[i, j, n, m] =
NIntegrate[current[x, y, i, j]*Aelem[x, y, n, m], {y, -Ly, Ly}, {x, -Lx, Lx}, Method -> {"LocalAdaptive"}]

And then these current functions have the shape

g[x_,y_,n_,m_]:=(2*Lx/(n*Pi))*Sin[n*Pi*(x + Lx)/(2*Lx)]*Sin[m*Pi*(y + Ly)/(2*Ly)]
current[x_,y_,n_,m_]:=D[g[x, y, n, m], y]
Helem[1,1,1,1]//AbsoluteTiming

I'm trying out some different integration strategies, and different ways to define the current function, to try to get the best solution. Running Helem[1,1,1,1]//AbsoluteTiming has some problems: 1. Sometimes the evaluation just takes for ever 2. I get the " NIntegrate::inumr error message: "The integrand has evaluated to non-numerical values (...)"" message, so I use ?NumericQ, but then I have problems with defining the current as the derivative of g; 3. Sometimes I get a really weird error message:

Throw::nocatch: Uncaught Throw[Holonomic`DifferentialRootReduceDump`y[NIntegrate`LevinRuleDump`x]+(Holonomic`DifferentialRootReduceDump`y^\[Prime]\[Prime])[NIntegrate`LevinRuleDump`x],NIntegrate`LevinRuleDump`FastLookupHolonomicDifferentialEquation] returned to top level.

I really need to be able to perform the integral as fast as possible, I need the function H for many combinations of (i,j,n,m), and I feel like I am shooting in the dark about the integration strategies and all. Does anyone have any input on how to improve this and make the calculations faster? Any input is appreciated!

Edit (25/06/2018): I am now trying to integrate all at once, so something like:

{Lx,Ly}={1/50, 1/50};
gnm[x_, y_, n_, m_] := (2*Lx/(n*Pi))*Sin[n*Pi*(x + Lx)/(2*Lx)]*Sin[m*Pi*(y + Ly)/(2*Ly)]
current[x_, y_, n_, m_] = D[gnm[x, y, n, m], y]; 
atonce[i_, j_, n_, m_] := atonce[i, j, n, m] =
NIntegrate[current[x, y, i, j]*current[x1, y1, n, m]/Sqrt[(x - x1)^2 + (y - y1)^2]
, {y, -Lx, Lx}, {x, -Lx, Lx}, {y1, -Lx, Lx}, {x1, -Lx, Lx}]
atonce[1,2,1,2]//AbsoluteTiming

And this is a little insane, whenever I try a different integration strategy (GlobalAdaptive, DuffyCoordinates, DoubleExponential), increase the MaxRecursion or the WorkingPrecision, the result becomes very different. I've tried to put both these parameters as high as 300, but I keep getting the error message "NIntegrate failed to converge to prescribed accuracy after X recursive bissections", or equivalent. I am considering trying something like integrating first on dx1dy1, somehow defining "x" and "y" as singularities, and only then integrate on dxdy, but I am unsure how to do this.

Edit 26/06/2018: Corrected the gnm function