From 2b58c3f7f3628e5315b52514053e80672221fd06 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 19 Dec 2020 19:44:04 -0600 Subject: [PATCH] add padua nodes --- examples/plot-padua-nodes.py | 65 +++++++++++++++++++++++ modepy/__init__.py | 7 ++- modepy/nodes.py | 100 +++++++++++++++++++++++++++++++++++ 3 files changed, 170 insertions(+), 2 deletions(-) create mode 100644 examples/plot-padua-nodes.py diff --git a/examples/plot-padua-nodes.py b/examples/plot-padua-nodes.py new file mode 100644 index 0000000..5190836 --- /dev/null +++ b/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) diff --git a/modepy/__init__.py b/modepy/__init__.py index fe800f8..1b86eea 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -21,6 +21,7 @@ THE SOFTWARE. """ +from pytools import MovedFunctionDeprecationWrapper from modepy.matrices import ( diff_matrices, @@ -58,6 +59,8 @@ legendre_gauss_lobatto_tensor_product_nodes, legendre_gauss_tensor_product_nodes, node_tuples_for_space, + padua_jacobi_nodes, + padua_nodes, random_nodes_for_shape, tensor_product_nodes, warp_and_blend_nodes, @@ -164,6 +167,8 @@ "nodal_quad_mass_matrix_for_face", "node_tuples_for_space", "orthonormal_basis_for_space", + "padua_jacobi_nodes", + "padua_nodes", "quadrature_for_space", "random_nodes_for_shape", "resampling_matrix", @@ -178,7 +183,5 @@ "warp_and_blend_nodes", ] -from pytools import MovedFunctionDeprecationWrapper - get_warp_and_blend_nodes = MovedFunctionDeprecationWrapper(warp_and_blend_nodes) diff --git a/modepy/nodes.py b/modepy/nodes.py index 84a8838..85be59d 100644 --- a/modepy/nodes.py +++ b/modepy/nodes.py @@ -28,6 +28,9 @@ .. autofunction:: tensor_product_nodes .. autofunction:: legendre_gauss_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, " \ @@ -425,6 +428,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 `_. + + 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