LEBREF2D: a 2D mesh-refinement toolbox

A short documentation

Leonardo Rocchi
September, 2018

LEBREF2D is a MATLAB package for the mesh-refinement of 2D polygonal domains based

on an efficient usage of MATLAB built-in functions as well as vectorization features. The
package includes functions for uniform mesh-refinements and local mesh-refinements
based on the longest edge bisection (LEB) algorithm which is a variant of the newest vertex
bisection (NVB) algorithm. LEBREF2D can be easily incorporated into, e.g., an adaptive
Finite Element algorithm for solving partial differential equations.

Contents
1 Data structures 2

2 Longest Edge Bisection 4

3 Detail-grid 5

4 Solving PDEs 6

5 List of functions 7

6 Useful commands 10

References 11
1 Data structures
Consider a two-dimensional domain D ⊂ ‘2 , with a regular polygonal boundary denoted by
D . Let T be a conforming mesh (or triangulation) of the domain D , i.e.,
(i) T = {K1 , . . . , KM } is a finite set of M > 0 triangles (elements) Ki := conv{x j , xm , xn } such that
|Ki | > 0, i = 1, . . . , M. Here, x j , xm , and xn are the nodes (or vertices) of Ki and | · | denotes
the area;
(ii) the triangulation T covers the closure D of D ;
(iii) the intersection of two elements is either empty, a common node, or a common edge
(an edge is naturally defined as the one-dimensional segment connecting two nodes, i.e.,
conv{xm , xn });
(iv) a boundary edge (i.e., those edges lying on D ) belongs to one element only.
Figure 1 shows three examples of non-conforming meshes where, in particular, property (iii)
above is not satisfied.
Let N be the number of nodes (including the boundary ones), i.e., N := #N , where N :=
{x1 , . . . , xN } denotes the set of nodes. Also, let M bet the number of elements of the mesh, i.e.,
M := #T . The mesh is represented by a structure MESH which contains five fields, coord, elem,
int, bnd, elbnd:
• MESH.coord is a N ×2 matrix containing the physical coordinates of all nodes of the mesh.
The i-th row of MESH.coord stores the coordinates of the i-th node xi = (xi , yi ) ∈ D :
MESH.coord(i,:) = [xi yi ].

• MESH.elem is a M × 3 matrix containing the elements nodes’ numbers. That is, the „-th
element K„ = {xi , x j , xk } ∈ T with nodes xi , x j , xk ∈ N is stored in
MESH.elem(„,:) = [i j k],
with the nodes counted in counterclockwise order.
• MESH.int is a ND × 1 vector containing the interior nodes’ numbers (here, ND denotes the
number of interior nodes). For example, if the interior nodes are x3 , x9 , x15 , and x21 , then
MESH.int = [3 9 15 21]T .

• MESH.bnd is a ND × 1 vector containing the boundary nodes’ numbers (here, ND denotes
the number of boundary nodes). For example, if the boundary nodes are the x2 , x10 , x19 ,
then
MESH.bnd = [2 10 19]T .
Notice that N = ND + ND and that {MESH.int, MESH.bnd} = {1, . . . , N}.

K
K′ K K′

Figure 1. Example of non conforming triangulations. Intersection of two elements which is not a node (left). In-
tersection of two elements which is not an edge (middle). Hanging nodes of a mesh represented by white dots
(right).

2
 7 8 9
1 .coord .elem .bnd .elbnd .int
1 −1 −1 1 5 4 1 1 1 1 1 3 1 5
7 6
2 0 −1 2 1 2 5 2 2 2 2 1
8 5 3 1 −1 3 5 2 3 3 3 3 3 3

0 4
5
6
4 −1 0 4 3 6 5 4 4 4 4 1
5 0 0 5 5 6 9 5 6 5 5 3
1 4 6 1 0 6 9 8 5 6 7 6 6 1

2 3 7 −1 1 7 5 8 7 7 8 7 7 3
8 0 1 8 7 4 5 8 9 8 8 1
−1
1 2 3 9 1 1

−1 0 1

Figure 2. Mesh T of the square domain D = (−1, 1)2 with M = 8 elements (circled numbers in black) specified by
the MESH.elem matrix and N = 9 vertices (numbers in blue) specified by MESH.coord matrix. Interior and boundary
vertices are stored in MESH.int and MESH.bnd vectors, respectively.

