Lecture 6 (d)

Nodes vs elements

Unstructured grids

• Most problems have complex domains.
• Use unstructured grids.
• Use Finite Elements.
• Move from nodes as key to elements.

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

Matrix form

Shape functions \(N_A(x)\) take PDE to

\[ \sum_A K_{AB} \psi_A = F_B \, . \]

\(\psi_A\) are nodal values of \(\psi\).

\(K_{AB}\) from all overlapping shape functions at node.

Gets too complex when

• grid is not regular;
• move to higher dimensions.

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.

In the earlier examples we noted that row \(A\) of the stiffness matrix \(K_{AB}\) is linked to node \(A\) of the grid, and hence to the basis function \(N_A\). We then compute the entries in this row by looking at every overlap integral with all basis functions, including itself.

However, this can only be done simply and systematically when the grid is even. In that case all the results are the same, except for when boundaries are involved. In a general unstructured grid we have to redo the calculation for every entry, which is painful.

Element approach

Focus on one element only.

\[ x \in [x_A, x_{A+1}] \to \xi \in [-1, 1]. \]

Two shape functions in element,

\[ N_a = \tfrac{1}{2} (1 + \xi_a \xi), \quad a = 1, 2 \, . \]

These contribute to element stiffness

\[ \begin{aligned} k^{(e)}_{ab} &= \int_{x_A}^{x_{A+1}} \mathrm{d}x \, \partial_x N_a \partial_x N_b \\ &= \int_{-1}^1 \mathrm{d}\xi (\partial_\xi x)^{-1} \partial_\xi N_a \partial_\xi N_b \\ &= \frac{(-1)^{a+b}}{x_{A+1} - x_A} \, . \end{aligned} \]

Rather than re-computing things for every element separately, we can change coordinates in order to map to a reference element. By convention this has coordinates that run from -1 to 1, because we are going to integrate over the element, and Gauss quadrature uses that standard coordinate range for numerical integration.

Within a single element there are two shape functions that are non-zero (in one dimension). We can therefore compute the contribution of the overlap of these shape functions, within this element, giving a two-by-two matrix. That includes the overlap of each shape function with itself, and with the other. In the reference coordinates this integration is very easy. We need to include the determinant of the coordinate transformation in order to get back to the original coordinates. When using linear elements the result is directly given by the size of the element in the original coordinates.

Assembly

All \(k^{e}_{ab}\) add up to \(K_{AB}\). Need \(\{a\}\to\{A\}\) map.

No formula: set location matrix LM s.t.

\[ K_{LM(a, e), LM(b, e)} = K_{LM(a, e), LM(b, e)} + k^{(e)}_{ab} \, . \]

In 1d

• \(LM(1,e)\) is left node of element \(e\)
• \(LM(2,e)\) is right node of element \(e\).

The full domain is made up of lots of elements, not just the one single element. More importantly, the linear system is a set of simultaneous equations for the basis coefficients, which (at least in the linear case) are the values of the functions at the nodes. We therefore need to map the local stiffness matrix labels to the global stiffness matrix.

We also need to note that the local matrix does not contain all the information in the global matrix. It contains all the contributions from within a single element. As, in one dimension, most nodes are next to two elements, we would expect a given row of the global matrix to contain contributions from two local stiffness matrices.

This process is called the assembly of the global problem. We start by setting the global matrix to zero. We then loop over every element, finding its local stiffness matrix, and then adding its entries to the global entries associated with the nodes at the end of the element.

We need a mapping from element numbers to node numbers to do this. This is usually called the location matrix. In one dimension it is usually simple, as all elements are adjacent. In higher dimensions we have to think of things as truly unstructured, so the node numbers in “close” elements could be very different.

The location matrix is something that we have to set up when the grid is set. There’s no mathematics here; it’s just connecting points. In one dimension the easiest thing to do is to count elements from the left, and link to node numbers counted from the left.

One crucial point is that the location matrix should say which nodes are to be ignored. As we noted in the last lecture, the node at the left boundary does not need to be computed as its value is set by the Dirichlet boundary condition. To encode this, any node where a Dirichlet boundary condition is imposed is given the value -1.

Force vector

Write \(S(x) = \sum_a S_a N_a(x)\) in element.

\[ \begin{aligned} f^{(e)}_b &= \int_0^1 \mathrm{d}x N_b(x) S(x) \\ &= \sum_a \int_{-1}^1 \mathrm{d}\xi (\partial_\xi x) N_b N_a S_a \\ & \simeq \frac{x_{A+1} - x_A}{6} \begin{pmatrix} 2 S_1 + S_2 \\ S_1 + 2 S_2 \end{pmatrix} \end{aligned} \]

