add padua nodes

Author: alexfikl

examples/plot-padua-nodes.py

@@ -0,0 +1,65 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+import modepy as mp
+
+
+def _first_padua_curve(n, t):
+    return np.vstack([-np.cos((n + 1) * t), -np.cos(n * t)])
+
+
+def _second_padua_curve(n, t):
+    return np.vstack([-np.cos(n * t), -np.cos((n + 1) * t)])
+
+
+def _third_padua_curve(n, t):
+    return np.vstack([np.cos((n + 1) * t), np.cos(n * t)])
+
+
+def _fourth_padua_curve(n, t):
+    return np.vstack([np.cos(n * t), np.cos((n + 1) * t)])
+
+
+def plot_padua_nodes(alpha, beta, order, family):
+    if family == "first":
+        curve_fn = _first_padua_curve
+    elif family == "second":
+        curve_fn = _second_padua_curve
+    elif family == "third":
+        curve_fn = _third_padua_curve
+    elif family == "fourth":
+        curve_fn = _fourth_padua_curve
+
+    t = np.linspace(0, np.pi, 512)
+    curve = curve_fn(order, t)
+    nodes = mp.padua_jacobi_nodes(alpha, beta, order, family)
+    assert nodes.shape[1] == ((order + 1) * (order + 2) // 2)
+
+    fig = plt.figure()
+    ax = fig.gca()
+    ax.grid()
+
+    ax.plot(curve[0], curve[1])
+    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])
+    ax.set_aspect("equal")
+    fig.savefig(f"padua_nodes_order_{order:02d}_family_{family}")
+    plt.show()
+
+
+if __name__ == "__main__":
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--order", type=int, default=5)
+    parser.add_argument("--family", default="first",
+            choices=["first", "second", "third", "fourth"])
+    parser.add_argument("--alpha", type=float, default=-0.5)
+    parser.add_argument("--beta", type=float, default=-0.5)
+    args = parser.parse_args()
+
+    plot_padua_nodes(args.alpha, args.beta, args.order, args.family)

modepy/__init__.py

@@ -21,6 +21,7 @@
 THE SOFTWARE.
 """
 
+from pytools import MovedFunctionDeprecationWrapper
 
 from modepy.matrices import (
     diff_matrix_permutation, differentiation_matrices, inverse_mass_matrix,
@@ -37,7 +38,8 @@
 from modepy.nodes import (
     edge_clustered_nodes_for_space, equidistant_nodes, equispaced_nodes_for_space,
     legendre_gauss_lobatto_tensor_product_nodes, node_tuples_for_space,
-    random_nodes_for_shape, tensor_product_nodes, warp_and_blend_nodes)
+    padua_jacobi_nodes, padua_nodes, random_nodes_for_shape, tensor_product_nodes,
+    warp_and_blend_nodes)
 from modepy.quadrature import (
     LegendreGaussTensorProductQuadrature, Quadrature, QuadratureRuleUnavailable,
     TensorProductQuadrature, ZeroDimensionalQuadrature, quadrature_for_space)
@@ -82,6 +84,7 @@
 
         "equidistant_nodes", "warp_and_blend_nodes",
         "tensor_product_nodes", "legendre_gauss_lobatto_tensor_product_nodes",
+        "padua_jacobi_nodes", "padua_nodes",
         "node_tuples_for_space",
         "edge_clustered_nodes_for_space", "equispaced_nodes_for_space",
         "random_nodes_for_shape",
@@ -107,9 +110,6 @@
         "WitherdenVincentQuadrature",
         ]
 
-from pytools import MovedFunctionDeprecationWrapper
-
-
 get_simplex_onb = MovedFunctionDeprecationWrapper(simplex_onb)
 get_grad_simplex_onb = MovedFunctionDeprecationWrapper(grad_simplex_onb)
 get_warp_and_blend_nodes = MovedFunctionDeprecationWrapper(warp_and_blend_nodes)

modepy/nodes.py

@@ -27,6 +27,9 @@
 
 .. autofunction:: tensor_product_nodes
 .. autofunction:: legendre_gauss_lobatto_tensor_product_nodes
+
+.. autofunction:: padua_jacobi_nodes
+.. autofunction:: padua_nodes
 """
 
 __copyright__ = "Copyright (C) 2009, 2010, 2013 Andreas Kloeckner, " \
