]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a class that implements reordering of grids whenever this is
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 27 Sep 2000 15:26:17 +0000 (15:26 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 27 Sep 2000 15:26:17 +0000 (15:26 +0000)
necessary to fulfill the assumptions of the triangulation class.

git-svn-id: https://svn.dealii.org/trunk@3351 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_reordering.h [new file with mode: 0644]
deal.II/deal.II/source/grid/grid_reordering.cc [new file with mode: 0644]

diff --git a/deal.II/deal.II/include/grid/grid_reordering.h b/deal.II/deal.II/include/grid/grid_reordering.h
new file mode 100644 (file)
index 0000000..67f7c8c
--- /dev/null
@@ -0,0 +1,448 @@
+//----------------------------  grid_reordering.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  grid_reordering.h  ---------------------------
+#ifndef __deal2__grid_reordering_h
+#define __deal2__grid_reordering_h
+
+
+#include <map>
+#include <vector>
+#include <grid/tria.h>
+
+
+/**
+ * Class declaring some dimension dependent numbers which are needed
+ * for the grid reordering class.
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+template <int dim>
+class GridReorderingInfo
+{};
+
+
+/**
+ * Class declaring some dimension dependent numbers which are needed
+ * for the grid reordering class. This is the specialization for the
+ * 2d case.
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+template <>
+class GridReorderingInfo<2>
+{
+  public:
+                                    /**
+                                     * Number of possible valid
+                                     * orientations of a cell. They
+                                     * are the state in which it was
+                                     * delivered and three possible
+                                     * rotations in counter-clockwise
+                                     * sense, thus a total of four.
+                                     */
+    static const unsigned int rotational_states_of_cells = 4;
+
+                                    /**
+                                     * Number of possible
+                                     * orientations of a face in
+                                     * 2d. It is the face and the
+                                     * face with vertices exchanged,
+                                     * thus two.
+                                     */
+    static const unsigned int rotational_states_of_faces = 2;
+};
+
+
+
+/**
+ * Class declaring some dimension dependent numbers which are needed
+ * for the grid reordering class. This is the specialization for the
+ * 2d case.
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+template <>
+class GridReorderingInfo<3>
+{
+  public:
+                                    /**
+                                     * ???
+                                     */
+    static const unsigned int rotational_states_of_cells = static_cast<unsigned int>(-1);
+    static const unsigned int rotational_states_of_faces = 4;
+};
+
+
+  
+
+
+/**
+ * @author Wolfgang Bangerth, 2000
+ */
+template <int dim>
+class GridReordering : private GridReorderingInfo<dim>
+{
+  public:
+                                    /**
+                                     * This is the main function,
+                                     * doing what is announced in the
+                                     * general documentation of this
+                                     * class.
+                                     */
+    static void reorder_cells (vector<CellData<dim> > &original_cells);
+    
+  private:
+
+                                    /**
+                                     * Forward declarations of local
+                                     * classes.
+                                     */
+    class Cell;
+    class Face;
+    class FaceData;
+    
+                                    /**
+                                     * Class that describes the
+                                     * properties of cells beyond
+                                     * what is provided by the data
+                                     * that is available from the
+                                     * calling functions of this
+                                     * class. In particular, several
+                                     * fields are available that
+                                     * describe connections of cells
+                                     * to faces and to
+                                     * neighbors. These fields are
+                                     * filled in a first pass before
+                                     * the actual reoordering starts,
+                                     * as they are needed for the
+                                     * latter purpose.
+                                     *
+                                     * Since this class is derived
+                                     * from the @ref{CellData} class,
+                                     * it also contains all the
+                                     * information available
+                                     * beforehand.
+                                     *
+                                     * @author Wolfgang Bangerth, 2000
+                                     */
+    struct Cell : public CellData<dim>
+    {
+                                        /**
+                                         * Value to be used if a
+                                         * neighbor does not exist,
+                                         * i.e. if the cell is at the
+                                         * boundary of the domain
+                                         * with a certain face.
+                                         */
+       static const unsigned int invalid_neighbor = static_cast<unsigned int>(-1);
+
+                                        /**
+                                         * Pointers to the faces of
+                                         * this cell and their
+                                         * rotations. If the first
+                                         * index is zero, then the
+                                         * faces denote the faces in
+                                         * their standard direction
+                                         * with respect to the
+                                         * ordering in this cell. If
+                                         * it is nonzero, then they
+                                         * denote the faces that
+                                         * would be needed if the
+                                         * cell were rotate so often.
+                                         */
+       typename map<Face,FaceData>::iterator
+       faces[GridReorderingInfo<dim>::rotational_states_of_cells][GeometryInfo<dim>::faces_per_cell];
+
+                                        /**
+                                         * Cell indices of the
+                                         * neighbors of this cell in
+                                         * the global array of cells.
+                                         */
+       unsigned int neighbors[GeometryInfo<dim>::faces_per_cell];
+
+                                        /**
+                                         * The index of this cell in
+                                         * the global array of cells.
+                                         */
+       unsigned int cell_no;
+
+                                        /**
+                                         * If we fail to insert this
+                                         * cell, then we have to
+                                         * track back. We could track
+                                         * back right to the previous
+                                         * cell, but we can do better
+                                         * than that, by tracking
+                                         * back to the cell indicated
+                                         * by this field. Which value
+                                         * it has is described in the
+                                         * documentation of the
+                                         * @ref{GridReordering}
+                                         * class.
+                                         */
+       unsigned int track_back_to_cell;
+
+                                        /**
+                                         * Default
+                                         * constructor. Invalidate
+                                         * all data.
+                                         */
+       Cell ();
+       
+                                        /**
+                                         * Constructor that copies
+                                         * the data of an object of
+                                         * the base class and
+                                         * requires to be given the
+                                         * index of this cell in the
+                                         * global array.
+                                         */
+       Cell (const CellData<dim> &cd,
+             const unsigned int   cell_no);
+       
+                                        /**
+                                         * Count the existing neighbors
+                                         * of this cell.
+                                         */
+       unsigned int count_neighbors () const;
+
+                                        /**
+                                         * Insert the faces of the
+                                         * present cell into the map
+                                         * of all faces. This
+                                         * function inserts them in
+                                         * all orientations possible
+                                         * if the given cell is
+                                         * rotated. The function also
+                                         * takes care to fill in the
+                                         * @p{adjacent_cells} field
+                                         * of the inserted faces.
+                                         */
+       void insert_faces (map<Face,FaceData > &global_faces);
+
+                                        /**
+                                         * Find out the neighbors of the
+                                         * given cell by looking at the
+                                         * @p{adjacent_cells} field of
+                                         * the faces of this cell. Store
+                                         * the neighbor indices in the
+                                         * present object.
+                                         */
+       void fix_cell_neighbors ();
+
+                                        /**
+                                         * Compute back to which cell we
+                                         * have to backtrack in case we
+                                         * can't insert this cell in any
+                                         * orientation into the already
+                                         * existing part of the
+                                         * triangulation. The method of
+                                         * how to determine the point to
+                                         * which we have to backtrack is
+                                         * described in the documentation
+                                         * of the @ref{GridReordering}
+                                         * class.
+                                         */
+       void find_backtracking_point ();
+
+                                        /**
+                                         * Find out whether the cell
+                                         * could be inserted into the
+                                         * already existing part of
+                                         * the triangulation with
+                                         * orientation given by the
+                                         * parameter, by checking
+                                         * that no face would be
+                                         * inserted twice in
+                                         * different orientations.
+                                         */
+       bool check_consistency (const unsigned int rot) const;
+
+                                        /**
+                                         * Tag the faces of this cell in
+                                         * the given orientation as used
+                                         * by this cell.
+                                         */
+       void mark_faces_used (const unsigned int rot);
+
+                                        /**
+                                         * Remove the use tags on the
+                                         * faces of this cell in the
+                                         * given orientation by this
+                                         * cell. Tags may remain if
+                                         * there is another cell that
+                                         * uses a given face.
+                                         */
+       void mark_faces_unused (const unsigned int rot);
+    };
+
+
+                                    /**
+                                     * Structure describing a face of
+                                     * a cell. This class is used as
+                                     * key in a map storing all faces
+                                     * possible in a triangulation,
+                                     * i.e. all faces between cells
+                                     * in all possible orientations.
+                                     *
+                                     * @author Wolfgang Bangerth, 2000
+                                     */
+    struct Face
+    {
+                                        /**
+                                         * Indices of the vertices of
+                                         * this face.
+                                         */
+       unsigned int vertices[GeometryInfo<dim>::vertices_per_face];
+
+                                        /**
+                                         * Comparison operator. Use
+                                         * the vertex indices as
+                                         * primary, secondary,
+                                         * ... criteria for
+                                         * comparison.
+                                         */
+       bool operator < (const Face &face) const;
+    };
+
+
+                                    /**
+                                     * Class describing some data to
+                                     * be stored on each
+                                     * face. Objects of this type are
+                                     * used as values in a map
+                                     * containing all faces.
+                                     *
+                                     * @author Wolfgang Bangerth, 2000
+                                     */
+    struct FaceData
+    {
+                                        /**
+                                         * Value denoting
+                                         * non-existing adjacent
+                                         * cells of this face,
+                                         * i.e. when the face is at
+                                         * the boundary of the domain
+                                         * and has only one adjacent
+                                         * cell.
+                                         */
+       static const unsigned int invalid_adjacent_cell = static_cast<unsigned int>(-1);
+
+                                        /**
+                                         * Pointers to the same face
+                                         * but in all other
+                                         * orientations. Storing
+                                         * these pointers makes it
+                                         * much easier to find out
+                                         * whether a given faces has
+                                         * already been used in
+                                         * another direction, thus
+                                         * forbidding the present
+                                         * face to be used.
+                                         */
+       typename map<Face,FaceData >::const_iterator
+       reverse_faces[GridReorderingInfo<dim>::rotational_states_of_faces-1];
+
+                                        /**
+                                         * Indices of the one or two
+                                         * adjacent cells of this
+                                         * face in the global array
+                                         * of cells.
+                                         */
+       unsigned int adjacent_cells[2];
+
+                                        /**
+                                         * Number of cells presently
+                                         * using this face in the
+                                         * orientation represented by
+                                         * this object. May be zero,
+                                         * one, or two.
+                                         */
+       unsigned int use_count;
+
+                                        /**
+                                         * Default constructor.
+                                         */
+       FaceData ();
+    };
+
+
+                                    /**
+                                     * If we couldn't insert a cell
+                                     * into the already existing part
+                                     * of the mesh, then we need to
+                                     * track back a while. This
+                                     * function does so, given the
+                                     * array of cells, the stack of
+                                     * rotation states, and the
+                                     * position to which backtracking
+                                     * shall take place.
+                                     *
+                                     * In some cases, it is possible
+                                     * that from the place where we
+                                     * backtracked to, there is no
+                                     * more possibility to orient a
+                                     * cell. Then we will have to
+                                     * backtrack further until we
+                                     * come to a place where further
+                                     * work is possible; this
+                                     * recursive backtracking is also
+                                     * done by this function,
+                                     * although it is not implemented
+                                     * as recursive calls but rather
+                                     * as eliminated tail-recursion.
+                                     */
+    static void track_back (vector<Cell>                               &cells,
+                           stack<unsigned int, vector<unsigned int> > &rotation_states,
+                           unsigned int                                track_back_to_cell);
+
+                                    /**
+                                     * This is the main function that
+                                     * does the main work. It is
+                                     * called by the
+                                     * @p{reorder_cells} function
+                                     * after all the preparations
+                                     * have been completed and
+                                     * operates on the @p{cells}
+                                     * array. After a way to reorder
+                                     * the cells has been found, the
+                                     * @p{original_cells} are reorder
+                                     * accordingly, where the
+                                     * @p{new_cell_numbers} array is
+                                     * needed to find the connection
+                                     * between original cells and
+                                     * presorted cells.
+                                     */
+    static void find_reordering (vector<Cell>               &cells,
+                                vector<CellData<dim> >     &original_cells,
+                                const vector<unsigned int> &new_cell_numbers);
+
+                                    /**
+                                     * Preorder the incoming cells by
+                                     * some kind of Cuthill-McKee
+                                     * algorithm. The reason for the
+                                     * need to do so is described in
+                                     * the general documentation.
+                                     *
+                                     * Return a vector in which for
+                                     * each old cell the new index is
+                                     * stored.
+                                     */
+    static
+    vector<unsigned int>
+    presort_cells (vector<Cell>       &cells,
+                  map<Face,FaceData> &faces);
+};
+
+
+
+#endif
diff --git a/deal.II/deal.II/source/grid/grid_reordering.cc b/deal.II/deal.II/source/grid/grid_reordering.cc
new file mode 100644 (file)
index 0000000..a5e9c9d
--- /dev/null
@@ -0,0 +1,805 @@
+//----------------------------  grid_reordering.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  grid_reordering.cc  ---------------------------
+
+
+#include <algorithm>
+#include <stack>
+
+#include <grid/grid_reordering.h>
+
+
+
+
+template <int dim>
+GridReordering<dim>::Cell::Cell () :
+               cell_no (invalid_neighbor)
+{
+  for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+    neighbors[i] = invalid_neighbor;
+};
+
+
+
+template <int dim>
+GridReordering<dim>::Cell::Cell (const CellData<dim> &cd,
+                                const unsigned int   cell_no) :
+               CellData<dim> (cd), cell_no(cell_no)
+{
+  for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+    neighbors[i] = invalid_neighbor;
+};
+
+
+
+template <int dim>
+inline
+unsigned int GridReordering<dim>::Cell::count_neighbors () const
+{
+  unsigned int n = 0;
+  for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+    if (neighbors[i] != invalid_neighbor)
+      ++n;
+  return n;
+};
+
+
+
+template <int dim>
+void
+GridReordering<dim>::Cell::insert_faces (map<Face,FaceData> &/*global_faces*/)
+{
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <>
+void
+GridReordering<2>::Cell::insert_faces (map<Face,FaceData> &global_faces)
+{
+  const unsigned int dim = 2;
+
+                                  // first compute index numbers for
+                                  // the faces in usual order as
+                                  // defined by the order of vertices
+                                  // in the cell object
+  Face new_faces[GeometryInfo<dim>::faces_per_cell]
+    = { { {vertices[0], vertices[1]} },
+        { {vertices[1], vertices[2]} },
+       { {vertices[3], vertices[2]} },
+       { {vertices[0], vertices[3]} } };
+
+                                  // then insert them into the global
+                                  // list and store iterators to
+                                  // them. note that if the face
+                                  // already exists, then the stored
+                                  // data is not touched.
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    faces[0][face] = global_faces.insert (make_pair(new_faces[face],
+                                                   FaceData())).first;
+
+
+                                  // then for each of the faces also
+                                  // insert the reverse form and
+                                  // store pointers to them. note
+                                  // that the rotational state in
+                                  // which all faces are reverted is
+                                  // `2'
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    {
+      swap (new_faces[face].vertices[0],
+           new_faces[face].vertices[1]);
+      faces[2][face] = global_faces.insert (make_pair(new_faces[face],
+                                                     FaceData())).first;
+    };
+
+                                  // then finally fill in rotational
+                                  // states 1 and 3 of the cell. the
+                                  // faces of these states can be
+                                  // obtained from states 0 and 2
+  faces[1][0] = faces[2][0];
+  faces[1][1] = faces[0][1];
+  faces[1][2] = faces[2][2];
+  faces[1][3] = faces[0][3];
+  
+  faces[3][0] = faces[0][0];
+  faces[3][1] = faces[2][1];
+  faces[3][2] = faces[0][2];
+  faces[3][3] = faces[2][3];
+  
+
+                                  // finally fill the crosslink and
+                                  // other fields of the new
+                                  // entries. note that since
+                                  // rotational states 0 and 2 of the
+                                  // cell are exactly reverted, we
+                                  // only have to operate on the face
+                                  // pointers of these two states to
+                                  // reach all possible faces and
+                                  // permutations thereof
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    {
+      if (faces[0][face]->second.adjacent_cells[0] ==
+         FaceData::invalid_adjacent_cell)
+       {
+                                          // face had not been
+                                          // inserted by previous
+                                          // cells, since first
+                                          // adjacent cell is still
+                                          // untouched. provide
+                                          // xlinks to rotated faces
+         faces[0][face]->second.reverse_faces[0] = faces[2][face];
+         faces[2][face]->second.reverse_faces[0] = faces[0][face];
+
+                                          // and insert this cell as
+                                          // adjacent_cell of the faces
+         faces[0][face]->second.adjacent_cells[0] = cell_no;
+         faces[2][face]->second.adjacent_cells[0] = cell_no;
+       }
+      else
+       {
+                                          // face had already been
+                                          // inserted. make sure that
+                                          // it was in the same way:
+         Assert (faces[0][face]->second.reverse_faces[0] == faces[2][face],
+                 ExcInternalError());    
+         Assert (faces[2][face]->second.reverse_faces[0] == faces[0][face],
+                 ExcInternalError());
+
+                                          // now insert ourselves as
+                                          // second
+                                          // adjacent_cell. the
+                                          // respective slots must
+                                          // necessarily be empty
+                                          // still
+         Assert (faces[0][face]->second.adjacent_cells[1] ==
+                 FaceData::invalid_adjacent_cell,
+                 ExcInternalError());
+         Assert (faces[2][face]->second.adjacent_cells[1] ==
+                 FaceData::invalid_adjacent_cell,
+                 ExcInternalError());
+         faces[0][face]->second.adjacent_cells[1] = cell_no;
+         faces[2][face]->second.adjacent_cells[1] = cell_no;
+       };
+    };
+};
+
+
+
+template <int dim>
+void GridReordering<dim>::Cell::fix_cell_neighbors ()
+{
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    {
+                                      // first assert that the
+                                      // neighborship info of all
+                                      // versions of the same face is
+                                      // identical
+      for (unsigned int rot=1; rot<rotational_states_of_faces; ++rot)
+       for (unsigned int adjacent_cell=0; adjacent_cell<2; ++adjacent_cell)
+         Assert (faces[rot][face]->second.adjacent_cells[adjacent_cell]
+                 ==
+                 faces[0][face]->second.adjacent_cells[adjacent_cell],
+                 ExcInternalError());
+      
+
+                                      // then insert the neighbor
+                                      // behind this face as neighbor
+                                      // of the present cell. note
+                                      // that it is not relevant to
+                                      // which permutation of a face
+                                      // we refer. note that it might
+                                      // well be that some of the
+                                      // neighbor indices are
+                                      // FaceData::invalid_adjacent_cell
+      if (faces[0][face]->second.adjacent_cells[0] == cell_no)
+       neighbors[face] = faces[0][face]->second.adjacent_cells[1];
+      else
+       neighbors[face] = faces[0][face]->second.adjacent_cells[0];
+    };
+};
+
+
+
+template <int dim>
+void GridReordering<dim>::Cell::find_backtracking_point ()
+{
+                                  // we know what neighbors we have,
+                                  // we can determine the neighbor
+                                  // with the maximal cell_no that is
+                                  // smaller than that of the present
+                                  // cell. we need this information
+                                  // in the backtracking process and
+                                  // don't want to compute it every
+                                  // time again
+  track_back_to_cell = FaceData::invalid_adjacent_cell;
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    if ((neighbors[face] != FaceData::invalid_adjacent_cell)
+       &&
+       (neighbors[face] < cell_no)
+       &&
+       ((neighbors[face] > track_back_to_cell)
+        ||
+        (track_back_to_cell == FaceData::invalid_adjacent_cell)))
+      track_back_to_cell = neighbors[face];
+
+                                  // if this cell had no neighbors
+                                  // with lower cell numbers, we
+                                  // still need to know what cell to
+                                  // track back to in case some
+                                  // higher cell than the present one
+                                  // failed to coexist with the
+                                  // existing part of the mesh
+                                  // irrespective of the rotation
+                                  // state of this present cell. we
+                                  // then simply track back to the
+                                  // cell before this one, lacking a
+                                  // better alternative. this does,
+                                  // of course, not hold for cell 0,
+                                  // from which we should never be
+                                  // forced to track back
+  if (track_back_to_cell == 0)
+    track_back_to_cell = 0;
+  else
+    if (track_back_to_cell == FaceData::invalid_adjacent_cell)
+      track_back_to_cell = cell_no-1;
+};
+
+
+
+template <int dim>
+inline
+bool GridReordering<dim>::Cell::check_consistency (const unsigned int /*rot*/) const
+{
+                                  // might be that we can use the
+                                  // specialization for dim==2 after
+                                  // some modification, but haven't
+                                  // thought about that yet
+  Assert (false, ExcNotImplemented());
+  return false;
+};
+
+
+
+template <>
+inline
+bool GridReordering<2>::Cell::check_consistency (const unsigned int rot) const
+{
+  const unsigned int dim = 2;
+                                  // make sure that for each face of
+                                  // the cell the permuted faces are
+                                  // not already in use, as that
+                                  // would make the cell disallowed
+  for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+    {
+      const FaceData &face = faces[rot][face_no]->second;
+
+      for (unsigned int face_rot=0; face_rot<rotational_states_of_faces-1; ++face_rot)
+       {
+         const FaceData &reverse_face = face.reverse_faces[face_rot]->second;
+         if (reverse_face.use_count != 0)
+           return false;
+       };
+    };
+
+                                  // no conflicts found
+  return true;
+};
+
+
+template <int dim>
+inline
+void GridReordering<dim>::Cell::mark_faces_used (const unsigned int /*rot*/)
+{
+                                  // might be that we can use the
+                                  // specialization for dim==2 after
+                                  // some modification, but haven't
+                                  // thought about that yet. in
+                                  // particular, I don't know whether
+                                  // we have to treat edges in 3d
+                                  // specially
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <int dim>
+inline
+void GridReordering<dim>::Cell::mark_faces_unused (const unsigned int /*rot*/)
+{
+                                  // might be that we can use the
+                                  // specialization for dim==2 after
+                                  // some modification, but haven't
+                                  // thought about that yet. in
+                                  // particular, I don't know whether
+                                  // we have to treat edges in 3d
+                                  // specially
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+template <>
+inline
+void GridReordering<2>::Cell::mark_faces_used (const unsigned int rot)
+{
+  const unsigned int dim=2;
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    {
+      Assert (faces[rot][face]->second.use_count < 2,
+             ExcInternalError());
+      ++faces[rot][face]->second.use_count;
+    };
+};
+
+
+
+template <>
+inline
+void GridReordering<2>::Cell::mark_faces_unused (const unsigned int rot)
+{
+  const unsigned int dim=2;
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    {
+      Assert (faces[rot][face]->second.use_count > 0,
+             ExcInternalError());
+      --faces[rot][face]->second.use_count;
+    };
+};
+
+
+
+template <int dim>
+bool GridReordering<dim>::Face::operator < (const Face &face) const
+{
+  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+    {
+                                      // if vertex index is smaller,
+                                      // then comparison is true
+      if (vertices[v] < face.vertices[v])
+       return true;
+      else
+                                        // if vertex index is greater,
+                                        // then comparison is false
+       if (vertices[v] > face.vertices[v])
+         return false;
+                                      // if indices are equal, then test
+                                      // next index
+    };
+
+                                  // if all indices are equal:
+  return false;
+};
+
+
+
+template <int dim>
+GridReordering<dim>::FaceData::FaceData () :
+               use_count (0)
+{
+  adjacent_cells[0] = adjacent_cells[1] = invalid_adjacent_cell;
+};
+
+
+
+
+
+
+template <int dim>
+inline
+void GridReordering<dim>::track_back (vector<Cell<dim> >                         &cells,
+                                     stack<unsigned int, vector<unsigned int> > &rotation_states,
+                                     unsigned int                                track_back_to_cell)
+{
+  top_of_function:
+  
+  Assert (track_back_to_cell > 0, ExcInternalError());
+
+  unsigned int last_rotation_state;
+  for (unsigned int cell_no=rotation_states.size()-1; cell_no>=track_back_to_cell; --cell_no)
+    {
+                                      // store rotation state of
+                                      // topmost cell, as we will
+                                      // have to advance that by one
+      last_rotation_state = rotation_states.top();
+      
+                                      // first mark faces of that
+                                      // cell as no more used
+      cells[cell_no].mark_faces_unused (last_rotation_state);
+
+                                      // then pop state from
+                                      // stack
+      rotation_states.pop();
+    };
+
+                                  // now we will have to find out
+                                  // whether we can try the last cell
+                                  // we have popped from the stack in
+                                  // another rotation state, or will
+                                  // have to backtrack further:
+  if (last_rotation_state < rotational_states_of_cells-1)
+    {
+                                      // possible. push that state to
+                                      // the stack and leave
+      rotation_states.push (last_rotation_state+1);
+      return;
+    }
+  else
+    {
+                                      // last cell can't be rotated
+                                      // further. go on with
+                                      // backtracking
+      const typename vector<Cell>::iterator
+       try_cell = cells.begin() + rotation_states.size();
+//        cout << "Further backtracking from " << rotation_states.size();
+//        cout << ". Neighbors are ";
+//        for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+//     cout << (int)try_cell->neighbors[i] << ' ';
+//        cout << "Will track back to cell " << try_cell->track_back_to_cell << endl;
+
+      track_back_to_cell = try_cell->track_back_to_cell;
+      Assert (track_back_to_cell > 0, ExcInternalError());
+
+                                      // track further back. this
+                                      // could be done by recursive
+                                      // calls of this function,
+                                      // which in this case would
+                                      // represent a tail-recursion
+                                      // as there is nothing more to
+                                      // be done after calling the
+                                      // function recursively, but we
+                                      // prefer to write down the
+                                      // tail-recursion by hand using
+                                      // a goto, since the compiler
+                                      // seems to have problems to
+                                      // rewrite the tail recursion
+                                      // as a goto.
+      goto top_of_function;
+    };
+};
+
+
+
+template <int dim>
+void GridReordering<dim>::find_reordering (vector<Cell<dim> >         &cells,
+                                          vector<CellData<dim> >     &original_cells,
+                                          const vector<unsigned int> &new_cell_numbers)
+{
+  const unsigned int n_cells = cells.size();
+  
+                                  // stack of value indicating that
+                                  // the nth cell needs to be rotated
+                                  // so-and-so often, where n is the
+                                  // position on the stack
+  stack<unsigned int, vector<unsigned int> > rotation_states;
+
+                                  // for the first cell, the
+                                  // rotational state can never be
+                                  // important, since we can rotate
+                                  // all other cells
+                                  // accordingly. therefore preset
+                                  // the rotation state of the first
+                                  // cell
+  rotation_states.push (0);
+  cells[0].mark_faces_used (rotation_states.top());
+  
+  while (true)
+    {
+                                      // if all cells have a coherent
+                                      // orientation, then we can
+                                      // exit the main loop
+      if (rotation_states.size() == n_cells)
+       break;
+      
+                                      // try to push back another
+                                      // cell in orientation zero
+      rotation_states.push (0);
+
+                                      // check whether the present
+                                      // cell in the present
+                                      // orientation is valid
+      check_topmost_cell:
+
+      static unsigned int max_size = 0;
+      if (rotation_states.size() > max_size)
+       {
+         if (max_size % 10 == 0)
+           cout << "New max size " << rotation_states.size() << endl;
+         max_size = rotation_states.size();
+       };
+      
+      const typename vector<Cell>::iterator
+       try_cell = cells.begin() + rotation_states.size()-1;
+      if (try_cell->check_consistency (rotation_states.top()))
+       {
+                                          // yes, works, we found a
+                                          // way of how to add the
+                                          // present cell to the
+                                          // existing cells without
+                                          // violating any ordering
+                                          // constraints. now mark
+                                          // the respective faces as
+                                          // used and go on with the
+                                          // next cell
+         try_cell->mark_faces_used (rotation_states.top());
+         
+//       cout << "Added cell " << try_cell->cell_no
+//            << " in rotation " << rotation_states.top() << endl;
+                                          // go on with next cell
+         continue;
+       }
+      else
+       {
+                                          // no, doesn't work. see if
+                                          // we can rotate the top
+                                          // cell so that it works
+         if (rotation_states.top()+1 < rotational_states_of_cells)
+           {
+                                              // yes, can be
+                                              // done. then do so and
+                                              // check again
+             ++rotation_states.top();
+             goto check_topmost_cell;
+           }
+         else
+           {
+                                              // no, no more
+                                              // orientation of the
+                                              // top cell possible,
+                                              // we have to backtrack
+                                              // some way
+//           cout << "Failure with cell " << rotation_states.size()-1;
+//           cout << ". Neighbors are ";
+//           for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+//             cout << (int)try_cell->neighbors[i] << ' ';
+//           cout << "Will track back to cell " << try_cell->track_back_to_cell << endl;
+
+                                              // first pop rotational
+                                              // state of top cell,
+                                              // since for that no
+                                              // faces have been
+                                              // marked as used yet
+             rotation_states.pop();
+                                              // then track back
+             track_back (cells, rotation_states, try_cell->track_back_to_cell);
+                                              // and go on by
+                                              // checking the now
+                                              // topmost cell
+             goto check_topmost_cell;
+           };
+       };
+    };
+
+
+
+                                  // rotate the cells according to
+                                  // the results we have found. since
+                                  // we operate on a stack, we do the
+                                  // rotations from the back of the
+                                  // array to the front
+  while (rotation_states.size() != 0)
+    {
+      const unsigned int
+       new_cell_number = rotation_states.size()-1;
+      const unsigned int
+       old_cell_number = find (new_cell_numbers.begin(),
+                               new_cell_numbers.end(),
+                               new_cell_number) - new_cell_numbers.begin();
+      Assert (old_cell_number < cells.size(), ExcInternalError());
+
+      original_cells[old_cell_number].rotate (rotation_states.top());
+      rotation_states.pop ();
+    };
+};
+
+
+
+template <int dim>
+vector<unsigned int>
+GridReordering<dim>::presort_cells (vector<Cell<dim> > &cells,
+                                   map<Face,FaceData> &faces)
+{
+                                  // first find the cell with the
+                                  // least neighbors
+  unsigned int min_neighbors           = cells[0].count_neighbors();
+  unsigned int cell_with_min_neighbors = 0;
+  for (unsigned int i=1; i<cells.size(); ++i)
+    if (min_neighbors > cells[i].count_neighbors())
+      {  
+       min_neighbors = cells[i].count_neighbors();
+       cell_with_min_neighbors = i;
+       if (min_neighbors == 1)
+                                          // better is not possible
+         break;
+      };
+
+                                  // have an array into which we
+                                  // insert the new cells numbers of
+                                  // each cell
+  const unsigned int invalid_cell_number = static_cast<unsigned int>(-1);
+  vector<unsigned int> new_cell_numbers (cells.size(), invalid_cell_number);
+
+                                  // and have an array of the next
+                                  // cells to be numbered (old numbers)
+  vector<unsigned int> next_round_cells (1, cell_with_min_neighbors);
+
+  unsigned int next_free_new_number = 0;
+  
+                                  // while there are still cells to
+                                  // be renumbered:
+  while (next_round_cells.size() != 0)
+    {
+      for (unsigned int i=0; i<next_round_cells.size(); ++i)
+       {
+         Assert (new_cell_numbers[next_round_cells[i]] == invalid_cell_number,
+                 ExcInternalError());
+         
+         new_cell_numbers[next_round_cells[i]] = next_free_new_number;
+         ++next_free_new_number;
+       };
+
+                                      // for the next round, find all
+                                      // neighbors of the cells of
+                                      // this round which have not
+                                      // yet been renumbered
+      vector<unsigned int> new_next_round_cells;
+      for (unsigned int i=0; i<next_round_cells.size(); ++i)
+       for (unsigned int n=0; n<GeometryInfo<dim>::faces_per_cell; ++n)
+         if (cells[next_round_cells[i]].neighbors[n] != Cell<dim>::invalid_neighbor)
+           if (new_cell_numbers[cells[next_round_cells[i]].neighbors[n]]
+               == invalid_cell_number)
+             new_next_round_cells.push_back (cells[next_round_cells[i]].neighbors[n]);
+
+           
+      
+                                      // if no more cells have been
+                                      // found, then we must have
+                                      // renumbered all cells already
+      if (new_next_round_cells.size() == 0)
+       Assert (next_free_new_number == cells.size(), ExcInternalError());
+
+                                      // eliminate duplicates from
+                                      // the new_next_round_cells
+                                      // array. note that a cell
+                                      // which is entered into this
+                                      // array might have been
+                                      // entered more than once since
+                                      // it might be a neighbor of
+                                      // more than one cell of the
+                                      // present round
+                                      //
+                                      // in order to eliminate
+                                      // duplicates, we first sort
+                                      // tha array and then copy over
+                                      // only unique elements to the
+                                      // next_round_cells array,
+                                      // which is needed for the next
+                                      // loop iteration anyway
+      sort (new_next_round_cells.begin(), new_next_round_cells.end());
+      next_round_cells.clear ();
+      unique_copy (new_next_round_cells.begin(), new_next_round_cells.end(),
+                  back_inserter(next_round_cells));
+    };
+  Assert (find (new_cell_numbers.begin(), new_cell_numbers.end(), invalid_cell_number)
+         ==
+         new_cell_numbers.end(),
+         ExcInternalError());
+
+                                  // now that we know in which order
+                                  // to sort the cells, do so:
+  vector<Cell<dim> > new_cells (cells.size());
+  for (unsigned int i=0; i<cells.size(); ++i)
+    new_cells[new_cell_numbers[i]] = cells[i];
+                                  // then switch old and new array
+  swap (cells, new_cells);
+  
+                                  // now we still have to convert all
+                                  // old cell numbers to new cells
+                                  // numbers
+  for (unsigned int c=0; c<cells.size(); ++c)
+    {
+      cells[c].cell_no = new_cell_numbers[cells[c].cell_no];
+      Assert (cells[c].cell_no == c, ExcInternalError());
+
+      for (unsigned int n=0; n<GeometryInfo<dim>::faces_per_cell; ++n)
+       cells[c].neighbors[n] = new_cell_numbers[cells[c].neighbors[n]];
+    };
+
+  for (typename map<Face,FaceData>::iterator i=faces.begin(); i!=faces.end(); ++i)
+    for (unsigned int k=0; k<2; ++k)
+      if (i->second.adjacent_cells[k] != FaceData::invalid_adjacent_cell)
+       i->second.adjacent_cells[k] = new_cell_numbers[i->second.adjacent_cells[k]];
+
+  return new_cell_numbers;
+};
+
+                     
+
+template <int dim>
+void GridReordering<dim>::reorder_cells (vector<CellData<dim> > &original_cells)
+{
+                                  // the present function might
+                                  // actually work in 3d, but the
+                                  // ones it calls probably not. will
+                                  // have to think about what needs
+                                  // to be changed. however, rather
+                                  // than killing the program when
+                                  // calling functions that may not
+                                  // work, kill it here if the
+                                  // dimension is not appropriate
+  Assert (dim==2, ExcNotImplemented());
+  
+                                  // we need more information than
+                                  // provided by the input parameter,
+                                  // in particular we need
+                                  // neighborship relations between
+                                  // cells. therefore copy over the
+                                  // old cells to another class that
+                                  // provides space to these
+                                  // informations
+  vector<Cell<dim> > cells;
+  cells.reserve (original_cells.size());
+  for (unsigned int i=0; i<original_cells.size(); ++i)
+    cells.push_back (Cell<dim>(original_cells[i], i));
+  
+                                  // first generate all the faces
+                                  // possible, i.e. in each possible
+                                  // direction and rotational state
+  map<Face,FaceData> faces;
+  for (unsigned int cell_no=0; cell_no<cells.size(); ++cell_no)
+    cells[cell_no].insert_faces (faces);
+
+                                  // after all faces have been filled
+                                  // and the faces have indices of
+                                  // their neighbors, we may also
+                                  // insert the neighbor indices into
+                                  // the cells themselves
+  for (unsigned int cell_no=0; cell_no<cells.size(); ++cell_no)
+    cells[cell_no].fix_cell_neighbors ();
+
+
+                                  // do a preordering step in order
+                                  // to make further backtracking
+                                  // more local
+  const vector<unsigned int>
+    new_cell_numbers = presort_cells (cells, faces);
+
+                                  // finally do some preliminary work
+                                  // to make backtracking simpler
+                                  // later
+  for (unsigned int cell_no=0; cell_no<cells.size(); ++cell_no)
+    cells[cell_no].find_backtracking_point ();
+  
+                                  // now do the main work
+  find_reordering (cells, original_cells, new_cell_numbers);
+};
+
+
+
+template <>
+void GridReordering<1>::reorder_cells (vector<CellData<1> > &)
+{
+                                  // there should not be much to do
+                                  // in 1d...
+};
+
+
+
+
+// explicit instantiations. only require the main function, it should
+// then claim whatever templates it needs
+template
+void
+GridReordering<deal_II_dimension>::
+reorder_cells (vector<CellData<deal_II_dimension> > &);
+

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.