# 5. Functions in finite element spaces¶

Recall that the general form of a function in a finite element space is:

Where the \(\phi_i(x)\) are now the global basis functions achieved by stitching together the local basis functions defined by the finite element.

## 5.1. A python implementation of functions in finite element spaces¶

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

Imperial students can also watch this video on Panopto

The `Function`

class provides a
simple implementation of function storage. The input is a
`FunctionSpace`

which defines the
mesh and finite element to be employed, to which the
`Function`

adds an array of degree of
freedom values, one for each node in the
`FunctionSpace`

.

## 5.2. Interpolating values into finite element spaces¶

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

Imperial students can also watch this video on Panopto

Suppose we have a function \(g(x): \mathbb{R}^n \rightarrow \mathbb{R}\) which we wish to approximate as a function \(f(x)\) in some finite element space \(V\). In other words, we want to find the \(f_i\) such that:

The simplest way to do this is to *interpolate* \(g(x)\) onto \(V\). In
other words, we evaluate:

where \(n_i\) is the node associated with \(\phi_i\) as considered over the
entire mesh (*globally*). Since we are only concerned with point
evaluation nodes, this is equivalent to:

where \(x_i\) is the coordinate vector of the point defining the node
\(n_i\). This looks straightforward, however the \(x_i\) are the *global*
node points, and so far we have only defined the node points in
*local* coordinates on the reference element.

### 5.2.1. Changing coordinates between reference and physical space¶

We’ll refer to coordinates on the global mesh as being in *physical
space* while those on the reference element are in *local
space*. We’ll use case to distinguish local and global objects, so
local coordinates will be written as \(X\) and global coordinates as
\(x\). The key observation is that within each cell, the global
coordinates are the linear interpolation of the global coordinate
values at the cell vertices. In other words, if \(\{\Psi_j\}\) is the
local basis for the **linear** lagrange elements on the reference cell and
\(\hat{x}_j\) are the corresponding global vertex locations on a cell \(c\)
then:

Remember that we know the location of the nodes in local coordinates,
and we have the
`tabulate()`

method to
evaluate all the basis functions of an element at a known set of
points. So if we write:

where \(\{X_i\}\) are the local node points of our finite element, then the corresponding global node points \(\{x_i\}\) are given by

where \(\hat{x}\) is the \((\dim+1, \dim)\) array whose rows are the current
element vertex coordinates. Specifically, \(x\) is the \((\textrm{nodes}, \dim)\) array whose
rows are the global coordinates of the nodes in the current
element. We can then apply \(g()\) to each row of \(x\) in turn and record
the result as the `Function`

value
for that node.

Hint

The observant reader will notice that this algorithm is inefficient because the function values at nodes on the boundaries of elements are evaluated more than once. This can be avoided with a little tedious bookkeeping but we will not concern ourselves with that here.

### 5.2.2. Looking up cell coordinates and values¶

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

Imperial students can also watch this video on Panopto

In the previous section we used the vertex coordinates of a cell to
find the node coordinates, and then we calculated
`Function`

values at those
points. The coordinates are stored in a single long list associated
with the `Mesh`

, and the
`Function`

contains a single long
list of values. We need to use *indirect addressing* to access these
values. This is best illustrated using some Python code.

Suppose `f`

is a `Function`

.
For brevity, we write `fs = f.function_space`

, the
`FunctionSpace`

associated with
`f`

. Now, we first need a linear element and a corresponding
`FunctionSpace`

:

```
cg1 = fe_utils.LagrangeElement(fs.mesh.cell, 1)
cg1fs = fe_utils.FunctionSpace(fs.mesh, cg1)
```

Then the vertex indices of cell number `c`

in the correct order for the linear Lagrange element are:

```
cg1fs.cell_nodes[c, :]
```

and therefore the set of coordinate vectors for the vertices of
element `c`

are:

```
fs.mesh.vertex_coords[cg1fs.cell_nodes[c, :], :]
```

That is, the `cg1fs.cell_nodes`

array is used to look up the right
vertex coordinates. By a similar process we can access the values
associated with the nodes of element `c`

:

```
f.values[fs.cell_nodes[c, :]]
```

### 5.2.3. A Python implementation of interpolation¶

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

Imperial students can also watch this video on Panopto

Putting together the change of coordinates with the right indirect
addressing, we can provide the
`Function`

