Variational formulations with finite elements¶

We shall now take the ideas from previous chapters and put together such that we can solve PDEs using the flexible finite element basis functions. This is quite a machinery with many details, but the chapter is mostly an assembly of concepts and details we have already met.

Computing with finite elements¶

The purpose of this section is to demonstrate in detail how the finite element method can then be applied to the model problem

$$-u''(x) = 2,\quad x\in (0,L),\ u(0)=u(L)=0,$$

with variational formulation

$$(u',v') = (2,v)\quad\forall v\in V{\thinspace .}$$

Any $v\in V$ must obey $v(0)=v(L)=0$ because of the Dirichlet conditions on $u$.

Finite element mesh and basis functions¶

We introduce a finite element mesh with $N_e$ cells, all with length $h$, and number the cells from left to right. Choosing P1 elements, there are two nodes per cell, and the coordinates of the nodes become

$$x_{i} = i h,\quad h=L/N_e,\quad i=0,\ldots,N_n-1=N_e{\thinspace .}$$

Any node $i$ is associated with a finite element basis function ${\varphi}_i(x)$. When approximating a given function $f$ by a finite element function $u$, we expand $u$ using finite element basis functions associated with all nodes in the mesh. The parameter $N$, which counts the unknowns from 0 to $N$, is then equal to $N_n-1$ such that the total number of unknowns, $N+1$, is the total number of nodes. However, when solving differential equations we will often have $N<N_n-1$ because of Dirichlet boundary conditions. The reason is simple: we know what $u$ are at some (here two) nodes, and the number of unknown parameters is naturally reduced.

In our case with homogeneous Dirichlet boundary conditions, we do not need any boundary function $B(x)$, so we can work with the expansion

$$\begin{equation} u(x) = \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x){\thinspace .} \label{fem:deq:1D:fem1:ex:u} \tag{64} \end{equation}$$

Because of the boundary conditions, we must demand ${\psi}_i(0)={\psi}_i(L)=0$, $i\in{\mathcal{I}_s}$. When ${\psi}_i$ for all $i=0,\ldots,N$ is to be selected among the finite element basis functions ${\varphi}_j$, $j=0,\ldots,N_n-1$, we have to avoid using ${\varphi}_j$ functions that do not vanish at $x_{0}=0$ and $x_{N_n-1}=L$. However, all ${\varphi}_j$ vanish at these two nodes for $j=1,\ldots,N_n-2$. Only basis functions associated with the end nodes, ${\varphi}_0$ and ${\varphi}_{N_n-1}$, violate the boundary conditions of our differential equation. Therefore, we select the basis functions ${\varphi}_i$ to be the set of finite element basis functions associated with all the interior nodes in the mesh:

$${\psi}_i={\varphi}_{i+1},\quad i=0,\ldots,N{\thinspace .}$$

The $i$ index runs over all the unknowns $c_i$ in the expansion for $u$, and in this case $N=N_n-3$.

In the general case, and in particular on domains in higher dimensions, the nodes are not necessarily numbered from left to right, so we introduce a mapping from the node numbering, or more precisely the degree of freedom numbering, to the numbering of the unknowns in the final equation system. These unknowns take on the numbers $0,\ldots,N$. Unknown number $j$ in the linear system corresponds to degree of freedom number $\nu (j)$, $j\in{\mathcal{I}_s}$. We can then write

$${\psi}_i={\varphi}_{\nu(i)},\quad i=0,\ldots,N{\thinspace .}$$

With a regular numbering as in the present example, $\nu(j) = j+1$, $j=0,\ldots,N=N_n-3$.

Computation in the global physical domain¶

We shall first perform a computation in the $x$ coordinate system because the integrals can be easily computed here by simple, visual, geometric considerations. This is called a global approach since we work in the $x$ coordinate system and compute integrals on the global domain $[0,L]$.

The entries in the coefficient matrix and right-hand side are

$$A_{i,j}=\int_0^L{\psi}_i'(x){\psi}_j'(x) {\, \mathrm{d}x},\quad b_i=\int_0^L2{\psi}_i(x) {\, \mathrm{d}x}, \quad i,j\in{\mathcal{I}_s}{\thinspace .}$$

Expressed in terms of finite element basis functions ${\varphi}_i$ we get the alternative expressions

$$A_{i,j}=\int_0^L{\varphi}_{i+1}'(x){\varphi}_{j+1}'(x) {\, \mathrm{d}x},\quad b_i=\int_0^L2{\varphi}_{i+1}(x) {\, \mathrm{d}x},\quad i,j\in{\mathcal{I}_s}{\thinspace .}$$

For the following calculations the subscripts on the finite element basis functions are more conveniently written as $i$ and $j$ instead of $i+1$ and $j+1$, so our notation becomes

$$A_{i-1,j-1}=\int_0^L{\varphi}_{i}'(x){\varphi}_{j}'(x) {\, \mathrm{d}x},\quad b_{i-1}=\int_0^L2{\varphi}_{i}(x) {\, \mathrm{d}x},$$

where the $i$ and $j$ indices run as $i,j=1,\ldots,N+1=N_n-2$.

The ${\varphi}_i(x)$ function is a hat function with peak at $x=x_{i}$ and a linear variation in $[x_{i-1},x_{i}]$ and $[x_{i},x_{i+1}]$. The derivative is $1/h$ to the left of $x_{i}$ and $-1/h$ to the right, or more formally,

$$\begin{equation} {\varphi}_i'(x) = \left\lbrace\begin{array}{ll} 0, & x < x_{i-1},\\ h^{-1}, & x_{i-1} \leq x < x_{i},\\ -h^{-1}, & x_{i} \leq x < x_{i+1},\\ 0, & x\geq x_{i+1} \end{array} \right. \label{fem:approx:fe:Dphi:1:formula2} \tag{65} \end{equation}$$

Figure shows ${\varphi}_2'(x)$ and ${\varphi}_3'(x)$.

Illustration of the derivative of piecewise linear basis functions associated with nodes in cell 2.

We realize that ${\varphi}_i'$ and ${\varphi}_j'$ has no overlap, and hence their product vanishes, unless $i$ and $j$ are nodes belonging to the same cell. The only nonzero contributions to the coefficient matrix are therefore

$$\begin{align*} A_{i-1,i-2} &=\int_0^L{\varphi}_i'(x) {\varphi}_{i-1}'(x) {\, \mathrm{d}x},\\ A_{i-1,i-1}&=\int_0^L{\varphi}_{i}'(x)^2 {\, \mathrm{d}x}, \\ A_{i-1,i}&=\int_0^L{\varphi}_{i}'(x){\varphi}_{i+1}'(x) {\, \mathrm{d}x}, \end{align*}$$

for $i=1,\ldots,N+1$, but for $i=1$, $A_{i-1,i-2}$ is not defined, and for $i=N+1$, $A_{i-1,i}$ is not defined.

From Figure, we see that ${\varphi}_{i-1}'(x)$ and ${\varphi}_i'(x)$ have overlap of one cell $\Omega^{(i-1)}=[x_{i-1},x_{i}]$ and that their product then is $-1/h^{2}$. The integrand is constant and therefore $A_{i-1,i-2}=-h^{-2}h=-h^{-1}$. A similar reasoning can be applied to $A_{i-1,i}$, which also becomes $-h^{-1}$. The integral of ${\varphi}_i'(x)^2$ gets contributions from two cells, $\Omega^{(i-1)}=[x_{i-1},x_{i}]$ and $\Omega^{(i)}=[x_{i},x_{i+1}]$, but ${\varphi}_i'(x)^2=h^{-2}$ in both cells, and the length of the integration interval is $2h$ so we get $A_{i-1,i-1}=2h^{-1}$.

The right-hand side involves an integral of $2{\varphi}_i(x)$, $i=1,\ldots,N_n-2$, which is just the area under a hat function of height 1 and width $2h$, i.e., equal to $h$. Hence, $b_{i-1}=2h$.

To summarize the linear system, we switch from $i$ to $i+1$ such that we can write

$$A_{i,i-1}=A_{i,i+1}=-h^{-1},\quad A_{i,i}=2h^{-1},\quad b_i = 2h{\thinspace .}$$

The equation system to be solved only involves the unknowns $c_i$ for $i\in{\mathcal{I}_s}$. With our numbering of unknowns and nodes, we have that $c_i$ equals $u(x_{i+1})$. The complete matrix system then takes the following form:

$$\begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} 1 & -1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & -1 & 1 \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} 2h \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ 2h \end{array} \right) \label{fem:deq:1D:ex1:Ab:glob} \tag{66} \end{equation}$$

Comparison with a finite difference discretization¶

A typical row in the matrix system (66) can be written as

$$\begin{equation} -\frac{1}{h}c_{i-1} + \frac{2}{h}c_{i} - \frac{1}{h}c_{i+1} = 2h{\thinspace .} \label{fem:deq:1D:fem:ex1:c} \tag{67} \end{equation}$$

Let us introduce the notation $u_j$ for the value of $u$ at node $j$: $u_j=u(x_{j})$, since we have the interpretation $u(x_{j})=\sum_jc_j{\varphi}(x_{j})=\sum_j c_j\delta_{ij}=c_j$. The unknowns $c_0,\ldots,c_N$ are $u_1,\ldots,u_{N_n-2}$. Shifting $i$ with $i+1$ in (67) and inserting $u_i = c_{i-1}$, we get

$$\begin{equation} -\frac{1}{h}u_{i-1} + \frac{2}{h}u_{i} - \frac{1}{h}u_{i+1} = 2h, \label{fem:deq:1D:fem:ex1} \tag{68} \end{equation}$$

A finite difference discretization of $-u''(x)=2$ by a centered, second-order finite difference approximation $u''(x_i)\approx [D_x D_x u]_i$ with $\Delta x = h$ yields

$$\begin{equation} -\frac{u_{i-1} - 2u_{i} + u_{i+1}}{h^2} = 2, \label{_auto35} \tag{69} \end{equation}$$

which is, in fact, equal to (68) if (68) is divided by $h$. Therefore, the finite difference and the finite element method are equivalent in this simple test problem.

Sometimes a finite element method generates the finite difference equations on a uniform mesh, and sometimes the finite element method generates equations that are different. The differences are modest, but may influence the numerical quality of the solution significantly, especially in time-dependent problems. It depends on the problem at hand whether a finite element discretization is more or less accurate than a corresponding finite difference discretization.

Cellwise computations¶

Software for finite element computations normally employs the cell by cell computational procedure where an element matrix and vector are calculated for each cell and assembled in the global linear system. Let us go through the details of this type of algorithm.

All integrals are mapped to the local reference coordinate system $X\in [-1,1]$. In the present case, the matrix entries contain derivatives with respect to $x$,

$$A_{i-1,j-1}^{(e)}=\int_{\Omega^{(e)}} {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x} = \int_{-1}^1 \frac{d}{dx}{\tilde{\varphi}}_r(X)\frac{d}{dx}{\tilde{\varphi}}_s(X) \frac{h}{2} {\, \mathrm{d}X},$$

where the global degree of freedom $i$ is related to the local degree of freedom $r$ through $i=q(e,r)$. Similarly, $j=q(e,s)$. The local degrees of freedom run as $r,s=0,1$ for a P1 element.

The integral for the element matrix¶

There are simple formulas for the basis functions ${\tilde{\varphi}}_r(X)$ as functions of $X$. However, we now need to find the derivative of ${\tilde{\varphi}}_r(X)$ with respect to $x$. Given