@@ -397,6 +400,103 @@ def legendre_gauss_lobatto_tensor_product_nodes(dims: int, n: int) -> np.ndarray
 # }}}
 
 
+# {{{ Padua nodes
+
+def _make_padua_grid_nodes(
+        alpha: float, beta: float, order: int
+        ) -> Tuple[np.ndarray, np.ndarray]:
+    from modepy.quadrature.jacobi_gauss import jacobi_gauss_lobatto_nodes
+    mu = jacobi_gauss_lobatto_nodes(alpha, beta, order)
+    eta = jacobi_gauss_lobatto_nodes(alpha, beta, order + 1)
+
+    return mu, eta
+
+
+def _make_padua_jacobi_nodes(
+        mu: np.ndarray, eta: np.ndarray, odd_or_even: int
+        ) -> np.ndarray:
+    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)
+
+    return nodes[:, indices % 2 == odd_or_even].reshape(2, -1)
+
+
+def _first_padua_jacobi_nodes(alpha: float, beta: float, order: int) -> np.ndarray:
+    mu, eta = _make_padua_grid_nodes(alpha, beta, order)
+    return _make_padua_jacobi_nodes(mu, eta, 0)
+
+
+def _second_padua_jacobi_nodes(alpha: float, beta: float, order: int) -> np.ndarray:
+    # 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 _third_padua_jacobi_nodes(alpha: float, beta: float, order: int) -> np.ndarray:
+    # 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: float, beta: float, order: int) -> np.ndarray:
+    # 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: float, beta: float, order: int,
+        family: str = "first") -> np.ndarray:
+    r"""Generalized Padua-Jacobi nodes.
+
+    The Padua-Jacobi nodes are constructed from an interlaced grid of
+    standard Jacobi-Gauss-Lobatto nodes, making use of
+    :func:`~modepy.quadrature.jacobi_gauss.jacobi_gauss_lobatto_nodes`.
+    This construction is detailed in
+
+        M. Briani, A. Sommariva, M. Vianello,
+        *Computing Fekete and Lebesgue Points: Simplex, Square, Disk*,
+        Journal of Computational and Applied Mathematics, Vol. 236,
+        pp. 2477--2486, 2012, `DOI <http://dx.doi.org/10.1016/j.cam.2011.12.006>`_.
+
+    The values of the parameters :math:`(\alpha, \beta)` can have an effect
+    on the Lebesgue constant of the resulting set, but all of them have
+    optimal growth of :math:`\mathcal{O}(\log^2 n)`.
+
+    The Padua-Jacobi nodes are not rotationally symmetric.
+
+    :arg family: one of the four families of Padua-Jacobi nodes. The three
+        additional families are :math:`90^\circ` rotations of the first one.
+    """
+
+    if family == "first":
+        nodes = _first_padua_jacobi_nodes(alpha, beta, order)
+    elif family == "second":
+        nodes = _second_padua_jacobi_nodes(alpha, beta, order)
+    elif family == "third":
+        nodes = _third_padua_jacobi_nodes(alpha, beta, order)
+    elif family == "fourth":
+        nodes = _fourth_padua_jacobi_nodes(alpha, beta, order)
+    else:
+        raise ValueError(f"unknown Padua-Jacobi node family: '{family}'")
+
+    return nodes
+
+
+def padua_nodes(order: int, family: str = "first") -> np.ndarray:
+    r"""Standard Padua nodes.
+
+    Padua nodes are Padua-Jacobi nodes with :math:`\alpha = \beta = -0.5`,
+    i.e. they are constructed from the Chebyshev-Gauss-Lobatto nodes.
+    """
+
+    return padua_jacobi_nodes(-0.5, -0.5, order, family=family)
+
+# }}}
+
+
 # {{{ space-based interface
 
 @singledispatch