Better Gauss quadrature methods work for higher accuracy bases.

Assembly: \[ F_{LM(b, e)} = F_{LM(b, e)} + f^{(e)}_{b} \, . \]

In addition to the assembly of the stiffness matrix, we also have to assemble the force vector. First we need to compute it on the reference element. To do this we could analytically calculate the integral of the source term over the element. However, for complex sources this can be impractical.

Instead, we can use numerical quadrature to approximate the integral. In general the shape functions allow us to evaluate any function, when represented in terms of this basis, at any point inside the domain. So we compute the source term at the nodes at the boundary of the element in order to compute the basis coefficients, and from that we can evaluate \(S(\xi)\) wherever we need. Two point Gauss quadrature combined with linear elements gives accurate enough results for these basis elements. When using more complex elements higher order Gauss quadrature can be used without increasing cost excessively.

Boundaries

Weak form includes \([ w(x) \partial_x \psi(x) ]_0^1\).

• Neumann

\[ \partial_x \psi|_{x=1} = \beta \quad \implies \quad F_C = F_C + \beta \]

where \(x_C = 1\).

• Dirichlet

\[ \psi_{x=0} = \alpha \quad \implies \quad F_C = F_C - \alpha k^{(e)}_{12} \]

where \(e\) is the element at boundary, \(x_C\) in \(e\), not at boundary.

Earlier we threw away the boundary term at the analytical level and hand-waved away the impact at the numerical level. For homogeneous boundaries, where the value of the solution and its derivative are set to zero at Dirichlet and Neumann boundaries respectively, it turns out that this is the right thing to do.

At Neumann boundary conditions the value of the derivative appears directly in the formula, as seen from the boundary term that we threw away. This appears only at the boundary node. This can be “moved directly into the force vector”, as it is effectively a constant term applied at a single point.

At Dirichlet boundaries the node is not included in the calculation as its value is already known. However, the value of the node will impact on its neighbour. The details of how this is done is outlined in the notes, but the key point is that again the force vector should be modified, but it’s not the entry at the boundary that’s modified (as that is not included in the calculation), but the other node in the element that touches that Dirichlet boundary. The influence of the value couples through the shape functions, and this scales the boundary value by an entry in the stiffness matrix. When moving to higher dimensions, every node in the element that touches the boundary needs modifying.

Example

Ne = 10
nodes = np.array([0, *sorted(np.random.rand(Ne-1)), 1])
LM = np.zeros((2, Ne), dtype=np.int64)
for e in range(Ne):
    if e==0:
        LM[0, e] = -1
        LM[1, e] = 0
    else:
        LM[0, e] = LM[1, e-1]
        LM[1, e] = LM[0, e] + 1
K = np.zeros((Ne, Ne))
F = np.zeros((Ne,))
for e in range(Ne):
    k_e = stiffness(nodes[e:e+2])
    f_e = force(nodes[e:e+2], S)
    for a in range(2):
        A = LM[a, e]
        for b in range(2):
            B = LM[b, e]
            if (A >= 0) and (B >= 0):
                K[A, B] += k_e[a, b]
        if (A >= 0):
            F[A] += f_e[a]
    # Modify force vector for Dirichlet BC
    if e == 0:
        F[0] -= alpha * k_e[1, 0]
# Modify force vector for Neumann BC
F[-1] += beta
# Solve
Psi_A = np.zeros_like(nodes)
Psi_A[0] = alpha
Psi_A[1:] = np.linalg.solve(K, F)

Here is some example code to illustrate how much work is needed.

The first part of the code is setting up the grid. We construct a grid using randomly spaced nodes inside the domain. This defines the elements, which we then convert into the location matrix, mapping between node and element numbers. Note that we have already fixed the boundary types (the left boundary is Dirichlet and therefore set to -1 in the mapping, whilst all other nodes including the Neumann boundary point are included in the map).

The second part of the code sets up the global stiffness matrix and force vector. It loops over every element. Within each element the local stiffness matrix and force vector are computed, which needs the locations of the nodes bounding the element. We then loop over all the entries, using the location matrix to add these into the global arrays.

Next, we modify the force vector for the boundary conditions. The only reason to do the Dirichlet correction inside the loop is to pick up the correct local stiffness matrix.

Finally we can solve the linear system.

Now, this simplified code will scale to moderate numbers of elements before the matrix uses too much memory for a reasonable machine. Shifting away from simple arrays to sparse matrix structures is needed for practical algorithms.

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.
• Work on elements, not nodes: applies to higher dimensions.

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.