Mesh Generation: Application to Finite Elements

By Pascal Frey and Paul Louis George

The aim of the second edition of this book is to provide a comprehensive survey of the different algorithms and data structures useful for triangulation and meshing construction. In addition, several aspects are given full coverage, such as mesh modification tools, mesh evaluation criteria, mesh optimization, adaptive mesh construction and parallel meshing techniques.
This new edition has been comprehensively updated and also includes a new chapter on mobile or deformable meshes.

• Language: English
• Publisher: Wiley
• Release date: Mar 11, 2013
• ISBN: 9781118623824

Book preview

Chapter 1

General Definitions

Before going further, it seems important to clarify the terminology and to provide some basic definitions together with some notions of general interest. First, we define the covering-up of a bounded domain, then we present the notion of a triangulation before introducing a particular triangulation, namely the well-known Delaunay triangulation.

A domain covering-up simply corresponds to the naive meaning of this word and the term may be taken at face value. On the other hand, a triangulation is a specific covering-up that has certain specific properties. Triangulation problems concern the construction, of a covering-up of the convex hull of a given set of points. In general, a triangulation is a set of simplices, triangles in two dimensions, tetrahedra in three dimensions, with certain properties. If, in addition to a set of vertices, the boundary of a domain (more precisely a discretization of this boundary whose vertices are in the above set) is specified or, simply if any set of required edges (faces) is provided, we encounter a problem of constrained triangulation. In this case, the expected triangulation of the convex hull must contain these required items.

In contrast, the notion of a mesh may now be specified. Given a domain, namely defined by a discretization of its boundary, the problem comes down to constructing a triangulation that accurately matches this specific domain. In a way, we are dealing with a constrained triangulation but, now, we no longer face a convex hull problem and, moreover, the mesh elements are not necessarily simplices.

After having established triangulation and mesh definitions, some other aspects are discussed, including a suitable element definition (as an element is the basic component of both a triangulation and a mesh), finite element definition as well as mesh data structure definition which are the fundamental ingredients of any further processing (such as using a finite element method). In addition, we introduce some definitions related to certain data structures which are widely used in mesh construction and mesh optimization processes. To conclude, we propose measures of mesh quality and of mesh optimality.

Obviously this chapter cannot claim to be exhaustive. In fact, more specific ideas will be introduced and discussed as required throughout the book.

1.1 Covering-up and triangulation

If is a finite set of points in d (d = 2 or d = 3), the convex hull of , denoted as Conv(S), defines a domain Ω in d. Let K be a simplex¹ (triangle or tetrahedron according to d, always considered as a connected and closed set). Then a covering up r of Ω by means of simplices corresponds to the following definition:

Definition 1.1 r is a simplicial covering-up of Ω if the following conditions hold

• (HO) The set of element vertices in r is exactly .

• (H1) , where K is a simplex.

• (H2) The interior of every element K in r is non empty.

• (H3) The intersection of the interior of two elements is an empty set.

Here is a natural definition. With respect to condition (H1) (where while not strictly necessary, we restrict ourselves to simplicial elements), one can see that Ω is the open set corresponding to the domain that means, in particular, that . Condition (H2) is not strictly necessary to define a covering-up, but KeTr it is nevertheless practical with respect to the context and, thus, will be assumed. Condition (H3) means that element overlapping is proscribed.

Similarly, we will consider conforming coverings-up, referred to as triangulations.

Definition 1.2 r is a conforming triangulation or simply a triangulation of Ω if r is a covering-up following Definition (1.1) and if in addition, the following condition holds:

• (H4) the intersection of two elements in r is either reduced to

– the empty set or to

– a vertex, an edge or a face (for d = 3).

More generally, in d dimensions, such an intersection must be a k-face² , for k = − 1,…, d − 1, d being the spatial dimension

Figure 1.1 Conformal triangles (left-hand side) and non-conformal triangles (right-hand side). Note the vertex located on one edge in this case.

Remark 1.1 For the moment, we are not concerned with the existence and possibly uniqueness of such a triangulation for a given set of points. Nevertheless, a theorem of existence will be provided below and, based on some specific assumptions, the particular case of a Delaunay triangulation will be described.

Euler characteristics. The Euler formula, and its extensions, the Dehn-Sommerville relationships, relate the number of k-faces (k = 0, …,d − 1) in a triangulation of Ω. Such formula can be used to check the topological validity of a given mesh or also for other purposes, such as the determination of the genus of a surface.

Definition 1.3 The Euler characteristics of a triangulation r, is the alterned summation:

(1.1)

where nk = 0,.., d denotes the number of the k-faces in the triangulation.

When the triangulation is homotopic to the topological ball, its characteristic is 1. If the triangulation is homeomorphic to the topological sphere, its Euler characteristic is 1 + (−l)d. In two dimensions, the following relation holds:

where nv, ne and nt are respectively the number of vertices, edges and triangles in the triangulation, c corresponds to the number of connected components of the boundary of Ω. More precisely, if the triangulation includes no hole, then nv − ne + nt = 1. In three dimensions, the above formula becomes:

where nf is the number of faces, nt the number of tets and g stands for the genus of the surface (i.e., the number of holes) of the triangulation. Thus, a triangulation of a closed surface is such that nv − ne + nf = 2.

Delaunay triangulation. Among the different possible types of triangulations, the Delaunay triangulation is of great interest. Let us recall that is a set (a cloud) of points (sites) and that Ω is Conv(S), the convex hull of .

Definition 1.4 r is the Delaunay triangulation of Ω if the open discs (balls) circumscribed to any of its elements does not contain any vertex of S.

Figure 1.2 The empty sphere criterion is violated, the disc of K encloses the point P. Similarly, the circumdisc of the triangle with vertex P includes the vertex of triangle K opposite the common edge (the criterion is symmetric for any pair of adjacent elements).

