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 2D) or tetrahedra (in 3D) 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 onebyone from an existing front of elements towards the interior until the region is filled. Advancing front methods generally create highquality meshes close to the domain boundaries but can have difficulties in regions where advancing fronts merge. (2) Octreebased 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 Delaunaybased methods [7, 23, 8, 26] robust and efficient. However, in 3D they can generate very poorly shaped tetrahedra with four almost coplanar vertex nodes. These socalled ’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], 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 higherquality meshes. A widely used open access communitycode for 2D mesh generation is Triangle [25], however there is no 3D version of this mesh generator. DistMesh [22] is an elegant and simple springbased 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 wellshaped elements to improve the performance of iterative methods such as the conjugate gradient method [27]. Frequently used mesh generators in 3D 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 socalled ”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 meshgeneration 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 subregions 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 builtin 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 guidemesh 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 meshgeneration 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 welldefined rectangular 2D region (Section 2). In Section 3 we show how a 2D 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 3D spherical shell mesh that we are using to solve for mantle flow.
2 2D Rectangular work flow
This mesh generation algorithm has its simplest form as a program to create a 2D rectangular mesh with an embedded high resolution subregion. The white and yellow boxes in Fig. 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 () and coarse () 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 springlength function that is defined on a socalled ’guidemesh’. This approach is very similar to the background grid approach created by [20]. The generation of a refined rectangular mesh using the guidemesh approach involves the following steps. First, create a (coarse) mesh to serve as a guidemesh with only a small number of nodes defining the boundaries of the domain and the internal boundaries of the embedded high resolution and transition subregions. Second, create the design function , for each node of the guidemesh. This function defines the desired length for the springs around those points. Third, the function , is evaluated at the midpoint of all springs using linear Finite Element shape functions. We find that a coarse guidemesh is a simple and flexible way to control nodal spacing during the generation of a Finite Element mesh. Fig. 2a shows the guidemesh for a rectangular mesh example whose parameters are listed in Table 1. Red and blue dots represent nodes in the guidemesh with defined and , respectively. The red region represents the refined region of the mesh with spring length approximately equal to . The green region defines the transition region where the length of the springs smoothly varies from to . The blue region represents the coarse region of the mesh with a apporximate spring length of .
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 and 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 and 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 and . The initial placement for boundary and interior nodes in the transition region is created using 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 quasiregular mesh of triangles within the refined and coarse regions, with a poorly structured transition region between them is created (Fig. 2b). Fig. 2c shows a zoom of the initial mesh with the guidemesh also shown.
Step 2: Springbased 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 spring is described by Hooke’s law
(1) 
where is the force acting at each end of spring, is the stiffness of the spring, and is the distance the spring is stretched or compressed from its equilibrium length . Forces and nodal positions are expressed in , coordinates in 2D (Fig. 3a). Because Hooke’s law is formulated along the spring direction it is necessary to introduce the axis as the local 1D reference system to solve for the nodal positions. Hooke’s law for each spring in the local 1D reference system is given by equationparentequation
(2a)  
(2b) 
where and are the force and position of the ends of the spring given by the subscripts 1 and 2, respectively.
Writing equations Eqs. 2b and 2a in matrix form, and moving the force terms to the left hand side yields
(3) 
In order to solve for the nodal positions in 2D, a change from local coordinates (, ; , ) to global coordinates (, ; , ) is needed. This change of coordinates is described in matrix form as
(4) 
where is the angle of the axis measured from the axis in the counterclockwise direction (Fig. 3a). Applying equation Eq. 4 to equation Eq. 3 (see SM1. Derivation of equation Eq. 5 for further details), equation Eq. 3 becomes
(5) 
where and . Equation Eq. 5 can be written in the matrix form as
(6) 
where is the stiffness matrix, is the nodal displacement vector, is the element force vector and is the forceterm created by the fact that the springs would have zeroforce at their desired length. Because the system of equations is solved for its equilibrium steady state, . 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) 
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  or axis. We achieve this by setting one of its or 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  or 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 (Fig. 4). Node 2 is free to slide along the tilted segment (yellow dashed line in Fig. 4) since defines the boundary constraint. The boundary condition is imposed by a rotation of coordinate system for node 2 given by the transformation matrix that relates global coordinates to local coordinates by
(8) 
Applying the transformation matrix to the stiffness matrix and force vector
(9) 
(10) 
the new system of equations is given by
(11) 
which is solved for . When desired, the original global coordinates are recovered through the transformation matrix
(12) 
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 poorlyshaped 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 guidemesh approach described above (see Fig. 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 2D), volumes (in 3D), or the radius of its inscribed and circumscribed circles/spheres, see e.g., [12, 27]. Here we use a normalized quality factor, which in 2D is given by
(13) 
where is the radius of the element’s inscribed circle and is the radius of its circumscribed circle. and can be expressed as
(14) 
(15) 
where , and 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 (Fig. 1).
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 (Fig. 1). This is given by
(16) 
where is the actual spring length, is the desired spring length and is the total number of springs in the mesh. Nodes are added or rejected (see below) if . When 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 for 2D 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 subregion. The logic for adding or rejecting nodes is based on the relative length change of the springs connecting nodes
(17) 
indicating whether springs are stretched or compressed with respect to their desired lengths. A new node is created at the midpoint of those springs with , i.e., springs stretched more than 50% greater than their desired length. One node at the end of a spring is rejected when , 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 subiteration within the main iteration in which nodal positions are found. Subiterations are performed until the percentage of springs with in the subiteration is higher than in the subiteration . 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.,
(18) 
where are the new coordinates of the interior node, is the number of elements surrounding the interior node and are the barycentre coordinates of the th surrounding element. Fig. 13 shows an example of smoothing positions of interior nodes for a 2D mesh.
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 meshgeneration recipe. As an example, we show the results for a rectangular box with an embedded highresolution subregion (code available in SM3. Code for Rectangular mesh generation). 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. Fig. 1). Fig. 5a shows the final mesh (top) and a zoom around the left boundary of the refined region (bottom) for the iteration 8 (see Fig. 14 for iterations 0 (initial mesh) and 1). The final mesh has 22000 nodes forming 43000 triangles (Table 2) with an edgelength factor . 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 Fig. 5b) and the mean quality factor for all elements is 0.99 (blue line in Fig. 5b). Only 0.12% of the triangles have a quality factor lower than 0.6 (green line in Fig. 5b). Fig. 5c shows the fraction of elements as a function of quality factor for the final mesh.
Symbol  Meaning 





