In this paper, we study the efficient numerical integration of functions with sharp
gradients and cusps. An adaptive integration algorithm is presented that systematically
improves the accuracy of the integration of a set of functions. The algorithm is based on a
divide and conquer strategy and is independent of the location of the sharp gradient or
cusp. The error analysis reveals that for a $C^0$ function (derivative-discontinuity at a
point), a rate of convergence of $n+1$ is obtained in $R^n$. Two applications of the
adaptive integration scheme are studied. First, we use the adaptive quadratures for the
integration of the regularized Heaviside function---a strongly localized function that is
used for modeling sharp gradients. Then, the adaptive quadratures are employed in the
enriched finite element solution of the all-electron Coulomb problem in crystalline
diamond. The source term and enrichment functions of this problem have sharp gradients and
cusps at the nuclei. We show that the optimal rate of convergence is obtained with only a
marginal increase in the number of integration points with respect to the pure finite
element solution with the same number of elements. The adaptive integration scheme is
simple, robust, and directly applicable to any generalized finite element method employing
enrichments with sharp local variations or cusps in $n$-dimensional parallelepiped
elements.