This criterion, the so-called empty sphere criterion or Delaunay criterion, means that all open balls associated with all elements do not contain any vertex, a closed ball containing the vertices of the element under consideration only. This is the main characterization of the Delaunay triangulation. The Delaunay criterion leads to several other characteristics of any Delaunay triangulation. Figure 1.2 shows an example of an element K which does not meet the Delaunay criterion.

A basic theoretical issue follows.

Theorem 1.1 There exists a unique Delaunay triangulation of a set of points.

The proof is evident by involving the duality with the Voronoï’ diagram associated with the set of points (cf. Chapter 7). The existence is then immediate and the uniqueness is achieved as the points are assumed in general position³ if one wishes to have a simplicial triangulation. Otherwise, the following remark holds.

Remark 1.2 In the case of more than three co-circular (resp. four co-spherical) points, a circle (resp. sphere) exists enclosing these points. If the related disk (resp. ball) is empty, the Delaunay triangulation exists but contains non-simplicial elements such as polygons (resp. polyhedra).

Hence, the uniqueness holds if non-simplicial elements are allowed while if the latter are subdivided by means of simplices, several solutions can be found. Nevertheless, while it may be excessive, we will continue to speak of the Delaunay triangulation by observing that all any partitions of a non-simplicial element are equivalent after swapping⁴ a k-face.

A brief digression. The notion of a Voronoï diagram (though it had yet to be called as such!) first appeared in the work of the French philosopher R. Descartes (1596-1650) who introduced this notion in 1644 in his Principia Philosophiae, which aimed to give a mathematical description of the arrangement of matter in the solar system. In 1850, G. Dirichlet (1805-1859) studied this idea in two and three dimensions and this diagram came to be called the Dirichlet tessellation [Dirichlet-1850]. However, its definitive name came after M.G. Voronoï (1868− 1908), who generalized these results in d dimensions [Voronoï–1908].

Nature provides numerous examples of arrangements and quasi-regular paving which bear a strange resemblance to Voronoï diagrams. Figure 1.3 illustrates some of these typical arrangements⁵.

Constrained triangulation. Provided a set of points and, in addition, a set of edges (resp. edges and faces in three dimensions), an important problem is to ensure the existence of these edges (resp. these edges and faces) in a triangulation. In the following, Const denotes a set of such entities.

Definition 1.5 r is a constrained triangulation of Ω for Const if all and any element of Const is an entity of r.

In particular, a constrained triangulation⁶ can satisfy the Delaunay criterion locally, except in some neighborhood of the constraints.

Remark 1.3 As above, provided a set of points and a constraint, we are not concerned here with the existence of a solution triangulation.

Figure 1.3 Top, the wings of a dragonfly (doc. A. LeBéon) show an alveolar structure apparently close to a Voronoï diagram (left-hand side) and one of the more representative examples of regular paving (consisting of hexagonal cells) is that of a bee’s nest (right-hand side). Bottom, two examples of natural arrangements. Left-hand side: the basaltic rock site of the Giant’s Causeway, Co Antrim, Northern Ireland (photo credit: John Hinde Ltd.). Right-hand side, desert region of Atacama (Chile), the drying earth forms patterns close to Voronoï cells.

1.2 Mesh, mesh element, finite element mesh

Now we turn to a different problem. Let Ω be a closed bounded domain in ² or ³. The question is how to construct a conforming triangulation of this domain. Such a triangulation will be referred to as a mesh of Ω and will be denoted by r or h for reasons that will be made clear in the following. Thus,

Definition 1.6 h is a mesh of Ω if

• (H1)

• (H2) The interior of every element K in h is non-empty.

• (H3) The intersection of the interior of two elements is empty.

Condition (H2) is clearly not verified for a beam element for instance. Condition (H3) avoids element overlapping. In contrast to the definition of a triangulation, Condition (HO) is no longer assumed, which means that the vertices are not, in general, given a priori (see hereafter) and, in (H1), the K’s are not necessarily simplices.

Most computational schemes using a mesh as a spatial support assume that this mesh is conforming (although, this property is not strictly necessary for some solution methods).

Definition 1.7 r is a conformal mesh of Ω if Definition (1.6) holds and

• (H4) the intersection of two elements in h is either the empty set, a vertex, an edge or a face (d = 3).

Clearly, the set of definitions related to a triangulation is again met. There is a fundamental difference between a triangulation and a mesh. A triangulation is a covering-up of the convex hull of a given set of points which, in general, is composed of simplicial elements. A mesh is a covering-up of a given domain defined, in most of the applications, via a given discretization of its boundary, this covering-up being composed of possibly non simplicial elements. On the other hand, at least two new problems occur, namely:

• the respect or enforcement, in some sense, of the boundary of the domain so that the triangulation is a constrained triangulation,

• the necessity of constructing the set of points which will define the vertices of the mesh. Usually the boundary points of the given boundary discretization are given as sole input and field points must be explicitly created.

Remark 1.4 For a boundary discretization defining a domain, the existence of a mesh conforming to this discretization holds in two dimensions but is still, at least from a computer point of view, a delicate question in three dimensions.

Remark 1.5 In the finite element method, the meshes⁷ are generally denoted by h, where the index h of the notation refers to the diameters of the elements in the mesh, these quantities being used in error bound theorems.

As previously mentioned, a mesh can be composed of elements of different geometric natures. A mesh consists of a finite number of segments in one dimension, segments, triangles and quadrilaterals (quads for short) in two dimensions and the above elements, tetrahedra (tets), pentahedra and hexahedra (hexes) in three dimensions. The mesh elements must generally satisfy some specific properties depending on the application involved.

Meshes can be classified into three main classes according to their connectivity

Definition 1.8 The connectivity of a mesh is the definition of the connection between its vertices.

Then, following this definition

Definition 1.9 A mesh is called structured (resp. unstructured) if its connectivity is of the finite difference type (resp. any other type).

