FunctionSpace: enable quad_scheme and docs

Author: pbrubeck

docs/source/element_list.py

@@ -1,5 +1,4 @@
 from finat.ufl.elementlist import ufl_elements
-# ~ from ufl.finiteelement.elementlist import ufl_elements
 from finat.element_factory import supported_elements
 import csv
 

docs/source/extruded-meshes.rst

@@ -304,7 +304,7 @@ supports a more general syntax:
 
     V = FunctionSpace(mesh, element)
 
-where ``element`` is a UFL :py:class:`~ufl.finiteelement.finiteelement.FiniteElement` object. This
+where ``element`` is a UFL :py:class:`~finat.ufl.finiteelement.FiniteElement` object. This
 requires generation and manipulation of ``FiniteElement`` objects.
 
 Geometrically, an extruded mesh cell is the *product* of a base, "horizontal",
@@ -319,7 +319,7 @@ The ``TensorProductElement`` operator
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 To create an element compatible with an extruded mesh, one should use
-the :py:class:`~ufl.finiteelement.tensorproductelement.TensorProductElement`
+the :py:class:`~finat.ufl.tensorproductelement.TensorProductElement`
 operator. For example,
 
 .. code-block:: python3
@@ -359,7 +359,7 @@ but is piecewise constant horizontally.
 
 A more complicated element, like a Mini horizontal element with linear
 variation in the vertical direction, may be built using the
-:py:class:`~ufl.finiteelement.enrichedelement.EnrichedElement` functionality
+:py:class:`~finat.ufl.enrichedelement.EnrichedElement` functionality
 in either of the following ways:
 
 .. code-block:: python3
@@ -392,7 +392,7 @@ The ``HDivElement`` and ``HCurlElement`` operators
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 For moderately complicated vector-valued elements,
-:py:class:`~ufl.finiteelement.tensorproductelement.TensorProductElement`
+:py:class:`~finat.ufl.tensorproductelement.TensorProductElement`
 does not give enough information to unambiguously produce the desired
 space. As an example, consider the lowest-order *Raviart-Thomas* element on a
 quadrilateral. The degrees of freedom live on the facets, and consist of
@@ -417,7 +417,7 @@ However, this is only scalar-valued. There are two natural vector-valued
 elements that can be generated from this: one of them preserves tangential
 continuity between elements, and the other preserves normal continuity
 between elements. To obtain the Raviart-Thomas element, we must use the
-:py:class:`~ufl.finiteelement.hdivcurl.HDivElement` operator:
+:py:class:`~finat.ufl.hdivcurl.HDivElement` operator:
 
 .. code-block:: python3
 
@@ -433,7 +433,7 @@ between elements. To obtain the Raviart-Thomas element, we must use the
   :align: center
 
   The RT quadrilateral element, requiring the use
-  of :py:class:`~ufl.finiteelement.hdivcurl.HDivElement`
+  of :py:class:`~finat.ufl.hdivcurl.HDivElement`
 
 Another reason to use these operators is when expanding a vector into a
 higher dimensional space. Consider the lowest-order Nedelec element of the
@@ -453,7 +453,7 @@ the product of this with a continuous element on an interval:
 
 This element still only has two components. To expand this into a
 three-dimensional curl-conforming element, we must use the
-:py:class:`~ufl.finiteelement.hdivcurl.HCurlElement` operator; the syntax is:
+:py:class:`~finat.ufl.hdivcurl.HCurlElement` operator; the syntax is:
 
 .. code-block:: python3
 

docs/source/manual.rst

@@ -19,6 +19,7 @@ Manual
    duals
    interpolation
    point-evaluation
+   quadrature
    external_operators
    adjoint
    visualisation

docs/source/quadrature.rst

