make padua nodes vary slowest in the first dimension

examples/plot-padua-nodes.py

@@ -40,7 +40,9 @@ def plot_padua_nodes(alpha, beta, order, family):
     ax.grid()
 
     ax.plot(curve[0], curve[1])
-    ax.plot(nodes[0], nodes[1], "ko")
+    for i, xi in enumerate(nodes.T):
+        ax.plot(nodes[0], nodes[1], "ko", markersize=8)
+        ax.text(*xi, str(i), color="k", fontsize=24, fontweight="bold")
 
     ax.set_xlim([-1, 1])
     ax.set_ylim([-1, 1])

modepy/nodes.py

@@ -370,50 +370,44 @@ def legendre_gauss_lobatto_tensor_product_nodes(dims, n):
 
 # {{{ Padua nodes
 
-def _first_padua_jacobi_nodes(alpha, beta, order):
+def _make_padua_grid_nodes(alpha, beta, order):
     from modepy.quadrature.jacobi_gauss import jacobi_gauss_lobatto_nodes
-    mu = jacobi_gauss_lobatto_nodes(alpha, beta, order)[::-1]
-    eta = jacobi_gauss_lobatto_nodes(alpha, beta, order + 1)[::-1]
+    mu = jacobi_gauss_lobatto_nodes(alpha, beta, order)
+    eta = jacobi_gauss_lobatto_nodes(alpha, beta, order + 1)
 
-    nodes = np.stack(np.meshgrid(mu, eta))
-    return np.hstack([
-        nodes[:, ::2, 1::2].reshape(2, -1),
-        nodes[:, 1::2, ::2].reshape(2, -1)
-        ])
+    return mu, eta
 
 
-def _second_padua_jacobi_nodes(alpha, beta, order):
-    rot = np.array([
-        [np.cos(np.pi/2), -np.sin(np.pi/2)],
-        [np.sin(np.pi/2), np.cos(np.pi/2)],
-        ])
-    if order % 2 == 0:
-        rot = -rot
+def _make_padua_jacobi_nodes(mu, eta, odd_or_even):
+    nodes = np.stack(np.meshgrid(mu, eta, indexing="ij"))
+    indices = np.sum(
+            np.meshgrid(np.arange(mu.size), np.arange(eta.size), indexing="ij"),
+            axis=0)
 
-    nodes = _first_padua_jacobi_nodes(alpha, beta, order)
-    return rot @ nodes
+    return nodes[:, indices % 2 == odd_or_even].reshape(2, -1)
 
 
-def _third_padua_jacobi_nodes(alpha, beta, order):
-    rot = np.array([
-        [np.cos(np.pi), -np.sin(np.pi)],
-        [np.sin(np.pi), np.cos(np.pi)],
-        ])
+def _first_padua_jacobi_nodes(alpha, beta, order):
+    mu, eta = _make_padua_grid_nodes(alpha, beta, order)
+    return _make_padua_jacobi_nodes(mu, eta, 0)
 
-    nodes = _first_padua_jacobi_nodes(alpha, beta, order)
-    return rot @ nodes
 
+def _second_padua_jacobi_nodes(alpha, beta, order):
+    # NOTE: these are just "rotated" by pi/2 from the first family
+    mu, eta = _make_padua_grid_nodes(alpha, beta, order)
+    return _make_padua_jacobi_nodes(eta, mu, 0)
 
-def _fourth_padua_jacobi_nodes(alpha, beta, order):
-    rot = np.array([
-        [np.cos(np.pi/2), -np.sin(np.pi/2)],
-        [np.sin(np.pi/2), np.cos(np.pi/2)],
-        ])
-    if order % 2 == 1:
-        rot = -rot
 
-    nodes = _first_padua_jacobi_nodes(alpha, beta, order)
-    return rot @ nodes
+def _third_padua_jacobi_nodes(alpha, beta, order):
+    # NOTE: these are just "rotated" by pi from the first family
+    mu, eta = _make_padua_grid_nodes(alpha, beta, order)
+    return _make_padua_jacobi_nodes(mu, eta, 1)
+
+
+def _fourth_padua_jacobi_nodes(alpha, beta, order):
+    # NOTE: these are just "rotated" by 2 pi/3 from the first family
+    mu, eta = _make_padua_grid_nodes(alpha, beta, order)
+    return _make_padua_jacobi_nodes(eta, mu, 1)
 
 
 def padua_jacobi_nodes(alpha, beta, order, family="first"):