A structured mesh can be termed as a grid⁸. In two dimensions, a grid element is a quadrilateral while, in three dimensions, a grid consists of hexahedra. The connectivity between nodes is of the type (i, j, k), i.e., assuming the indices of a given node, the node with indices (i, j, k) has the node with indices ((i − 1), j, k) as its left neighbor and that with indices ((i + 1), j, k) as its right neighbor; this kind of mesh is convenient for geometries for which such properties are suitable, i.e., for generalized quadrilateral or hexahedral configurations.

Remark 1.6 Peculiar meshes other than quad or hex meshes could have a structured connectivity. For instance, one can consider a classical grid of quads where each of them are subdivided into two triangles using the same subdivision pattern.

Such a mesh is usually composed of triangles (tetrahedra) but can also be a set of quadrilaterals (hexahedra) or, more generally, a combination of elements of a different geometric nature. Note that quad or hex unstructured meshes are such that the internal vertices may be shared by more than 4 (8) elements (unlike the case of structured meshes).

For completeness, we introduce two more definitions.

Definition 1.10 A mesh is said to be mixed if it includes some elements of a different geometric nature.

Definition 1.11 A mesh is said to be hybrid if it includes some elements with a different spatial dimension.

A mixed mesh, in two dimensions, is composed of triangles and quads. A hybrid mesh, again in two dimensions, is clearly a mixed mesh but, for instance, includes some triangles together with some segments.

To complete this classification, a mesh may be manifold or not. This point concerns only surface meshes.

Definition 1.12 A (conformal) surface mesh is called manifold if its internal edges are shared by exactly two elements or only one element in the case of a boundary edge for an open surface.

Otherwise, the surface mesh is said to be non-manifold. This is the case of surface meshes which include stiffeners or which have two or more connected components.

Mesh element

The elements are the basic components of a mesh. An element is defined by its geometric nature (triangle, quadrilateral, etc.) and a list of vertices. This list, enriched with some conventions (see hereafter), allows the complete definition of an element, including the definition of its edges and faces (in three dimensions).

Definition 1.13 The connectivity of a mesh element is the definition of the connections between the vertices at the element level.

This connectivity, the local equivalent of the mesh connectivity, makes the description of the topology of the element possible

Definition 1.14 The topology of a mesh element is a definition of the relationships between its faces, edges and vertices.

Triangle connectivity and topology. For convenient purposes, the (local) numbering of vertices and edges is pre-defined in such a way that some properties are implicitly induced⁹. This definition is only a convention leading to implicit properties. In particular, a ordered numbering of the vertices enables us to compute the surface area of a triangle with a positive, or directional, sense. It also allows us to evaluate directional normals for each edge.

In the case of a triangle with connectivity [1,2,3], the first vertex (1) having been chosen, the numbering of the others is deduced counterclockwise. Then the topology can be well defined by means of the edge definition:

– edge [1] runs from vertex (1) to vertex (2),

– edge [2] : (2) → (3),

– edge [3] : (3) → (1),

or alternatively,

– edge [1] is opposite vertex (1), it runs from vertex (2) to vertex (3),

– edge [2] : (3) → (1),

– edge [3] : (1) → (2).

Once a topology has been chosen, all mesh elements must conform to this rule. Such an implicit definition will be a source of simplicity hereafter, avoiding explicit definitions at the element level during the computational step, as mentioned earlier.

Usual element connectivities and topologies. Elements other than triangles are now defined in terms of the two above definitions.

• The segment: [1,2], (1)→ (2).

• The quadrilateral: [1, 2, 3, 4] with a numbering as for the triangle,

edge [1] : (1) → (2) edge [2] : (2) → (3)

edge [3] : (3) → (4) edge [4] : (4) → (1)

Figure 1.4 Local vertex numbering of segment, triangle and quadrilateral, given the first vertex index.

Figure 1.5: Tetrahedron, pentahedron and hexahedron

• The tetrahedron¹⁰ : [1,2,3,4] with assumed to be positive with, for the edges:

and, for the faces:

• The pentahedron: [1,2,3,4,5,6] with assumed to be positive, with, for the edges:

and, for the faces:

• The hexahedron: [1, 2, 3, 4, 5, 6, 7, 8], assumed to be positive, with

and, for the faces:

Other types such as pyramid may be defined. Actually, this type of element allows for some flexibility in mixed meshes. This is the case when structured hex meshes must be combined with unstructured tet meshes.

Finite element mesh

So far, we have considered meshes as geometric entities. Now we turn to the notion of finite element meshes since we are mainly interested in finite element computation as the mesh has a great deal of importance for this application. As will be discussed in Chapter 20, finite elements will be constructed based on the mesh element. To this end, it will be necessary to properly define the nodes, degrees of freedom, interpolation schemes, etc. so as to define the required structures (stiffness matrix, right-hand side, etc.) in order to compute the solution to the problem at hand.

Let us briefly recall for those not so familiar with a finite element style computation, that the classical scheme of such calculus includes the following steps (for the simple case of a linear system to solve):

• a definition of the computational domain,

• a mesh construction step whose purpose is to complete the list of the (geometric) elements of this mesh,

• an interpolation step which constructs the finite elements from the mesh elements,

• a matrix and right-hand side construction step to complete the system corresponding to the discretization of the initial equations, based on the element connectivity,

• a solution step which computes the solution of the above system.

This being clarified, we can proceed by giving some definitions related to the finite element meshes.

Node definition. The finite elements will be associated with the mesh elements at the computational step. For the moment, a finite element is a geometric element supplied with a list of nodes.

Definition 1.15 A node is a point supporting one or several unknowns or degrees of freedom (dof).

The nodes are defined according to the interpolation used in the computation. For a given geometric element, several finite elements may be exhibited as a function of the interpolation step. The simplest finite element is the Lagrange P¹ finite element whose nodes are the element vertices. A Lagrange P² finite element includes as nodes the element’s vertices and a point on each of its edges (in the usual case, these nodes are the edge midpoints). Other finite elements may involve several nodes for each edge, nodes located on faces or inside the element while the element vertices may be nodes or not (see Chapter 20).

