Get 2D Gauss–Legendre quadrature points and weights using NumPy

Question

NumPy provides the np.polynomial.legendre.leggauss() function to compute the sample points and weights for Gauss-Legendre quadrature. I'm trying to use this function to get the 2D points and weights for a quadrilateral. I am able to use the function as shown here:

def gauss2D(n):
    x, w = np.polynomial.legendre.leggauss(n)

    weights = []
    gauss_pts = []

    for i in range(n):
        for j in range(n):
            wts = w[i] * w[j]
            weights.append(wts)

            g = [x[i], x[j]]
            gauss_pts.append(g)

    return np.array(weights), np.array(gauss_pts)

Which outputs the following for n=2 integration points:

weights
 [1. 1. 1. 1.]

points
 [[-0.57735027 -0.57735027]
 [-0.57735027  0.57735027]
 [ 0.57735027 -0.57735027]
 [ 0.57735027  0.57735027]]

If it's possible, I would like to get rid of the for-loops by using NumPy arrays but my attempt (see function below) does not generate the correct results.

def gauss2Dnumpy(n):
    x, w = np.polynomial.legendre.leggauss(n)

    weights = np.concatenate((w, w))

    gauss_pts = np.zeros((n*2, n))

    for i in range(n*2):
        for j in range(n):
            gauss_pts[i] = x[j], x[j]

    return weights, gauss_pts

Is there a way to accomplish this without the for-loops?

Answer 1

This would be:

x, w = np.polynomial.legendre.leggauss(n)

gauss_pts = np.array(np.meshgrid(x,x,indexing='ij')).reshape(2,-1).T
weights = (w*w[:,None]).ravel()