Thank you fedja, that's exactly what I was looking for. Now I will try to generalize this idea for greater dimensions.

$Q_1 = [0,\frac{1}{2}]\times[0,\frac{1}{2}], Q_2 = [\frac{1}{2},\frac{3}{4}]\times[0,\frac{1}{4}], Q_3=[0,\frac{1}{4}]\times[\frac{1}{2},\frac{3}{4}], Q_4=[\frac{3}{4},\frac{7}{8}]\times[0,\frac{1}{8}]...$ and so on.

Maybe I did not explain well the problem, I'm not looking for an algorithm like the one given by the Riemann sum because in that case when you have finish with the n^2 cuboids (following your program) and you want a better approximation then you have to start again (with a slightly different version of the algorithm) and construct for example (n+1)^2 new cuboids (which are obviously not disjoint of the ones constructed before). I'm looking more for a generalization of this construction in $\mathbb{R}^2$: To fill the triangle with vertices (0,0), (0,1) and (1,0) you can put:

Gerhard, actually I was thinking more in a "real world" programmable algorithm, but thank you anyway. And fedja, how did you arrive to this conclusion? What constructions do you have in mind?