A variety of new and old ideas have been combined to develop a variational algorithm for generating an integration mesh for use in electronic-structure calculations. While our interests are in deriving an efficient yet accurate mesh for local-orbital calculations, the method which we introduce is independent of basis set, very efficient, and allows systematic refinements. The method incorporates a new algorithm for integrating over arbitrary interstitial regions and an improved atomic-spheres integration technique (errors decrease exponentially with the number of radial points). The mesh points and volume elements are optimized to ensure that a large class of integrals are evaluated accurately. We employ two new ideas in optimizing the mesh. First, we introduce a set of one-dimensional Gaussian-quadrature schemes depending continuously on a variational parameter. The introduction of this variational parameter leads to substantially more efficient integration meshes compared to traditional Gaussian-quadrature schemes. Second, we show how force algorithms may be used to refine the resulting meshes. The accuracy and efficiency of our scheme is assessed by presenting charge-density integrations on a variety of large isolated clusters. The number of mesh points required using this scheme is as much as an order of magnitude less than what is required with conventional special-points algorithms. For example, for the tetrahedral ground state of a seventeen-atom C5H12 crystal fragment, a total of 6081 points is needed to integrate the total charge and total energy to six-decimal-place accuracy.