The inclusion of domain (point) sources into a three dimensional boundary element method while solving the Helmholtz equation is described. The method is fully desingularized which allows for the use of higher order quadratic elements on the surfaces of the problem with ease. The effect of the monopole sources ends up on the right hand side of the resulting matrix system. Several carefully chosen examples are shown, such as sources near and within a concentric spherical core-shell scatterer as a verification case, a curved focusing surface and a multi-scale acoustic lens.