In this paper, we consider a general decomposition of any convex polytope P \subset \BbbR n into a set of subpolytopes \Omega i and develop methods for approximating a definite integral of a given function f over P when, rather than its values at some points, a number of integrals of f over the faces of \Omega i are only available. We present several new families of extended integration formulas that contain such integrals and provide in a special case of our result the multivariate analogues of midpoint, trapezoidal, Hammer, and Simpson rules. The paper also presents the best possible explicit constants for their approximation errors. Here we succeed in finding the connection between minimization of the global error estimate and construction of centroidal Voronoi tessellations of a given polytope with special density function depending on properties of the integrand. In the case of integrands with strong singularities, it leads to essential reduction of the error. These ideas were extended to a more general case, in which the domain is not necessary polytope and is not necessary convex. We perform numerical tests with integrands having steep gradients which allow the comparison of the new cubature formulas and show their accuracy and rates of convergence.