$${\tilde{\varphi}}_0(X)=\frac{1}{2}(1-X),\quad{\tilde{\varphi}}_1(X)=\frac{1}{2}(1+X),$$

we can easily compute $d{\tilde{\varphi}}_r/ dX$:

$$\frac{d{\tilde{\varphi}}_0}{dX} = -\frac{1}{2},\quad \frac{d{\tilde{\varphi}}_1}{dX} = \frac{1}{2}{\thinspace .}$$

From the chain rule,

$$\begin{equation} \frac{d{\tilde{\varphi}}_r}{dx} = \frac{d{\tilde{\varphi}}_r}{dX}\frac{dX}{dx} = \frac{2}{h}\frac{d{\tilde{\varphi}}_r}{dX}{\thinspace .} \label{_auto36} \tag{70} \end{equation}$$

The transformed integral is then

$$A_{i-1,j-1}^{(e)}=\int_{\Omega^{(e)}} {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x} = \int_{-1}^1 \frac{2}{h}\frac{d{\tilde{\varphi}}_r}{dX}\frac{2}{h}\frac{d{\tilde{\varphi}}_s}{dX} \frac{h}{2} {\, \mathrm{d}X} {\thinspace .}$$

The integral for the element vector¶

The right-hand side is transformed according to

$$b_{i-1}^{(e)} = \int_{\Omega^{(e)}} 2{\varphi}_i(x) {\, \mathrm{d}x} = \int_{-1}^12{\tilde{\varphi}}_r(X)\frac{h}{2} {\, \mathrm{d}X},\quad i=q(e,r),\ r=0,1 {\thinspace .}$$

Detailed calculations of the element matrix and vector¶

Specifically for P1 elements we arrive at the following calculations for the element matrix entries:

$$\begin{align*} \tilde A_{0,0}^{(e)} &= \int_{-1}^1\frac{2}{h}\left(-\frac{1}{2}\right) \frac{2}{h}\left(-\frac{1}{2}\right)\frac{h}{2} {\, \mathrm{d}X} = \frac{1}{h}\\ \tilde A_{0,1}^{(e)} &= \int_{-1}^1\frac{2}{h}\left(-\frac{1}{2}\right) \frac{2}{h}\left(\frac{1}{2}\right)\frac{h}{2} {\, \mathrm{d}X} = -\frac{1}{h}\\ \tilde A_{1,0}^{(e)} &= \int_{-1}^1\frac{2}{h}\left(\frac{1}{2}\right) \frac{2}{h}\left(-\frac{1}{2}\right)\frac{h}{2} {\, \mathrm{d}X} = -\frac{1}{h}\\ \tilde A_{1,1}^{(e)} &= \int_{-1}^1\frac{2}{h}\left(\frac{1}{2}\right) \frac{2}{h}\left(\frac{1}{2}\right)\frac{h}{2} {\, \mathrm{d}X} = \frac{1}{h} \end{align*}$$

The element vector entries become

$$\begin{align*} \tilde b_0^{(e)} &= \int_{-1}^12\frac{1}{2}(1-X)\frac{h}{2} {\, \mathrm{d}X} = h\\ \tilde b_1^{(e)} &= \int_{-1}^12\frac{1}{2}(1+X)\frac{h}{2} {\, \mathrm{d}X} = h{\thinspace .} \end{align*}$$

Expressing these entries in matrix and vector notation, we have

$$\begin{equation} \tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{rr} 1 & -1\\ -1 & 1 \end{array}\right),\quad \tilde b^{(e)} = h\left(\begin{array}{c} 1\\ 1 \end{array}\right){\thinspace .} \label{fem:deq:1D:ex1:Ab:elm} \tag{71} \end{equation}$$

Contributions from the first and last cell¶

The first and last cell involve only one unknown and one basis function because of the Dirichlet boundary conditions at the first and last node. The element matrix therefore becomes a $1\times 1$ matrix and there is only one entry in the element vector. On cell 0, only ${\psi}_0={\varphi}_1$ is involved, corresponding to integration with ${\tilde{\varphi}}_1$. On cell $N_e$, only ${\psi}_N={\varphi}_{N_n-2}$ is involved, corresponding to integration with ${\tilde{\varphi}}_0$. We then get the special end-cell contributions

$$\begin{equation} \tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{r} 1 \end{array}\right),\quad \tilde b^{(e)} = h\left(\begin{array}{c} 1 \end{array}\right), \label{fem:deq:1D:ex1:Ab:elm:ends} \tag{72} \end{equation}$$

for $e=0$ and $e=N_e$. In these cells, we have only one degree of freedom, not two as in the interior cells.

Assembly¶

The next step is to assemble the contributions from the various cells. The assembly of an element matrix and vector into the global matrix and right-hand side can be expressed as

$$A_{q(e,r),q(e,s)} = A_{q(e,r),q(e,s)} + \tilde A^{(e)}_{r,s},\quad b_{q(e,r)} = b_{q(e,r)} + \tilde b^{(e)}_{r},\quad$$

for $r$ and $s$ running over all local degrees of freedom in cell $e$.

To make the assembly algorithm more precise, it is convenient to set up Python data structures and a code snippet for carrying out all details of the algorithm. For a mesh of four equal-sized P1 elements and $L=2$ we have

The total number of degrees of freedom is 3, being the function values at the internal 3 nodes where $u$ is unknown. In cell 0 we have global degree of freedom 0, the next cell has $u$ unknown at its two nodes, which become global degrees of freedom 0 and 1, and so forth according to the dof_map list. The mathematical $q(e,r)$ quantity is nothing but the dof_map list.

Assume all element matrices are stored in a list Ae such that Ae[e][i,j] is $\tilde A_{i,j}^{(e)}$. A corresponding list for the element vectors is named be, where be[e][r] is $\tilde b_r^{(e)}$. A Python code snippet illustrates all details of the assembly algorithm:

The general case with N_e P1 elements of length h has

Carrying out the assembly results in a linear system that is identical to (66), which is not surprising, since the procedures is mathematically equivalent to the calculations in the physical domain.

So far, our technique for computing the matrix system have assumed that $u(0)=u(L)=0$. The next section deals with the extension to nonzero Dirichlet conditions.

Boundary conditions: specified nonzero value¶

We have to take special actions to incorporate nonzero Dirichlet conditions, such as $u(L)=D$, into the computational procedures. The present section outlines alternative, yet mathematically equivalent, methods.

General construction of a boundary function¶

In previous sections, we introduced a boundary function $B(x)$ to deal with nonzero Dirichlet boundary conditions for $u$. The construction of such a function is not always trivial, especially not in multiple dimensions. However, a simple and general constructive idea exists when the basis functions have the property

$${\varphi}_i(x_{j}) = \delta_{ij},\quad \delta_{ij} = \left\lbrace\begin{array}{ll} 1, & i=j,\\ 0, & i\neq j, \end{array}\right.$$

where $x_{j}$ is a boundary point. Examples on such functions are the Lagrange interpolating polynomials and finite element functions.

Suppose now that $u$ has Dirichlet boundary conditions at nodes with numbers $i\in{I_b}$. For example, ${I_b} = \{0,N_n-1\}$ in a 1D mesh with node numbering from left to right and Dirichlet conditions at the end nodes $i=0$ and $i=N_n-1$. Let $U_i$ be the corresponding prescribed values of $u(x_{i})$. We can then, in general, use

$$\begin{equation} B(x) = \sum_{j\in{I_b}} U_j{\varphi}_j(x){\thinspace .} \label{_auto37} \tag{73} \end{equation}$$

It is easy to verify that $B(x_{i})= \sum_{j\in{I_b}} U_j{\varphi}_j(x_{i})=U_i$.

The unknown function can then be written as

$$\begin{equation} u(x) = \sum_{j\in{I_b}} U_j{\varphi}_j(x) + \sum_{j\in{\mathcal{I}_s}}c_j{\varphi}_{\nu(j)}, \label{_auto38} \tag{74} \end{equation}$$

where $\nu(j)$ maps unknown number $j$ in the equation system to node $\nu(j)$, ${I_b}$ is the set of indices corresponding to basis functions associated with nodes where Dirichlet conditions apply, and ${\mathcal{I}_s}$ is the set of indices used to number the unknowns from zero to $N$. We can easily show that with this $u$, a Dirichlet condition $u(x_{k})=U_k$ is fulfilled:

$$u(x_{k}) = \sum_{j\in{I_b}} U_j\underbrace{{\varphi}_j(x_k)}_{\neq 0\, \Leftrightarrow\,j=k} + \sum_{j\in{\mathcal{I}_s}} c_j\underbrace{{\varphi}_{\nu(j)}(x_{k})}_{=0,\ k\not\in{\mathcal{I}_s}} = U_k$$

Some examples will further clarify the notation. With a regular left-to-right numbering of nodes in a mesh with P1 elements, and Dirichlet conditions at $x=0$, we use finite element basis functions associated with the nodes $1, 2, \ldots, N_n-1$, implying that $\nu(j)=j+1$, $j=0,\ldots,N$, where $N=N_n-2$. Consider a particular mesh:

The expansion associated with this mesh becomes

$$u(x) = U_0{\varphi}_0(x) + c_0{\varphi}_1(x) + c_1{\varphi}_2(x) + \cdots + c_4{\varphi}_5(x){\thinspace .}$$

Switching to the more standard case of left-to-right numbering and boundary conditions $u(0)=C$, $u(L)=D$, we have $N=N_n-3$ and

$$\begin{align*} u(x) &= C{\varphi}_0 + D{\varphi}_{N_n-1} + \sum_{j\in{\mathcal{I}_s}} c_j{\varphi}_{j+1}\\ &= C{\varphi}_0 + D{\varphi}_{N_n} + c_0{\varphi}_1 + c_1{\varphi}_2 +\cdots + c_N{\varphi}_{N_n-2}{\thinspace .} \end{align*}$$

Finite element meshes in non-trivial 2D and 3D geometries usually leads to an irregular cell and node numbering. Let us therefore take a look at an irregular numbering in 1D:

Say we in this mesh have Dirichlet conditions on the left-most and right-most node, with numbers 3 and 1, respectively. We can number the unknowns at the interior nodes as we want, e.g., from left to right, resulting in $\nu(0)=0$, $\nu(1)=4$, $\nu(2)=5$, $\nu(3)=2$. This gives

$$B(x) = U_3{\varphi}_3(x) + U_1{\varphi}_1(x),$$

and

$$u(x) = B(x) + \sum_{j=0}^3 c_j{\varphi}_{\nu(j)} = U_3{\varphi}_3 + U_1{\varphi}_1 + c_0{\varphi}_0 + c_1{\varphi}_4 + c_2{\varphi}_5 + c_3{\varphi}_2{\thinspace .}$$

The idea of constructing $B$, described here, generalizes almost trivially to 2D and 3D problems: $B=\sum_{j\in{I_b}}U_j{\varphi}_j$, where ${I_b}$ is the index set containing the numbers of all the nodes on the boundaries where Dirichlet values are prescribed.

Example on computing with a finite element-based boundary function¶

Let us see how the model problem $-u''=2$, $u(0)=C$, $u(L)=D$, is affected by a $B(x)$ to incorporate boundary values. Inserting the expression

$$u(x) = B(x) + \sum_{j\in{\mathcal{I}_s}}c_j{\psi}_j(x)$$

in $-(u'',{\psi}_i)=(f,{\psi}_i)$ and integrating by parts results in a linear system with

