@acer 

Maybe Intel oneAPI with its new OpenMP offload to GPU? I don't know yet how it behaves but it looks promising, and I believe it was the point of Intel switching to these new compilers.

@vv It's doing exactly what I was looking for. I'll give it a try. Thanks!

After another attempt, it works with:

int(piecewise(x^2+y^2<=1,1,0),[x=-1..1,y=-1..1]);

Initially I didn't even think Maple would yield a different result than with Heaviside, but apparently it does.

I'm still interested to know if there are better ways.

(answering as my comment to acer's answer was removed as well as his answer, and it might be useful to the reader passing by)

I think I got it. The problem is not so much n, as I don't see how it would be evalf'ed before entering u. The problem is 1.

If I understand correctly, in u(n)$n=1..10, if there are quotes around u(n), Maple will build u(1),u(2)... as expected, and evaluate with an integer n.

Without quotes, it's evalf'ed first, but since n is not yet known I thought nothing happens. Wrong! The expression is evalf'ed as much as possible, which means n remains a symbol, and 1 is replaced by an sfloat, and when n is replaced, there must be some type promotion from integer to sfloat, and everything is evaluated as if n was passed as an sfloat.

Good to know.

@acer 

I think I got it. The problem is not so much n, as I don't see how it would be evalf'ed before entering u. The problem is 1.

If I understand correctly, in u(n)$n=1..10, if there are quotes around u(n), Maple will build u(1),u(2)... as expected, and evaluate with an integer n.

Without quotes, it's evalf'ed first, but since n is not yet known I thought nothing happens. Wrong! The expression is evalf'ed as much as possible, which means n remains a symbol, and 1 is replaced by an sfloat, and when n is replaced, there must be some type promotion from integer to sfloat, and everything is evaluated as if n was passed as an sfloat.

Good to know.