Lecture 6 (c)

Finite elements in one dimension

Unstructured grids

• Most problems have complex domains.

• Boundary conditions crucial.

• Little “structure” to problem.

• What use is “structure” in the grid?

• Use unstructured grids.

• Use Finite Elements.

Most “real” applications happen on complex domains, which require complex grids.

Function bases

Use shape or indicator functions

\[ N_A(x) = \begin{cases} 1 & x = x_A \\ 0 & x = x_j, \quad j \ne A \end{cases} \]

as basis,

\[ \phi = \sum_A c_A(t) N_A(x) \, . \]

Many possibilities; linear case shown.

Here the shape functions are acting as the basis functions. Using indicator functions means that at the grid point, or grid node, \(A\), where \(x = x_A\), the function value is exactly the coefficient. So we immediately have \(c_A = \phi(x_A)\). This works not just for the solution but for any function that we want to represent on the unstructured grid.

This choice of basis (or shape, or indicator) functions is not unique: the functions can do whatever they like between nodes as long as they have the correct behaviour at nodes. The plot shows linear basis functions and will be our focus. However, increasing the order of accuracy can be done by using more complex basis functions.

Weak form

Focus on time independent case

\[ 0 = \partial_{xx} \psi + S \quad \psi(0) = 0 \quad \partial_x \psi_{x=1} = 0 \, . \]

Switch to weak form: multiply by smooth \(w(x)\), integrate:

\[ \begin{split} 0 = \left[ w(x) \partial_x \psi(x) \right]_0^1 - \int_0^1 \partial_x \psi(x) \partial_x w(x) \, \text{d}x + \\ \int_0^1 w(x) S(x) \, \text{d}x. \end{split} \]

Use basis functions for \(w, \psi, S\): linear system!

Through the method of lines we can separate out time and space steps. So, for simplicity, we’ll concentrate on the solution of elliptic equations to start. Once we know how to form finite element approximations of spatial derivatives on unstructured grids in higher dimensions, then we can add time dependence and solve using the semi-discrete approach.

The advection diffusion equation simplifies down by first choosing the comoving frame (so the advection velocity goes to zero) and then looking for stationary solutions (so the time derivative drops out). We’re left with the diffusion term driven by a source. We’re going to impose boundary conditions of Dirichlet type at the left (no pollutant is allowed at the left boundary) and Neumann type at the right (where the pollutant escapes).

Dealing with two derivatives can be annoying, especially when working with linear shape functions. Equally, we saw in the spectral case that we minimise the residual by projecting along the basis functions. In finite element terms, we can combine this at the analytic level by writing the PDE in weak form. To do this, pick an arbitrary, smooth, weight function \(w(x)\), multiply the PDE by the weight function, and then integrate over the domain. We can integrate by parts to shift one derivative from the solution \(\psi\) onto the weight function \(w\). As the weight function is smooth, this gives us no problems whatsoever.

At this point we can approximate everything using our basis functions. We have the freedom, in principle, to use different basis functions to approximate different terms. However, the simplest (Galerkin) approach is to use the same approximation for all functions. When we do that we will get a linear system for the coefficients.

Discrete form

Use the basis functions: \[ \begin{aligned} && 0 & = - \int_0^1 \partial_x \psi(x) \partial_x w(x) \, \text{d}x + \int_0^1 w(x) S(x) \, \text{d}x \\ \to && & -\sum_{A, B} \psi_A w_B \int_0^1 \partial_x N_A(x) \partial_x N_B(x) + \sum_{B} w_B \int_0^1 N_B(x) S(x) \, . \end{aligned} \]

The weight function is arbitrary: \[ \begin{aligned} 0 &= \sum_B w_B \left\{ K_{AB} \psi_A - F_B \right\} \, , \\ K_{AB} &= \int_0^1 \partial_x N_A(x) \partial_x N_B(x) \, , & F_B &= \int_0^1 N_B(x) S(x) \, . \end{aligned} \]

\(K\) mass or stiffness matrix, \(F\) force vector.

We drop the boundary term for now: we will come back to it later.

We have substituted in the shape functions, using them to approximate the solution and the weight function. We could use them to approximate the source function, but for now it is not needed. Having substituted in we have a double sum over the coefficients approximating the solution and weight function. We then rearrange the sums to isolate the term involving the coefficients of the weight function.

Now, the weight function is acting as a test function: it is arbitrary (and smooth). That means the equation must hold for any non-zero value of the weight function, and hence any non-zero value of the \(w_B\) coefficients. Therefore the term in curly brackets must vanish.

