For Hankel transform, one can classify the methods into four major groups:
- Numerical quadrature-based.
- Fourier-based ones.
- Asymptotic expansion of Bessel into sines and cosines.
- Projection-slice methods.
The following paper gives a nice overview of this methods (types of methods):
According to this paper (p. 3, Sec. 4):
The main disadvantage of the trapezoidal rule is its poor computational efficiency. ... the trapezoidal method was found to be almost always as reliable as any other method tested. ...we found it useful as a benchmark to test more efficient algorithms.
So, two directions are possible:
- Either find a more efficient numerical quadrature rule.
- Follow along with the Hankel transform-like direction.
In fact, this integral is very similar to the $n$th order Hankel transform. However, I see as a major complication here that $f(r)$ is not only extremely expensive, but oscillatory, and, potentially, highly oscillatory. In the given formulation, $f(r)$ is almost a black-box with some known properties which are actually bad. A lot of the standard techniques for evaluation Hankel transforms use the fact that $f(r)$ is eventually monotonic.
Because of that, I would look into the direction of Filon-type quadrature, that is used for highly-oscillatory integrals. But I think you would have to use some knowledge about the behaviour of $f(r)$ to make it work. That is regarding the direction 1. The following reference (as well as the one for the next part) will be helpful for general understanding:
For direction 2, I would advise trying to develop some form of a projection-slice methods. I am personally not aware of a ready method of this type developed for your integral.
The following references might be useful (for the Filon quadrature, as well):
Another interesting approach worth trying (me stumbling upon this paper was the main motivation to write this answer) is to apply Mellin transform. However, it might not be worth it if it required too many evaluations of $f(r)$ in the first place: