Actually, it defaults to the Euclidean norm (vecnorm) for arrays, and you can pass a different norm if desired via the norm keyword. It only falls back to elementwise isapprox if the norm is Inf or NaN. This is all documented. You don’t have to guess at the behavior.
We shouldn’t warn people against using this for arrays — if anything, we should warn people against implementing their own approximate equality for arrays, because experience suggests they will probably get it wrong. (e.g. NumPy got it wrong: their isclose function is only elementwise, has a nonzero default atol, and even has the property that isclose(x,y) is inequivalent to isclose(y,x)!)
I don’t think it should throw a DomainError for x ≈ 0. This is perfectly well defined, it just means x == 0 since the definition is norm(x-y) ≤ rtol*max(norm(x),norm(y)) for atol==0. It would be very weird for 0 ≈ 0 to throw an error rather than returning the correct answer (true), and this would be an especially unfortunate behavior for elementwise isapprox comparisons.