basicnumerics.com Blog

Lets say we have a rectangular mesh, perhaps created with our half mesh tool, and we want to imprint upon it some polygon shape. For instance, we may want to model the temperature or fluid flow in such a polygon and need a mesh in that shape. By imprinting the polygon onto the mesh, we can create a simple stair step mesh that can be used for modeling.

What is the criteria for determining if a cell is inside the polygon or not? Clearly any cell that is fully inside the polygon would be considered “in”. What about one where one of the vertices is outside? Or two are outside? For the figure above, the imprinting follows a simple rule: if the cell center is in the polygon, then the cell is inside.

Shooting Rays

To create our stair step mesh, we need to cycle through the cells in the rectangular mesh and check each cell center to see if its inside. The classic way to determine if a point lies inside a polygon is to shoot a ray from this point and see how many times it crosses the polygon. If its an even number, the point is outside. An odd number of crossings means its inside. The Wikipedia has a nice article on this here.

Hitting a vertex – A rare Event?

There is a subtlety with this scheme. Your ray may strike one of the vertices of the polygon. In that case, you will likely get double counting. It will appear that your ray has crossed both polygon segments N and N+1. As a result, you may conclude exactly the wrong result. My initial reaction to this situation was: it could happen, but will be extremely rare. Wrong! Both meshes and polygons are created by humans. And humans tend to like round numbers. It is much more likely you will have a vertex point at (1,1) than at (1.0001, 0.9999). Similarly for the cell positions. In my own testing, I hit this problem right away.

The solution is found in how we think about a segment. Rather than thinking of a segment as two points and the line that connects them, we can think of it as just one point and the line. The figure shows the idea.

This can be implemented with a judicious use of the > and >= operators.

Implementations

A quick search of the internet reveals many implementations of this odd-even crossing approach. Bob Stein’s article compares a few implementations, looking at lines of code, number of if statements and so on. The routines he looks at range from 44 to 91 lines of code and include from 2 to 19 if statements, and 0 to 2 divides. When I first implemented this, some years ago, my code ran 70 lines and included 7 if statements and 3 divides. In my defense, my code also checked how close a point was to a segment and included the point if it was within some tolerance. Still, it seems a bit bulky compared to some of the implementations Stein looked at. He was justifiably proud of his implementation, which is only 44 long with no divides.

A bit more looking on the web reveals the incredibly tiny 12 line routine of W. Randolph Franklin. It contains only one if and one divide. I tried it out for a couple of cases and it does seem to capture the essense of the method without a single wasted character. Here is the Franklin’s routine:

int pnpoly(int npol, float *xp, float *yp, float x, float y)
{
  int i, j, c = 0;
  for (i = 0, j = npol-1; i < npol; j = i++) {
    if ((((yp[i]<=y) && (y<yp[j])) ||
         ((yp[j]<=y) && (y<yp[i]))) &&
        (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))

      c = !c;
  }
  return c;
}

(Note this code is Copyright © 1970-2003, Wm. Randolph Franklin )

Stair Step Mesh

With a routine like this, we can easily create a stairstep mesh. You can envision a program something like this:

int main()
{
    MESH_TYPE *pMesh;
    CELL_TYPE *pCell;
    int i;

    //Make a triangle
    float xp[10] = {0., 1.0, 1.0};
    float yp[10] = {0., 0., 1.0};

    //Create mesh
    pMesh = createMesh(0., 0., 1.0, 1.0, 0.1, 0.1);

    //Determine which cells in mesh are 'live'
    FOR_EACH(i, pCell, *(pMesh->cells))
    {
        pCell->live = pnpoly(3, xp, yp, pCell->x, pCell->y,)
    }
    return 0;
}

For this to work, we need to add a few more fields to the cell structure defined previously. Specifically, we need the cell center coordinates and a flag stating whether the cell is live or dead.

Some additional work would be required to actually use this stair step mesh in a CFD code. For instance, it is not efficient to keep dead cells in our cell list. Every time we loop through the cells, we will need an if statement:

. Also, we have not identified the boundary faces or made any effort to transfer the boundary conditions from the polygon segments to the edge faces of the stair step mesh. These are tasks for another day.

