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 8bb5510..62541cc 100644 --- a/modepy/__init__.py +++ b/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) diff --git a/modepy/nodes.py b/modepy/nodes.py index 7f0f72c..4e81e47 100644 --- a/modepy/nodes.py +++ b/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 `_. + + 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