Log In Sign Up

Unified Implementation of Adaptive Finite Element Methods in Matlab

We provide a unified implementation of the adaptive finite element methods for the Poisson equation in Matlab by combing the use of the software package varFEM. The design idea can be extended to other Galerkin-based methods, for example, the discontinuous Galerkin methods and the virtual element methods.


page 5

page 11


Hybridization and postprocessing in finite element exterior calculus

We hybridize the methods of finite element exterior calculus for the Hod...

Unfitted Trefftz discontinuous Galerkin methods for elliptic boundary value problems

We propose a new geometrically unfitted finite element method based on d...

Stable implementation of adaptive IGABEM in 2D in MATLAB

We report on our MATLAB program package IGABEM2D, which provides an easi...

Implementation of hp-adaptive discontinuous finite element methods in Dune-Fem

In this paper we describe generic algorithms and data structures for the...

Automatic symbolic computation for discontinuous Galerkin finite element methods

The implementation of discontinuous Galerkin finite element methods (DGF...

Discretizations of Stochastic Evolution Equations in Variational Approach Driven by Jump-Diffusion

Stochastic evolution equations with compensated Poisson noise are consid...

1 Finite element methods for the Poisson equation

In this section we briefly review the finite element method and introduce the numerical implementation in varFEM, a Matlab software package for the finite element method. For simplicity, we consider the Poisson equation with the homogeneous Dirichlet boundary conditions.

1.1 The finite element method for the Poisson equation

Let be a bounded Lipschitz domain in with polygonal boundary . Consider the Poisson equation


where is a given function. The continuous variational problem is to find such that


For the finite element discretization, we discuss the conforming Lagrange elements. Let be a shape regular triangulation. The generic element will be denoted as in the sequel. Define

where . The discrete problem is to seek satisfying


1.2 Data structures for triangular meshes

We adopt the data structures given in FEM [1]. All related data are stored in the Matlab structure Th, which is computed by using the subroutine FeMesh2d.m as

1Th = FeMesh2d(node,elem,bdStr);

where the basic data structures node and elem are generated by

1[node,elem] = squaremesh([x1 x2 y1 y2], h1, h2);

For clarity, we take a simple mesh shown in Fig. 1 as an example.

Fig. 1: Illustration of the data structures

The triangular meshes are represented by two basic data structures node and elem, where node is an matrix with the first and second columns contain - and -coordinates of the nodes in the mesh, and elem is an matrix recording the vertex indices of each element in a counterclockwise order, where N and NT are the numbers of the vertices and triangular elements. For the mesh given in Fig. 1,

In the current version, we only consider the -Lagrange finite element spaces with up to 3. In this case, there are two important data structures edge and elem2edge. In the matrix edge(1:NE,1:2), the first and second rows contain indices of the starting and ending points. The column is sorted in the way that for the -th edge, edge(k,1)¡edge(k,2) for . For the given triangulation,

The matrix elem2edge establishes the map of local index of edges in each triangle to its global index in matrix edge. By convention, we label three edges of a triangle such that the -th edge is opposite to the -th vertex. For the given mesh,

We refer the reader to for some detailed information.

To deal with boundary integrals, we first extract the boundary edges from edge and store them in matrix bdEdge. In the input of FeMesh2d, the string bdStr is used to indicate the interested boundary part in bdEdge. For example, for the unit square ,

  • bdStr = ’x==1’ divides bdEdge into two parts: bdEdgeType{1} gives the boundary edges on , and bdEdgeType{2} stores the remaining part.

  • bdStr = {’x==1’,’y==0’} separates the boundary data structure bdEdge into three parts: bdEdgeType{1} and bdEdgeType{2} give the boundary edges on and , respectively, and bdEdgeType{3} stores the remaining part.

  • bdStr = [] implies that bdEdgeType{1} = bdEdge.

We also use bdEdgeIdxType to record the index in matrix edge, and bdNodeIdxType to store the nodes for respective boundary parts. Note that we determine the boundary of interest by the coordinates of the midpoint of the edge, so ’x==1’ can also be replaced by a statement like ’x¿0.99’.

1.3 Implementation of the FEM in varFEM

FreeFEM is a popular 2D and 3D partial differential equations (PDE) solver based on finite element methods

[2], which has been used by thousands of researchers across the world. The highlight is that the programming language is consistent with the variational formulation of the underlying PDEs, referred to as the variational formulation based programming in [3], where we have developed an FEM package in a similar way of FreeFEM using the language of Matlab, named varFEM. The similarity here only refers to the programming style of the main or test script, not to the internal architecture of the software.

Consider the Poisson equation with Dirichlet boundary condition on the unit square. The exact solution is given by

The PDE data is generated by pde = Poissondata_afem(). The function file is simply given as follows.

1function uh = varPoisson(Th,pde,Vh,quadOrder)
3%% Assemble stiffness matrix
4Coef  = 1;
5Test  = ’v.grad’;

In the above code, the structure pde stores the information of the PDE, including the exact solution pde.uexact, the gradient pde.Du, etc. We set up the triple (Coef,Test,Trial) for the coefficients, test functions and trial functions in variational form, respectively. It is obvious that v.grad is for and v.val is for itself. The routine assem2d.m computes the stiffness matrix corresponding to the bilinear form on the two-dimensional region, i.e.

where are the global shape functions of the finite element space Vh. The integral of the bilinear form, as , is approximated by using the Gaussian quadrature formula with quadOrder being the order of accuracy.

We remark that the Coef has three forms:

  1. A function handle or a constant.

  2. The numerical degrees of freedom of a finite element function.

  3. A coefficient matrix CoefMat resulting from the numerical integration.

In the computation, the first two forms in fact will be transformed to the third one. Given a function , where , the coefficient matrix is in the following form


Here, are the quadrature points on the element .

We display the numerical result in Fig. 2 for the uniform triangular mesh with generated by

1[node,elem] = squaremesh([0 1 0 1], 1/50, 1/50);

Fig. 2: Numerical and exact solutions with

2 Adaptive finite element methods for the Poisson equation

In this section we briefly introduce the ingredients of the adaptive finite element method, and provide the overall structure of the implementation.

2.1 The ingredients of the adaptive FEM

For the conforming FEM, one can establish the following residual based a-posteriori error estimate

where and


are referred to as the global and local error indicators, respectively. Here, is the interior residual, and

is the jump term, where

is the unit outer normal vector of the target element

. Note that for an edge on the domain boundary , we assume .

Standard adaptive algorithms based on the local mesh refinement can be written as loops of the form

Given an initial subdivision , to get from we first solve the FEM problem under consideration to get the numerical solution on . The error is then estimated by using , and the a posteriori error bound . The local error bound is used to mark a subset of elements in for refinement. The marked triangles and possible more neighboring elements are refined in such a way that the subdivision meets certain conditions, for example, the resulting triangular mesh is still shape regular. The above procedures are included in the test script with an overview given as follows

1for k = 1:maxIt
2    % Step 1: SOLVE
3    Th = FeMesh2d(node,elem,bdStr);
4    uh = varPoisson(Th,pde,Vh,quadOrder);
5    % Step 2: ESTIMATE
6    eta = Poisson_indicator(Th,uh,pde,Vh,quadOrder);
7    % Step 3: MARK
8    elemMarked = mark(elem,eta,theta);
9    % Step 4: REFINE
10    [node,elem] = bisect(node,elem,elemMarked);

Three important modules are involved: the local error indicator in Step 2, the marking algorithm in Step 3 and the local refinement algorithm in Step 4. We employ the Dörfler marking strategy to select the subset of elements and then use the newest vertex method to refine the mesh. Note that the subroutines mark.m and bisect.m are extracted from FEM [1] with some minor modifications. In this article, we are only concerned with the computation of the error indicator. For the later two steps, we refer the reader to [1] for details on the implementation.

2.2 The unified implementation of the error indicator

2.2.1 The computation of the elementwise residuals

With the help of varFEM, the first term in (2.1) can be simply computed as

1%% elementwise residuals
2fc = interp2dMat(pde.f,’.val’,Th,Vh,quadOrder);

In the above code, elemIh stores all the local error indicators ; The function interp2dMat is used to generate the coefficient matrix (see (1.3)). It is obvious that the coefficient matrix of is (fc+uxxc+uyyc).^2, where fc, uxxc and uyyc are the coefficient matrices of , and , respectively.

In what follows, we focus on the unified implementation of the jump term or jump integral . We remark that any other type of jump terms can be easily adapted or designed accordingly as will be seen.

2.2.2 The elementwise interior and exterior indices of the quadrature points

The integral over is calculated by using the one-dimensional Gaussian numerical integration formula

where and are the quadrature weights and points on , and is the number of the quadrature points. Note that the endpoints of are not included.

Fig. 3: Illustration of the quadrature points
  1. Indexing rule for the quadrature points. In Fig. 3 the direction of the edge is specified by the arrow. In the implementation, the direction is determined by the data structure edge which satisfies edge(k,1)¡edge(k,2). For the first edge in edge, the quadrature points will be numbered by or for short when restricted to the left triangle. Similarly, the quadrature points on the -th edge will be numbered by

    The right-hand restrictions are accordingly numbered as

  2. Elementwise sign matrix. To characterize the direction of an edge on an element , one can introduce the elementwise sign matrix

    1% sign of elementwise edges
    2sgnelem = sign([elem(:,3)-elem(:,2), elem(:,1)-elem(:,3), elem(:,2)-elem(:,1)]);

    In some cases, it is better to restore the positive sign for edges on the boundary of the domain, or one can set the signs to be zero for later use:

    1E = false(NE,1); E(bdEdgeIdx) = 1; sgnbd = E(elem2edge);
    2sgnelem(sgnbd) = 0;

    Here, the data structures elem2edge and bdEdgeIdx have been introduced in Subsect. 1.2.

  3. Elementwise interior indices. According to the indexing rule, one easily obtains the interior indices of the first sides of all triangles:

    1% first side
    2e1 = elem2edge(:,1); sgn1 = sgnelem(:,1);
    3id = repmat(1:ng, NT,1);
    4id(sgn1<0,:) = repmat((ng:-1:1)+NE*ng, sum(sgn1<0), 1);
    5elemQuade1 = id + (e1-1)*ng;

    Note that for the edges with positive and zero signs the natural indices are , while for those with negative signs the natural indices are reversed. One can similarly introduce the interior indices of other sides, and hence give the elementwise interior indices

    1elemQuadM = [elemQuade1,elemQuade2,elemQuade3];

    where the letter M stands for “Minus”.

  4. Elementwise exterior indices. To compute the jump, we also need to introduce the elementwise exterior indices elemQuadP, where P is for “Plus”. Given an edge, assume that the interior indices of the quadrature points are , and the exterior indices are . Since , one just needs to subtract for those indices in elemQuadM greater than , and add to those indices less than .

    1index = ( elemQuadM > ng*NE ) ;
    2elemQuadP = elemQuadM + (-ng*NE)*index + ng*NE*(~index);

    Obviously, for edges on the domain boundary, the exterior indices are greater than .

2.2.3 The elementwise interior and exterior evaluations

Given , we now determine the interior evaluations elemuhM and the exterior evaluations elemuhP of corresponding to elemQuadM and elemQuadP. To this end, one can compute the evaluations of the basis functions at the quadrature points along the boundary . Given a triangle , let , and be the barycentric coordinate functions. Since are usually used to construct the basis functions for the FEMs, one can first specify the evaluations of and the derivatives and at the quadrature points.

Given some points on the triangle , for example, the quadrature points for the integration of the bilinear forms. In varFEM, the evaluations of are stored in the following form

At this time, one can choose as the quadrature points along the boundary . Let . The quadrature points on are denoted by . In this case, , and one can specify the values by using the 1-D quadrature points . The Gaussian quadrature points and weights and on are given by quadpts1.m:

1[lambda1d,weight1d] = quadpts1(4);  ng = length(weight1d);
2[~,id] = sort(lambda1d(:,1));
3lambda1d = lambda1d(id,:); weight1d = weight1d(id);

Here we use sort to guarantee . Note that

and all the row sums are 1, i.e. . By the definition of the barycentric coordinate functions, the coordinates on the three sides of are:

  • 1-th side: ;

  • 2nd side: ;

  • 3rd side: .

Therefore the points can be given as

1function [lambdaBd,weightBd] = quadptsBd(order)
2%%Gauss quadrature points along the boundary of triangles
4% quadrature on [0,1]
5[lambda1d,weight1d] = quadpts1(order);  ng = length(weight1d);
6[~,id] = sort(lambda1d(:,1));
7lambda1d = lambda1d(id,:); weight1d = weight1d(id);
8% quadrature on each side
9lambdae1 = [zeros(ng,1), lambda1d(:,2), lambda1d(:,1)];
10lambdae2 = [lambda1d(:,1), zeros(ng,1), lambda1d(:,2)];
11lambdae3 = 1 - lambdae1 - lambdae2;
12% quadrature along the boundary of each element
13lambdaBd = [lambdae1; lambdae2; lambdae3];
14weightBd = repmat(weight1d,1,3);

With these quadrature points one can construct the evaluations of the basis functions. The details are omitted. Please refer to Base2dBd.m for implementation details. For the element, the basis functions are , and the elementwise interior evaluations of are given as follows

1elemuhM = zeros(NT,3*ng);
2for p = 1:3*ng
3    % interpolation at the p-th quadrture point
4    for i = 1:length(phi)
5        base = phi{i};
6        elemuhM(:,p) = elemuhM(:,p) + uh(elem2dof(:,i)).*base(:,p);
7    end

For other types of finite elements, one just needs to modify the basis functions base. For evaluations of the derivatives, for example, , one needs to change base to the partial derivative of the basis functions. The elementwise exterior evaluations of can be computed as

1uhI = zeros(2*NE*ng,1);
2uhI(elemQuadM) = elemuhM;
3elemuhP = uhI(elemQuadP);

We remark that the exterior evaluations on the domain boundary are zero since the corresponding indices in elemQuadM are less than or equal to .

The above discussions are summarized in a subroutine:

1function [elemuhM,elemuhP] = elem2edgeInterp(wStr,Th,uh,Vh,quadOrder)

One can obtain the evaluations of , and by setting wStr as .val, .dx and .dy, respectively. For convenience, we also return the elementwise unit normal vectors on the quadrature points:

1%% elementwise unit normal vectors of edges on quadrature points
2if nargout==4
3    rep = ones(1,ng);
4    z1 = node(elem(:,1),:);
5    z2 = node(elem(:,2),:);
6    z3 = node(elem(:,3),:);
7    e1 = z2-z3;   e1 = e1./vecnorm(e1,2,2); % -e1
8    e2 = z3-z1;   e2 = e2./vecnorm(e2,2,2);
9    e3 = z1-z2;   e3 = e3./vecnorm(e3,2,2);
10    elemQuadnx = -[e1(:,2*rep),e2(:,2*rep),e3(:,2*rep)];
11    elemQuadny = [e1(:,rep),e2(:,rep),e3(:,rep)];

2.2.4 The computation of the jump integral

With the above preparations, we are able to compute the jump integral

  • The first step is to compute the elementwise interior and exterior evaluations of and :

    1%% elementwise interior and exterior evaluations at quadrature points
    2[elemuhxM,elemuhxP,elemnx,elemny] = elem2edgeInterp(’.dx’,Th,uh,Vh,quadOrder);
  • The elementwise jumps are

    1elem2Jumpx = elemuhxM - elemuhxP;
    2elem2Jumpy = elemuhyM - elemuhyP;
  • The jump integral can be computed by looping of triangle sides:

    1elemJump = zeros(NT,1);
    2[~,weight1d] = quadpts1(quadOrder);
    3ng = length(weight1d);
    4for i = 1:3 % loop of triangle sides
    5    hei = he(elem2edge(:,i));
    6    id = (1:ng)+(i-1)*ng;
    7    cei = hei;
    8    neix = elemnx(:,id);
    9    neiy = elemny(:,id);
    10    Jumpnx = elem2Jumpx(:,id).*neix;
    11    Jumpny = elem2Jumpy(:,id).*neiy;
    12    Jumpn = (Jumpnx+Jumpny).^2;
    13    elemJump = elemJump + cei.*hei.*(Jumpn*weight1d(:));
  • We finally get the local error indicators

    1%% Local error indicator
    2eta = (abs(elemRes) + elemJump).^(1/2);

    Here we add abs since the third-order quadrature rule has negative weights (???).

3 Numerical example

Consider the example given in Subsect. 1.3. We employ the Dörfler marking strategy with parameter to select the subset of elements for refinement.

Fig. 4: The initial and the final adapted meshes. (a) The initial mesh; (b) After 20 refinement steps; (c) After 30 refinement steps

The initial mesh and the final adapted meshes after 20 and 30 refinement steps are presented in Fig. 4 (a-c), respectively, where the element is used. We also plot the adaptive approximation in Fig. 5, which almost coincides with the exact solution. The convergence rates of the error estimators and the errors in norm are shown in Fig. 6, from which we observe the optimal rates of convergence as predicted by the theory. The full code is available from varFEM package ( The subroutine PoissonVEM_indicator.m is used to compute the local indicator and the test script is main_Poisson_afem.m.

Fig. 5: The exact and numerical solutions
Fig. 6: The convergence rates of the error estimators and the errors in norm

4 Concluding remarks

In this paper, a unified implementation of the adaptive finite element methods was presented for the model problem, which combines the use of the Matlab software package varFEM. The design idea can be extended to other types of Galerkin methods.