  We present a novel method to perform numerical integration over curved
polyhedra enclosed by high-order parametric surfaces. Such a polyhedron is
first decomposed into a set of triangular and/or rectangular pyramids, whose
certain faces correspond to the given parametric surfaces. Each pyramid serves
as an integration cell with a geometric mapping from a standard parent domain
(e.g., a unit cube), where the tensor-product Gauss quadrature is adopted. As
no constraint is imposed on the decomposition, certain resulting pyramids may
intersect with themselves, and thus their geometric mappings may present
negative Jacobian values. We call such cells the folded cells and refer to the
corresponding decomposition as a folded decomposition. We show that folded
cells do not cause any issues in practice as they are only used to numerically
compute certain integrals of interest. The same idea can be applied to planar
curved polygons as well. We demonstrate both theoretically and numerically that
folded cells can retain the same accuracy as the cells with strictly positive
Jacobians. On the other hand, folded cells allow for a much easier and much
more flexible decomposition for general curved polyhedra, on which one can
robustly compute integrals. In the end, we show that folded cells can flexibly
and robustly accommodate real-world complex geometries by presenting several
examples in the context of immersed isogeometric analysis, where involved sharp
features can be well respected in generating integration cells.
