GENERATION OF UNSTRUCTURED MESHES IN 2-D, 3-D, AND
SPHERICAL GEOMETRIES WITH EMBEDDED HIGH
RESOLUTION SUB-REGIONS∗
JORGE M. TARAMÓN† , JASON P. MORGAN† ‡ , CHAO SHI‡ , AND JÖRG

arXiv:1711.06333v1 [cs.CG] 14 Nov 2017

HASENCLEVER§
Abstract. We present 2-D, 3-D, and spherical mesh generators for the Finite Element Method
(FEM) using triangular and tetrahedral elements. The mesh nodes are treated as if they were linked
by virtual springs that obey Hooke’s law. Given the desired length for the springs, the FEM is used to
solve for the optimal nodal positions for the static equilibrium of this spring system. A ’guide-mesh’
approach allows the user to create embedded high resolution sub-regions within a coarser mesh. The
method converges rapidly. For example, in 3-D, the algorithm is able to refine a specific region within
an unstructured tetrahedral spherical shell so that the edge-length factor l0r /l0c = 1/33 within a few
iterations, where l0r and l0c are the desired spring length for elements inside the refined and coarse
regions respectively. One use for this type of mesh is to model regional problems as a fine region
within a global mesh that has no fictitious boundaries, at only a small additional computational cost.
The algorithm also includes routines to locally improve the quality of the mesh and to avoid badly
shaped ’slivers-like’ tetrahedra.
Key words. Finite Element Method, Unstructured tetrahedral mesh, Embedded high resolution
sub-region
AMS subject classifications. 65D18, 68U01, 68W05

1. Introduction. Mesh generation and (adaptive) refinement are essential ingredients for computational modelling in various scientific and industrial fields. A
particular design metric or goal is the quality of the generated mesh, because lowquality meshes can potentially lead to larger numerical approximation errors. A highquality mesh would consist of triangles (in 2-D) or tetrahedra (in 3-D) that have
aspect ratios near 1, i.e. their sides should have similar lengths. The techniques to
generate meshes can be crudely classified into three groups: (1) The advancing front
method [20, 24, 9, 15] starts from the boundary of the domain. New elements are
created one-by-one from an existing front of elements towards the interior until the
region is filled. Advancing front methods generally create high-quality meshes close
to the domain boundaries but can have difficulties in regions where advancing fronts
merge. (2) Octree-based methods [21, 18, 16] produce graded meshes through recursive subdivision of the domain. The simplicity of these methods makes them very
efficient. However, poorly shaped elements can be introduced near region boundaries.
(3) Delaunay Triangulation ensures that the circumcircle/circumsphere associated to
each triangle/tetrahedron does not contain any other point in its interior. This feature makes Delaunay-based methods [7, 23, 8, 26] robust and efficient. However, in
3-D they can generate very poorly shaped tetrahedra with four almost coplanar vertex
nodes. These so-called ’sliver’ elements have a volume near zero. Several techniques to
remove slivers have been proposed [6, 19, 5] although some slivers near the boundaries
can typically persist [13].
Current mesh generation algorithms oriented to engineering such as Gmsh [14],
∗ Submitted

to the editors 14Nov2017.
Funding: This work was funded by the COMPASS Consortium.
† Department of Earth Sciences, Royal Holloway University of London,
(Jorge.TaramonGomez.2014@live.rhul.ac.uk).
‡ EAS Dept., Snee Hall, Cornell University, Ithaca, NY, USA
§ Institute of Geophysics, Hamburg University, Hamburg, Germany
1

Egham,

UK

2

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