Once the node location has been established, it is a case of defining a local numbering for the nodes of the finite element. This task is trivial when the only element vertices are the nodes as the node numbering follows the vertex numbering. If the nodes are defined elsewhere, the local numbering must be well defined. It can be either implicitly defined as for the vertex or explicitly defined in some cases where an implicit definition is not possible, for instance, when the number of nodes varies from one edge to the other as it does for some finite elements. If we consider a Lagrangian P² triangle, it is common to define the three first nodes as the three vertices and then to define as fourth node the node located along the first edge of the triangle (and so on for the other nodes). See Chapter 20) and Chapter 17) for node (re)numbering issues.

Physical attributes. At the solution step, the finite elements are the support of various computations or specific treatments. The mesh, through its elements must contain information making it possible the selection of a set of elements, a set of faces, set of edges or a set of nodes, in other words, making possible any processing concerning these entities, in particular to take into account the loads, flows, pressures, boundary conditions, graphic requirements at the time the solutions must be displayed, etc., related to the problem considered. This will allow to carry out the adequate assignment (i.e., associate the relevant value with such or such an item) or the proper computation of the useful integrals over particular entities.

It is then convenient to associate a physical attribute with all the mesh entities (elements, faces, edges, nodes). This task can be carried on in many manners and, at the computer level, can be implemented in various different ways.

Geometric attributes. For similar reasons, a geometric attribute (provided in some way) proves to be useful for some operations such as the proper definition of the nodes in the case of a finite element other than P¹ or the definition of a subdomain. Thus, it is also convenient to associate a geometric attribute with all the element entities (faces, edges, nodes).

1.3 Mesh data structures

A data structure is simply, but not only, a way of organizing a set of given information (values). A mesh data structure is a way to store all the values that will be useful for a further processing. Following this idea, two categories of structures can be distinguished. One is typically used when constructing the mesh. In other words, such a structure, referred to as the internal mesh data structure, is the organization of the values describing the mesh which are required inside the mesh generation method chosen for this task. On the other hand, the structure referred to as the output mesh data structure is the mesh organization used at the computational step, i.e., outside the mesher. Thus, the values stored in it as well as the way in which they are organized can slightly differ from those in the previous structure. The possible differences between the two mesh structures depend both on the type of the mesh generation method and the solution method that are used. However, for the output structure, it could be desirable to have a universal structure or, at least, a (presumably) generic structure, more suitable for data exchanges.

The internal mesh data structure

The key idea is to define as simple a structure as possible which is well suited to the problem. In some sense, it is the minimal amount of information needed for the application. Various reasons can justify this policy. Among them, a given mesh generation method, due to the algorithm used and the nature of the meshes it is capable of processing, may require a given data organization that is different from another method. Hence, a universal structure is not, a priori, a good solution, as it would be unnecessarily complex in certain cases. Thus, for a given method, one has to find what is needed specifically and what is not strictly required. In this respect, efficiency as well as memory occupation reasons can be invoqued to justify such or such a choice.

The internal mesh structure is only used within the algorithm and, when the mesh is created, the information stored within this structure is transformed to complete the output data structure which is the natural and unique link with the other computational steps.

Some values and data items are specific to one structure or the other. Some others must be included in both. The point is to make this requirement precise and to define the structure accordingly. In the next chapter, we will give some indication about what such an internal structure could be.

The output mesh data structure

Defining a suitable, general purpose and reasonably simple output data structure for mesh storage is not a trivial task. A number of issues must be addressed in order to fulfill various requirements. These include:

• the definition of the nodes, when defining the finite elements at the so-called interpolation step,

• the definition of boundary conditions and loads, the computation and the assembly of the relevant matrices and right-hand sides,

• the solution of the resulting system(s),

• the visualization of the results

• and many others types of processing related to the nature of the problem to solve.

In this respect, it is clear that a classic computational scheme involving a mesh generation step, an interpolation step, the computation of the system to be solved and the solution step is less demanding than an adaptive scheme requiring a loop of such operations where, for instance, the mesh (and the geometry itself, in some cases) is modified at each iteration step of the loop. Thus, a mesh data structure must be defined in such a way as to provide easy access to:

• the vertex coordinates

• the element vertices,

• the physical attributes of the mesh entities,

• the visualization of the results

• the geometrical attributes of these entities,

in order to make the previous computational requirements possible.

These basic principles being stated, it is beyond the scope of this book to discuss further what the ideal data structure could be. Nevertheless, it must be observed that, at this time, existing norms do not give a satisfactory answer to the question of knowing what a mesh data structure should be. On the other hand, in Chapter 2 and Chapter 20, one can gain some indications about data structures when seen from a more abstract point of view or from a purely practical point of view.

The .mesh data structure

Flexibility and versatility have been the major concerns when designing the following mesh data structure, proposed in [George, Borouchaki-1997], Chapter 10. With no surprise, this mesh data structure is named with the suffix .mesh.

Notations. The terms in policy font are file items. The blanks, <> and tabs are item separators. The comments start with the character "#" and end at the end of the line, unless if they are in a string. The comments are placed between the fields.

The notation ( …, i=l,n ) stands for an implicit DO loop.

The syntactic entities are field names, integer values (I) , (double) floating values (R) , strings (C∗) (up to 1024 characters) being placed between . The blanks and < < new lines>> are significant when used between quotes and to use a quote in a string, one has to type it twice "

Booleans (B) : 0 for false and any other value for true (1 in general). Numbers, for instance, a vertex number, is denoted by @Vertex.

