Multivariate Orthogonal Polynomial Generation

I'm trying to apply the stochastic galerkin method to partial differential equation with multiple uniform random coefficients. I'm puzzled as to how to extend the corresponding orthogonal (legendre) polynomial basis into higher dimensions. Is there a systematic method for doing this?

Is there is a paper/book that discusses an algorithm to generate orthogonal polynomials (hermite, laguerre, jacobi,...) into higher dimensional spaces ($d\ge 2$) ?

Suppose that you know the orthogonal polynomial basis in a single dimension $(x)$ of each degree $i$ up to some desired order $K$. That is, we know $$p_0(x),p_1(x),...p_i(x),...,p_K(x)$$
To extend this into a two dimensions $(x,y)$, we need only consider the product between 1D polynomials in (x) and (y) and collect only the products whose total degree is less than or equal to $K$. That is, we collect the term $p_i(x)\cdot p_j(y)$ if $i+j\leq K$.
To extend this further into any arbitrary dimension $d$, simply created one dimensional orthogonal polynomial basis in each direction $1,2,...d$. Then, to obtain the multidimensional basis functions, collect mutual products of these bases $$p_{i_1}(x_1)\cdot p_{i_2}(x_2)\cdot ... \cdot p_{i_d}(x_d)$$
such that $i_1+i_2+...i_d\leq K$.