Chapter 9
Assembly
Context: Deriving the discretized equations in matrix form can be tedious. This chapter describes how to derive them from the properties of the individual elements. This leads to the concept of a per-element stiffness matrix. The global system matrix can be constructing by assembling the individual per-element matrices. This splits the complexity of the discretization into two steps: Determining the element matrices and assembling them.
9.1 Shape functions
The finite element method is often formulated not with the help of basis functions but with the help of so-called shape functions. The reason for this is that the shape functions provide a more intuitive approach to the interpolation rule behind the finite elements, thus allowing for easier extension to multiple dimensions. In the one-dimensional case, the shape functions may appear more complicated than the approach outlined above, but for the formulation of the finite element method in several dimensions, this will be the simpler approach.
We first have to talk about what an element is. The tent functions lead to a linear interpolation between the nodes. The areas between the nodes are called elements. In one dimension, these are one-dimensional intervals, in higher dimensions the elements can take on complex forms. A linear basis function, as introduced here, is centered on the node and non-zero on two elements.
Instead of defining the interpolation in terms of the basis functions, we can also demand that within an element the function values on the knots be interpolated in a given form, here initially linearly. For the \(n\)-th element between the nodes \(n\) and \(n+1\), i.e. \(x\in \Omega ^{(n)}=[x_n,x_{n+1}]\), only the basis functions \(\varphi _n\) and \(\varphi _{n+1}\) contribute, since all other basis functions vanish on this interval. Here, \(\Omega ^{(n)}\) denotes the domain of the \(n\)th element. Thus, in this element, the function \(\Phi _N(x)\) has the form \begin {equation} \begin {split} \phi ^{(n)}(x) =& a_n \varphi _n(x) + a_{n+1}\varphi _{n+1}(x) \\ =& a_n \frac {x_{n+1}-x}{x_{n+1} - x_n} + a_{n+1}\frac {x-x_n}{x_{n+1} - x_n} \\ =& a_n N^{(n)}_{1}(\xi ^{(n)}(x)) + a_{n+1} N^{(n)}_{2}(\xi ^{(n)}(x)), \label {eq:basistoform} \end {split} \end {equation} with \(\xi ^{(n)}(x)=(x-x_n)/(x_{n+1}-x_n)\) and \(x\in \Omega ^{(n)}\). Here and in the following, superscripts \(x^{(n)}\) denote elements and subscripts \(x_n\) denote nodes. The functions \begin {equation} N_{1}^{(n)}(\xi ) = 1-\xi \quad \text {and}\quad N_{2}^{(n)}(\xi ) = \xi \label {eq:linel1d} \end {equation} with \(\xi \in [0,1]\) are called shape functions and \(\xi ^{(n)}(x)\) is a rescaling function that becomes \(0\) at the left edge of the \(n\)th element and \(1\) at the right edge of the element. This decouples the size of the element from the interpolation rule (the “shape” of the element).
Note: There is exactly one element less than there are nodes in 1D, which we have so far also denoted by the index \(n\). In the one-dimensional case, the relationship between the global node index and the element index is trivial; in higher dimensions, keeping track of the indices becomes complicated. The combination of the local index of the node within an element (here \(1\) for the left node and \(2\) for the right node) and the element index yields the global node index (here \(n\)). We will use capital Latin letters \(I, J,K\) to indicate local node indices.
Note: The interpolation rule depends on the element index \(n\) because the elements can have different sizes. There is a shape function for each node of the element, which are referred to here as \(N_{1}^{(n)}\) on the left and \(N_{2}^{(n)}\) on the right. The set of shape functions of an element determine the element type. Equations \eqref{eq:linel1d} define a linear element. In principle, the element types can be different for each element. However, for the basis set we have used so far, the left and right shape functions (and thus the element type) are identical for all elements.
Note: The basis functions \(\varphi _n(x)\) are defined globally, i.e. they live on the entire simulation domain \(\Omega \) (but vanish in sections of it). The shape function \(N_I^{(n)}(x)\) is defined only on the individual element \(n\), i.e. they live on \(\Omega _n\). Here \(I\) denotes the node within the element, i.e. it is a local note index. However, we can of course express the basis functions using the shape functions by collecting the functions that have the corresponding expansion coefficient \(a_n\) as a prefactor. In the one-dimensional case, we obtain \begin {equation} \varphi _n(x) = N^{(n)}_{1}(x) + N^{(n-1)}_{2}(x) \label {eq:formtobasis} \end {equation} with the compact notation \(N^{(n)}_{I}(x)=N^{(n)}_{I}(\xi ^{(n)}(x))\). Equation \eqref{eq:formtobasis} is to be interpreted such that the shape functions vanish if the argument is not in the element, i.e., for \(x\not \in \Omega ^{(n)}\). In two or three dimensions, it is usually easier to work with the shape functions than with the basis functions. Conversely, the shape functions are ultimately the part of the basis functions that lives on the elements.
Shape functions are useful because they can be used to write the approximated solution as a sum over elements, i.e. in the one-dimensional case \begin {equation} \Phi _N(x) = \sum _{n=1}^{N} \phi ^{(n)}(x) = \sum _{n=1}^{N} \left ( a_n N^{(n)}_{1}(x) + a_{n+1} N^{(n)}_{2}(x) \right ). \end {equation} For a general PDE, \(R = \mathcal {L} u_N - f\), the Galerkin condition becomes \begin {equation} (\varphi _k, R) = (N^{(k)}_{1} + N^{(k-1)}_{2}, R) = 0 \label {eq:galerkinform} \end {equation} with \begin {align} (N^{(k)}_{I}, R) =& \sum _{n=1}^{N} \left ( a_n (N^{(k)}_{I}, \mathcal {L} N^{(n)}_{1}) + a_{n+1} (N^{(k)}_{I}, \mathcal {L} N^{(n)}_{2}) \right ) - (N^{(k)}_{I}, f) \label {eq:formprod1} \\ =& a_k (N^{(k)}_{I}, \mathcal {L} N^{(k)}_{1}) + a_{k+1} (N^{(k)}_{I}, \mathcal {L} N^{(k)}_{2}) - (N^{(k)}_{I}, f) \end {align}
where the sum in Eq. \eqref{eq:formprod1} disappears because the shape functions on different elements are trivially orthogonal. (The shape function itself is nonzero only on a single element.) Here and in the following, capital letters \(I\) denote local node indices, while small letters \(i\) denote global node indices.
Motivated by this equation, we define an element matrix (or element stiffness matrix) for element \(n\) as \begin {equation} K_{IJ}^{(n)} = (N_{I}^{(n)}, \mathcal {L} N_{J}^{(n)}) \end {equation} and \begin {equation} f_I^{(n)} = (N_I^{(n)}, f). \end {equation} In this notation, we obtain \begin {equation} (N_I^{(n)}, R) = \sum _{J=1}^2 K_{IJ}^{(n)} a_{J}^{(n)} - f_I^{(n)} \end {equation} where \(J\) runs over the nodes within the element \(n\) and thus (in one dimension) \(n+J-1\) is the global node index of the corresponding left or right node, i.e. \(a_J^{(n)}=a_{n+J-1}\). The Galerkin condition Eq. \eqref{eq:galerkinform} thus becomes \begin {equation} \begin {split} 0=(\varphi _k, R) =& (N_1,R) + (N_2,R) \\ =& \sum _{J=1}^2 K_{1J}^{(k)} a_{J}^{(k)} - f_1^{(k)} + \sum _{J=1}^2 K_{2J}^{(k-1)} a_{J}^{(k-1)} - f_2^{(k-1)} \end {split} \label {eq:galerkinform2} \end {equation} This condition corresponds to a row of the system matrix and the right-hand side. Row \(k\) of the system matrix, which corresponds to node \(k\), thus has a contribution from element \(k\) (whose left node is node \(k\)) and element \(k-1\) (whose right node is node \(k\)). This process of constructing or “assembling” the global system matrix from the element matrices is often referred to as assembly.
Note: The element matrix is constructed by evaluating \begin {equation} (N_I^{(n)}, R). \end {equation} Since the shape function is nonzero only on element \(n\), we can restrict the integration inside the scalar product to just this element.
9.2 Assembling the system matrix
As an example, we reformulate the example problem (Poisson equation) here again with shape functions. This is a process that requires three steps:
-
Element matrix and load vector: Since all elements are identical, the individual element matrices are also identical. First, we apply the partial integration to reduce the condition of differentiability. We get the expression \begin {equation} K_{IJ}^{(n)} = \left (N^{(n)}_{I}, \frac {\dif ^2 N^{(n)}_{J}}{\dif x^2} \right ) = -\left (\frac {\dif N^{(n)}_{I}}{\dif x}, \frac {\dif N^{(n)}_{J}}{\dif x}\right ), \end {equation} which only contains first derivatives of the shape functions. These derivatives are constants because the shape functions are linear, Eq. \eqref{eq:linel1d}, \begin {align} \frac {\dif N^{(n)}_{1}}{\dif x} &= -\frac {1}{\Delta x^{(n)}} \\ \frac {\dif N^{(n)}_{2}}{\dif x} &= \frac {1}{\Delta x^{(n)}}, \end {align}
where \(\Delta x^{(n)}\) is the size of element \(n\). This gives the element stiffness matrix \begin {equation} \t {K}^{(n)} = \frac {1}{\Delta x^{(n)}} \begin {pmatrix} -1 & 1 \\ 1 & -1 \end {pmatrix}, \label {eq:elmat1d} \end {equation} which is identical for all elements if they have the same size and use the same shape functions. The right-hand side for the elements becomes \begin {equation} \v {f}^{(n)} = \begin {pmatrix} (N^{(n)}_{1}, \rho )/\varepsilon \\ (N^{(n)}_{2}, \rho )/\varepsilon \end {pmatrix}. \end {equation} This right-hand side can be different for the different elements if \(\rho \) varies spatially.
- Assembly: The rules for assembling the system matrix follow from the Galerkin condition, Eq. \eqref{eq:galerkinform2}. To do this, the \(2\times 2\) element matrix must be expanded to the \(6\times 6\) system matrix and summed: rows and columns of the element matrix correspond to element nodes and these must now be mapped to the corresponding global nodes in the system matrix. In our one-dimensional case, this is trivial. We get \begin {equation} \t {K}\Delta x = \underbrace { \begin {pmatrix} -1 & 1 & 0 & 0 & 0 & 0 \\ 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end {pmatrix} }_{\text {Element}\,n=1} + \underbrace { \begin {pmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -1 & 1 & 0 & 0 & 0 \\ 0 & 1 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end {pmatrix} }_{\text {Element}\,n=2} + \cdots \end {equation} and thus \begin {equation} \t {K} = \frac {1}{\Delta x}\begin {pmatrix} -1 & 1 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 1 & -1 \\ \end {pmatrix}. \label {eq:sysmatform} \end {equation} where we have assumed that all elements have the same size \(\Delta x=\Delta x^{(n)}\). The advantage of the formulation with shape functions is that we only have to calculate the element matrix once for an element type, and then you can simply assemble the system matrix from it. The assembly of the right-hand side follows analogously, but is easier since it is a vector.
- Boundary conditions: We have not yet discussed boundary conditions in the context of the shape functions. Here, the rows of the system matrix and the load vector must be replaced by the corresponding boundary condition. The matrix \(\t {K}\) from Eq. \eqref{eq:sysmatform} is singular without the corresponding boundary conditions.
9.3 Nonuniform one-dimensional grids
Let us construct system matrix for a system consisting of three elements with different sizes \(\Delta x^{(1)}\), \(\Delta x^{(2)}\) and \(\Delta x^{(3)}\). The element matrices for the operator \(\mathcal {L}=\dif ^2/\dif x^2\) are given by Eq. \eqref{eq:elmat1d}. We simply sum the three element matrices, yielding \begin {equation} \begin {split} \t {K} &= \begin {pmatrix} -\frac {1}{\Delta x^{(1)}} & \frac {1}{\Delta x^{(1)}} & 0 & 0 \\ \frac {1}{\Delta x^{(1)}} & -\frac {1}{\Delta x^{(1)}} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end {pmatrix} + \begin {pmatrix} 0 & 0 & 0 & 0 \\ 0 & -\frac {1}{\Delta x^{(2)}} & \frac {1}{\Delta x^{(2)}} & 0 \\ 0 & \frac {1}{\Delta x^{(2)}} & -\frac {1}{\Delta x^{(2)}} & 0 \\ 0 & 0 & 0 & 0 \end {pmatrix} + \begin {pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & -\frac {1}{\Delta x^{(3)}} & \frac {1}{\Delta x^{(3)}} \\ 0 & 0 & \frac {1}{\Delta x^{(3)}} & -\frac {1}{\Delta x^{(3)}} \end {pmatrix} \\ &= \begin {pmatrix} -\frac {1}{\Delta x^{(1)}} & \frac {1}{\Delta x^{(1)}} & 0 & 0 \\ \frac {1}{\Delta x^{(1)}} & -\frac {1}{\Delta x^{(1)}}-\frac {1}{\Delta x^{(2)}} & \frac {1}{\Delta x^{(2)}} & 0 \\ 0 & \frac {1}{\Delta x^{(2)}} & -\frac {1}{\Delta x^{(2)}}-\frac {1}{\Delta x^{(3)}} & \frac {1}{\Delta x^{(3)}} \\ 0 & 0 & \frac {1}{\Delta x^{(3)}} & -\frac {1}{\Delta x^{(3)}} \end {pmatrix}. \end {split} \end {equation} This expression would have been more difficult to derive from using the basis functions directly.
9.4 Element matrices
We have in this chapter discussed the discretization of the second derivative \(\mathcal {L}=\dif ^2/\dif x^2\). The element matrix for zero and first derivatives in one-dimensions with linear shape functions are summarized in Tab. 9.1. Deriving the explicit expression in this table requires solving the integrals behind the scalar products shown in the third column. Those integrals are at most over quadratic functions and hence trivial to carry out.
| Name | Operator \(\mathcal {L}\) | \(K_{IJ}\) | \(\t {K}\) |
| Mass matrix | \(1\) | \(\left (N_I,N_J\right )\) | \(\Delta x \begin {pmatrix} 1/3 & 1/6 \\ 1/6 & 1/3 \end {pmatrix}\) |
| \(\frac {\dif }{\dif x}\) | \(\left (N_I,\frac {\dif N_J}{\dif x}\right )\) | \(\begin {pmatrix} -1/2 & 1/2 \\ -1/2 & 1/2 \end {pmatrix}\) | |
| Laplacian matrix | \(\frac {\dif ^2}{\dif x^2}\) | \(-\left (\frac {\dif N_I}{\dif x},\frac {\dif N_J}{\dif x}\right )\) | \(\frac {1}{\Delta x} \begin {pmatrix} -1 & 1 \\ 1 & -1 \end {pmatrix}\) |
9.5 Implementation
Example: As an example of a differential equation, we now discuss the numerical solution of the (linearized) Poisson-Boltzmann (PB) equation. For two identical species with opposite charge, the PB equation is given by \begin {equation} \begin {split} \nabla ^2 \Phi &= - \frac {c^\infty }{\varepsilon } \left [ q_+ \exp \left (-\frac {q_+ \Phi }{k_B T}\right ) + q_- \exp \left (-\frac {q_- \Phi }{k_B T}\right ) \right ] \\ &= \frac {2\rho _0}{\varepsilon } \sinh \left ( \frac {|e| \Phi }{k_B T} \right ) \end {split} \end {equation} where \(\rho _0=|e|c^\infty \) is a reference charge density and \(q_+=|e|\) and \(q_-=-|e|\) are the ionic charges of the two species. Since for small \(x\) we have \(\sinh x\approx x\), the linearized version of the PB equation is \begin {equation} \nabla ^2\Phi =\Phi /\lambda ^2 \end {equation} with the Debye length \(\lambda =\sqrt {\varepsilon k_B T/(2|e|\rho _0)}\).
We here show how to implement the solution of the linearized Poisson-Boltzmann equation with the residual \begin {equation} R=\nabla ^2\Phi -\Phi /\lambda ^2 \end {equation} with the concepts described in this chapter. This equation employs two operators, the Laplacian and the unit operator, which means each element has a contribution from the Laplacian and mass matrices (see Tab. 9.1). The overall element matrix of this differential equation is therefore given by \begin {equation} \t {K}^{(n)} = \t {K}^{(n)}_\text {Laplacian} - \t {K}^{(n)}_\text {mass} / \lambda ^2. \end {equation} We implement these matrices as individual functions, that take the width \(\Delta x\) of each element as input:
1def laplacian_1d(width): 2 """ 3 Return the Laplacian element matrix for a linear 1D 4 element of width ‘width‘. 5 """ 6 return np.array([[-1, 1], [1, -1]]) / width 7 8def mass_1d(width): 9 """ 10 Returns the mass element matrix for a linear 1D 11 element of width ‘width‘. 12 """ 13 return np.array([[2, 1], [1, 2]]) * width / 6 14 15def pb_1d(width, lam): 16 """ 17 Returns the Poisson-Boltzmann element matrix for a linear 18 1D element of width ‘width‘ and a Debye length ‘lam‘. 19 """ 20 return laplacian_1d(width) - mass_1d(width) / lam
Given the nodal positions, we can now compute an array containing the matrices for each element:
1lam = 0.5 # Debye length 2x_g = np.array([0, 0.1, 0.2, 0.5, 1.0, 2.0]) # Node positions 3widths_e = np.diff(x_g) # Element widths 4K_ell = np.array([pb_1d_element(width, lam) # El. matrices 5 for width in widths_e])
What is now left is to assemble these into a global stiffness matrix \(\t {K}\). We write a utility function that achieves this for the one-dimensional case:
1def assemble_1d(K_ell): 2 """ 3 Assemble the global matrix from the element matrices. 4 """ 5 e, l, _ = K_ell.shape # Nb elements, nb element nodes 6 K_gg = np.zeros((e+1, e+1)) 7 for i in range(e): # Loop over all elements and sum 8 K_gg[i:i+l, i:i+l] += K_ell[i] 9 return K_gg
Running this code for the above node positions yields the system matrix with values \begin {equation} \t {K} = \begin {pmatrix} -10.07 & 9.97 & 0.00 & 0.00 & 0.00 & 0.00\\ 9.97 & -20.13 & 9.97 & 0.00 & 0.00 & 0.00\\ 0.00 & 9.97 & -13.60 & 3.23 & 0.00 & 0.00\\ 0.00 & 0.00 & 3.23 & -5.87 & 1.83 & 0.00\\ 0.00 & 0.00 & 0.00 & 1.83 & -4.00 & 0.67\\ 0.00 & 0.00 & 0.00 & 0.00 & 0.67 & -1.67 \end {pmatrix}. \end {equation}
The code above uses a notation where the suffixes indicate what an array dimensions represents. The suffix _e indicates and element index, suffix _l a local (per-element) node and suffix _g denotes the index of the global node. Suffixes are combined to show overall dimensions of an array. For example, variable K_gg contains the system matrix \(\t {K}\) and therefore needs two global node indices. Futhermore, K_ell contains element matrices for each element in the system and therefore needs the three indices, one for the element and the other two for the matrix dimensions.