The entities of number type (assuming that the numbering starts from 1) are the vertex numbers @Vertex, the edge numbers @Edge, the triangle numbers @Tria, the quadrilateral numbers @Quad, the tetrahedron numbers @Tetra, the pentahedron numbers @Penta, the hexahedron numbers @Hexa, the numbers of a vertex in the appropriate support (described later), @Vertex supp , the numbers of a support edge @Edge supp , the numbers of a support triangle @Tria supp , the numbers of a support quadrilateral @Quad supp , the numbers of a support tetrahedron @Tetra supp , the numbers of a support pentahedron @Penta supp , the numbers of a support hexahedron @Hexa supp .

In addition, Refϕi denotes a number based on a physical attribute.

Description in extenso. The data structure first includes a string identifying the release and then the various fields that can be used.

• IncludeFile (C*) filename

• BoundingBox ( (R) Mini (R) Maxi , i=1, dim )

• End

Remark 1.7 In the following, we give some comments concerning the different fields. At first, it may be seen that some fields are mandatory while some others are optional¹¹.

The comments and remarks about the fields are given according to their introduction order in the above description. Info contained in the data structure are supposed to be sufficient to make any processing possible, if the mesh tool processing the data has not the relevant capability, the corresponding info is simply ignored.

The string MeshVersionFormatted indicates the release identificator and the type of the file. MeshVersionUnf ormatted is an alternative case for this field.

The edge table, Edges, includes only, a priori, the edges with a significant reference number Refϕ).

The elements are given with respect to their geometric nature (triangle, quadrilateral, etc.). In this way, when several types of elements coexist in the mesh, it is not required to manage a table of element types.

The sub-domains are defined using one edge in two dimensions or one face in three dimensions combined with the orientation information, (Orientationi), indicating on which side of this entity the sub-domain lies. The sub-domain number is Refϕs.

A corner point, Corners (for a support type structure), is a point where there is only a C° continuity between the edges sharing the point. Thus, a corner is necessarily an existing mesh vertex, listed in the Vertices list.

A ridge is an edge where there is a C° continuity between the faces sharing it. Thus, a ridge is necessarily a mesh edge, listed in the Edges list.

The required vertices, RequiredVertices, are the vertices of the support that must be present in the mesh as element vertices. Similarly, some edges or (triangular or quadrilateral) faces can be tagged as required entities.

The tangent vector to an edge, Tangent At Edges, gives the tangent vector (with respect to the surface) for this edge at the indicated endpoint. Giving the tangent vector of an edge by means of the tangent vector at a point enables us to deal with the case where several edges (boundary lines) emanate from a point.

The normal at a vertex, NormalAtVert ices, gives the (unit) vector normal the surface at this vertex.

The normal at a vertex of a triangle (a face), NormalAtTrianglesVertices, gives the normal vector at the vertex of the specified triangle. The difference between the normal at vertices and the normal at a triangle vertex allows for local discontinuities between neighbouring triangles.

The corner threshold, AngleOf CornerBound, is a value which enables us to determine the continuity type between two edges or two faces that was not clearly defined or not explicitely specified

The mesh vertices are related to the type of support in which they exist. There are two categories of supports, a geometric support and a current mesh

When the support is of a geometric nature, Geometry, and is defined by a file, it gives the relationships between the vertices, boundary edges and boundary faces of the current mesh with the geometric entities. Thus, a mesh vertex can be identical to a geometric vertex, a mesh edge can have a geometric edge as support and, in three dimensions, a face can have a geometric face as support. These relationships allow us to classify the entities of the current mesh with respect to an entity defining the domain geometry, this information will be particularly useful when constructing finite elements of higher order.

When the support is a (usual) mesh by itself, MeshSupport Of Vert ices, and is defined by a file, it gives the relationships between the current mesh and the above mesh. A vertex of the current mesh belongs to an entity¹² of the support mesh. This information may be relevant when interpolating or transferring a solution from one mesh to another, in an adaptive iterative process for instance.

Hence, in an iterative computational process, the support for the mesh at a given iteration step is the mesh of the previous step. In this way, we indicate that a vertex, i, of the current mesh

• is identical with a vertex of the support

• lies on an edge of the support at abcissa u,

• falls within a triangle of the support, u, v being the coordinates in the reference element,

• falls within a quadrilateral of the support, u, v (idem),

• falls within a tetrahedron of the support, u,v,w (idem),

• falls within a pentahedron of the support, u,v,w (idem),

• falls within an hexahedron of the support, u,v,w (idem).

A vertex not in this table is considered to be a free vertex. The relationships defined in this way enable us to know the location of a vertex using the reference element related to the support entity which includes this vertex. Using the reference element to come to the current element, requires using one of the following relations according to the geometric type of the element (see Figures 1.4 and Figure 1.5 for the numbering convention):

• for an edge with endpoints k1 and k2

• for a triangle with vertices kl, l = 1,3

• for a quad with vertices kl, l = 1,4

• for a tetrahedron with vertices kl, l = 1,4

• for a pentahedron with vertices kl, l = 1,6

• for an hex with vertices kl, l = 1,8

Remark 1.8 In principle, this information is naturally known by the mesh generation algorithm and is relatively easy to obtain. Moreover, when simplicial elements are used, the barycentric coordinates are trivial to obtain and thus do not strictly need to be stored.

Crack definition is the purpose of three fields, CrackedEdges, CrackedTriangles and CrackedQuadrilaterals; we specify then that an edge (a face, respectively) is identical in terms of geometry to another edge (face).

The three fields, EquivalentEdges, EquivalentTriangles and EquivalentQuadrilaterals indicate that two edges (resp. faces) must be meshed the same way (for instance, in periodic meshes).

A comment about the meaning of the physical reference numbers is provided in the field PhysicsRef erence.

It is possible to include a file in the data structure, IncludeFile. This inclusion will be made without ensuring any compatibility.

For some applications, it is useful to know the size of the domain, i.e., the extrema of its point coordinates. This will be stored in the field BoundingBox

