The stencil module consists of two parts:
In Directions.h, all possible directions are defined as an enumeration. Each direction has several properties. For example: the inverse direction, the offset in x,y,z direction, etc. These properties are stored as constant arrays which can be indexed using the direction enum value:
The only context in which direction enums should be used as integers, are above property arrays. Do not use them as indices to the fourth coordinate of fields!
An alternative implementation would have been to store directions as classes which hold these properties (cx, inverseDir, ...) as members. This was not done since this prevents some compile-time optimizations of the compiler. For best possible performance, these compile-time optimizations are, however, essential! Above code snippets can all be evaluated at compile time, since all arrays are declared constant.
With all possible directions and their properties defined, a stencil class now only has to pick a subset of these directions. The generation of the stencil classes are automated by the python script generate.py. The script uses the file Stencil.in.h as template and substitutes strings according to the stencil definition. Which stencils are generated is specified at the beginning of the script:
Defines a D-dimensional stencil info with Q directions. A stencil is defined by picking a subset of all possible directions. These directions are listed in the dir array. So when iterating a stencil, one iterates the dir array which contains the directions. To get the properties for these directions, the global arrays of Directions.h can be used:
Since all these arrays are constant, the lookup can be resolved at compile time and no performance overhead should be generated.
For more convenient iteration, use the provided iterator. The following code is equivalent to the code example above:
In LBM simulation, one typically stores a PDF for every direction contained in a stencil. To write code that is independent of the stencil, the stencil type is usually passed as template parameter. The size of the f coordinate (= fourth dimension of the field) equals the number of directions contained in the stencil and can be obtained by Stencil::Size (or simply Stencil::Q).
The reverse mapping is also possible: When one has a stencil::Direction enum at hand and needs to determine the index of this direction, the Stencil::idx array has to be used. Do not use the enum directly as index: field( x, y, z, W ) is not correct in the general case! Typically you want field( x, y, z, Stencil::idx[W] )!
There are special cases (for example D3Q19 or D3Q27) where both variants are coincidentally the same.