class with a
`interpolate()`

method which
interpolates a user-provided function onto the
`Function`

.

Read and understand the
`interpolate()`

method. Use
`plot_sin_function`

to investigate interpolating different
functions onto finite element spaces at differering resolutions and
polynomial degrees.

Hint

There is no implementation work associated with this exercise, but
the programming constructs used in
`interpolate()`

will be
needed when you implement integration.

## 5.3. Integration¶

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

Imperial students can also watch this video on Panopto

We now come to one of the fundamental operations in the finite element
method: integrating a `Function`

over
the domain. The full finite element method actually requires the
integration of expressions of unknown test and trial functions, but we
will start with the more straightforward case of integrating a single,
known, `Function`

over a domain
\(\Omega\):

where \(\mathrm{d}x\) should be understood as being the volume measure with the correct dimension for the domain and \(V\) is some finite element space over \(\Omega\). We can express this integral as a sum of integrals over individual cells:

So we have in fact reduced the integration problem to the problem of integrating \(f\) over each cell. In a previous part of the module we implemented quadrature rules which enable us to integrate over specified reference cells. If we can express the integral over some arbitrary cell \(c\) as an integral over a reference cell \(c_0\) then we are done. In fact this simply requires us to employ the change of variables formula for integration:

where \(|J|\) is the absolute value of the determinant of the Jacobian matrix. \(J\) is given by:

Hint

We will generally adopt the convention of using Greek letters to indicate indices in spatial dimensions, while we will use Roman letters in the sequence \(i,j,\ldots\) for basis function indices. We will continue to use \(q\) for the index over the quadrature points.

Evaluating (5.11) depends on having an expression for \(x\) in terms of \(X\). Fortunately, (5.5) is exactly this expression, and applying the usual rule for differentiating functions in finite element spaces produces:

where \(\{\Psi_j\}\) is once again the degree 1 Lagrange basis and \(\{\tilde{x}_j\}\) are the coordinates of the corresponding vertices of cell \(c\). The presence of \(X\) in (5.12) implies that the Jacobian varies spatially across the reference cell. However since \(\{\Psi_j\}\) is the degree 1 Lagrange basis, the gradients of the basis functions are constant over the cell and so it does not matter at which point in the cell the Jacobian is evaluated. For example we might choose to evaluate the Jacobian at the cell origin \(X=0\).

Hint

When using simplices with curved sides, and on all but the simplest quadrilateral or hexahedral meshes, the change of coordinates will not be affine. In that case, to preserve full accuracy it will be necessary to compute the Jacobian at every quadrature point. However, non-affine coordinate transforms are beyond the scope of this course.

### 5.3.1. Expressing the function in the finite element basis¶

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

Imperial students can also watch this video on Panopto

Let \(\{\Phi_i(X)\}\) be a **local** basis for \(V\) on the reference element
\(c_0\). Then our integral becomes:

where \(F\) is the vector of global coefficient values of \(f\), and \(M\) is the cell node map.

### 5.3.2. Numerical quadrature¶

The actual evaluation of the integral will employ the quadrature rules we discussed in a previous section. Let \(\{X_q\}, \{w_q\}\) be a quadrature rule of sufficient degree of precision that the quadrature is exact. Then:

### 5.3.3. Implementing integration¶

Use (5.12) to implement the
`jacobian()`

method of
`Mesh`

. `test/test_09_jacobian.py`

is
available for you to test your results.

Hint

The \(\nabla_\beta\Psi_j(X)\) factor in (5.12) is the same
for every cell in the mesh. You could make your implementation more
efficient by precalculating this term in the `__init__()`

method of `Mesh`

.

Use (5.9) and (5.14) to implement
`integrate()`

.
`test/test_10_integrate_function.py`

may be used to test your
implementation.

Hint

Your method will need to:

Construct a suitable

`QuadratureRule`

.`tabulate()`

the basis functions at each quadrature point.Visit each cell in turn.

Construct the

`jacobian()`

for that cell and take the absolute value of its determinant (`numpy.absolute`

and`numpy.linalg.det()`

will be useful here).Sum all of the arrays you have constructed over the correct indices to a contribution to the integral (

`numpy.einsum()`

may be useful for this).

Hint

You might choose to read ahead before implementing
`integrate()`

, since the
`errornorm()`

function is very similar and may provide a useful
template for your work.