1. Numerical quadrature¶
A video recording of the following material is available here.
Imperial students can also watch this video on Panopto
The core computational operation with which we are concerned in the finite element method is the integration of a function over a known reference element. It’s no big surprise, therefore, that this operation will be at the heart of our finite element implementation.
The usual way to efficiently evaluate arbitrary integrals numerically is numerical quadrature. This basic idea will already be familiar to you from undergraduate maths (or maybe even high school calculus) as it’s the generalisation of the trapezoidal rule and Simpson’s rule for integration.
The core idea of quadrature is that the integral of a function \(f(X)\) over an element \(e\) can be approximated as a weighted sum of function values evaluated at particular points:
we term the set \(\{X_q\}\) the set of quadrature points and the corresponding set \(\{w_q\}\) the set of quadrature weights. A set of quadrature points and their corresponding quadrature weights together comprise a quadrature rule for \(e\). For an arbitrary function \(f\), quadrature is only an approximation to the integral. The global truncation error in this approximation is invariably of the form \(O(h^n)\) where \(h\) is the diameter of the element.
If \(f\) is a polynomial in \(X\) with degree \(p\) such that \(p\leq n2\) then it is easy to show that integration using a quadrature rule of degree \(n\) results in exactly zero error.
The degree of precision of a quadrature rule is the largest \(p\) such that the quadrature rule integrates all polynomials of degree \(p\) without error.
1.1. Exact and incomplete quadrature¶
In the finite element method, integrands are very frequently polynomial. If the quadrature rule employed for a particular interval has a sufficiently high degree of precision such that there is no quadrature error in the integration, we refer to the quadrature as exact or complete. In any other case we refer to the quadrature as incomplete.
Typically, higher degree quadrature rules have more quadrature points than lower degree rules. This results in a tradeoff between the accuracy of the quadrature rule and the number of function evaluations, and hence the computational cost, of an integration using that rule. Complete quadrature results in lower errors, but if the error due to incomplete quadrature is small compared with other errors in the simulation, particularly compared with the discretisation error, then incomplete quadrature may be advantageous.
1.2. Examples in one dimension¶
A video recording of the following material is available here.
Imperial students can also watch this video on Panopto
We noted above that a few one dimensional quadrature rules are commonly taught in introductory integration courses. The first of these is the midpoint rule:
In other words, an approximation to the integral of \(f\) over an interval can be calculated by multiplying the value of \(f\) at the midpoint of the interval by the length of the interval. This amounts to approximating the function over the integral by a constant value.
If we improve our approximation of \(f\) to a straight line over the interval, then we arrive at the trapezoidal (or trapezium) rule:
while if we employ a quadratic function then we arrive at Simpson’s rule:
1.3. Reference cells¶
A video recording of the following material is available here.
Imperial students can also watch this video on Panopto
As a practical matter, we wish to write down quadrature rules as arrays of numbers, independent of \(h\). In order to achieve this, we will write the quadrature rules for a single, reference cell. When we wish to actually integrate a function over cell, we will change coordinates to the reference cell. We will return to the mechanics of this process later, but for now it means that we need only consider quadrature rules on the reference cells we choose.
A commonly employed one dimensional reference cell is the unit interval \([0,1]\), and that is the one we shall adopt here (the other popular alternative is the interval \([1, 1]\), which some prefer due to its symmetry about the origin).
In two dimensions, the cells employed most commonly are triangles and quadrilaterals. For simplicity, in this course we will only consider implementing the finite element method on triangles. The choice of a reference interval implies a natural choice of reference triangle. For the unit interval the natural correspondence is with the triangle with vertices \([(0,0), (1,0), (0,1)]\), though different choices of vertex numbering are possible.
1.3.1. Reference cell topology¶
A cell is composed of topological entities, that is to say vertices, edges, faces and so forth. The topology of the cell is given by the connectivity of its entities, for example which vertices make up each edge. It is useful to define some terms to describe the cell topology:
The dimension of a cell is the maximal dimension of the topological entities that make up the cell.
A topological entity of codimension \(n\) is a topological entity of dimension \(dn\) where \(d\) is the dimension of the cell.
Armed with these definitions we are able to define names for topological entities of various dimension and codimension:
entity name 
dimension 
codimension 

vertex 
0 

edge 
1 

face 
2 

facet 
1 

