Summary: The solution of elliptic partial differential equations on composite overlapping grids using multigrid is discussed. An approach is described that provides a fast and memory efficient scheme for the solution of boundary value problems in complex geometries. The key aspects of the new scheme are an automatic coarse grid generation algorithm, an adaptive smoothing technique for adjusting residuals on different component grids, and the use of local smoothing near interpolation boundaries. Other important features include optimizations for Cartesian component grids, the use of over-relaxed red-black smoothers, and the generation of coarse grid operators through Galerkin averaging. Numerical results in two and three dimensions show that very good multigrid convergence rates can be obtained for Dirichlet, Neumann, and mixed boundary conditions. A comparison to Krylov-based solvers shows that the multigrid solver can be much faster and require significantly less memory.