Having done much of my programming in Basic, I have been sweetly oblivious to the importance of case in naming variables. In Basic, NumNodes is the same as numnodes, numNodes and NUMNODES. As a result, you can type any combination that suits you and all is well. Your eye doesn’t even notice the differences after awhile. That all changes when you program in C. I am translating some meshing code from PowerBasic to C for my next blog entry and the gcc compiler is giving me tons of errors like:

HalfMesher.c: In function `createMesh':
HalfMesher.c:242: error: `OldNF' undeclared (first use in this function)
HalfMesher.c:254: error: structure has no member named `Faces'

And that is after going through the code line by line looking for such things!

CompoundNames

This got me thinking about naming conventions for variable names and functions. There is an amazing amount written on the topic and feverous religious battles have been fought. Still, a few patterns emerged. In the C world, there is a tendency to use all lower case for simple variable and function names. For compound names, like CellVolume, there are several options:

• Run together: cellvolume
• Underscores: cell_volume
• Lower mixed case: cellVolume
• Upper mixed case: CellVolume

Mixed case is one of several names for capitalizing the first letter of each word. (Others are camel case, medial capitals, InterCaps.) In lower mixed case, this is applied only to the words after the first. In the list above, the first two combining methods seem the most common for C programmers, followed by lower and finally upper mixed case.

I have to say that I don’t care for the run together approach because it can be ambiguous for some word combinations. For instance consider the variable cellslip. It could be interpreted as cell_slip or cells_lip. As you read the code, you need to stop and decipher this variable name. That may only take a moment, but it is an extra burden on the reader. The underscore method addresses this issue nicely, but adds to the bulk of variables, and can clutter the code. Mixed case gives you the benefits of word separation, without the extra bulk. Because C programmers have tended to prefer lower case, lower mixed case seems like the best approach for compound names.

Other Issues

It turns out there are many more issues that come up with name conventions. For instance, terseness. The C examples in Kernighan and Ritchie are filled with one and two letter variable names. Their code is compact, but it can be hard to read. At the other extreme, there is the natural language approach of Kari Laitinen. Here is a snippet of code taken from his site:

cout <<  "\n Please, type in two integers separated "
        <<  "with spaces:  " ;

   cin  >>  first_integer_from_keyboard
        >>  second_integer_from_keyboard ;

   sum_of_two_integers  =  first_integer_from_keyboard  +
                           second_integer_from_keyboard  ;

It is very readable, but that is partly because it is also very simple. If you are solving complex equations, with lots of variables in them, your code will be highly cluttered if you follow this approach. Which is easier to read:

y = a*x*x + b*x + c;

or

Answer_to_problem_7 = First_Coefficient*x_variable*x_variable
                    + Second_Coefficient*x_variable
                    + Third_Coefficient;

The example above introduces the importance of idioms. Even very terse names can convey meaning if they follow a common idiom. Any mathematically aware persion will recognize the above equation as that of a second order polynomial. For this case single letter names are fine. Similarly, using i and j as index names for loops is a common idiom in programming, and so is fine too.

What is the optimum length for variable names? Steve McConnell discusses this in Code Complete. He references a 1990 study by Gorla, Benander and Benander that suggests names in the range of 10 to 16 characters where the easiest to debug.

Another important point is that care is required with global variables. Local variables can be understood in the context of the function they reside in. Global variables are generally without context. As you read code, they pop up randomly, and you need to understand them. Thus, their names must be descriptive enough to standalone, within reason.

Naming Convention

Here are the rules I have come up with for myself when programming in C. Following these should both make my life easier, and hopefully make my code easier to read.

• Single word variable names are lower case: face
• Compound variable names use lower mixed case: theLongName
• Functions follow same convention as variables
• Macros are upper case: MAX
• Typedefs are upper case and end in _TYPE for complex types: FACE_TYPE
• Names for lists or collections typically end in ’s’: cells
• Global variables start with g and are more detailed: gDiscretizationScheme
• If it helps readability, pointers can be prefixed with p: pFace
• Abbreviations are to be minimized.
• Idioms are used to shorten names where it makes sense