The string End indicates the end of the data struture wherever it is encountered.

To conclude with this description, one must consider the data structure together with a suite of vizualisation and manipulation (reading, writing, converting) tools to make the things coherent and consistent altogether.

1.4 Control space and neighborhood space

Control space

It is useful to introduce the notion of a control space for various purposes, as will be seen in the chapters devoted to adapted meshes. Indeed, this space serves to determine the current background.

For the sake of simplicity, an ideal control space is simply the specification of a function H(x, y) defined at any point P(x , y) of ² in a two-dimensional problem. In other words, function H is defined analytically and is used to specify the size and the directional features that must be conformed to by the mesh elements anywhere in the space. However, from a practical point of view, the control spaces will not be ideal in the sense that above function H is actually provided only in a discrete manner. Indeed, formally speaking, a control space can be defined as follows, [George-1991]:

Definition 1.16 (Δ, H) is a control space for the mesh T of a given domain Ω if:

• Ω ⊆ Δ where Δ covers the domain Ω

• a function H(P, ) is associated with every point P ∈ Δ, where is the direction of the disc ¹ (or the sphere ² in three dimensions):

Thus a control space includes two related ingredients. First a covering triangulation Δ, is defined. Next, a function H is defined, whose support is a covering triangulation A. This pair allows for the specification of some properties or criteria to which the elements of the mesh should conform.

In terms of geometry, Δ is an arbitrary covering triangulation. For example, it can be one of the following types:

type 1: a regular partition, such as a finite difference type,

type 2: a quadtree or octree-based partition (see Chapters 2 and 5),

type 3: an arbitrary user-constructed mesh or the current mesh, for instance, in an iterative process, the last computed mesh which then serves as a background mesh.

In addition to this partitioning aspect, (Δ, H) contains, by means of the function H, the global information related to different aspects. One could be the geometry of the domain, the other could be the physics of the problem. These values allow us to determine whether the mesh T, under construction, conforms to the function everywhere.

To construct H, one can consider one of the following approaches

• compute, from the data, the local stepsizes h (the value h being the desired distance between two points) related to the given points. A generalized interpolation then enables us to obtain H. This process is purely geometric in the sense that it relies on geometric data properties: boundary edge lengths, etc.;

• manually define the value of H for every element of Δ. A desired size is then given everywhere in the space for isotropic control, or the desired sizes according to specific directions are given for anisotropic control;

• manually specify H by giving its value for each element of the covering triangulation constructed for this purpose (we return here to a type 3 space as introduced above);

• use the cell sizes (in the above type 2 case), where this size is used to encode the value of H. This then leads to the construction of the (Δ, H) space so as to satisfy this requirement;

• define H from an a posteriori error estimate. We are then in an iterative adaptive process. A mesh T is constructed, the corresponding solution is computed and the error estimate analyzes this solution so as to complete H. Then Δ is set to T and the pair (T, H) forms the control space used to govern the new mesh construction (cf. Chapter 21).

For each of the different types, this definition results in one or the other control space types. In what follows, we will show how to create the internal points in accordance to the specifications contained in this space (for a mesh construction method) or to optimize a mesh (for a mesh optimization method).

Definition 1.17 When the geometric locus of points P + H(P, ) is a circle (resp. a sphere in three dimensions), with P in Δ and d varying, the control space is isotropic. When this locus is an ellipse (resp. an ellipsoid), the control space is anisotropic.

Only these two cases will be discussed hereafter, leading to the definition of the metric maps which are used to govern the mesh construction process.

Remark 1.9 The popular notion of a background mesh, extensively used for instance, in adaptive mesh construction processes, is nothing more than a control space. In this case the covering triangulation is the mesh of the previous iteration step while the function results from the analysis of the current solution by means of an a posteriori error estimate.

Remark 1.10 The above function H is strictly related to a metric tensor field as will be seen later. In this context the discrete meaning of H(P, ) will be seen as a d × d symmetric definite matrix denoted, at point P, by M(P).

Note that some authors, while using different approaches, return in fact to this notion of a control space. For instance, source points can be introduced and this information is used to complete an equivalent of the above H function, [Lohner-1989].

Neighborhood space

Close to the previous idea of control space is the notion of a neighborhood space. Such a space acts to help any neighboring information or queries. For instance, this is a way to facilitate the search for all the vertices located in some regions as such a query could be of great interest for various purposes

The neighborhood space structure is similar to that of a control space. It includes two components, a spatial support together with some information related to this support. The spatial support may be defined as above with, in this case, a special emphasis on simple structures (grid, tree structure such as a binary tree or a quadtree in two dimensions) in which localization problems are evident. Along with this support, the space includes some information of a topological nature. For instance, for a grid, one encodes in its cells the fact that a point, an edge, etc., exists or not. Thus, when such a cell is queried, we know whether a point, an edge, etc., exists within a certain neighborhood. This is a simple example of what could be encoded in this kind of space.

1.5 Mesh quality and mesh optimality

The purpose of constructing a mesh is not simply a question of creating a mesh (a covering-up) of the domain of interest but to obtain a good quality mesh. Therefore, the immediate question is: what is a good quality mesh? or, similarly: how is it possible to define an optimality criterion?.

In response to these questions, the literature contains a number of tentative definitions which are more or less naive (and, in some cases, fanciful to the point of eccentricity !). In this respect, let us mention some widely held, although not necessarily true, ideas. For instance, a nice looking mesh is a good mesh. This raises the problem of what nice looking actually means. A good mesh is a mesh whose elements have the same size and conform to a nice aspect ratio. In addition, some peremptory assertions can be found such as a Delaunay mesh is optimal. The list goes on, but our aim here is to find some reasonable criteria that are well suited to qualifying a mesh.