GiD (https://www.gidhome.com) or TetGen [28] are based on the methods described
above. Variational methods [1] rely on energy minimization to optimize the mesh
during the generation procedure in order to create higher-quality meshes. A widely
used open access community-code for 2-D mesh generation is Triangle [25], however
there is no 3-D version of this mesh generator. DistMesh [22] is an elegant and simple
spring-based method that allows the user to create 2D and 3D unstructured meshes
based on the distance from any point to the boundary of the domain. However this
algorithm is often slow, requiring many steps to converge,
Any ’good’ mesh should be able to meet the following requirements [3]: (1) It
conforms to the boundary; (2) It is fine enough in those regions where the problem
to be solved demands higher accuracy; (3) Its total number of elements is as small
as possible to reduce the size of the problem and the computational costs to solve
it; (4) It has well-shaped elements to improve the performance of iterative methods
such as the conjugate gradient method [27]. Frequently used mesh generators in 3-D
geodynamic problems are the ones included in the ASPECT [17], Rhea [4] and Fluidity
[11] codes. ASPECT and Rhea are written in C++ with adaptive mesh refinement
(AMR). However their regular hexahedral elements create so-called ”hanging nodes”
in regions where the resolution changes and cannot be directly applied to create wellformed tetrahedral elements. Fluidity is another example of AMR for a tetrahedral
mesh. However it has very limited mesh generation capabilities, and in this context
mesh-generation should not be confused with mesh adaptivity.
Here we present a new unstructured mesh generator that is based on a finite
element implementation of the DistMesh approach using virtual springs between nodes
and solving for the equilibrium positions of the nodes. We modify the Distmesh
solution procedure to directly solve for static equilibrium. Our method is considerably
faster than the DistMesh code. It also allows the user to create tetrahedral meshes
without hanging nodes. The user can also create embedded high resolution sub-regions
within a global coarse mesh. This approach becomes very useful when the goal is to
create a mesh that minimizes the number of fictitious internal boundaries, within a
computational problem.
A key design goal is the generation of a Delaunay mesh using a built-in MATLAB
triangulation function called ’delaunay’. Throughout the algorithm, this ’delaunay’
function is called to generate the spring connectivity matrix that relates nodes to
triangles or tetrahedra. We have also developed and tested techniques for adding or
rejecting nodes in regions where the mesh resolution is too high or too low respectively.
A smooth variation in the element size between high resolution and low resolution
regions is achieved by using a guide-mesh approach. These local operations improve
the quality of the relatively few poorly shaped elements that can result from the
ficticious spring algorithm to determine good nodal locations. The mesh-generation
code is written in vectorized MATLAB, and can be easily used within the MATLAB
working environment.
We will present this approach first in its simplest form for making a mesh in a
well-defined rectangular 2-D region (Section 2). In Section 3 we show how a 2-D
cylindrical annulus mesh can be generated with small modifications to the previous
rectangular mesh generator algorithm. In Section 4 we present the modifications
needed to create the 3-D spherical shell mesh that we are using to solve for mantle
flow.
2. 2-D Rectangular work flow. This mesh generation algorithm has its simplest form as a program to create a 2-D rectangular mesh with an embedded high

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

3

Fig. 1. Flow chart for the mesh generator iterative process. Yellow, orange and green boxes
represent the routines exclusively used for creating 2-D rectangular meshes, 2-D cylindrical annulus
meshes and 3-D spherical shell meshes, respectively. White boxes represent the shared routines to
all mesh generators. µ is the mean of the misfit spring lengths (equation (16)) and q is the quality
factor of the elements (equations (13) and (31) for triangular and tetrahedral elements respectively).
Tolerance parameters µt , qt and q̄t are listed in Table 1.

resolution sub-region. The white and yellow boxes in Figure 1 show the flowchart
that describes this algorithm.
Step 1: Definition of preferred nodal distances and initial placement
of the nodes. The first step in this recipe is to define the preferred nodal distances
within the refined (l0r ) and coarse (l0c ) regions as well as the dimensions of the regions.
In order to avoid poor quality elements, an appropriate smooth transition for the mesh
refinement should be specified. Here we choose a preferred spring-length function that
is defined on a so-called ’guide-mesh’. This approach is very similar to the background
grid approach created by [20]. The generation of a refined rectangular mesh using the
guide-mesh approach involves the following steps. First, create a (coarse) mesh to
serve as a guide-mesh with only a small number of nodes defining the boundaries of
the domain and the internal boundaries of the embedded high resolution and transition
sub-regions. Second, create the design function l0 (x, y) for each node of the guidemesh. This function defines the desired length for the springs around those points.
Third, the function l0 (x, y) is evaluated at the midpoint of all springs using linear
Finite Element shape functions. We find that a coarse guide-mesh is a simple and
flexible way to control nodal spacing during the generation of a Finite Element mesh.
Figure 2a shows the guide-mesh for a rectangular mesh example whose parameters are
listed in Table 1. Red and blue dots represent nodes in the guide-mesh with defined

4

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

Fig. 2. (a) Guide-mesh defined by a few nodes in Cartesian coordinates for a rectangular mesh.
The parameters for this mesh are listed in Table 1. Each node is assigned a value for the desired
spring length, being l0r for red dots and l0c for blue dots. The length of the springs within the refined
region (in red) is approximately equal to l0r . The length of the springs within the transition region
(in green) varies smoothly from l0r to l0c . The length of springs within the coarse region (in blue)
is approximately equal to l0c . (b) Initial guess for the rectangular mesh. (c) Zoom around the left
boundary of the refined region for the initial guess (yellow line in (b)). The guide-mesh defining
refined (red) and transition (green) regions is shown in white dashed lines.

l0r and l0c , respectively. The red region represents the refined region of the mesh
with spring length approximately equal to l0r . The green region defines the transition
region where the length of the springs smoothly varies from l0r to l0c . The blue region
represents the coarse region of the mesh with a apporximate spring length of l0c .
The next step is to create a starting guess for the locations of the nodes. Computational work is reduced considerably with a good initial guess for the density of the
nodes. Nodes on the boundary and within the domain are created taking into account
both the location of the refined region and the desired springs length for elements inside the refined and coarse regions. Boundary nodes in the refined and coarse regions
are created using l0r and l0c respectively for the spacing between the nodes. The
interior nodes within the refined and coarse regions are created using a circle packing
lattice with radius equal to l0r /2 and l0c /2 respectively. This fills each region with
an equilateral triangular tiling. In the transition region the size of the elements is
expected to change smoothly between l0r and l0c . The initial placement for boundary
and interior nodes in the transition region is created using l0r as explained above.
After this step, the rejection method described in [22] is used to discard points and
create a ’balanced’ intitial distribution of nodes. After performing a Delaunay triangulation, a quasi-regular mesh of triangles within the refined and coarse regions, with
a poorly structured transition region between them is created (Figure 2b). Figure 2c
shows a zoom of the initial mesh with the guide-mesh also shown.
Step 2: Spring-based solver. Inspired by [22], to generate an unstructured
mesh we link the future locations of finite element nodes with virtual elastic springs.
The spring length is used to define the desired nodal distance within any mesh region,
i.e., short springs lead to mesh regions with higher resolution and longer springs lead
to lower resolution mesh regions. Nodal positions are solved for so that the global
network of virtual springs is in static equilibrium. The behaviour of each ficticious

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

5

Fig. 3. (a) Virtual spring in the 2-D space. Both global reference system (X, Y ) and local
reference system (X 0 , Y 0 ) are shown. (b) Virtual spring in the 3-D space. Both global reference
system (X, Y , Z) and local reference system (X 0 , Y 0 , Z 0 ) are shown. Grey dots represent two nodes
linked by the virtual spring. Red arrows represent the forces acting at each end of the spring.

spring is described by Hooke’s law
F = −kδs

(1)

where F is the force acting at each end of spring, k is the stiffness of the spring, and
δs is the distance the spring is stretched or compressed from its equilibrium length
l0 . Forces and nodal positions are expressed in x, y coordinates in 2-D (Figure 3a).
Because Hooke’s law is formulated along the spring direction it is necessary to introduce the X 0 axis as the local 1-D reference system to solve for the nodal positions.
Hooke’s law for each spring in the local 1-D reference system is given by
(2a)

f1 0 =

(2b)

f2 0 = −kδs = −k(x2 0 − x1 0 − l0 )

kδs =

k(x2 0 − x1 0 − l0 )

where f 0 and x0 are the force and position of the ends of the spring given by the
subscripts 1 and 2, respectively. Writing equations (2a) and (2b) in matrix form, and
moving the force terms to the left hand side yields
 0 






f1
−1
1
0
−1
1
x1 0
(3)
+k
=k
f2 0
1 −1
l0
1 −1
x2 0
In order to solve for the nodal positions in 2-D, a change from local coordinates (x1 0 ,
0; x2 0 , 0) to global coordinates (x1 , y1 ; x2 , y2 ) is needed. This change of coordinates
is described in matrix form as


cos α sin α
0
0
(4)
R2D =
0
0
cos α sin α
where α is the angle of the X 0 axis measured from the X axis in the counterclockwise
direction (Figure 3a). Applying equation (4) to equation (3) (see section 5 for further
details), equation (3) becomes


 



−cα 2 −sα cα
cα 2
sα cα
x1
f1,x
cα
 −sα cα −sα 2

 



sα cα
sα 2 
  y1  =  f1,y  + kl0  sα 
(5) k 
2
 cα 2






sα cα
−cα
−sα cα
z1
f2,x
−cα 
2
2
sα cα
sα
−sα cα −sα
x2
f2,y
−sα

6

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

Fig. 4. Implementation of boundary conditions along a straight tilted segment (yellow dashed
line) for one triangle. A rotation is needed for the node 2 in order to pass from the global reference
system (X, Y ) to the local reference system (X 0 , Y 0 ) where y2 0 = 0 is the constrained boundary
condition.

where sα ≡ sin α and cα ≡ cos α. Equation (5) can be written in the matrix form as
(6)

Kx = f + fl0

where K is the stiffness matrix, x is the nodal displacement vector, f is the element
force vector and fl0 is the force-term created by the fact that the springs would
have zero-force at their desired length. Because the system of equations is solved
for its equilibrium steady state, f = 0. A vectorized ’blocking’ technique based
on the MATLAB methodology described in the MILAMIN code [10] is employed to
speed up the assembly of the stiffness matrix. The solution to this problem is the
’optimal’ position of each node obtained from the inversion of the system of static
force equilibrium equations
(7)

x = K −1 fl0

Straight line Boundary Conditions. Boundary conditions are necessary to
constrain the mesh to the desired domain boundaries, and to differentiate between
boundary and interior nodes. In the simple case of a rectangular mesh, a boundary
node is free to slide along a domain edges parallel to the X- or Y -axis. We achieve
this by setting one of its yi or xi values to be fixed and letting the other value vary
so that the node is free to move along the boundary segment. In the case of a
general line that is not parallel to the X- or Y -axes, this requires a transformation
from global coordinates to a new local coordinate system in which the constraint
direction is parallel to a local coordinate axis. In other words, the new local axes
have to be parallel to and perpendicular to the boundary segment. For simplicity,
the mathematical implementation is shown for one triangle (Figure 4). Node 2 is
free to slide along the tilted segment (yellow dashed line in Figure 4) since y2 0 = 0
defines the boundary constraint. The boundary condition is imposed by a rotation
of coordinate system for node 2 given by the transformation matrix T that relates

7

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

global coordinates x to local coordinates x0 by

 
x1
1
 y1  
1

 
 x2  
cos α2 − sin α2

 
(8)
 y2  = 
sin α2
cos α2

 
 x3  
y3
| {z } |
{z
x

T










1
1

}|

x1
y1
x2 0
0
x3
y3
{z
x0









}

Applying the transformation matrix to the stiffness matrix and force vector
(9)

K 0 = T T KT

(10)

fl0 0 = T T fl0

the new system of equations is given by
(11)

K 0 x0 = fl0 0

which is solved for x0 . When desired, the original global coordinates are recovered
through the transformation matrix
(12)

x = T x0

Step 3: Mesh refinement. In this algorithm we refine a mesh by decreasing
the element size in the region of interest. One common issue in the refinement process
arises from the size contrast between large and small elements within a short spatial
interval so that poorly-shaped elements with short and long edges may form. In order
to mitigate this issue a transition region surrounding the refined region is defined
using the guide-mesh approach described above (see Figure 2a).
Quality factor for triangles. The ’quality’ of a mesh is determined by assessing the quality of its individual elements. This usually involves measures of angles,
edge lengths, areas (in 2-D), volumes (in 3-D), or the radius of its inscribed and circumscribed circles/spheres, see e.g., [12, 27]. Here we use a normalized quality factor,
which in 2-D is given by
(13)

q2D =

2rc
Rc

where rc is the radius of the element’s inscribed circle and Rc is the radius of its
circumscribed circle. Rc and rc can be expressed as
r
1 (b + c − a)(c + a − b)(a + b − c)
(14)
rc =
2
a+b+c
(15)

abc
Rc = p
(a + b + c)(b + c − a)(c + a − b)(a + b − c)

where a, b and c are the side lengths of the triangle. A fair criteria to evaluate the
quality of a mesh is to provide the minimum and mean values of the quality factor,
cf. [1]. Here both are used as control parameters to determine when the iterative
algorithm has reached the desired mesh quality tolerances (Figure 1).

8

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

Step 4: Local mesh improvements. So far the above algorithm would only
move nodes within the domain to meet the desired spring lengths/internodal distances.
However, in general we do not know a priori how many nodes are needed for a mesh.
Therefore we use algorithms to locally add and remove nodes where the spacing is too
loose or tight in the equilibrium configuration. After solving for nodal positions, we
check if the mesh has reached the expected nodal density by determining the mean
of the misfit in spring lengths (Figure 1). This is given by
(16)

µ=

N
1 X li − l0i
N i=1
l0i

where l is the actual spring length, l0 is the desired spring length and N is the total
number of springs in the mesh. Nodes are added or rejected (see below) if µ ≥ µt .
When µ < µt the expected nodal density is achieved and element shape improvements
(see below) are applied to obtain higher quality elements. After some experimentation
we found it appropriate to use 0.02 < µt < 0.05 for 2-D meshes.
Add/reject nodes. In the iterative process of mesh generation the possibility
to either add or reject nodes plays an important local role. This feature is especially
relevant when the goal is to create a global coarse mesh with an embedded high
resolution sub-region. The logic for adding or rejecting nodes is based on the relative
length change of the springs connecting nodes
(17)

=

l − l0
l0

indicating whether springs are stretched ( > 0) or compressed ( < 0) with respect
to their desired lengths. A new node is created at the midpoint of those springs
with  > 0.5, i.e., springs stretched more than 50% greater than their desired length.
One node at the end of a spring is rejected when  < −0.5, i.e., springs compressed
more than 50% below their desired length. In order to save computational time,
the add/reject nodes routine is called as a sub-iteration within the main iteration in
which nodal positions are found. Sub-iterations are performed until the percentage
of springs with || > 0.5 in the sub-iteration j + 1 is higher than in the sub-iteration
j. This implementation is especially useful when a large fraction of nodes need to be
either added or rejected within a particular region of the mesh, e.g., when a relatively
poor initial guess is used.
Smooth positions of the interior nodes. Good quality meshes are directly
related to the generation of isotropic elements [1]. A Laplacian smoothing criteria, cf.
[9], is used to improve the shape of poorly shaped elements, i.e., to make elements as
close to a equilateral triangles or regular tetrahedra as possible. This method is only
applied to interior nodes. The routine repositions interior nodes towards the mean of
the barycentres of their surrounding elements, i.e.,
N
P

(18)

xs =

xbi

i=1

N

where xs are the new coordinates of the interior node, N is the number of elements
surrounding the interior node and xbi are the barycentre coordinates of the i-th
surrounding element. Figure 13 shows an example of smoothing positions of interior
nodes for a 2-D mesh.

9

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS
Table 1
Mesh Parameters.

Symbol

Meaning

d
l
ri
ro
x0
z0
θ0
φ0
r0

Depth
Length
Inner radius
Outer radius
x-coordinate centre of refined region
z-coordinate centre of refined region
Colatitude centre of refined region
Longitude centre of refined region
Radial distance centre of refined region
Desired spring length for elements
inside the coarse region
Desired spring length for elements
inside the refined region
Transition region depth
Transition region length
Transition region width
Refined region depth
Refined region length
Refined region width
Tolerance for minimum quality factor
Tolerance for mean quality factor
Tolerance for mean misfit spring
length

l0

c

l0

r

dt
lt
wt
dr
lr
wr
qt
q̄t
µt

Rectangular
box

Cylindrical
annulus

Spherical
shell

2900 km
40000 km
0 km
0 km
-

3471 km
6371 km
90◦
6371 km

3471 km
6371 km
90◦
90◦
6371 km

1500 km

2000 km

2000 km

7.5 km

10 km

60 km

2900 km
8000 km
300 km
3333 km
0.45

2900 km
8000 km
300 km
3333 km
0.30

2900 km
6800 km
9600 km
300 km
2200 km
5000 km
0.23

0.89
0.025

0.93
0.04

0.86
0.11

Example: Rectangular mesh with an embedded high resolution region. Several tests have been performed with the above implementations in order to
demonstrate the robustness of this mesh-generation recipe. As an example, we show
the results for a rectangular box with an embedded high-resolution sub-region (code
available in section 5). The input parameters that control the algorithm are listed
in Table 1. The algorithm created the mesh in 9 s (all tests in this study have been
performed using MATLAB R2015a (8.5.0.197613) on a 3.2 GHz Intel Core i5 (MacOSX 10.12.5) with 24 GB of 1600 MHz DDR3 memory) after eight outermost loop
iterations (cf. Figure 1). Figure 5a shows the final mesh (top) and a zoom around
the left boundary of the refined region (bottom) for the iteration 8 (see Figure 14
for iterations 0 (initial mesh) and 1). The final mesh has 22000 nodes forming 43000
triangles (Table 2) with an edge-length factor l0r /l0c = 1/200. The percentage of
triangles within the coarse, transition and refined regions is 0.3%, 6.3% and 93.4%
respectively. The lowest quality factor for an element is 0.51 (red line in Figure 5b)
and the mean quality factor for all elements is 0.99 (blue line in Figure 5b). Only
0.12% of the triangles have a quality factor lower than 0.6 (green line in Figure 5b).
Figure 5c shows the fraction of elements as a function of quality factor for the final
mesh.

10

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

Fig. 5. (a) Final mesh (top) for a rectangular box with an embedded high resolution sub-region
and a zoom around the left boundary of the refined region (bottom). (b) Minimum quality factor (red
line), mean quality factor for all elements (blue line) and percentage of elements having a quality
factor lower than 0.6% (green line) as a function of iteration number. (c) Histogram of the fraction
of elements as a function of quality factor for the final mesh.
Table 2
Information on example meshes.

Mesh
Rectangular
box
Cylindrical
annulus
Spherical
shell

nodes

elements

time (s)

iterations

time per
node

time per
element

22000

43000

9

8

4.1 · 10−4

2.1 · 10−4

12000

23000

17

5

1.4 · 10−3

7.4 · 10−4

27000

150000

224

10

8.3 · 10−3

1.5 · 10−3

3. 2-D Cylindrical annulus work flow. The algorithm presented above needs
to be slightly modified to generate a cylindrical annulus mesh. The white and orange
boxes in Figure 1 show the flowchart that describes this modified algorithm. Since
the general algorithm is the same, in this section we only discuss the parts that differ
from the rectangular mesh generator described previously.
Cylindrical annulus guide-mesh. The generation of a refined cylindrical annulus mesh using the guide-mesh involves the same steps as for a rectangular mesh
except that the function l0 (x, y) becomes l0 (θ, r). In this case the guide-mesh is a
coarse cylindrical annulus mesh defined in polar coordinates. Figure 6a shows the
guide-mesh (white dashed lines) defining the refined (red), transition (green) and
coarse (blue) regions and the parameters are listed in Table 1. Red and blue dots represent l0 r and l0 c respectively. The initial triangulation is shown in black solid lines.
Figure 6c shows a zoom of the guide-mesh defined in polar coordinates. Green dots
represent the points where the function l0 (θ, r) is interpolated. The use of a guide-

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

11

Fig. 6. (a) Guide-mesh (white dashed lines) defined by a few nodes (red and blue dots represent
l0 r and l0 c respectively) in polar coordinates for a cylindrical annulus mesh (initial guess is shown
in black solid lines). Red, green and blue colours represent the refined, transition and coarse regions
respectively. (b) Guide-mesh defined in Cartesian coordinates. Same colours as in (a). (c) Zoom
around an edge of the transition region in polar coordinates. The function l0 (θ, r) can be interpolated
at green dots with maximum precision since both boundaries – the cylindrical annulus mesh and its
guide-mesh – are overlapping. (d) Zoom around an edge of the transition region in Cartesian
coordinates. The function l0 (x, y) cannot be interpolated at magenta dots since they lay outside of
the outer boundary of a Cartesian guide-mesh. The precision of the interpolated l0 values at yellow
dots is reduced since both boundaries – the cylindrical annulus mesh and its guide-mesh – do not
overlap.

mesh defined in polar coordinates (white dashed lines in Figure 6a and Figure 6c)
instead of Cartesian coordinates (white dashed lines in Figure 6b and Figure 6d)
takes advantage of higher precision when l0 values are interpolated in points both
close and on the boundaries (green dots in Figure 6c). This is because the shapes of
the outer and inner boundaries of any cylindrical annulus mesh defined in Cartesian
coordinates is not perfectly circular (Figure 6b). Therefore, it may occur that some
boundary points (magenta dots in Figure 6d) may lay outside of the boundaries of
a Cartesian guide-mesh (which can be a very coarse mesh) preventing accurate interpolation for the desired length at those points. Furthermore, the fact that both
boundaries – the cylindrical annulus mesh and its guide-mesh – would not overlap in
a Cartesian geometry would reduce the precision of the interpolated l0 values (yellow
dots in Figure 6d).

12

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

Fig. 7. (a) Conceptual diagram for circular boundary conditions. The motion of boundary
nodes is first restricted to be along the tangent line to the circle. Then they are ’pulled back’ to the
circle by projecting in the radial direction. (b) Implementation of circular boundary conditions for
one triangle. A rotation is needed for the node 2 in order to pass from the global reference system
(X, Y ) to the local surface-parallel reference system (X 0 , Y 0 ) where y2 0 = |r| is the constrained
boundary condition.

Circular Boundary Conditions. Boundary conditions for a cylindrical annulus mesh are a generalization to the treatment for a straight-sided boundary linesegment. We denote the inner and outer boundaries Σ of the cylindrical annulus
mesh as radii r = rinner and r = router respectively. Ω is the interior region confined
between both boundaries. A useful boundary condition is to prescribe nodes on Σ
that are free to move along the circular boundary. This nodal motion is generated by
two independent steps (Figure 7a): 1) The node is allowed to move along the tangent
line to the circle at its current location, and 2) the node is place onto the circle by projecting its new location in the radial direction. This approximation assumes that the
radial distance needed to put the node back onto the circle is small compared to the
distance moved along the tangent line. For simplicity, the mathematical implementation is presented here only for one triangle (Figure 7b). The boundary condition for
node 2 is that it slides along its tangent line (dashed line in Figure 7b) since y2 0 = |r|,
where r is the radial distance from the centre of the cylindrical annulus mesh to the
boundary. The boundary condition is imposed by a rotation of the coordinate system
for node 2 given by the transformation matrix T that relates global coordinates x
with local coordinates x0 (local surface-parallel reference system (X 0 , Y 0 ) in green in
Figure 7b) by

 


