See also the recent thread: Numerical integration over 3D domain of vector-valued function
You can use HCubature.jl, Cubature.jl, nested quadgk calls (made more efficient by IteratedIntegration.jl), MeshIntegrals.jl, Cuba.jl, or explicit tensor products of fixed-order 1d quadrature rules (which has the tradeoff off not being adaptive, but for smooth integrands can be very efficient), for example.
As discussed in e.g. the linked thread, if you have a singularity in you integrand, then if at all possible you want to do a transformation that removes the singularity (or choose a quadrature rule that builds that singularity in).