5 4
1

3 4
6 2 2 3 2 3
2
1 3
1 3
1
2 1
1 2 3

Figure 3. Local enumeration of the edges of a conforming mesh T with four elements; numbers in blues are vertices’
global numbers of T whereas circled numbers are elements’ numbers. Coloured numbers denote local elements’
edges. The longest edges of the elements are the second ones. Note that the i-th edge of an element may not be the
the i-th edge of a neighbour element (i = 1, 2, 3). For instance, the edge conv{x5 , x2 } is the 2-nd edge of element K2
as well as the 3-rd edge of element K3 .

• MESH.elbnd is a MD × 2 matrix containing the boundary elements and the (local) num-
ber of the corresponding edge lying on the boundary (here, MD denotes the number of
boundary elements). A boundary element, is an element which has at least one (or at most
two) edge lying on the boundary. For example, if the first 4 boundary elements of a mesh
are K6 , K10 , K22 , and K34 , with the 1-st, 3-th, 2-nd, and 3-th edges on the boundary,
respectively, then
 
 6 1 
10 3 
 
MESH.elbnd = 22 2  .
 
34 3 
 
... ...
 

Figure 2 shows a visual representation of the data structure MESH for the default mesh of
the square domain D = (−1, 1)2 .
Two important features/assumptions of LEBREF2D are the following:

(A1) (enumeration of edges). Given the element K = {xi , x j , xk } ∈ T , with the 1-st node being
xi , the 2-nd node bei ngx j , and 3-rd node being xk , counted counterclockwise, we assume
that

– the 1-st edge of K is the one in front of xi , i.e., conv{x j , xk },

3
 (a) (b) (c) (d)

Figure 4. NVB bisections. (a) One, (b)-(c) two, and (d) three bisections of the edges of an element in T . Double lines
denote the reference edges, black dots denote the edges to be bisected, and white dots denote the newest vertices.
Reference edges are always bisected first.

2 3
4 4
1 2 3 3 2
1 1 2 3 2 3
Figure 5. Similarity classes which occur using iterated longest edge bisections on a two-dimensional element.

– the 2-nd edge of K is the one in front of x j , i.e., conv{xk , xi },

– the 3-rd edge of K is the one in front of xk , i.e., conv{xi , x j }.

(A2) (longest edge). The local enumeration is such that the longest edge of the generic element
K„ is the second one (i.e., conv{xk , xi }).

See Figure 3 for an example.

2 Longest Edge Bisection

The Longest Edge Bisection (LEB) algorithm [5] is a special case of the Newest Vertex Bisection
algorithm [2, 4, 7].
Let K ∈ T„ be an element of the mesh T„ at the „-th iteration (in other words, after 0 refine-
ments). The NVB algorithm chooses a reference edge of K. In the LEB algorithm, the reference
edge is chosen to be the longest edge. Given a set of elements M„ ⊆ T„ to be refined, successive
NVB refinements of T„ are as follows:

(a) for each elements K ∈ M„ , the midpoint xL of the longest edge is connected with the
vertex of K in front of the longest edge. Then, xL becomes a new vertex;

(b) such bisection (for all K ∈ M„ ) produces two new elements (the children of K) whose
reference edges are the edges in front of the newest vertex xL ; see Figure 4(a);

Usually, one iteration of the following steps above will results in hanging nodes. An hanging
node appears in the interior of those edges E = K1 ∩ K2 , with K1 , K2 ∈ T„ , for which the refine-
ment of only one between T1 and T2 involves the bisection of E. Therefore, additional bisections
have to be performed in order to obtain a refined conforming mesh. This is typically referred
to as completion, or mesh-closure, step.
Figure 4(left) shows one iteration of steps (a)-(b). If there is an hanging node of K, we
then assume that the edge on which the hanging node lies has to be bisected together with the
longest edge of K which needs to be bisected first (see Figure 4(middle-left and middle-right)).
If there are two hanging nodes, then all edges have to be bisected, that is, three bisections are

4
 8 8 8
7 7
10 9 10 9

2 1 2 1
6 6 6
3 5 3 5

4 4

