5. Functions in finite element spaces¶

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

(5.1)¶\[f(x) = \sum_i f_i \phi_i(x)\]

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¶

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¶

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:

(5.2)¶\[\sum_i f_i \phi_i(x) \approx g(x)\]

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

(5.3)¶\[f_i = n_i(g(x))\]

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:

(5.4)¶\[f_i = g(x_i)\]

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:

(5.5)¶\[x = \sum_j \hat{x}_j \Psi_j(X) \quad \forall x \in c.\]

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:

(5.6)¶\[A_{i,j} = \Psi_j(X_i)\]

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

(5.7)¶\[x = A\cdot \hat{x}\]

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¶

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¶

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.

Exercise 5.37

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¶

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\):

(5.8)¶\[\int_\Omega f \mathrm{d} x \quad f \in V\]

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:

(5.9)¶\[\int_\Omega f \mathrm{d} x = \sum_{c\in\Omega} \int_c f \mathrm{d} x.\]

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:

(5.10)¶\[\int_{c} f(x) \mathrm{d} x = \int_{c_0} f(X) |J|\mathrm{d} X\]

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

(5.11)¶\[J_{\alpha\beta} = \frac{\partial x_\alpha}{\partial X_\beta}.\]

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:

(5.12)¶\[J_{\alpha\beta} = \sum_j (\tilde{x}_j)_\alpha \nabla_\beta\Psi_j(X)\]

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¶

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

(5.13)¶\[\int_c f(x)\mathrm{d}x = \int_{c_0} \sum_i F(M(c,i))\,\Phi_i(X)\, |J|\,\mathrm{d} X\]

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.14)¶\[\int_c f(x)\mathrm{d}x = \sum_q \sum_i F(M(c,i))\,\Phi_i(X_q)\, |J|\,w_q\]

5.3.3. Implementing integration¶

Exercise 5.38

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.

Exercise 5.39

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.