Some Notes on Delaunay Triangulations

This is just a sketch of some tools that are available for computing Delaunay triangulations of 2D pointsets. In some cases, there are similar tools that can be used for pointsets in 3D and higher dimensions. We are most concerned about

In more advanced cases, we may "constrain" the problem, by requiring that a particular set of edges be included as part of the triangulation, or by specifying that the region we work in is not the convex hull of the points, but some smaller region. The boundaries of this smaller region will be line segments that begin and end on elements of the pointset. The resulting region might also include holes.

Definition of a Triangulation

You would surely know a triangulation if you saw one. You have a set of points P, and you connect them with lots of non-intersecting lines so the only shapes you have are triangles.

However, your computer won't be able to see it that way. And if we need to answer a technical question, we may need to have a more precise definition.

We begin with a set of points P in the two dimensional plane. We may reasonably assume that these points are distinct. In some cases, we may want to assume that no three points are collinear, that is, that they do not lie along a line.

Given the set P, if we connect some or all of them by some straight line segments which never intersect each other, then we have created a planar subdivision with vertices P. It divides the plane up into regions, with the edges forming boundaries.

For most planar subdivisions, it is easy to add another edge. However, no matter what the form of our set P, there is always at least one maximal planar subdivision, to which no more edges can be added. (Actually, there are usually many different maximal planar subdivisions.)

It turns out that if a subdivision is maximal, then all the (finite) regions must be triangles; so we now define a triangulation of a point set P to be a maximal planar subdivision whose vertices are the point set P.

Facts about any Triangulation

(Number of triangles): The number of triangles and edges in any triangulation can be determined from two other values, namely NP, the number of points in P, and NB, the number of those points that lie on the convex hull of the pointset. If we rule out the case where all the points in P are collinear, then any triangulation has 2*NP-2-NB triangles and 3*NP-3-NB edges. There are NB edges on the boundary, and hence 3*NP-3-2*NB interior edges.

Note, for example, in the following simple graphic, that there are NP=9 points, and that there are NB=7 of them on the boundary. This immediately implies that there are 2*9-2-7=2 interior nodes, and 3*9-3-7=17 edges, of which 7 are boundary edges, and 10 interior.

      |\  |  /|\
      | \ | / | \
      |  \|/  |  \
       \  |  /|  /
        \ | / | /
         \|/  |/

For constrained problems, this result may still apply. For instance, if we have only constrained the problem by specifying a "tighter" boundary, but have not introduced any "holes" in the region, then the formula still applies, with NB now counting the number of points that lie along the user-specified boundary line.

For example, let us suppose that we have constrained the previous problem by requiring that the boundary look as follows:

      |        \
      |         \
      |          \
      X---X   X  X
          |      /
          |     /
          |    /
There are still NP=9 points, and but now NB=8 of them on the boundary. There are 2*9-2-8=1 interior nodes, and 3*9-3-8=16 edges, of which 8 are boundary edges, and 8 interior.
      |\  |  /|\
      | \ | / | \
      |  \|/  |  \
          |  /|  /
          | / | /
          |/  |/

A Delaunay Triangulation

Given a set of points P, there are usually a vast number of triangulations possible. For each triangulation T of these points, we can assign a value Theta(T), the minimum angle measure. Theta(T) is defined to be the smallest angle in all the triangles that make up the triangulation. We may now define the Delaunay triangulation of a set of points P to be any triangulation T for which Theta(T) achieves its maximum value.

Facts and Properties about the Delaunay Triangulation

(Maximum Theta): we repeat the defining property of a Delaunay triangulation: that it maximizes the value of Theta(T), the minimum angle in the triangles. This fact alone is enough to recommend it to people who create meshes for solving partial differential equations. Assuming that you must use a given set of points, then the Delaunay triangulation is the best in this sense.

(Not unique): The Delaunay triangulation need not be unique. For the simplest example, take P to be the four points that form a square. There are two possible triangulations, they have the same value of Theta(T), so they are both Delaunay triangulations. For general point sets, however, it is usually the case that the Delaunay trianulation is unique.

(Empty circumcircle): every (non-degenerate) triangle has a unique circumcircle, that is, a triangle which passes through the three vertices. Every triangle of a Delaunay triangulation satisfies the empty circumcircle criterion. That is, for each Delaunay triangle, the circumcircle does not contain any other points of the pointset. (In some special ambiguous cases, other points might lie just on the circumcircle, but they will never lie within the circumcircle.)

