  As part of our development of a computer code to perform 3D `constrained
evolution' of Einstein's equations in 3+1 form, we discuss issues regarding the
efficient solution of elliptic equations on domains containing holes (i.e.,
excised regions), via the multigrid method. We consider as a test case the
Poisson equation with a nonlinear term added, as a means of illustrating the
principles involved, and move to a "real world" 3-dimensional problem which is
the solution of the conformally flat Hamiltonian constraint with Dirichlet and
Robin boundary conditions. Using our vertex-centered multigrid code, we
demonstrate globally second-order-accurate solutions of elliptic equations over
domains containing holes, in two and three spatial dimensions. Keys to the
success of this method are the choice of the restriction operator near the
holes and definition of the location of the inner boundary. In some cases (e.g.
two holes in two dimensions), more and more smoothing may be required as the
mesh spacing decreases to zero; however for the resolutions currently of
interest to many numerical relativists, it is feasible to maintain second order
convergence by concentrating smoothing (spatially) where it is needed most.
This paper, and our publicly available source code, are intended to serve as
semi-pedagogical guides for those who may wish to implement similar schemes.
