Fix support for zero-dimensional tensor product discretizations

modepy/__init__.py

@@ -28,12 +28,8 @@
     nodal_quad_mass_matrix_for_face, resampling_matrix, vandermonde)
 from modepy.modes import (
     Basis, BasisNotOrthonormal, TensorProductBasis, basis_for_space, grad_jacobi,
-    grad_legendre_tensor_product_basis, grad_simplex_best_available_basis,
-    grad_simplex_monomial_basis, grad_simplex_onb, grad_tensor_product_basis, jacobi,
-    legendre_tensor_product_basis, monomial_basis_for_space,
-    orthonormal_basis_for_space, simplex_best_available_basis,
-    simplex_monomial_basis, simplex_monomial_basis_with_mode_ids, simplex_onb,
-    simplex_onb_with_mode_ids, symbolicize_function, tensor_product_basis)
+    jacobi, monomial_basis_for_space, orthonormal_basis_for_space,
+    symbolicize_function)
 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,
@@ -112,6 +108,4 @@
 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/modes.py

@@ -27,12 +27,11 @@
 from typing import (
     TYPE_CHECKING, Callable, Hashable, Iterable, List, Optional, Sequence, Tuple,
     TypeVar, Union)
-from warnings import warn
 
 import numpy as np
 
 from modepy.shapes import Shape, Simplex, TensorProductShape
-from modepy.spaces import PN, QN, FunctionSpace, TensorProductSpace
+from modepy.spaces import PN, FunctionSpace, TensorProductSpace
 
 
 if TYPE_CHECKING:
@@ -465,193 +464,6 @@ def diff_monomial(r, o):
 # }}}
 
 