A preliminary and obvious remark helps us in this discussion. The mesh is constructed for a specific purpose: to solve a given problem. Therefore, the true quality or optimality problem is related to the solution that can be computed with the mesh as a support. In this respect, it makes sense to claim that the mesh quality is good if the resulting solution quality is good. As a consequence and following the same approach, optimality is obtained if the mesh size is minimal (in some sense, for instance, if the number of nodes and vertices is minimal) resulting in a minimal cost when computing the problem solution. This characterization is then related to the nature of the problem and, moreover, implicitly assumes that a suitable way to qualify the solution quality is given.

In a numerical simulation by means of finite elements, an error estimate is the sole judge. A nice mesh is one which leads to the best possible numerical accuracy, i.e., to a minimal bound for the approximation error. For other cases, for instance, in a surface visualization problem based on a mesh of the surface of interest, the quality must be measured with regard to the approximation quality or the graphic aspect of the surface mesh with respect to the real surface. For other problems, the quality aspect may be related to different objectives.

At first, this computational process is based on an approximation of domain Ω where the problem is formulated. Thus, a first condition regarding the mesh quality is to make this approximation precise. In fact, we want to solve the given problem in domain Ω and not in an approximate domain different from the real domain. This being established (that results in conditions about the boundary elements of the mesh leading to a good boundary approximation), the mesh quality is related to the solution and thus to the nature of the problem under investigation. Using an error estimate enables us to see the quality of the mesh and, more precisely, to know if its elements conform to a desirable size, a nice directional aspect.

When such an appreciation is not available, i.e., when no error estimate is used, one can only guess in advance the mesh quality notion and it becomes then important to translate this evaluation in terms of a quality function about the mesh elements. For a problem of an isotropic nature, considering the equilateral triangle as an optimal element is a natural choice, while for a quad the square is a priori optimal. When no sizing information (about the element size or the desired edge length) is provided, the only reasonable behavior is to take advantage of the boundary (assumed to be appropriately sized) and to complete regular elements whose size follows at best the size of the boundary items. A finer discretization in a boundary region results in smaller elements, while a coarser discretization may lead to larger elements. Then, these sizing features will be used to find a reasonable size value for the elements located at some distance from the domain boundaries. In particular, a progressive gradation of the sizes is often a desirable feature when considering the elements in a given neighborhood.

When a size map is provided via a control space, we will see that a unit length in the metric associated with this control space is the targeted value. This leads us to the notion of a unit mesh as the good mesh we like to create.

Definition 1.18 A unit mesh is a mesh with unit edge lengths.

Comparing the edge length with the metric specification as defined in the control space provides the expected actual values. For the moment, it is enough to say that, in an isotropic context, a unit length in a metric map specifying a value h actually gives a length h in the usual sense and that a similar notion is also valid in an anisotropic context.

Remark 1.11 A triangle with unit edge lengths is necessarily a good quality triangle. This is not the case of a tetrahedron that may have an inconsistent volume (quite small) while its edges are (about) unit edges.

After this remark about simplicial elements and due to the other various geometries that can be envisaged, edge length control must be coupled with other criteria to make sure that a good quality is obtained.

Notice that the unit value is related to the underlying metric map and that the latter could be the combination of a map of a geometric nature (due to the domain geometry) with one or several other maps of a physical nature depending on the behavior of the problem under investigation.

To conclude this brief overview of the main concepts of mesh generation, let us indicate that numerous detailled discussions will be provided in the following chapters to clarify all the notions introduced in this chapter as well as practical advices to deal with them in computational applications.

---

