Archive for the ‘Meshing’ Category

Algebraic Mesher in Python

Monday, July 14th, 2008

This year I have tried to make it a habit to track how much time I spend on numerical programming. I’d like to put in 10 hours or so a week and this lets me know if I am doing that. Anyway, I went back to my records from February when I was developing the Algebraic Mesher in C. It took me 20.9 hours to complete, with about half of that time spent in a frustrating effort to track down numerous memory problems. The promise of scripting languages is they allow you to develop programs much more efficiently, at the price of slower run time performance. So I decided to create the mesher in Python and see how long it took me.

Method

Four months had elapsed since I finished the C version and I couldn’t recall how I had programmed it. All I remembered is it was a miserable experience. Thus it seemed that I could create a Python version from scratch without undue influence from the previous effort. To get started, I reviewed the equations linking the ID numbers for vertices, faces and cells from earlier, but nothing more. I didn’t look at the C code or my follow up blog entry on the experience.

So I opened up the IDLE editor and looked at a perfectly blank page. I wasn’t sure where to start. Then I thought: I need a cell object. So I made a very simple one:

class cellclass:
    def __init__(self,CID):
        self.CID = CID

All it does is store the ID number of the cell. I made a list of cells and stored some numbers. So far so good. I then made a vertex object, and created a list of vertices:

verts=[]
for i in range(NumRows+1):
    for j in range(NumCols+1):
        VID = (NumCols+1)*i + j
        x = j*delx; y = i*dely
        temvert = vertclass(VID,x,y)
        verts.append(temvert)

In this way I added features incrementally, in baby steps, without really any planning. Proceeding along in this fashion I had the mesh generator finished in 1.5 hours. I couldn’t believe how fast it went. The code wasn’t especially pretty – nearly everything is lumped in one routine – but it seemed to run fine.

Testing

I wanted to make sure I tested the code as thoroughly as I had tested the C version, so I opened up the C code at this point, reviewed the tests there, and replicated them, one by one in Python. With the C version, every test revealed new bugs which had to be laboriously tracked down and fixed. With the Python version, the code passed the first four tests the first time. The fifth test found a bug where I was writing multiple copies of vertices to the face objects. I fixed that in a couple of minutes.

In a little more than 2 hours I had a working mesher. And it was fun to write. The code just seemed to flow out of my finger tips (and I am not especially skilled at Python). Yes, this code probably should be rewritten into a more modular style. And perhaps it could be made more ‘Pythonic’. So maybe another hour or two are needed to really polish it up. But there is something extremely satisfying about getting a rough version working with so little effort.

One major disappointment with the C version of the program was the memory problems. I had spent some time developing a generic list structure in C, so I could make my programming in that language less painful. The list and FOR_EACH macro were supposed to make C more like Python. The thought was good, but the experience was not.

Stats

The table shows the time and lines of code needed for each version of the code. From Prechelt’s study, we expect Python to be 3.8x faster and 3.6x shorter than C/C++. Roughly the speed up in development time mirrors that in lines of code. In my case, the speed up was nearly 9x, while the needed lines of code was about 2x less. This may be due to the inordinate amount debugging time the C code required. In any case the general trend holds: Python is a lot faster to develop in. 

.            Time            --Lines of code --
Lang      (hrs)   Mesher   Test Code  Total
-------------------------------------------------------------------------------------
 C         20.9    299      163         462
Python      2.4    124       89         213
-------------------------------------------------------------------------------------
Ratio       8.7x   2.4x     1.8        2.2x

Run Times

The big downside of fast development time is slow execution time. Python does all kinds of nice things for you (like memory management), but that comes at a price. To measure this, I had both the C and Python meshers create a million cell mesh and timed them. Here are the results:

.            Time      RAM
Lang      (sec)     (GB)
------------------------------------------------------------------------------------
 C            6.9   0.366
Python      344     1.33
------------------------------------------------------------------------------------
Ratio       50x     3.6x

This test was run on an Intel Core2 Duo (T7300) at 2 GHz, with 2GB of RAM. I think much of the time, for both programs, went into memory allocation. This was corroborated when I tried compiling the C version with various optimizations and got no change in the timing. There may some ways to allocate memory in bigger chunks, and thus speed up the code.

There are tools for speeding up the execution of Python programs. Tools like Psyco and Pyrex can offer dramatic improvements, at least under some conditions. That might be a good experiment for next time.

File: PyAlgMesher.py

Posted in Meshing, Programming, Python | No Comments »

Algebraic Mesher

Thursday, March 13th, 2008

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.

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

Cell Data Structure

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.

Posted in Data Structures, Meshing, Programming | No Comments »

Meshing by the Numbers

Saturday, March 8th, 2008