x1
1
x1
 y1  
  y1 
1

 


 x2  
  x2 0 
cos
θ
sin
θ
2
2






(19)
 y2  = 
  |r| 
− sin θ2 cos θ2



 
 x3  
  x3 
1
y3
1
y3
| {z } |
{z
} | {z }
x

T

x0

where θ2 is the angle of the node 2 measured from the Y axis in the clockwise direction.
After applying the transformation matrix to the stiffness matrix and force vector
(20)

K 0 = T T KT

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

13

Fig. 8. (a) Final mesh for a cylindrical annulus with an embedded high resolution sub-region.
(b) Zoom around an edge of the refined region. (c) Minimum quality factor (red line), mean quality
factor for all elements (blue line) and percentage of elements having a quality factor lower than
0.6% (green line) as a function of iteration number. (d) Histogram of the fraction of elements as a
function of quality factor for the final mesh.

(21)

fl0 0 = T T fl0

the new system of equations is given by
(22)

K 0 x0 = fl0 0

which is then solved for x0 . Global coordinates are recovered through the transformation matrix
(23)

x = T x0

Add/reject nodes in cylindrical annulus meshes. The routine to add or
reject nodes for a cylindrical annulus mesh works like the one explained above for
a rectangular mesh. The only difference appears when a new node is added on a
boundary spring. In this case, the new boundary node needs to be projected onto the
surface along the radial direction.
Example: Cylindrical annulus mesh with an embedded high resolution
region. We show the results for a cylindrical annulus mesh with an embedded highresolution sub-region (code available in section 5). The input generation parameters
are listed in Table 1. The algorithm created the mesh in 17 s after 5 iterations.
Figure 8a shows the final mesh (top) and a zoom around an edge of the refined region
(bottom) for iteration 5 (see Figure 15 for iterations 0 (initial mesh) and 1). The final