Depth  2900 km      
Length  40000 km      
Inner radius    3471 km  3471 km  
Outer radius    6371 km  6371 km  
xcoordinate centre of refined region  0 km      
zcoordinate centre of refined region  0 km      
Colatitude centre of refined region    
Longitude centre of refined region      
Radial distance centre of refined region    6371 km  6371 km  
Desired spring length for elements inside the coarse region  1500 km  2000 km  2000 km  
Desired spring length for elements inside the refined region  7.5 km  10 km  60 km  
Transition region depth  2900 km  2900 km  2900 km  
Transition region length  8000 km  8000 km  6800 km  
Transition region width      9600 km  
Refined region depth  300 km  300 km  300 km  
Refined region length  3333 km  3333 km  2200 km  
Refined region width      5000 km  
Tolerance for minimum quality factor  0.45  0.30  0.23  
Tolerance for mean quality factor  0.89  0.93  0.86  
Tolerance for mean misfit spring length  0.025  0.04  0.11 
Mesh  nodes  elements  time (s)  iterations 





22000  43000  9  8  

12000  23000  17  5  

27000  150000  224  10 
3 2D 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 Fig. 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 guidemesh
The generation of a refined cylindrical annulus mesh using the guidemesh involves the same steps as for a rectangular mesh except that the function , becomes , . In this case the guidemesh is a coarse cylindrical annulus mesh defined in polar coordinates. Fig. 6a shows the guidemesh (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 and respectively. The initial triangulation is shown in black solid lines. Fig. 6c shows a zoom of the guidemesh defined in polar coordinates. Green dots represent the points where the function ,
is interpolated. The use of a guidemesh defined in polar coordinates (white dashed lines in
Fig. 6a and Fig. 6c) instead of Cartesian coordinates (white dashed lines in Fig. 6b and Fig. 6d) takes advantage of higher precision when values are interpolated in points both close and on the boundaries (green dots in Fig. 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 (Fig. 6b). Therefore, it may occur that some boundary points (magenta dots in Fig. 6d) may lay outside of the boundaries of a Cartesian guidemesh (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 guidemesh – would not overlap in a Cartesian geometry would reduce the precision of the interpolated values (yellow dots in Fig. 6d).Circular Boundary Conditions
Boundary conditions for a cylindrical annulus mesh are a generalization to the treatment for a straightsided boundary linesegment. We denote the inner and outer boundaries of the cylindrical annulus mesh as radii and 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 (Fig. 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 (Fig. 7b). The boundary condition for node 2 is that it slides along its tangent line (dashed line in Figure 7b) since , where 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 that relates global coordinates with local coordinates (local surfaceparallel reference system (, ) in green in Fig. 7b) by
(19) 
where is the angle of the node 2 measured from the axis in the clockwise direction. After applying the transformation matrix to the stiffness matrix and force vector
(20) 
(21) 
the new system of equations is given by
(22) 
which is then solved for . Global coordinates are recovered through the transformation matrix
(23) 
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 subregion (code available in SM4. Code for Cylindrical annulus mesh generation). The input generation parameters are listed in Table 1. The algorithm created the mesh in 17 s after 5 iterations. Fig. 8a shows the final mesh (top) and a zoom around an edge of the refined region (bottom) for iteration 5 (see Fig. 15 for iterations 0 (initial mesh) and 1). The final mesh has 12000 nodes forming 23000 triangular elements (Table 2) with an edgelength factor . 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 Fig. 8b) and the mean quality factor of all elements is 0.98 (blue line in Fig. 8b). Only 0.13% of the triangles have a quality factor lower than 0.6 (green line in Fig. 8b). Fig. 8c shows the fraction of elements as a function of their quality factor for the final mesh.
4 3D Spherical shell work flow
The algorithm presented above was developed as an intermediate step towards the generation of 3D spherical shell meshes that include an embedded high resolution subregion. The white and green backgrounds in Fig. 1 show the flowchart that describes the 3D 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 3D
The boundary nodes in the refined and coarse regions are created by recursively splitting an initial dodecahedron according to and
respectively. This gives a uniform distribution of equilateral triangles on the spherical surface. In contrast to equilateral triangles in 2D, 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 closepacking of equal spheres with radii equal to
and respectively. The initial placement for boundary and interior nodes in the transition region is created using as explained above. Then the rejection method described in [22] is used to discard points and create a weighted distribution of nodes.Spherical shell guidemesh
The generation of a refined spherical shell mesh using the guidemesh involves steps similar to those described above except that the preferred length function , is now , , . In this case the guidemesh is a coarse spherical shell mesh defined in spherical coordinates (Figure 9a).
Springbased solver in 3D
The springbased solver described above naturally extends to 3D. Forces and nodal positions are expressed in , and coordinates (Fig. 3b). In order to solve for nodal positions in 3D, a change from local coordinates (, , ; , , ) to global coordinates (, , ; , , ) is needed. This change of coordinates consists of a 3D rotation described by the rotation matrix
(24) 
where and are angles equivalents to latitude and longitude, respectively (Fig. 3b). Applying equation Eq. 24 to equation Eq. 3 (see SM2. Derivation of equation Eq. 25 for details), equation Eq. 3 becomes
(25) 
where , , and . The system of equations is solved as described above (see equation Eq. 7).
Spherical Boundary Conditions
For 3D applications, we currently focus on developing unstructured spherical meshes. Using a notation similar to that for 2D circular boundary conditions, we denote by the inner and outer boundaries of the spherical shell with radii and 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 (Fig. 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 (Fig. 10b). Node 2 is free to slide along the tangent plane since the boundary condition is , where 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 axis by an angle , which is the longitude of node 2 (local reference system in blue in Fig. 10b). The second rotation is around the axis by an angle , which is the colatitude for node 2 (local reference system in green in Fig. 10b). The complete rotation is given by the transformation matrix that relates global coordinates with local coordinates as follows
(26) 
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) 
(28) 
the new system of equations is given by
(29) 
which is solved for . Global Cartesian coordinates are recovered through the transformation matrix
(30) 
Quality factor for tetrahedra
The 3D quality factor for a tetrahedron is defined by
(31) 
where is the radius of the tetrahedron’s inscribed sphere and is the radius of its circumscribed sphere. and are given by
(32) 
(33) 
where , and are vectors pointing from one node, O, to the three other nodes of the tetrahedron , and respectively (Figure 11
a). 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 that depends on tetrahedral volume and the lengths of its edges. Computing and for the same mesh gives differences of up to 0.1 for the worst element (Fig. 11b). The quality factor that we choose to use is a more restrictive aspect ratio than the shape factor measure .Element shape improvements
In 3D, even when the expected nodal density is achieved () 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 although this can vary from 0.1 to 0.2 depending on the degree of mesh refinement. The value of for 2D meshes is smaller than for 3D meshes due to the shape compactness that can be achieved on a 2D planar surface.
Methods based on swapping edges or faces to improve element quality can possibly generate nonDelaunay 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 3D 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 , where , we select the spring with the maximum distortion, i.e. . If , 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 . 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. Fig. 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 . 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 (Fig. 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 highresolution subregion (code available in SM5. Code for Spherical shell mesh generation). 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 guidemesh 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 coremantle boundary and the Earth’s surface (Figure 9b). The smallest tetrahedra with quasiuniform size lie inside the high resolution region (red tesseroid in Fig. 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 Fig. 12a for a cross section of the final mesh). Fig. 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 . The fraction of tetrahedra within the coarse, transition and refined regions is 0.8%, 20.0% and 79.2% respectively (see Fig. 18). The worst quality factor for an element is 0.23 (red line in Fig. 12c) and the mean of the quality factor for all elements is 0.87 (blue line in Fig. 12c). Only 1% of the tetrahedra have a quality factor lower than 0.4. Fig. 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 2D, 3D, and spherical geometries that can contain embedded high resolution subregions. While we do not discuss the recipe for the (simpler) generation of a Cartesian 3D mesh, only small modifications to the 3D 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 springlike system of preferred nodal positions. Straight line, circular and spherical boundary conditions are imposed to constrain the shape of the mesh. We use a guidemesh 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 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. Taramón.
References
 [1] P. Alliez, D. CohenSteiner, 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/S00220000(05)800595.
 [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. ACMSIAM Symp. Discret. algorithms, 33 (2002), pp. 137–146, https://doi.org/10.1137/S0097539703418808.
 [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, GuaranteedQuality Triangular Meshes, tech. report, Department of Computer Science, Cornell University, Ithaca, New York, 1989.
 [8] L. P. Chew, GuaranteedQuality 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: MATLABbased 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. Labbé, 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 threedimensional finite element mesh generator with builtin preand postprocessing 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, Octreebased reasonablequality 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.1365246X.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 wellshaped Delaunay meshed in 3D, in 12th Annu. ACMSIAM Symp. Discret. algorithms, Washington, D. C., 2001, pp. 28–37.
 [20] R. Löhner and P. Parikh, Generation of threedimensional unstructured grids by the advancingfront 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 2Dimensional Mesh Generation, J. Algorithms, 18 (1995), pp. 548–585, https://doi.org/10.1006/jagm.1995.1021.
 [24] J. Schöberl, An advancing front 2D/3Dmesh 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 DelaunayBased Quality Tetrahedral Mesh Generator, AMC Trans. Math. Softw., 41 (2015), p. 36, https://doi.org/10.1145/2629697.
Supplementary Materials
SM1. Derivation of equation Eq. 5
The 2D development of equation Eq. 3, rewritten here for convenience
(34) 
is given by two steps. First, develop the right hand side of equation Eq. 34 by writing local coordinates as a function of global coordinates (see Fig. 3a)
(35) 
where and . Second, express the global coordinates of the force vector as a function of local coordinates (see Fig. 3a)
(36) 
Combining equations Eqs. 35 and 34 gives
(37) 
Substituting equation Eq. 37 into equation Eq. 36 and reordering gives
(38) 
which is equivalent to equation Eq. 5.
SM2. Derivation of equation Eq. 25
The 3D development of equation Eq. 3, rewritten here for convenience
(39) 
also involves two steps. First, develop the right hand side of equation (39) by writing local coordinates as a function of global coordinates (see Fig. 3b)
(40) 
where , , and . Second, express the global coordinates of the force vector as a function of local coordinates (see Fig. 3b)
(41) 
Combining equations Eqs. 40 and 39 gives
(42) 
Substituting equation Eq. 42 into equation Eq. 41 and reordering gives
(43) 
which is equivalent to equation Eq. 25.
Comments
There are no comments yet.