I will present a new method of generating high order quadrature formulas for domains with piecewise smooth curved boundaries and for surfaces, developed jointly with Bruno Degli Esposti. The main idea is to use the divergence theorem thatconnects two  integrals, one over the domain and another over its boundary, and require exactness of the quadrature for the difference of the two integrals applied to the divergence and the normal projection, respectively, of the elements of a finite dimensional approximation space. This way, in contrast to traditional methods, there is no need to precompute the  integrals of the basis functions or to partition the domain into smooth images of simple elements like cubes or simplices, which are challenging tasks for complicated domains.