basicnumerics.com Blog

In looking for an alternative to the half edge meshing scheme, I hit upon the idea of using algebraic relations between cells, faces and vertices in order to build a mesh. Those relationships were described in the previous post. Here I want to describe how I implemented this in C.

Data Structures

Like my Half Edge mesher, the MESH_TYPE data structure consists of lists of all the vertices, faces and cells in the mesh.

struct mesh_type {       //contains all mesh info
    ALIST cells;   //cell list
    ALIST faces;   //face list
    ALIST verts;   //vert list
};

ALIST is the generic array I had defined earlier. It consists of a small data structure that includes a pointer to an array on the heap containing the actual data. In this way, the array can be resized and bounds checking can be done. The figure shows the mesh data structure.

Up to this point, the structures are the same as in the Half Mesh program. Where things change is in the data structures for the cells, faces and vertices. For instance a face structure looks like

struct face_type {
    ALIST verts;  //list of ptrs to face verts
    ALIST cells;  //list of ptrs to faces 1 or 2 cells
};

Each face in the list of faces contains lists of pointers to elements in these three master lists. Trying to draw all the relationships can get confusing. This pictures gives an idea for the faces of a cell

Each cell in the master list contains a list of four face pointers. These pointers point to the four faces in the master list where this cell’s faces reside. Imagine how confusing that picture would be if I were to try and draw all the connections between cells, vertices and faces?

Accessing the Data

Suppose you want to average the x and y coordinates all the vertices of a cell, to make sure it matches the cell center. (This is one of the tests I used to check the meshing code.) This can be done easily using the new improved FOR_EACH macro from the generic list include file:

double totX, totY,dx,dy;
VERT_TYPE **ppVert;
FOR_EACH(pCell, pMesh->cells)
{
    totX = 0; totY = 0;
    FOR_EACH(ppVert, pCell->verts)
    {
        pVert = *ppVert;
        totX += pVert->x[0];
        totY += pVert->x[1];
    }
    //Find diff betw cell center and vert ave
    dx = pCell->x[0] - totX/4;
    dy = pCell->x[1] - totY/4;

The outer loop is over the cells in the master list of cells for the mesh. The way the FOR_EACH works, is it puts a pointer to the current cell into pCell on each pass through the loop. This allows us to dereference the fields of the cell. It is the inner loop, over the cell vertices, where confusion can set in. The vertex list in a cell structure holds pointers to the actual vertices residing in the master vertex list. Since FOR_EACH returns a pointer to whatever is in the list, it returns a pointer to a pointer. That is why I defined the variable ppVert. In order to gain access to the x and y position of the vertex, we dereference ppVert, to get a pointer to the vertex. From their we can get the position. These dereferencing subtleties were not at all clear to me as I wrote this. In addition, I hadn’t always allocated memory for my lists. The result was the memory debugging nightmare I mentioned previously.

Assessment

The entire motivation for this project was to generate a meshing program that had some of the elegance of the half edge mesh, while being better suited for CFD applications. Was it successful? It is hard for me to say at this point. In terms of code length, the two implementations are nearly identical. Thus one is not a lot clumsier than the other. One thing I find troubling about the algebraic meshing approach is you can’t look at the code and tell if its right. When you see a line like

pVert = getItem(verts,CID+i+1);  //SE vert

how do you know that pVert is really the south east vertex of cell CID? You need to sketch out the grid and work out the relationships. Also, I am still reeling a bit from all the memory problems I had. The half edge mesh routine came together a lot more easily.

All in all, I would say the jury is still out on this one. But now that I have created it, I hope to use it to do some work with Immersed Boundaries.

Download

If you want to play with the code to generate meshes, you can download it here. Put the files ArrayType.inc and AlgMesher.c in the same folder and compile and run:

gcc -Wall AlgMesher.c
./a.exe

When main runs, it will call a test function that tests five different aspects of a simple 100×100 mesh. It writes out the results to a text file. In this way, I tried to make the code relatively bug free. Enjoy.