(a) (b) (c)

(d) (e)

Figure 6. (a) Initial mesh T with 10 elements. (b) Marked elements M (in orange). (c) Refinement of the marked
elements may introduce hanging nodes. Elements sharing the introduced hanging nodes are also marked for re-
finement. (d) Refinement of the new marked elements. (e) Final refined mesh (total refined elements in green).

performed as shown in Figure 4(right). This way, an element can generate 2, 3, or 4 children
according to the number of edges bisected. Obviously, the procedure terminates after a finite
number of bisections due to the finite number of edges of the mesh. Notice that the longest
edge bisection iterations lead to only 4 similarity classes of elements (see Figure 5).
Figure 6 shows how the LEB iterations are performed in LEBREF2D: (a) It is given an ex-
ample mesh T with 10 elements. (b) Initially 3 elements are marked for refinement. These
elements are M = {K3 , K5 , K9 } ⊆ T . (c) Their longest edges are bisected. This introduces two
hanging nodes, the red midpoints of the edges shared by elements K9 and K8 , and K5 and K6 ,
respectively. In order to eliminate the hanging nodes, elements K8 , K6 are marked for refine-
ment. (d) Refinement of K6 does not lead to further hanging nodes. Instead, the bisection of the
longest edge of K8 introduces one more hanging node on the edge shared by K7 . (e) As before,
also K7 is marked for refinement. Once K7 is refined, no further hanging nodes are introduced.
The mesh-refinement is then terminated. Notice that although #M = 3, the final number of
elements that have been refined is 6.

3 Detail-grid
We first define what a uniform refinement is. We say that Tb is a uniform refinement of a given
mesh T if all edges of T are bisected, i.e., the refinement of each element produces 4 children;
in particular #Tb = 4 · #T .
By detail-grid we refer to as the “mesh” of edge-midpoints that should be introduced, or
added, to the current mesh T in order to obtain the uniformly refined mesh Tb. That is, given
a mesh stored in MESHX, its associated detail-grid is the data structure MESHY containing infor-
mation about all edge-midpoints coordinates, position, and numbers.
The terminology comes from the Finite Element Method. Suppose that X = X(T ) is the

5
 6
1
.coord .elem .bnd
9 7 1 0.50 0.25 1 1 2 8 1 4
0.75
4 2 0.25 0.25 2 4 6 1 2 5

3 5 3 0.25 0.75 3 5 2 3 3 6
0.5 4
4 0.75 0.00 4 3 7 9 4 7

3 5 0.00 0.25 5 8
2
0.25 5 6 6 0.75 0.25 6 9
1 1 2 7 0.25 0.75
8 0.25 0.00
0
1 8 2 4 3 9 0.00 0.75

0 0.25 0.5 0.75 1

Figure 7. A detail grid associated with a given mesh T (in blue) with 4 elements. Representation of the associated
MESHY structure’s fields storing information about all edge-midpoints of T . The detail grid is returned by the
function detailgrid.

