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. Implementation can be an $i$FEM project.

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.

%% 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;