1 Introduction
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 timeconsuming than the problem solution itself.
In the case of threedimensional 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 vertexcentered 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.
[10]).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 [7]. One of the pioneering work in this direction is due to Hermeline [5] and to George (see e.g. [3]). Several celebrated work followed included or quoted in publications such as [2] and [4], 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 [6], among many others.
It would be difficult to be exhaustive about the stateoftheart 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 spherelike 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 quasiuniformity or uniform regularity (see e.g.[1]
), is sufficient to guarantee accuracy improvement of the discretization method as the mesh is refined. But more than this, quasiuniformity 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.
[1]). 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 subdomain having a non negligible measure with respect to its own measure. Then any point in the interior of this subdomain 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 threedimensional counterpart of a oneparameter triangular mesh generation procedure studied in [8] for star shaped twodimensional 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. [9]).
Likewise its twodimensional analog [8], 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 oneparameter 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 macrotetrahedra.
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 threedimensional 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 macrotetrahedra, 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 macrotetrahedron 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 macrotetrahedron
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 macrotetrahedron can be written as follows:
(1) 
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 macrotetrahedra with three plane faces and one curved face contained in , where the triple subscript is defined as above. Notice that each one of these macrotetrahedra 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 semiaxis. The third edge of anyone of such trihedra is the bisector of one of the two quadrants formed by the above positive semiaxis and another positive coordinate semiaxis (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
axis .
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) macrotetrahedra 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 macrotetrahedron .
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 macrotetrahedron . 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 onetoone 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 macrotetrahedra 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 twodimensional star shaped domains proposed in [8], 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 byproduct of the numbering of the octants and macrotetrahedra 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 vector
oriented 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:(2) 
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 easytocompute 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.
Remark 1
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 roundoff 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. [11]) we choose both the ratio given by
(3) 
and the minimum of the normalized JoeLiu parameter over all mesh elements . Referring to [11], and denoting by the six edges of for , our normalization consists of taking the square root of the usual value of this parameter (cf. [11]), that is,
(4) 
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 quasiuniformity property of thus generated families of meshes. On the other hand the
metrics based on the JoeLiu parameter can indicate that a family of meshes is (shape) regular in the sense of [1], but by no means whether or not it is
quasiuniform.
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 cigarshaped 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
(5) 
for a parameter . Notice that if the value of the polar radius ranges between and within rather small subdomains 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.
10  

6,000  
0.717640  0.717640  0.717641  0.717640  0.717640  
0.824084  0.819960  0.818755  0.818185  0.817809 
1.0  0.8  0.6  0.4  0.2  0.1  

for  0.717640  0.749023  0.573403  0.390018  0.201354  0.112168  
for  0.717640  0.748662  0.572133  0.387183  0.194082  0.097374  
for  0.824084  0.779466  0.721367  0.659456  0.511977  0.300437  
for  0.817809  0.773979  0.719027  0.659455  0.511977  0.290470 
0.0  0.1  0.2  0.3  0.4  

for  0.717640  0.664304  0.441041  0.287242  0.181349  
for  0.717640  0.663084  0.439356  0.285817  0.180264  
for  0.824084  0.586850  0.298063  0.185235  0.128325  
for  0.817809  0.113584  0.052925  0.032887  0.022714 
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 twodimensional case for the triangulation procedure studied in [8]. Actually in that paper indications are given on how to remedy eventual ”insideout 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 byproduct of our mesh generation procedure is a family of quasiuniform 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 [10], 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 [10] 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/20085.
References
 [1] P.G. Ciarlet. The Finite Element Method for Elliptic Problems. North Holland, Amsterdam, 1978.
 [2] P.J. Frey & P.L George. Mesh generation: Application to Finite Elements, Hermes Science Publishing, Oxford, UK, 2000.
 [3] P.L. George. Automatic Mesh Generation, Wiley, 1991.
 [4] 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), 187218.
 [5] F. Hermeline, Triangulation automatique d’un polyèdre en dimension n. R.A.I.R.O. Numerical Analysis, 163 (1982), 211242.
 [6] 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), 975995.
 [7] D.S.H. Lo, Finite Element Mesh Generation, CRC Press, Taylor & Francis group, 2015.
 [8] V. Ruas. Automatic generation of triangular finite element meshes. Computer and Mathematics with Applications, 5 (1979) 125–140.
 [9] V. Ruas. Numerical Methods for Partial Differential Equations. An Introduction, Wiley, 2016.
 [10] V. Ruas. Methods of arbitrary optimal order with tetrahedral finiteelement meshes forming polyhedral approximations of curved domains, to appear.
 [11] 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.379386, 2007.
Comments
There are no comments yet.