1Let us briefly recall the definition of a d-simplex: we consider d + 1 points aj = (aij =1 ∈ d, 1 ≤ j ≤ d + 1, , not all in the same hyper-plane, i.e., such that the matrix of order d + 1 :

is invertible. D-simplex K whose vertices are the aj is the convex hull of these points aj. Every point x in dwith Cartesian coordinates xi is fully specified by the data of d + 1 scalar values λj = λj (x) that are solutions of the linear system:

whose matrix is A. The λj ( x ) are the barycentric coordinates of point x with respect to the points aj.

2A (−l)-face is the empty set, a 0-face is a vertex, a 1-face is an edge, a k-face is in fact a k-simplex with k < d, d being the spatial dimension.A (−l)-face is the empty set, a 0-face is a vertex, a 1-face is an edge, a k-face is in fact a k-simplex with k < d, d being the spatial dimension.

3A set of points is said to be in general position if there is no configuration of more than three points that are co-circular (more than four co-spherical points) such that the corresponding open disk (ball) is empty.

4A 2-face swap (flip) consists of replacing the diagonal of the convex quadrilateral made up of two adjacent triangles by the alternate configuration, see Chapter 18 for the precise definition.

5Given a set of geometric objects, an arrangement is a covering-up of the space by means of the regions (cells) formed by the given objects and their (potential) intersections.

6Whereas a constrained Delaunay triangulation in two dimensions is a triangulation which satisfies the empty sphere criterion, where a open ball can contain a vertex in the case where the latter is not seen, due to a constrained edge, by all the vertices of the considered element. In other words, a constrained entity exists which separates the above vertices and the others.

7It should be noted that people with a finite element background use the term triangulation and use the term mesh synonymously.

8Note that some authors use the term grid to refer to any kind of mesh whatever its connectivity.

9Given a vertex numbering (index) based on an implicit definition results in implicit definitions for both the edges and the faces, thus avoiding an explicit definition of these entities at the element level, which would be not unique and memory consuming.

10Similarly to the triangle, an alternative definition also suits well where face [i] is opposite vertex (i). Actually, the latter convention leads to greater simplicity.

11In this way, it will be possible to add some fields that are not yet defined at the time such or such capability must be made available.

12For a boundary element, a projection will be needed to obtain the desired location.

Chapter 2

Basic Structures and Algorithms

The aim of this chapter is to introduce a variety of data structures and to show how they can be used profitably in a mesh generation context. To this end, some basic as well as more sophisticated data structures are recalled together with some algorithms of greater or lesser complexity. The discussion is then developed by means of various application examples related to situations extensively used in a meshing context.

While people having a computer science background may be familiar with these basic notions, we would nevertheless like to address this topic here in order to allow people less directly concerned to gain some knowledge of this (vast) topic, notably applied mathematicians or numericians. Moreover, in the mesh generation context, specific applications and uses of the classical data structures lead to specific situations and merit some comments.

The literature about data structures and algorithms is quite abundant. Among the usual references, the textbooks would include [Ahoet al. 1983], [Wirth-1986], [Sedgewick-1988], [Cormenet al 1990], [Samet-1990] and [Gonnetet al 1991] as well as the unchallengeable [Knuth-1998a], not to mention many others that can also be consulted.

The complexity of an algorithm, both in terms of the number of operations and of the memory resource allocated, is analyzed from a theoretical point of view. However, we will show that specific theoretical results obtained inad-hoc academic situations must be slightly nuanced. Indeed, numerous assumptions like the points must be in general position in a triangulation problem or all operations involved in a given numerical process have the same cost or again are exactly computed are unlikely to be what we meet in real life. Nevertheless, despite these remarks, theoretical results allow for a good understanding of some difficulties and will help us to find appropriate solutions.

Therefore, after having described the theoretical point of view in the first sections of this chapter, we give some indications and remarks about the most commonly encountered difficulties in realistic applications. After that, we turn to some application examples to illustrate how to benefit from the theoretical material.

The first section introduces the general problem from an academic point of view. The second section presents the most commonly used elementary data structures (array, list, stack, etc.). The third section deals with complexity problems for a given algorithm. Section four analyzes the sorting and searching techniques and introduces three main paradigms used in many methods. Data structures in one and two dimensions are discussed in sections five and six, while topological data structures are mentioned in section seven. Sections eight and nine deal respectively with the notions of robustness and optimality of an implementation. The last section proposes several practical application examples.

2.1. Why use data structures

As an introduction, we look at a naive algorithm that can be used to construct a triangulation. Let consider a set of points each of which contained (to simplify even more) in a single initial triangle. The algorithm consists of finding, for each and any point P, the triangle K enclosing it and then to subdivide K into three new triangles by connecting P to the three vertices of K:

• For all P ∈

– Find triangle K containing point P,

– Subdivide this triangle into three.

• End

While very simple, this algorithm raises several questions. Among these, we simply mention the need to define the concept of a triangulation and how to represent it, for instance, by using adjacency relationships between the triangles. Another question is related to the quick identification of the triangle containing the point P to be inserted. Should we examine all triangles of the current triangulation or take into account the fact that any triangle is obtained by subdividing its parent?

This simple example gives some indications on how to proceed and what to know to implement such an algorithm (simple as it may be). This is not indeed restricted to defining the operations required to code this algorithm, but also to finding the data structure(s) adapted to the problem, in such a way as to define a Program. According to [Wirth-1986], we have the following paradigm:

Algorithm + Data Structures = Program.

There is obviously a close link between an algorithm and the data structures it uses. Usually, the more complex a data structure, the simpler the algorithm will be, although the simplicity of the algorithm is generally altered during the data structure update. For triangulation (meshing) algorithms in particular, a rich data structure allows useful data to be stored and retrieved thus simplifying the task of the algorithm, but on the other hand, any modification of the mesh induces a set of modifications and updates of the data structure. On the other hand, a simple data structure is efficient to update but forces the algorithm to perform explicitly a set of operations to recreate some data each time needed. As an example, a data structure that keeps only the adjacency relationships between the triangles of a mesh provides instantaneously the neighbors of a given triangle but requires an algorithm to identify all the triangles sharing a common vertex. However, if this information is stored, it can be retrieved immediately, provided it has been updated as the mesh evolves.

Notice that a (mesh) data structure contains integer values (numbers, indices, etc.) as well as real values (triplets of coordinates of the vertices in three dimensions, for instance). In Section 2.2, we will describe data structures allowing this kind of information to be stored.

A data structure being fixed, we discuss the behavior of the algorithm. Therefore, we recall, in Section 2.3, several fundamental notions about the complexity, allowing to analyze the efficiency of standard algorithms and basic data structures. In Section 2.4, we describe a series of algorithms, based on one of the computer science paradigms,namely, Divide and Conquer. We give some examples of methods for searching and sorting, some of which we describe, such as the insertion sorting technique, the quicksort and bucket sorting as well as binary searching (dichotomy) or interpolation methods. Section 2.5 discusses the manipulation of entities of dimension one (integers). To this end, we look at:

• general data structures allowing to store, retrieve or analyze sets of objects,

• structures allowing a selective access to some entities already stored. The access can be performed according to several criteria of selection. We find here, for instance, the approaches where the smallest item (in some sense), the first or the last recorded, the neighbor(s) of a given item, etc. is sought. Here we will find the data structures like stack (LIFO), queue (FIFO), priority queues, array with sorting and binary searching trees.

• data structures like dictionaries that can provide answers to questions like does this item exist? and allow items to be inserted or suppressed. We will find here BST and hash coding techniques.

In Section 2.6, we discuss how to use data structures in two and three dimensions for fast storing and retrieving of items such as points, segments (edges) or polygons. Section 2.7 is devoted to the computer implementation of topological data. After this overview of basic data structures and algorithms, we discuss robustness problems inherent to any implementation of a mathematical expression in a computer. The degree of the problems and the notion of predicate are then analyzed as well as the cost in terms of the number of operations and of memory requirements (Sections 2.8 and Sections 2.9). To conclude, we mention some applications where the previously described material can be used, in the specific context of the development of mesh generation and modification algorithms (Sections 2.10).

2.2 Elementary structures

In this section, we describe tables (arrays), pointers, lists, stacks and queues. These structures are briefly introduced below on deliberately simple examples.

Table or array

The table or the array is most certainly