Engineering applications frequently require the numerical solution of elliptic boundary value problems in irregularly shaped domains. For smooth problems, spectral element methods have proved very successful, since they can accommodate fairly complicated geometries while retaining a rapid rate of convergence. Geometric singularities, however, often give rise to singular solutions. The accuracy of the spectral element methods is then degraded, and they offer no apparent advantage over low-order finite element methods. In many cases, however, the singular structure of the solution is known, and its form may be exploited by the spectral element method. Among the Various ways of doing so (through supplementary basis functions, eigenfunction expansions, and graded meshes), the method of auxiliary mapping proves to be particularly effective. For certain simple cases, the problem is transformed to a coordinate system in which the solution is analytic, and exponential convergence is recovered. Even when this is not possible, the singularity is usually much weaker after mapping, so that other treatments are more effective in the new coordinate system. In this paper, we study different ways of treating singularities, and in particular, the method of auxiliary mapping coupled with the use of supplementary basis functions. Error estimates are presented explaining why the combined approach is more effective, and these estimates are confirmed through a number of numerical experiments for the Laplace, Poisson, and Helmholtz equations, (C) 1995 Academic Press, Inc.