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
(1.1) 
where is a given function. The continuous variational problem is to find such that
where
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) 
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
where the basic data structures node and elem are generated by
For clarity, we take a simple mesh shown in Fig. 1 as an example.
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 https://www.math.uci.edu/~chenlong/ifemdoc/mesh/auxstructuredoc.html 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.
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 twodimensional 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:

A function handle or a constant.

The numerical degrees of freedom of a finite element function.

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
(1.3) 
Here, are the quadrature points on the element .
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 aposteriori error estimate
where and
(2.1) 
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
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
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 onedimensional 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.

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 righthand restrictions are accordingly numbered as

Elementwise sign matrix. To characterize the direction of an edge on an element , one can introduce the elementwise sign matrix
1% sign of elementwise edges2sgnelem = 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:
Here, the data structures elem2edge and bdEdgeIdx have been introduced in Subsect. 1.2.

Elementwise interior indices. According to the indexing rule, one easily obtains the interior indices of the first sides of all triangles:
1% first side2e1 = 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 + (e11)*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”.

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 .
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 1D quadrature points . The Gaussian quadrature points and weights and on are given by quadpts1.m:
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:

1th side: ;

2nd side: ;

3rd side: .
Therefore the points can be given as
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
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
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:
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:
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 points2[elemuhxM,elemuhxP,elemnx,elemny] = elem2edgeInterp(’.dx’,Th,uh,Vh,quadOrder);3[elemuhyM,elemuhyP]␣=␣elem2edgeInterp(’.dy’,Th,uh,Vh,quadOrder); 
The elementwise jumps are

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 sides5 hei = he(elem2edge(:,i));6 id = (1:ng)+(i1)*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(:));14end 
We finally get the local error indicators
Here we add abs since the thirdorder 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.
The initial mesh and the final adapted meshes after 20 and 30 refinement steps are presented in Fig. 4 (ac), 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 (https://github.com/Terenceyuyue/varFEM). The subroutine PoissonVEM_indicator.m is used to compute the local indicator and the test script is main_Poisson_afem.m.
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.
References
 [1] L. Chen. iFEM: an integrated finite element method package in MATLAB. Technical report, University of California at Irvine, 2009.
 [2] F. Hecht. New development in freefem++. J. Numer. Math., 20(34):251–265, 2012.
 [3] Y. Yu. varFEM: variational formulation based programming for finite element methods in Matlab. arXiv:2206.06918, 2022.