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 low-quality meshes can potentially lead to larger numerical approximation errors. A high-quality 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 .
Current mesh generation algorithms oriented to engineering such as Gmsh , GiD (https://www.gidhome.com) or TetGen  are based on the methods described above. Variational methods  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 , however there is no 3-D version of this mesh generator. DistMesh  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 : (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 . Frequently used mesh generators in 3-D geodynamic problems are the ones included in the ASPECT , Rhea  and Fluidity  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 well-formed 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 resolution sub-region. 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 spring-length function that is defined on a so-called ’guide-mesh’. This approach is very similar to the background grid approach created by . 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 , for each node of the guide-mesh. 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 guide-mesh is a simple and flexible way to control nodal spacing during the generation of a Finite Element mesh. Fig. 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 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  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 (Fig. 2b). Fig. 2c shows a zoom of the initial mesh with the guide-mesh also shown.
Step 2: Spring-based solver
Inspired by , 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
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 2-D (Fig. 3a). Because Hooke’s law is formulated along the spring direction it is necessary to introduce the 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 equationparentequation
where and are the force and position of the ends of the spring given by the subscripts 1 and 2, respectively.
In order to solve for the nodal positions in 2-D, a change from local coordinates (, ; , ) to global coordinates (, ; , ) is needed. This change of coordinates is described in matrix form as
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
where and . Equation Eq. 5 can be written in the matrix form as
where is the stiffness matrix, is the nodal displacement vector, is the element force vector and 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, . A vectorized ’blocking’ technique based on the MATLAB methodology described in the MILAMIN code  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
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
Applying the transformation matrix to the stiffness matrix and force vector
the new system of equations is given by
which is solved for . When desired, the original global coordinates are recovered through the transformation matrix
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 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 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
where is the radius of the element’s inscribed circle and is the radius of its circumscribed circle. and can be expressed as
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. . 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
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 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
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 sub-iteration within the main iteration in which nodal positions are found. Sub-iterations are performed until the percentage of springs with in the sub-iteration is higher than in the sub-iteration . 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 . A Laplacian smoothing criteria, cf. , 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.,
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 2-D 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 mesh-generation recipe. As an example, we show the results for a rectangular box with an embedded high-resolution sub-region (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 (18.104.22.168613) 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 edge-length 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.
|Inner radius||-||3471 km||3471 km|
|Outer radius||-||6371 km||6371 km|
|x-coordinate centre of refined region||0 km||-||-|
|z-coordinate 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|
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 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 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 , becomes , . In this case the guide-mesh is a coarse cylindrical annulus mesh defined in polar coordinates. Fig. 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 and respectively. The initial triangulation is shown in black solid lines. Fig. 6c shows a zoom of the guide-mesh defined in polar coordinates. Green dots represent the points where the function ,
is interpolated. The use of a guide-mesh defined in polar coordinates (white dashed lines inFig. 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 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 values (yellow dots in Fig. 6d).
Circular Boundary Conditions
Boundary conditions for a cylindrical annulus mesh are a generalization to the treatment for a straight-sided boundary line-segment. 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 surface-parallel reference system (, ) in green in Fig. 7b) by
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
the new system of equations is given by
which is then solved for . Global coordinates are recovered through the transformation matrix
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 high-resolution sub-region (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 edge-length 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 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 Fig. 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 and
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 toand 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  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 , is now , , . In this case the guide-mesh is a coarse spherical shell mesh defined in spherical coordinates (Figure 9a).
Spring-based solver in 3-D
The spring-based solver described above naturally extends to 3-D. Forces and nodal positions are expressed in , and coordinates (Fig. 3b). In order to solve for nodal positions in 3-D, a change from local coordinates (, , ; , , ) to global coordinates (, , ; , , ) is needed. This change of coordinates consists of a 3-D rotation described by the rotation matrix
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
where , , and . The system of equations is solved as described above (see equation Eq. 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 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
This transformation matrix contains a and angle for each node on the spherical boundary. Applying the transformation matrix to stiffness matrix and force vector
the new system of equations is given by
which is solved for . Global Cartesian coordinates are recovered through the transformation matrix
Quality factor for tetrahedra
The 3-D quality factor for a tetrahedron is defined by
where is the radius of the tetrahedron’s inscribed sphere and is the radius of its circumscribed sphere. and are given by
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, 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 3-D, 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 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  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 well-spaced 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 high-resolution sub-region (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 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 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 edge-length 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.
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 convergence. This approach could also be extended to be used as part or all of an adaptive mesh refinement routine.
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.
-  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.
-  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.
-  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.
-  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.
-  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.
-  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.
-  L. P. Chew, Guaranteed-Quality Triangular Meshes, tech. report, Department of Computer Science, Cornell University, Ithaca, New York, 1989.
-  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.
-  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.
-  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.
-  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.
-  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.
-  H. Edelsbrunner and D. Guoy, An Experimental Study of Sliver Exudation, Eng. Comput., 18 (2002), pp. 229–240, https://doi.org/10.1007/s003660200020.
-  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.
-  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.
-  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.
-  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.
-  F. Labelle and J. R. Shewchuk, Isosurface stuffing, ACM Trans. Graph., 26 (2007), p. 57, https://doi.org/10.1145/1276377.1276448.
-  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.
-  R. Lö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.
-  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.
-  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.
-  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.
-  J. Schö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.
-  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.
-  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.
-  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.
-  H. Si, TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator, AMC Trans. Math. Softw., 41 (2015), p. 36, https://doi.org/10.1145/2629697.
SM1. Derivation of equation Eq. 5
The 2-D development of equation Eq. 3, rewritten here for convenience
where and . Second, express the global coordinates of the force vector as a function of local coordinates (see Fig. 3a)
which is equivalent to equation Eq. 5.
SM2. Derivation of equation Eq. 25
The 3-D development of equation Eq. 3, rewritten here for convenience
where , , and . Second, express the global coordinates of the force vector as a function of local coordinates (see Fig. 3b)
which is equivalent to equation Eq. 25.