We saw earlier that the half edge data structure is a very elegant way to create a mesh. However, as I pointed out recently, this approach does have a couple of drawbacks for general CFD applications. Are there other elegant approaches worth considering? One interesting possibility is to use the algebraic relationships between vertices, faces and cells to generate the connections between them. If the mesh is rectangular, and you number the components systematically, you can work out these relationships.

Cell and Vertex Numbering

Cells & Vertices

Lets say we number our cells from left to right, starting with the bottom row, as in the figure. Any given cell can be expressed in terms of its column number and row number:

CID = N_cols*i + j

The Cell ID (CID) is found by multiplying the number of columns in the grid (N_cols) by the current row number i, and adding the current column number j. We are following the C convention of starting our numbering from zero. We can number the vertices in the same way, starting in the lower left corner. Without too much trouble, we can see that if you know the ID of the cell, you can work out the ID of the south west vertex of that cell:

VID_SW = CID + i

where VID_SW is the ID number of the southwest vertex. The equation works out this way because any row of N cells has N+1 vertices.

The ID for the southeast vertex is just one more bigger than the southwest:

VID_SE = CID + i + 1

To get the index for the northern vertices, we need to add a row of vertices to the index of the corresponding southern vertex:

VID_NW = CID + i + (N_Cols + 1)
VID_NE = CID + i + 1 + (N_Cols + 1)

Face Numbering

Faces

The numbering scheme for the faces is a bit trickier. You can just name them sequentially starting from the south west corner, but I think its a bit easier if you create two lists: a list of horizontal faces (i.e., the north and south faces) and a list of vertical faces (i.e., the east and west faces). The figure shows the idea, with the horizontal faces numbered in black and the verticals in blue. With this approach, the southern faces have the same numbering as the cells, thus

FID_S = CID

The northern face ID numbers can be found by taking the southern face ID and adding the number of columns in the grid:

FID_N = CID + N_cols

For the east and west faces, each row has one more face than there are cells. This is the same as the vertex numbering, thus

FID_W = CID + i
FID_E = CID + i + 1

In practice, we don’t really want to have two lists of faces. It is more convenient if we can loop through a single list. This can be achieved by appending the east-west list onto the north-south list. Once that is done, the indices for the east and west faces will need to be increased by the number of horizontal faces:

FID_W = N_HFaces + CID + i
FID_E = N_HFaces + CID + i + 1

where

N_HFaces = N_Cells + N_cols

.

Implementation

I have written a C program to implement these ideas, but it was a torture. I made heavy use of my generic list structure and in the process hit a bunch snags. C is notorious for giving you the power to shoot yourself in the foot, and I dare say that I don’t have many toes left. I made a ton of memory mistakes and it took me more than 10 hours of debugging to find and correct them all. In the process, I made a number of improvements to the generic list routines (including writing a bunch of documentation), so hopefully this won’t happen again. I’ll discuss the implementation more in the next post.

Posted in Meshing, Programming | No Comments »

Rethinking the Half Edge Mesh

Saturday, March 1st, 2008

In earlier discussions (here and here), I extolled the benefits of the half edge mesh. It is an elegant way of compactly storing 2d mesh information. Moreover, it allows you to easily do adjacency tests (e.g., what are all the edges touching this vertex?). However, upon further reflection, two weaknesses have occurred to me:

  1. It is lacking for the finite volume method.
  2. Its not extendable to 3d

I don’t know that these are show stopper considerations, but they do give one pause. Lets look at them in more detail.

Finite Volume

The key tenet of the finite volume method is the balance of fluxes leaving and entering each cell face. This is what makes this approach intrinsically conservative, unlike the finite element method.

At some point in the execution of a finite volume code, we loop over the faces in the domain and calculate the fluxes crossing those faces. Because each pair of cells share one face, we know that the material leaving cell A must flow into cell B. But the half edge data structure doesn’t have faces; it has half faces. There is no single entity that connects both cells. So how do we achieve this conservation? The simplest way is probably to create a full face data structure that each half face points to. We can then calculate the fluxes on this face as normal. While this adds additional complexity to the scheme and reduces its elegance, it should be perfectly doable.

3d

What really makes me doubt the viability of the half edge approach is in generalizing it to 3d. In 2d, the vertices, edges and cells are connected together in a simple fashion. For instance, the edges around a cell are connected, one to the next, in a chain. This makes it easy to create a linked list of edges to go around the cell.

Face linkage in 2d, 3d

In 3d, each cell is surrounded by six faces. There are many ways to connect these faces together in a linked list. The figure shows one of them. That is workable, but how do you get an adjacency test out of it? For instance, lets say you want to know all the faces that surround a vertex. I can’t figure out any clever way of looping through the half faces to achieve this.

Next Steps

I am working on another mesh generator that uses algebraic relationships between cells, faces and vertices to establish the connections between these entities. It has a certain elegance to it as well. Unfortunately, in implementing it, I hit all kinds of memory and pointer problems that took awhile to work out.

Posted in Data Structures, Meshing | No Comments »