$$A_{i,j} = \int_0^L {\psi}_i'(x){\psi}_j'(x) {\, \mathrm{d}x},\quad b_i = \int_0^L (f(x){\psi}_i(x) - B'(x){\psi}_i'(x)) {\, \mathrm{d}x}{\thinspace .}$$

We choose ${\psi}_i={\varphi}_{i+1}$, $i=0,\ldots,N=N_n-3$ if the node numbering is from left to right. (Later we also need the assumption that cells too are numbered from left to right.) The boundary function becomes

$$B(x) = C{\varphi}_0(x) + D{\varphi}_{N_n-1}(x){\thinspace .}$$

The expansion for $u(x)$ is

$$u(x) = B(x) + \sum_{j\in{\mathcal{I}_s}} c_j{\varphi}_{j+1}(x){\thinspace .}$$

We can write the matrix and right-hand side entries as

$$\begin{align*} A_{i-1,j-1} &= \int_0^L {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x},\\ b_{i-1} &= \int_0^L (f(x){\varphi}_i'(x) - (C{\varphi}_{0}'(x) + D{\varphi}_{N_n-1}'(x)){\varphi}_i'(x) ){\, \mathrm{d}x}, \end{align*}$$

for $i,j = 1,\ldots,N+1=N_n-2$. Note that we have here used $B'=C{\varphi}_0' + D{\varphi}_{N_n-1}'$.

Computations in physical coordinates¶

Most of the terms in the linear system have already been computed so we concentrate on the new contribution from the boundary function. The integral $C\int_0^L {\varphi}_{0}'(x)){\varphi}_i'(x) {\, \mathrm{d}x}$, associated with the Dirichlet condition in $x=0$, can only get a nonzero contribution from the first cell, $\Omega^{(0)}=[x_{0},x_{1}]$ since ${\varphi}_{0}'(x)=0$ on all other cells. Moreover, ${\varphi}_{0}'(x){\varphi}_i'(x) {\, \mathrm{d}x} \neq 0$ only for $i=0$ and $i=1$ (but node $i=0$ is excluded from the formulation), since ${\varphi}_{i}=0$ on the first cell if $i>1$. With a similar reasoning we realize that $D\int_0^L {\varphi}_{N_n-1}'(x)){\varphi}_i'(x) {\, \mathrm{d}x}$ can only get a nonzero contribution from the last cell.

$$\begin{align*} \int_0^L {\varphi}_{0}'(x){\varphi}_{1}'(x) {\, \mathrm{d}x} &= (-\frac{1}{h})\cdot\frac{1}{h}\cdot h = -\frac{1}{h},\\ \int_0^L {\varphi}_{N_n-1}'(x){\varphi}_{N_n-2}'(x) {\, \mathrm{d}x} &= \frac{1}{h}\cdot(-\frac{1}{h})\cdot h = -\frac{1}{h}{\thinspace .} \end{align*}$$

With these expressions we get

$$b_0 = \int_0^Lf(x){\varphi}_1{\, \mathrm{d}x} - C(-\frac{1}{h}),\quad b_N = \int_0^L f(x){\varphi}_{N_n-2}{\, \mathrm{d}x} - D(-\frac{1}{h}){\thinspace .}$$

Cellwise computations on the reference element¶

As an equivalent alternative, we now turn to cellwise computations. The element matrices and vectors are calculated as in the section Cellwise computations, so we concentrate on the impact of the new term involving $B(x)$. This new term, $B'=C{\varphi}_0' + D{\varphi}_{N_n-1}'$, vanishes on all cells except for $e=0$ and $e=N_e$. Over the first cell ($e=0$) the $B'(x)$ function in local coordinates reads

$$\frac{dB}{dx} = C\frac{2}{h}\frac{d{\tilde{\varphi}}_0}{dX},$$

while over the last cell ($e=N_e$) it looks like

$$\frac{dB}{dx} = D\frac{2}{h}\frac{d{\tilde{\varphi}}_1}{dX}{\thinspace .}$$

For an arbitrary interior cell, we have the formula

$$\tilde b_r^{(e)} = \int_{-1}^1 f(x(X)){\tilde{\varphi}}_r(X)\frac{h}{2}{\, \mathrm{d}X},$$

for an entry in the local element vector. In the first cell, the value at local node 0 is known so only the value at local node 1 is unknown. The associated element vector entry becomes

$$\begin{align*} \tilde b_0^{(1)} &= \int_{-1}^1 \left(f{\tilde{\varphi}}_1 - C\frac{2}{h}\frac{d{\tilde{\varphi}}_0}{dX}\frac{2}{h}\frac{d{\tilde{\varphi}}_1}{dX}\right) \frac{h}{2} {\, \mathrm{d}X} \\ & = \frac{h}{2} 2\int_{-1}^1 {\tilde{\varphi}}_1 {\, \mathrm{d}X} - C\frac{2}{h}(-\frac{1}{2})\frac{2}{h}\frac{1}{2}\frac{h}{2}\cdot 2 = h + C\frac{1}{h}{\thinspace .} \end{align*}$$

The value at local node 1 in the last cell is known so the element vector here is

$$\begin{align*} \tilde b_0^{N_e} &= \int_{-1}^1 \left(f{\tilde{\varphi}}_0 - D\frac{2}{h}\frac{d{\tilde{\varphi}}_1}{dX}\frac{2}{h}\frac{d{\tilde{\varphi}}_0}{dX}\right) \frac{h}{2} {\, \mathrm{d}X}\\ & = \frac{h}{2} 2\int_{-1}^1 {\tilde{\varphi}}_0 {\, \mathrm{d}X} - D\frac{2}{h}\frac{1}{2}\frac{2}{h}(-\frac{1}{2})\frac{h}{2}\cdot 2 = h + D\frac{1}{h}{\thinspace .} \end{align*}$$

The contributions from the $B(x)$ function to the global right-hand side vector becomes $C/h$ for $b_0$ and $D/h$ for $b_N$, exactly as we computed in the physical domain.

Modification of the linear system¶

From an implementational point of view, there is a convenient alternative to adding the $B(x)$ function and using only the basis functions associated with nodes where $u$ is truly unknown. Instead of seeking

$$\begin{equation} u(x) = \sum_{j\in{I_b}} U_j{\varphi}_j(x) + \sum_{j\in{\mathcal{I}_s}}c_j{\varphi}_{\nu(j)}(x), \label{fem:deq:1D:fem:essBC:Bfunc:modsys:utrad} \tag{75} \end{equation}$$

we use the sum over all degrees of freedom, including the known boundary values:

$$\begin{equation} u(x) = \sum_{j\in{\mathcal{I}_s}}c_j{\varphi}_j(x){\thinspace .} \label{fem:deq:1D:fem:essBC:Bfunc:modsys:uall} \tag{76} \end{equation}$$

Note that the collections of unknowns $\left\{ {c}_i \right\}_{i\in{\mathcal{I}_s}}$ in (75) and (76) are different. The index set ${\mathcal{I}_s}=\{0,\ldots,N\}$ always goes to $N$, and the number of unknowns is $N+1$, but in (75) the unknowns correspond to nodes where $u$ is not known, while in (76) the unknowns cover $u$ values at all the nodes. So, if the index set ${I_b}$ contains $N_b$ node numbers where $u$ is prescribed, we have that $N=N_n-N_b$ in (75) and $N=N_n$ in (76).

The idea is to compute the entries in the linear system as if no Dirichlet values are prescribed. Afterwards, we modify the linear system to ensure that the known $c_j$ values are incorporated.

A potential problem arises for the boundary term $[u'v]_0^L$ from the integration by parts: imagining no Dirichlet conditions means that we no longer require $v=0$ at Dirichlet points, and the boundary term is then nonzero at these points. However, when we modify the linear system, we will erase whatever the contribution from $[u'v]_0^L$ should be at the Dirichlet points in the right-hand side of the linear system. We can therefore safely forget $[u'v]_0^L$ at any point where a Dirichlet condition applies.

Computations in the physical system¶

Let us redo the computations in the example in the section General construction of a boundary function. We solve $-u''=2$ with $u(0)=0$ and $u(L)=D$. The expressions for $A_{i,j}$ and $b_i$ are the same, but the numbering is different as the numbering of unknowns and nodes now coincide:

$$A_{i,j} = \int_0^L {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x},\quad b_{i} = \int_0^L f(x){\varphi}_i(x) {\, \mathrm{d}x},$$

for $i,j = 0,\ldots,N=N_n-1$. The integrals involving basis functions corresponding to interior mesh nodes, $i,j=1,\ldots,N_n-2$, are obviously the same as before. We concentrate on the contributions from ${\varphi}_0$ and ${\varphi}_{N_n-1}$:

$$\begin{align*} A_{0,0} &= \int_0^L ({\varphi}_0')^2{\, \mathrm{d}x} = \int_{0}^{x_{1}} = ({\varphi}_0')^2{\, \mathrm{d}x} \frac{1}{h},\\ A_{0,1} &= \int_0^L {\varphi}_0'{\varphi}_1'{\, \mathrm{d}x} = \int_{0}^{x_{1}} {\varphi}_0'{\varphi}_1'{\, \mathrm{d}x} = -\frac{1}{h},\\ A_{N,N} &= \int_0^L ({\varphi}_N')^2{\, \mathrm{d}x} = \int_{x_{N_n-2}}^{x_{N_n-1}} ({\varphi}_N')^2{\, \mathrm{d}x} = \frac{1}{h},\\ A_{N,N-1} &= \int_0^L {\varphi}_N'{\varphi}_{N-1}'{\, \mathrm{d}x} =\int_{x_{N_n-2}}^{x_{N_n-1}} {\varphi}_N'{\varphi}_{N-1}'{\, \mathrm{d}x} = -\frac{1}{h}{\thinspace .} \end{align*}$$

The new terms on the right-hand side are also those involving ${\varphi}_0$ and ${\varphi}_{N_n-1}$:

$$\begin{align*} b_0 &= \int_0^L 2{\varphi}_0(x) {\, \mathrm{d}x} = \int_0^{x_{1}} 2{\varphi}_0(x){\, \mathrm{d}x} = h,\\ b_N &= \int_0^L 2{\varphi}_{N_n-1}{\, \mathrm{d}x} = \int_{x_{N_n-2}}^{x_{N_n-1}} 2{\varphi}_{N_n-1}{\, \mathrm{d}x} = h{\thinspace .} \end{align*}$$

The complete matrix system, involving all degrees of freedom, takes the form

$$\begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} 1 & -1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & -1 & 1 \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} h \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h\\ h \end{array} \right) \label{fem:deq:1D:ex1:Ab:glob2} \tag{77} \end{equation}$$

Incorporation of Dirichlet values can now be done by replacing the first and last equation by the very simple equations $c_0=0$ and $c_N=D$, respectively. Note that the factor $1/h$ in front of the matrix then requires a factor $h$ to be introduce appropriately on the diagonal in the first and last row of the matrix.

$$\begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} h & 0 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 0 & h \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} 0 \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h\\ D \end{array} \right) \label{fem:deq:1D:ex1:Ab:glob3} \tag{78} \end{equation}$$