space of first-order piecewise-linear basis functions (hat functions) associated with the nodes
of the mesh T and let Y = Y(T ) be the space of first-order piecewise-linear basis functions
associated with the edge midpoints of T (i.e., these hat functions attain the value 1 at edge-
midpoints and 0 at mesh nodes). The space Y is called the detail space (with respect to X). In
particular, the enriched space X b Tb) = X ⊕ Y is associated with the uniform refinement Tb
b = X(
of T .
Figure 7 shows an example of a detail-grid for a given mesh T (in blue) with 4 elements. In
green, it is shown the uniform refinement which would introduces the edge-midpoints. Here,
the refinement is done by 3 longest edge bisection iterations on each element (see Figure 4(d)).
Notice that edge-midpoints have their own enumeration which is independent by the enumer-
ation of the nodes. In particular, this means that the edge-midpoints’ numbers are also the
(global) numbers of the edges of the mesh. Moreover, the edge-midpoints of an element K are
stored so that the j-th local midpoint of K is the one lying on the j-th edge of K, j = 1, 2, 3; see
MESHY.elem in Figure 7 and cf. Figure 3.
The function detailgrid (see Section 5) takes in input a mesh MESHX and returns the asso-
ciated MESHY, which is a data structure composed of three field, coord, elem, and bnd, which
store the coordinates of all midpoints, the element-mapping matrix with respect to the edge,
and the midpoints lying on the boundary, respectively.
The function detailgrid is used by the adaptive refinement routines in order to store
information about those midpoints that will become the new nodes of the refined mesh.

4 Solving PDEs
The local refinement routines can be particularly useful if LEBREF2D is incorporated in, e.g.,
an adaptive Finite Element code for solving Partial Differential Equations (PDEs). In such a
case, at current iteration „, one computes a discrete solution u„ to the problem under consid-
eration on a given underlying mesh T„ of the domain D , and then estimates the local error
contributions from each element of the mesh using some strategy (e.g., some a posteriori error
estimation technique, see [1]). This way, a set of marked elements M„ ⊆ T„ can be determined
and a new mesh T„+1 is obtained by refining this set of marked elements.
Although this is not the primal goal of the package, the routines of LEBREF2D can be easily
adapted and used in the setting of adaptive Finite Elements codes [6, 3].

6
5 List of functions
In the following, we list the functions of LEBREF2D:

• [MESHX] = adjustunstructmesh(MESHX)
If the input MESHX contains an unstructured mesh, the function changes the local order
of the elements in MESHX.elem so that the longest edges of all elements are the second
ones (a generic unstructured mesh may not satisfy this condition which is required by
the refinement routine; see assumption (A2)). More info in help adjustunstructmesh.

• [MESHX.elem] = bisection(MMele,MMedge,markedge,MESHX.elem,MESHY.elem);
Main routine performing the bisection of the marked elements. In particular, MMele is
the overall set of marked elements, including the extra elements which have to be refined
in order to avoid hanging nodes (see Figure 6), and MMedge is the associated set of overall
marked edges which will be bisected.
The output MESHX.elem is the element mapping matrix of the new refined mesh.

• [MESH] = crackdomain(reflev,slith)
Structured crack domain D = (−1, 1)2 \(−1, 0)×{0} grid generation (square domain (−1, 1)2
with a crack along the line (−1, 0) × {0}).
reflev is an optional input which denotes the initial refinement level. This is a non nega-
tive integer (i.e., 0, 1, 2,. . . ) which indicates how many uniform refinements the default
mesh (for reflev=0) consists of. The default mesh for reflev=0 has 9 vertices and 8
elements.
slith denotes instead half of the maximum “height” of the crack. That is, the crack is
created by connecting two points, (−1, −slith) and (−1, +slith) with the origin (0, 0).
Default value slith = 0.01.

• [MESHY,edgelep] = detailgrid(MESHX,ixy)
Linear detail-grid generator. For an input MESHX which contains information about a
mesh T , the function returns the data structure MESHY associated with all edge-midpoints
that would be introduced by a uniform refinement of T ; cf. description in the Section 3
and Figure 7 to see how the data structure MESHY looks like.
The extra output edgelep is the edge-midpoint element-position NE × 4 matrix (here, NE
is the number of edges of the mesh). Suppose that its i-th row is as follows:

(i-th row) 3 4 3 1

This mean that the i-th (global) edge-midpoint lies on the (i-th) edge of the mesh in
common with elements K3 and K4 , and that it lies on the 3-rd local edge of K3 and on the
1-st local edge of K4 . For example, this is the case of the row number 3 for the edgelep
matrix (i.e., the 3-rd edge of the mesh) in the example showed in Figure 7.
Basically, the first two entries (1-st, 2-nd columns) indicate which elements the edge-
midpoints are shared by, while the second two entries (3-rd, 4-th columns) indicate which
local edges of these elements the edge-midpoint lie on.

7
 The full edgelep matrix for the example in Figure 7 is the following:

1 2 1 3
 
1 3 2 2
 
 
3 4 3 1
2 2 1 1
 
 
edgelep = 3 3 1 1

2 2 2 2
 
 
4 4 2 2
 
1 1 3 3
 
4 4 3 3
Notice that for boundary edge-midpoints, the 1-st and 2-nd columns as well as the 3-rd
and 4-th columns are repeated (in the example there are only 3 internal edge-midpoints
and 6 boundary edge-midpoints).
Finally, ixy is an optional input to decide whether or not the coordinates of the edge-
midpoints, i.e. the field MESHY.coord, have to be computed or not.
• [edgenodes] = edge2nodes(MESHX,edgenum,edgelep)
Given in input the edge number edgenum, it returns their two extremal nodes xk , xm , i.e.,
such that the edgenum-th edge of the mesh MESHX is conv{xk , xm }.
edgelep is an optional input to be used when the detail grid is not available; see detail-
grid.
• [MMele,MMedge] = getallmarkelem(Mele,MESHY.elem,edgelep)
Given in input the set Mele of marked elements for refinement, it returns the set of over-
all marked elements to be refined so that no hanging nodes will be introduced during
the mesh refinement step. Then, the output MMele is the union of the set Mele and all
the additional elements containing hanging nodes that would be introduced during the
refinement of Mele (see example in Figure 6).
The second output MMedge is the set of overall marked edges associated with MMele that
have to be bisected.
Note that no real mesh refinement is performed in this function.
• lebdemo(domain)
Demo performing automatic adaptive mesh refinements of 4 available domains accord-
ing to the optional input domain:
– domain = 1 for structured square domain D = (−1, 1)2 ; see squaredomain;
– domain = 2 for structured the L-shaped domain D = (−1, 1)2 \ (−1, 0]2 ;
see lshapedomain;
– domain = 3 for unstructured L-shaped domain D = (−1, 1)2 \ (−1, 0]2 ;
see lshapedomainunstruct;
– domain = 4 for structured crack domain D = (−1, 1)2 \ (−1, 0) × {0}; see crackdomain.
Default domain = 1. Throughout iterations, the number of marked elements are chosen
randomly using a certain threshold parameter markpar ∈ (0, 1] and the MATLAB built-in
function randperm:
Mele = randperm(nel,ceil(markpar*nel))’.
Here, nel is the number of current elements of the mesh. Notice that if markpar = 1, all
elements are marked.

8
• [MESHX,MMele] = lebmeshref(MESHX,Mele)
Main routine for the mesh refinement step. It takes in input the current mesh T„ stored
in MESHX and the set Mele of marked elements M„ ⊆ T„ and returns the new refined mesh
T„+1 as well as the set MMele of overall elements that have been refined from T„ → T„+1 .

• [MESH] = lshapedomainunstruct
Unstructured L-shaped domain D = (−1, 1)2 \ (−1, 0]2 grid generation;

• [MESH] = lshapedomain(reflev)
Structured L-shaped domain D = (−1, 1)2 \ (−1, 0]2 grid generation;
reflev is an optional input which denotes the initial refinement level; cf. description in
crackdomain.

• plotmarkedelem(MESH,Mele,plotelem,plotvtx)
Plots an input mesh and a set of marked elements Mele with colour. Default colour is
orange.
plotelem and plotvtx are optional inputs to plot also the elements’ numbers and the
nodes’ numbers, respectively.

• plotmesh(MESH,titleplot,plotelem,plotvtx)
Plots the input mesh. titleplot is an optional input for the title of the plot.
plotelem and plotvtx are also optional inputs; see description in plotmarkedelem.

• plotmeshxandy(MESHX,MESHY,reftype,plotelem,plotvtx)
Plots the mesh MESHX and its associated detail-grid MESHY.
reftype decides the type of uniform refinements the detail-grid is associated with. They
are the red (reftype = 1) and bisec3 (reftype = 2) refinement.
Figure 8 shows an example. The red uniform refinement refines an element by connecting
the 3 edge-midpoints. This introduces a “central” child in the middle of the element (Fig-
ure 8(a)). The bisec3 uniform refinement is the result of 3 longest edge bisection iterations
on each element (Figure 8(c)). The plots on the right show the resulted uniformly refined
meshes if either red or bisec3 refinement is done (Figure 8(b) and (d), respectively).
Remind that a detail-grid MESHY does not provide information about the type of bisections
to be done. It only store information about midpoints. However, a mesh resulted by a
bisec3 uniform refinement, as in Figure 8(d), is what naturally follows by longest edge
bisection iterations.
plotelem and plotvtx are optional inputs; see description in plotmarkedelem.

• [MESH] = squaredomain(reflev)
Structured square domain D = (−1, 1)2 grid generation.
reflev is an optional input which denotes the initial refinement level; cf. description in
crackdomain.

• [MESHREF] = unimeshref(MESH,reftype,iplot)
Given an input mesh MESH, it returns the uniformly refined mesh MESHREF. The input
reftype decides the refinement-type, red (reftype = 1) or bisec3 (reftype = 2); see
Figure 8.
iplot is an optional input to plot the uniformly refined mesh.

9
 7 12 8 15 9 7 24 8 25 9
27 26 22 21
7 6 29 28 24 19
16 8 6 5 11 18 20 21 22 23
8 5 32 25 23 20
30 31 17 18
4 7 5 4 6 4 17 5 19 6
2 1 15 14
1 4 4 7 9 16
9 1 2 3 14 11 12 14 15 16
2 3 3 8 12 13
5 6 10 11
1 13 2 10 3 1 10 2 13 3
(a) (b)

7 12 8 15 9 7 24 8 25 9
28 27 22 21
7 6 29 26 23 20
16 8 6 5 11 18 20 21 22 23
8 5 30 25 24 19
31 32 17 18
4 7 5 4 6 4 17 5 19 6
2 1 16 15
1 4 3 8 9 14
9 1 2 3 14 11 12 14 15 16
2 3 4 7 10 13
5 6 11 12
1 13 2 10 3 1 10 2 13 3
(c) (d)

Figure 8. Two detail-grids for a given mesh T and associated uniformly refined meshes Tb. Blue numbers denotes
the nodes of T and black numbers the elements. (a)-(b) Detail-grid (edge-midpoints in red) corresponding to a red
uniform refinement; (c)-(d) Detail-grid (edge-midpoints in green) corresponding to a bisec3 uniform refinement.

6 Useful commands
• number of elements and number of nodes of the mesh, respectively:

nel = size(MESH.elem,1);
nvtx = size(MESH.coord,1);

• recover coordinates of internal and boundary nodes, respectively:

xyint = MESH.coord(MESH.int,:);
xybd = MESH.coord(MESH.bnd,:);

• get the nodes of the k-th element:

elnodes = MESH.elem(k,:);

• get all elements sharing the node xk :

find( sum( ismember(MESH.elem, k ), 2) ) ;

• number of edge-midpoints (and then edges) of the mesh (using MESHY):

nedg = size(MESHY.coord,1);

10
 • get the midpoints of the k-th element (using MESHY):

midpel = MESHY.elem(k,:);

• recover the two elements sharing the k-th edge of the mesh (using MESHY):

find( sum( ismember(MESHY.elem, k ), 2) );

• recover the 3 edge-neighbors of the i-th element Ki , i.e., the elements which share at most
one edge with Ki (using MESHY and edgelep):

neigh = edgelep(MESHY.elem(i,:), 1:2);

neigh = neigh(neigh ∼= i );

If Ki is a boundary element, then #neigh = 2 instead of 3.

References
[1] M. Ainsworth and J. T. Oden, A posteriori error estimation in finite element analysis, Pure
and Applied Mathematics (New York), Wiley, 2000.

[2] E. Bänsch, Local mesh refinement in 2 and 3 dimensions, IMPACT Comput. Sci. Engrg., 3
(1991), pp. 181–191.

[3] A. Bespalov and L. Rocchi, Stochastic T-IFISS, February 2019. Available online at
http://web.mat.bham.ac.uk/A.Bespalov/software/index.html#stoch tifiss.

[4] I. Kossaczky, A recursive approach to local mesh refinement in two and three dimensions, J.
Comput. Appl. Math., 55 (1995), pp. 275–288.

[5] M. C. Rivara, Mesh refinement processes based on the generalized bisection of simplices, SIAM
J. Numer. Anal., 21 (1984), pp. 604–613.

[6] D. J. Silvester, A. Bespalov, Q. Liao, and L. Rocchi, Triangular IFISS (T-IFISS) version
1.2. Available online at https://personalpages.manchester.ac.uk/staff/david.silvester/-
ifiss/tifiss.html, December 2018.

[7] R. Stevenson, The completion of locally refined simplicial partitions created by bisection, Math.
Comp., 77 (2008), pp. 227–241.

11