Refactor shapes

doc/nodes.rst

@@ -1,39 +1,4 @@
 Interpolation Nodes
 ===================
 
-Simplices
-^^^^^^^^^
-
-Transformations between coordinate systems
-------------------------------------------
-
-.. currentmodule:: modepy.tools
-
-All of these expect and return arrays of shape *(dims, npts)*.
-
-.. autofunction:: equilateral_to_unit
-.. autofunction:: barycentric_to_unit
-.. autofunction:: unit_to_barycentric
-.. autofunction:: barycentric_to_equilateral
-
-Node sets for interpolation
----------------------------
-
-.. currentmodule:: modepy
-
-.. autofunction:: equidistant_nodes
-.. autofunction:: warp_and_blend_nodes
-
-Also see :class:`modepy.VioreanuRokhlinSimplexQuadrature` if nodes on the
-boundary are not required.
-
-Hypercubes
-^^^^^^^^^^
-
-Node sets for interpolation
----------------------------
-
-.. currentmodule:: modepy
-
-.. autofunction:: tensor_product_nodes
-.. autofunction:: legendre_gauss_lobatto_tensor_product_nodes
+.. automodule:: modepy.nodes

doc/quadrature.rst

@@ -6,10 +6,6 @@ Base classes
 
 .. automodule:: modepy.quadrature
 
-.. currentmodule:: modepy
-
-.. autoclass:: Quadrature
-
 Jacobi-Gauss quadrature in one dimension
 ----------------------------------------
 

doc/tools.rst

@@ -24,11 +24,6 @@ Modal decay/residual
 
 .. automodule:: modepy.modal_decay
 
-Interpolation quality
----------------------
-
-.. currentmodule:: modepy.tools
-
-.. autofunction:: estimate_lebesgue_constant
+.. automodule:: modepy.tools
 
 .. vim: sw=4

modepy/__init__.py

