Diagonal TP mass matrices and related changes

README.rst

@@ -1,5 +1,5 @@
-modepy: Basis Functions and Node Sets for Interpolation
-=======================================================
+modepy: Basis Functions, Node Sets, Quadratures
+===============================================
 
 .. image:: https://gitlab.tiker.net/inducer/modepy/badges/main/pipeline.svg
     :alt: Gitlab Build Status
@@ -18,12 +18,41 @@ modepy: Basis Functions and Node Sets for Interpolation
 simplices (i.e. segments, triangles and tetrahedra) and tensor products of
 simplices (i.e. squares, cubes, prisms, etc.). These are a key building block
 for high-order unstructured discretizations, as often used in a finite
-element context. It closely follows the approach taken in the book
+element context. Features include:
+
+- Support for simplex and tensor product elements in any dimension.
+- Orthogonal bases:
+    - Jacobi polynomials with derivatives
+    - Orthogonal polynomials for simplices up to 3D and tensor product elements
+      and their derivatives.
+    - All bases permit symbolic evaluation, for code generation.
+- Access to numerous quadrature rules:
+    - Jacobi-Gauss, Jacobi-Gauss-Lobatto in 1D
+      (includes Legendre, Chebyshev, ultraspherical, Gegenbauer)
+    - Clenshaw-Curtis and Fejér in 1D
+    - Grundmann-Möller on the simplex
+    - Xiao-Gimbutas on the simplex
+    - Vioreanu-Rokhlin on the simplex
+    - Jaśkowiec-Sukumar on the tetrahedron
+    - Witherden-Vincent on the hypercube
+    - Generic tensor products built on the above, e.g. for prisms and hypercubes
+- Matrices for FEM, usable across all element types:
+    - generalized Vandermonde,
+    - mass matrices (including lumped diagonal),
+    - face mass matrices,
+    - differentiation matrices, and
+    - resampling matrices.
+- Objects to represent 'element shape' and 'function space',
+  generic node/mode/quadrature retrieval based on them.
+
+Its roots closely followed the approach taken in the book
 
   Hesthaven, Jan S., and Tim Warburton. "Nodal Discontinuous Galerkin Methods:
   Algorithms, Analysis, and Applications". 1st ed. Springer, 2007.
   `Book web page <http://nudg.org>`_
 
+but much has been added beyond that basic functionality.
+
 Resources:
 
 * `documentation <http://documen.tician.de/modepy>`_

doc/quadrature.rst

@@ -11,27 +11,6 @@ Jacobi-Gauss quadrature in one dimension
 
 .. automodule:: modepy.quadrature.jacobi_gauss
 
-.. currentmodule:: modepy
-
-.. autoclass:: JacobiGaussQuadrature
-    :members:
-    :show-inheritance:
-
-.. autoclass:: LegendreGaussQuadrature
-    :show-inheritance:
-
-.. autoclass:: ChebyshevGaussQuadrature
-    :show-inheritance:
-
-.. autoclass:: GaussGegenbauerQuadrature
-    :show-inheritance:
-
-.. currentmodule:: modepy.quadrature.jacobi_gauss
-
-.. autofunction:: jacobi_gauss_lobatto_nodes
-
-.. autofunction:: legendre_gauss_lobatto_nodes
-
 Clenshaw-Curtis and Fejér quadrature in one dimension
 -----------------------------------------------------
 

modepy/__init__.py

@@ -27,25 +27,28 @@
     inverse_mass_matrix, mass_matrix, modal_mass_matrix_for_face,
     modal_quad_mass_matrix, modal_quad_mass_matrix_for_face, multi_vandermonde,
     nodal_mass_matrix_for_face, nodal_quad_mass_matrix,