Note that because we do not require ${\varphi}_i(0)=0$ and ${\varphi}_i(L)=0$, $i\in{\mathcal{I}_s}$, the boundary term $[u'v]_0^L$, in principle, gives contributions $u'(0){\varphi}_0(0)$ to $b_0$ and $u'(L){\varphi}_N(L)$ to $b_N$ ($u'{\varphi}_i$ vanishes for $x=0$ or $x=L$ for $i=1,\ldots,N-1$). Nevertheless, we erase these contributions in $b_0$ and $b_N$ and insert boundary values instead. This argument shows why we can drop computing $[u'v]_0^L$ at Dirichlet nodes when we implement the Dirichlet values by modifying the linear system.

Symmetric modification of the linear system¶

The original matrix system (66) is symmetric, but the modifications in (78) destroy this symmetry. Our described modification will in general destroy an initial symmetry in the matrix system. This is not a particular computational disadvantage for tridiagonal systems arising in 1D problems, but may be more serious in 2D and 3D problems when the systems are large and exploiting symmetry can be important for halving the storage demands and speeding up computations. Methods for solving symmetric matrix are also usually more stable and efficient than those for non-symmetric systems. Therefore, an alternative modification which preserves symmetry is attractive.

One can formulate a general algorithm for incorporating a Dirichlet condition in a symmetric way. Let $c_k$ be a coefficient corresponding to a known value $u(x_{k}) = U_k$. We want to replace equation $k$ in the system by $c_k=U_k$, i.e., insert zeroes in row number $k$ in the coefficient matrix, set 1 on the diagonal, and replace $b_k$ by $U_k$. A symmetry-preserving modification consists in first subtracting column number $k$ in the coefficient matrix, i.e., $A_{i,k}$ for $i\in{\mathcal{I}_s}$, times the boundary value $U_k$, from the right-hand side: $b_i \leftarrow b_i - A_{i,k}U_k$, $i=0,\ldots,N$. Then we put zeroes in both row number $k$ and column number $k$ in the coefficient matrix, and finally set $b_k=U_k$. The steps in algorithmic form becomes

1. $b_i \leftarrow b_i - A_{i,k}U_k$ for $i\in{\mathcal{I}_s}$

2. $A_{i,k} = A_{k,i} = 0$ for $i\in{\mathcal{I}_s}$

3. $A_{k,k}=1$

4. $b_i = U_k$

This modification goes as follows for the specific linear system written out in (77) in the section Modification of the linear system. First we subtract the first column in the coefficient matrix, times the boundary value, from the right-hand side. Because $c_0=0$, this subtraction has no effect. Then we subtract the last column, times the boundary value $D$, from the right-hand side. This action results in $b_{N-1}=2h+D/h$ and $b_N=h-2D/h$. Thereafter, we place zeros in the first and last row and column in the coefficient matrix and 1 on the two corresponding diagonal entries. Finally, we set $b_0=0$ and $b_N=D$. The result becomes

$$\begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} h & 0 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & 0 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 0 & h \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} 0 \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h +D/h\\ D \end{array} \right) \label{fem:deq:1D:ex1:Ab:glob3:symm} \tag{79} \end{equation}$$

Modification of the element matrix and vector¶

The modifications of the global linear system can alternatively be done for the element matrix and vector. Let us perform the associated calculations in the computational example where the element matrix and vector is given by (71). The modifications are needed in cells where one of the degrees of freedom is known. In the present example, this means the first and last cell. We compute the element matrix and vector as if there were no Dirichlet conditions. The boundary term $[u'v]_0^L$ is simply forgotten at nodes that have Dirichlet conditions because the modification of the element vector will anyway erase the contribution from the boundary term. In the first cell, local degree of freedom number 0 is known and the modification becomes

$$\begin{equation} \tilde A^{(0)} = A = \frac{1}{h}\left(\begin{array}{rr} h & 0\\ -1 & 1 \end{array}\right),\quad \tilde b^{(0)} = \left(\begin{array}{c} 0\\ h \end{array}\right){\thinspace .} \label{fem:deq:1D:ex1:Ab:elm:bc:0} \tag{80} \end{equation}$$

In the last cell we set

$$\begin{equation} \tilde A^{(N_e)} = A = \frac{1}{h}\left(\begin{array}{rr} 1 & -1\\ 0 & h \end{array}\right),\quad \tilde b^{(N_e)} = \left(\begin{array}{c} h\\ D \end{array}\right){\thinspace .} \label{fem:deq:1D:ex1:Ab:elm:bc:N} \tag{81} \end{equation}$$

We can also perform the symmetric modification. This operation affects only the last cell with a nonzero Dirichlet condition. The algorithm is the same as for the global linear system, resulting in

$$\begin{equation} \tilde A^{(N_e)} = A = \frac{1}{h}\left(\begin{array}{rr} 1 & 0\\ 0 & h \end{array}\right),\quad \tilde b^{(N_e)} = \left(\begin{array}{c} h + D/h\\ D \end{array}\right){\thinspace .} \label{fem:deq:1D:ex1:Ab:elm:bc:N:symm} \tag{82} \end{equation}$$

The reader is encouraged to assemble the element matrices and vectors and check that the result coincides with the system (79).

Boundary conditions: specified derivative¶

Suppose our model problem $-u''(x)=f(x)$ features the boundary conditions $u'(0)=C$ and $u(L)=D$. As already indicated, the former condition can be incorporated through the boundary term that arises from integration by parts. The details of this method will now be illustrated in the context of finite element basis functions.

The variational formulation¶

Starting with the Galerkin method,

$$\int_0^L(u''(x)+f(x)){\psi}_i(x) {\, \mathrm{d}x} = 0,\quad i\in{\mathcal{I}_s},$$

integrating $u''{\psi}_i$ by parts results in

$$\int_0^Lu'(x)'{\psi}_i'(x) {\, \mathrm{d}x} -(u'(L){\psi}_i(L) - u'(0){\psi}_i(0)) = \int_0^L f(x){\psi}_i(x) {\, \mathrm{d}x}, \quad i\in{\mathcal{I}_s}{\thinspace .}$$

The first boundary term, $u'(L){\psi}_i(L)$, vanishes because $u(L)=D$. The second boundary term, $u'(0){\psi}_i(0)$, can be used to implement the condition $u'(0)=C$, provided ${\psi}_i(0)\neq 0$ for some $i$ (but with finite elements we fortunately have ${\psi}_0(0)=1$). The variational form of the differential equation then becomes

$$\int_0^Lu'(x){\varphi}_i'(x) {\, \mathrm{d}x} + C{\varphi}_i(0) = \int_0^L f(x){\varphi}_i(x) {\, \mathrm{d}x},\quad i\in{\mathcal{I}_s}{\thinspace .} \label{fem:deq:1D:BC:nat:varform} \tag{83}$$

Boundary term vanishes because of the test functions¶

At points where $u$ is known we may require ${\psi}_i$ to vanish. Here, $u(L)=D$ and then ${\psi}_i(L)=0$, $i\in{\mathcal{I}_s}$. Obviously, the boundary term $u'(L){\psi}_i(L)$ then vanishes.

The set of basis functions $\left\{ {{\psi}}_i \right\}_{i\in{\mathcal{I}_s}}$ contains, in this case, all the finite element basis functions on the mesh, except the one that is 1 at $x=L$. The basis function that is left out is used in a boundary function $B(x)$ instead. With a left-to-right numbering, ${\psi}_i = {\varphi}_i$, $i=0,\ldots,N_n-2$, and $B(x)=D{\varphi}_{N_n-1}$:

$$u(x) = D{\varphi}_{N_n-1}(x) + \sum_{j=0}^{N=N_n-2} c_j{\varphi}_j(x){\thinspace .}$$

Inserting this expansion for $u$ in the variational form (83) leads to the linear system

$$\begin{equation} \sum_{j=0}^{N}\left( \int_0^L {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x} \right)c_j = \int_0^L\left(f(x){\varphi}_i(x) -D{\varphi}_{N_n-1}'(x){\varphi}_i'(x)\right) {\, \mathrm{d}x} - C{\varphi}_i(0), \label{fem:deq:1D:natBC} \tag{84} \end{equation}$$

for $i=0,\ldots,N=N_n-2$.

Boundary term vanishes because of linear system modifications¶

We may, as an alternative to the approach in the previous section, use a basis $\left\{ {{\psi}}_i \right\}_{i\in{\mathcal{I}_s}}$ which contains all the finite element functions on the mesh: ${\psi}_i={\varphi}_i$, $i=0,\ldots,N_n-1=N$. In this case, $u'(L){\psi}_i(L)=u'(L){\varphi}_i(L)\neq 0$ for the $i$ corresponding to the boundary node at $x=L$ (where ${\varphi}_i=1$). The number of this node is $i=N_n-1=N$ if a left-to-right numbering of nodes is utilized.

However, even though $u'(L){\varphi}_{N_n-1}(L)\neq 0$, we do not need to compute this term. For $i<N_n-1$ we realize that ${\varphi}_i(L)=0$. The only nonzero contribution to the right-hand side comes from $i=N$ ($b_N$). Without a boundary function we must implement the condition $u(L)=D$ by the equivalent statement $c_N=D$ and modify the linear system accordingly. This modification will erase the last row and replace $b_N$ by another value. Any attempt to compute the boundary term $u'(L){\varphi}_{N_n-1}(L)$ and store it in $b_N$ will be lost. Therefore, we can safely forget about boundary terms corresponding to Dirichlet boundary conditions also when we use the methods from the section Modification of the linear system or the section Symmetric modification of the linear system.

The expansion for $u$ reads

$$u(x) = \sum_{j\in{\mathcal{I}_s}} c_j{\varphi}_j(x){\thinspace .}$$

Insertion in the variational form (83) leads to the linear system

$$\begin{equation} \sum_{j\in{\mathcal{I}_s}}\left( \int_0^L {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x} \right)c_j = \int_0^L\left(f(x){\varphi}_i(x)\right) {\, \mathrm{d}x} - C{\varphi}_i(0),\quad i\in{\mathcal{I}_s} {\thinspace .} \label{fem:deq:1D:natBC2} \tag{85} \end{equation}$$

After having computed the system, we replace the last row by $c_N=D$, either straightforwardly as in the section Modification of the linear system or in a symmetric fashion as in the section Symmetric modification of the linear system. These modifications can also be performed in the element matrix and vector for the right-most cell.

Direct computation of the global linear system¶

We now turn to actual computations with P1 finite elements. The focus is on how the linear system and the element matrices and vectors are modified by the condition $u'(0)=C$.

Consider first the approach where Dirichlet conditions are incorporated by a $B(x)$ function and the known degree of freedom $C_{N_n-1}$ is left out of the linear system (see the section Boundary term vanishes because of the test functions). The relevant formula for the linear system is given by (84). There are three differences compared to the extensively computed case where $u(0)=0$ in the sections Computation in the global physical domain and Cellwise computations. First, because we do not have a Dirichlet condition at the left boundary, we need to extend the linear system (66) with an equation associated with the node $x_{0}=0$. According to the section Modification of the linear system, this extension consists of including $A_{0,0}=1/h$, $A_{0,1}=-1/h$, and $b_0=h$. For $i>0$ we have $A_{i,i}=2/h$, $A_{i-1,i}=A_{i,i+1}=-1/h$. Second, we need to include the extra term $-C{\varphi}_i(0)$ on the right-hand side. Since all ${\varphi}_i(0)=0$ for $i=1,\ldots,N$, this term reduces to $-C{\varphi}_0(0)=-C$ and affects only the first equation ($i=0$). We simply add $-C$ to $b_0$ such that $b_0=h - C$. Third, the boundary term $-\int_0^L D{\varphi}_{N_n-1}(x){\varphi}_i{\, \mathrm{d}x}$ must be computed. Since $i=0,\ldots,N=N_n-2$, this integral can only get a nonzero contribution with $i=N_n-2$ over the last cell. The result becomes $-Dh/6$. The resulting linear system can be summarized in the form

$$\begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} 1 & -1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & -1 & 2 \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} h - C \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ 2h - Dh/6 \end{array} \right){\thinspace .} \label{fem:deq:1D:BC:nat:Aub:system} \tag{86} \end{equation}$$

