This paper introduces and studies a new class of multidimensional numerical integration, which we call " strongly positive definite cubature formulas ". We establish, among others, a characterization theorem providing necessary and sufficient conditions for the approximation error (based on such cubature formulas) to be bounded by the approximation error of the quadratic function. This result is derived as a consequence of two characterization results, which are of independent interest, for linear functionals obtained in a more general seeting. Thus, this paper extends some result previously reported in [2, 3] when convexity in the classical sense is only assumed. We also show that the centroidal Voronoi Tesselations provide an efficient way for constructing a class of optimal of cubature formulas. Numerical results for the two-dimensional test functions are given to illustrate the efficiency of our resulting cubature formulas