We present an efficient method to solve the Poisson equation in embedded boundary (EB) domains. The original problem is divided into an inhomogeneous problem without the effects of EB and a homogeneous problem that imposes the effects of EB. The inhomogeneous problem is efficiently solved through a geometric multi-grid (GMG) solver and the homogenous problem is solved through a boundary element method (BEM) utilizing the free space Green’s function. Our method is robust and can handle sharp geometric features without any special treatment. Analytical expressions are presented for the boundary and the domain integrals in BEM to reduce the computational cost and integration error relative to numerical quadratures. Furthermore, a fast multipole method (FMM) is employed to evaluate the boundary integrals in BEM and reduce the computational complexity of BEM. Our method inherits the complementary advantages of both GMG and FMM and presents an efficient alternative with linear computational complexity even for problems involving complex geometries. We observe that the overall computational cost is an order of magnitude lower compared with a stand-alone FMM and is similar to that of an ideal GMG solver.