Next we consider the technique where we modify the linear system to incorporate Dirichlet conditions (see the section Boundary term vanishes because of linear system modifications). Now $N=N_n-1$. The two differences from the case above is that the $-\int_0^LD{\varphi}_{N_n-1}{\varphi}_i{\, \mathrm{d}x}$ term is left out of the right-hand side and an extra last row associated with the node $x_{N_n-1}=L$ where the Dirichlet condition applies is appended to the system. This last row is anyway replaced by the condition $c_N=D$ or this condition can be incorporated in a symmetric fashion. Using the simplest, former approach gives

$$\begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} 1 & -1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & -1 & 2 & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 0 & h \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} h - C \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h \\ D \end{array} \right){\thinspace .} \label{fem:deq:1D:BC:nat:Aub:system:mod} \tag{87} \end{equation}$$

Cellwise computations¶

Now we compute with one element at a time, working in the reference coordinate system $X\in [-1,1]$. We need to see how the $u'(0)=C$ condition affects the element matrix and vector. The extra term $-C{\varphi}_i(0)$ in the variational formulation only affects the element vector in the first cell. On the reference cell, $-C{\varphi}_i(0)$ is transformed to $-C{\tilde{\varphi}}_r(-1)$, where $r$ counts local degrees of freedom. We have ${\tilde{\varphi}}_0(-1)=1$ and ${\tilde{\varphi}}_1(-1)=0$ so we are left with the contribution $-C{\tilde{\varphi}}_0(-1)=-C$ to $\tilde b^{(0)}_0$:

$$\begin{equation} \tilde A^{(0)} = A = \frac{1}{h}\left(\begin{array}{rr} 1 & 1\\ -1 & 1 \end{array}\right),\quad \tilde b^{(0)} = \left(\begin{array}{c} h - C\\ h \end{array}\right){\thinspace .} \label{fem:deq:1D:ex1:Ab:elm:bc:nat} \tag{88} \end{equation}$$

No other element matrices or vectors are affected by the $-C{\varphi}_i(0)$ boundary term.

There are two alternative ways of incorporating the Dirichlet condition. Following the section Boundary term vanishes because of the test functions, we get a $1\times 1$ element matrix in the last cell and an element vector with an extra term containing $D$:

$$\begin{equation} \tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{r} 1 \end{array}\right),\quad \tilde b^{(e)} = h\left(\begin{array}{c} 1 - D/6 \end{array}\right), \label{fem:deq:1D:ex1:Ab:elm:ends2} \tag{89} \end{equation}$$

Alternatively, we include the degree of freedom at the node with $u$ specified. The element matrix and vector must then be modified to constrain the $\tilde c_1 = c_N$ value at local node $r=1$:

$$\begin{equation} \tilde A^{(N_e)} = A = \frac{1}{h}\left(\begin{array}{rr} 1 & 1\\ 0 & h \end{array}\right),\quad \tilde b^{(N_e)} = \left(\begin{array}{c} h\\ D \end{array}\right){\thinspace .} \label{fem:deq:1D:ex1:Ab:elm:bc:nat:mod} \tag{90} \end{equation}$$

Implementation of finite element algorithms¶

At this point, it is sensible to create a program with symbolic calculations to perform all the steps in the computational machinery, both for automating the work and for documenting the complete algorithms. As we have seen, there are quite many details involved with finite element computations and incorporation of boundary conditions. An implementation will also act as a structured summary of all these details.

Extensions of the code for approximation¶

Implementation of the finite element algorithms for differential equations follows closely the algorithm for approximation of functions. The new additional ingredients are

1. other types of integrands (as implied by the variational formulation)

2. additional boundary terms in the variational formulation for Neumann boundary conditions

3. modification of element matrices and vectors due to Dirichlet boundary conditions

Point 1 and 2 can be taken care of by letting the user supply functions defining the integrands and boundary terms on the left- and right-hand side of the equation system:

• Integrand on the left-hand side: ilhs(e, phi, r, s, X, x, h)

• Integrand on the right-hand side: irhs(e, phi, r, X, x, h)

• Boundary term on the left-hand side: blhs (e, phi, r, s, X, x, h)

• Boundary term on the right-hand side: brhs (e, phi, r, s, X, x, h)

Here, phi is a dictionary where phi[q] holds a list of the derivatives of order q of the basis functions with respect to the physical coordinate $x$. The derivatives are available as Python functions of X. For example, phi[0][r](X) means ${\tilde{\varphi}}_r(X)$, and phi[1][s](X, h) means $d{\tilde{\varphi}}_s (X)/dx$ (we refer to the file fe1D.py for details regarding the function basis that computes the phi dictionary). The r and s arguments in the above functions correspond to the index in the integrand contribution from an integration point to $\tilde A^{(e)}_{r,s}$ and $\tilde b^{(e)}_r$. The variables e and h are the current element number and the length of the cell, respectively. Specific examples below will make it clear how to construct these Python functions.

Given a mesh represented by vertices, cells, and dof_map as explained before, we can write a pseudo Python code to list all the steps in the computational algorithm for finite element solution of a differential equation.

Having implemented this function, the user only has supply the ilhs, irhs, blhs, and brhs functions that specify the PDE and boundary conditions. The rest of the implementation forms a generic computational engine. The big finite element packages Diffpack and FEniCS are structured exactly the same way. A sample implementation of ilhs for the 1D Poisson problem is:

which returns $d{\tilde{\varphi}}_i (X)/dx d{\tilde{\varphi}}_j (X)/dx$. Reducing the amount of code the user has to supply and making the code as close as possible to the mathematical formulation makes the software user-friendly and easy to debug.

A complete function finite_element1D_naive for the 1D finite algorithm above, is found in the file fe1D.py. The term "naive" refers to a version of the algorithm where we use a standard dense square matrix as global matrix A. The implementation also has a verbose mode for printing out the element matrices and vectors as they are computed. Below is the complete function without the print statements. You should study in detail since it contains all the steps in the finite element algorithm.

The timing object is a dictionary holding the CPU spent on computing A and the CPU time spent on solving the linear system. (We have left out the timing statements.)

Utilizing a sparse matrix¶

A potential efficiency problem with the finite_element1D_naive function is that it uses dense $(N+1)\times (N+1)$ matrices, while we know that only $2d+1$ diagonals around the main diagonal are different from zero. Switching to a sparse matrix is very easy. Using the DOK (dictionary of keys) format, we declare A as

Assignments or in-place arithmetics are done as for a dense matrix,

    A[i,j] += term
    A[i,j]  = term

but only the index pairs (i,j) we have used in assignments or in-place arithmetics are actually stored. A tailored solution algorithm is needed. The most reliable is sparse Gaussian elimination. SciPy gives access to the UMFPACK algorithm for this purpose:

The declaration of A and the solve statement are the only changes needed in the finite_element1D_naive to utilize sparse matrices. The resulting modification is found in the function finite_element1D.

Application to our model problem¶

Let us demonstrate the finite element software on

$$-u''(x)=f(x),\quad x\in (0,L),\quad u'(0)=C,\ u(L)=D{\thinspace .}$$

This problem can be analytically solved by the model2 function from the section "Simple model problems and their solutions" in previous chapter. Let $f(x)=x^2$. Calling model2(x**2, L, C, D) gives

$$u(x) = D + C(x-L) + \frac{1}{12}(L^4 - x^4)$$

The variational formulation reads

$$(u', v) = (x^2, v) - Cv(0){\thinspace .}$$

The entries in the element matrix and vector, which we need to set up the ilhs, irhs, blhs, and brhs functions, becomes

$$\begin{align*} A^{(e)}_{r,s} &= \int_{-1}^1 \frac{d{\tilde{\varphi}}_r}{dx}\frac{{\tilde{\varphi}}_s}{dx}(\det J{\, \mathrm{d}X}),\\ b^{(e)} &= \int_{-1}^1 x^2{\tilde{\varphi}}_r\det J{\, \mathrm{d}X} - C{\tilde{\varphi}}_r(-1)I(e,0), \end{align*}$$

where $I(e)$ is an indicator function: $I(e,q)=1$ if $e=q$, otherwise $I(e)=0$. We use this indicator function to formulate that the boundary term $Cv(0)$, which in the local element coordinate system becomes $C{\tilde{\varphi}}_r(-1)$, is only included for the element $e=0$.

The functions for specifying the element matrix and vector entries must contain the integrand, but without the $\det J{\, \mathrm{d}X}$ term as this term is taken care of by the quadrature loop, and the derivatives $d{\tilde{\varphi}}_r(X)/dx$ with respect to the physical $x$ coordinates are contained in phi[1][r](X), computed by the function basis.

We can then make the call to finite_element1D_naive or finite_element1D to solve the problem with two P1 elements:

It remains to plot the solution (with high resolution in each element). To this end, we use the u_glob function imported from fe1D, which imports it from fe_approx1D_numit (the u_glob function in fe_approx1D.py works with elements and nodes, while u_glob in fe_approx1D_numint works with cells, vertices, and dof_map):

The result is shown in Figure. We see that the solution using P1 elements is exact at the nodes, but feature considerable discrepancy between the nodes. Exercise 10: Investigate exact finite element solutions asks you to explore this problem further using other $m$ and $d$ values.

Finite element and exact solution using two cells.

Variational formulations in 2D and 3D¶

The major difference between deriving variational formulations in 2D and 3D compared to 1D is the rule for integrating by parts. The cells have shapes different from an interval, so basis functions look a bit different, and there is a technical difference in actually calculating the integrals over cells. Otherwise, going to 2D and 3D is not a big step from 1D. All the fundamental ideas still apply.

Integration by parts¶

A typical second-order term in a PDE may be written in dimension-independent notation as

$$\nabla^2 u \quad\hbox{or}\quad \nabla\cdot\left( {\alpha}(\boldsymbol{x})\nabla u\right) {\thinspace .}$$

The explicit forms in a 2D problem become

$$\nabla^2 u = \nabla\cdot\nabla u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2},$$

and

