# 2. Constructing finite elements¶

##
A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

At the core of the finite element method is the representation of finite-dimensional function spaces over elements. This concept was formalised by [Cia02]:

A *finite element* is a triple \((K, P, N)\) in which \(K\) is a cell,
\(P\) is a space of functions \(K\rightarrow\mathbb{R}^n\) and \(N\), the
set of *nodes*, is a basis for \(P^*\), the dual space to \(P\).

Note that this definition includes a basis for \(P^*\), but not a basis for \(P\). It turns out to be most convenient to specify the set of nodes for an element, and then derive an appropriate basis for \(P\) from that. In particular:

Let \(N = \{\phi^*_j\}\) be a basis for \(P^*\). A *nodal
basis*, \(\{\phi_i\}\) for \(P\) is a basis for \(P\)
with the property that \(\phi^*_j(\phi_i) = \delta_{ij}\).

## 2.1. A worked example¶

##
A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

To illustrate the construction of a nodal basis, let’s consider the linear polynomials on a triangle. We first need to define our reference cell. The obvious choice is the triangle with vertices \(\{(0,0), (1,0), (0,1)\}\)

Functions in this space have the form \(a + bx + cy\). So the function space has three unknown parameters, and its basis (and dual basis) will therefore have three members. In order to ensure the correct continuity between elements, the dual basis we need to use is the evaluation of the function at each of the cell vertices. That is:

We know that \(\phi_i((x,y))\) has the form \(a_i + b_ix + c_iy\) so now we can use the definition of the nodal basis to determine the unknown coefficients:

So for \(\phi_0\) we have:

Which has solution \(\phi_0 = 1 - x - y\). We can write the equations for all the basis functions at once as a single matrix equation:

By which we establish that the full basis is given by:

## 2.2. Types of node¶

We have just encountered nodes given by the evaluation of the function at a given point. Other forms of functional are also suitable for use as finite element nodes. Examples include the integral of the function over the cell or some sub-entity and the evaluation of the gradient of the function at some point. For some vector-valued function spaces, the nodes may be given by the evaluation of the components of the function normal or tangent to the boundary of the cell at some point.

In this course we will only consider point evaluation nodes. The implementation of several other forms of node are covered in [Kir04].

## 2.3. The Lagrange element nodes¶

The number of coefficients of a degree \(p\) polynomial in \(d\) dimensions is given by the combination \(\binom{p+d}{d}\). The simplest set of nodes which we can employ is simply to place these nodes in a regular grid over the reference cell. Given the classical relationship between binomial coefficients and Pascal’s triangle (and between trinomial coefficients and Pascal’s pyramid), it is unsurprising that this produces the correct number of nodes.

The set of equally spaced points of degree \(p\) on the triangle is:

The finite elements with this set of nodes are called the *equispaced
Lagrange* elements and are the most commonly used elements for
relatively low order computations.

While this is the simplest node ordering to construct, when we come to build finite element spaces over a whole computational mesh in Section 4, it will be much more straightforward if the nodes are numbered in topological order. That is to say, the lowest numbered nodes are those associated with the vertices, followed by those associated with the edges and finally, in two dimensions, those associated with the cell. For reasons that will become apparent when we consider the continuity of finite element spaces, the nodes associated with the edges need to be in edge orientation order. That is to say, the node number increases as one moves along the edge in the direction of the arrow. In two dimensions, the ordering of nodes in the cell interior is arbitrary.

Note

At higher order the equispaced Lagrange basis is poorly conditioned and creates unwanted oscillations in the solutions. However for this course Lagrange elements will be sufficient.

Implement
`lagrange_points()`

. Make sure your
algorithm also works for one-dimensional elements. Some basic tests
for your code are to be found in
`test/test_02_lagrange_points.py`

. You can also test your lagrange
points on the triangle by running:

```
plot_lagrange_points degree
```

Where `degree`

is the degree of the points to plot.

## 2.4. Solving for basis functions¶

##
A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

The matrix in (2.3) is a *generalised Vandermonde* [1]
matrix . Given a list of points \((x_i,y_i) \in \mathbb{R}^2, 0\leq i< m\)
the corresponding degree \(n\) generalised Vandermonde matrix is given by:

If we construct the Vandermonde matrix for the nodes of a finite element, then the equation for the complete set of basis function polynomial coefficients is:

where the \(j\)-th column of \(C\) contains the polynomial coefficients of the basis function corresponding to the \(j\)-th node. For (2.8) to be well-posed, there must be a number of nodes equal to the number of coefficients of a degree \(n\) polynomial. If this is the case, then it follows immediately that:

The same process applies to the construction of basis functions for elements in one or three dimensions, except that the Vandermonde matrix must be modified to exclude powers of \(y\) (in one dimension) or to include powers of \(z\).

Note

Here we employ a monomial basis to represent polynomial spaces: any polynomial is given as a linear sum of monomials such as \(x\), \(xy\) or \(x^2\). This basis becomes increasingly ill-conditioned at higher order, so it may be advantageous to employ a different basis in the construction of the Vandermonde matrix. See [Kir04] for an example.

Use (2.7) to implement
`vandermonde_matrix()`

. Think
carefully about how to loop over each row to construct the correct
powers of \(x\) and \(y\). For the purposes of this exercise you should
ignore the `grad`

argument.

Tests for this function are in `test/test_03_vandermonde_matrix.py`

Hint

You can use numpy array operations to construct whole columns of the matrix at once.

## 2.5. Implementing finite elements in Python¶

##
A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

The Ciarlet triple \((K, P, N)\) also provides a
good abstraction for the implementation of software objects
corresponding to finite elements. In our case \(K\) will be a
`ReferenceCell`