-    nodal_quad_mass_matrix_for_face, resampling_matrix, vandermonde)
+    nodal_quad_mass_matrix_for_face, resampling_matrix,
+    spectral_diag_nodal_mass_matrix, vandermonde)
 from modepy.modes import (
     Basis, BasisNotOrthonormal, TensorProductBasis, basis_for_space, grad_jacobi,
-    jacobi, monomial_basis_for_space, orthonormal_basis_for_space,
+    jacobi, monomial_basis_for_space, orthonormal_basis_for_space, scaled_jacobi,
     symbolicize_function)
 from modepy.nodes import (
     edge_clustered_nodes_for_space, equidistant_nodes, equispaced_nodes_for_space,
     legendre_gauss_lobatto_tensor_product_nodes, legendre_gauss_tensor_product_nodes,
     node_tuples_for_space, random_nodes_for_shape, tensor_product_nodes,
     warp_and_blend_nodes)
 from modepy.quadrature import (
+    LegendreGaussLobattoTensorProductQuadrature,
     LegendreGaussTensorProductQuadrature, Quadrature, QuadratureRuleUnavailable,
     TensorProductQuadrature, ZeroDimensionalQuadrature, quadrature_for_space)
 from modepy.quadrature.clenshaw_curtis import (
     ClenshawCurtisQuadrature, FejerQuadrature)
 from modepy.quadrature.grundmann_moeller import GrundmannMoellerSimplexQuadrature
 from modepy.quadrature.jacobi_gauss import (
-    ChebyshevGaussQuadrature, GaussGegenbauerQuadrature, JacobiGaussQuadrature,
-    LegendreGaussQuadrature)
+    ChebyshevGaussQuadrature, GaussGegenbauerQuadrature,
+    JacobiGaussLobattoQuadrature, JacobiGaussQuadrature,
+    LegendreGaussLobattoQuadrature, LegendreGaussQuadrature)
 from modepy.quadrature.jaskowiec_sukumar import JaskowiecSukumarQuadrature
 from modepy.quadrature.vioreanu_rokhlin import VioreanuRokhlinSimplexQuadrature
 from modepy.quadrature.witherden_vincent import WitherdenVincentQuadrature
@@ -68,7 +71,7 @@
 
         "FunctionSpace", "TensorProductSpace", "PN", "QN", "space_for_shape",
 
-        "jacobi", "grad_jacobi",
+        "jacobi", "grad_jacobi", "scaled_jacobi",
         "symbolicize_function",
         "Basis", "BasisNotOrthonormal", "TensorProductBasis",
         "basis_for_space", "orthonormal_basis_for_space", "monomial_basis_for_space",
@@ -86,18 +89,23 @@
         "diff_matrix_permutation",
         "inverse_mass_matrix", "mass_matrix",
         "modal_quad_mass_matrix", "nodal_quad_mass_matrix",
+        "spectral_diag_nodal_mass_matrix",
         "modal_mass_matrix_for_face", "nodal_mass_matrix_for_face",
         "modal_quad_mass_matrix_for_face",
         "nodal_quad_mass_matrix_for_face",
 
         "Quadrature", "QuadratureRuleUnavailable",
         "ZeroDimensionalQuadrature",
         "TensorProductQuadrature", "LegendreGaussTensorProductQuadrature",
+        "LegendreGaussLobattoTensorProductQuadrature",
         "quadrature_for_space",
 
         "JacobiGaussQuadrature", "LegendreGaussQuadrature",
         "GaussLegendreQuadrature", "ChebyshevGaussQuadrature",
         "GaussGegenbauerQuadrature",
+        "JacobiGaussLobattoQuadrature",
+        "LegendreGaussLobattoQuadrature",
+
         "XiaoGimbutasSimplexQuadrature", "GrundmannMoellerSimplexQuadrature",
         "VioreanuRokhlinSimplexQuadrature",
         "ClenshawCurtisQuadrature", "FejerQuadrature",

modepy/matrices.py

@@ -28,7 +28,7 @@
 import numpy.linalg as la
 
 from modepy.modes import Basis, BasisNotOrthonormal
-from modepy.quadrature import Quadrature
+from modepy.quadrature import Quadrature, TensorProductQuadrature
 from modepy.shapes import Face
 
 