-# {{{ DEPRECATED dimension-independent interface for simplices
-
-def simplex_onb_with_mode_ids(dims, n):
-    """Return a list of orthonormal basis functions in dimension *dims* of maximal
-    total degree *n*.
-
-    :returns: a tuple ``(mode_ids, basis)``, where *basis* is a class:`tuple`
-        of functions, each of  which accepts arrays of shape *(dims, npts)* and
-        return the function values as an array of size *npts*.  'Scalar'
-        evaluation, by passing just one vector of length *dims*, is also supported.
-        *mode_ids* is a tuple of the same length as *basis*, where each entry
-        is a tuple of integers describing the order of the mode along the coordinate
-        axes.
-
-    See the following publications:
-
-    * |proriol-ref|
-    * |koornwinder-ref|
-    * |dubiner-ref|
-
-    .. versionadded:: 2018.1
-    """
-    warn("simplex_onb_with_mode_ids is deprecated. "
-            "Use orthonormal_basis_for_space instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    if dims == 1:
-        mode_ids = tuple(range(n+1))
-        return mode_ids, tuple(partial(jacobi, 0, 0, i) for i in mode_ids)
-    else:
-        b = _SimplexONB(dims, n)
-        return b.mode_ids, b.functions
-
-
-def simplex_onb(dims, n):
-    """Return a list of orthonormal basis functions in dimension *dims* of maximal
-    total degree *n*.
-
-    :returns: a :class:`tuple` of functions, each of  which
-        accepts arrays of shape *(dims, npts)*
-        and return the function values as an array of size *npts*.
-        'Scalar' evaluation, by passing just one vector of length *dims*,
-        is also supported.
-
-    See the following publications:
-
-    * |proriol-ref|
-    * |koornwinder-ref|
-    * |dubiner-ref|
-
-    .. versionchanged:: 2013.2
-
-        Made return value a tuple, to make bases hashable.
-    """
-    warn("simplex_onb is deprecated. "
-            "Use orthonormal_basis_for_space instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    mode_ids, basis = simplex_onb_with_mode_ids(dims, n)
-    return basis
-
-
-def grad_simplex_onb(dims, n):
-    """Return the gradients of the functions returned by :func:`simplex_onb`.
-
-    :returns: a :class:`tuple` of functions, each of  which
-        accepts arrays of shape *(dims, npts)*
-        and returns a :class:`tuple` of length *dims* containing
-        the derivatives along each axis as an array of size *npts*.
-        'Scalar' evaluation, by passing just one vector of length *dims*,
-        is also supported.
-
-    See the following publications:
-
-    * |proriol-ref|
-    * |koornwinder-ref|
-    * |dubiner-ref|
-
-    .. versionchanged:: 2013.2
-
-        Made return value a tuple, to make bases hashable.
-    """
-    warn("grad_simplex_onb is deprecated. "
-            "Use orthonormal_basis_for_space instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    from pytools import (
-        generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam)
-
-    if dims == 1:
-        return tuple(partial(grad_jacobi, 0, 0, i) for i in range(n+1))
-    elif dims == 2:
-        return tuple(partial(grad_pkdo_2d, order) for order in gnitstam(n, dims))
-    elif dims == 3:
-        return tuple(partial(grad_pkdo_3d, order) for order in gnitstam(n, dims))
-    else:
-        raise NotImplementedError("%d-dimensional bases" % dims)
-
-
-def simplex_monomial_basis_with_mode_ids(dims, n):
-    """Return a list of monomial basis functions in dimension *dims* of maximal
-    total degree *n*.
-
-    :returns: a tuple ``(mode_ids, basis)``, where *basis* is a class:`tuple`
-        of functions, each of  which accepts arrays of shape *(dims, npts)* and
-        return the function values as an array of size *npts*.  'Scalar'
-        evaluation, by passing just one vector of length *dims*, is also supported.
-        *mode_ids* is a tuple of the same length as *basis*, where each entry
-        is a tuple of integers describing the order of the mode along the coordinate
-        axes.
-
-    .. versionadded:: 2018.1
-    """
-    warn("simplex_monomial_basis_with_mode_ids is deprecated. "
-            "Use monomial_basis_for_space instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    from modepy.nodes import node_tuples_for_space
-    mode_ids = node_tuples_for_space(PN(dims, n))
-
-    return mode_ids, tuple(partial(monomial, order) for order in mode_ids)
-
-
-def simplex_monomial_basis(dims, n):
-    """Return a list of monomial basis functions in dimension *dims* of maximal
-    total degree *n*.
-
-    :returns: a :class:`tuple` of functions, each of  which
-        accepts arrays of shape *(dims, npts)*
-        and return the function values as an array of size *npts*.
-        'Scalar' evaluation, by passing just one vector of length *dims*,
-        is also supported.
-
-    .. versionadded:: 2016.1
-    """
-    mode_ids, basis = simplex_monomial_basis_with_mode_ids(dims, n)
-    return basis
-
-
-def grad_simplex_monomial_basis(dims, n):
-    """Return the gradients of the functions returned by
-    :func:`simplex_monomial_basis`.
-
-    :returns: a :class:`tuple` of functions, each of  which
-        accepts arrays of shape *(dims, npts)*
-        and returns a :class:`tuple` of length *dims* containing
-        the derivatives along each axis as an array of size *npts*.
-        'Scalar' evaluation, by passing just one vector of length *dims*,
-        is also supported.
-
-    .. versionadded:: 2016.1
-    """
-
-    warn("grad_simplex_monomial_basis_with_mode_ids is deprecated. "
-            "Use monomial_basis_for_space instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    from pytools import (
-        generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam)
-    return tuple(partial(grad_monomial, order) for order in gnitstam(n, dims))
-
-
-def simplex_best_available_basis(dims, n):
-    warn("simplex_best_available_basis is deprecated. "
-            "Use basis_for_space instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    return basis_for_space(PN(dims, n), Simplex(dims)).functions
-
-
-def grad_simplex_best_available_basis(dims, n):
-    warn("grad_simplex_best_available_basis is deprecated. "
-            "Use basis_for_space instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    return basis_for_space(PN(dims, n), Simplex(dims)).gradients
-
-# }}}
-
-
 # {{{ tensor product basis helpers
 
 class _TensorProductBasisFunction:
@@ -694,7 +506,17 @@ def __call__(self, x):
         assert x.shape[0] == self.ndim
 
         n = 0
-        result = 1
+
+        if x.shape[1:]:
+            # This ensures the result has the right shape even
+            # if we're in zero dimensions.
+            # The dtype gymnastics are intended as a way to coerce the
+            # dtype to be floating point without changing the accuracy.
+            result = np.ones(x.shape[1:], dtype=(x+np.float32(1)).dtype)
+        else:
+            # Likely we're evaluating symbolically.
+            result = 1
+
         for d, function in zip(self.dims_per_function, self.functions):
             result *= function(x[n:n + d])
             n += d
@@ -799,82 +621,6 @@ def __repr__(self):
 # }}}
 
 
-# {{{ DEPRECATED dimension-independent basis getters
-
-def tensor_product_basis(dims, basis_1d):
-    """Adapt any iterable *basis_1d* of 1D basis functions into a *dims*-dimensional
-    tensor product basis.
-
-    :returns: a tuple of callables representing a *dims*-dimensional basis
-
-    .. versionadded:: 2017.1
-    """
-    warn("tensor_product_basis is deprecated. "
-            "Use TensorProductBasis instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    from modepy.nodes import node_tuples_for_space
-    mode_ids = node_tuples_for_space(QN(dims, len(basis_1d) - 1))
-
-    return tuple(
-            _TensorProductBasisFunction(
-                order, [basis_1d[i] for i in order],
-                dims_per_function=(1,) * len(order))
-            for order in mode_ids)
-
-
-def grad_tensor_product_basis(dims, basis_1d, grad_basis_1d):
-    """Provides derivatives for each of the basis functions generated by
-    :func:`tensor_product_basis`.
-
-    :returns: a :class:`tuple` of callables, where each one returns a
-        *dims*-dimensional :class:`tuple`, one for each derivative.
-
-    .. versionadded:: 2020.2
-    """
-    warn("grad_tensor_product_basis is deprecated. "
-            "Use TensorProductBasis instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    from pytools import wandering_element
-
-    from modepy.nodes import node_tuples_for_space
-    mode_ids = node_tuples_for_space(QN(dims, len(basis_1d) - 1))
-
-    func = (basis_1d, grad_basis_1d)
-    return tuple(
-            _TensorProductGradientBasisFunction(order, [
-                [func[i][k] for i, k in zip(iderivative, order)]
-                for iderivative in wandering_element(dims)
-                ], dims_per_function=(1,) * len(order))
-            for order in mode_ids)
-
-
-def legendre_tensor_product_basis(dims, order):
-    warn("legendre_tensor_product_basis is deprecated. "
-            "Use orthonormal_basis_for_space instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    basis = [partial(jacobi, 0, 0, n) for n in range(order + 1)]
-    return tensor_product_basis(dims, basis)
-
-
-def grad_legendre_tensor_product_basis(dims, order):
-    warn("grad_legendre_tensor_product_basis is deprecated. "
-            "Use orthonormal_basis_for_space instead. "
-            "This function will go away in 2022.",
-            DeprecationWarning, stacklevel=2)
-
-    basis = [partial(jacobi, 0, 0, n) for n in range(order + 1)]
-    grad_basis = [partial(grad_jacobi, 0, 0, n) for n in range(order + 1)]
-    return grad_tensor_product_basis(dims, basis, grad_basis)
-
-# }}}
-
-
 # {{{ conversion to symbolic
 
 def symbolicize_function(

modepy/nodes.py

@@ -386,7 +386,10 @@ def tensor_product_nodes(
         result[d:d+node_dims] = nodes.reshape(nodes.shape + (1,)*(dims-node_dims-d))
         d += node_dims
 
-    return result.reshape((dims, -1), order="F").copy(order="C")
+    if dims == 0:
+        return result.reshape(dims, 1)
+    else:
+        return result.reshape((dims, -1), order="F").copy(order="C")
 
 
 def legendre_gauss_lobatto_tensor_product_nodes(dims: int, n: int) -> np.ndarray:

modepy/test/test_modes.py

@@ -308,6 +308,18 @@ def test_modal_coeffs_by_projection(dim):
 # }}}
 
 
+def test_tp_0d():
+    basis = mp.TensorProductBasis([])
+    nodes = np.zeros((0, 15))
+    for f in basis.functions:
+        res = f(nodes)
+        assert res.shape == (15,)
+
+    for f in basis.gradients:
+        res = f(nodes)
+        assert len(res) == 0
+
+
 # You can test individual routines by typing
 # $ python test_modes.py 'test_routine()'
 

modepy/test/test_nodes.py

@@ -283,6 +283,19 @@ def test_random_nodes_for_tensor_product():
 # }}}
 
 
+def test_tp_0d():
+    import modepy as mp
+    shape = mp.Hypercube(0)
+    space = mp.QN(0, 5)
+
+    for node_func in [
+        mp.equispaced_nodes_for_space,
+        mp.edge_clustered_nodes_for_space,
+    ]:
+        nodes = node_func(space, shape)
+        assert nodes.shape == (0, 1)
+
+
 # You can test individual routines by typing
 # $ python test_nodes.py 'test_routine()'
 

modepy/tools.py

@@ -316,7 +316,7 @@ def plot_element_values(
 
     if resample_n is not None:
         import modepy as mp
-        basis = mp.simplex_onb(dims, n)
+        basis = mp.orthonormal_basis_for_space(mp.PN(dims, n), mp.Simplex(dims))
         fine_nodes = mp.equidistant_nodes(dims, resample_n)
 
         values = np.dot(mp.resampling_matrix(basis, fine_nodes, nodes), values)