Reference documentation for deal.II version GIT relicensing837gaf73e64026 20240615 22:40:02+00:00

A module describing how deal.II consistently orients Triangulation objects. More...
Classes  
struct  GridTools::PeriodicFacePair< CellIterator > 
class  ReferenceCell 
class  internal::TriangulationImplementation::TriaObjectsOrientations 
Functions  
void  TriaAccessor< structdim, dim, spacedim >::set_combined_face_orientation (const unsigned int face, const unsigned char combined_orientation) const 
Querying or modifying topological information  
template<int dim>  
void  GridTools::consistently_order_cells (std::vector< CellData< dim > > &cells) 
Orientation of subobjects  
unsigned char  TriaAccessor< structdim, dim, spacedim >::combined_face_orientation (const unsigned int face) const 
bool  TriaAccessor< structdim, dim, spacedim >::face_orientation (const unsigned int face) const 
bool  TriaAccessor< structdim, dim, spacedim >::face_flip (const unsigned int face) const 
bool  TriaAccessor< structdim, dim, spacedim >::face_rotation (const unsigned int face) const 
bool  TriaAccessor< structdim, dim, spacedim >::line_orientation (const unsigned int line) const 
Orientation of subobjects  
static unsigned char  TriaAccessor< 0, dim, spacedim >::combined_face_orientation (const unsigned int face) 
Always return 0.  
static bool  TriaAccessor< 0, dim, spacedim >::face_orientation (const unsigned int face) 
Always return false.  
static bool  TriaAccessor< 0, dim, spacedim >::face_flip (const unsigned int face) 
Always return false.  
static bool  TriaAccessor< 0, dim, spacedim >::face_rotation (const unsigned int face) 
Always return false.  
static bool  TriaAccessor< 0, dim, spacedim >::line_orientation (const unsigned int line) 
Always return false.  
Orientation of subobjects  
static unsigned char  TriaAccessor< 0, 1, spacedim >::combined_face_orientation (const unsigned int face) 
Always return 0.  
static bool  TriaAccessor< 0, 1, spacedim >::face_orientation (const unsigned int face) 
Always return false.  
static bool  TriaAccessor< 0, 1, spacedim >::face_flip (const unsigned int face) 
Always return false.  
static bool  TriaAccessor< 0, 1, spacedim >::face_rotation (const unsigned int face) 
Always return false.  
static bool  TriaAccessor< 0, 1, spacedim >::line_orientation (const unsigned int line) 
Always return false.  
A module describing how deal.II consistently orients Triangulation objects.
A Triangulation is built not just of cells but also of lowerdimensional objects. In particular, a line in 2D is a face which may be shared by two cells, whereas in 3D a line may be shared by an arbitrary number of cells.
Lines (i.e., faces in 2D) are implicitly defined by the vertex numbering of the cell on which they are defined (i.e., there is no separate line or face object). This information is encoded in various ReferenceCell functions. For example, the first face of a triangle with nodes {0, 1, 2} is {0, 1}. Similarly, the first face of a triangle with nodes {1, 0, 3} is {1, 0}. By itself, this would be inconsistent, since the same line would be defined in twice: once by each of the two adjacent cells, and they might not agree whether the line should be {0, 1} or {1, 0}. To solve this problem deal.II also stores an unsigned char
on each cell for each of the cell's faces; this value encodes the orientation of the line as seen from this cell. In this particular case the orientation of the {1, 0} line from the {0, 1, 2} cell is ReferenceCell::default_combined_face_orientation() whereas the orientation of the {0, 1} line from the {1, 0, 3} cell is ReferenceCell::reversed_combined_line_orientation().
In deal.II, we express the orientation of an object with three booleans: orientation, rotate, and flip. The default values for these are true, false, and false. These values are typically encoded or decoded from or to a single unsigned char
by the internal::combined_face_orientation() and internal::split_face_orientation() functions.
For a quadrilateral, these values correspond to
true
is the default orientation and false
means vertices 1 and 2 are swapped.For a triangle, these values correspond to
true
is the default orientation and false
means vertices 1 and 2 are swapped.Here, 'clockwise' is relative to the vector defined by the cross product of two lines adjacent to the zeroth vertex in their standard orientation (which, e.g., points into the hexahedron for face 0 but out of the hexahedron for face 1).
For triangles, to enable indexing from the combined orientation, we do not consider fliprotate or fliporientrotate as those cases are equivalent, respectively, to the identity operation or the orientation = true
case as fliprotate is equal to the identity operation. As a consequence, there are only six valid orientations for triangles as faces of tetrahedra.
Purely quadrilateral meshes are a special case, since deal.II will (with the exception of faces which are neighbors across periodic boundaries) consistently orient purely quadrilateral meshes. Hence, in this case, the orientation of all lines will be ReferenceCell::default_combined_face_orientation(). See [1] for more information on this algorithm.
For example, in two dimensions, a quad consists of four lines which have a direction, which is by definition as follows:
* 2>3 *   * ^ ^ *   * 0>1 *
Now, two adjacent cells must have a vertex numbering such that the direction of the common side is the same. For example, the following two quads
* 345 *    * 012 *
may be characterized by the vertex numbers {0, 1, 3, 4} and {1, 2, 4, 5}, since the middle line would be {1, 4} when viewed from both cells. The numbering {0, 1, 3, 4} and {5, 4, 2, 1} would not be allowed, since the left quad would give the common line the direction {1, 4} whereas the right one would use {4, 1}, leading to an inconsistency.
As a sidenote, we remark that if one adopts the idea that having directions of faces is useful, then the orientation of the four faces of a cell as shown above is almost necessary. In particular, it is not possible to orient them such that they represent a (counter)clockwise sense, since then we couldn't already find a valid orientation of the following patch of three cells:
* o * / \ * o o *  \ /  *  o  *    * ooo *
(The reader is asked to try to find a conforming choice of line directions; it will soon be obvious that there can't exists such a thing, even if we allow that there might be cells with clockwise and counterclockwise orientation of the lines at the same time.)
One might argue that the definition of unique directions for faces and edges, and the definition of directions relative to the cells they bound, is a misfeature of deal.II. In fact, it makes reading in grids created by mesh generators rather difficult, as they usually don't follow these conventions when generating their output. On the other hand, there are good reasons to introduce such conventions, as they can make programming much simpler in many cases, leading to an increase in speed of some computations as one can avoid expensive checks in many places because the orientation of faces is known by assumption that it is guaranteed by the triangulation.
As a last question for this first section: is it guaranteed that such orientations of faces always exist for a given subdivision of a domain into cells? The linear complexity algorithm described below for 2d also proves that the answer is yes for 2d. For 3d, the answer is no (which also underlines that using such orientations might be an – unfortunately uncurable – misfeature of deal.II). A simple counterexample in 3d illustrates this: take a string of 3d cells and bend it together to a torus. Since opposing lines in a cell need to have the same direction, there is a simple ordering for them, for example all lines radially outward, tangentially clockwise, and axially upward. However, if before joining the two ends of the string of cells, the string is twisted by 180 degrees, then no such orientation is possible any more, as can easily be checked. In effect, some meshes could not be used in deal.II. In order to overcome this problem, the face_rotation
, face_flip
and line_orientation
flags have been introduced. With these, it is possible to treat all purely hexahedral meshes. However, in order to reduce the effect of possible bugs, it should still be tried to reorder a grid. Only if this procedure fails, the original connectivity information should be used.
As noted, reordering the vertex lists of cells such that the resulting grid is not a trivial problem. In particular, it is often not sufficient to only look at the neighborhood of a cell that cannot be added to a set of other cells without violating the requirements stated above. We will show two examples where this is obvious.
The first such example is the following, which we will call the `‘four cells at the end’' because of the four cells that close of the right end of a row of three vertical cells each (in the following picture we only show one such column of three cells at the left, but we will indicate what happens if we prolong this list):
* 91011 *   /  * 678  *     * 345  *   \  * 012 *
Assume that you had numbered the vertices in the cells at the left boundary in a way, that the following line directions are induced:
* 9>1011 * ^ ^ /  * 6>78  * ^ ^   * 3>45  * ^ ^ \  * 0>12 *
(This could for example be done by using the indices {0, 1, 3, 4}, {3, 4, 6, 7}, {6, 7, 9, 10} for the three cells). Now, you will not find a way of giving indices for the right cells, without introducing either ambiguity for one line or other, or without violating that within each cells, there must be one vertex from which both lines are directed away and the opposite one to which both adjacent lines point to.
The solution in this case is to renumber one of the three left cells, e.g. by reverting the sense of the line between vertices 7 and 10 by numbering the top left cell by {9, 6, 10, 7}:
* 9>1011 * v v /  * 6>78  * ^ ^   * 3>45  * ^ ^ \  * 0>12 *
The point here is the following: assume we wanted to prolong the grid to the left like this:
* oooooo *      /  * oooooo  *        * oooooo  *      \  * oooooo *
Then we run into the same problem as above if we order the cells at the left uniformly, thus forcing us to revert the ordering of one cell (the one which we could order as {9, 6, 7, 10} above). However, since opposite lines have to have the same direction, this in turn would force us to rotate the cell left of it, and then the one left to that, and so on until we reach the left end of the grid. This is therefore an example we have to track back right until the first column of three cells to find a consistent ordering, if we had initially ordered them uniformly.
As a second example, consider the following simple grid, where the order in which the cells are numbered is important:
* 32oo ... o76 *        *  0  N  N1  ...  2  1  *        * 01oo ... o45 *
We have here only indicated the numbers of the vertices that are relevant. Assume that the user had given the cells 0 and 1 by the vertex indices {0, 1, 3, 2} and {6, 7, 5, 4}. Then, if we follow this orientation, the grid after creating the lines for these two cells would look like this:
* 3>2oo ... o7<6 *        * ^ 0 ^ N  N1  ...  2 v 1 v *        * 0>1oo ... o4<5 *
Now, since opposite lines must point in the same direction, we can only add the cells 2 through N1 to cells 1 such that all vertical lines point down. Then, however, we cannot add cell N in any direction, as it would have two opposite lines that do not point in the same direction. We would have to rotate either cell 0 or 1 in order to be able to add all the other cells such that the requirements of deal.II triangulations are met.
These two examples demonstrate that if we have added a certain number of cells in some orientation of faces and can't add the next one without introducing faces that had already been added in another direction, then it might not be sufficient to only rotate cells in the neighborhood of the cell that we failed to add. It might be necessary to go back a long way and rotate cells that have been entered long ago.
From the examples above, it is obvious that if we encounter a cell that cannot be added to the cells which have already been entered, we can not usually point to a cell that is the culprit and that must be entered in a different orientation. Furthermore, even if we knew which cell, there might be large number of cells that would then cease to fit into the grid and which we would have to find a different orientation as well (in the second example above, if we rotated cell 1, then we would have to rotate the cells 1 through N1 as well).
A brute force approach to this problem is the following: if cell N can't be added, then try to rotate cell N1. If we can't rotate cell N1 any more, then try to rotate cell N2 and try to add cell N with all orientations of cell N1. And so on. Algorithmically, we can visualize this by a tree structure, where node N has as many children as there are possible orientations of node N+1 (in two space dimensions, there are four orientations in which each cell can be constructed from its four vertices; for example, if the vertex indices are {0 1 3 2}, then the four possibilities would be {0, 1, 3, 2}, {1, 3, 2, 0}, {3, 2, 0, 1}, and {2, 0, 1, 3}. When adding one cell after the other, we traverse this tree in a depthfirst (preorder) fashion. When we encounter that one path from the root (cell 0) to a leaf (the last cell) is not allowed (i.e. that the orientations of the cells which are encoded in the path through the tree does not lead to a valid triangulation), we have to track back and try another path through the tree.
In practice, of course, we do not follow each path to a final node and then find out whether a path leads to a valid triangulation, but rather use an inductive argument: if for all previously added cells the triangulation is a valid one, then we can find out whether a path through the tree can yield a valid triangulation by checking whether entering the present cell would introduce any faces that have a nonunique direction; if that is so, then we can stop following all paths below this point and track back immediately.
Nevertheless, it is already obvious that the tree has \(4^N\) leaves in two space dimensions, since each of the \(N\) cells can be added in four orientations. Most of these nodes can be discarded rapidly, since firstly the orientation of the first cell is irrelevant, and secondly if we add one cell that has a neighbor that has already been added, then there are already only two possible orientations left, so the total number of checks we have to make until we find a valid way is significantly smaller than \(4^N\). However, the algorithm is still exponential in time and linear in memory (we only have to store the information for the present path in form of a stack of orientations of cells that have already been added).
In fact, the two examples above show that the exponential estimate is not a pessimistic one: we indeed have to track back to one of the very first cells there to find a way to add all cells in a consistent fashion.
This discouraging situation is greatly improved by the fact that we have an alternative algorithm for 2d that is always linear in runtime (discovered and implemented by Michael Anderson of TICAM, University of Texas, in 2003), and that for 3d we can find an algorithm that in practice is usually only roughly linear in time and memory. We will describe these algorithms in the following. A full description and theoretical analysis is given in [1] .
The algorithm uses the fact that opposite faces of a cell need to have the same orientation. So you start with one arbitrary line, choose an orientation. Then the orientation of the opposite face is already fixed. Then go to the two cells across the two faces we have fixed: for them, one face is fixed, so we can also fix the opposite face. Go on with doing so. Eventually, we have done this for a string of cells. Then take one of the nonfixed faces of a cell which has already two fixed faces and do all this again.
In more detail, the algorithm is best illustrated using an example. We consider the mesh below:
* 91011 *   / *   /  *   /  * 678  *     *     *     * 345  *   \  *   \  *   \ * 012 *
First a cell is chosen ( (0,1,3,4) in this case). A single side of the cell is oriented arbitrarily (3>4). This choice of orientation is then propagated through the mesh, across sides and elements. (0>1), (6>7) and (9>10). The involves edgehopping and face hopping, giving a path through the mesh shown in dots.
* 9>1011 *  .  / *  .  /  *  .  /  * 6>78  *  .    *  .    *  .    * 3>45  *  .  \  *  X  \  *  .  \ * 0>12 *
This is then repeated for the other sides of the chosen element, orienting more sides of the mesh.
* 9>1011 *   / * v.....v.......V  *   /.  * 6>78 .  *    .  *    .  *    .  * 3>45 .  *   \.  * ^..X..^.......^  *   \ * 0>12 *
Once an element has been completely oriented it need not be considered further. These elements are filled with o's in the diagrams. We then move to the next element.
* 9>10>11 *  ooo  . / * v ooo v . V  *  ooo  . /  * 6>7>8  *   .   *   .   *   .   * 3>4>5  *  ooo  . \  * ^ ooo ^ X ^  *  ooo  . \ * 0>1>2 *
Repeating this gives
* 9>10>11 *  ooo  oooooo / * v ooo v ooooo V  *  ooo  oooo /  * 6>7>8  *     * ^.....^..X..^...^ *     * 3>4>5  *  ooo  oooo \  * ^ ooo ^ ooooo ^  *  ooo  oooooo \ * 0>1>2 *
and the final oriented mesh is
* 9>10>11 *   / * v v V  *   /  * 6>7>8  *     * ^ ^ ^ ^ *     * 3>4>5  *   \  * ^ ^ ^  *   \ * 0>1>2 *
It is obvious that this algorithm has linear runtime, since it only ever touches each face exactly once.
The algorithm just described in the twodimensional case is implemented for both 2d and (in generalized form) for 3d in this class. The 3d case uses sheets instead of strings of cells to work on. If a grid is orientable, then the algorithm is able to do its work in linear time; if it is not orientable, then it aborts in linear time as well.
Both algorithms are described in the paper "On orienting edges of unstructured two and threedimensional meshes", R. Agelek, M. Anderson, W. Bangerth, W. L. Barth, ACM Transactions on Mathematical Software, vol. 44, article 5, 2017. A preprint is available as arxiv 1512.02137.
Prior to the implementation of the algorithms described above (originally implemented by Michael Anderson in 2002, and reimplemented by Wolfgang Bangerth in 2016 based on the work in [1]), we used a branchandcut algorithm initially implemented in 2000 by Wolfgang Bangerth. Although it is no longer used, here is how it works, and why it doesn't always work for large meshes since its runtime can be exponential in bad cases.
The first observation is that although there are counterexamples, problems are usually local. For example, in the second example mentioned above, if we had numbered the cells in a way that neighboring cells have similar cell numbers, then the amount of backtracking needed is greatly reduced. Therefore, in the implementation of the algorithm, the first step is to renumber the cells in a CuthillMcKee fashion: start with the cell with the least number of neighbors and assign to it the cell number zero. Then find all neighbors of this cell and assign to them consecutive further numbers. Then find their neighbors that have not yet been numbered and assign to them numbers, and so on. Graphically, this represents finding zones of cells consecutively further away from the initial cells and number them in this frontmarching way. This already greatly improves locality of problems and consequently reduced the necessary amount of backtracking.
The second point is that we can use some methods to prune the tree, which usually lead to a valid orientation of all cells very quickly.
The first such method is based on the observation that if we fail to insert one cell with number N, then this may not be due to cell N1 unless N1 is a direct neighbor of N. The reason is obvious: the chosen orientation of cell M could only affect the possibilities to add cell N if either it were a direct neighbor or if there were a sequence of cells that were added after M and that connected cells M and N. Clearly, for M=N1, the latter cannot be the case. Conversely, if we fail to add cell N, then it is not necessary to track back to cell N1, but we can track back to the neighbor of N with the largest cell index and which has already been added.
Unfortunately, this method can fail to yield a valid path through the tree if not applied with care. Consider the following situation, initially extracted from a mesh of 950 cells generated automatically by the program BAMG (this program usually generates meshes that are quite badly balanced, often have many – sometimes 10 or more – neighbors of one vertex, and exposed several problems in the initial algorithm; note also that the example is in 2d where we now have the much better algorithm described above, but the same observations also apply to 3d):
* 131415 *  \   *  \ 4  5  *  \   *  121011 *     *    7  *     *  3  89 *     *    6  *     * 4567 *     *  2  1  0  *     * 0123 *
Note that there is a hole in the middle. Assume now that the user described the first cell 0 by the vertex numbers {2, 3, 6, 7}, and cell 5 by {15, 14, 11, 10}, and assume that cells 1, 2, 3, and 4 are numbered such that 5 can be added in initial rotation. All other cells are numbered in the usual way, i.e. starting at the bottom left and counting counterclockwise. Given this description of cells, the algorithm will start with cell zero and add one cell after the other, up until the sixth one. Then the situation will be the following:
* 13>14<15 *  \   *  > 4 v 5 v *  \   *  12>10<11 *     * ^   7  *     *  3 ^ 8>9 *     *   ^ 6 ^ *     * 4>5>6>7 *     * ^ 2 ^ 1 ^ 0 ^ *     * 0>1>2>3 *
Coming now to cell 7, we see that the two opposite lines at its top and bottom have different directions; we will therefore find no orientation of cell 7 in which it can be added without violation of the consistency of the triangulation. According to the rule stated above, we track back to the neighbor with greatest index, which is cell 6, but since its bottom line is to the right, its top line must be to the right as well, so we won't be able to find an orientation of cell 6 such that 7 will fit into the triangulation. Then, if we have finished all possible orientations of cell 6, we track back to the neighbor of 6 with the largest index and which has been added already. This would be cell 0. However, we know that the orientation of cell 0 can't be important, so we conclude that there is no possible way to orient all the lines of the given cells such that they satisfy the requirements of deal.II triangulations. We know that this can't be, so it results in an exception be thrown.
The bottom line of this example is that when we looked at all possible orientations of cell 6, we couldn't find one such that cell 7 could be added, and then decided to track back to cell 0. We did not even attempt to turn cell 5, after which it would be simple to add cell 7. Thus, the algorithm described above has to be modified: we are only allowed to track back to that neighbor that has already been added, with the largest cell index, if we fail to add a cell in any orientation. If we track back further because we have exhausted all possible orientations but could add the cell (i.e. we track back since another cell, further down the road couldn't be added, irrespective of the orientation of the cell which we are presently considering), then we are not allowed to track back to one of its neighbors, but have to track back only one cell index.
The second method to prune the tree is that usually we cannot add a new cell since the orientation of one of its neighbors that have already been added is wrong. Thus, if we may try to rotate one of the neighbors (of course making sure that rotating that neighbor does not violate the consistency of the triangulation) in order to allow the present cell to be added.
While the first method could be explained in terms of backtracking in the tree of orientations more than one step at once, turning a neighbor means jumping to a totally different place in the tree. For both methods, one can find arguments that they will never miss a path that is valid and only skip paths that are invalid anyway.
These two methods have proven extremely efficient. We have been able to read very large grids (several ten thousands of cells) without the need to track back much. In particular, the time to find an ordering of the cells was found to be mostly linear in the number of cells, and the time to reorder them is usually much smaller (for example by one order of magnitude) than the time needed to read the data from a file, and also to actually generate the triangulation from this data using the Triangulation::create_triangulation() function.
void GridTools::consistently_order_cells  (  std::vector< CellData< dim > > &  cells  ) 
Given a vector of CellData objects describing a mesh, reorder their vertices so that all lines are consistently oriented.
The expectations on orientation and a discussion of this function are available in the reordering module.
cells  The array of CellData objects that describe the mesh's topology. 
Definition at line 1615 of file grid_tools_topology.cc.
unsigned char TriaAccessor< structdim, dim, spacedim >::combined_face_orientation  (  const unsigned int  face  )  const 
Return an integer representation that uniquely encodes the orientation, flip, and rotation of a face
.

private 
Set the combined face orientation (i.e., the integer that uniquely encodes the orientation, flip, and rotation). This function is only implemented for objects which have faces, i.e., for structdim == dim.

static 
Always return 0.

static 
Always return 0.
bool TriaAccessor< structdim, dim, spacedim >::face_orientation  (  const unsigned int  face  )  const 
Return whether the face with index face
has its normal pointing in the standard direction (true
) or whether it is the opposite (false
). Which is the standard direction is documented with the GeometryInfo class. In 1d and 2d, this is always true
, but in 3d it may be different, see the respective discussion in the documentation of the GeometryInfo class.
This function is really only for internal use in the library unless you absolutely know what this is all about.
bool TriaAccessor< structdim, dim, spacedim >::face_flip  (  const unsigned int  face  )  const 
Return whether the face with index face
is rotated by 180 degrees (true
) or not (false
). In 1d and 2d, this is always false
, but in 3d it may be different, see the respective discussion in the documentation of the GeometryInfo class.
This function is really only for internal use in the library unless you absolutely know what this is all about.
bool TriaAccessor< structdim, dim, spacedim >::face_rotation  (  const unsigned int  face  )  const 
Return whether the face with index face
is rotated by 90 degrees (true
) or not (false
). In 1d and 2d, this is always false
, but in 3d it may be different, see the respective discussion in the documentation of the GeometryInfo class.
This function is really only for internal use in the library unless you absolutely know what this is all about.
bool TriaAccessor< structdim, dim, spacedim >::line_orientation  (  const unsigned int  line  )  const 
Return whether the line with index line
is oriented in standard direction. true
indicates, that the line is oriented from vertex 0 to vertex 1, whereas it is the other way around otherwise. In 1d and 2d, this is always true
, but in 3d it may be different, see the respective discussion in the documentation of the GeometryInfo class.
This function is really only for internal use in the library unless you absolutely know what this is all about.
This function queries ReferenceCell::standard_vs_true_line_orientation().

static 
Always return false.

static 
Always return false.

static 
Always return false.

static 
Always return false.

static 
Always return false.

static 
Always return false.

static 
Always return false.

static 
Always return false.