@@ -62,6 +62,7 @@
 .. autofunction:: nodal_mass_matrix_for_face
 .. autofunction:: modal_quad_mass_matrix_for_face
 .. autofunction:: nodal_quad_mass_matrix_for_face
+.. autofunction:: spectral_diag_nodal_mass_matrix
 
 Differentiation is also convenient to express by using :math:`V^{-1}` to
 obtain modal values and then using a Vandermonde matrix for the derivatives
@@ -396,6 +397,27 @@ def nodal_quad_mass_matrix(
                     modal_quad_mass_matrix(quadrature, test_functions))
 
 
+def spectral_diag_nodal_mass_matrix(
+            quadrature: TensorProductQuadrature
+        ) -> np.ndarray:
+    """Return the diagonal mass matrix for use in the spectral element method.
+    This mass matrix is exact for Lagrange polynomials with respect to
+    Gauss-Legendre (GL) nodes, using GL nodal degrees of freedom.
+    It is approximate for Lagrange polynomials with respect to the
+    Gauss-Lobatto-Legendre (GLL) nodes, using GLL nodal degrees of freedom.
+
+    Returns the vector of diagonal entries.
+
+    .. versionadded:: 2024.2
+    """
+    if not isinstance(quadrature, TensorProductQuadrature):
+        raise ValueError("only applicable to tensor product discretizations")
+    if not all(quad.dim == 1 for quad in quadrature.quadratures):
+        raise ValueError("constituent quadratures of TP quadrature must be 1D")
+
+    return quadrature.weights
+
+
 def modal_mass_matrix_for_face(
             face: Face, face_quad: Quadrature,
             trial_functions: Sequence[Callable[[np.ndarray], np.ndarray]],

modepy/modes.py

@@ -66,6 +66,7 @@
 
 .. autofunction:: jacobi
 .. autofunction:: grad_jacobi
+.. autofunction:: scaled_jacobi
 
 Conversion to Symbolic
 ----------------------
@@ -209,6 +210,56 @@ def grad_jacobi(alpha: float, beta: float, n: int, x: RealValueT) -> RealValueT:
     else:
         return math.sqrt(n*(n+alpha+beta+1)) * jacobi(alpha+1, beta+1, n-1, x)
 
+
+def binom(x, y):
+    # FIXME This may overflow quickly.
+    # mpmath has clever 'gammaprod' pole handling in case that's ever needed:
+    # https://github.com/mpmath/mpmath/blob/75a2ed37c4f2c576a9d01d360ee4c94ead57c7ff/mpmath/functions/factorials.py#L61-L65
+    from math import gamma
+    return gamma(x + 1) / (gamma(y + 1) * gamma(x - y + 1))
+
+
+def scaled_jacobi(alpha: float, beta: float, n: int, x: RealValueT) -> RealValueT:
+    r"""
+    Same as :func:`jacobi`, but scaled so that
+
+    .. math::
+
+        P_n^{(\alpha, \beta)}(1) = \binom{ n+\alpha \\ n},
+
+    and
+
+    .. math::
+
+        P_n^{(\alpha, \beta)}(-1) = (-1)^n \binom{ n+\beta \\ n}.
+
+    This is the more conventional definition of the Jacobi polynomials.
+    """
+    from math import factorial, gamma, sqrt
+
+    eps = 100 * np.finfo((np.float32(0) + np.asarray(x)).dtype).eps
+    if n == 0 and abs(alpha + beta + 1) < eps:
+        # maxima claims limit(1/(x*gamma(x)),x, 0) == 1
+        scaling = (
+            # 2**(alpha + beta + 1) \approx 1
+            gamma(n + alpha + 1)
+            * gamma(n + beta + 1)
+            / factorial(n)
+        )
+    else:
+        # source:
+        # https://en.wikipedia.org/w/index.php?title=Jacobi_polynomials&oldid=1219118998
+        scaling = (
+            2**(alpha + beta + 1)
+            / (2*n + alpha + beta + 1)
+            * gamma(n + alpha + 1)
+            * gamma(n + beta + 1)
+            / gamma(n + alpha + beta + 1)
+            / factorial(n)
+        )
+
+    return sqrt(scaling) * jacobi(alpha, beta, n, x)
+
 # }}}
 
 

modepy/nodes.py

@@ -395,20 +395,22 @@ def tensor_product_nodes(
 
 def legendre_gauss_tensor_product_nodes(dims: int, n: int) -> np.ndarray:
     """
-    :arg n: the number of points in one dimension.
+    :arg n: the degree of polynomial exactly interpolated by the nodes.
+        The one-dimensional base quadrature has *n+1* nodes.
 
     .. versionadded:: 2024.2
     """
     from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature
-    gl_quad = LegendreGaussQuadrature(n)
+    gl_quad = LegendreGaussQuadrature(n, force_dim_axis=True)
     gl_nodes = gl_quad.nodes.reshape(1, -1)
 
     return tensor_product_nodes(dims, gl_nodes)
 
 
 def legendre_gauss_lobatto_tensor_product_nodes(dims: int, n: int) -> np.ndarray:
     """
-    :arg n: the number of points in one dimension.
+    :arg n: the degree of polynomial exactly interpolated by the nodes.
+        The one-dimensional base quadrature has *n+1* nodes.
     """
     from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes
     return tensor_product_nodes(dims,

modepy/quadrature/__init__.py

@@ -35,7 +35,7 @@
 """
 
 from functools import singledispatch
-from typing import Callable, Iterable, Optional
+from typing import Callable, Iterable, Optional, Sequence
 
 import numpy as np
 
@@ -143,9 +143,15 @@ def __init__(self, quad: Quadrature, left: float, right: float) -> None:
 class TensorProductQuadrature(Quadrature):
     r"""A tensor product quadrature of one-dimensional :class:`Quadrature`\ s.
 
+    .. autoattribute:: quadratures
     .. automethod:: __init__
     """
 
+    quadratures: Sequence[Quadrature]
+    """The lower-dimensional quadratures from which the tensor product quadrature is
+    composed.
+    """
+
     def __init__(self, quads: Iterable[Quadrature]) -> None:
         """
         :arg quad: a iterable of :class:`Quadrature` objects for one-dimensional
@@ -154,7 +160,7 @@ def __init__(self, quads: Iterable[Quadrature]) -> None:
 
         from modepy.nodes import tensor_product_nodes
 
-        quads = list(quads)
+        quads = tuple(quads)
         x = tensor_product_nodes([quad.nodes for quad in quads])
         w = np.prod(tensor_product_nodes([quad.weights for quad in quads]), axis=0)
         assert w.size == x.shape[1]
@@ -167,6 +173,8 @@ def __init__(self, quads: Iterable[Quadrature]) -> None:
 
         super().__init__(x, w, exact_to=exact_to)
 
+        self.quadratures = quads
+
 
 class LegendreGaussTensorProductQuadrature(TensorProductQuadrature):
     """A tensor product using only :class:`~modepy.LegendreGaussQuadrature`
@@ -182,6 +190,20 @@ def __init__(self,
             ] * dims)
 
 
+class LegendreGaussLobattoTensorProductQuadrature(TensorProductQuadrature):
+    """A tensor product using only :class:`~modepy.LegendreGaussLobattoQuadrature`
+    one-dimenisonal rules.
+    """
+
+    def __init__(self,
+                 N: int, dims: int,  # noqa: N803
+                 backend: Optional[str] = None) -> None:
+        from modepy.quadrature.jacobi_gauss import LegendreGaussLobattoQuadrature
+        super().__init__([
+            LegendreGaussLobattoQuadrature(N, backend=backend, force_dim_axis=True)
+            ] * dims)
+
+
 # {{{ quadrature
 
 @singledispatch