14

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

Fig. 9. (a) Guide-mesh defined by a few nodes (red and blue dots represent l0r and l0c respectively) in spherical coordinates for a spherical shell. The length of the springs within the refined
region (red) is approximately equal to l0r . The length of the springs within the transition region
(green) smoothly varies from l0r to l0c . Outside the transition region the length of the springs is
approximately equal to l0c . (b) Model domain representing a 3-D spherical shell with an embedded
high resolution sub-region.

mesh has 12000 nodes forming 23000 triangular elements (Table 2) with an edge-length
factor l0r /l0c = 1/200. The percentage of triangles within the coarse, transition and
refined regions is 0.2%, 6.0% and 93.8% respectively. The worst quality factor for an
element is 0.44 (red line in Figure 8b) and the mean quality factor of all elements
is 0.98 (blue line in Figure 8b). Only 0.13% of the triangles have a quality factor
lower than 0.6 (green line in Figure 8b). Figure 8c shows the fraction of elements as
a function of their quality factor for the final mesh.
4. 3-D Spherical shell work flow. The algorithm presented above was developed as an intermediate step towards the generation of 3-D spherical shell meshes that
include an embedded high resolution sub-region. The white and green backgrounds in
Figure 1 show the flowchart that describes the 3-D spherical algorithm. In this section
we discuss those parts of the algorithm that differ from the cylindrical annulus mesh
generator.
Initial placement of the nodes in 3-D. The boundary nodes in the refined and
coarse regions are created by recursively splitting an initial dodecahedron according
to l0r and l0c respectively. This gives a uniform distribution of equilateral triangles
on the spherical surface. In contrast to equilateral triangles in 2-D, which are able to
fill up the plane, regular tetrahedra do not fill up the entire space. However, there
do exist some compact lattices, e.g., the hexagonal close packing (hcp) lattice, that
create a distribution of nodes that leads to well shaped tetrahedra. The interior nodes
within the refined and coarse regions are created by a close-packing of equal spheres
with radii equal to l0r /2 and l0c /2 respectively. The initial placement for boundary
and interior nodes in the transition region is created using l0r as explained above.
Then the rejection method described in [22] is used to discard points and create a
weighted distribution of nodes.
Spherical shell guide-mesh. The generation of a refined spherical shell mesh
using the guide-mesh involves steps similar to those described above except that the
preferred length function l0 (x, y) is now l0 (θ, φ, r). In this case the guide-mesh is a
coarse spherical shell mesh defined in spherical coordinates (Figure 9a).

