In the framework of the numerical solution of boundary value problems by the finite element method or the finite volume method, mesh generation plays a fundamental
role. It has even become a crucial issue in contemporary techniques for numerical simulation such as adaptivity, for which the generation of meshes is sometimes more time-consuming than the problem solution itself.
In the case of three-dimensional problems widespread numerical techniques of the kind are based on partitions of equation’s spatial domain into tetrahedra, by virtue of their flexibility to fit irregular shapes. Moreover the geometry of a tetrahedron conforms very well to simple algebra for both methods. This is the case for instance of linear finite element or vertex-centered finite volume schemes, in which the approximation of a curved domain by a polyhedron equal to the union of mesh tetrahedra does no harm in terms of accuracy. Even in the case of higher order methods this kind of geometrical approximation is acceptable, provided a suitable boundary condition interpolation is employed (see e.g.).
For all those reasons high quality mesh generation is a vast subject, to which an increasing number of respected specialists are steadily contributing. Most devoted themselves to the development of procedures for tetrahedral mesh generation as general as possible. A good survey on this topic can be found in . One of the pioneering work in this direction is due to Hermeline  and to George (see e.g. ). Several celebrated work followed included or quoted in publications such as  and , to name just a few. In the framework of both discretization methods under consideration, the construction of Delaunay tesselations is very important. Therefore several authors contributed in this direction. In this respect we could quote for instance , among many others.
It would be difficult to be exhaustive about the state-of-the-art of tetrahedral mesh generation in a single article. Our point here is that, for obvious reasons, the use of a very general procedure is not so suitable to generate tetrahedral meshes with the best properties in practical terms, when the domains has a more particular shape. If we take the example of a sphere, it is clear that a procedure especially designed for its shape would be preferable to a general one based on a triangular mesh of domain’s surface, like most mesh generators use. The present contribution lies precisely in an extension of such a framework, for we deal here with a very simple procedure to generate high quality tetrahedral meshes of sphere-like domains, or equivalently spheroidal domains. Here quality means that optimal node numbering is achieved without any complex algorithm, for a mesh can be generated by inputting only a single integer parameter defining its degree of refinement. It also means that the mesh tetrahedra have approximately the same shape and volume, as long as the domain is not too distorted as compared to a sphere. Besides the usual stability and consistency requirements, the latter property, known as quasi-uniformity or uniform regularity (see e.g.
), is sufficient to guarantee accuracy improvement of the discretization method as the mesh is refined. But more than this, quasi-uniformity is the condition under which important tools of the mathematical analysis of a numerical method for partial differential equations apply. This is for instance the case of inverse inequalities for Sobolev norms (cf.). Moreover it is always handy to use simple procedures, that nevertheless attempt to distribute mesh elements in such a way that smaller elements are naturally assigned to domain’s narrowest zones. This is the case of the one proposed in this paper, even though only one parameter determines the construction of the partition.
As one can infer from the above introduction, the main limitation of our mesh generation method is the fact that it requires that the domain be not very irregular. More specifically we confine ourselves to the case where its boundary can be expressed in spherical coordinates for a suitable origin located in its interior. Actually in order to guarantee that the already mentioned regularity properties will hold, it is advisable to further require that the domain is star shaped with respect to all points of a sub-domain having a non negligible measure with respect to its own measure. Then any point in the interior of this sub-domain can be taken as the origin of spherical coordinates. For most practical geometries our method is designed for, the best choice of the origin is obvious, such as in the case of a sphere or an ellipsoid.
As we should point out, the method to be described hereafter is a non trivial three-dimensional counterpart of a one-parameter triangular mesh generation procedure studied in  for star shaped two-dimensional domains. Likewise the boundary of the domains it applies to can be expressed in polar coordinates with origin at a suitable point in its interior. Very nice meshes of disks, ellipses, among less classical domains encountered in practical applications have been generated with such a procedure in several author’s work (see e.g. ).
Likewise its two-dimensional analog , the mesh generation method considered in this article is particularly suitable to check the order of a new discretization method, in case the equation to solve is posed in a curved domain. This is because mesh successive refinement is very easy to carry out, and a roughly uniform sequence of meshes is thus generated, as seen in the examples given in the sequel.
An outline of the paper is as follows. In Section 2 we describe our one-parameter tetrahedrization procedure, by defining its vertices together with the way thay are linked together in order to form the final partition. In Section 3 we complete this description by specifying the steps allowing for the practical calculation of the vertex coordinates. In Section 4 we construct meshes of some star shaped domains in order to exemplify our tetrahedrization procedure. In particular we observe numerically mesh quality in terms of both refinement and domain distortion. Finally we conclude in Section 5 with a few remarks.
2 Partition Description
To begin with we consider a modification of the usual partition of a unit cube into tetrahedra, based on its subdivision into macro-tetrahedra.
Taking the origin of a system of cartesian coordinates , , to be the center of , the corresponding axes are chosen parallel
to its edges. Those axes subdivide into eight equal cubes, each one of them corresponding to an octant of the three-dimensional space.
Next we denote these eight cubes and octants by and , respectively, where is triple subscript such that for , being any non zero value of the -th coordinate in .
Now referring to Figure 1 we take as a model octant with for the purpose of this description. In doing so let be the diagonal of which is a half diagonal of , , , be the diagonals of the faces of intercepting at , and , , be the diagonals of the faces of intercepting at the end of . , , , , , subdivide into six equal macro-tetrahedra, which we denote by , where is another triple subscript corresponding to a permutation of . More precisely the position of illustrated in Figure 1 is such that in each point of this macro-tetrahedron we have . Now given an integer parameter , , we subdivide into equal cubes. Next we bring together the vertices of those cubes located in the interior and faces of each macro-tetrahedron by segments parallel to its six edges. In this manner a partition of into equal tetrahedra is generated.
Finally, taking the half diagonals of as a starting point, we proceed in the same manner for the other seven octants using symmetry, thereby generating a partition of the unit cube into equal tetrahedra. Notice that the cartesian coordinates of vertices of all tetrahedra of such a partition are of the form where the s are integers in the interval . Furthermore it is possible to number the vertices of the partition in a structured manner from one through . More precisely, for example, we can number the vertices located
on the face given by one after the other from up to , in the usual way for squares, as shown in Figure 2 for the -th face.
For the later convenience we point out that the coordinates of the vertices belonging to macro-tetrahedron can be written as follows:
Now let be a star shaped domain of with boundary assumed to be of the -class. Such an assumption is not mandatory, and is aimed at simplifying the presentation. is defined by an equation of the form in spherical coordinates with origin conveniently chosen in the interior of . is the radial coordinate, is the azimuthal angle (or longitude) and (or , as some authors prefer) is the polar angle (or colatitude). We shall generate a partition of this domain into tetrahedra by a method
entirely analogous to the one we just described for the unit cube . The idea is to transform cartesian coordinates into spherical coordinates
in a specific way for each one of the trihedra can be subdivided into, corresponding to the tetrahedra .
To begin with, here again we first subdivide into eight octants defined by the cartesian axes, the latter being also associated with the spherical coordinates with the same origin . Akin to the case of the cube we denote by the subset of contained in the octant , and take as a model a partition of with defined as follows:
First of all we observe that is characterized by and , and set and . Next, referring to Figure 3, we subdivide into six disjoint subsets quite abusively called macro-tetrahedra with three plane faces and one curved face contained in , where the triple subscript is defined as above. Notice that each one of these macro-tetrahedra correspond to the intersection with of one of the six trihedra with vertex , having an edge aligned with the line given by and (i.e. the with the segment in Figure 3), a second edge being a positive (cartesian) coordinate semi-axis. The third edge of anyone of such trihedra is the bisector of one of the two quadrants formed by the above positive semi-axis and another positive coordinate semi-axis (in Figure 3 the curved thetrahdra are separated by the curved dashed lines and the segment ).
Now let , and be the three vertices of located on
. Let also be the angular spherical coordinates of for . As a reference we
take and for all and choose to be the vertex located on
Next we consider homothetic transformations of with origin and ratio and let be its boundary, for . For each the vertices of the partition are the points , for integers and with and defined in the following manner :
First of all we set for every . Then for a given and we denote by the angle with vertex at the origin whose edges contain points and different from , respectively. For mere convenience we call the ”left end” and the ”right end” of the angle and denote by the point given by for every . For instance, an illustration of the points , , for and is supplied in Figure 4. Let also and be the intersection with of the polar radii that subdivide the angles and into equal angles in the same plane, respectively, for . These points are numbered from through from angle’s left end to angle’s right end. The points are the intersections with of the polar radii that subdivide into equal angles (in the same plane) numbered from through , from angle’s left end to angle’s right end. An illustration of the points is given in Figure 4 for and .
Finally we construct a partition of the entire domain by application of the principle we have just described in a symmetrically analogous manner to the other seven spatial octants. This means that for each octant we define six (curved) macro-tetrahedra in such a way that the axis contains a straight edge of for each , and a face of the same is contained in the plane . Similarly, and , is the point of whose angular spherical coordinates are given respectively by:
while is the point of located on the axis . Then the vertices of the tetrahedra of the final partition are determined in the same manner as for the macro-tetrahedron .
In Figure 4 the vertices of the partition belonging to with located on , are illustrated for .
Once the vertices of the partition are known there are different possibilities to define the final tetrahedrization of within each curved macro-tetrahedron . We will chose the following one that ensures mesh compatibility on the interfaces of the s, as seen hereafter. First of all we refer to the already described tetrahedrization of the unit cube . Recalling the expressions (1) of the vertex coordinates for that partition, we can immediately establish a one-to-one correspondence between them and the above defined vertices of the intended tetrahedrization of . More specifically this means that corresponds to the point of the unit cube whose cartesian coordinates are given by (1). It follows that, if we assign to the s the same number as its counterpart in we can generate the tetrahedra in the partition of , by simply defining their edges as the segments whose ends carry the same pair of vertex numbers as for the edges of the tetrahedra in the partition of the unit cube.
It is clear that in the above manner we construct a tetrahedrization of consisting of elements. These tetrahedra can obviously be numbered in the same way as for the unit cube , i.e., the number of each tetrahedron in the partition of is the same as the number of an element in the partition of , whenever the numbers of their four vertices coincide.
To conclude we observe that the faces of the tetrahedra contained in the plane interfaces of two contiguous macro-tetrahedra form a triangulation of a plane sector with angle equal to . It turns out that such triangulation coincides with the one constructed by the procedure for two-dimensional star shaped domains proposed in , as illustrated in Figure 5 for .
3 Determining vertex coordinates
It is possible to set up a method for calculating the cartesian coordinates of every vertex of the partition into tetrahedra of a spheroidal domain described in the previous section, given its number , with . This is a simple by-product of the numbering of the octants and macro-tetrahedra advocated therein, together with the integer superscripts , , , which can be associated with the three spherical coordinates as seen below.
First of all we determine the three integers , , with for , that fulfill . In doing so the values of , and are given by
where . Next setting for , is determined by ordering the s in such a way that .
Now all that is left to do is to subdivide the angles in the way described in Section 2, to obtain its cartesian coordinates according to the following recipe: Let and be two points whose angular spherical coordinates are and , respectively, be the measure of and , for two integers and satisfying and . In practice we will have either and or and , for and . Now the components ,,
of the unit vectororiented like the polar radius in the plane of that subdivides this angle into two angles (in the same plane) contiguous to and , with the complementary measures and respectively, satisfy the following equations:
The two first equations express the fact that the point is located on the surfaces of two cones with vertex at the origin, axes and , and apertures equal to and , respectively. Noticing that there is only one point located simultaneously on the surface of the unit ball centered at the origin and on the surfaces of both cones, system (2) has a unique easy-to-compute solution. Finally, once the components of the unit vector in the direction are determined, where is a generic notation for the vertex whose coordinates we are calculating, we can compute associated spherical coordinates . Then from the radial coordinate of given by , we can immediately determine its cartesian coordinates.
In practice it is not necessary to solve (2) as a non linear system of algebraic equations. This is because it necessarily has a unique solution. Therefore, after elimination of two unknowns using the first two equations of (2), we come up with a quadratic equation for the remaining unknown , which may be either , or . Disregarding round-off errors the (unique) solution of this equation must be and we are done.
4 Quality assessment
In this section we assess the quality of the meshes generated by the procedure described in the previous sections for some representative spheroidal domains. More precisely we use two types of data to work this out. The first one allows to check, on the basis of two different metrics, the quality of meshes with a fixed , of domains with decreasing aspect ratios, in such a way that the mesh parameter also remains fixed. A second type of data refers to a given domain, for whose meshes we observe the evolution of the same metrics as above, as increases at the same rate as . Denoting by the volume of a mesh tetrahedron , among the most used metrics (see e.g. ) we choose both the ratio given by
and the minimum of the normalized Joe-Liu parameter over all mesh elements . Referring to , and denoting by the six edges of for , our normalization consists of taking the square root of the usual value of this parameter (cf. ), that is,
Before starting the evaluation of our mesh generation procedure it is important to list some facts about the above metrics.
First of all corresponds to a uniform mesh, such as the one of the unit cube described in Section 2. It should also be noted that and in case is equilateral. The volume ratio criterion will enable us to exhibit (or not) the quasi-uniformity property of thus generated families of meshes. On the other hand the metrics based on the Joe-Liu parameter can indicate that a family of meshes is (shape) regular in the sense of , but by no means whether or not it is quasi-uniform.
4.1 Mesh quality for ellipsoids with decreasing aspect ratios
Here we consider to be the ellipsoid given by . We will take , and more particularly we will let
the value of this parameter decrease in such a way that we will gradually switch from a sphere for , to a cigar-shaped domain with . Owing to symmetry the mesh will be generated only in the octant .
The least to be expected of the procedure being checked is that it engenders nicely regular meshes of a sphere. Thus to begin with we take , and display in Table 1 the evolution of the metrics and as increases. The number of tetrahedra in each mesh equal to is also supplied. The figures clearly indicate that the meshes behave roughly like uniform meshes of a unit cube, for which both metrics are invariant with . Moreover the mesh elements are not so different from each other, since both their volumes and their shapes are rather close, as indicated by and respectively.
As for ellipsoids, we display in Table 2 the evolution of metrics and for equal to , , , , and for two different meshes, namely, for and . Good news here is the low sensitivity to mesh refinement of both metrics. As for the shapes and volumes a steady but relatively moderate deterioration of rates is observed, as the aspect ratio decreases. However this is more than natural, taking into account the significant variation of domain’s shape.
4.2 Mesh quality for domains with pronounced concavities
In these experiments is the domain given by
for a parameter . Notice that if the value of the polar radius ranges between and within rather small sub-domains of . Whatever the case, for , has boundary concavities that become sharper as increases. Akin to Subsection 4.1 and for the same reason, only the octant will be taken into account in the mesh generation process.
In Table 3 we present the same type of results as in Table 2. However, in contrast to the latter case, now indicates a clear degeneracy of element shapes, as the domain becomes more distorted, i.e. as increases. This effect is amplified by the discrepancy between values of this parameter as the mesh is refined. Nevertheless a rather stable behavior of parameter can be observed.
5 Final comments
Some problems may arise when using the mesh generation procedure described in this article, in case the function defining the boundary of the domain has large local Lipschitz constants with respect to the spherical coordinates and . A similar situation may happen in the two-dimensional case for the triangulation procedure studied in . Actually in that paper indications are given on how to remedy eventual ”inside-out turning” of elements, which may occur in such cases. However we will not further elaborate on those issues here, since anyway it is not advisable to mesh too distorted domains using our method.
The unknown numbering issue for a discretization method to be used in connection with the tetrahedrization proposed in this work has been examined as well. As one can easily guess, it is possible to conceive rather simple algorithms for optimal unknown or node numbering using the analogies with the unit cube. However for the sake of brevity we skip details.
As one can easily infer from the description and the examples given in the previous sections, a natural by-product of our mesh generation procedure is a family of quasi-uniform triangulations of the surface of spheroidal domains indexed by the mesh parameter , in case it is not too distorted. This kind of mesh is very useful in shell modeling and in CAD, among other applications.
Local mesh refinement is possible with the procedure proposed in this paper. It suffices to start from a locally refined mesh of the unit cube, and then map the resulting vertex coordinates into the true curved domain in the way prescribed in Sections 2 and 3. However our method is not well adapted to such refinements because after all it only generates structured meshes in a certain sense. This means that local refinement necessarily impacts zones far away from the one where it is necessary, likewise finite difference grids. Moreover the number of mesh elements and nodes must remain constant for a given value of the mesh integer parameter , and therefore local refinement necessarily implies mesh coarsening away from the refined zone. As a corollary, our mesh generation method is unsuitable to adaptivity techniques. Nevertheless it is certainly very useful whenever one is dealing with problems having a smooth solution in a curved domain not so irregular. In this case the user can take the best advantage of these features, by avoiding low quality meshes that might result from general meshing algorithms.
The procedure studied in this paper was first proposed by the author in two papers quoted in , published in the 80’s. One of them written in Portuguese appeared in Revista Brasileira de Computação; the other one was its abridged translation into English published in a conference proceedings. However, to the best of author’s knowledge, this procedure was implemented for the first time in  and had not been the object of any assessment prior to the present work.
Acknowledgment: The author is grateful to CNPq for the financial support through grant 307996/2008-5.
-  P.G. Ciarlet. The Finite Element Method for Elliptic Problems. North Holland, Amsterdam, 1978.
-  P.J. Frey & P.L George. Mesh generation: Application to Finite Elements, Hermes Science Publishing, Oxford, UK, 2000.
-  P.L. George. Automatic Mesh Generation, Wiley, 1991.
-  P.L. George, F. Hecht & E. Saltel. Fully automatic mesh generator for 3D domains of any shape. Impact of Computers in Science and Engineering, 2 (1990), 187-218.
-  F. Hermeline, Triangulation automatique d’un polyèdre en dimension n. R.A.I.R.O. Numerical Analysis, 16-3 (1982), 211-242.
-  F. Hermeline & P.L. George. Delaunay’s Mesh of a Convex Polyhedron in Dimension d. Application to Arbitrary Polyhedra Int. J. Num. Meths. Engin., 33 (1992), 975-995.
-  D.S.H. Lo, Finite Element Mesh Generation, CRC Press, Taylor & Francis group, 2015.
-  V. Ruas. Automatic generation of triangular finite element meshes. Computer and Mathematics with Applications, 5 (1979) 125–140.
-  V. Ruas. Numerical Methods for Partial Differential Equations. An Introduction, Wiley, 2016.
-  V. Ruas. Methods of arbitrary optimal order with tetrahedral finite-element meshes forming polyhedral approximations of curved domains, to appear.
-  Y. Zhang, T.J.R. Hughes & C.L. Bajaj. Automatic 3D Mesh Generation for a Domain with Multiple Materials. Proc. 16th Int. Meshing Roundtable. M.L. Brewer & D. Marcum eds., p.379-386, 2007.