Archive for the ‘BFlow’ Category

Half Edge Mesh Implementation

Saturday, February 2nd, 2008

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.

MeshSteps

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.

Faces around a vertex

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.

Half Edge Mesh – Intro

Saturday, January 26th, 2008

When I created the mesh generator for BFlow, I used a straightforward but unimaginative set of data structures for storage. A mesh is made up of cells, faces, edges and vertices. Unstructured meshes have arbitrary relationships between these (unlike a structured mesh), and thus those relationships must be stored.

In a 2d, a mesh looks like this:

mesh

My BFlow data structures to hold the information for a mesh look like this:

TYPE CellType node(4) AS LONG '1=SW,2=SE,3=NE,4=NW (1st 4 are the corners) Face(4) AS LONG '1=East, 2=west, 3=North, 4=South nb(4) AS LONG 'Neighbor cells. 1=East, 2=west, 3=North, 4=South END TYPE TYPE FaceType node(2) AS LONG 'ptrs to two nodes of face cell(2) AS LONG 'ptrs to neighbor cells END TYPE TYPE VertType face(4) AS LONG 'pointer to touching faces, 1=East,2-West,3=North,4=South cell(4) AS LONG 'pointer to touching cells, 1=SW,2=SE,3=NE,4=NW END TYPE

(These are Basic declarations, but they would look very similar in C.) You can see that everything is stored. A cell can have four faces and four nodes. Pointers to them are stored in that cell structure. A vertex can be touched by four faces. Pointers to all those faces are stored for each vertex.

There is a very elegant alternative storage scheme called the Half Edge Mesh. Amazingly little information is stored, but it provides everything you need! The key idea is that rather than store edges (or faces in 2d), you store half edges (half faces). While an ordinary face is defined by its two end points, a half face is directed and defined by just a single point. The picture shows the idea.

mesh2

This probably looks like a step in the wrong direction, if we are trying to simplify things. Yet, this change is really helpful. The data structure for a half face looks like:

struct face_type { VERT_TYPE *vert; //vertex at the end of the half-face FACE_TYPE *pair; //Oppositely oriented adjacent half-face CELL_TYPE *cell; //cell that this face belongs to FACE_TYPE *nextF; //next half-face around the cell };

The vertex and cell structures are given below:

struct vert_type { double x; double y; FACE_TYPE *face; //Ptr to a touching face }; struct cell_type { FACE_TYPE *face; //south face of cell };

The cell structure seems especially spartan, at least compared with what I started with. How do you access the faces of the cell? Or its vertices? It turns out to be very easy:

FACE_TYPE* face = cell->face; do { // ...some face calcs... face = face->nextF; } while (face != cell->face);

The cell faces are effectively organized into a circular linked list. Thus the above code works equally well for a quad mesh as a polyhedral mesh.

Accessing all of the vertices around the cell is nearly as easy:

VERT_TYPE *vert; FACE_TYPE *face = cell->face; do { vert = face->vert; // ...some vert calcs... face = face->next; } while (face != cell->face);

Now lets try something trickier. Say you needed to loop over all of the cells surrounding a vertex? This could come up in a CFD model, for instance, when interpolating a field value from the cell centers to the vertices. Here is the code:

CELL_TYPE *cell; FACE_TYPE* face = vert->face; do { cell = face->cell; // ...cell calcs... face = face->nextF->pair; } while (face != vert->face);

This takes advantage of the fact that pairs of faces are oriented in opposite directions. In the figure below, face F1 is pointed upward. Its pair face, F2, is pointed downward. Each face stores a pointer to the vertex that it is pointing to. F1 stores a pointer to V1, while F2 stores a pointer to V2.

Mesh3

You can trace out logic by picking one of the faces pointing toward vertex V2, then following the steps. It is very elegant.

Compared to my original scheme, the Half Edge approach is both more efficient (in terms of memory use) and more flexible (since its not limited to quad meshes). In my next entry, I will explore this in more detail.

Why Basic?

Tuesday, December 18th, 2007

You may wonder why I chose the Basic computer language for writing the BFlow fluid library in. Basic isn’t exactly a mainstream language in computational modeling. The first version of Numerical Recipes supported it, both in text and CD-ROM, but later editions have not. A recent check of SourceForge found 11726 projects listed under Science and Engineering. Of those, only 30 were in Basic. Here is the breakdown:

Lang. Num Projs --------------------- Fortran 224 ( 1.9%) C 2477 (21%) C++ 3665 (31%) Python 1103 ( 9.4%) Basic 30 ( 0.3%) --------------------- Total 11726

I didn’t check all languages and some projects use more than one language. That is why the numbers don’t total up.

In addition to being little used in science, Basic suffers from a PR problem. The common view is that Basic is…well, basic. It is a toy language suitable for people unable or incapable of mastering “real” languages like C, C++ or Java. On top of all that, Basic is typically an interpreted language and thus is dog slow.

Given that dismal preamble, how did I end up choosing this language? It was historical. I generally use Excel for data analysis and quick calculations. Around 1995, Excel version 5 was released. This version introduced Visual Basic for Applications (VBA) to replace the Excel macro language. I found it incredibly easy to code in VBA. The code just flowed. And debugging was a joy. You just hold the cursor over any variable in scope and its value would pop up. In a very interactive way, you could write programs.

I began writing a flow solver in VBA. I knew it would be slow, but I wanted to experiment with some ideas I had, and thought it would be cool to have Excel solve some flow problems. Eventually I got the thing working, and sure enough, it was very slow. For many problems it was intolerably slow.

PowerBasic

In 2001, I looked around for alternatives. I didn’t want to rewrite everything, so I was hoping for a Basic compiler that was similar to VBA. I considered Microsoft Visual Basic, but Microsoft products tend to be big and bloated and overly complex. It turns out there are many other Basics out there. I came across the PowerBasic compiler. They emphasized the speed of their code and the easy porting from Visual Basic. So I sent in my $200 and gave it a try. In short order, I had my code ported over and was enjoying a 20x speedup.

I have been using PowerBasic since 2001 and find it to be a simple, robust compiler. The debugger isn’t as nice as the VBA debugger, but I can live with it. As implemented by PowerBasic, this language is quite powerful, having pretty much all the features found in C: arrays, structures, wide range of looping constructs, pointers, function pointers. It offers full access to the Windows API, for creating Windows applications. It comes with tools for creating user interfaces, doing TCP communication, talking to COM objects, inline assembler and many other things. All in all it is very capable.

I have found a few things that are unattractive for scientific computing. The compiler is optimized more for file size than speed. With care, you can make programs as fast as C or Fortran. With less care, you code can easily be 2x – 5x slower than an optimized C or Fortran program. Another problem: PowerBasic doesn’t allow forward declaration of user defined types (structs). This makes it hard to create several structures that reference each other.

In this blog, I will mostly use C for my experiments. It is widely used and understood and is fast. Over time, I may migrate BFlow to C if it makes sense to do so.