15

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

Spring-based solver in 3-D. The spring-based solver described above naturally
extends to 3-D. Forces and nodal positions are expressed in x, y and z coordinates
(Figure 3b). In order to solve for nodal positions in 3-D, a change from local coordinates (x1 0 , 0, 0; x2 0 , 0, 0) to global coordinates (x1 , y1 , z1 ; x2 , y2 , z2 ) is needed. This
change of coordinates consists of a 3-D rotation described by the rotation matrix

(24)

R3D =

cos α cos β
0

cos α sin β
0

sin α
0

0
cos α cos β

0
cos α sin β

0
sin α



where α and β are angles equivalents to latitude and longitude, respectively (Figure 3b). Applying equation (24) to equation (3) (see section 5 for details), equation
(3) becomes
−cα 2 sβ cβ −sα cα cβ
−cα 2 cβ 2
−cα 2 sβ 2
−sα cα sβ
 −cα 2 sβ cβ

−s
c
c
−s
−sα 2

α α β
α cα s β
k
s α cα cβ
cα 2 sβ cβ
 cα 2 cβ 2
 c 2s c
s α cα s β
cα 2 s β 2
α β β
s α cα cβ
s α cα s β
sα 2




f1,x
cα cβ
 f1,y 
 cα s β 




sα
 f



=  1,z  + kl0 

 f2,x 
 −cα cβ 
 f
 −c s 

α β
2,y
f2,z
−sα


(25)

cα 2 cβ 2
cα 2 s β cβ
s α cα cβ
−cα 2 cβ 2
−cα 2 sβ cβ
−sα cα cβ

cα 2 sβ cβ
cα 2 s β 2
sα cα sβ
−cα 2 sβ cβ
−cα 2 sβ 2
−sα cα sβ

s α cα cβ
s α cα s β
sα 2
−sα cα cβ
−sα cα sβ
−sα 2









x1
y1
z1
x2
y2
z2









where sα ≡ sin α, cα ≡ cos α, sβ ≡ sin β and cβ ≡ cos β. The system of equations is
solved as described above (see equation (7)).
Spherical Boundary Conditions. For 3-D applications, we currently focus on
developing unstructured spherical meshes. Using a notation similar to that for 2-D
circular boundary conditions, we denote by Σ the inner and outer boundaries of the
spherical shell with radii r = rinner and r = router respectively. Ω is the interior
region between the boundaries. A useful boundary condition consists in prescribing
boundary nodes that are free to slide along the local tangent plane to the spherical
surface. Nodal sliding is generated in two independent steps (Figure 10a): 1) The
node is allowed to move along the local tangent plane to the sphere, and 2) the
node is returned to the sphere’s surface by projecting in the radial direction. This
approximation assumes that the radial distance needed to pull the node back to the
surface of the sphere is small compared to the distance moved along the tangent plane.
For simplicity, the mathematical implementation of the spherical boundary conditions
is presented here only for one tetrahedron (Figure 10b). Node 2 is free to slide along
the tangent plane since the boundary condition is z2 00 = |r|, where r is the radial
distance from the centre of the sphere to the surface. This boundary condition is
imposed by two rotations of the coordinate system for node 2. The first rotation is
around the Z axis by an angle φ2 , which is the longitude of node 2 (local reference
system (X 0 , Y 0 , Z 0 ) in blue in Figure 10b). The second rotation is around the Y 0 axis
by an angle θ2 , which is the colatitude for node 2 (local reference system (X 00 , Y 00 , Z 00 )
in green in Figure 10b). The complete rotation is given by the transformation matrix

16

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

Fig. 10. (a) Conceptual diagram for spherical boundary conditions. The motion of boundary
nodes is first restricted to be along the tangent plane to the sphere. Then, they are ’pulled back’ to
the sphere’s surface by projecting in the radial direction. (b) Implementation of spherical boundary
conditions for one tetrahedron. Two rotations are needed for node 2 to pass from the global reference
system (X, Y , Z) to the local reference system (X 00 , Y 00 , Z 00 ), where z2 00 = |r| is the boundary
condition.