(Acyclicity): Given v, an arbitrary viewpoint in the plane, we say that triangle t1 lies in front of triangle t2 if there is a ray starting at v that passes first through t1 and then through t2. In a Delaunay triangulation, this induces a partial ordering of the triangles. (Surprisingly, this is not the case with a general triangulation. It is possible in that case for, say, t1 to be in front of t2 which is in front of t3 which is in front of t1!) This induced ordering is important, for instance, in computer graphics, where the "Painter's Algorithm" is used to draw a collection of objects correctly. You simply draw the "far away" ones first. If a near and far object overlap in the visual field, then the near one will show up, because it will have been drawn after the far away one.

This property is also useful when trying to determine which triangle contains an arbitrary point in the plane. You can start in any triangle whatsoever, and consider the straight line from the center of that triangle towards the point in question. For a Delaunay triangulation, you can always reach the point by stepping from your current triangle to the neighbor that shares the edge through which the line passes.

(Voronoi Dual): The dual graph of the Delaunay triangulation of the pointset P is the Voronoi diagram of P, and vice versa.

(Number of triangles): since the Delaunay triangulation is a triangulation, we have the same result, that there are 2*NP-2-NB triangles and 3*NP-3-NB edges in the triangulation.

DIAMOND: A sample dataset

To keep things simple, we'll only be concerned with 2D Delaunay triangulations, and we'll concentrate on a single small dataset of points whose coordinates are contained in the unit square. Here is a typewriter plot of the data:

        2  .  .  .  .  .  .  .  .  .  9
        .  .  .  .  .  .  .  .  .  .  .
        .  .  .  .  .  .  .  .  .  .  .
        .  .  .  .  .  .  .  .  .  .  .
        .  .  .  4  .  .  .  .  .  .  .
        .  .  3  .  5  .  7  .  .  .  .
        .  .  .  .  .  .  .  .  .  .  .
        .  .  .  .  .  .  6  .  .  .  .
        .  .  .  .  .  .  .  .  .  .  .
        .  .  .  .  .  .  .  .  .  .  .
        1  .  .  .  .  .  .  .  .  .  8

I've stored a copy of the coordinates of the points in the file ../data/xy/diamond_02_00009.xy. Here's what it looks like:

        #  diamond_02_00009.xy
        #  This set of data was set up simply to study the computation
        #  and reporting of the vertices of Voronoi cells.
        #  In particular, the data was selected so that the point
        #    ( 0.4, 0.5 )
        #  has a Voronoi cell with vertices:
        #    ( 0.3, 0.3 ), ( 0.5, 0.4 ), ( 0.5, 0.7 ), ( 0.3, 0.6 ).
          0.0  0.0
          0.0  1.0
          0.2  0.5
          0.3  0.6
          0.4  0.5
          0.6  0.3
          0.6  0.5
          1.0  0.0
          1.0  1.0

One useful thing about this dataset is that I can tell you EXACTLY the shape of the Voronoi polygon around the point labeled 5. It consists of the four points (0.3,0.2), (0.5,0.4), (0.5,0.7), and (0.3,0.5), which form a sort of diamond around the point. This fact should help later on, when we are trying to understand the output of some of the programs.

What Data Do We Want?

In the simplest cases, the input to the problem is simply the number NP and coordinates of the set of points P. For a constrained problem, it may be necessary to supply somewhow the number NC of edges that must be included and the pairs of nodes that form them, or the sequences of nodes that are connected to form one or more connected curves that are the boundary of the region.

For output, it is often enough for us to return triangle_num, the number of triangles, and a list of the triples of points that form the triangles of the Delaunay triangulation. When convenient, we may refer to this as the triangle array, and think of it as having 3 rows and triangle_num columns.

It is customary to list the vertices in a consistent orientation. In mathematical geometry, this ordering is usually such that the vertices are given in counterclockwise orientation. Even then, there are three possible orderings for the vertices, and it may be helpful to choose the ordering in which the vertex of lowest index is listed first. Finally, to produce a unique description of a given triangulation, we should also lexically order the list of triangles. It is also worth noting that FORTRAN indices begin at 1, while C/C++ indices start at 0. This affects the values of the triangle array, since they are actually a list of indices.

