A new method is described for constructing geomagnetic flux coordinates, constrained by magnetospheric and ionospheric boundary surfaces. The technique is especially useful for computing boundary value solutions in a geometry defined by a realistic, three-dimensional geomagnetic field with boundary conditions specified on surfaces oriented arbitrarily relative to the background magnetic field. An adaptable algorithm for calculating the metric tensor for a general class of geomagnetic flux coordinate systems is presented. Application to ultralow-frequency wave propagation in the magnetosphere is illustrated by solving the equations of linear, one-fluid magnetohydrodynamics for the driven field line resonance in a dipolar geomagnetic field bounded by a spherical ionosphere. A novel diagnostic for energy flow in the driven resonance problem shows that MHD wave power flows radially inward as a compressional wave from the boundary driver to the resonant flux surface, is diverted azimuthally at the resonant surface, and becomes field-aligned in the resonant surface by coupling to the magnetically guided Alfv¿n wave at locations where the shear mode is in azimuthal phase quadrature with the compressional driver.