  We present a systematic computational framework for generating positive
quadrature rules in multiple dimensions on general geometries. A direct
moment-matching formulation that enforces exact integration on polynomial
subspaces yields nonlinear conditions and geometric constraints on nodes and
weights. We use penalty methods to address the geometric constraints, and
subsequently solve a quadratic minimization problem via the Gauss-Newton
method. Our analysis provides guidance on requisite sizes of quadrature rules
for a given polynomial subspace, and furnishes useful user-end stability bounds
on error in the quadrature rule in the case when the polynomial moment
conditions are violated by a small amount due to, e.g., finite precision
limitations or stagnation of the optimization procedure. We present several
numerical examples investigating optimal low-degree quadrature rules, Lebesgue
constants, and 100-dimensional quadrature. Our capstone examples compare our
quadrature approach to popular alternatives, such as sparse grids and
quasi-Monte Carlo methods, for problems in linear elasticity and topology
optimization.