@@ -0,0 +1,117 @@
+.. only:: html
+
+  .. contents::
+
+Quadrature rules
+================
+
+To numerically compute the integrals in the variational formulation,
+a quadrature rule is required.
+By default, Firedrake obtains a quadrature rule by estimating the polynomial
+degree of the integrands within a :py:class:`ufl.Form`. Sometimes
+this estimate might be quite large, and a warning like this one will be raised:
+
+.. code-block::
+
+   tsfc:WARNING Estimated quadrature degree 13 more than tenfold greater than any argument/coefficient degree (max 1)
+
+For integrals with very complicated nonlinearities, the estimated quadrature
+degree might be in the hundreds or thousands, rendering the integration
+prohibitively expensive, or leading to segfaults.
+
+Specifying the quadrature rule in the variational formulation
+-------------------------------------------------------------
+
+To manually override the default, the
+quadrature degree can be prescribed on each integral :py:class:`~ufl.measure.Measure`,
+
+.. code-block:: python3
+
+   inner(sin(u)**4, v) * dx(degree=4)
+
+Setting ``degree=4`` means that the quadrature rule will be exact only for integrands
+of total polynomial degree up to 4. This, of course, will introduce a greater numerical error than the default.
+
+For integrals that do not specify a quadrature degree, the default may be keyed as
+``"quadrature_degree"`` in the ``form_compiler_parameters`` dictionary passed on to
+:py:func:`~.solve`, :py:func:`~.project`, or :py:class:`~.NonlinearVariationalProblem`.
+
+.. code-block:: python3
+
+   F = inner(grad(u), grad(v))*dx(degree=0) + inner(exp(u), v)*dx - inner(1, v)*dx
+
+   solve(F == 0, u, form_compiler_parameters={"quadrature_degree": 4})
+
+In the example above, only the integrals with unspecified quadrature degree
+will be computed on a quadrature rule that exactly integrates polynomials of
+the degree set in ``form_compiler_parameters``.
+
+Another way to specify the quadrature rule is through the ``scheme`` keyword. This could be
+either a :py:class:`~finat.quadrature.QuadratureRule`, or a string. Supported string values
+are ``"default"``, ``"canonical"``, and ``"KMV"``. For more details see
+:py:func:`~FIAT.quadrature_schemes.create_quadrature`.  
+
+Lumped quadrature schemes
+-------------------------
+
+Spectral elements, such as Gauss-Legendre-Lobatto and `KMV`_, may be used with
+lumped quadrature schemes to produce a diagonal mass matrix.
+
+.. literalinclude:: ../../tests/firedrake/regression/test_quadrature_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_lump_scheme 1]
+   :end-before: [test_lump_scheme 2]
+
+.. Note::
+
+   To obtain the lumped mass matrix with ``scheme="KMV"``,
+   the ``degree`` argument should match the degree of the :py:class:`~finat.ufl.finiteelement.FiniteElement`.
+
+The Quadrature space
+--------------------
+
+It is possible to define a finite element :py:class:`~.Function` on a quadrature rule.
+The ``"Quadrature"`` and ``"Boundary Quadrature"`` spaces are useful to
+interpolate data at quadrature points on cell interiors and cell boundaries,
+respectively.
+
+.. literalinclude:: ../../tests/firedrake/regression/test_quadrature_manual.py
+   :language: python3
+   :dedent:
+   :start-after: [test_quadrature_space 1]
+   :end-before: [test_quadrature_space 2]
+
+The ``quad_scheme`` keyword argument again may be either
+:py:class:`~finat.quadrature.QuadratureRule` or a string.
+If a :py:class:`~.Function` in the ``"Quadrature"`` space appears within an
+integral, Firedrake will automatically select the quadrature rule that corresponds
+to ``dx(degree=quad_degree, scheme=quad_scheme)`` to match the one associated
+with the quadrature space.
+
+.. _element_quad_scheme:
+
+Specifying the quadrature for integral-type degrees of freedom
+--------------------------------------------------------------
+
+Finite element spaces with :ref:`integral-type degrees of freedom <element_variants>`
+support different quadrature rules.
+These are selected by passing a string to the ``"quad_scheme"`` keyword argument of
+the :py:class:`~finat.ufl.finiteelement.FiniteElement` or
+:py:func:`~.FunctionSpace` constructors. For example, to construct a
+Crouzeix-Raviart space with degrees of freedom consisting of integrals along the edges
+computed from a 2-point average at the endpoints, one can set ``quad_scheme="KMV"``:
+
+.. code-block:: python3
+
+    fe = FiniteElement("Crouzeix-Raviart", triangle, 1, variant="integral", quad_scheme="KMV")
+
+.. Note::
+
+    Finite elements with integral-type degrees of freedom only accept string
+    values for ``quad_scheme``, since it does not make sense to specify a
+    concrete :py:class:`~finat.quadrature.QuadratureRule` when the degrees of
+    freedom are defined on cell entities of different dimensions.
+
+
+.. _KMV : https://defelement.org/elements/kong-mulder-veldhuizen.html

docs/source/variational-problems.rst

@@ -180,7 +180,7 @@ family such as ``"Raviart-Thomas"``. Secondly, you may use the
 family, which gives a vector-valued space where each component is
 identical to the appropriate scalar-valued
 :py:class:`~.FunctionSpace`.  Thirdly, you can create a