Under these conventions, the triangle array for the diamond data set is:

           1     3     2
           1     6     3
           1     8     6
           2     3     4
           2     4     9
           3     5     4
           3     6     5
           4     5     7
           4     7     9
           5     6     7
           6     8     7
           7     8     9

For more sophisticated uses, it may be helpful to have another list, which, for any side of any triangle, gives the index of the triangle that shares that side, or returns a value of -1, perhaps, if there is no such triangle. (A value of -1 indicates that that side of the triangle forms part of the boundary).

PLOT_POINTS - A Skimpy but Exact Picture

Let's use the program PLOT_POINTS. This is an interactive program. We'll keep the interaction short though:

        plot_points diamond_02_00009.xy
If we didn't include the input file name on the command line, we could have indicated it directly to the program by a command like
        read diamond_02_00009.xy
The delaunay command tells the program to compute the Delaunay triangulation and save it to an Encapsulated PostScript file. I converted this to a JPEG file using GSView.

MATHEMATICA'S DelaunayTriangulation Command

The Mathematica program includes a discrete mathematics package that can compute Delaunay triangulations. The following commands set up our data. (Note that, when working with Mathematica, we have to avoid names containing an underscore because Mathematica uses the underscore character for other purposes!)

We begin by clearing the names of the variables we intend to use, read in the appropriate package, define our pointset, and make a plot of the points.

        Clear[ points ]
        Needs[ "DiscreteMath`ComputationalGeometry`" ]
        gxy = {  {0.0,0.0},  {0.0,1.0},  {0.2,0.5}, {0.3,0.6}, 
                  {0.4,0.5},  {0.6,0.3},  {0.6,0.5}, {1.0,0.0}, 
                  {1.0,1.0} } 
        ListPlot [ gxy, PlotJoined->False ]

Now we can compute the convex hull of the pointset. If we save the output, which is a list of indices, then we can also draw a plot of the convex hull:

        c_list = ConvexHull [ gxy ]

Now we invoke the DelaunayTriangulation command, which returns a "vertex adjacency list". A vertex adjacency list is a a list of lists. Item 1 of the list is an ordered pair containing the index of the node 1, and a list of the nodes to which it is connected. Item 2 reports the same things for node 2, and so on.

        v_adj = DelaunayTriangulation [ gxy ]

It should be clear that, given a vertex adjacency list, and an understanding of how to extract information from it, we could draw the Delaunay triangulation line by line. Simply consider each item in the vertex adjacency list, which has the form

        { i, { neighbors of node i} }
and draw a line from the node to each of its neighbors.

Although we don't deserve this good luck, it turns out that we can easily get a plot of the Delaunay triangulation without having to understand all these details, because Mathematica has a command which does all the work for us:

        PlanarGraphPlot [ gxy ]

Mathematica also includes a command that, given a pointset and a vertex adjacency list that represents a triangulation, returns TRUE or FALSE, depending on whether the triangulation is a Delaunay triangulation:

        DelaunayTriangulationQ [ gxy, v_adj ]

Thus, depending on our needs, we know the Mathematica commands necessary to make the various kinds of plots, or get the list of Delaunay triangles for our given pointset.


The MATLAB program is a very handy interactive program. It includes a number of computational geometry routines, and in particular, a delaunay command that computes the Delaunay triangulation. (It also includes the DELAUNAYN command for higher dimensions, and a very interesting search command using the Delaunay triangulation.)

With an eye to future work, let's see if we can get the data in the DIAMOND file into MATLAB. One way to do this is to use the MATLAB TABLE_READ routine. If we've copied the file table_read.m, we can tell MATLAB to get our file as follows:

        p = table_read ( 'diamond_02_00009.xy' );
Then we can compute the Delauany triangulation with the command
        triangle = delaunay ( p(1,:), p(2,:) );

This command returns a triangle array. For convenience, we also display the vertex sorted, and then triangle sorted versions:

           7  9  8  => 7  8  9 => 1  3  2
           6  8  1     1  6  8    1  3  6
           6  7  8     6  7  8    1  6  8
           3  2  1     1  3  2    2  4  3
           3  6  1     1  3  6    2  4  9
           4  3  2     2  4  3    3  5  4
           4  9  2     2  4  9    3  6  5
           4  7  9     4  7  9    4  7  5
           5  4  7     4  7  5    4  7  9
           5  4  3     3  5  4    5  6  7
           5  6  7     5  6  7    6  7  8
           5  3  6     3  6  5    7  8  9
