  A non-singular formulation of the boundary integral method (BIM) is presented
for the Laplace equation whereby the well-known singularities that arise from
the fundamental solution are eliminated analytically. A key advantage of this
approach is that numerical errors that arise due to the proximity of nodes
located on osculating boundaries are suppressed. This is particularly relevant
in multi-scale problems where high accuracy is required without undue increase
in computational cost when the spacing between boundaries become much smaller
than their characteristic dimensions. The elimination of the singularities
means that standard quadrature can be used to evaluate the surface integrals
and this results in about 60% savings in coding effort. The new formulation
also affords a numerically robust way to calculate the potential close to the
boundaries. Detailed implementations of this approach are illustrated with
problems involving osculating boundaries, 2D domains with corners and a wave
drag problem in a 3D semi-infinite domain. The explicit formulation of problems
with axial symmetry is also given.
