From ee88fce027f4842acc2e27ad55def74d18ce2f25 Mon Sep 17 00:00:00 2001 From: David Wells Date: Mon, 19 Apr 2021 15:30:11 -0400 Subject: [PATCH] Move the reordering text into a module. --- doc/doxygen/headers/grid.h | 6 +- doc/doxygen/headers/reordering.h | 615 +++++++++++++++++++++++++ include/deal.II/grid/grid_reordering.h | 595 +----------------------- 3 files changed, 620 insertions(+), 596 deletions(-) create mode 100644 doc/doxygen/headers/reordering.h diff --git a/doc/doxygen/headers/grid.h b/doc/doxygen/headers/grid.h index 7e321f3750..e0cb8ed79a 100644 --- a/doc/doxygen/headers/grid.h +++ b/doc/doxygen/headers/grid.h @@ -15,10 +15,10 @@ /** - * @defgroup grid Grid classes + * @defgroup grid Grids and Triangulations * - * This module groups classes that have to do with the topology and - * geometry of meshes. A mesh can be thought of as a collection of cells; + * This module groups functions and classes that have to do with the topology + * and geometry of meshes. A mesh can be thought of as a collection of cells; * if the mesh has been refined (possibly in an adaptive way), then * this collection is grouped into a hierarchy of refinement * levels. In addition to cells, the geometric objects that make up a diff --git a/doc/doxygen/headers/reordering.h b/doc/doxygen/headers/reordering.h new file mode 100644 index 0000000000..84d280cfa7 --- /dev/null +++ b/doc/doxygen/headers/reordering.h @@ -0,0 +1,615 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2002 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/** + * @defgroup reordering Grid reordering and cell orientation + * + * @brief A module describing how deal.II consistently orients Triangulation + * objects. + * + * This class reorders the vertices of cells such that they meet the standard + * requirements of the Triangulation class when creating grids, i.e. all lines + * have a unique orientation with respect to all neighboring cells. This class + * is mainly used when reading in grids from files and converting them to + * deal.II triangulations. + * + * @note In contrast to the rest of the deal.II library, by default this class + * uses the old deal.II numbering scheme, which was used up to deal.II version + * 5.2 (but the main function of this class takes a flag that specifies + * whether it should do an implicit conversion from the new to the old format + * before doing its work, and then back again after reordering). In this old + * format, the vertex and face ordering in 2d is assumed to be + * @verbatim + * 2 + * 3--->---2 + * | | + * 3^ ^1 + * | | + * 0--->---1 + * 0 + * @endverbatim + * the vertices in 3d: + * @verbatim + * 7-------6 7-------6 + * /| | / /| + * / | | / / | + * / | | / / | + * 3 | | 3-------2 | + * | 4-------5 | | 5 + * | / / | | / + * | / / | | / + * |/ / | |/ + * 0-------1 0-------1 + * @endverbatim + * and the faces in 3d: + * @verbatim + * *-------* *-------* + * /| | / /| + * / | 1 | / 4 / | + * / | | / / | + * * | | *-------* | + * | 5 *-------* | | 3 * + * | / / | | / + * | / 2 / | 0 | / + * |/ / | |/ + * *-------* *-------* + * @endverbatim + * After calling the GridReordering::reorder_cells() function the CellData is + * still in this old numbering scheme. Hence, for creating a Triangulation + * based on the resulting CellData the + * Triangulation::create_triangulation_compatibility() (and not the + * Triangulation::create_triangulation()) function must be used. For a + * typical use of the reorder_cells() function see the implementation of the + * GridIn read_*() functions. + * + * + *

Statement of problems

+ * + * Triangulations in deal.II have a special structure, in that there are not + * only cells, but also faces, and in 3d also edges, that are objects of their + * own right. Faces and edges have unique orientations, and they have a + * specified orientation also with respect to the cells that are adjacent. + * Thus, a line that separates two cells in two space dimensions does not only + * have a direction, but it must also have a well-defined orientation with + * respect to the other lines bounding the two quadrilaterals adjacent to the + * first line. Likewise definitions hold for three dimensional cells and the + * objects (lines, quads) that separate them. + * + * For example, in two dimensions, a quad consists of four lines which have a + * direction, which is by definition as follows: + * @verbatim + * 3-->--2 + * | | + * ^ ^ + * | | + * 0-->--1 + * @endverbatim + * 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 + * @verbatim + * 3---4---5 + * | | | + * 0---1---2 + * @endverbatim + * may be characterised by the vertex numbers (0 1 4 3) and (1 2 + * 5 4), since the middle line would get the direction 1->4 when + * viewed from both cells. The numbering (0 1 4 3) and (5 4 1 + * 2) would not be allowed, since the left quad would give the common + * line the direction 1->4, while the right one would want to use + * 4->1, leading to an ambiguity. + * + * 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: + * @verbatim + * o + * / \ + * o o + * | \ / | + * | o | + * | | | + * o---o---o + * @endverbatim + * (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. + * + * The purpose of this class is now to find an ordering for a given set of + * cells such that the generated triangulation satisfies all the requirements + * stated above. To this end, we will first show some examples why this is a + * difficult problem, and then develop algorithms that finds such a + * reordering. Note that the algorithm operates on a set of CellData objects + * that are used to describe a mesh to the triangulation class. These objects + * are, for example, generated by the GridIn class, when reading in grids from + * input files. + * + * 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 counter-example 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. + * + * + *

Examples of problems

+ * + * 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): + * @verbatim + * 9---10-----11 + * | | / | + * 6---7---8 | + * | | | | + * 3---4---5 | + * | | \ | + * 0---1------2 + * @endverbatim + * Assume that you had numbered the vertices in the cells at the left boundary + * in a way, that the following line directions are induced: + * @verbatim + * 9->-10-----11 + * ^ ^ / | + * 6->-7---8 | + * ^ ^ | | + * 3->-4---5 | + * ^ ^ \ | + * 0->-1------2 + * @endverbatim + * (This could for example be done by using the indices (0 1 4 3), + * (3 4 7 6), (6 7 10 9) 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 7 10): + * @verbatim + * 9->-10-----11 + * v v / | + * 6->-7---8 | + * ^ ^ | | + * 3->-4---5 | + * ^ ^ \ | + * 0->-1------2 + * @endverbatim + * + * The point here is the following: assume we wanted to prolong the grid to + * the left like this: + * @verbatim + * o---o---o---o---o------o + * | | | | | / | + * o---o---o---o---o---o | + * | | | | | | | + * o---o---o---o---o---o | + * | | | | | \ | + * o---o---o---o---o------o + * @endverbatim + * 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: + * @verbatim + * 3-----2-----o-----o ... o-----7-----6 + * | | | | | | | + * | 0 | N | N-1 | ... | 2 | 1 | + * | | | | | | | + * 0-----1-----o-----o ... o-----4-----5 + * @endverbatim + * 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 2 3 and 6 7 4 5. Then, if we follow this orientation, + * the grid after creating the lines for these two cells would look like this: + * @verbatim + * 3-->--2-----o-----o ... o-----7--<--6 + * | | | | | | | + * ^ 0 ^ N | N-1 | ... | 2 v 1 v + * | | | | | | | + * 0-->--1-----o-----o ... o-----4--<--5 + * @endverbatim + * Now, since opposite lines must point in the same direction, we can only add + * the cells 2 through N-1 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. + * + * + *

Solution

+ * + * 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 N-1 as well). + * + * A brute force approach to this problem is the following: if cell N can't be + * added, then try to rotate cell N-1. If we can't rotate cell N-1 any more, + * then try to rotate cell N-2 and try to add cell N with all orientations of + * cell N-1. 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 2 3), then the four + * possibilities would be (0 1 2 3), (1 2 3 0), (2 3 0 + * 1), and (3 0 1 2)). When adding one cell after the other, we + * traverse this tree in a depth-first (pre-order) 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 + * pessimized 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. + * + * + *

The 2d linear complexity algorithm

+ * + * 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 + * non-fixed 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: + * @verbatim + * 9------10-------11 + * | | /| + * | | / | + * | | / | + * 6------7-----8 | + * | | | | + * | | | | + * | | | | + * 3------4-----5 | + * | | \ | + * | | \ | + * | | \| + * 0------1---------2 + * @endverbatim + * First a cell is chosen ( (0,1,4,3) 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 edge-hopping and face hopping, giving a path through + * the mesh shown in dots. + * @verbatim + * 9-->--10-------11 + * | . | /| + * | . | / | + * | . | / | + * 6-->--7-----8 | + * | . | | | + * | . | | | + * | . | | | + * 3-->--4-----5 | + * | . | \ | + * | X | \ | + * | . | \| + * 0-->--1---------2 + * @endverbatim + * This is then repeated for the other sides of the chosen element, orienting + * more sides of the mesh. + * @verbatim + * 9-->--10-------11 + * | | /| + * v.....v.......V | + * | | /. | + * 6-->--7-----8 . | + * | | | . | + * | | | . | + * | | | . | + * 3-->--4-----5 . | + * | | \. | + * ^..X..^.......^ | + * | | \| + * 0-->--1---------2 + * @endverbatim + * 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. + * @verbatim + * 9-->--10->-----11 + * | ooo | . /| + * v ooo v . V | + * | ooo | . / | + * 6-->--7-->--8 | + * | | . | | + * | | . | | + * | | . | | + * 3-->--4-->--5 | + * | ooo | . \ | + * ^ ooo ^ X ^ | + * | ooo | . \| + * 0-->--1-->------2 + * @endverbatim + * Repeating this gives + * @verbatim + * 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 + * @endverbatim + * and the final oriented mesh is + * @verbatim + * 9-->--10->-----11 + * | | /| + * v v V | + * | | / | + * 6-->--7-->--8 | + * | | | | + * ^ ^ ^ ^ + * | | | | + * 3-->--4-->--5 | + * | | \ | + * ^ ^ ^ | + * | | \| + * 0-->--1-->-------2 + * @endverbatim + * It is obvious that this algorithm has linear run-time, since it only ever + * touches each face exactly once. + * + * The algorithm just described in the two-dimensional 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 three-dimensional 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. + * + * + *

For the curious

+ * + * Prior to the implementation of the algorithms described above (originally + * implemented by Michael Anderson in 2002, and re-implemented by Wolfgang + * Bangerth in 2016), we used a branch-and-cut 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 run-time 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 Cuthill-McKee 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 front-marching 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 N-1 unless N-1 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=N-1, the latter + * cannot be the case. Conversely, if we fail to add cell N, then it is not + * necessary to track back to cell N-1, 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): + * @verbatim + * 13----------14----15 + * | \ | | + * | \ 4 | 5 | + * | \ | | + * | 12-----10----11 + * | | | | + * | | | 7 | + * | | | | + * | 3 | 8-----9 + * | | | | + * | | | 6 | + * | | | | + * 4-----5-----6-----7 + * | | | | + * | 2 | 1 | 0 | + * | | | | + * 0-----1-----2-----3 + * @endverbatim + * 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 7 6, and cell 5 by + * 15 14 10 11, 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: + * @verbatim + * 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 + * @endverbatim + * 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. + * + * @ingroup grid + */ diff --git a/include/deal.II/grid/grid_reordering.h b/include/deal.II/grid/grid_reordering.h index 841608bf89..9a9386a889 100644 --- a/include/deal.II/grid/grid_reordering.h +++ b/include/deal.II/grid/grid_reordering.h @@ -34,599 +34,8 @@ DeclExceptionMsg(ExcMeshNotOrientable, /** - * This class reorders the vertices of cells such that they meet the standard - * requirements of the Triangulation class when creating grids, i.e. all lines - * have a unique orientation with respect to all neighboring cells. This class - * is mainly used when reading in grids from files and converting them to - * deal.II triangulations. - * - * @note In contrast to the rest of the deal.II library, by default this class - * uses the old deal.II numbering scheme, which was used up to deal.II version - * 5.2 (but the main function of this class takes a flag that specifies - * whether it should do an implicit conversion from the new to the old format - * before doing its work, and then back again after reordering). In this old - * format, the vertex and face ordering in 2d is assumed to be - * @verbatim - * 2 - * 3--->---2 - * | | - * 3^ ^1 - * | | - * 0--->---1 - * 0 - * @endverbatim - * the vertices in 3d: - * @verbatim - * 7-------6 7-------6 - * /| | / /| - * / | | / / | - * / | | / / | - * 3 | | 3-------2 | - * | 4-------5 | | 5 - * | / / | | / - * | / / | | / - * |/ / | |/ - * 0-------1 0-------1 - * @endverbatim - * and the faces in 3d: - * @verbatim - * *-------* *-------* - * /| | / /| - * / | 1 | / 4 / | - * / | | / / | - * * | | *-------* | - * | 5 *-------* | | 3 * - * | / / | | / - * | / 2 / | 0 | / - * |/ / | |/ - * *-------* *-------* - * @endverbatim - * After calling the GridReordering::reorder_cells() function the CellData is - * still in this old numbering scheme. Hence, for creating a Triangulation - * based on the resulting CellData the - * Triangulation::create_triangulation_compatibility() (and not the - * Triangulation::create_triangulation()) function must be used. For a - * typical use of the reorder_cells() function see the implementation of the - * GridIn read_*() functions. - * - * - *

Statement of problems

- * - * Triangulations in deal.II have a special structure, in that there are not - * only cells, but also faces, and in 3d also edges, that are objects of their - * own right. Faces and edges have unique orientations, and they have a - * specified orientation also with respect to the cells that are adjacent. - * Thus, a line that separates two cells in two space dimensions does not only - * have a direction, but it must also have a well-defined orientation with - * respect to the other lines bounding the two quadrilaterals adjacent to the - * first line. Likewise definitions hold for three dimensional cells and the - * objects (lines, quads) that separate them. - * - * For example, in two dimensions, a quad consists of four lines which have a - * direction, which is by definition as follows: - * @verbatim - * 3-->--2 - * | | - * ^ ^ - * | | - * 0-->--1 - * @endverbatim - * 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 - * @verbatim - * 3---4---5 - * | | | - * 0---1---2 - * @endverbatim - * may be characterised by the vertex numbers (0 1 4 3) and (1 2 - * 5 4), since the middle line would get the direction 1->4 when - * viewed from both cells. The numbering (0 1 4 3) and (5 4 1 - * 2) would not be allowed, since the left quad would give the common - * line the direction 1->4, while the right one would want to use - * 4->1, leading to an ambiguity. - * - * 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: - * @verbatim - * o - * / \ - * o o - * | \ / | - * | o | - * | | | - * o---o---o - * @endverbatim - * (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. - * - * The purpose of this class is now to find an ordering for a given set of - * cells such that the generated triangulation satisfies all the requirements - * stated above. To this end, we will first show some examples why this is a - * difficult problem, and then develop algorithms that finds such a - * reordering. Note that the algorithm operates on a set of CellData objects - * that are used to describe a mesh to the triangulation class. These objects - * are, for example, generated by the GridIn class, when reading in grids from - * input files. - * - * 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 counter-example 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. - * - * - *

Examples of problems

- * - * 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): - * @verbatim - * 9---10-----11 - * | | / | - * 6---7---8 | - * | | | | - * 3---4---5 | - * | | \ | - * 0---1------2 - * @endverbatim - * Assume that you had numbered the vertices in the cells at the left boundary - * in a way, that the following line directions are induced: - * @verbatim - * 9->-10-----11 - * ^ ^ / | - * 6->-7---8 | - * ^ ^ | | - * 3->-4---5 | - * ^ ^ \ | - * 0->-1------2 - * @endverbatim - * (This could for example be done by using the indices (0 1 4 3), - * (3 4 7 6), (6 7 10 9) 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 7 10): - * @verbatim - * 9->-10-----11 - * v v / | - * 6->-7---8 | - * ^ ^ | | - * 3->-4---5 | - * ^ ^ \ | - * 0->-1------2 - * @endverbatim - * - * The point here is the following: assume we wanted to prolong the grid to - * the left like this: - * @verbatim - * o---o---o---o---o------o - * | | | | | / | - * o---o---o---o---o---o | - * | | | | | | | - * o---o---o---o---o---o | - * | | | | | \ | - * o---o---o---o---o------o - * @endverbatim - * 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: - * @verbatim - * 3-----2-----o-----o ... o-----7-----6 - * | | | | | | | - * | 0 | N | N-1 | ... | 2 | 1 | - * | | | | | | | - * 0-----1-----o-----o ... o-----4-----5 - * @endverbatim - * 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 2 3 and 6 7 4 5. Then, if we follow this orientation, - * the grid after creating the lines for these two cells would look like this: - * @verbatim - * 3-->--2-----o-----o ... o-----7--<--6 - * | | | | | | | - * ^ 0 ^ N | N-1 | ... | 2 v 1 v - * | | | | | | | - * 0-->--1-----o-----o ... o-----4--<--5 - * @endverbatim - * Now, since opposite lines must point in the same direction, we can only add - * the cells 2 through N-1 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. - * - * - *

Solution

- * - * 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 N-1 as well). - * - * A brute force approach to this problem is the following: if cell N can't be - * added, then try to rotate cell N-1. If we can't rotate cell N-1 any more, - * then try to rotate cell N-2 and try to add cell N with all orientations of - * cell N-1. 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 2 3), then the four - * possibilities would be (0 1 2 3), (1 2 3 0), (2 3 0 - * 1), and (3 0 1 2)). When adding one cell after the other, we - * traverse this tree in a depth-first (pre-order) 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 - * pessimized 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. - * - * - *

The 2d linear complexity algorithm

- * - * 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 - * non-fixed 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: - * @verbatim - * 9------10-------11 - * | | /| - * | | / | - * | | / | - * 6------7-----8 | - * | | | | - * | | | | - * | | | | - * 3------4-----5 | - * | | \ | - * | | \ | - * | | \| - * 0------1---------2 - * @endverbatim - * First a cell is chosen ( (0,1,4,3) 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 edge-hopping and face hopping, giving a path through - * the mesh shown in dots. - * @verbatim - * 9-->--10-------11 - * | . | /| - * | . | / | - * | . | / | - * 6-->--7-----8 | - * | . | | | - * | . | | | - * | . | | | - * 3-->--4-----5 | - * | . | \ | - * | X | \ | - * | . | \| - * 0-->--1---------2 - * @endverbatim - * This is then repeated for the other sides of the chosen element, orienting - * more sides of the mesh. - * @verbatim - * 9-->--10-------11 - * | | /| - * v.....v.......V | - * | | /. | - * 6-->--7-----8 . | - * | | | . | - * | | | . | - * | | | . | - * 3-->--4-----5 . | - * | | \. | - * ^..X..^.......^ | - * | | \| - * 0-->--1---------2 - * @endverbatim - * 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. - * @verbatim - * 9-->--10->-----11 - * | ooo | . /| - * v ooo v . V | - * | ooo | . / | - * 6-->--7-->--8 | - * | | . | | - * | | . | | - * | | . | | - * 3-->--4-->--5 | - * | ooo | . \ | - * ^ ooo ^ X ^ | - * | ooo | . \| - * 0-->--1-->------2 - * @endverbatim - * Repeating this gives - * @verbatim - * 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 - * @endverbatim - * and the final oriented mesh is - * @verbatim - * 9-->--10->-----11 - * | | /| - * v v V | - * | | / | - * 6-->--7-->--8 | - * | | | | - * ^ ^ ^ ^ - * | | | | - * 3-->--4-->--5 | - * | | \ | - * ^ ^ ^ | - * | | \| - * 0-->--1-->-------2 - * @endverbatim - * It is obvious that this algorithm has linear run-time, since it only ever - * touches each face exactly once. - * - * The algorithm just described in the two-dimensional 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 three-dimensional 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. - * - * - *

For the curious

- * - * Prior to the implementation of the algorithms described above (originally - * implemented by Michael Anderson in 2002, and re-implemented by Wolfgang - * Bangerth in 2016), we used a branch-and-cut 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 run-time 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 Cuthill-McKee 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 front-marching 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 N-1 unless N-1 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=N-1, the latter - * cannot be the case. Conversely, if we fail to add cell N, then it is not - * necessary to track back to cell N-1, 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): - * @verbatim - * 13----------14----15 - * | \ | | - * | \ 4 | 5 | - * | \ | | - * | 12-----10----11 - * | | | | - * | | | 7 | - * | | | | - * | 3 | 8-----9 - * | | | | - * | | | 6 | - * | | | | - * 4-----5-----6-----7 - * | | | | - * | 2 | 1 | 0 | - * | | | | - * 0-----1-----2-----3 - * @endverbatim - * 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 7 6, and cell 5 by - * 15 14 10 11, 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: - * @verbatim - * 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 - * @endverbatim - * 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. - * - * @ingroup grid + * A class implementing various grid reordering algorithms. For more information + * see the @ref reordering "reordering module". */ template class GridReordering -- 2.39.5