T that relates global coordinates x with local coordinates x00 as follows

(26)

 x
1
 y1
 z
 1
 x
 2
 y2

 z2

 x3

 y3

 z3

 x4

y4
z4
| {z



1





































}

=

x1
  y1
 z
1

  x 00
 2
  y2 00

  |r|

  x3

  y3

  z3

  x4

y4
1
z4
}|
{z


1
1
cos φ2 cos θ2 − sin φ2 cos φ2 sin θ2
sin φ2 cos θ2
cos θ2 sin φ2 sin θ2
− sin θ2
0
cos θ2
1
1
1
1
1
{z

|

x

T

x00



















}

This transformation matrix contains a θ and φ angle for each node on the spherical
boundary. Applying the transformation matrix to stiffness matrix and force vector
(27)

K 00 = T T KT

(28)

fl0 00 = T T fl0

the new system of equations is given by
(29)

K 00 x00 = fl0 00

which is solved for x00 . Global Cartesian coordinates are recovered through the transformation matrix
(30)

x = T x00

Quality factor for tetrahedra. The 3-D quality factor for a tetrahedron is
defined by
(31)

q3D =

3rs
Rs

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

17

Fig. 11. (a) Tetrahedron with vertices OABC. R and r are the radius of the circumscribed
and inscribed spheres respectively. (b) Number of tetrahedra as a function of the quality factor q3D
(green) and the shape measure s (red) for the same mesh.

where rs is the radius of the tetrahedron’s inscribed sphere and Rs is the radius of its
circumscribed sphere. Rs and rs are given by
(32)

(33)

rs =

|a · (b × c) |
(|a × b| + |b × c| + |c × a| + | (a × b) + (b × c) + (c × a) |)

Rs =

|a2 · (b × c) + b2 · (c × a) + c2 · (a × b) |
2|a · (b × c) |

where a, b and c are vectors pointing from one node, O, to the three other nodes
of the tetrahedron A, B and C respectively (Figure 11a). This quality factor is
normalized to be 0 for degenerate tetrahedra and 1 for regular tetrahedra. Note that
different definitions for normalized aspect ratios can lead to different estimators for
the global quality of a mesh. For example, [2] define a shape measure s that depends
on tetrahedral volume and the lengths of its edges. Computing q3D and s for the same
mesh gives differences of up to 0.1 for the worst element (Figure 11b). The quality
factor q3D that we choose to use is a more restrictive aspect ratio than the shape
factor measure s.
Element shape improvements. In 3-D, even when the expected nodal density
is achieved (µ < µt ) by adding or rejecting nodes, a considerable number of poorly
shaped tetrahedra can still persist. Local improvements are needed to ensure that the
mesh is robust enough to perform optimal FEM calculations. After some experimentation, we found it appropriate to use µt = 0.11 although this can vary from 0.1 to
0.2 depending on the degree of mesh refinement. The value of µt for 2-D meshes is
smaller than for 3-D meshes due to the shape compactness that can be achieved on a
2-D planar surface.
Methods based on swapping edges or faces to improve element quality can possibly
generate non-Delaunay triangulations, which will cause problems in algorithms that
rely on a mesh created by a Delaunay triangulation (e.g. point search algorithms).
Hence, as an alternative and in addition to smoothing the position of interior nodes, we
recommend two additional operations to improve the quality of tetrahedral elements.
Improvement of badly shaped tetrahedra. Unstructured 3-D meshes are
composed of irregular tetrahedra. Some may be quite poor in terms of their shape
and quality factor (see [6] for a complete categorization of badly shaped tetrahedra).
The first improvement for tetrahedral shapes acts locally and only modifies one node
of each badly shaped tetrahedron. For each badly shaped tetrahedron, identified by

18

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

q < qbad , where 0.2 ≤ qbad ≤ 0.3, we select the spring with the maximum distortion,
i.e. max(||). If  > 0, a new node is created in the midpoint of the selected spring,
while a node at one end of the selected spring is removed if  < 0. A new connectivity
is then created by another Delaunay triangulation. The new connectivity is only
modified in the surroundings of nodes that have been added or removed, keeping the
rest of the connectivity to be the same as the old triangulation. Figure 16 illustrates a
simple example that improves badly shaped tetrahedra when meshing the unit cube.
Removing slivers. Slivers are degenerate tetrahedra whose vertices are wellspaced and near the equator of their circumsphere, hence their quality factor and
enclosed volume are close to zero. We define a sliver as a tetrahedron with q < 0.1.
Our routine for removing slivers is purely geometrical, i.e., it does not take into
account the actual or desired length of the springs. The four vertices of each sliver
are replaced by the three mesh points of the best potential triangle that can be
generated from all permutations of its vertices and potential new nodes created at the
midpoints of its springs (Figure 17). Delaunay triangulation is called afterwards to
create the connectivity matrix around the changed nodes.
Example: Spherical shell mesh with an embedded high resolution region. We show the results for a spherical shell mesh with an embedded high-resolution
sub-region (code available in section 5). The input mesh parameters are listed in Table 1. We recommend to set the point around which the refined region is created
far from the polar axis since the guide-mesh can have difficulties in interpolating the
desired spring lengths near the polar axis.
For this example, the domain of the mesh is a spherical shell whose boundaries
represent the core-mantle boundary and the Earth’s surface (Figure 9b). The smallest
tetrahedra with quasi-uniform size lie inside the high resolution region (red tesseroid
in Figure 9b). This region is embedded within a coarser global mesh. A transition
region (green tesseroid in Figure 9b) guarantees a gradual change in tetrahedral size
from the high resolution region to the coarse region. The algorithm created the mesh
in 224 s after 10 iterations (see Figure 12a for a cross section of the final mesh).
Figure 12b shows a detail of the mesh around the northern boundary of the refined
region. The mesh has 27000 nodes forming 150000 tetrahedra (Table 2) with an edgelength factor l0r /l0c = 1/33. The fraction of tetrahedra within the coarse, transition
and refined regions is 0.8%, 20.0% and 79.2% respectively (see Figure 18). The worst
quality factor for an element is 0.23 (red line in Figure 12c) and the mean of the quality
factor for all elements is 0.87 (blue line in Figure 12c). Only 1% of the tetrahedra
have a quality factor lower than 0.4. Figure 12d shows the fraction of elements as a
function of their quality factor for the final mesh.
5. Summary. We have developed the tools for generating unstructured meshes
in 2-D, 3-D, and spherical geometries that can contain embedded high resolution
sub-regions. While we do not discuss the recipe for the (simpler) generation of a
Cartesian 3-D mesh, only small modifications to the 3-D spherical code are needed
to assign boundary points to lie along small sets of linear boundary edges and planar
boundary surfaces. The algorithm employs the FEM to solve for the optimal nodal
positions of a spring-like system of preferred nodal positions. Straight line, circular
and spherical boundary conditions are imposed to constrain the shape of the mesh.
We use a guide-mesh approach to smoothly refine the mesh around regions of interest.
Methods for achieving the expected nodal density and improving the element shape
and quality are also introduced to give robustness to the mesh. These allow it to make
Finite Element meshes capable of higher computational accuracy and faster iterative

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

19

Fig. 12. (a) Cross section of the final mesh with an embedded high resolution sub-region after
refinement using the guide-mesh. (b) Zoom around the boundary of the refined region. (c) Minimum
quality factor (red line), mean quality factor for all elements (blue line) and fraction of elements
having a quality factor lower than 0.4% (green line) as a function of iteration number. (d) Histogram
of the fraction of elements as a function of quality factor for the final mesh.

convergence. This approach could also be extended to be used as part or all of an
adaptive mesh refinement routine.
Acknowledgments. We thank Cornell University for supporting the initial work
on this problem by Morgan and Shi, and the COMPASS Consortium for the Ph.D.
support for Jorge M. Taramón.
REFERENCES
[1] P. Alliez, D. Cohen-Steiner, M. Yvinec, and M. Desbrun, Variational tetrahedral meshing,
ACM Trans. Graph., 24 (2005), p. 617, https://doi.org/10.1145/1073204.1073238.
[2] A. Anderson, X. Zheng, and V. Cristini, Adaptive unstructured volume remeshing - I: The
method, J. Comput. Phys., 208 (2005), pp. 616–625, https://doi.org/10.1016/j.jcp.2005.02.
023.
[3] M. Bern, D. Eppstein, and J. Gilbert, Provably good mesh generation, J. Comput. Syst.
Sci., 48 (1994), pp. 384–409, https://doi.org/10.1016/S0022-0000(05)80059-5.
[4] C. Burstedde, O. Ghattas, M. Gurnis, G. Stadler, Eh Tan, T. Tu, L. C. Wilcox, and
S. Zhong, Scalable adaptive mantle convection simulation on petascale supercomputers, in
2008 SC - Int. Conf. High Perform. Comput. Networking, Storage Anal., 2008, pp. 1–15,
https://doi.org/10.1109/SC.2008.5214248.
[5] S. Cheng and T. Dey, Quality meshing with weighted Delaunay refinement, Proc. Thirteen.
Annu. ACM-SIAM Symp. Discret. algorithms, 33 (2002), pp. 137–146, https://doi.org/10.
1137/S0097539703418808.

20

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

[6] S.-W. Cheng, T. K. Dey, H. Edelsbrunner, M. A. Facello, and S.-H. Teng, Silver exudation, J. ACM, 47 (2000), pp. 883–904, https://doi.org/10.1145/355483.355487.
[7] L. P. Chew, Guaranteed-Quality Triangular Meshes, tech. report, Department of Computer
Science, Cornell University, Ithaca, New York, 1989.
[8] L. P. Chew, Guaranteed-Quality Delaunay Meshing in 3D (short version), Proc. Thirteen. Annu. Symp. Comput. Geom., (1997), pp. 391–393, https://doi.org/10.1145/262839.
263018.
[9] W.-Y. Choi, D.-Y. Kwak, I.-H. Son, and Y.-T. Im, Tetrahedral mesh generation based on
advancing front technique and optimization scheme, Int. J. Numer. Methods Eng., 58
(2003), pp. 1857–1872, https://doi.org/10.1002/nme.840.
[10] M. Dabrowski, M. Krotkiewski, and D. W. Schmid, MILAMIN: MATLAB-based finite
element method solver for large problems, Geochem. Geophys. Geosyst., 9 (2008), https:
//doi.org/10.1029/2007GC001719.
[11] D. R. Davies, C. R. Wilson, and S. C. Kramer, Fluidity: A fully unstructured anisotropic
adaptive mesh computational modeling framework for geodynamics, Geochem. Geophys.
Geosyst., 12 (2011), https://doi.org/10.1029/2011GC003551.
[12] J. Dompierre, P. Labbé, F. Guibault, and R. Camarero, Proposal of benchmarks for 3D
unstructured tetrahedral mesh optimization, in 7th Int. Meshing Roundtable, 1998, pp. 525–
537.
[13] H. Edelsbrunner and D. Guoy, An Experimental Study of Sliver Exudation, Eng. Comput.,
18 (2002), pp. 229–240, https://doi.org/10.1007/s003660200020.
[14] C. Geuzaine and J.-F. Remacle, Gmsh: a three-dimensional finite element mesh generator
with built-in pre-and post-processing facilities, Int. J. Numer. Methods Eng. 79(11), 0
(2009), pp. 1309–1331, https://doi.org/10.1002/nme.2579.
[15] Y. Ito, A. M. Shih, and B. K. Soni, Reliable Isotropic Tetrahedral Mesh Generation Based
on an Advancing Front Method, in 13th Int. Meshing Roundtable, 2004, pp. 95–105.
[16] Y. Ito, A. M. Shih, and B. K. Soni, Octree-based reasonable-quality hexahedral mesh generation using a new set of refinement templates, Int. J. Numer. Methods Eng., 77 (2009),
pp. 1809–1833, https://doi.org/10.1002/nme.2470.
[17] M. Kronbichler, T. Heister, and W. Bangerth, High accuracy mantle convection simulation through modern numerical methods, Geophys. J. Int., 191 (2012), pp. 12–29,
https://doi.org/10.1111/j.1365-246X.2012.05609.x.
[18] F. Labelle and J. R. Shewchuk, Isosurface stuffing, ACM Trans. Graph., 26 (2007), p. 57,
https://doi.org/10.1145/1276377.1276448.
[19] X.-Y. Li and S.-H. Teng, Generating well-shaped Delaunay meshed in 3D, in 12th Annu.
ACM-SIAM Symp. Discret. algorithms, Washington, D. C., 2001, pp. 28–37.
[20] R. Löhner and P. Parikh, Generation of three-dimensional unstructured grids by the
advancing-front method, Int. J. Numer. Methods Fluids, 8 (1988), pp. 1135–1149, https:
//doi.org/10.1002/fld.1650081003.
[21] S. A. Mitchell and S. Vavasis, Quality mesh generation in three dimensions, in Proc. Eighth
Annu. Symp. Comput. Geom. ACM, 1992, pp. 212–221, https://doi.org/10.1145/142675.
142720.
[22] P.-O. Persson and G. Strang, A Simple Mesh Generator in MATLAB, SIAM Rev., 46
(2004), pp. 329–345, https://doi.org/10.1137/S0036144503429121.
[23] J. Ruppert, A Delaunay Refinement Algorithm for Quality 2-Dimensional Mesh Generation,
J. Algorithms, 18 (1995), pp. 548–585, https://doi.org/10.1006/jagm.1995.1021.
[24] J. Schöberl, An advancing front 2D/3D-mesh generator based on abstract rules, Comput.
Vis. Sci., 1 (1997), pp. 41–52, https://doi.org/10.1007/s007910050004.
[25] J. R. Shewchuk, Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator, in Appl. Comput. Geom. Towar. Geom. Eng., M. C. Lin and D. Manocha,
eds., vol. 1148, Springer Berlin Heidelberg, 1996, pp. 203–222, https://doi.org/10.1007/
BFb0014497.
[26] J. R. Shewchuk, Tetrahedral mesh generation by Delaunay refinement, in 14th Annu. Symp.
Comput. Geom. - SCG ’98, 1998, pp. 86–95, https://doi.org/10.1145/276884.276894.
[27] J. R. Shewchuk, What is a Good Linear Element? Interpolation, Conditioning, and Quality
Measures, in Elev. Int. Meshing Roundtable, 2002, pp. 115–126, https://doi.org/10.1.1.68.
8538.
[28] H. Si, TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator, AMC Trans. Math.
Softw., 41 (2015), p. 36, https://doi.org/10.1145/2629697.

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

21

SUPPLEMENTARY MATERIALS.
SM1. Derivation of equation (5). The 2-D development of equation (3),
rewritten here for convenience
 0 
h
i

h
i 0 
f1
−1
1
0
−1
1
x1
(34)
+
k
=
k
0
f
1 −1
l
1 −1
x 0
2

0

2

is given by two steps. First, develop the right hand side of equation (34) by writing
local coordinates as a function of global coordinates (see Figure 3a)

k

−1
1



1
−1

x1 0
x2 0



x2 0 − x1 0
−(x2 0 − x1 0 )
 
 
(x2 − x1 )cα + (y2 − y1 )sα 
=k
− (x2 − x1 )cα + (y2 − y1 )sα



−1
1
x1 c α + y1 c α
=k
1 −1
x2 c α + y2 c α


=k

(35)



=k

−1
1



1
−1

cos α
0

sin α
0

0
cos α

0
sin α




x1
 y1 


 x2 
y2

where sα ≡ sin α and cα ≡ cos α. Second, express the global coordinates of the force
vector as a function of local coordinates (see Figure 3a)

 

f1,x
cα
0
 0 
f
s
0
α
1,y

=
 f1 0
(36)
f
0
c
f
α

2,x

f2,y

0

2

sα

Combining equations (34) and (35) gives
(37)


x1
 0 
h
ih
i
h
i

f1
−1
1
cos α sin α
0
0
1
0
 y1  − k −1
=k
f2 0
1 −1
0
0
cos α sin α
x2
1 −1
l0
y2

Substituting equation (37) into equation (36) and reordering gives




cα
0
x1
h
ih
i
s
0
−1
1
cos
α
sin
α
0
0
y
 1 
k  0α c 
1 −1
0
0
cos α sin α
x2
α
0

sα


(38)

=

y2

f1,x
f1,y
f2,x
f2,y





+k

cα
sα
0
0

0
0
cα
sα


h


−1
1

1
−1

i

0
l0



which is equivalent to equation (5).
SM2. Derivation of equation (25). The 3-D development of equation (3),
rewritten here for convenience
 0 
h
i

h
i 0 
f1
−1
1
0
−1
1
x1
(39)
+
k
=
k
0
f
1 −1
l
1 −1
x 0
2

0

2

22

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

also involves two steps. First, develop the right hand side of equation (39) by writing
local coordinates as a function of global coordinates (see Figure 3b)

k

−1
1

1
−1



x1 0
x2 0



=k

= k

(40)


=k


=k


x2 0 − x1 0
0
0
−(x2 − x1 )



(x2 − x1 )cβ + (y2 − y1 )sβ cα + (z2 − z1 )sα



− (x2 − x1 )cβ + (y2 − y1 )sβ cα + (z2 − z1 )sα


−1
1
x1 c α c β + y1 c α s β + z 1 s α
x2 c α c β + y2 c α s β + z 2 s α
1 −1

x1
y1



−1
1
cα cβ cα s β s α 0
0
0  z1

1 −1
0
0
0 c α c β c α s β s α  x2
 y
2
z2












where sα ≡ sin α, cα ≡ cos α, sβ ≡ sin β and cβ ≡ cos β. Second, express the global
coordinates of the force vector as a function of local coordinates (see Figure 3b)
 f
1,x
 f1,y
 f1,z
 f2,x


(41)

f2,y
f2,z

 c c
α β
  cα s β
  sα
= 0
 


0
0

0
0
0
cα cβ
cα s β
sα


 0 
 f1
 f2 0


Combining equations (39) and (40) gives
(42)


f1 0
f2 0

 x
1
i  y1

h
ih
−1
1
cα cβ cα s β s α 0
0
0  z1
=k
1 −1
0
0
0 cα cβ cα s β s α 
 x2
y2
z2



h
i


−1
1
0
−k
1 −1
l0


Substituting equation (42) into equation (41) and reordering gives
 c c

 x
0
α β
1
c
s
0
α
β

h
ih
i  y1
 s
0  −1
1
cα cβ cα s β s α 0
0
0  z1
k  0α c c 
1 −1
0
0
0 cα cβ cα s β s α 
α β 

 x2
0
0

(43)

cα s β
sα

 f
1,x
 f1,y
 f
=  f1,z
 2,x
f2,y
f2,z

y2
z2

 c c
α β

 cα s β

 s
 + k  0α




0
0

0
0
0
cα cβ
cα s β
sα








h
i

 −1
1
0

1 −1
l0


which is equivalent to equation (25).
SM3. Code for Rectangular mesh generation. Code to reproduce the
example shown in Figure 5 and Figure 14

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

23

SM4. Code for Cylindrical annulus mesh generation. Code to reproduce
the example shown in Figure 8 and Figure 15
SM5. Code for Spherical shell mesh generation. Code to reproduce the
example shown in Figure 12 and Figure 18
SM6. Additional Figures.

24

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

Fig. 13. (a) Initial 2-D mesh. (b) Mesh after applying the Laplacian correction to smooth
positions of its interior nodes. Blue points are the barycentres of the triangles. Green and black
crosses are the nodal positions before and after smoothing, respectively. Red arrows indicate the
motions of interior nodes.

Fig. 14. (a) Initial mesh (top) for a rectangular box with an embedded high resolution subregion and a zoom around the left boundary of the refined region (bottom). (b) Mesh (top) and zoom
(bottom) after the first iteration.

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

25

Fig. 15. (a) Initial mesh (top) for a cylindrical annulus with an embedded high resolution subregion and a zoom around an edge of the refined region (bottom). (b) Mesh (top) and zoom (bottom)
after the first iteration.

26

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

Fig. 16. (a) Initial mesh with badly shaped tetrahedra (in blue). Rejected nodes in red. (b)
Badly shaped tetrahedra. (c) Mesh after improving badly shaped tetrahedra contains no badly shaped
tetrahedra. (d) Fraction of tetrahedra for a given quality factor for both before (dashed line) and
after (solid line) local improvements to the shape of badly shaped tetrahedra. The minimum quality
factor for the initial mesh is 0.04 and for the final mesh is 0.39.

MESH GENERATION WITH EMBEDDED HI-RES SUB-REGIONS

27

Fig. 17. Removing a sliver (represented by black lines and dashed grey line for hidden edge).
Possible triangles (grey and green colours) created from permutations of the vertices and midpoints
of the edges of a sliver. Black, red and green points represent unaltered, removed and added nodes,
respectively. qtri is the quality factor for each triangle. The four vertices of the sliver are replaced
by the three mesh points of the potential triangle with the best quality factor (green colour).

28

J. M. TARAMÓN, J. P. MORGAN, C. SHI, AND J. HASENCLEVER

Fig. 18. (a) Tetrahedra within the coarse region. (b) Tetrahedra within the transition region.
(c) Tetrahedra within the refined region.