@@ -22,28 +22,43 @@
 """
 
 
+from modepy.shapes import (
+        Shape, Face, Simplex, Hypercube,
+
+        biunit_vertices_for_shape, faces_for_shape
+        )
 from modepy.modes import (
         jacobi, grad_jacobi,
         simplex_onb, grad_simplex_onb, simplex_onb_with_mode_ids,
         simplex_monomial_basis, grad_simplex_monomial_basis,
         simplex_monomial_basis_with_mode_ids,
         simplex_best_available_basis, grad_simplex_best_available_basis,
         tensor_product_basis, grad_tensor_product_basis,
-        legendre_tensor_product_basis, grad_legendre_tensor_product_basis)
+        legendre_tensor_product_basis, grad_legendre_tensor_product_basis,
+
+        Basis, BasisNotOrthonormal,
+        basis_for_shape, orthonormal_basis_for_shape, monomial_basis_for_shape)
 from modepy.nodes import (
         equidistant_nodes, warp_and_blend_nodes,
-        tensor_product_nodes, legendre_gauss_lobatto_tensor_product_nodes)
+        tensor_product_nodes, legendre_gauss_lobatto_tensor_product_nodes,
+
+        node_count_for_shape, node_tuples_for_shape,
+        equispaced_nodes_for_shape, edge_clustered_nodes_for_shape,
+        random_nodes_for_shape)
 from modepy.matrices import (vandermonde,
         resampling_matrix, differentiation_matrices,
         diff_matrix_permutation,
         inverse_mass_matrix, mass_matrix,
-        modal_face_mass_matrix, nodal_face_mass_matrix)
+        modal_face_mass_matrix, nodal_face_mass_matrix,
+        modal_mass_matrix_for_face, nodal_mass_matrix_for_face)
 from modepy.quadrature import (
         Quadrature, QuadratureRuleUnavailable,
-        TensorProductQuadrature, LegendreGaussTensorProductQuadrature)
+        TensorProductQuadrature, LegendreGaussTensorProductQuadrature,
+        quadrature_for_shape)
 from modepy.quadrature.jacobi_gauss import (
         JacobiGaussQuadrature, LegendreGaussQuadrature, ChebyshevGaussQuadrature,
-        GaussGegenbauerQuadrature)
+        GaussGegenbauerQuadrature,
+        )
 from modepy.quadrature.xiao_gimbutas import XiaoGimbutasSimplexQuadrature
 from modepy.quadrature.vioreanu_rokhlin import VioreanuRokhlinSimplexQuadrature
 from modepy.quadrature.grundmann_moeller import GrundmannMoellerSimplexQuadrature
@@ -58,23 +73,35 @@
 __all__ = [
         "__version__",
 
+        "Shape", "Face", "Simplex", "Hypercube",
+        "biunit_vertices_for_shape", "faces_for_shape",
+
         "jacobi", "grad_jacobi",
         "simplex_onb", "grad_simplex_onb", "simplex_onb_with_mode_ids",
         "simplex_monomial_basis", "grad_simplex_monomial_basis",
         "simplex_monomial_basis_with_mode_ids",
         "simplex_best_available_basis", "grad_simplex_best_available_basis",
         "tensor_product_basis", "grad_tensor_product_basis",
         "legendre_tensor_product_basis", "grad_legendre_tensor_product_basis",
+        "Basis", "BasisNotOrthonormal",
+        "basis_for_shape", "orthonormal_basis_for_shape", "monomial_basis_for_shape",
 
         "equidistant_nodes", "warp_and_blend_nodes",
         "tensor_product_nodes", "legendre_gauss_lobatto_tensor_product_nodes",
+        "node_count_for_shape", "node_tuples_for_shape",
+        "edge_clustered_nodes_for_shape", "equispaced_nodes_for_shape",
+        "random_nodes_for_shape",
 
         "vandermonde", "resampling_matrix", "differentiation_matrices",
         "diff_matrix_permutation",
-        "inverse_mass_matrix", "mass_matrix", "modal_face_mass_matrix",
-        "nodal_face_mass_matrix",
+        "inverse_mass_matrix", "mass_matrix",
+        "modal_face_mass_matrix", "nodal_face_mass_matrix",
+        "modal_mass_matrix_for_face", "nodal_mass_matrix_for_face",
 
         "Quadrature", "QuadratureRuleUnavailable",
+        "TensorProductQuadrature", "LegendreGaussTensorProductQuadrature",
+        "quadrature_for_shape",
+
         "JacobiGaussQuadrature", "LegendreGaussQuadrature",
         "GaussLegendreQuadrature", "ChebyshevGaussQuadrature",
         "GaussGegenbauerQuadrature",
@@ -83,7 +110,6 @@
         "ClenshawCurtisQuadrature", "FejerQuadrature",
 
         "WitherdenVincentQuadrature",
-        "TensorProductQuadrature", "LegendreGaussTensorProductQuadrature",
         ]
 
 from pytools import MovedFunctionDeprecationWrapper

modepy/matrices.py

@@ -21,9 +21,12 @@
 """
 
 
+from warnings import warn
 import numpy as np
 import numpy.linalg as la
 
+import modepy.shapes as shp
+
 
 __doc__ = r"""
 .. currentmodule:: modepy
@@ -48,12 +51,9 @@
 where :math:`(\phi_i)_i` is the basis of functions underlying :math:`V`.
 
 .. autofunction:: inverse_mass_matrix
-
 .. autofunction:: mass_matrix
-
-.. autofunction:: modal_face_mass_matrix
-
-.. autofunction:: nodal_face_mass_matrix
+.. autofunction:: modal_mass_matrix_for_face
+.. autofunction:: nodal_mass_matrix_for_face
 
 Differentiation is also convenient to express by using :math:`V^{-1}` to
 obtain modal values and then using a Vandermonde matrix for the derivatives
@@ -75,7 +75,7 @@ def vandermonde(functions, nodes):
 
     *functions* are allowed to return :class:`tuple` instances.
     In this case, a tuple of matrices is returned--i.e. this function
-    works directly on :func:`modepy.grad_simplex_onb` and returns
+    works directly on :func:`modepy.Basis.gradients` and returns
     a tuple of matrices.
     """
 