Note that in this example, the vertices are not consistently listed in clockwise or counterclockwise order; the vertices in a triangle are not ordered; the triangles are not ordered; vertex indexing begins at 1.

Of course, our first desire is to see the triangulation, but to do that, we need to massage the data a little, by pasting on a fourth column that is a copy of the first one:

        [ ntri, ncol ] = size ( tri );
        tri(:,4) = tri(:,1)
Now we can plot the points, and then add lines for each triangle:
        scatter ( x, y )
        for i = 1 : ntri
          line ( x(tri(i,:)), y(tri(i,:)) )

The result can be exported directly from MATLAB to a JPEG file, and looks something like this:


For serious applications (including datasets in 3D or higher dimension), MATLAB includes a more general function called DELAUNAYN.

The DELAUNAYN command returns the Delaunay tessellation of a pointset in M dimensions. The form of the command is

        t = delaunayn ( x )
where x is an N (how many) by M (spatial dimension) array of points. The output t is an NUMT by M+1 array, where each row is a set of M+1 indices of points in x which form a Delaunay simplex. (Remember that in 2D, we get groups of 3 points that make a triangle, in 3D we get 4 points that form a tetrahedron.)

The CONVEXHULLN command returns the indices of a pointset that form its convex hull:

        k = convexhulln ( x )
where x is an N (how many) by M (spatial dimension) array of points. The output k is a list of NUMB indices of the points that lie on the convex hull. (The documentation is ambiguous, but I believe this is correct.) You can also get the volume of the convex hull using the form of the command
        [ k, v ] = convexhulln ( x )

Given a tessellation t of points x, the TSEARCHN command returns that simplex in t which encloses the point xi, or NaN if the point lies outside the convex hull.

        index = tsearchn ( x, t, xi )

Given a tessellation t of points x, the DSEARCHN command returns the index of the closest point in x to the point xi. It is also possible to have xi be a vector of points.

        index = dsearchn ( x, t, xi )
Ignoring the cost of setting up the tessellation, this approach is much faster than the brute force method of simply computing the distance of xi to each point x and taking the closest. The fixed cost of setting up the tessellation can be paid off quickly if there are many points xi, and if the spatial dimension is not too large.

Fortune's SWEEP2

Steve Fortune wrote a C program called SWEEP2 which can produce the Voronoi diagram or Delaunay triangulation of a set of points. To compute the Delaunay triangulation, the point set must be in a file, say points.txt, with each line containing the coordinates of a single point. Issuing the command

sweep2 -t < points.txt > points_del.txt
will create an output file containing the Delaunay triangulation. Each line of the output file is a triple of indices that indicate the points included in one triangle. The indexing starts at 0.

For the diamond dataset, the sweep2 program produced an output file containing the triangle data. We have listed the original data, then incremented the indices by 1, then reordered the vertices, then sorted the triangles.

        0 5 7  =>  1 6 8 => 1 6 8 => 1 2 3
        2 5 0      3 6 1    1 3 6    1 3 6
        2 4 5      3 5 6    3 5 6    1 6 8
        4 6 5      5 7 6    5 7 6    2 4 3
        3 4 2      4 5 3    3 4 5    2 9 4
        6 7 5      7 8 6    6 7 8    3 4 5
        3 6 4      4 7 5    4 7 5    3 5 6
        1 3 2      2 4 3    2 4 3    4 7 5
        8 7 6      9 8 7    7 9 8    4 9 7
        0 1 2      1 2 3    1 2 3    5 7 6
        3 8 6      4 9 7    4 9 7    6 7 8
        1 8 3      2 9 4    2 9 4    7 9 8
We remark that the original data is in clockwise order; the vertices are not sorted; the triangles are not sorted; vertex indexing begins at 0.

I wrote a program, SWEEP2_DELAUNAY_EPS which reads the text files of the point locations and the Delaunay triangles computed by SWEEP2 and returns an EPS file containing an image of the triangulation. For the diamond dataset, the result is:

Barry Joe's GEOMPACK

