# iFEM project: implementing the hierarchical basis for linear finite elements

This note walks through the construction of the stiffness matrix for the finite element approximation of Poisson equation in 2D using the hierarchical basis functions for linear elements.

## Setting

For simplicity the model problem is the pure diffusion problem with homogeneous Dirichlet boundary condition in a polygonal/polyhedral bounded domain $\Omega\subset \mathbb{R}^d$ with $d=2,3$. The construction can be extended for the diffusion reaction problem and/or with other types of boundary conditions very naturally.

\left\{ \begin{aligned} -\nabla\!\cdot\!(A \nabla u) &= f, \; &\text{ in }\, \Omega, \\ u &=0, \; &\text{ on }\,\partial \Omega. \end{aligned} \right. \tag{1}\label{eq:model}

The weak formulation is for problem \eqref{eq:model} is

\left\{ \begin{aligned} & \text{Find } \, u\in H^1_0(\Omega)\, \text{ such that } \\ & \mathcal{A}(u,v) = f(v), \quad \forall v\in H^1_0(\Omega), \end{aligned} \right. \tag{2}\label{eq:pb-w}

where

$\mathcal{A}(u,v) := \int_{\Omega} \nabla u\cdot \nabla v, \quad \text{ and }\; f(v) := \int_{\Omega} fv .$

Suppose we have a triangulation $\mathcal{T}(h)$ of $\Omega$, with mesh size

$h = \max_{K\in \mathcal{T}(h)} \operatorname{diam} K.$

The finite element approximation problem is then replacing the test and trial function spaces $H^1_0(\Omega)$ by the linear Lagrange element space

$V_h := \{v|_K \in P^1(K): v\in C^0(\Omega), \text{and } v|_{\partial \Omega} = 0 \}.$

The stiffness matrix

$A \in M^{N \times N}(\mathbb{R})$ of linear Lagrange element space using the conventional nodal basis function associated with each vertex is then

$A_{ij} = \mathcal{A}(\lambda_j,\lambda_i) = \int_{\Omega} \nabla \lambda_j \cdot \nabla \lambda_i,$

where the $i$-th row corresponds to the $i$-th nodal basis test function, the $j$-th column corresponds to the $j$-th nodal basis for the trial function space, and $N = \dim (\mathcal{N}_h) = \dim V_h$ which is the number of vertices in $\mathcal{T}(h)$.

## Hierarchical basis for Lagrange elements

Now suppose we have two triangulations $\mathcal{T}_H$ and $\mathcal{T}(h)$ with the finer one being refined from the coarser one. The refining procedure can be either longest edge bisection (number of elements doubles), or uniformly refining by connecting the midpoints of each edge (number of elements quadruples).

Coarse and fine triangulation’s corresponding Lagrange elements spaces are $V_H$ and $V_h$ respectively, with $V_H \subset V_h$. For convenience the nodal basis function for the coarse and fine triangulations are denoted differently with a subscript, i.e., $\lambda_{H,i}$ is the $i$-th nodal basis for $V_H$. Let the set of nodal basis functions be $\mathcal{V}_H$ and $\mathcal{V}_h$ respectively. Then denote the index set of the nodal basis function for $V_H$ and $V_h$ as

\begin{aligned} &\Lambda_H = \{ i \in \mathbb{Z}^+: \lambda_{H,i} \in \mathcal{V}_H \text{ is a nodal basis of } V_H \}, \\ \text{ and } & \; \Lambda_h = \{ i \in \mathbb{Z}^+: \lambda_{h,i} \in \mathcal{V}_h \text{ is a nodal basis of } V_h \}, \end{aligned}

so that

$V_H = \operatorname{span}_{i\in \Lambda_H}{\lambda_{H,i}} = \operatorname{span} \mathcal{V}_H \;\text{ and }\; V_h = \operatorname{span}_{i\in \Lambda_h}{\lambda_{h,i}} = \operatorname{span} \mathcal{V}_h.$

Now suppose the indexing of the vertices in $\mathcal{T}_H$ and $\mathcal{T}(h)$ is in a way such that the newly created vertices in $\mathcal{T}(h)$ is labeled as the $(N_H+1)$-th vertex to the $N_h$-th vertex. The enumeration of vertices of $\mathcal{T}_H$ with indices in $\Lambda_H$ is

$\boldsymbol{z}_1, \boldsymbol{z}_2, \dots, \boldsymbol{z}_{N_H}.$

Then the enumeration of its counterpart $\Lambda_h$ of $\mathcal{T}(h)$ is

$\boldsymbol{z}_1, \boldsymbol{z}_2, \dots, \boldsymbol{z}_{N_H}, \boldsymbol{z}_{N_H+1}, \dots, \boldsymbol{z}_{N_h}.\tag{3}\label{eq:ind-h}$

The first $N_H$ vertices in $\mathcal{T}(h)$ inherits from the coarser triangulation.

The hierarchical basis set for linear Lagrange elements on the fine triangulation $\mathcal{T}(h)$ is then:

$\mathcal{S}_h := \{ \lambda_{H,i}\}_{i\in \Lambda_H}\cup \{ \lambda_{h,i}\}_{i\in \Lambda_h\backslash \Lambda_H}. \tag{4}\label{eq:b-tg}$

And the hierarchical Lagrange elements space on $\mathcal{T}(h)$ is then:

$S_h := \Bigl(\operatorname{span}_{i\in \Lambda_H}{\lambda_{H,i}} \Bigr) \cup \Bigl(\operatorname{span}_{i\in \Lambda_h\backslash \Lambda_H}{\lambda_{h,i}}\Bigr) = \operatorname{span} \mathcal{S}_H$
If $\mathcal{T}(h)$ is refined from $\mathcal{T}_H$ using longest edge bisection, then Lagrange element spaces spanned by the nodal basis functions $\mathcal{V}_h$ and the hierarchical basis functions $\mathcal{S}_h$ coincides, i.e.: $$S_h = V_h.$$
It is straightforward to verify that, using the enumeration \eqref{eq:ind-h}, $\lambda_{H,i}$ associated vertex $\boldsymbol{z}_i$ can be represented by the $\{\lambda_{h,j}\}$ associated with $\{\boldsymbol{z}_j\}$'s including $\boldsymbol{z}_i$ and the newly created vertices around $\boldsymbol{z}_i$: $$\lambda_{H,i} = \lambda_{h,i} + \frac{1}{2}\sum_{j\in \mathcal{N}_{h,\boldsymbol{z}_i} } \lambda_{h,j},$$ where $\mathcal{N}_{h,\boldsymbol{z}_i}$ is the set of the vertices neighboring $\boldsymbol{z}_i$ on the fine triangulation $$\mathcal{N}_{h,\boldsymbol{z}_i} := \{\boldsymbol{z}_j \in \mathcal{N}_h: \boldsymbol{z}_i \boldsymbol{z}_j \text{ is an edge of }\mathcal{T}(h) \}.$$

## Multigrid projections to construct matrices

Denote the stiffness matrices using nodal basis for $V_H$ and $V_h$ as $A^H$ and $A_h$ respectively. Denote the stiffness matrices using hierarchical basis for $S_h$ as $\widehat{A}^h$.

To construct $\widehat{A}^h$, we need a few more notations and operators usually used in multigrid language. The prolongation operator or natural inclusion $\mathcal{P}_{H \to h}$, then can be defined as

$\mathcal{P}_{H\to h}: V_H \to V_h,\; w_H \mapsto w_H.$

Any coarse piecewise linear continuous function can be written as follows:

$w_H = \sum_{i=1}^{N_H} w_H(\boldsymbol{z}_i)\lambda_i = \mathbf{w}_H\cdot \boldsymbol{\lambda}_H,\; \text{ where } \;\boldsymbol{\lambda}_H = (\lambda_{H,1},\dots,\lambda_{H,N_H})^T.$

Now the prolongation operator has the following matrix form:

$\mathcal{P}_{H\to h} w_H = \mathcal{P}_{H\to h} \Big(\mathbf{w}_H\cdot \boldsymbol{\lambda}_H\Big) = (P_H \mathbf{w}_H)^T \boldsymbol{\lambda}_h,$

where

$P_H \in M_{N_h\times N_H}(\mathbb{R}) :\mathbb{R}^{N_H} \to \mathbb{R}^{N_h}$

is the matrix representation of the prolongation operator.

The stiffness matrix for the basis set \eqref{eq:b-tg} can be then computed by evaluating the following three parts. The first part is for

$\{ \lambda_{H,i}\}_{i\in \Lambda_H}$

by directly passing the stiffness matrix $A_H$ from the coarser triangulation. The second part is for

$\{ \lambda_{h,i}\}_{i\in \Lambda_h\backslash \Lambda_H}$

by taking the block from $A_h$ accounting for the nodal basis function associated with the newly created vertices. The third part are for the cross terms.

For the first part, the integral on the coarser triangulation can be naturally computed on the finer triangulation since $V_H\subset V_h$: for $1\leq i,j \leq N_H$

$\mathcal{A}(\lambda_{H,j}, \lambda_{H,i}) = \mathcal{A}\Bigl(\mathcal{P}_{H\to h}\lambda_{H,j}, \mathcal{P}_{H\to h}\lambda_{H,i}\Bigr),$

and this implies by matrix representation

$A_H = (P_H)^T A_h P_H.$

Similarly, the cross term then is: for $1\leq i \leq N_H$, and $N_H+1\leq j\leq N_h$,

$\mathcal{A}(\lambda_{h,j}, \lambda_{H,i}) = \mathcal{A}\Bigl(\lambda_{h,j}, \mathcal{P}_{H\to h}\lambda_{H,i}\Bigr).$

Lastly, the lower right block consists the Lagrange basis functions for the newly created vertices on the fine mesh $N_H+1\leq i,j \leq N_h$, and its matrix representation is

$Q^T A_h Q := \begin{pmatrix}0 & I \end{pmatrix} A_h \begin{pmatrix} 0 \\ I \end{pmatrix},$

where $I = I_{N_h-N_H}$ is the identity matrix.

To sum up, the $\widehat{A}^h$ consists four blocks:

$\label{eq:mat-h} \widehat{A}^h = \begin{pmatrix} (P_H)^TA_h P_H & (Q^T A_h P_H)^T \\ Q^T A_h P_H & Q^T A_h Q \end{pmatrix}.$

Now if we have $n$ levels, and the stiffness matrix $\widehat{A}^{h_n}$ for the hierarchical basis is know, then the stiffness matrix $\widehat{A}^{h_{n+1}}$ is just

$\label{eq:mat-hn} \widehat{A}^{h_{n+1}} = \begin{pmatrix} \widehat{A}^{h_n} & * \\ * & (Q^{n})^T {A}^{h_{n+1}} Q^{n} \end{pmatrix},$

where

$(Q_{n})^T := \begin{pmatrix}0 & I_{N_{h_{n+1}} - N_{h_{n}}} \end{pmatrix},$

the off-diagonal blocks are a bit complicated to be represented in a compact matrix form. Here $(Q^{n})^T {A}^{h_{n+1}} Q^{n}$ represent the lower right block from $A^{h_{n+1}}$ that accounts for the basis of the newly created vertices.

For 2 levels of triangulation ($\mathcal{T}({h_i})$ with $i=1,2$), and we want to obtain the stiffness matrix for level 3 triangulation $\mathcal{T}({h_3})$, we have:

$\widehat{A}^{h_{3}} = \begin{pmatrix} A^{h_1} & * & (Q_3^T A^{h_{3}} P_{h_2} P_{h_1})^T \\ * & Q_{2}^T {A}^{h_{2}} Q_{2} & (Q_3^T A^{h_{3}} P_{h_2} Q_2)^T \\ Q_3^T A^{h_{3}} P_{h_2} P_{h_1} & Q_3^T A^{h_{3}} P_{h_2} Q_2 & Q_{3}^T A^{h_{3}} Q_{3} \end{pmatrix}.\tag{5} \label{eq:mat-h3}$

For example, the $Q_3^T A^{h_{3}} P_{h_2} Q_2$ is taking the $(N_{h_2}+1)$-th row to the $N_{h_3}$-th row, and $(N_{h_1}+1)$-th column to the $N_{h_2}$-th column from $A^{h_{3}} P_{h_2}$. This blocks accounts for the $\mathcal{A}(\lambda_{h_2,j}, \lambda_{h_3,i})$ for $N_{h_1}+1\leq j\leq N_{h_2}$, which corresponds to the trial functions on $\mathcal{T}({h_2})$, and $N_{h_2}+1\leq i\leq N_{h_3}$, which corresponds to the test functions on
$\mathcal{T}({h_3})$.

To sum up, the construction is exploiting the existing multigrid data structure.

## MATLAB implementation using iFEM

The following implementation assumes that we have iFEM by Long Chen in our MATLAB’s $PATH. The functions we have used are: • bisect: bisecting the uniform mesh to get a criss-cross type mesh. • HBstructure: the hierarchical structure, a newly added node is neighboring to which parents. • transferoperator: interpolations from coarse to fine mesh; and restrictions from fine to coarse mesh. • getH1error: computing$H^1\$-error of the approximation to the solution.
• Poisson: computing a reference solution in non-hierarchical basis, and assembling matrices in a single level.

With the help of these functions, this is a good exercise to test our skills in writing for loops which heavy interactions between loops.

Tags:

Categories:

Updated: