electromagnetic reciprocity. The cost of computation involved in naive Galerkin implementations is typically quite severe. One must loop over each pair of elements (so we get n2 interactions) and for each pair of elements we loop through Gauss points in the elements producing a multiplicative factor proportional to the number of Gauss-points squared. Also, the function evaluations required are typically quite expensive, involving trigonometric/hyperbolic function calls. Nonetheless, the principal source of the computational cost is this double-loop over elements producing a fully populated matrix. The Green's functions, or fundamental solutions, are often problematic to integrate as they are based on a solution of the system equations subject to a singularity load (e.g. the electrical field arising from a point charge). Integrating such singular fields is not easy. For simple element geometries (e.g. planar triangles) analytical integration can be used. For more general elements, it is possible to design purely numerical schemes that adapt to the singularity, but at great computational cost. Of course, when source point and target element (where the integration is done) are far-apart, the local gradient surrounding the point need not be quantified exactly and it becomes possible to integrate easily due to the smooth decay of the fundamental solution. It is this feature that is typically employed in schemes designed to accelerate boundary element problem calculations. Derivation of closed-form Green's functions is of particular interest in boundary element method, especially in electromagnetics. Specifically in the analysis of layered media, derivation of spatial-domain Green's function necessitates the inversion of analytically-derivable spectral-domain Green's function through Sommerfeld path integral. This integral can not be evaluated analytically and its numerical integration is costly due to its oscillatory and slowly-converging behaviour. For a robust analysis, spatial Green's functions are approximated as complex exponentials with methods such as Prony's method or generalized pencil of function, and the integral is evaluated with Sommerfeld identity. This method is known as discrete complex image method. Comparison to other methods The boundary element method is often more efficient than other methods, including finite elements, in terms of computational resources for problems where there is a small surface/volume ratio. Conceptually, it works by constructing a "mesh" over the modelled surface. However, for many problems boundary element methods are significantly less efficient than volume-discretisation methods (finite element method, finite difference method, finite volume method). A good example of application of the boundary element method is efficient calculation of