]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Documentation.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 1 Oct 2000 22:49:38 +0000 (22:49 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 1 Oct 2000 22:49:38 +0000 (22:49 +0000)
git-svn-id: https://svn.dealii.org/trunk@3363 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_reordering.h
deal.II/deal.II/include/grid/tria.h

index 02dad5b89cf163d44f05916263c3585d6470d5dc..289c1f0aa2f4bddb5eb62d3bf616a3c870ec2d58 100644 (file)
@@ -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{<dim>::create_triangulation} function.
+ *
  * @author Wolfgang Bangerth, 2000
  */
 template <int dim>
index 009ff436c30dd4274f32281f0de03cf3465bdef7..2134e36d24063b69d44e6c4e27feaa82b589c18a 100644 (file)
@@ -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{<dim>::read_*} functions. The appropriate function to be
- *        called is @p{Triangulation<dim>::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<dim>::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<dim>::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<dim>::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<dim>,
                                      * @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<Point<dim> >    &vertices,
                                       const vector<CellData<dim> > &cells,
@@ -3002,15 +3017,18 @@ class Triangulation : public TriaDimensionInfo<dim>,
     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 ();
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.