The term in curly brackets re-writes the earlier equation, isolating any pieces involving integrals over shape functions in specific variables. The term multiplying the coefficients of the unknown function we want to find, the \(\psi_A\) coefficients, involves two shape functions, so therefore has two indices. This is naturally written as a matrix. It is usually referred to as the mass or stiffness matrix, depending on the field. The term involving the source function \(S\) can, in principle, be computed directly using a single shape function, as \(S\) is known everywhere. Therefore it is written as a vector, called the force vector.

This combination then gives us a linear system. The matrix \(K\) can be computed explicitly ahead of time, as soon as we have chosen what the shape functions \(N\) are. In time dependent problems we would usually choose time independent shape functions, so even in those cases it only needs computing once. The force vector can be computed from the shape function and the source term. In time dependent problems where the source depends on time this term may need recomputing at each timestep.

Important to note at this point: the \(A, B\) indexes label the nodes. It is best to think of these as labels rather than components, in order to emphasize that in the general case the nodes will be unstructured.

Linear shape functions, even grid

\[ \begin{aligned} N_A(x) & = \begin{cases} (x - x_{A-1}) / \Delta x & x_{A-1} \le x \le x_A \\ (x_{A+1} - x) / \Delta x & x_{A} \le x \le x_{A+1} \\ 0 & \text{otherwise} \end{cases} \\ \partial_x N_A(x) &= \begin{cases} 1 / \Delta x & x_{A-1} \le x \le x_A \\ -1 / \Delta x & x_{A} \le x \le x_{A+1} \\ 0 & \text{otherwise} \end{cases} \\ K_{AB} &= \int_0^1 \partial_x N_A(x) \partial_x N_B(x) \\ & \sim \frac{1}{\Delta x} \begin{pmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 1 \end{pmatrix} \end{aligned} \]

The simplest grid would be evenly spaced. The simplest choice of shape functions would be linear. As we need the shape functions to be one at the associated node and zero at all others, we end up with “tent functions”.

When we plug in the calculation for the derivatives we see that, in this case, they are all constants and proportional to the inverse grid spacing. When we compute the stiffness matrix we multiply two of these derivatives together, getting something that is either zero (if the shape functions do not overlap here) or proportional to the inverse grid spacing squared. Finally, we integrate over the whole domain. However, the individual shape functions are only non-zero on a width of two grid spacings. Therefore the final result will be proportional to the inverse grid spacing again.

The final matrix only has three rows. This is because the boundary condition fixes the value of the unknown at the left boundary, so we do not need to included that entry. In the interior we get terms that look roughly like the coefficients for second order centred in space finite differencing of the second derivative - only the factors of the grid spacing are different. The final row looks slightly different. This is the influence of the boundary condition, which is implicit in the entries.

Example

Source \(S = 2\).

Exact solution \(\psi = 2x - x^2\).

Force vector \(F_B = \Delta x \, (2, 2, \dots, 1)^T\).

Ne = 4; dx = 1 / Ne
psi = np.zeros(Ne+1)
K = np.diag(2*np.ones(Ne))
K -= np.diag(np.ones(Ne-1), -1)
K -= np.diag(np.ones(Ne-1), +1)
K[-1, -1] = 1
K /= dx
F = 2*np.ones(Ne)
F[-1] = 1
F *= dx
psi[1:] = np.linalg.solve(K, F)

This is a very simple example. Remember that we need the boundary conditions to hold where the solution vanishes at the left boundary \(x=0\) and its derivative vanishes at the right boundary \(x=1\). Choosing this solution gives a constant source term.

With the constant source term we can analytically compute the force vector. Again, the final term differs because of the boundary.

We can then set up the code to solve it. First we set up the number of elements, which is here quite small. We then set up the stiffness matrix by setting its diagonal entries - the matrix is tridiagonal, so this is rapid. We have to adjust the entry next to the boundary. We then set up the force vector, following the analytic calculation.

Finally we can solve the linear system. This immediately gives the solutions on the nodes bordering the elements. The first node is not set as its value is fixed by the boundary condition. As it’s zero, there’s no need to modify the result further.

Summary

• Finite elements use function basis.
• Basis not orthogonal: overlap leads to stiffness matrix.
• Use weak form to shift derivatives to test function.
• Galerkin: single set of basis functions for everything.
• Get linear system.

We will want to build methods that take as much advantage of spectral convergence as possible, but are still robust enough to work on complex grids.