Barry Joe wrote several versions of a FORTRAN package called GEOMPACK, which computes the Voronoi diagram or Delaunay triangulation using fairly simple data types. Most versions of the package contain many more routines than you need. The fundamental routine is called either DTRIS or DTRIS2 or RTRIS or RTRIS2; the D or R was used to distinguish double and single precision versions, and the 2 was used after Barry Joe developed 3D versions of the code as well.

I have a version of this program, which I converted to FORTRAN90, in the directory ../f_src/geompack/geompack.html

I have a second version of a subset of this program, which I converted to C++, in the directory ../cpp_src/geompack/geompack.html While it is much less extensive than the FORTRAN version, it includes the DTRIS2 routine, and so can do a Delaunay triangulation.

The DTRIS2 code computes the Delaunay triangulation directly from the coordinates of a point set.

Start with g_num points with coordinates in g_xy. To generate the Delaunay triangulation, we call

        call dtris2 ( g_num, g_xy, triangle_num, triangle, tnbr )

Warning! It is important to note that DTRIS2, as written, sorts the point coordinates, and on exit, returns the sorted point coordinates. This can be very disruptive if you are not aware of it! I find this fact so unpleasant that, in many cases, I have modified DTRIS2 so that it works with an index sort vector instead, producing the same output as before, except that the point coordinates are not modified. Before I made this change, I repeatedly encountered severe problems, such as making a plot using the old node coordinates, but with triangulation data that assumed the node coordinates had been shuffled by the sorting. The result is a meaningless spaghetti plot! And this hardens my belief that a subroutine's arguments should be INPUT or OUTPUT but NOT both!

Interpreting the output of DTRIS2

The number of Delaunay triangles is returned as triangle_num; for our sample problem, this value is triangle_num=12.

DTRIS2 returns a triangle array, which we can vertex sort and then triangle sort for convenience:

         2  1  3  =>  1  3  2 => 1  3  2
         3  1  6      1  6  3    1  6  3
         2  3  4      2  3  4    1  8  6
         4  3  5      3  5  4    2  3  4
         7  4  5      4  5  7    2  4  9
         5  3  6      3  6  5    3  5  4
         7  5  6      5  6  7    3  6  5
         9  4  7      4  7  9    4  5  7
         6  1  8      1  8  6    4  7  9
         7  6  8      6  8  7    5  6  7
         7  8  9      7  8  9    6  8  7
         2  4  9      2  4  9    7  8  9
We remark that this DTRIS2 triangle array has vertices listed in counterclockwise order; vertices are not ordered; triangles are not ordered; vertex indexing begins at 1.

The array tnbr indicates, for each Delaunay triangle the index of the neighboring triangle along each side. However, some triangles don't have a neighbor on a particular side. In such cases, the corresponding entry is a negative value:

         1    -1   2   3
         2     1   9   6
         3     1   4  12
         4     3   6   5
         5     8   4   7
         6     4   2   7
         7     5   6  10
         8    12   5  11
         9     2  -2  10
        10     7   9  11
        11    10  -3   8
        12     3   8  -4