-:py:class:`~ufl.classes.VectorElement` directly (which is itself
+:py:class:`~finat.ufl.mixedelement.VectorElement` directly (which is itself
 *vector-valued* and pass that to the :py:func:`~.FunctionSpace`
 constructor).
 
@@ -275,9 +275,11 @@ Firedrake supports the use of the following finite elements.
     :file: element_list.csv
 
 In addition, the
-:py:class:`~ufl.finiteelement.tensorproductelement.TensorProductElement`
+:py:class:`~finat.ufl.tensorproductelement.TensorProductElement`
 operator can be used to create product elements on extruded meshes.
 
+.. _element_variants:
+
 Element variants
 ~~~~~~~~~~~~~~~~
 
@@ -288,15 +290,49 @@ For discontinuous elements these are the Gauss-Legendre points, and for
 continuous elements these are the Gauss-Lobatto-Legendre points.
 For CG and DG spaces on simplices, Firedrake offers both equispaced points and
 the better conditioned recursive Legendre points from :cite:`Isaac2020` via the
-`recursivenodes`_ module. These are selected by passing `variant="equispaced"`
-or `variant="spectral"` to the :py:class:`~ufl.classes.FiniteElement` or
+`recursivenodes`_ module. These are selected by passing ``variant="equispaced"``
+or ``variant="spectral"`` to the :py:class:`~finat.ufl.finiteelement.FiniteElement` or
 :py:func:`~.FunctionSpace` constructors. For example:
 
 .. code-block:: python3
 
     fe = FiniteElement("RTCE", quadrilateral, 2, variant="equispaced")
 
-The default is the spectral variant.
+For CG and DG spaces, and most tensor-product elements, the default is ``variant="spectral"``.
+
+Integral-type degrees of freedom are also available through ``variant="integral"``. These
+enable orthogonal bases for DG on any cell type:
+
+.. code-block:: python3
+
+    V = FunctionSpace(mesh, "DG", 1, variant="integral")
+
+In this case, the degrees of freedom are integral moments against a basis of polynomials,
+which need to be computed with a quadrature rule. The degree of accuracy of the quadrature
+can be increased also through the variant option. For instance, ``variant="integral(4)"``
+will use a quadrature that exactly evaluates the degrees of freedom on polynomials of
+4 degrees higher than that of the space.
+Furthermore, the quadrature rule can be specified via the ``"quad_scheme"`` keyword argument.
+For more details, see the manual section on :ref:`quadrature rules <element_quad_scheme>`.
+
+Integral-type degrees of freedom are useful to nearly preserve curl-free or
+divergence-free functions when interpolating into the finite element space.
+H(curl) and H(div) elements, such as Nédélec, Brezzi-Douglas-Marini, and Raviart-Thomas,
+have ``variant="integral"`` as default. These elements also support ``variant="point"``,
+where the degrees of freedom are vector components evaluated at equispaced points. 
+
+Macroelements can be constructed by splitting each cell into subcells and tiling
+the element on each subcell, without refining the mesh. For example,
+the continuous Lagrange space where each simplex is subdivided
+by connecting each of its vertices to the barycenter is constructed as
+
+.. code-block:: python3
+
+    V = FunctionSpace(mesh, "CG", 1, variant="alfeld")
+
+The space P1-iso-P2 may be similarly constructed with ``variant="iso(2)"``.
+Macroelements also support integral-type degrees of freedom,
+for example ``variant="alfeld,integral(2)"``.
 
 
 Expressing a variational problem
@@ -685,7 +721,7 @@ More complicated forms
 
 UFL is a fully-fledged language for expressing variational problems,
 and hence has operators for all appropriate vector calculus operations
-along with special support for discontinuous galerkin methods in the
+along with special support for discontinuous Galerkin methods in the
 form of symbolic expressions for facet averages and jumps.  For an
 introduction to these concepts we refer the user to the `UFL manual
 <UFL_package_>`_ as well as the :ref:`Firedrake tutorials

firedrake/cofunction.py

@@ -344,7 +344,7 @@ def interpolate(self,
             block on the Pyadjoint tape.
         **kwargs
             Any extra kwargs are passed on to the interpolate function.
-            For details see `firedrake.interpolation.interpolate`.
+            For details see :func:`firedrake.interpolation.interpolate`.
 
         Returns
         -------

firedrake/function.py

@@ -373,7 +373,7 @@ def interpolate(self,
             block on the Pyadjoint tape.
         **kwargs
             Any extra kwargs are passed on to the interpolate function.
-            For details see `firedrake.interpolation.interpolate`.
+            For details see :func:`firedrake.interpolation.interpolate`.
 
         Returns
         -------