diff --git a/.basedpyright/baseline.json b/.basedpyright/baseline.json index f4c0af4..ea4b8f0 100644 --- a/.basedpyright/baseline.json +++ b/.basedpyright/baseline.json @@ -438,6 +438,192 @@ } } ], + "./examples/plot-padua-nodes.py": [ + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 10, + "endColumn": 20, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 4, + "endColumn": 11, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 4, + "endColumn": 11, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 12, + "endColumn": 20, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 22, + "endColumn": 30, + "lineCount": 1 + } + }, + { + "code": "reportUnknownVariableType", + "range": { + "startColumn": 11, + "endColumn": 13, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 27, + "endColumn": 34, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 27, + "endColumn": 34, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 8, + "endColumn": 15, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 16, + "endColumn": 24, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 26, + "endColumn": 34, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 8, + "endColumn": 15, + "lineCount": 1 + } + }, + { + "code": "reportCallIssue", + "range": { + "startColumn": 21, + "endColumn": 27, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 4, + "endColumn": 15, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 4, + "endColumn": 12, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 21, + "endColumn": 31, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 21, + "endColumn": 31, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 33, + "endColumn": 42, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 33, + "endColumn": 42, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 44, + "endColumn": 54, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 44, + "endColumn": 54, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 56, + "endColumn": 67, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 56, + "endColumn": 67, + "lineCount": 1 + } + } + ], "./examples/plot-quadrature-nodes.py": [ { "code": "reportMissingTypeStubs", @@ -3795,6 +3981,22 @@ "lineCount": 1 } }, + { + "code": "reportAny", + "range": { + "startColumn": 4, + "endColumn": 11, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 11, + "endColumn": 55, + "lineCount": 1 + } + }, { "code": "reportUnusedParameter", "range": { diff --git a/examples/plot-padua-nodes.py b/examples/plot-padua-nodes.py new file mode 100644 index 0000000..a92c97d --- /dev/null +++ b/examples/plot-padua-nodes.py @@ -0,0 +1,75 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +import matplotlib.pyplot as plt +import numpy as np + +import modepy as mp + + +if TYPE_CHECKING: + from modepy.typing import ArrayF + + +def _first_padua_curve(n: int, t: ArrayF) -> ArrayF: + return np.vstack([-np.cos((n + 1) * t), -np.cos(n * t)]) + + +def _second_padua_curve(n: int, t: ArrayF) -> ArrayF: + return np.vstack([-np.cos(n * t), -np.cos((n + 1) * t)]) + + +def _third_padua_curve(n: int, t: ArrayF) -> ArrayF: + return np.vstack([np.cos((n + 1) * t), np.cos(n * t)]) + + +def _fourth_padua_curve(n: int, t: ArrayF) -> ArrayF: + return np.vstack([np.cos(n * t), np.cos((n + 1) * t)]) + + +def plot_padua_nodes(alpha: float, beta: float, order: int, family: str) -> None: + 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 + else: + raise ValueError(f"Unknown Padua curve family: '{family}'") + + 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 6f22b14..6beaf0a 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -61,6 +61,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, @@ -168,6 +170,8 @@ "nodal_quadrature_test_matrix", "node_tuples_for_space", "orthonormal_basis_for_space", + "padua_jacobi_nodes", + "padua_nodes", "quadrature_for_space", "random_nodes_for_shape", "resampling_matrix", diff --git a/modepy/nodes.py b/modepy/nodes.py index 434564a..4df654b 100644 --- a/modepy/nodes.py +++ b/modepy/nodes.py @@ -32,6 +32,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 """ from __future__ import annotations @@ -445,6 +448,101 @@ def legendre_gauss_lobatto_tensor_product_nodes(dims: int, n: int) -> ArrayF: # }}} +# {{{ Padua nodes + +def _make_padua_grid_nodes( + alpha: float, beta: float, order: int + ) -> tuple[ArrayF, ArrayF]: + 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: ArrayF, eta: ArrayF, odd_or_even: int) -> ArrayF: + 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) -> ArrayF: + 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) -> ArrayF: + # 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) -> ArrayF: + # 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) -> ArrayF: + # 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") -> ArrayF: + 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") -> ArrayF: + 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