# 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\]## 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 *i*FEM

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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221

%% Solving the diffusion equation test problem on a square [0,1]^2
% DIFFUSION2DHR solves the Poisson equation using direct method
% with linear hierarchical basis in multi-level
% this an un-incapsulated file to view all data structures
%
% Equation: $-div(a grad u)=f$
% Domain: $\Omega = (0,1)^2$
% Dirichlet data: $u = 0$ on $\partial \Omega$
% Finite element spaces: hierarchical linear Lagrange element
% Requirements of packages in MATLAB
% 1. Computational Geometry, TriRep class
% 2. Linear Algebra, mldivide
% 3. Visualization, rendering patch object
% 4. Version >= 2013a
% Requirements for external pkgs
% Long Chen's iFEM
%
% See also delaunayTriangulation, mldivide.
% S Cao 2015
%% set up
clear; close all force; clc;
help DIFFUSION2DHR
maxLvl = input(['Enter number of max levels\n',...
'(>=2, default = 12, recommend >=5) \n',...
'hit enter for the default value:']);
if isempty(maxLvl) || (~isempty(maxLvl) && ~isnumeric(maxLvl)) || maxLvl<=1
maxLvl = 12; % default number of levels
elseif ~isinteger(maxLvl)
maxLvl = ceil(maxLvl);
end
maxNnode = 1e5;
%% coarsest possible criss-cross type mesh
% initial mesh
% 4-----3
% | \ / |
% | 5 |
% | / \ |
% 1-----2
nodeH = [-1,-1; 1,-1; 1,1; -1,1; 0,0];
elemH = [5,2,3; 5,4,1; 5,1,2; 5,3,4];
NH = size(nodeH,1);
% #nodes at coarse
elemH = label(nodeH,elemH);
%% set up multilevel datastructure
node = cell(maxLvl,1);
Nnode = zeros(maxLvl,1); % number of vertices
elem = cell(maxLvl,1);
Nelem = zeros(maxLvl,1); % number of elements
dof = cell(maxLvl,1); % degrees of freedom excluding the boundary nodes
nodeNew = cell(maxLvl,1); % newly created dof at each lvl
Nnew = zeros(maxLvl,1); % Nnew(lvl) = length(nodeNew{lvl})
NL = cell(maxLvl,1);
HB = cell(maxLvl,1);
% AHr is the hierarchical stiffness matrix after each refinement
AHr = cell(maxLvl,1);
% Ah is the stiffness matrix on current level after each refinement
Ah = cell(maxLvl,1);
% PAHh is the transfer matrix of all previous hierarchical spaces to current lvl
PAHh = cell(maxLvl-1,1);
bHr = cell(maxLvl,1);
% bHr is the hierarchical test bases' RHS
uHr = cell(maxLvl,1);
% uHr is the hierarchical element approximation
uHr2h = cell(maxLvl,1);
% uHr2h is the hierarchical element approximation interpolated to Lagrange
uh = cell(maxLvl,1);
% uh is the ordinary lagrange linear element approximation
bh = cell(maxLvl,1);
% bh is the lagrange test bases' RHS at each level
Pro = cell(maxLvl,1); %Pro{lvl-1}: V_{lvl-1} to V_{lvl}
Res = cell(maxLvl,1); %Res{lvl}: V_{lvl} to V_{lvl-1}
ProH = cell(maxLvl,1); % ProH{lvl}: coarsest to V_{lvl}
ResH = cell(maxLvl,1); % ResH{lvl}: V_{lvl} to coarsest
Proh = cell(maxLvl,1); % Proh{k}: V_{k} to current finest lvl
Resh = cell(maxLvl,1); % Resh{k}: current finest lvl to V_{k}
etaL = cell(maxLvl,1); % etaL{k}: error indicator for Lagrange element at lvl=k
etaH = cell(maxLvl,1); % etaH{k}: error indicator for Hierarchical basis at lvl=k
etaLTotal = zeros(maxLvl,1);
etaHTotal = zeros(maxLvl,1);
errL = zeros(maxLvl,1);
errH = zeros(maxLvl,1);
%% Right hand side data, diffusion coefficient, boundary condition
uH1seminorm = sqrt(2*pi^2); % \|\nabla u\|_{L^2}
u = @(p) sin(pi*p(:,1)).*sin(pi*p(:,2)); % u is true solution
f = @(p) 2*pi^2*u(p); % f is the right hand side data
gradu = @(p) pi*[cos(pi*p(:,1)).*sin(pi*p(:,2)) ...
sin(pi*p(:,1)).*cos(pi*p(:,2))];
a = @(p) ones(size(p,1),1); % diffusion constant is 1 here
pde = struct('f',f,'u',u,'g_D',0,'Du',gradu,'a',a);
%% generating the coarse matrices use the subroutine from iFEM
lvl=1;
node{lvl} = nodeH; elem{1} = elemH;
Nnode(lvl) = size(node{lvl},1);
Nelem(lvl) = size(elem{lvl},1);
uh{lvl} = zeros(Nnode(lvl),1);
[uh{lvl},~,eqn] = Poisson(node{lvl},elem{lvl},pde,[]);
uHr{lvl} = uh{lvl};
uHr2h{lvl} = uh{lvl};
AHr{lvl} = eqn.Lap; % stiffness (un-modified version)
bHr{lvl} = eqn.b;
Ah{lvl} = AHr{lvl}; % coarsest level hierarchical = lagrange
bh{lvl} = bHr{lvl};
dof{lvl} = eqn.freeNode;
nodeNew{lvl} = 1:Nnode(lvl);
NL{lvl}= [0 Nnode(lvl)];
fprintf('\nLevel is %d with #DoF=%2d. \n',lvl,Nnode(lvl));
%% compute the error at lvl=1
errL(lvl) = getH1error(node{lvl},elem{lvl},pde.Du,uh{lvl});
errH(lvl) = errL(lvl); % same error for hierarchical or Lagrange
%% the iterative refinement loop
for lvl = 2:maxLvl
%% refine from lvl-1 (known) to lvl
[node{lvl},elem{lvl},~,HB{lvl}] = bisect(node{lvl-1},elem{lvl-1});
Nnode(lvl) = size(node{lvl},1);
Nelem(lvl) = size(elem{lvl},1);
fprintf('\nLevel is %d with #DoF=%2d. \n',lvl,Nnode(lvl));
fprintf('Rel-err in | |_{H^1} is %6g in previous level. \n', errL(lvl-1)/uH1seminorm);
%% find the lagrange basis matrix
[uh{lvl},~,eqn] = Poisson(node{lvl},elem{lvl},pde);
Ah{lvl} = eqn.Lap; % stiffness
bh{lvl} = eqn.b;
dof{lvl} = eqn.freeNode;
%% find the hierarchical structure between lvl-1 and lvl
NL{lvl} = [0 Nnode(lvl-1) Nnode(lvl)];
N0 = size(elem{lvl-1});
[HB_tg, ~, ~] = HBstructure(elem{lvl},N0);
[Pro_temp,Res_temp] = transferoperator(HB_tg,NL{lvl});
Pro{lvl-1} = Pro_temp{1};
Res{lvl} = Res_temp{2};
if lvl-1 == 1 % when there are only two levels, the operators needs special treatments
ProH{lvl-1} = Pro{lvl-1};
ResH{lvl} = Res{lvl};
else
ProH{lvl-1} = Pro{lvl-1}*ProH{lvl-2};
ResH{lvl} = ResH{lvl-1}*Res{lvl};
end
Proh{lvl-1} = Pro{lvl-1};
Resh{lvl} = Res{lvl};
for j = 1:lvl-2
Proh{j} = Pro{lvl-1}*Proh{j};
Resh{j+1} = Resh{j+1}*Res{lvl};
end
%% generating the hierarchical matrices for two levels
% Nnode(lvl-1)+1:Nnode(lvl) are newly created nodes (DoF) at lvl
nodeNew{lvl} = Nnode(lvl-1)+1:Nnode(lvl);
Nnew(lvl) = length(nodeNew{lvl});
PAHh{lvl-1} = zeros(Nnew(lvl),Nnode(lvl-1));
% PAh: inner prod matrix of \int \nabla Prolongation of phi_{H,i} \cdot \nabla phi_{h,j}
if lvl == 2
PAh_tmp = Ah{lvl}*Pro{lvl-1};
PAHh{lvl-1} = PAh_tmp(nodeNew{lvl},:);
else
% jj = 1
PAh_tmp = Ah{lvl}*Proh{1};
PAHh{1} = PAh_tmp(nodeNew{lvl},:);
for jj = 2:lvl-1
PAh_tmp = Ah{lvl}*Proh{jj};
PAHh{jj} = [PAHh{jj-1} PAh_tmp(nodeNew{lvl},nodeNew{jj})];
end
end
%% compute the stiffness matrix entry between levels
AHr{lvl} = [AHr{lvl-1} PAHh{lvl-1}'; ...
PAHh{lvl-1}, Ah{lvl}(nodeNew{lvl},nodeNew{lvl})];
%% compute the right handside bHr
% each time after refine, the RHS is computed using finer representation
if lvl == 2
bHr{lvl} = [bHr{lvl-1}; bh{lvl}(nodeNew{lvl})];
else
bHr{1} = Proh{1}'*bh{lvl};
for jj = 2:lvl-1
bHh_tmp = Proh{jj}'*bh{lvl};
bHr{jj} = [bHr{jj-1}; bHh_tmp(nodeNew{jj})];
end
bHr{lvl} = [bHr{lvl-1}; bh{lvl}(nodeNew{lvl})];
end
%% solve the hierarchical element linear ssystem using A\b
uHr{lvl} = zeros(Nnode(lvl),1);
tic;
uHr{lvl}(dof{lvl}) = AHr{lvl}(dof{lvl},dof{lvl})\bHr{lvl}(dof{lvl});
htime=toc;
fprintf('MATLAB A\b hierarchical linear system time = %6.4g s\n',htime);
%% interpolate the hierarchical solution to the current level to verify
uHr2h{lvl} = zeros(Nnode(lvl),1);
% interpolate every part of the current level hierarchical solution
% from bases of lvl=1 to (lvl-1) onto the current lvl
uHr2h{lvl} = Proh{1}*uHr{lvl}(nodeNew{1});
for k = 2:lvl-1
uHr2h{lvl} = uHr2h{lvl} + Proh{k}*[zeros(Nnode(k-1),1);uHr{lvl}(nodeNew{k})];
end
uHr2h{lvl} = uHr2h{lvl} + [zeros(Nnode(lvl-1),1); uHr{lvl}(nodeNew{lvl})];
%% compute the error and estimator for uh and hierarchical uh at lvl
errL(lvl) = getH1error(node{lvl},elem{lvl},pde.Du,uh{lvl});
errH(lvl) = getH1error(node{lvl},elem{lvl},pde.Du,uHr2h{lvl});
%% stopping criteria for
if (size(node{lvl},1)>maxNnode); break; end
end
%% plot convergence
N = Nnode(1:lvl);
errL = errL(1:lvl)/uH1seminorm;
errH = errH(1:lvl)/uH1seminorm;

## Comments