$$\nabla\cdot\left( a(\boldsymbol{x})\nabla u\right) = \frac{\partial}{\partial x}\left( {\alpha}(x,y)\frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left( {\alpha}(x,y)\frac{\partial u}{\partial y}\right) {\thinspace .}$$

We shall continue with the latter operator as the former arises from just setting ${\alpha} =1$.

The integration by parts formula for $\int\nabla\cdot({\alpha}\nabla)$.

The general rule for integrating by parts is often referred to as Green's first identity:

$$\begin{equation} -\int_{\Omega} \nabla\cdot ({\alpha}(\boldsymbol{x})\nabla u) v{\, \mathrm{d}x} = \int_{\Omega} {\alpha}(\boldsymbol{x})\nabla u\cdot\nabla v {\, \mathrm{d}x} - \int_{\partial\Omega} a\frac{\partial u}{\partial n} v {\, \mathrm{d}s}, \label{fem:deq:2D:int:by:parts} \tag{91} \end{equation}$$

where $\partial\Omega$ is the boundary of $\Omega$ and $\partial u/\partial n = \boldsymbol{n}\cdot\nabla u$ is the derivative of $u$ in the outward normal direction, $\boldsymbol{n}$ being an outward unit normal to $\partial\Omega$. The integrals $\int_\Omega (){\, \mathrm{d}x}$ are area integrals in 2D and volume integrals in 3D, while $\int_{\partial\Omega} (){\, \mathrm{d}s}$ is a line integral in 2D and a surface integral in 3D.

It will be convenient to divide the boundary into two parts:

• $\partial\Omega_N$, where we have Neumann conditions $-a\frac{\partial u}{\partial n} = g$, and

• $\partial\Omega_D$, where we have Dirichlet conditions $u = u_0$.

The test functions $v$ are (as usual) required to vanish on $\partial\Omega_D$.

Example on a multi-dimensional variational problem¶

Here is a quite general, stationary, linear PDE arising in many problems:

$$\begin{equation} \label{varform:conv:diff:pde:pre} \tag{92} \boldsymbol{v}\cdot\nabla u + \beta u = \nabla\cdot\left( {\alpha}\nabla u\right) + f, \quad\boldsymbol{x}\in\Omega, \end{equation}$$

$$\begin{equation} \label{varform:conv:diff:bc1:pre} \tag{93} u = u_0,\quad\boldsymbol{x}\in\partial\Omega_D, \end{equation}$$

$$\begin{equation} \label{varform:conv:diff:bc2:pre} \tag{94} -{\alpha}\frac{\partial u}{\partial n} = g,\quad\boldsymbol{x}\in\partial\Omega_N {\thinspace .} \end{equation}$$

The vector field $\boldsymbol{v}$ and the scalar functions $a$, $\alpha$, $f$, $u_0$, and $g$ may vary with the spatial coordinate $\boldsymbol{x}$ and must be known.

Such a second-order PDE needs exactly one boundary condition at each point of the boundary, so $\partial\Omega_N\cup\partial\Omega_D$ must be the complete boundary $\partial\Omega$.

Assume that the boundary function $u_0(\boldsymbol{x})$ is defined for all $\boldsymbol{x}\in\Omega$. The unknown function can then be expanded as

$$u = B + \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j,\quad B = u_0 {\thinspace .}$$

As long as any ${\psi}_j=0$ on $\partial\Omega_D$, we realize that $u=u_0$ on $\partial\Omega_D$.

The variational formula is obtained from Galerkin's method, which technically means multiplying the PDE by a test function $v$ and integrating over $\Omega$:

$$\int_{\Omega} (\boldsymbol{v}\cdot\nabla u + \beta u)v{\, \mathrm{d}x} = \int_{\Omega} \nabla\cdot\left( {\alpha}\nabla u\right){\, \mathrm{d}x} + \int_{\Omega}fv {\, \mathrm{d}x} {\thinspace .}$$

The second-order term is integrated by parts, according to the formula (91):

$$\int_{\Omega} \nabla\cdot\left( {\alpha}\nabla u\right)v {\, \mathrm{d}x} = -\int_{\Omega} {\alpha}\nabla u\cdot\nabla v{\, \mathrm{d}x} + \int_{\partial\Omega} {\alpha}\frac{\partial u}{\partial n} v{\, \mathrm{d}s} {\thinspace .}$$

Galerkin's method therefore leads to

$$\int_{\Omega} (\boldsymbol{v}\cdot\nabla u + \beta u)v{\, \mathrm{d}x} = -\int_{\Omega} {\alpha}\nabla u\cdot\nabla v{\, \mathrm{d}x} + \int_{\partial\Omega} {\alpha}\frac{\partial u}{\partial n} v{\, \mathrm{d}s} + \int_{\Omega} fv {\, \mathrm{d}x} {\thinspace .}$$

The boundary term can be developed further by noticing that $v\neq 0$ only on $\partial\Omega_N$,

$$\int_{\partial\Omega} {\alpha}\frac{\partial u}{\partial n} v{\, \mathrm{d}s} = \int_{\partial\Omega_N} {\alpha}\frac{\partial u}{\partial n} v{\, \mathrm{d}s},$$

and that on $\partial\Omega_N$, we have the condition $a\frac{\partial u}{\partial n}=-g$, so the term becomes

$$-\int_{\partial\Omega_N} gv{\, \mathrm{d}s}{\thinspace .}$$

The final variational form is then

$$\int_{\Omega} (\boldsymbol{v}\cdot\nabla u + \beta u)v{\, \mathrm{d}x} = -\int_{\Omega} {\alpha}\nabla u\cdot\nabla v {\, \mathrm{d}x} - \int_{\partial\Omega_N} g v{\, \mathrm{d}s} + \int_{\Omega} fv {\, \mathrm{d}x} {\thinspace .}$$

Instead of using the integral signs, we may use the inner product notation:

$$(\boldsymbol{v}\cdot\nabla u, v) + (\beta u,v) = - ({\alpha}\nabla u,\nabla v) - (g,v)_{N} + (f,v) {\thinspace .}$$

The subscript $\,{}_N$ in $(g,v)_{N}$ is a notation for a line or surface integral over $\partial\Omega_N$, while $(\cdot,\cdot)$ is the area/volume integral over $\Omega$.

We can derive explicit expressions for the linear system for $\left\{ {c}_j \right\}_{j\in{\mathcal{I}_s}}$ that arises from the variational formulation. Inserting the $u$ expansion results in

$$\begin{align*} \sum_{j\in{\mathcal{I}_s}} ((\boldsymbol{v}\cdot\nabla {\psi}_j, {\psi}_i) &+ (\beta {\psi}_j ,{\psi}_i) + ({\alpha}\nabla {\psi}_j,\nabla {\psi}_i))c_j = \\ & (g,{\psi}_i)_{N} + (f,{\psi}_i) - (\boldsymbol{v}\cdot\nabla u_0, {\psi}_i) + (\beta u_0 ,{\psi}_i) + ({\alpha}\nabla u_0,\nabla {\psi}_i) {\thinspace .} \end{align*}$$

This is a linear system with matrix entries

$$A_{i,j} = (\boldsymbol{v}\cdot\nabla {\psi}_j, {\psi}_i) + (\beta {\psi}_j ,{\psi}_i) + ({\alpha}\nabla {\psi}_j,\nabla {\psi}_i)$$

and right-hand side entries

$$b_i = (g,{\psi}_i)_{N} + (f,{\psi}_i) - (\boldsymbol{v}\cdot\nabla u_0, {\psi}_i) + (\beta u_0 ,{\psi}_i) + ({\alpha}\nabla u_0,\nabla {\psi}_i),$$

for $i,j\in{\mathcal{I}_s}$.

In the finite element method, we usually express $u_0$ in terms of basis functions and restrict $i$ and $j$ to run over the degrees of freedom that are not prescribed as Dirichlet conditions. However, we can also keep all the $\left\{ {c}_j \right\}_{j\in{\mathcal{I}_s}}$ as unknowns, drop the $u_0$ in the expansion for $u$, and incorporate all the known $c_j$ values in the linear system. This has been explained in detail in the 1D case, and the technique is the same for 2D and 3D problems.

Transformation to a reference cell in 2D and 3D¶

The real power of the finite element method first becomes evident when we want to solve partial differential equations posed on two- and three-dimensional domains of non-trivial geometric shape. As in 1D, the domain $\Omega$ is divided into $N_e$ non-overlapping cells. The elements have simple shapes: triangles and quadrilaterals are popular in 2D, while tetrahedra and box-shapes elements dominate in 3D. The finite element basis functions ${\varphi}_i$ are, as in 1D, polynomials over each cell. The integrals in the variational formulation are, as in 1D, split into contributions from each cell, and these contributions are calculated by mapping a physical cell, expressed in physical coordinates $\boldsymbol{x}$, to a reference cell in a local coordinate system $\boldsymbol{X}$. This mapping will now be explained in detail.

We consider an integral of the type

$$\begin{equation} \int_{{\Omega}^{(e)}} {\alpha}(\boldsymbol{x})\nabla{\varphi}_i\cdot\nabla{\varphi}_j{\, \mathrm{d}x}, \label{_auto39} \tag{95} \end{equation}$$

where the ${\varphi}_i$ functions are finite element basis functions in 2D or 3D, defined in the physical domain. Suppose we want to calculate this integral over a reference cell, denoted by $\tilde\Omega^r$, in a coordinate system with coordinates $\boldsymbol{X} = (X_0, X_1)$ (2D) or $\boldsymbol{X} = (X_0, X_1, X_2)$ (3D). The mapping between a point $\boldsymbol{X}$ in the reference coordinate system and the corresponding point $\boldsymbol{x}$ in the physical coordinate system is given by a vector relation $\boldsymbol{x}(\boldsymbol{X})$. The corresponding Jacobian, $J$, of this mapping has entries

$$J_{i,j}=\frac{\partial x_j}{\partial X_i}{\thinspace .}$$

The change of variables requires ${\, \mathrm{d}x}$ to be replaced by $\det J{\, \mathrm{d}X}$. The derivatives in the $\nabla$ operator in the variational form are with respect to $\boldsymbol{x}$, which we may denote by $\nabla_{\boldsymbol{x}}$. The ${\varphi}_i(\boldsymbol{x})$ functions in the integral are replaced by local basis functions ${\tilde{\varphi}}_r(\boldsymbol{X})$ so the integral features $\nabla_{\boldsymbol{x}}{\tilde{\varphi}}_r(\boldsymbol{X})$. We readily have $\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r(\boldsymbol{X})$ from formulas for the basis functions in the reference cell, but the desired quantity $\nabla_{\boldsymbol{x}}{\tilde{\varphi}}_r(\boldsymbol{X})$ requires some efforts to compute. All the details are provided below.

Let $i=q(e,r)$ and consider two space dimensions. By the chain rule,

$$\frac{\partial {\tilde{\varphi}}_r}{\partial X} = \frac{\partial {\varphi}_i}{\partial X} = \frac{\partial {\varphi}_i}{\partial x}\frac{\partial x}{\partial X} + \frac{\partial {\varphi}_i}{\partial y}\frac{\partial y}{\partial X},$$

and

$$\frac{\partial {\tilde{\varphi}}_r}{\partial Y} = \frac{\partial {\varphi}_i}{\partial Y} = \frac{\partial {\varphi}_i}{\partial x}\frac{\partial x}{\partial Y} + \frac{\partial {\varphi}_i}{\partial y}\frac{\partial y}{\partial Y} {\thinspace .}$$

We can write these two equations as a vector equation

$$\left[\begin{array}{c} \frac{\partial {\tilde{\varphi}}_r}{\partial X}\\ \frac{\partial {\tilde{\varphi}}_r}{\partial Y} \end{array}\right] = \left[\begin{array}{cc} \frac{\partial x}{\partial X} & \frac{\partial y}{\partial X}\\ \frac{\partial x}{\partial Y} & \frac{\partial y}{\partial Y} \end{array}\right] \left[\begin{array}{c} \frac{\partial {\varphi}_i}{\partial x}\\ \frac{\partial {\varphi}_i}{\partial y} \end{array}\right]$$

Identifying

$$\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r = \left[\begin{array}{c} \frac{\partial {\tilde{\varphi}}_r}{\partial X}\\ \frac{\partial {\tilde{\varphi}}_r}{\partial Y} \end{array}\right], \quad J = \left[\begin{array}{cc} \frac{\partial x}{\partial X} & \frac{\partial y}{\partial X}\\ \frac{\partial x}{\partial Y} & \frac{\partial y}{\partial Y} \end{array}\right], \quad \nabla_{\boldsymbol{x}}{\varphi}_r = \left[\begin{array}{c} \frac{\partial {\varphi}_i}{\partial x}\\ \frac{\partial {\varphi}_i}{\partial y} \end{array}\right],$$

we have the relation

$$\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r = J\cdot\nabla_{\boldsymbol{x}}{\varphi}_i,$$

which we can solve with respect to $\nabla_{\boldsymbol{x}}{\varphi}_i$:

$$\begin{equation} \nabla_{\boldsymbol{x}}{\varphi}_i = J^{-1}\cdot\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r{\thinspace .} \label{_auto40} \tag{96} \end{equation}$$

On the reference cell, ${\varphi}_i(\boldsymbol{x}) = {\tilde{\varphi}}_r(\boldsymbol{X})$, so

$$\begin{equation} \nabla_{\boldsymbol{x}}{\tilde{\varphi}}_r(\boldsymbol{X}) = J^{-1}(\boldsymbol{X})\cdot\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r(\boldsymbol{X}){\thinspace .} \label{_auto41} \tag{97} \end{equation}$$

This means that we have the following transformation of the integral in the physical domain to its counterpart over the reference cell:

$$\begin{equation} \int_{\Omega^{(e)}} {\alpha}(\boldsymbol{x})\nabla_{\boldsymbol{x}}{\varphi}_i\cdot\nabla_{\boldsymbol{x}}{\varphi}_j{\, \mathrm{d}x} = \int_{\tilde\Omega^r} {\alpha}(\boldsymbol{x}(\boldsymbol{X}))(J^{-1}\cdot\nabla_{\boldsymbol{X}}{\tilde{\varphi}}_r)\cdot (J^{-1}\cdot\nabla{\tilde{\varphi}}_s)\det J{\, \mathrm{d}X} \label{_auto42} \tag{98} \end{equation}$$

Numerical integration¶

Integrals are normally computed by numerical integration rules. For multi-dimensional cells, various families of rules exist. All of them are similar to what is shown in 1D: $\int f {\, \mathrm{d}x}\approx \sum_jw_if(\boldsymbol{x}_j)$, where $w_j$ are weights and $\boldsymbol{x}_j$ are corresponding points.

The file numint.py contains the functions quadrature_for_triangles(n) and quadrature_for_tetrahedra(n), which returns lists of points and weights corresponding to integration rules with n points over the reference triangle with vertices $(0,0)$, $(1,0)$, $(0,1)$, and the reference tetrahedron with vertices $(0,0,0)$, $(1,0,0)$, $(0,1,0)$, $(0,0,1)$, respectively. For example, the first two rules for integration over a triangle have 1 and 3 points:

Rules with 1, 3, 4, and 7 points over the triangle will exactly integrate polynomials of degree 1, 2, 3, and 4, respectively. In 3D, rules with 1, 4, 5, and 11 points over the tetrahedron will exactly integrate polynomials of degree 1, 2, 3, and 4, respectively.

Convenient formulas for P1 elements in 2D¶

We shall now provide some formulas for piecewise linear ${\varphi}_i$ functions and their integrals in the physical coordinate system. These formulas make it convenient to compute with P1 elements without the need to work in the reference coordinate system and deal with mappings and Jacobians. A lot of computational and algorithmic details are hidden by this approach.

Let $\Omega^{(e)}$ be cell number $e$, and let the three vertices have global vertex numbers $I$, $J$, and $K$. The corresponding coordinates are $(x_{I},y_{I})$, $(x_{J},y_{J})$, and $(x_{K},y_{K})$. The basis function ${\varphi}_I$ over $\Omega^{(e)}$ have the explicit formula

$$\begin{equation} {\varphi}_I (x,y) = \frac{1}{2}\Delta \left( \alpha_I + \beta_Ix + \gamma_Iy\right), \label{fem:approx:fe:2D:phi:I} \tag{99} \end{equation}$$

where

$$\begin{equation} \alpha_I = x_{J}y_{K} - x_{K}y_{J}, \label{fem:approx:fe:2D:phi:alpha:I} \tag{100} \end{equation}$$

$$\begin{equation} \beta_I = y_{J} - y_{K}, \label{fem:approx:fe:2D:phi:beta:I} \tag{101} \end{equation}$$

$$\begin{equation} \gamma_I = x_{K} - x_{J}, \label{fem:approx:fe:2D:phi:gamma:I} \tag{102}, \end{equation}$$

and

$$\begin{equation} 2\Delta = \det\left(\begin{array}{rrr} 1 & x_{I} & y_{I} \\ 1 & x_{J} & y_{J} \\ 1 & x_{K} & y_{K} \end{array}\right) {\thinspace .} \label{fem:approx:fe:2D:phi:Delta} \tag{103} \end{equation}$$

The quantity $\Delta$ is the area of the cell.

The following formula is often convenient when computing element matrices and vectors:

$$\begin{equation} \int_{\Omega^{(e)}} {\varphi}_I^{p}{\varphi}_J^{q}{\varphi}_K^{r} dx dy = {p!q!r!\over (p+q+r+2)!}2\Delta \label{fem:approx:fe:2D:phi:integral} \tag{104} {\thinspace .} \end{equation}$$

(Note that the $q$ in this formula is not to be mixed with the $q(e,r)$ mapping of degrees of freedom.)

As an example, the element matrix entry $\int_{\Omega^{(e)}} {\varphi}_I{\varphi}_J{\, \mathrm{d}x}$ can be computed by setting $p=q=1$ and $r=0$, when $I\neq J$, yielding $\Delta/12$, and $p=2$ and $q=r=0$, when $I=J$, resulting in $\Delta/6$. We collect these numbers in a local element matrix:

$$\frac{\Delta}{12} \left[\begin{array}{ccc} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{array}\right]$$

The common element matrix entry $\int_{\Omega^{(e)}} \nabla{\varphi}_I\cdot\nabla{\varphi}_J{\, \mathrm{d}x}$, arising from a Laplace term $\nabla^2u$, can also easily be computed by the formulas above. We have

$$\nabla{\varphi}_I\cdot\nabla{\varphi}_J = \frac{\Delta^2}{4}(\beta_I\beta_J + \gamma_I\gamma_J) = \hbox{const},$$

so that the element matrix entry becomes $\frac{1}{4}\Delta^3(\beta_I\beta_J + \gamma_I\gamma_J)$.

From an implementational point of view, one will work with local vertex numbers $r=0,1,2$, parameterize the coefficients in the basis functions by $r$, and look up vertex coordinates through $q(e,r)$.

Similar formulas exist for integration of P1 elements in 3D.

Implementation in 2D and 3D via FEniCS¶

From a principle of view, we have seen that variational forms of the type: find $a(u,v)=L\ \forall v\in V$ (and even general nonlinear problems $F(u;v)=0$), can apply the computational machinery of introduced for the approximation problem $u=f$. We actually need two extensions only:

1. specify Dirichlet boundary conditions as part of $V$

2. incorporate Neumann flux boundary conditions in the variational form

The algorithms are all the same in any space dimension, we only need to choose the element type and associated integration rule. Once we know how to compute things in 1D, and made the computer code sufficiently flexible, the method and code should work for any variational form in any number of space dimensions! This fact is exactly the idea behind the FEniCS finite element software.

Therefore, if we know how to set up an approximation problem in any dimension in FEniCS, and know how to derive variational forms in higher dimensions, we are (in principle!) very close to solving a PDE problem in FEniCS. We shall now solve a quite general 1D/2D/3D Poisson problem in FEniCS. There is much more to FEniCS than what is shown in this example, but it illustrates the fact that when we go beyond 1D, there exists software which leverage the full power of the finite element method as a method for solving any problem on any mesh in any number of space dimensions.

Mathematical problem¶

The following model describes the pressure $u$ in the flow around a bore hole of radius $a$ in a porous medium. If the hole is long in the vertical direction (z-direction) then it is natural to assume that the vertical changes are small and $u_z\approx\mbox{constant}$. Therefore, we can model it by a 2D domain in the cross section.

$$\begin{equation} \nabla\cdot \left( {\alpha}\nabla u\right) = 0, \quad a < ||\boldsymbol{x}|| < b, \label{fem:varform:fenics:problem:pde} \tag{105} \end{equation}$$

$$\begin{equation} u(\boldsymbol{x}) = U_a, \quad ||\boldsymbol{x}|| = a, \label{fem:varform:fenics:problem:ua} \tag{106} \end{equation}$$

$$\begin{equation} u(\boldsymbol{x}) = U_b \quad ||\boldsymbol{x}|| = b{\thinspace .} \label{fem:varform:fenics:problem:ub} \tag{107} \end{equation}$$

That is, we have a hollow circular 2D domain with inner radius $a$ and outer radius $b$. The pressure is known on these two boundaries, so this is a pure Dirichlet problem.

Symmetry¶

The first thing we should observe is that the problem is radially symmetric, so we can change to polar coordinates and obtain a 1D problem in the radial direction:

$$(r{\alpha} u')' = 0,\quad u(a)=U_a, u(b)=U_b{\thinspace .}$$

This is not very exciting beyond being able to find an analytical solution and compute the true error of a finite element approximation.

However, many software packages solve problems in Cartesian coordinates, and FEniCS basically do this, so we want to take advantage of symmetry in Cartesian coordinates and reformulate the problem in a smaller domain.

Looking at the domain as a cake with a hole, any piece of the cake will be a candidate for a reduced-size domain. The solution is symmetric about any line $\theta = \hbox{const}$ in polar coordinates, so at such lines we have the symmetry boundary condition $\partial u/\partial n=0$, i.e., a homogeneous Neumann condition. In Figure we have plotted a possible mesh of cells as triangles, here with dense refinement toward the bore hole, because we know the solution will decay most rapidly toward the origin. This mesh is a piece of the cake with four sides: Dirichlet conditions on the inner and outer boundary, named $\Gamma_{D_a}$ and $\Gamma_{D_b}$, and $\partial u/\partial n=0$ on the two other sides, named $\Gamma_N$. In this particular example, the arc of the piece of the cake is 45 degrees, but any value of the arc will work.

Mesh of a hollow cylinder, with refinement and utilizing symmetry.

The boundary problem can then be expressed as

$$\begin{equation} \nabla\cdot \left( {\alpha}\nabla u\right) = 0, \quad \boldsymbol{x}\in\Omega, \label{fem:varform:fenics:problem:pde2} \tag{108} \end{equation}$$

$$\begin{equation} u(\boldsymbol{x}) = U_a, \quad \boldsymbol{x}\in\Gamma_{D_a}, \label{fem:varform:fenics:problem:Ga} \tag{109} \end{equation}$$

$$\begin{equation} u(\boldsymbol{x}) = U_b, \quad \boldsymbol{x}\in\Gamma_{D_b}, \label{fem:varform:fenics:problem:Gb} \tag{110} \end{equation}$$

$$\begin{equation} \frac{\partial u}{\partial n} =0,\quad \boldsymbol{x}\in\Gamma_N{\thinspace .} \label{fem:varform:fenics:problem:GN} \tag{111} \end{equation}$$

Variational formulation¶

To obtain the variational formulation, we multiply the PDE by a test function $v$ and integrate the second-order derivatives by part:

$$\begin{align*} \int_\Omega \nabla\cdot ({\alpha}\nabla u) v {\, \mathrm{d}x} &= -\int_\Omega {\alpha}\nabla u\cdot\nabla v{\, \mathrm{d}x} + \int_{\Gamma_N}{\alpha} \frac{\partial u}{\partial n}v{\, \mathrm{d}s}\\ &= -\int_\Omega {\alpha}\nabla u\cdot\nabla v{\, \mathrm{d}x},\quad \forall v\in V\\ \end{align*}$$

We are left with a problem of the form: find $u$ such that $a(u,v)=L(v)\ \forall v\in V$, with

$$\begin{equation} a(u,v) = \int_\Omega {\alpha}\nabla u\cdot\nabla v{\, \mathrm{d}x}, \label{_auto43} \tag{112} \end{equation}$$

$$\begin{equation} L(v) = \int_\Omega 0v{\, \mathrm{d}x} {\thinspace .} \label{_auto44} \tag{113} \end{equation}$$

We write the integrand as $0v{\, \mathrm{d}x}$ even though $L=0$, because it is necessary in FEniCS to specify $L$ as a linear form (i.e., a test function and some form of integration need to be present) and not the number zero. The Dirichlet conditions make a nonzero solution.

The FEniCS solver¶

We suppose that we have a function make_mesh that can make the mesh for us. More details about this function will be provided later. A next step is then to define proper Dirichlet conditions. This might seem a bit complicated, but we introduce markers at the boundary for marking the Dirichlet boundaries. The inner boundary has marker 1 and the outer has marker 2. In this way, we can recognize the nodes that are on the boundary. It is usually a part of the mesh making process to compute both the mesh and its markers, so make_mesh returns a Mesh object as mesh and a MeshFunction object markers. Setting Dirichlet conditions in the solver is then a matter of introducing DirichletBC objects, one for each part of the boundary marked by markers, and then we collect all such Dirichlet conditions in a list that is used by the assembly process to incorporate the Dirichlet conditions in the linear system. The code goes like this:

Here, u_a and u_b are constants (floats) set by the user. In general, anything that can be evaluated pointwise can be used, such as Expression, Function, and Constant. The next step is to define the variational problem and solve it:

In order to avoid L=0 (L equal to the float zero), we have to tell FEniCS that is a linear form, so zero must be specified as Constant(0).

Note that everything is the same as for the approximation problem in the section fe:approx:fenics, except for the Dirichlet conditions and the formulas for a and L. FEniCS has, of course, access to very efficient solution methods, so we could add arguments to the solve call to apply state-of-the-art iterative methods and preconditioners for large-scale problems. However, for this little 2D case a standard sparse Gaussian elimination, as implied by solve(a = L, u, bcs) is a sufficient approach.

Finally, we can save the solution to file for using professional visualization software and, if desired, add a quick plotting using the built-in FEniCS tool plot:

(The u.rename call is just for getting a more readable title in the plot.)

The above statements are collected in a function solver in the file borehole_fenics.py:

Be careful with name clashes!

It is easy when coding mathematics to use variable names that correspond to one-letter names in the mathematics. For example, in the mathematics of this problem there are two $a$ variables: the radius of the inner boundary and the bilinear form in the variational formulation. Using a for the inner boundary in solver does not work: it is quickly overwritten by the bilinear form. We therefore have to introduce x_a. Long variable names are to be preferred for safe programming, though short names corresponding to the mathematics are often nicer.

Making the mesh¶

The hardest part of a finite element problem is very often to make the mesh. This is particularly the case in large industrial projects, but also often academic projects quickly lead to time-consuming work with constructing finite element meshes. In the present example we create the mesh for the symmetric problem by deforming an originally rectangular mesh. The rectangular mesh is made by the FEniCS object RectangleMesh on $[a,b]\times [0,1]$. Therefore, we stretch the mesh towards the left before we bend the rectangle onto to "a piece of cake". Figure shows an example on the resulting mesh. The stretching gives us refinement in the radial direction because we expect the variations to be quite large in this direction, but uniform in $\theta$ direction.

We first make the rectangle and set boundary markers here for the inner and outer boundaries (since these are particularly simple: $x=a$ and $x=b$). Here is how we make the rectangular mesh from lower left corner $(a,0)$ to upper left corner $(b,1)$ with nr quadrilateral cells in $x$ direction (later to become the radial direction) and nt quadrilateral cells in the $y$ direction:

Each quadrilateral cell is divided into two triangles with right or left going diagonals, or four triangles using both diagonals. These choices of producing triangles from rectangles are named right, left, and crossed. Recall that FEniCS can only work with cells of triangular shape only and where the sides are straight. This means that we need a good resolution in $\theta$ direction to represent a circular boundary. With isoparametric elements, it is easier to get a higher-order polynomial approximation of curved boundaries.

We must then mark the boundaries for boundary conditions. Since we do not need to do anything with the homogeneous Neumann conditions, we can just mark the inner and outer boundary of the hole cylinder. This is very easy to do as long as the mesh is a rectangle since then the specifications of the boundaries are $x=a$ and $x=b$. The relevant FEniCS code requires the user to define a subclass of SubDomain and implement a function inside for indicating whether a given point x is on the desired boundary or not:

Finite element mesh for a porous medium outside a bore hole (hollow cylinder).

With the instances inner and outer we fill a marker object, called MeshFunction in FEniCS. For this purpose we must introduce our own convention of numbering boundaries: here we use 1 for all points on the inner boundary and 2 for the outer boundary, while all other points are marked by 0. The solver applies the markers object to set the right Dirichlet boundary conditions.

The next step is to deform the mesh. Given coordinates $x$, we can map these onto a stretched coordinate $\bar x$ by

$$\begin{equation} \bar x = a + (b-a)\left({x-a\over b-a}\right)^s, \label{_auto45} \tag{114} \end{equation}$$

where $s$ is a parameter that controls the amount of stretching. The formula above gives a stretching towards $x=a$, while the next one stretches the coordinates towards $x=b$:

$$\begin{equation} \bar x = a + (b-a)\left({x-a\over b-a}\right)^{1/s}{\thinspace .} \label{_auto46} \tag{115} \end{equation}$$

The code shown later shows the details of mapping coordinates in a FEniCS mesh object.

The final step is to map the rectangle to a part of a hollow cylinder. Mathematically, a point $(x,y)$ in the rectangle is mapped onto $\bar x, \bar y)$ in our final geometry:

$$\hat x = \bar x\cos (\Theta \bar y),\quad \hat y = \bar x\sin (\Theta \bar y){\thinspace .}$$

The relevant FEniCS code becomes

Fortunately, the solver is independent of all the details of the mesh making. We could also have used the mesh tool mshr in FEniCS, but with our approach here we have full control of the refinement towards the hole.

Solving a problem¶

We assume that ${\alpha}$ is constant. Before solving such a specific problem, it can be wise to scale the problem since it often reduces the amount of input data in the model. Here, the variation in $u$ is typically $|u_a-u_b|$ so we use that as characteristic pressure. The coordinates may be naturally scaled by the bore hole radius, so we have new, scaled variables

$$\bar u = \frac{u-u_a}{u_a-u_b},\quad \bar x = \frac{x}{a},\quad \bar y = \frac{y}{a}{\thinspace .}$$

Now, we expect $\bar u\in [0,1]$, which is a goal of scaling. Inserting this in the problem gives the PDE

$$\nabla^2 \bar u = 0$$

in a domain with inner radius 1 and $\bar u=0$, and outer radius

$$\beta = \frac{a}{b},$$

with $\bar u = 1$. Our solver can solve this problem by setting alpha=1, u_a=1, and u_b=0. We see that the dimensionless parameter $\beta$ goes to the mesh and not to the solver. Figure shows a solution for $\beta =2$ on a mesh with $4\cdot 20\cdot 20 = 1600$ triangles, 25 degree opening, and P1 elements. Switching to higher-order, say P3, is a matter of changing the degree parameter that goes to the function V in the solver:

The complete code is found in borehole_fenics.py. All fluids flow in the same way as long as the geometry is the same!

Solution for (scaled) fluid pressure around a bore hole in a porous medium.

How can we solve a 3D version of this problem? Then we would make a long cylinder. The assumption is that nothing changes in the third direction, so $\partial/\partial z=0$. This means that the cross sections at the end of the cylinder have homogeneous Neumann conditions $\partial u/\partial n=0$. Therefore, nothing changes in the variational form. Actually, all we have to do is to a generate a 3D box and use the same stretching and mapping to make the cylinder, and run the solver without changes!

Summary¶

• When approximating $f$ by $u = \sum_j c_j{\varphi}_j$, the least squares method and the Galerkin/projection method give the same result. The interpolation/collocation method is simpler and yields different (mostly inferior) results.

• Fourier series expansion can be viewed as a least squares or Galerkin approximation procedure with sine and cosine functions.

• Basis functions should optimally be orthogonal or almost orthogonal, because this makes the coefficient matrix become diagonal or sparse and results in little round-off errors when solving the linear system. One way to obtain almost orthogonal basis functions is to localize the basis by making the basis functions that have local support in only a few cells of a mesh. This is utilized in the finite element method.

• Finite element basis functions are piecewise polynomials, normally with discontinuous derivatives at the cell boundaries. The basis functions overlap very little, leading to stable and efficient numerics involving sparse matrices.

• To use the finite element method for differential equations, we use the Galerkin method or the method of weighted residuals to arrive at a variational form. Technically, the differential equation is multiplied by a test function and integrated over the domain. Second-order derivatives are integrated by parts to allow for typical finite element basis functions that have discontinuous derivatives.

• The least squares method is not much used for finite element solution of differential equations of second order, because it then involves second-order derivatives which cause trouble for basis functions with discontinuous derivatives.

• We have worked with two common finite element terminologies and associated data structures (both are much used, especially the first one, while the other is more general):

a. elements, nodes, and mapping between local and global node numbers

b. an extended element concept consisting of cell, vertices, degrees of freedom, local basis functions, geometry mapping, and mapping between local and global degrees of freedom

• The meaning of the word "element" is multi-fold: the geometry of a finite element (also known as a cell), the geometry and its basis functions, or all information listed under point 2 above.

• One normally computes integrals in the finite element method element by element (cell by cell), either in a local reference coordinate system or directly in the physical domain.

• The advantage of working in the reference coordinate system is that the mathematical expressions for the basis functions depend on the element type only, not the geometry of that element in the physical domain. The disadvantage is that a mapping must be used, and derivatives must be transformed from reference to physical coordinates.

• Element contributions to the global linear system are collected in an element matrix and vector, which must be assembled into the global system using the degree of freedom mapping (dof_map) or the node numbering mapping (elements), depending on which terminology that is used.

• Dirichlet conditions, involving prescribed values of $u$ at the boundary, are mathematically taken care of via a boundary function that takes on the right Dirichlet values, while the basis functions vanish at such boundaries. The finite element method features a general expression for the boundary function. In software implementations, it is easier to drop the boundary function and the requirement that the basis functions must vanish on Dirichlet boundaries and instead manipulate the global matrix system (or the element matrix and vector) such that the Dirichlet values are imposed on the unknown parameters.

• Neumann conditions, involving prescribed values of the derivative (or flux) of $u$, are incorporated in boundary terms arising from integrating terms with second-order derivatives by part. Forgetting to account for the boundary terms implies the condition $\partial u/\partial n=0$ at parts of the boundary where no Dirichlet condition is set.