diff --git a/docs/source/element_list.py b/docs/source/element_list.py index 6e8c6ab981..4671bfcf1b 100644 --- a/docs/source/element_list.py +++ b/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 diff --git a/docs/source/extruded-meshes.rst b/docs/source/extruded-meshes.rst index 363df0a6cc..b37532577c 100644 --- a/docs/source/extruded-meshes.rst +++ b/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 diff --git a/docs/source/manual.rst b/docs/source/manual.rst index 4cbf29c126..d7e7e655ff 100644 --- a/docs/source/manual.rst +++ b/docs/source/manual.rst @@ -19,6 +19,7 @@ Manual duals interpolation point-evaluation + quadrature external_operators adjoint visualisation diff --git a/docs/source/quadrature.rst b/docs/source/quadrature.rst new file mode 100644 index 0000000000..02f0bb7ec6 --- /dev/null +++ b/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 ` +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 diff --git a/docs/source/variational-problems.rst b/docs/source/variational-problems.rst index 8c38e7271d..e164f6c620 100644 --- a/docs/source/variational-problems.rst +++ b/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 `. + +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 `_ as well as the :ref:`Firedrake tutorials diff --git a/firedrake/cofunction.py b/firedrake/cofunction.py index b975ca56d1..b08743c4a6 100644 --- a/firedrake/cofunction.py +++ b/firedrake/cofunction.py @@ -324,7 +324,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 ------- diff --git a/firedrake/function.py b/firedrake/function.py index a6b658f356..d37e62b91e 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -375,7 +375,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 ------- diff --git a/firedrake/functionspace.py b/firedrake/functionspace.py index cefa3bf9a4..617e7f6780 100644 --- a/firedrake/functionspace.py +++ b/firedrake/functionspace.py @@ -20,7 +20,7 @@ @PETSc.Log.EventDecorator() -def make_scalar_element(mesh, family, degree, vfamily, vdegree, variant): +def make_scalar_element(mesh, family, degree, vfamily, vdegree, variant, quad_scheme): """Build a scalar :class:`finat.ufl.finiteelement.FiniteElement`. Parameters @@ -31,14 +31,16 @@ def make_scalar_element(mesh, family, degree, vfamily, vdegree, variant): The finite element family. degree : The degree of the finite element. - variant : - The variant of the finite element. vfamily : The finite element in the vertical dimension (extruded meshes only). vdegree : The degree of the element in the vertical dimension (extruded meshes only). + variant : + The variant of the finite element. + quad_scheme : + The quadrature scheme used to evaluate integral-type degrees of freedom. Notes ----- @@ -59,20 +61,24 @@ def make_scalar_element(mesh, family, degree, vfamily, vdegree, variant): and vfamily is not None and vdegree is not None: la = finat.ufl.FiniteElement(family, cell=cell.sub_cells()[0], - degree=degree, variant=variant) + degree=degree, + variant=variant, + quad_scheme=quad_scheme) # If second element was passed in, use it lb = finat.ufl.FiniteElement(vfamily, cell=ufl.interval, - degree=vdegree, variant=variant) + degree=vdegree, + variant=variant, + quad_scheme=quad_scheme) # Now make the TensorProductElement return finat.ufl.TensorProductElement(la, lb) else: - return finat.ufl.FiniteElement(family, cell=cell, degree=degree, variant=variant) + return finat.ufl.FiniteElement(family, cell=cell, degree=degree, variant=variant, quad_scheme=quad_scheme) @PETSc.Log.EventDecorator("CreateFunctionSpace") def FunctionSpace(mesh, family, degree=None, name=None, - vfamily=None, vdegree=None, variant=None): + vfamily=None, vdegree=None, variant=None, quad_scheme=None): """Create a :class:`.FunctionSpace`. Parameters @@ -93,6 +99,10 @@ def FunctionSpace(mesh, family, degree=None, name=None, meshes only). variant : The variant of the finite element. + For more details see the :ref:`manual section on element variants `. + quad_scheme : + The quadrature scheme used to evaluate integral-type degrees of freedom. + For more details see the :ref:`manual section on quadrature schemes `. Notes ----- @@ -101,13 +111,13 @@ def FunctionSpace(mesh, family, degree=None, name=None, are ignored and the appropriate :class:`.FunctionSpace` is returned. """ - element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant) + element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant, quad_scheme) return impl.WithGeometry.make_function_space(mesh, element, name=name) @PETSc.Log.EventDecorator() def DualSpace(mesh, family, degree=None, name=None, - vfamily=None, vdegree=None, variant=None): + vfamily=None, vdegree=None, variant=None, quad_scheme=None): """Create a :class:`.FunctionSpace`. Parameters @@ -128,6 +138,10 @@ def DualSpace(mesh, family, degree=None, name=None, meshes only). variant : The variant of the finite element. + For more details see the :ref:`manual section on element variants `. + quad_scheme : + The quadrature scheme used to evaluate integral-type degrees of freedom. + For more details see the :ref:`manual section on quadrature schemes `. Notes ----- @@ -136,13 +150,13 @@ def DualSpace(mesh, family, degree=None, name=None, other arguments are ignored and the appropriate :class:`.FunctionSpace` is returned. """ - element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant) + element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant, quad_scheme) return impl.FiredrakeDualSpace.make_function_space(mesh, element, name=name) @PETSc.Log.EventDecorator() -def VectorFunctionSpace(mesh, family, degree=None, dim=None, - name=None, vfamily=None, vdegree=None, variant=None): +def VectorFunctionSpace(mesh, family, degree=None, dim=None, name=None, + vfamily=None, vdegree=None, variant=None, quad_scheme=None): """Create a rank-1 :class:`.FunctionSpace`. Parameters @@ -166,6 +180,10 @@ def VectorFunctionSpace(mesh, family, degree=None, dim=None, meshes only). variant : The variant of the finite element. + For more details see the :ref:`manual section on element variants `. + quad_scheme : + The quadrature scheme used to evaluate integral-type degrees of freedom. + For more details see the :ref:`manual section on quadrature schemes `. Notes ----- @@ -178,10 +196,10 @@ def VectorFunctionSpace(mesh, family, degree=None, dim=None, pass it to :class:`.FunctionSpace` directly instead. """ - sub_element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant) + sub_element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant, quad_scheme) if dim is None: dim = mesh.geometric_dimension() - if not isinstance(dim, numbers.Integral) and dim > 0: + if not (isinstance(dim, numbers.Integral) and dim > 0): raise ValueError(f"Can't make VectorFunctionSpace with dim={dim}") element = finat.ufl.VectorElement(sub_element, dim=dim) return FunctionSpace(mesh, element, name=name) @@ -190,7 +208,7 @@ def VectorFunctionSpace(mesh, family, degree=None, dim=None, @PETSc.Log.EventDecorator() def TensorFunctionSpace(mesh, family, degree=None, shape=None, symmetry=None, name=None, vfamily=None, - vdegree=None, variant=None): + vdegree=None, variant=None, quad_scheme=None): """Create a rank-2 FunctionSpace. Parameters @@ -217,6 +235,10 @@ def TensorFunctionSpace(mesh, family, degree=None, shape=None, meshes only). variant : The variant of the finite element. + For more details see the :ref:`manual section on element variants `. + quad_scheme : + The quadrature scheme used to evaluate integral-type degrees of freedom. + For more details see the :ref:`manual section on quadrature schemes `. Notes ----- @@ -230,7 +252,7 @@ def TensorFunctionSpace(mesh, family, degree=None, shape=None, `FunctionSpace` directly instead. """ - sub_element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant) + sub_element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant, quad_scheme) if shape is None: shape = (mesh.geometric_dimension(),) * 2 element = finat.ufl.TensorElement(sub_element, shape=shape, symmetry=symmetry) diff --git a/tests/firedrake/regression/test_quadrature_manual.py b/tests/firedrake/regression/test_quadrature_manual.py new file mode 100644 index 0000000000..6f340f274a --- /dev/null +++ b/tests/firedrake/regression/test_quadrature_manual.py @@ -0,0 +1,41 @@ +from firedrake import * +import pytest +import numpy as np + + +@pytest.fixture +def mesh(): + return UnitSquareMesh(2, 2) + + +def test_lump_scheme(mesh): + # [test_lump_scheme 1] + degree = 3 + V = FunctionSpace(mesh, "KMV", degree) + u = TrialFunction(V) + v = TestFunction(V) + + a = inner(u, v) * dx(scheme="KMV", degree=degree) + # [test_lump_scheme 2] + + A = assemble(a) + d = assemble(inner(1, v)*dx) + with d.dat.vec as dvec: + Dmat = PETSc.Mat().createDiagonal(dvec) + + Bmat = A.petscmat.copy() + Bmat.axpy(-1, Dmat) + assert np.allclose(Bmat.norm(PETSc.NormType.FROBENIUS), 0) + + +@pytest.mark.parametrize("quad_scheme", (None, "default", "canonical")) +def test_quadrature_space(mesh, quad_scheme): + quad_degree = 4 + # [test_quadrature_space 1] + Q = FunctionSpace(mesh, "Quadrature", degree=quad_degree, quad_scheme=quad_scheme) + # [test_quadrature_space 2] + + f = sum(SpatialCoordinate(mesh)) ** 2 + q = Function(Q).interpolate(f) + + assert norm(q - f) < 1E-12