4. Function spaces: associating data with meshes

A finite element space over a mesh is constructed by associating a finite element with each cell of the mesh. We will refer to the basis functions of this finite element space as global basis functions, while those of the finite element itself we will refer to as local basis functions. We can establish the relationship between the finite element and each cell of the mesh by associating the nodes (and therefore the local basis functions) of the finite element with the topological entities of the mesh. This is a two stage process. First, we associate the nodes of the finite element with the local topological entities of the reference cell. This is often referred to as local numbering. Then we associate the correct number of degrees of freedom (i.e. number of basis functions) with each global mesh entity. This is the global numbering.

4.1. Local numbering and continuity

A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

Which nodes should be associated with which topological entities? The answer to this question depends on the degree of continuity required between adjacent cells. The nodes associated with topological entites on the boundaries of cells (the vertices in one dimension, the vertices and edges in two dimensions, and the vertices, edges and faces in three dimensions) are shared between cells. The basis functions associated with nodes on the cell boundary will therefore be continuous between the cells which share that boundary.

For the Lagrange element family, we require global \(C_0\) continuity. This implies that the basis functions are continuous everywhere. This has the following implications for the association of basis functions with local topological entites:

vertices

At the function vertices we can achieve continuity by requiring that there be a node associated with each mesh vertex. The basis function associated with that node will therefore be continuous. Since we have a nodal basis, all the other basis functions will vanish at the vertex so the global space will be continuous at this point.

edges

Where the finite element space has at least 2 dimensions we need to ensure continuity along edges. The restriction of a degree \(p\) polynomial over a \(d\)-dimensional cell to an edge of that cell will be a one dimensional degree \(p\) polynomial. To fully specify this polynomial along an edge requires \(p+1\) nodes. However there will already be two nodes associated with the vertices of the edge, so \(p-1\) additional nodes will be associated with the edge.

faces

For three-dimensional (tetrahedral) elements, the basis functions must also be continuous across faces. This requires that sufficient nodes lie on the face to fully specify a two dimensional degree \(p\) polynomial. However the vertices and edges of the face already have nodes associated with them, so the number of nodes required to be associated with the face itself is actually the number required to represent a degree \(p-2\) polynomial in two dimensions given by the combination \(\binom{p-1}{2}\).

This pattern holds more generally: for a \(C_0\) function space, the number of nodes which must be associated with a local topological entity of dimension \(d\) is \(\binom{p-1}{d}\).

Fig. 4.1 illustrates the association of nodes with reference entities for Lagrange elements on triangles.

Fig. 4.1 Association of nodes with reference entities for the degree 1, 2, and 3 equispaced Lagrange elements on triangles. Black nodes are associated with vertices, red nodes with edges and blue nodes with the cell (face). This is the same figure as Fig. 2.1.

4.2. Implementing local numbering

A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

Local numbering can be implemented by adding an additional data structure to the FiniteElement class. For each local entity this must record the local nodes associated with that entity. This can be achieved using a dictionary of dictionaries structure. For example employing the local numbering of nodes employed in Fig. 2.1, the entity_node dictionary for the degree three equispaced Lagrange element on a triangle is given by:

entity_node = {0: {0: [0],
                   1: [1],
                   2: [2]},
               1: {0: [3, 4],
                   1: [5, 6],
                   2: [7, 8]},
               2: {0: [9]}}

Note that the order of the nodes in each list is important: it must always consistently reflect the orientation of the relevant entity in order that all the cells which share that entity consistently interpret the nodes. In this case this has been achieved by listing the nodes in order given by the direction of the orientation of each edge.

Exercise 4.44

Extend the __init__() method of LagrangeElement so that it passes the correct entity_node dictionary to the FiniteElement it creates.

The test/test_08_entity_nodes.py script tests this functionality.

4.3. Global numbering

A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

Given a mesh and a finite element, the global numbering task is to uniquely associate the appropriate number of global node numbers with each global entity. One such numbering [1] is to allocate global numbers in ascending entity dimension order, and within each dimension in order of the index of each global topological entity. The formula for the first global node associated with entity \((d, i)\) is then:

(4.2)\[G(d, i) = \left(\sum_{\delta < d} N_\delta E_\delta\right) + iN_d\]

where \(N_d\) is the number of nodes which this finite element associates with a single entity of dimension \(d\), and \(E_d\) is the number of dimension \(d\) entities in the mesh. The full list of nodes associated with entity \((d, i)\) is therefore:

(4.3)\[[G(d, i), \ldots, G(d,i) + N_d - 1]\]

4.4. The cell-node map

A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

The primary use to which we wish to put the finite element spaces we are constructing is, naturally, the solution of finite element problems. The principle operation we will therefore need to support is integration over the mesh of mathematical expressions involving functions in finite element spaces. This will be accomplished by integrating over each cell in turn, and then summing over all cells. This means that a key operation we will need is to find the nodes associated with a given cell.

It is usual in finite element software to explicitly store the map from cells to adjacent nodes as a two-dimensional array with one row corresponding to each cell, and with columns corresponding to the local node numbers. The entries in this map will have the following values:

(4.4)\[M[c, e(\delta, \epsilon)] = [G(\delta, i), \ldots, G(\delta,i) + N_\delta - 1] \qquad\forall 0\leq\delta\leq\dim(c), \forall 0\leq\epsilon < \hat{E}_\delta\]

where:

(4.5)\[i = \operatorname{Adj}_{\dim(c), \delta}[c, \epsilon],\]

\(e(\delta, \epsilon)\) is the local entity-node list for this finite element for the \((\delta, \epsilon)\) local entity, \(\operatorname{Adj}\) has the meaning given under Python implementations of reference elements, \(\hat{E}_\delta\) is the number of dimension \(\delta\) entities in each cell, and \(G\) and \(N\) have the meanings given above. This algorithm requires a trivial extension to adjacency:

(4.6)\[\operatorname{Adj}_{\dim(c),\dim(c)}[c, 0] = c\]

Hint

In (4.4), notice that for each value of \(\delta\) and \(\epsilon\), \(e(\delta, \epsilon)\) is a vector of indices, so the equation sets the value of zero, one, or more defined entries in row \(c\) of \(M\) for each \(\delta\) and \(\epsilon\).

4.5. Implementing function spaces in Python

A video recording of the following material is available here.

Imperial students can also watch this video on Panopto

As noted above, a finite element space associates a mesh and a finite element, and contains (in some form) a global numbering of the nodes.

Exercise 4.45

Implement the __init__() method of fe_utils.function_spaces.FunctionSpace. The key operation is to set cell_nodes using (4.4).

You can plot the numbering you have created with the plot_function_space_nodes script. As usual, run the script passing the -h option to discover the required arguments.

Hint

Many of the terms in (4.4) are implemented in the objects in fe_utils. For example:

• \(\operatorname{Adj}_{\dim(c), \delta}\) is implemented by the adjacency() method of the Mesh.

• You have \(e(\delta, \epsilon)\) as entity_nodes. Note that in this case you need separate square brackets for each index:

  element.entity_nodes[delta][epsilon]
  

Hint

cell_nodes needs to be integer-valued. If you choose to use numpy.zeros() to create a matrix which you then populate with values, you need to explicitly specify that you want a matrix of integers. This can be achieved by passing the dtype argument to numpy.zeros(). For example numpy.zeros((nrows, ncols), dtype=int).

Footnotes

[1]

Many correct global numberings are possible, that presented here is simple and correct, but not optimal from the perspective of the memory layout of the resulting data.