In three dimensional boundary element analysis, computation of integrals is an important aspect since it governs the accuracy of the analysis and also because it usually takes the major part of the CPU time. The integrals which determine the influence matrices, the internal field and its gradients contain (nearly) singular kernels of order lIr a (0:= 1,2,3,4,.**) where r is the distance between the source point and the integration point on the boundary element. For planar elements, analytical integration may be possible 1,2,6. However, it is becoming increasingly important in practical boundary element codes to use curved elements, such as the isoparametric elements, to model general curved surfaces. Since analytical integration is not possible for general isoparametric curved elements, one has to rely on numerical integration. When the distance d between the source point and the element over which the integration is performed is sufficiently large compared to the element size (d> 1), the standard Gauss-Legendre quadrature formula 1,3 works efficiently. However, when the source is actually on the element (d=O), the kernel 1I~ becomes singular and the straight forward application of the Gauss-Legendre quadrature formula breaks down. These integrals will be called singular integrals. Singular integrals occur when calculating the diagonals of the influence matrices.