Chapter 8
Finite elements in one dimension
Context: In this chapter, some properties of finite elements are discussed using one-dimensional examples. In particular, we show how second-order PDEs can be discretized using linear elements and how boundary conditions can be incorporated into the finite elements.
8.1 Differentiability of the Basis Functions
The simplest case of a basis of finite elements in one dimension has already been introduced in the previous chapters. The basis functions are “tent” or “hat” functions, each of which is maximal (\(=1\)) on one node and then drops linearly to the two neighboring nodes. In contrast to the Fourier basis, this basis is not differentiable (in the classical sense, not even once). This actually means that we cannot naively use this basis as a series expansion for the ansatz function to solve a second-order (or even first-order) PDE in which derivatives occur. We will now show that a weak solution can be obtained in the context of a weak formulation that does not have to meet these requirements for differentiability. This means that these basis functions can be used for second-order PDEs after all.
As an example, we continue to consider the one-dimensional Poisson equation with the residual \begin {equation} R(x) = \frac {\dif ^2 \Phi }{\dif x^2} + \frac {\rho (x)}{\varepsilon }. \label {eq:strongresidual} \end {equation} In the weighted residual method, the scalar product of the residual with a test function \(v(x)\) must vanish, \((v,R)=0\). In this integral formulation, the rules of partial integration can now be used to transfer a derivative to the test function. This yields \begin {align} (v(x), R(x)) =& \int _a^b \dif x\, v(x) \left ( \frac {\dif ^2 \Phi }{\dif x^2} + \frac {\rho (x)}{\varepsilon }\right ) \label {eq:beforepartialint} \\ =& \int _a^b \dif x\, \left [ \frac {\dif }{\dif x} \left ( v(x) \frac {\dif \Phi }{\dif x} \right ) - \frac {\dif v}{\dif x} \frac {\dif \Phi }{\dif x} \right ] + \int _a^b \dif x\, v(x) \frac {\rho (x)}{\varepsilon } \\ =& \left . v(x) \frac {\dif \Phi }{\dif x} \right |_a^b - \int _a^b \dif x\, \frac {\dif v}{\dif x} \frac {\dif \Phi }{\dif x} + \int _a^b \dif x\, v(x) \frac {\rho (x)}{\varepsilon }. \label {eq:afterpartialint} \end {align}
Equation \eqref{eq:afterpartialint} now contains no second derivative. This means that both the test function \(v(x)\) and the solution \(\Phi (x)\) need only be differentiable once for all \(x\). We can therefore use linear basis functions to discretize a second-order PDE. In the following, we will denote the simulation domain as \(\Omega \). In the one-dimensional case discussed here, \(\Omega =[a,b]\).
Note: The linear basis functions are not even simply differentiable in the classical sense, because the left- and right-sided difference quotient differs at the kinks of the function. However, the so-called weak derivative \(f'(x)\) of the function \(f(x)\) exists if \begin {equation} \int _\Omega \dif x\, v(x) f'(x) = -\int _\Omega \dif x\, v'(x) f(x) \label {eq:weakderivative} \end {equation} for arbitrary (strongly differentiable) test functions \(v(x)\). The weak derivative is not unique; for example, at the kinks of the tent functions, you can assume \(0\) (or any other value) as the weak derivative. This works because this single point makes no contribution to the integral in Eq. \eqref{eq:weakderivative}. Just as the weak derivative is not unique, the weak solution of a differential equation is not unique.
In this learning material, we will not systematically distinguish between strong and weak derivatives but implicitly assume that we are operating with the weak derivative. We often perform intuitive calculations with weak derivatives, such as taking a step function as the derivative of the finite tents and a Dirac \(\delta \)-function as the derivative of the step function. This mathematically imprecise use of derivatives is common in engineering and physics and usually leads to correct results. Spaces of weakly differentiable functions are called Sobolev spaces. Sobolev spaces are always Hilbert spaces, since a scalar product is always needed for weak differentiability, see Eq. \eqref{eq:weakderivative}.
We have thus broadened our understanding of the solution of a PDE. While for the strong solution \(R(x)\equiv 0\), Eq. \eqref{eq:strongresidual}, the function \(\Phi (x)\) must be twice strongly differentiable, for the weak solution \((v,R)\equiv 0\), Eq. \eqref{eq:afterpartialint}, this function must only be once weakly differentiable. Of course, equation \eqref{eq:afterpartialint} must be fulfilled here for all singly- differentiable test functions \(v(x)\).
Note: In the collocation method, the “strong” requirement of differentiability twice is not removed. This can be seen, for example, from the fact that the Dirac \(\delta \) as a test function for the collocation method is not even differentiable in the weak sense. That is, the test functions of the collocation method are not in the set of test functions for which \((v,R)\equiv 0\) is required in the sense of the weak solution. This is the reason why the weighted residual method, and specifically the Galerkin method, has become so widespread.
8.2 Galerkin method
In the context of the weak formulation Eq. \eqref{eq:afterpartialint}, the Galerkin method is now applied. We write again \begin {equation} \Phi (x) \approx \Phi _N(x) = \sum _{n=0}^N a_n \varphi _n(x) \end {equation} as a (finite) series expansion with our finite elements, the tent functions \(\varphi _n(x)\). We first consider the case in which the test function vanishes on the boundary of the domain \(\Omega \), i.e. at \(x=a\) and \(x=b\); thus the first term in Eq. \eqref{eq:afterpartialint} vanishes. This term vanishes, for example, for periodic domains. (However, this term becomes important again when we talk about boundary conditions for non-periodic domains.)
The Galerkin condition thus becomes \begin {equation} (\varphi _k, R_N) = - \sum _n a_n \int _\Omega \dif x\, \frac {\dif \varphi _k}{\dif x} \frac {\dif \varphi _n}{\dif x} + \frac {1}{\varepsilon } \int _\Omega \dif x\, \varphi _k(x) \rho (x) = 0 \label {eq:fegaler1d} \end {equation} with unknown \(a_n\). The integrals in Eq. \eqref{eq:fegaler1d} are merely numbers; the equation is thus a system of coupled linear equations. With \begin {equation} K_{kn} = \int _\Omega \dif x\, \frac {\dif \varphi _k}{\dif x} \frac {\dif \varphi _n}{\dif x} and f_k = \frac {1}{\varepsilon } \int _\Omega \dif x\, \varphi _k(x) \rho (x) \end {equation} we thus obtain \begin {equation} \sum _n K_{kn} a_n = f_k, \label {eq:lineq} \end {equation} or in dyadic (matrix-vector) notation \begin {equation} \t {K}\cdot \v {a} = \v {f}. \label {eq:lineqdyad} \end {equation} Thus, the differential equation is transformed into an algebraic equation that can be solved numerically. The matrix \(\t {K}\) is called the system matrix or stiffness matrix. (The latter term comes from the application of finite elements in the context of structural mechanics.) The term \(\v {f}\) is often referred to as the “right hand side” (often abbreviated as “rhs”) or load vector (again from structural mechanics).
Note: In structural mechanics, which deals with the deformation of solids, \(\t {K}\) is something like a spring constant, \(\v {f}\) a force and \(\v {a}\) the displacements of the nodes, i.e. the distance to the undeformed state. This makes Eq. \eqref{eq:lineqdyad} something like the generalization of a spring law, Hooke’s law. For this reason, the symbols \(\t {K}\) and \(\v {f}\) are traditionally used. Structural mechanics is more complicated than most other numerical applications because the area on which the constitutive equations were discretized changes due to the deformation itself. This automatically leads to non-linear equations, so-called geometric non-linearities. In this case, \(\t {K}\) itself depends on \(\v {a}\).
We can now calculate the elements of the matrix \(\t {K}\) directly. For the finite element basis, the following applies \begin {equation} \frac {\dif \varphi _n(x)}{\dif x} = \left \{ \begin {array}{ll} \frac {1}{x_n - x_{n-1}} & \text {for}\; x\in [x_{n-1},x_n]\\ -\frac {1}{x_{n+1} - x_n} & \text {for}\; x\in [x_n,x_{n+1}] \\ 0 & \text {otherwise} \end {array} \right . \label {eq:finite-element-basis-derivative} \end {equation} and thus \begin {align} K_{nn} &= \frac {1}{x_n - x_{n-1}} + \frac {1}{x_{n+1} - x_n} \\ K_{n,n+1} &= -\frac {1}{x_{n+1} - x_n} \end {align}
and \(K_{kn} = 0\) for \(|n-k|>1\). The matrix \(\t {K}\) is therefore sparse, symmetrical and almost tridiagonal.
Example: For equidistant nodes with spacing \(\Delta x = x_{n+1}-x_n\) on a periodic domain, for example, we obtain \(6\) nodes (\(N=5\)): \begin {equation} \t {K} = \frac {1}{\Delta x} \begin {pmatrix} 2 & -1 & 0 & 0 & 0 & -1 \\ -1 & 2 & -1 & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 & 0 \\ 0 & 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & 0 & -1 & 2 & -1 \\ -1 & 0 & 0 & 0 & -1 & 2 \\ \end {pmatrix} \label {eq:systemmatrix1d-periodic-not-regular} \end {equation} Note that the \(-1\)s in the upper right and lower left corners (\(K_{0N}\) and \(K_{N0}\)) appear due to periodicity. The matrix is therefore not purely tridiagonal. The right-hand side \(f_k\) depends on the specific choice of the source term, i.e. the charge density \(\rho (x)\), and looks as follows: \begin {equation} \v {f} = \begin {pmatrix} (\varphi _0(x), \rho (x))/\varepsilon \\ (\varphi _1(x), \rho (x))/\varepsilon \\ (\varphi _2(x), \rho (x))/\varepsilon \\ (\varphi _3(x), \rho (x))/\varepsilon \\ (\varphi _4(x), \rho (x))/\varepsilon \\ (\varphi _5(x), \rho (x))/\varepsilon \end {pmatrix} \label {eq:rhs1d-periodic-not-regular} \end {equation}
In the course of this learning material, it has been completely neglected so far that boundary conditions are always needed to solve differential equations. Even this periodic case is not complete. In the Fourier representation, \(n=0\) (with wave vector \(q_0=0\)) represents the mean value of the Fourier series. The solution of the Poisson equation does not specify this mean value and it must therefore be given as an additional condition.
In the context of discretization using finite elements, this manifests itself in the fact that the system matrix \(\t {K}\), as written in e.g. Eq. \eqref{eq:systemmatrix1d-periodic-not-regular}, is not regular (also invertible or nonsingular). This can be seen, for example, from the fact that the determinant of \(\t {K}\) vanishes. In other words, the linear equations are not linearly independent. In the case discussed here, the rank of the matrix, i.e. the number of linearly independent rows, is exactly one less than its dimension. We can therefore remove one of these equations (or rows) and introduce the corresponding mean value condition instead. For the periodic case, this is \begin {equation} \int _\Omega \dif x\, \Phi _N(x) = \sum _n a_n \int _\Omega \dif x\, \varphi _n(x) = \sum _n a_n \Delta x = 0. \label {eq:feaverage} \end {equation} The regular system matrix then looks like this, \begin {equation} \t {K} = \frac {1}{\Delta x} \begin {pmatrix} 2 & -1 & 0 & 0 & 0 & -1 \\ -1 & 2 & -1 & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 & 0 \\ 0 & 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & 0 & -1 & 2 & -1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ \end {pmatrix} \label {eq:systemmatrix1d-periodic-regular} \end {equation} whereby \(f_N=0\) must now also apply. The right-hand side thus becomes \begin {equation} \v {f} = \begin {pmatrix} (\varphi _0(x), \rho (x))/\varepsilon \\ (\varphi _1(x), \rho (x))/\varepsilon \\ (\varphi _2(x), \rho (x))/\varepsilon \\ (\varphi _3(x), \rho (x))/\varepsilon \\ (\varphi _4(x), \rho (x))/\varepsilon \\ 0 \end {pmatrix}. \label {eq:rhs1d-periodic-regular} \end {equation} The last row corresponds to the mean value condition Eq. \eqref{eq:feaverage}, but any row could have been replaced by this condition.
8.3 Boundary Conditions
The mean value condition of the previous example is a special case of a boundary condition. In most cases, problems are treated on finite domains \(\Omega \) in which either the function value \(\Phi (x)\) or the derivative \(\dif \Phi /\dif x\) on the boundary \(\partial \Omega \) of the domain is given. (In our one-dimensional case, \(\partial \Omega =\{a,b\}\).) The first case is called a Dirichlet boundary condition, the second case is called a Neumann boundary condition. For differential equations of higher order than two, even higher derivatives could of course occur on the boundary. Combinations of Dirichlet and Neumann boundary conditions are also possible. Here we will only discuss pure Dirichlet and pure Neumann boundary conditions.
8.3.1 Dirichlet Boundary Conditions
In the one-dimensional case discussed here, the Dirichlet boundary conditions are \begin {equation} \Phi (a) \approx \Phi _N(a) = \Phi _a \quad \text {and}\quad \Phi (b) \approx \Phi _N(b) = \Phi _b. \end {equation} These conditions lead directly to \(a_0=\Phi _a\) and \(a_N=\Phi _b\). That is, the Dirichlet conditions directly fix the corresponding coefficients of the series expansion. Note that the conditions \((\varphi _0, R)=0\) and \((\varphi _N, R)=0\) have been implicitly taken from the system of equations. However, the basis functions \(\varphi _0(x)\) and \(\varphi _N(x)\) still appear in the series expansion \(\Phi _N(x)\).
Example: In our example, the system matrix with Dirichlet boundary conditions then becomes \begin {equation} \t {K} = \frac {1}{\Delta x} \begin {pmatrix} 1 & 0 & 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 & 0 & 1 \\ \end {pmatrix} \label {eq:systemmatrix1d-dirichlet} \end {equation} and \(f_0=\Phi _a/\Delta x\) and \(f_N=\Phi _b/\Delta x\): \begin {equation} \v {f} = \begin {pmatrix} \Phi _a/\Delta x \\ (\varphi _1(x), \rho (x))/\varepsilon \\ (\varphi _2(x), \rho (x))/\varepsilon \\ (\varphi _3(x), \rho (x))/\varepsilon \\ (\varphi _4(x), \rho (x))/\varepsilon \\ \Phi _b/\Delta x \end {pmatrix} \label {eq:rhs1d-dirichlet} \end {equation} This matrix \(\t {K}\) is regular.
The \(\Delta x\) only appears on the right side \(f_0\) and \(f_N\) because it is a prefactor in the system matrix. In an implementation, one would push \(\Delta x\) completely to the right side and it would thus completely disappear from the first and last line.
8.3.2 Neumann boundary conditions
The Neumann boundary conditions are \begin {equation} \left .\frac {\dif \Phi }{\dif x}\right |_a \approx \left .\frac {\dif \Phi _N}{\dif x}\right |_a = \Phi ^\prime _a \quad \text {and}\quad \left .\frac {\dif \Phi }{\dif x}\right |_b \approx \left .\frac {\dif \Phi _N}{\dif x}\right |_b = \Phi ^\prime _b. \end {equation} To include these conditions in our algebraic equation, we can no longer neglect the first term in Eq. \eqref{eq:afterpartialint}. This term is not relevant for the Dirichlet boundary conditions, since the basis functions \(\varphi _1(x)\) to \(\varphi _{N-1}(x)\) vanish on the boundary \(x=a,b\). The basis functions \(\varphi _0(x)\) and \(\varphi _N(x)\) do not vanish, but in the Dirichlet case they no longer appear in our set of test functions.
However, we now have to determine how to interpret the basis functions \(\varphi _0(x)\) and \(\varphi _N(x)\), which now extend beyond the boundary of the domain. A natural interpretation is to consider only the half of the “tent” that remains in the domain.
The Galerkin approach thus leads to \begin {equation} (\varphi _k, R_N) = \left . \varphi _k(x) \frac {\dif \Phi }{\dif x} \right |_a^b - \sum _n a_n \int _\Omega \dif x\, \frac {\dif \varphi _k}{\dif x} \frac {\dif \varphi _n}{\dif x} + \frac {1}{\varepsilon }(\varphi _k(x), \rho (x)) = 0. \end {equation} The additional term on the left-hand side only plays a role at \(k=0\) and \(k=N\). We obtain \begin {equation} (\varphi _0, R_N) = - \Phi ^\prime _a - \sum _n K_{0n} a_n + \frac {1}{\varepsilon }(\varphi _0(x), \rho (x)) \end {equation} with \begin {equation} K_{00} = \frac {1}{x_1-x_0} \quad \text {and}\quad K_{01} = -\frac {1}{x_1 - x_0} \end {equation} and all other \(K_{0n}=0\). We can write this again (see Eq. \eqref{eq:lineq}) with \begin {equation} f_0 = \frac {1}{\varepsilon }(\varphi _0(x), \rho (x)) - \Phi ^\prime _a \end {equation} A corresponding set of equations applies to the right boundary with \(k=N\).
Example: In our example, the system matrix for two Neumann boundary conditions then becomes \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:systemmatrix1d-neumann} \end {equation} Note that this system matrix is slightly different from the one for Dirichlet boundary conditions, Eq. This matrix is not regular and thus the problem with two Neumann boundary conditions is indeterminate. The reason for this is the same as in the periodic case: the Neumann boundary conditions do not fix the absolute value (mean) of the potential \(\Phi \). So you either need a Dirichlet boundary condition (left or right) or again the fixing of the mean value.
The right side becomes: \begin {equation} \v {f} = \begin {pmatrix} (\varphi _0(x), \rho (x))/\varepsilon - \Phi ^\prime _a \\ (\varphi _1(x), \rho (x))/\varepsilon \\ (\varphi _2(x), \rho (x))/\varepsilon \\ (\varphi _3(x), \rho (x))/\varepsilon \\ (\varphi _4(x), \rho (x))/\varepsilon \\ (\varphi _5(x), \rho (x))/\varepsilon + \Phi ^\prime _b \end {pmatrix}. \label {eq:rhs1d-neumann} \end {equation} Note that for the evaluation of \((\varphi _0, \rho )\) and \((\varphi _5, \rho )\) only half of the corresponding tent needs to be integrated.
8.4 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 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)}_{0}(\xi ^{(n)}(x)) + a_{n+1} N^{(n)}_{1}(\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_{0}^{(n)}(\xi ) = 1-\xi \quad \text {and}\quad N_{1}^{(n)}(\xi ) = \xi \label {eq:linel1e} \end {equation} with \(\xi \in [0,1]\) are called shape functions (also known as ansatz 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 “form”).
Note: Note that there is exactly one element less than nodes, 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 \(0\) for the left node and \(1\) for the right node) and the element index yields the global node index (here \(n\)).
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_{0}^{(n)}\) on the left and \(N_{1}^{(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(x)\) is defined only on the individual element \(n\), i.e. they live on \(\Omega _n\). Here \(i\) denotes the node within the element. 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)}_{0}(x) + N^{(n-1)}_{1}(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, see also Eq.
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=0}^{N-1} \phi ^{(n)}(x) = \sum _{n=0}^{N-1} \left ( a_n N^{(n)}_{0}(x) + a_{n+1} N^{(n)}_{1}(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)}_{0} + N^{(k-1)}_{1}, R) = 0 \label {eq:galerkinform} \end {equation} with \begin {align} (N^{(k)}_{I}, R) =& \sum _{n=0}^{N-1} \left ( a_n (N^{(k)}_{I}, \mathcal {L} N^{(n)}_{0}) + a_{n+1} (N^{(k)}_{I}, \mathcal {L} N^{(n)}_{1}) \right ) - (N^{(k)}_{I}, f) \label {eq:formprod1} \\ =& a_k (N^{(k)}_{I}, \mathcal {L} N^{(k)}_{0}) + a_{k+1} (N^{(k)}_{I}, \ Con{L} N^{(k)}_{1}) - (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. 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_I{}^{(n)}_{J} = (N_I{}^{(n)}, \mathcal {L} N_J{}^{(n)}) \end {equation} and \begin {equation} f_I = (N_I, f). \end {equation} In this notation, we obtain \begin {equation} (N_I, R) = \sum _{J=0}^1 K_I J a_{k+J} - f_I \end {equation} where \(J\) runs over the nodes within the element \(k\) and thus (in one dimension) \(k+J\) is the global node index of the corresponding left or right node. The Galerkin condition Eq. \eqref{eq:galerkinform} thus becomes \begin {equation} \begin {split} (\varphi _k, R) =& (N_0,R) + (N_1,R) \\ =& \sum _{I=0}^1 \left ( \sum _{J=0}^1 K_IJ a - f_I \right ) = 0. \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\)). The assembly of the system matrix is often referred to as ’assembly’.
Example: We reformulate the example problem (Poisson equation) here again in the context of the form functions. This is a process that requires three steps:
-
Element matrix and right side: Since all elements are identical, the individual element matrices are also identical. First, we apply the trick with partial integration again 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)}_{0}}{\dif x} &= -\frac {1}{\Delta x} \\ \frac {\dif N^{(n)}_{1}}{\dif x} &= \frac {1}{\Delta x} \end {align}
This gives the element stiffness matrix \begin {equation} \t {K}^{(n)} = \frac {1}{\Delta x} \begin {pmatrix} 1 & -1 \\ -1 & 1 \end {pmatrix}, \label {eq:elmat1d} \end {equation} which is identical for all elements. The right-hand side for the elements becomes \begin {equation} \t {f}^{(n)} = \begin {pmatrix} (N^{(n)}_{0}, \rho )/\varepsilon \\ (N^{(n)}_{1}, \rho )/\varepsilon \end {pmatrix}. \end {equation} This right-hand side can be different for the different elements if \(\rho \) varies spatially.
- Assembly of the system matrix and the right-hand side: 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=0} + \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=1} + \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} The advantage of the formulation with shape functions is that you 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 form functions. Here, the corresponding rows of the system matrix and the right-hand side must be replaced by the corresponding boundary condition. The matrix \(\t {K}\) from Eq. \eqref{eq:sysmatform} is not regular without the corresponding boundary conditions.