Numerical integration over 3D domain of vector-valued function

For adaptive integration, your allocations won’t be zero because of the data structures needed to keep track of the dynamic subdivision of the domain, though pre-allocation can help.

Your integrand has a singularity at the corner — you really want to do a transformations to eliminate this singularity to make the integration more efficient. Algorithms to do this are widely known in the integral-equation community (indeed, your integral looks much like a panel integral for a surface integral equation).

You might also want to lower the tolerance — by default it integrates to at least 8 digits, which might be more than you need.