cell 
0 
The cells can be polygons or polyhedra of any shape, however in this course we will restrict ourselves to intervals and triangles. The only other twodimensional cells frequently employed are quadrilaterals.
1.3.2. Reference cell entities¶
The topological entities of each dimension in a cell are distinguished by giving them unique numbers. We will identify topological entities by an index pair \((d, i)\) where \(i\) is the index of the entity within the set of \(d\)dimensional entities.
The particular choices of numbering we will use are shown in Fig. 1.5. The numbering is a matter of convention: that adopted here is that edges share the number of the opposite vertex. The orientation of the edges is also shown, this is always from the lower numbered vertex to the higher numbered one.
The ReferenceCell
class stores the
local topology of the reference cell. Read the source and ensure that you
understand the way in which this information is encoded.
The following animation of the numbering of the topological entities on the reference cell may help in understanding this.
1.3.3. Python implementations of reference elements¶
The ReferenceCell
class provides
Python objects encoding the geometry and topology of the reference
cell. At this stage, the relevant information is the dimension of the
reference cell and the list of vertices. The topology will become
important in the following chapters. The reference cells we will
require for this course are the
ReferenceInterval
and
ReferenceTriangle
.
1.4. Quadrature rules on reference elements¶
A video recording of the following material is available here.
Imperial students can also watch this video on Panopto
Having adopted a convention for the reference element, we can simply express quadrature rules as lists of quadrature points with corresponding quadrature weights. For example Simpson’s rule becomes:
We choose to write the quadrature points as 1tuples for consistency with the \(n\)dimensional case, in which the points will be \(n\)tuples.
The lowest order quadrature rule on the reference triangle is a single point quadrature:
This rule has a degree of precision of 1.
Hint
The weights of a quadrature rule always sum to the volume of the reference element. Why is this?
1.5. LegendreGauß quadrature in one dimension¶
The finite element method will result in integrands of different polynomial degrees, so it is convenient if we have access to quadrature rules of arbitrary degree on demand. In one dimension the LegendreGauß quadrature rules are a family of rules of arbitrary precision which we can employ for this purpose. Helpfully, numpy provides an implementation which we are able to adopt. The LegendreGauß quadrature rules are usually defined for the interval \([1, 1]\) so we need to change coordinates in order to arrive at a quadrature rule for our reference interval:
where \((\{X'_q\}, \{w'_q\})\) is the quadrature rule on the interval \([1, 1]\) and \((\{X_q\}, \{w_q\})\) is the rule on the unit interval.
LegendreGauß quadrature on the interval is optimal in the sense that it uses the minimum possible number of points for each degree of precision.
1.6. Extending LegendreGauß quadrature to two dimensions¶
We can form a unit square by taking the Cartesian product of two unit intervals: \((0, 1)\otimes (0, 1)\). Similarly, we can form a quadrature rule on a unit square by taking the product of two interval quadrature rules:
where \((X, w)\) is an interval quadrature rule. Furthermore, the degree of accuracy of \((X_\textrm{sq}, w_\textrm{sq})\) will be the same as that of the onedimensional rule.
However, we need a quadrature rule for the unit triangle. We can achieve this by treating the triangle as a square with a zero length edge. The Duffy transform maps the unit square to the unit triangle:
By composing the Duffy transform with (1.11) we can arrive at a quadrature rule for the triangle:
where \((X_v, w_v)\) is a reference interval quadrature rule with degree of precision \(n\) and \((X_h, w_h)\) is a reference interval quadrature rule with degree of precision \(n+1\). The combined quadrature rule \((X_\textrm{tri}, w_\textrm{tri})\) will then be \(n\). The additional degree of precision required for \((X_h, w_h)\) is because the Duffy transform effectively increases the polynomial degree of the integrand by one.
1.7. Implementing quadrature rules in Python¶
A video recording about how to do this exercise is available here.
Imperial students can also watch this video on Panopto
The fe_utils.quadrature
module provides the
QuadratureRule
class which records
quadrature points and weights for a given
ReferenceCell
. The
gauss_quadrature()
function creates
quadrature rules for a prescribed degree of precision and reference
cell.
The integrate()
method is
left unimplemented. Using (1.4), implement this method.
A test script for your method is provided in the test
directory
as test_01_integrate.py
. Run this script to test your code:
py.test test/test_01_integrate.py
from the Bash command line. Make sure you commit your modifications and push them to your fork of the course repository.
Hint
You can implement
integrate()
in one line
using a list
comprehension and numpy.dot()
.
Hint
Don’t forget to activate your Python venv!
A video recording of the solution to this exercise is available here.
Imperial students can also watch this video on Panopto