basicnumerics.com Blog

In a previous post, I discussed the elegance of the Half Edge data structure for storing mesh information. This time, we will implement such a scheme in C. The goal is simple: create a program that meshes a 2d rectangle with a structured grid, while storing the information in the Half Edge data structure.

If we know the length and width of the rectangle, and the cell size, we can compute the number of rows and columns in our mesh. We can then loop over these rows and columns, creating cells, faces and vertices and linking them together. As we create these mesh elements, we will store them in lists, using the generalized arrays we developed earlier.

We will build the mesh row by row. The first step is to create the south most row of vertices. As we construct a row of cells, we will be building them on top of the vertices of the previous row. So for the very first row cells, we must first create a row of verts. Once this is done, we enter a pair of nested for loops where we loop over the rows and columns of the mesh. Here is the basic idea

for (i=0; i<numRows; i++)
    {
        Prepare the row
        for (j=0; j<numCols; j++)
        {
           Set corner vertices of cell
           Create new cell
           Create faces of cell
           Make remaining connections
           Hook up the pairs of WF and NF faces
           Save a few items to help with next cell or next row
         }
        wrap up this row
    }

Once the first row of verts is created, we begin creating cells. For the first cell of a row, we must create a north west vertex. For all other cells on that row, the north east vertex of the previous cell becomes the north west vertex for the new cell. We then need only create one vertex per cell, specially the north east one. This is shown in Step 2 of the figure. With the four vertices of the cell defined, we can then create the four faces of the cell (step 3), then create the cell and link everything together. Its then on to the next cell. In this step by step fashion we build up the mesh.

A subtlety

There is a subtlety I didn’t mention in the previous post. There are 8 half faces that touch a vertex. How do you know which one to store in the vertex? For an interior vertex it doesn’t really matter, as long as its consistent with the scheme for cycling through the faces around the vertex. Specifically, if you store a face pointing toward the vertex, then you cycle through the faces with the command

face = face->pair->nextF;

If you pick a face that points toward the vertex, then you use

face = face->nextF->pair;

For vertices at the edge of the domain, greater care is needed. Not every face has a pair near the edge, so a command like

face = face->pair->nextF;

can lead to errors. Also, we need to pick the edge that allows us to go fully around the vertex. In the picture below, this can only be satisfied for vertex V by picking face F.

In general, we need to pick the most counter clockwise edge, relative to the vertex. (That is because our circular linked list of faces is organized in a counter clockwise direction.)

More on Lists

In writing this simple mesher, I made use of the generalized arrays developed earlier. In the process of using these, I came across various things that didn’t work that well, so I changed them. For instance, the listAppend function originally returned an index to the data array where the new data was written. It turned out to be more convenient to return a pointer to that data.

I also realized that the creation and destruction functions were overly restrictive. They assume the user has already created the list object and only needs the data array allocated. There are times however, when the user would like both allocated and deallocated for him. So there are now two sets of creation and destruction functions:

void createListData(ALIST *list, long maxEl, size_t elSize)
void destroyListData(ALIST *list)

void createList(ALIST **list, long maxEl, size_t elSize)
void destroyList(ALIST *list)

This provides greater flexibility, but also adds complexity. Speaking of complexity, the astute reader will note in the source code that I did not use these flexible arrays for the lower and oldNF arrays. Instead I used fixed arrays, dimensioned to 1000 elements. A flexible array would be perfect here, since it will auto-expand to however many rows are needed. But the extra hassle didn’t seem worth it for some something so simple.

If I had used these flexible arrays, setting and getting values would have looked like:

LL = *(VERT_TYPE *)getItem(lower, j);
setItem(lower, j) = UL;

instead of

LL = lower[j];
lower[j] = UL;

That is a serious loss of readability. When getting the item from the list, you have to remember that getItem returns a pointer to the data element, so you must dereference it. To be proper, you need to cast the pointer first. That makes for some pretty darn ugly code. I could use a macro to clean it up. Perhaps something like:

#define GET_ITEM(pList, index, datatype) *(datatype *)getItem(pList, index);

I will have to experiment with it.

Wrap Up

If you want to play with the code to generate a half edge mesh, you can download it here. Put the files ArrayType.inc nd HalfMesher.c in the same folder and compile and run:

gcc -Wall HalfMesher.c
./a.exe

When main runs, it will call a test function that tests different aspects of a simple 10×10 mesh. It writes out the results to a text file. This is how I made sure the code was relatively bug free.