@@ -110,7 +110,7 @@ def resampling_matrix(basis, new_nodes, old_nodes, least_squares_ok=False):
 
     :arg basis: A sequence of basis functions accepting
         arrays of shape *(dims, npts)*, like those returned by
-        :func:`modepy.simplex_onb`.
+        :func:`modepy.orthonormal_basis_for_shape`.
     :arg new_nodes: An array of shape *(dims, n_new_nodes)*
     :arg old_nodes: An array of shape *(dims, n_old_nodes)*
     :arg least_squares_ok: If *False*, then nodal values at *old_nodes*
@@ -148,8 +148,9 @@ def is_square(m):
                     np.dot(resample_vdm_new, la.pinv(vdm_old)),
                     order="C")
         else:
-            raise RuntimeError("number of input nodes and number "
-                    "of basis functions "
+            raise RuntimeError(
+                    f"number of input nodes ({old_nodes.shape[1]}) "
+                    f"and number of basis functions ({len(basis)}) "
                     "do not agree--perhaps use least_squares_ok")
 
 
@@ -160,10 +161,10 @@ def differentiation_matrices(basis, grad_basis, nodes, from_nodes=None):
 
     :arg basis: A sequence of basis functions accepting arrays
         of shape *(dims, npts)*,
-        like those returned by :func:`modepy.simplex_onb`.
+        like those returned by :func:`modepy.orthonormal_basis_for_shape`.
     :arg grad_basis: A sequence of functions returning the
         gradients of *basis*,
-        like those returned by :func:`modepy.grad_simplex_onb`.
+        like those in :attr:`modepy.Basis.gradients`.
     :arg nodes: An array of shape *(dims, n_nodes)*
     :arg from_nodes:  An array of shape *(dims, n_from_nodes)*.
         If *None*, assumed to be the same as *nodes*.
@@ -242,34 +243,68 @@ def mass_matrix(basis, nodes):
     return la.inv(inverse_mass_matrix(basis, nodes))
 
 
-def modal_face_mass_matrix(trial_basis, order, face_vertices,
-        test_basis=None, shape=None):
+def modal_mass_matrix_for_face(face, trial_functions, test_functions, order):
+    """
+    .. versionadded:: 2020.3
+    """
+
+    from modepy.quadrature import quadrature_for_shape
+    quad = quadrature_for_shape(face, order)
+
+    assert quad.exact_to > order*2
+    mapped_nodes = face.map_to_volume(quad.nodes)
+
+    result = np.empty((len(test_functions), len(trial_functions)))
+
+    for i, test_f in enumerate(test_functions):
+        test_vals = test_f(mapped_nodes)
+        for j, trial_f in enumerate(trial_functions):
+            result[i, j] = (test_vals*trial_f(quad.nodes)).dot(quad.weights)
+
+    return result
+
+
+def nodal_mass_matrix_for_face(face: shp.Face, trial_functions, test_functions,
+        volume_nodes, face_nodes, order):
+    """
+    .. versionadded :: 2020.3
+    """
+    face_vdm = vandermonde(trial_functions, face_nodes)
+    vol_vdm = vandermonde(test_functions, volume_nodes)
+
+    modal_fmm = modal_mass_matrix_for_face(
+            face, trial_functions, test_functions, order)
+    return la.inv(vol_vdm.T).dot(modal_fmm).dot(la.pinv(face_vdm))
+
+
+# {{{ deprecated junk
+
+def modal_face_mass_matrix(trial_basis, order, face_vertices, test_basis=None):
     """
     :arg face_vertices: an array of shape ``(dims, nvertices)``.
     :arg shape: a :class:`~modepy.shapes.Shape` that identifies the
-        reference face element.
+        reference volume element.
 
     .. versionadded :: 2016.1
-
-    .. versionchanged:: 2020.3
-
-        Added *shape* parameter and support for :math:`[-1, 1]^d` domains.
     """
 
+    warn("modal_face_mass_matrix is deprecated and will go away in 2022. "
+            "Use modal_mass_matrix_for_face instead.",
+            DeprecationWarning, stacklevel=2)
+
     if test_basis is None:
         test_basis = trial_basis
 
-    if shape is None:
-        from modepy.shapes import Simplex
-        shape = Simplex(face_vertices.shape[0])
+    vol_dims = face_vertices.shape[0]
+    face_shape = shp.Simplex(vol_dims - 1)
 
-    from modepy.shapes import get_face_map, get_quadrature
-    face = type(shape)(shape.dims - 1)
-    fmap = get_face_map(shape, face_vertices)
-    quad = get_quadrature(face, order)
+    from modepy.quadrature import quadrature_for_shape
+    quad = quadrature_for_shape(face_shape, order)
 
     assert quad.exact_to > order*2
-    mapped_nodes = fmap(quad.nodes)
+
+    from modepy.shapes import _simplex_face_to_vol_map
+    mapped_nodes = _simplex_face_to_vol_map(face_vertices, quad.nodes)
 
     nrows = len(test_basis)
     ncols = len(trial_basis)
@@ -278,43 +313,38 @@ def modal_face_mass_matrix(trial_basis, order, face_vertices,
     for i, test_f in enumerate(test_basis):
         test_vals = test_f(mapped_nodes)
         for j, trial_f in enumerate(trial_basis):
-            trial_vals = trial_f(mapped_nodes)
-
-            result[i, j] = (test_vals*trial_vals).dot(quad.weights)
+            result[i, j] = (test_vals*trial_f(mapped_nodes)).dot(quad.weights)
 
     return result
 
 
 def nodal_face_mass_matrix(trial_basis, volume_nodes, face_nodes, order,
-        face_vertices, test_basis=None, shape=None):
+        face_vertices, test_basis=None):
     """
     :arg face_vertices: an array of shape ``(dims, nvertices)``.
     :arg shape: a :class:`~modepy.shapes.Shape` that identifies the
         reference face element.
 
     .. versionadded :: 2016.1
-
-    .. versionchanged:: 2020.3
-
-        Added *shape* parameter and support for :math:`[-1, 1]^d` domains.
     """
 
+    warn("nodal_face_mass_matrix is deprecated and will go away in 2022. "
+            "Use nodal_mass_matrix_for_face instead.",
+            DeprecationWarning, stacklevel=2)
+
     if test_basis is None:
         test_basis = trial_basis
 
-    if shape is None:
-        from modepy.shapes import Simplex
-        shape = Simplex(face_vertices.shape[0])
-
-    from modepy.shapes import get_face_map
-    fmap = get_face_map(shape, face_vertices)
-    face_vdm = vandermonde(trial_basis, fmap(face_nodes))  # /!\ non-square
+    from modepy.shapes import _simplex_face_to_vol_map
+    face_vdm = vandermonde(
+            trial_basis,
+            _simplex_face_to_vol_map(face_vertices, face_nodes))
     vol_vdm = vandermonde(test_basis, volume_nodes)
 
     modal_fmm = modal_face_mass_matrix(
-            trial_basis, order, face_vertices,
-            test_basis=test_basis, shape=shape)
+            trial_basis, order, face_vertices, test_basis=test_basis)
     return la.inv(vol_vdm.T).dot(modal_fmm).dot(la.pinv(face_vdm))
 
+# }}}
 
 # vim: foldmethod=marker

modepy/modal_decay.py

@@ -25,7 +25,7 @@
 import numpy.linalg as la
 
 __doc__ = """Estimate the smoothness of a function represented in a basis
-returned by :func:`modepy.simplex_onb`.
+returned by :func:`modepy.orthonormal_basis_for_shape`.
 
 The method implemented in this module follows this article: