Node elimination is a numerical approach for obtaining cubature rules for the approximation of multivariate integrals over domains in Rn. Beginning with a known cubature, nodes are selected for elimination, and a new, more efficient rule is constructed by iteratively solving the moment equations. In this work, a new node elimination criterion is introduced that is based on linearization of the moment equations. In addition, a penalized iterative solver is introduced that ensures positivity of weights and interiority of nodes. We aim to construct a universal algorithm for convex polytopes that produces efficient cubature rules without any user intervention or parameter tuning, which is reflected in the implementation of our package gen-quad. Strategies for constructing the initial rules for various polytopes in several space dimensions are described. Highly efficient rules in four and higher dimensions are presented. The new rules are compared to those that are obtained by combining transformed tensor products of one dimensional quadrature rules, as well as with known analytically and numerically constructed cubature rules
