Nothing is really smooth on the computer — it is just a question of whether you tolerate non-smoothness on the level of roundoff errors (which is unavoidable) or on the level of the truncation error. If the latter is unacceptable, then I agree that you need a non-adaptive rule.
If you just want to do a fixed-order/non-adaptive tensor-product rule, you can implement that yourself easily. e.g. for d-dimensional integration on [a,b]^d given the points and weights (x,w) of a 1d rule on [a,b], you can just do:
function cubature(f, x::AbstractVector, w::AbstractVector{<:Real}, ::Val{d}) where {d}
(ax = axes(x,1)) == axes(w,1) || throw(DimensionMismatch())
integral = zero(eltype(w))
for i in Iterators.product(ntuple(_ -> ax, Val{d}())...)
point = getindex.((x,), i)
weight = prod(getindex.((w,), i))
integral += f(point...) * weight
end
return integral
end
For example, to integrate e^{x_1 x_2 x_3 } in [0,1]^3:
julia> using QuadGK
julia> x, w = gauss(10, 0, 1); # 10-point Gauss rule on [0,1]
julia> cubature((x1,x2,x3) -> exp(x1*x2*x3), x, w, Val(3))
1.1464990725286444
which is correct to 13 significant digits.
This should be AD-able as-is (with respect to parameters of the integrand) using e.g. ForwardDiff, assuming that your integrand f is AD-able.