Q:

Calculation of the EFIE integral

I need help computing the following integral:
$$
\int_{}\frac{(1+jk|\vec{r}-\vec{r}^\prime|)e^{-jk|\vec{r}-\vec{r}^\prime|}}{|\vec{r}-\vec{r}^\prime|}d\vec{r}^\prime
$$
in this integral $\vec{r}$ and $\vec{r}^\prime$ are the position vectors that correspond to the observation point and source point, respectively. 
This integral comes from the Electric-Field Integral Equation (EFIE) kernel.

A:

A classic paper for evaluation of the integrals commonly present in computational electromagnetics (EM) is:

D. R. Wilton, S. M. Rao, A. W. Glisson, D. H. Schaubert, O. M. Al-Bundak, and C. M. Butler, "Potential integrals for uniform and linear source distributions of polygonal and polyhedral domains," IEEE Trans. Antennas Propag., vol. AP-32, no. 3, pp. 276-281, Mar. 1984.

The challenge to evaluate the EFIE integrals is the required singularity treatment. If both $\mathbf{r}$ and $\mathbf{r}^\prime$ belong to the same patch in the discretized geometry, one cannot use just quadrature-based integration.
The idea is to subtract a singularity and integrate it analytically and the rest (already without a singularity), numerically. So, say you need to calculate the integral ($R=|\mathbf{r}-\mathbf{r'}|$):
$$
\int\limits_S\frac{e^{-jkR}}{R}ds'=\underbrace{\int\limits_S\frac{e^{-jkR}-1}{R}ds'}_{I_1}+\underbrace{\int\limits_S\frac{1}{R}ds'}_{I_2}
$$
Here, $S$ is the source patch and $\mathbf{r}$ is fixed.
Now, integral $I_1$ does not have a singularity and can be evaluated using, say, 2-D Gaussian quadrature (use barycentric coordinates and apply Gaussian quadrature over $\xi$, $\eta$ with a proper Jacobian)
Integral $I_2$ can be found analytically. The formula for a planar polygon $S$ is given in the referenced paper.
You still might want to use this technique of singularity subtraction even when observation and source patches are different, but relatively close together. Such singularity treatment will decrease the error coming from numerical integration and reduce the number of quadrature points required to achieve the accuracy for nearby interactions.
I assume, that right now you've discretized your EFIE using pulse-basis functions and you are testing your EFIE using pulses as well. That might be good enough for certain scenarios, but I will point out that the use of proper basis functions and testing procedure might be important later in your development process.
Note: 
there are other techniques to tackle the singularity in this Green's function, including special quadrature rules, but I personally prefer this approach.