. In this course we
will only implement finite element spaces consisting of complete
polynomial spaces so we will specify \(P\) by providing the maximum
degree of the polynomials in the space. Since we will only deal with
point evaluation nodes, we can represent \(N\) by a series of points at
which the evaluation should occur.

Implement the rest of the
`FiniteElement`

`__init__()`

method. You should construct a Vandermonde matrix for the nodes and
invert it to create the basis function coefs. Store these as
`self.basis_coefs`

.

Some basic tests of your implementation are in
`test/test_04_init_finite_element.py`

.

Hint

The `numpy.linalg.inv()`

function may be
used to invert the matrix.

## 2.6. Implementing the Lagrange Elements¶

##
A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

The `FiniteElement`

class implements
a general finite element object assuming we have provided the cell,
polynomial, degree and nodes. The
`LagrangeElement`

class is a
subclass of
`FiniteElement`

which will implement
the particular case of the equispaced Lagrange elements.

Implement the `__init__()`

method of
`LagrangeElement`

. Use
`lagrange_points()`

to obtain the
nodes. For the purpose of this exercise, you may ignore the
`entity_nodes`

argument.

**After** you have implemented
`tabulate()`

in the
next exercise, you can use
`plot_lagrange_basis_functions`

to visualise your
Lagrange basis functions.

## 2.7. Tabulating basis functions¶

##
A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

A core operation in the finite element method is integrating
expressions involving functions in finite element spaces. This is
usually accomplished using numerical quadrature. This means that we need to be able to evaluate the
basis functions at a set of quadrature points. The operation of
evaluating a set of basis functions at a set of points is called
*tabulation*.

Recall that the coefficients of the basis functions are defined with respect to the monomial basis in (2.9). To tabulate the basis functions at a particular set of points therefore requires that the monomial basis be evaluated at that set of points. In other words, the Vandermonde matrix needs to be evaluated at the quadrature points. Suppose we have a set of points \(\{X_i\}\) and a set of basis functions \(\{\phi_j\}\) with coefficents with respect to the monomial basis given by the matrix \(C\). Then the tabulation matrix is given by:

Implement `tabulate()`

.
You can use a Vandermonde matrix to evaluate the polynomial terms
and take the matrix product of this with the basis function
coefficients. The method should have at most two executable
lines. For the purposes of this exercise, ignore the `grad`

argument.

The test file `test/test_05_tabulate.py`

checks that tabulating the
nodes of a finite element produces the identity matrix.

## 2.8. Gradients of basis functions¶

##
A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

A function \(f\) defined over a single finite element with basis \(\{\phi_i\}\) is represented by a weighted sum of that basis:

In order to be able to represent and solve PDEs, we will naturally also have terms incorporating derivatives. Since the coefficients \(f_i\) are spatially constant, derivative operators pass through to apply to the basis functions:

This means that we will need to be able to evaluate the gradient of the basis functions at quadrature points. Recall once again that the basis functions are evaluated by multiplying the Vandermonde matrix evaluated at the relevant points by the matrix of basis function coefficients. Hence:

The last step follows because \(C\) is not a function of \(X\), so it passes through \(\nabla\). The effect of this is that evaluating the gradient of a function in a finite element field just requires the evaluation of the gradient of the Vandermonde matrix.

Extend `vandermonde_matrix()`

so that
setting `grad`

to `True`

produces a rank 3 generalised
Vandermonde tensor whose indices represent points, monomial basis function,
and gradient component respectively. That is:

In other words, each entry of \(V\) is replaced by a vector of the gradient of that polynomial term. For example, the entry \(x^2y^3\) would be replaced by the vector \([ 2xy^3, 3x^2y^2 ]\).

The `test/test_06_vandermonde_matrix_grad.py`

file has tests of this
extension. You should also ensure that you still pass
`test/test_03_vandermonde_matrix.py`

.

Hint

The `transpose()`

method of numpy arrays enables
generalised transposes swapping any dimensions.

Hint

At least one of the natural ways of implementing this function
results in a whole load of `nan`

values in the generalised
Vandermonde matrix. In this case, you might find
`numpy.nan_to_num()`

useful.

Extend `tabulate()`

to
pass the `grad`

argument through to
`vandermonde_matrix()`

. Then
generalise the matrix product in
`tabulate()`

so that
the result of this function (when `grad`

is true) is a rank 3
tensor:

where \(\mathbf{e}_0\ldots\mathbf{e}_{\dim -1}\) is the coordinate basis on the reference cell.

The `test/test_07_tabulate_grad.py`

script tests this
extension. Once again, make sure you still pass
`test/test_05_tabulate.py`

##
A video recording about this exercise is available here.

Imperial students can also watch this video on Panopto

Hint

The `numpy.einsum()`

function implements generalised tensor
contractions using Einstein summation notation. For
example:

```
A = numpy.einsum("ijk,jl->ilk", T, C)
```

is equivalent to \(A_{ilk} = \sum_j T_{ijk} C_{jl}\).

## 2.9. Interpolating functions to the finite element nodes¶

##
A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

Recall once again that a function can be represented on a single finite element as:

Since \(\{\phi_i\}\) is a nodal basis, it follows immediately that:

where \(\phi_i^*\) is the node associated with the basis function \(\phi_i\). Since we are only interested in nodes which are the point evaluation of their function input, we know that:

where \(X_i\) is the point associated with the \(i\)-th node.

Implement `interpolate()`

.

Once you have done this, you can use the script provided to plot functions of your choice interpolated onto any of the finite elements you can make:

```
plot_interpolate_lagrange "sin(2*pi*x[0])" 2 5
```

Hint

You can find help on the arguments to this function with:

```
plot_interpolate_lagrange -h
```

Footnotes