Remove example and trim down overview

Author: alexfikl

doc/paper/paper.md

@@ -146,143 +146,80 @@ to any array type (including GPU or other high-performance arrays) as long as
 
 # Overview
 
-The high-level concepts available in `modepy` are shapes (i.e the geometry of a
+The high-level concepts available in `modepy` are shapes (i.e. a reference
 domain), modes (i.e. the basis functions), and nodes (i.e. the degrees of
 freedom). These are implemented in a user-extensible fashion using the
 `singledispatch` mechanism, with inspiration taken from common idiomatic usage
-in Julia [@Bezanson2017]. Given a function space `space` on a `shape`, a set of
-nodal degrees of freedom can be obtained via
-`edge_clustered_nodes_for_space(space, shape)`, with an implementation selected
-based on the type of `space`. While the example shows that multiple dispatch
-could be beneficial, we have decided to forgo it to avoid potential soundness
-issues [@Julia2018].
+in Julia [@Bezanson2017]. 
 
 ## Shapes
 
 The geometry of a reference element is described in `modepy` by the `Shape`
 class. Built-in support exists for `Simplex` and `Hypercube` geometries,
 encompassing the commonly used interval, triangle, tetrahedron, quadrilateral,
 and hexahedral shapes (see \autoref{FigureSimplices}). `TensorProductShape` can
-be used to compose additional shapes, such as prisms, that are useful in
-specific applications and sometimes generated by meshing software (such as
-`gmsh` [@Geuzaine2009]). An example that constructs a prism and performs some
-standard operations on it is given in the next section.
+be used to compose additional shapes, that are useful in specific applications
+and sometimes generated by meshing software (such as `gmsh` [@Geuzaine2009]). 
 
 ![Domains corresponding to the one-, two-, and three-dimensional simplices.](images/simplices.png){#FigureSimplices width="80%"}
 
-To support a wide array of applications, shapes can be queried and manipulated
-in several ways. The `faces_for_shape` function can be used to return a list of
-faces, which are shapes themselves, that allow access to domain normals and
-mappings back to the reference element.
-
 ## Modes and Spaces
 
-To perform calculus operations, each reference element can be equipped with
-function spaces described by the `FunctionSpace` class.
-These represent a finite-dimensional space of functions $\phi_i: D \to
-\mathbb{R}$, where $D$ is the reference element domain, but make no specific
-choice of basis. The function spaces mirror the shape definitions and are
-intended to work in tandem. In fact, the `space_for_shape` function provides a
-default choice of space for each shape. Predefined choices include the `PN`
-space, containing polynomials of total degree at most `order`, and the `QN`
-space, containing polynomials of maximum degree at most `order`. As with
-shapes, these spaces can also be combined using `TensorProductSpace` to
-describe functional spaces for other domains or in higher dimensions.
-
-For calculations, each space can be equipped with a default basis using the
-`basis_for_space` function. `modepy` provides two special types of basis
-functions: orthonormal (accessed through `orthonormal_basis_for_space`) and
-monomial (accessed through `monomial_basis_for_space`). The basis is provided
-in the form of a `Basis` class, which allows access to the basis functions,
-their derivatives (per axis) and opaque `mode_ids`. As noted before, the basis
-functions and their derivatives can be evaluated directly or can be traced
-using the `pymbolic` library to aid in code generation.
+To perform calculus operations, each reference element can be equipped with a
+function space described by the `FunctionSpace` class. These represent a
+finite-dimensional space of functions $\phi_i: D \to \mathbb{R}$, where $D$ is
+the reference element domain, and no specific choice of basis.  Predefined
+choices include the `PN` space, containing polynomials of total degree at most
+`order`, and the `QN` space, containing polynomials of maximum degree at most
+`order`. As with shapes, these spaces can be combined using
+`TensorProductSpace`.
+
+The basis is provided in the form of a `Basis` class, which allows access to
+the basis functions, their derivatives (per axis) and opaque `mode_ids`. As
+noted before, the basis functions and their derivatives can be evaluated
+directly or can be traced using the `pymbolic` library to aid in code
+generation. Various standard basis functions are provided, such as the
+monomials, general Jacobi polynomials, and the
+Proriol-Koornwinder-Dubiner-Owens (PKDO) basis from [@Dubiner1991] (see
+\autoref{FigurePKDO}).
 
 ![PKDO basis functions for the triangle.](images/pkdo-2d.png){#FigurePKDO width="80%"}
 
-By default, the `Simplex` shape is equipped with a `PN` functional space. This
-function space uses the Proriol-Koornwinder-Dubiner-Owens (PKDO) basis (see
-\autoref{FigurePKDO}) as defined in [@Dubiner1991] for dimensions up to 3. The
-`Hypercube` is equipped with the `QN` functional space that uses a tensor
-product of standard Legendre polynomials. Additionally, general Jacobi
-polynomials are available and can be used as a basis for variants of `QN` with
-different weights, depending on the application.
-
 ## Nodes
 
 A final component in an FEM discretization (or 'Ciarlet triple' [@Brenner2007,
-Section 3.1]) is a set of 'degrees of freedom' ('DOFs'), that is, a vector of
-numbers (technically, functionals) that uniquely define a certain element in
-the span of a basis. The literature considers various types of degrees of
-freedom; the two most common are modal DOFs (i.e. basis coefficients) and nodal
-DOFs (i.e. function or derivative values at certain points). `modepy` provides
-easy access to both, while paying particular attention to their usability in
-the high-order setting. Due to Runge's phenomenon [@Trefethen2020], simple
-equispaced nodal DOFs rapidly become unusable with increasing polynomial order
-due to poor conditioning. As such, `modepy` offers access to many types of
-nodes with better properties for each reference element. There also exists an
-ability to construct one's own set of nodes given a set of functions in 2D via
-a procedure based on eigenvalues of multiplication operators [@Vioreanu2014].
-
-On each shape, three default families of nodes are provided: equispaced
-(accessed through `equispaced_nodes_for_space`), edge-clustered (accessed
-through `edge_clustered_nodes_for_space`) and randomly distributed (accessed
-through `random_nodes_for_space`).
-
-For more generality, `modepy` can directly interoperate with the
-`recursivenodes` library described in [@Isaac2020], which offers
-well-conditioned nodes on the simplex that obey a 'recursive' property, where
-nodes on a face of a higher-dimensional shape are exactly the lower-dimensional
-nodes. In addition, it offers many many tunable parameters (e.g. nodes on the
-boundary vs. not). On simplices, the "warp-and-blend" nodes [@Warburton2007]
-are also available for dimensions three and lower. On the hypercube, standard
-tensor product nodes are constructed from one-dimensional
-Legendre-Gauss-Lobatto nodes. As for the basis functions, the full family of
-Jacobi-Gauss and  Jacobi-Gauss-Lobatto quadrature points is available for
-specific applications.
-
-Taken together, the functionality described thus far provides comprehensive
-support for interpolation, that is, the evaluation of a polynomial defined by
-its nodal DOFs, at any point, and further the ability to obtain any derivative
-of the interpolating polynomial.
+Section 3.1]) is a set of 'degrees of freedom' ('DOFs') that uniquely define a
+certain function in the span of a basis. `modepy` provides easy access to modal
+DOFs (i.e. basis coefficients) and nodal DOFs (i.e. function or derivative
+values at a point). As equispaced nodal DOFs are not usable in the high-order
+setting [@Trefethen2020], `modepy` offers access to many types of nodes with
+better properties for each reference element. It can also construct custom
+nodes sets for a given set of basis functions in 2D via a procedure based on
+eigenvalues of multiplication operators [@Vioreanu2014].
+
+On simplices, the "warp-and-blend" nodes [@Warburton2007] are available, and on
+the hypercube, standard tensor product nodes are constructed from
+one-dimensional Legendre-Gauss(-Lobatto) nodes. `modepy` can also directly
+interoperate with the `recursivenodes` library described in [@Isaac2020], which
+offers additional well-conditioned nodes on the simplex.
 
 ## Quadrature
 
-While nodal degrees of freedom must be *unisolvent* or *interpolatory* (i.e.
-uniquely determine a function in a discrete function space), this requirement
-does not exist in the setting of quadrature, where the objective is to
-accurately approximate the value of an integral of a function on a shape. Using
-function values at points $\boldsymbol{x}_i$ ("nodes") and associated weights
-$w_i$, quadrature rules are generically expressed as
-$$
-\int_{D} f(\boldsymbol x) \mathrm{d}\boldsymbol{x} \approx
-    \sum_{k = 0}^N f(\boldsymbol{x}_k) w_k.
-$$
-
-Besides the standard building blocks consisting of domains, basis functions and
-node sets, `modepy` also offers a wide array of quadrature rules that can be
-used on each domain. The quadrature rules are provided as implementations of
-the `Quadrature` superclass.
-
-For the interval, a broad selection of quadrature rules are present:
-Clenshaw--Curtis, Fejér, and Jacobi-Gauss. The latter class includes the
-commonly used Chebyshev and Legendre types, among others. Lobatto (that is,
-endpoint-included) variants of these are also available. The line segment
-quadratures can be used as part of the `TensorProductQuadrature` class to
-define different higher-dimensional rules. Additionally, several
-state-of-the-art higher-dimensional rules are available:
+Besides the standard building blocks above, `modepy` also offers a wide array
+of quadrature rules that can be used on each reference element. The quadrature
+rules are provided as implementations of the `Quadrature` superclass. For the
+interval, a broad selection of quadrature rules are present: Clenshaw--Curtis,
+Fejér, and Jacobi-Gauss(-Lobatto). Additionally, several state-of-the-art
+higher-dimensional rules are available:
 
 * Grundmann--Möller [@Grundmann1978] rules for the $n$-simplex.
 * Vioreanu--Rokhlin [@Vioreanu2014] rules for the two- and
-  three-dimensional simplex (see \autoref{FigureQuadrature}). The quadrature
-  nodes from these rules are also suitable for interpolation.
+  three-dimensional simplex (see \autoref{FigureQuadrature}). 
 * Xiao--Gimbutas [@Xiao2010] rules for the two- and three-dimensional
   simplex.
 * Jáskowiec--Sukumar [@Jaskowiec2021] rules for the tetrahedron.
 * Witherden--Vincent [@Witherden2015] rules for the two- and
-  three-dimensional hypercube (see \autoref{FigureQuadrature}). These rules
-  have significantly fewer nodes than an equivalent tensor product rule, while
-  maintaining positive weights.
+  three-dimensional hypercube (see \autoref{FigureQuadrature}).
 
 ![(left) Vioreanu--Rokhlin quadrature points of order 11 and (right) Witherden--Vincent quadrature points of order 11.](images/quadrature_rules.png){#FigureQuadrature width="50%"}
 
@@ -296,84 +233,11 @@ constructing desired quadrature nodes on a given domain. This can be found in
 While `modepy` is sufficiently flexible to accommodate different applications,
 it also provides functionality specific to FEM needs. In particular, it allows
 constructing a variety of matrix operators that define necessary resampling
-operators and bilinear forms used in FEM codes.
-
-The Vandermonde matrix used in interpolation is provided by the `vandermonde`
-function. This can be used to obtain point-to-point interpolants using
-`resampling_matrix`. In a similar fashion, the `differentiation_matrices`
-function can be used to construct derivative operators along each reference
-axis.
-
-A specific bilinear form matrix can be obtained using
-`nodal_quadrature_bilinear_form_matrix`, which requires providing the usual
-test functions, trial functions and corresponding nodes for the input and
-output. This allows constructing very general families of bilinear forms with
-or without oversampling. One such example is the standard mass matrix, which
-can be obtained more easily using the `mass_matrix` function. If the provided
-functions are derivatives of the basis functions, standard weak forms of
-differentiation operators can also be obtained.
-
-# A Practical Example
-
-In the following example, we demonstrate how to use the library to define a
-non-default tensor product shape --- the prism --- and construct a first
-derivative in the $x$ direction on this reference element. To produce the
-derivative matrix operator, we go through all the objects described in the
-previous section: nodes, modes, and quadrature rules.
-
-```python
-import numpy as np
-import modepy as mp
-
-# Define the shape on which we will operate
-line = mp.Simplex(1)
-triangle = mp.Simplex(2)
-prism = mp.TensorProductShape((triangle, line))
-
-assert prism.dim == 3
-
-# Define a function space for the prism of order 12
-n = 12
-space = mp.TensorProductSpace((mp.PN(triangle.dim, n), mp.PN(line.dim, n)))
-
-assert space.order == n
-assert space.spatial_dim == 3
-assert space.space_dim == (n + 1) * (n + 1) * (n + 2) // 2
-
-# Define a basis function for the prism
-basis = mp.orthonormal_basis_for_space(space, prism)
-
-# Define a point set for the prism
-nodes = mp.edge_clustered_nodes_for_space(space, prism)
-
-# Define a quadrature rule for the prism
-quadrature = mp.TensorProductQuadrature([
-    mp.VioreanuRokhlinSimplexQuadrature(n, triangle.dim),
-    mp.LegendreGaussQuadrature(n),
-])
-
-# Define a bilinear form: weak derivative in the i=x direction
-i = 0
-weak_d = mp.nodal_quadrature_bilinear_form_matrix(
-    quadrature,
-    test_functions=basis.derivatives(i),
-    trial_functions=basis.functions,
-    nodal_interp_functions_test=basis.functions,
-    nodal_interp_functions_trial=basis.functions,
-    input_nodes=nodes,
-    output_nodes=nodes,
-)
-inv_mass = mp.inverse_mass_matrix(basis, nodes)
-
-# Compute derivative and check exact solution
-x = nodes[i]
-f = 1.0 - np.sin(x) ** 3
-f_ref = -3.0 * np.cos(x) * np.sin(x) ** 2
-f_approx = inv_mass @ weak_d.T @ f
-
-error = np.linalg.norm(f_approx - f_ref) / np.linalg.norm(f_ref)
-assert error < 1.0e-4
-```
+operators and bilinear forms used in FEM codes. A specific bilinear form matrix
+can be obtained using `nodal_quadrature_bilinear_form_matrix`, which requires
+providing the usual test functions, trial functions and corresponding nodes for
+the input and output. This allows constructing very general families of
+bilinear forms with or without oversampling. 
 
 # Acknowledgements