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.