From 6136ff40eaaa9c119698b9e4574b33b190c9c07c Mon Sep 17 00:00:00 2001 From: wolf Date: Sun, 1 Oct 2000 22:49:38 +0000 Subject: [PATCH] Documentation. git-svn-id: https://svn.dealii.org/trunk@3363 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/grid/grid_reordering.h | 376 ++++++++++++++++++ deal.II/deal.II/include/grid/tria.h | 70 ++-- 2 files changed, 420 insertions(+), 26 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_reordering.h b/deal.II/deal.II/include/grid/grid_reordering.h index 02dad5b89c..289c1f0aa2 100644 --- a/deal.II/deal.II/include/grid/grid_reordering.h +++ b/deal.II/deal.II/include/grid/grid_reordering.h @@ -86,6 +86,382 @@ class GridReorderingInfo<3> /** + * This class reorders the vertices of cells such that they meet the + * requirements of the @ref{Triangulation} class when creating + * grids. This class is mainly used when reading in grids from files + * and converting them to deal.II triangulations. + * + * + * @sect3{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: + * @begin{verbatim} + * 3-->--2 + * | | + * ^ ^ + * | | + * 0-->--1 + * @end{verbatim} + * 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 + * @begin{verbatim} + * 3---4---5 + * | | | + * 0---1---2 + * @end{verbatim} + * may be characterised by the vertex numbers @p{(0 1 4 3)} and + * @p{(1 2 5 4)}, since the middle line would get the direction @p{1->4} + * when viewed from both cells. The numbering @p{(0 1 4 3)} and + * @p{(5 4 1 2)} would not be allowed, since the left quad would give the + * common line the direction @p{1->4}, while the right one would want + * to use @p{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: + * @begin{verbatim} + * o + * / \ + * o o + * | \ / | + * | o | + * | | | + * o---o---o + * @end{verbatim} + * (The reader is aked 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 an + * algorithm that finds such a reordering. Note that the algorithm + * operates on a set of @ref{CellData} objects that are used to + * describe a mesh to the triangulation class. These objects are, for + * example, generated by the @ref{GridIn} class, when reading in grids + * from input files. + * + * + * @sect3{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): + * @begin{verbatim} + * 9---10-----11 + * | | / | + * 6---7---8 | + * | | | | + * 3---4---5 | + * | | \ | + * 0---1------2 + * @end{verbatim} + * Assume that you had numbered the vertices in the cells at the left boundary + * in a way, that the following line directions are induced: + * @begin{verbatim} + * 9->-10-----11 + * ^ ^ / | + * 6->-7---8 | + * ^ ^ | | + * 3->-4---5 | + * ^ ^ \ | + * 0->-1------2 + * @end{verbatim} + * (This could for example be done by using the indices @p{(0 1 4 3)}, @p{(3 4 7 6)}, + * @p{(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 @p{(9 6 7 10)}: + * @begin{verbatim} + * 9->-10-----11 + * v v / | + * 6->-7---8 | + * ^ ^ | | + * 3->-4---5 | + * ^ ^ \ | + * 0->-1------2 + * @end{verbatim} + * + * The point here is the following: assume we wanted to prolong the grid to + * the left like this: + * @begin{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 + * @end{verbatim} + * 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 @p{(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 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: + * @begin{verbatim} + * 3-----2-----o-----o ... o-----7-----6 + * | | | | | | | + * | 0 | N | N-1 | ... | 2 | 1 | + * | | | | | | | + * 0-----1-----o-----o ... o-----4-----5 + * @end{verbatim} + * 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 @p{0 1 2 3} and @p{6 7 4 5}. Then, if we follow this + * orientation, the grid after creating the lines for these two cells + * would look like this: + * @begin{verbatim} + * 3-->--2-----o-----o ... o-----7--<--6 + * | | | | | | | + * ^ 0 ^ N | N-1 | ... | 2 v 1 v + * | | | | | | | + * 0-->--1-----o-----o ... o-----4--<--5 + * @end{verbatim} + * 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 oeirntation of faces and can't add the next + * one without introducingfaces that had already been added in another + * direction, then it might not be sufficient to only rotate cells in + * the neighborhood of the 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. + * + * + * @sect3{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 oreintation. 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). + * + * The only solution to this problem seems to be the following: if + * cell N can't be added, the 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 indicaes are @p{(0 1 2 3)}, + * then the four possibilities would be @p{(0 1 2 3)}, @p{(1 2 3 0)}, + * @p{(2 3 0 1)}, and @p{(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 @p{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 + * @p{4**N}. However, an 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 geatly improved by the fact that we + * can find an algorithm that in practice is usually only roughly + * linear in time and memory. We will describe this algorithm in the + * following. + * + * The first observation is that although there are counterexamples, + * problems are usually local. For example, in the second example, if + * we had numbered the cells in a way that neighboring cells have + * similar cell numbers, then the amount pf 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 + * abvious: 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 + * 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 -- somtimes + * 10 or more -- neighbors of one vertex, and exposed several problems + * in the initial algorithm): + * @begin{verbatim} + * 12----13----14----15 + * | | | | + * | 5 | 7 | 6 | + * | | | | + * 8-----9-----10----11 + * | | | | + * | 3 | | 4 | + * | | | | + * 4-----5-----6-----7 + * | | | | + * | 0 | 1 | 2 | + * | | | | + * 0-----1-----2-----3 + * @end{verbatim} + */ +/* + * Note that there is a hole in the middle. Assume now that the user + * described the first cell 0 by the vertex numbers @p{0 1 5 4}, cell + * 5 by @p{12 8 9 13}, and cell 6 by @p{10 11 15 14}. All other cells + * are numbered in the usual way, i.e. starting at the bottom left and + * counting counterclockwise. Cell 5 therefore is the only one that + * does not follow this order; however, note that the bottom line of + * cell 5 given by this order of cell 5 does match with the top line + * of cell 4 in that orientation. 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: + * @begin{verbatim} + * 12->--13----14->--15 + * | | | | + * v 5 v 7 ^ 6 ^ + * | | | | + * 8-->--9-----10->--11 + * | | | | + * ^ 3 ^ ^ 4 ^ + * | | | | + * 4-->--5-->--6-->--7 + * | | | | + * ^ 0 ^ 1 ^ 2 ^ + * | | | | + * 0-->--1-->--2-->--3 + * @end{verbatim} + * + * Coming now to cell 7, we see that the two + * opposite lines to its left and right 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...... + * + * 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 backtrack 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 + * @ref{Triangulation}@p{::create_triangulation} function. + * * @author Wolfgang Bangerth, 2000 */ template diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index 009ff436c3..2134e36d24 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -683,18 +683,22 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * consists of a vector storing the indices of the vertices of this cell * in the vertex list. To see how this works, you can take a look at the * @ref{GridIn}@p{::read_*} functions. The appropriate function to be - * called is @p{Triangulation::create_triangulation (2)}. - * - * Creating the hierarchical information needed for this library from - * cells storing only vertex information can be quite a complex task. - * For example in 2d, we have to create lines between vertices (but only - * once, though there are two cells which link these two vertices) and - * we have to create neighborship information. Grids being read in - * should therefore not be too large, reading refined grids would be - * inefficient. Apart from the performance aspect, refined grids do not - * lend too well to multigrid algorithms, since solving on the coarsest - * level is expensive. It is wiser in any case to read in a grid as coarse - * as possible and then do the needed refinement steps. + * called is @p{Triangulation::create_triangulation ()}. + * + * Creating the hierarchical information needed for this + * library from cells storing only vertex information can be + * quite a complex task. For example in 2d, we have to create + * lines between vertices (but only once, though there are two + * cells which link these two vertices) and we have to create + * neighborship information. Grids being read in should + * therefore not be too large, reading refined grids would be + * inefficient (although there is technically no problem in + * reading grids with several 10.000 or 100.000 cells; the + * library can handle this without much problems). Apart from + * the performance aspect, refined grids do not lend too well + * to multigrid algorithms, since solving on the coarsest level + * is expensive. It is wiser in any case to read in a grid as + * coarse as possible and then do the needed refinement steps. * * It is your duty to guarantee that cells have the correct orientation. * To guarantee this, in the input vector keeping the cell list, the @@ -709,12 +713,16 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * quadrilaterals, i.e. two vertices interchanged; this results in * a wrong area element). * - * There are more subtle conditions which must be imposed upon the - * vertex numbering within cells. See the documentation for the - * @ref{GridIn} class for more details on this. They do not only - * hold for the data read from an UCD or any other input file, but - * also for the data passed to the - * @p{Triangulation::create_triangulation (2)} function. + * There are more subtle conditions which must be imposed upon + * the vertex numbering within cells. They do not only hold for + * the data read from an UCD or any other input file, but also + * for the data passed to the + * @p{Triangulation::create_triangulation ()} + * function. See the documentation for the @ref{GridIn} class + * for more details on this, and above all to the + * @ref{GridReordering} class that explains many of the + * problems and an algorithm to reorder cells such that they + * satisfy the conditions outlined above. * * @item Copying a triangulation: when computing on time dependant meshes * of when using adaptive refinement, you will often want to create a @@ -1767,6 +1775,13 @@ class Triangulation : public TriaDimensionInfo, * @p{virtual} to allow derived * classes to set up some data * structures as well. + * + * For conditions when this + * function can generate a valid + * triangulation, see the + * documentation of this class, + * and the @ref{GridIn} and + * @ref{GridReordering} class. */ virtual void create_triangulation (const vector > &vertices, const vector > &cells, @@ -3002,15 +3017,18 @@ class Triangulation : public TriaDimensionInfo, void fix_coarsen_flags (); /** - * Re-compute the number of lines, quads, - * etc. This function is called by - * @p{execute_{coarsening,refinement}} and - * by @p{create_triangulation} after the - * grid was changed. + * Re-compute the number of + * lines, quads, etc. This + * function is called by + * @p{execute_coarsening_and_refinement} + * and by + * @p{create_triangulation} after + * the grid was changed. * - * This function simply delegates to the - * functions below, which count - * only a certain class of objects. + * This function simply delegates + * to the functions below, which + * count only a certain class of + * objects. */ void update_number_cache (); -- 2.39.5