  A set of fully numerical algorithms for evaluating the four-dimensional
singular integrals arising from Galerkin surface integral equation methods over
conforming quadrilateral meshes is presented. This work is an extension of
DIRECTFN, which was recently developed for the case of triangular patches,
utilizing in a same fashion a series of coordinate transformations together
with appropriate integration re-orderings. The resulting formulas consist of
sufficiently smooth kernels and exhibit several favorable characteristics when
compared with the vast majority of the methods currently available. More
specifically, they can be applied---without modifications---to the following
challenging cases: 1) weakly and strongly singular kernels, 2) basis and
testing functions of arbitrary order, 3) planar and curvilinear patches, 4)
problem-specific Green functions (e.g. expressed in spectral integral form), 5)
spectral convergence to machine precision. Finally, we show that the overall
performance of the fully numerical schemes can be further improved by a
judicious choice of the integration order for each dimension.