(The actual negative values returned by DTRIS2 are different, and somewhat arbitrary. I've made them consecutive for simplicity.)

Versions of DTRIS2

As mentioned above, DTRIS2 and its underlying routines are available in the FORTRAN90 code GEOMPACK. That version still does the sorting of the point coordinates. I have also extracted these routines for use other packages, such as PLOT_POINTS program and . DELAUNAY_2D, where it is used to determine the triangulation of a set of points.

I also put a copy of the routines in the FORTRAN90 GEOMETRY library, along with many routines for printing out the arrays associated with the triangulation. The test code for GEOMETRY includes a setup of the DIAMOND dataset.

I have just recently worked out a C++ version of DTRIS2, made available in the C++ GEOMETRY library. This does not yet include the extra routines for printing out the triangulation arrays The test code includes a setup of the DIAMOND dataset.


Robert Renka has published a number of algorithms in ACM TOMS. In particular, TRIPACK (ACM TOMS Algorithm #751) computes the Delaunay triangulation of a set of points in the plane, in order to perform interpolation of scattered data in the plane. This is an interesting application of Delauanay triangulation in its own right. I have not looked at the details of how the Delaunay triangulation is computed, stored and used in this package.

Another algorithm published by Renka is STRIPACK, (ACM TOMS Algorithm #772). This is essentially the same as TRIPACK, except that the points lie on a sphere.

Shewchuk's TRIANGLE

TRIANGLE is Jonathan Shewchuk's C program for producing meshes, Delaunay triangulations and Voronoi diagrams. You can look at my local copy of TRIANGLE. The input file for the program, which specifies the location of the points, varies slightly from the pure list of coordinates we have used elsewhere. (In particular, it requires an explicit declaration of the spatial dimension and number of points.)

Given an appropriate input file, we invoke TRIANGLE with the "-p" option:

        triangle -p diamond_02_00009.poly
it quickly creates four files:

The file diamond_02_00009.1.ele corresponds to our data structure triangle. I indicate the vertex and triangle reorderings:

        12  3  0
           1       1  3  2 => 1  3  2 => 1  3  2
           2       4  3  5    3  5  4    1  6  3
           3       2  4  9    2  4  9    1  8  6
           4       4  2  3    2  3  4    2  3  4
           5       3  1  6    1  6  3    2  4  9
           6       6  8  7    6  8  7    3  5  4
           7       8  6  1    1  8  6    3  6  5
           8       7  9  4    4  7  9    4  5  7
           9       7  5  6    5  6  7    4  7  9
          10       9  7  8    7  8  9    5  6  7
          11       7  4  5    4  5  7    6  8  7
          12       6  5  3    3  6  5    7  8  9
        # Generated by triangle -v diamond_02_00009.node
We remark that the TRIANGLE vertices are counterclockwise oriented; vertices are not ordered; triangles are not ordered; vertex indexing begins at 1.

TRIANGLE includes another program called SHOWME which can be used to make an X-window display of the triangulation, or to save an image as a Encapsulated PostScript file.

To see the Delaunay triangulation we just created, it is only necessary to type

        showme diamond_02_00009.1.ele
I save the image as an EPS file, and converted that to JPEG using GSView.

Note that TRIANGLE has a built-in capability to handle constrained Delaunay triangulation, and holes.


The QVORONOI program is part of the QHULL package, which has a home page at

QVORONOI is an interactive program that reads an input file, computes the Delaunay triangulation, determines the Voronoi diagram, and either prints a summary, or outputs the information to a file, based on command line switches.

The format for the pointset coordinate file is fairly simple. Lines beginning with a nonnumeric character are treated as comments. The first noncomment line must contain the spatial dimension; the second line gives the number of points, and each subsequent line lists the coordinates of a single point. Thus, our standard diamond example file would be slightly modified from "xy" format to a new format we could call "qxy", which might look like this:

        #  diamond_02_00009.qxy
          0.0  0.0
          0.0  1.0
          0.2  0.5
          0.3  0.6
          0.4  0.5
          0.6  0.3
          0.6  0.5
          1.0  0.0
          1.0  1.0

The command line switch "i" requests that a triangle file be written containing the indices of the nodes that make up the Delaunay triangles. Thus, the command

        qvoronoi diamond_02_00009.qxy i TO diamond_del.txt
results in the creation of a file containing the Delaunay information. The contents of this file follow, along with versions of the data that have been incremented by 1, vertex sorted, and triangle sorted:
        7 5 0 => 8 6 1 => 1 8 6 => 1 3 2
        2 1 0    3 2 1    1 3 2    1 6 3
        5 2 0    6 3 1    1 6 3    1 8 6
        2 5 4    3 6 5    3 6 5    2 3 4
        8 6 7    9 7 8    7 8 9    2 4 9
        6 5 7    7 6 8    6 8 7    3 5 4
        5 6 4    6 7 5    5 6 7    3 6 5
        2 3 1    3 4 2    2 3 4    4 5 7
        3 2 4    4 3 5    3 5 4    4 7 9
        3 8 1    4 9 2    2 4 9    5 6 7
        6 3 4    7 4 5    4 5 7    6 8 7
        3 6 8    4 7 9    4 7 9    7 8 9
and we remark that in the QVORONOI output, vertices are listed in counterclockwise order; vertices are not ordered; triangles are not ordered; and vertex indexing begins at 0.

If you are working with qvoronoi and want to make plots, one option would be simply to use the sweep2_delaunay_eps program. In that case, you will need to delete the header lines from the node and triangle input files that report the dimension, the number of nodes and the number of triangles.

Considerations for Creating a Finite Element Mesh

When creating a mesh for a finite element problem, probably the most fundamental feature is the modeling of the boundary (and possibly internal walls) of the region. If we are free to choose where the points go, then we probably need to ensure that a good number of points end up exactly on the boundary, with special care taken to place points at corners and along curved boundaries. If the boundary does not coincide with the convex hull of the points, then we need to constrain our Delaunay calculation so that it does not attempt to generate extraneous exterior triangles, or worse yet, triangles that are partly in and partly outside the region. We also need to be careful about interior walls, and about holes, which are more troublesome to handle.

The Area of a Triangle in 2D

In the unconstrained case, and in the case where we specify the boundary, but that boundary is formed of straight line segments, and we have been careful to place a node at every vertex of this boundary, then our triangulation procedure will decompose the region into triangles. (This would not be true, for instance, if we did not have a point at every vertex, or our region had curved boundaries!)

In a finite element calculation, these triangles will play the role of elements. It is often necessary to compute the area of such an element, in order to approximate integrals of functions over the triangle. In more interesting cases, the elements may actually be more general polygons.

Whether the element is a triangle or polygon, the area of a polygon P with N vertices is easy to compute. We assume we are given the coordinates of the vertices and that the polygon is not degenerate (in particular, that the polygon is not folded over on itself.)

The simplest approach is to triangulate the polygon and add up the areas; if we aren't particular, we can choose vertex N as the base, and then consider the N-2 triangles formed by the base vertex plus consecutive pairs of vertices from (1,2) through (N-2,N-1):

        Area ( P ) = sum ( T in P ) area ( T )
where the area of a triangle is easily computed:
        Area ( T ) = 1/2 * ( x1 * ( y2 - y3 ) 
                           + x2 * ( y3 - y1 ) 
                           + x3 * ( y1 - y2 ) )
(In the formula for the triangle, the indices 1, 2, and 3 simply indicate the three vertices of the triangle, and don't refer to the original numbering of the vertices of the polygon.)

A second formula for polygonal area is simply

        Area ( P ) = 1/2 * sum ( 1 <= i <= N ) x(i) * ( y(i+1) - y(i-1) )
where we "wrap around" the indices when they exceed the legal range.

In MATLAB, if the coordinates of a polygon are stored in the column vectors p_x and p_y, then the area of the polygon can be found with the command

        area = polyarea ( p_x, p_y )

In some cases, instead of having the vertices of the Voronoi polygon around generator G, we might have the list of generators that are Delaunay neighbors. This is actually a fairly common situation. It is still possible to compute the area of the Voronoi polygon directly from this information, using a formula given in Spatial Tessellations [Okabe].

The GEOMETRY Library

The formulas for the area and centroid of a polygon in 2D, the volume and centroid of a polyhedron in 3D, the circumcenter of a triangle in 2D or tetrahedron in 3D, or the area, centroid and vertices of a Voronoi polygon in 2D given the Delaunay neighbors, are all implemented in the routines:

These routines are available in the GEOMETRY library, which can be had in a FORTRAN90 version, a C++ version, and a MATLAB version,


  1. Franz Aurenhammer,
    Voronoi diagrams - a study of a fundamental geometric data structure,
    ACM Computing Surveys,
    Volume 23, Number 3, pages 345-405, September 1991,
  2. C B Barber, D P Dobkin, H T Huhdanpaa,
    The Quickhull algorithm for convex hulls,
    ACM Transactions on Mathematical Software,
    December 1996.
  3. Marc deBerg, Marc van Kreveld, Mark Overmars, Otfried Schwarzkopf,
    Computational Geometry - Algorithms and Applications,
    Second Edition,
    Springer, 2000.
  4. Herbert Edelsbrunner,
    Geometry and Topology for Mesh Generation,
    Cambridge, 2001,
  5. Steve Fortune,
    A Sweepline Algorithm for Voronoi Diagrams,
    Algorithmica, Volume 2, pages 153-174, 1987.
  6. Barry Joe,
    GEOMPACK - a software package for the generation of meshes using geometric algorithms,
    Advances in Engineering Software,
    Volume 13, pages 325-331, 1991.
  7. Atsuyuki Okabe, Barry Boots, Kokichi Sugihara, Sung Nok Chiu,
    Spatial Tessellations: Concepts and Applications of Voronoi Diagrams,
    Second Edition,
    Wiley, 2000.
  8. Joseph ORourke,
    Computational Geometry,
    Cambridge, 1998,

You can return to the HTML page.

Last revised on 12 September 2005.