]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Check in first version of Mike's fabulous 3d reordering algorithm.
authormike <mike@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Sep 2003 22:11:56 +0000 (22:11 +0000)
committermike <mike@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Sep 2003 22:11:56 +0000 (22:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@8011 0785d39b-7218-0410-832d-ea1e28bc413d

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

index e839f9b72e0c7b9940d9de0777f93bdab2b05e66..7a71c0468b0380aabb4a6bcf95a1f916bf95cdd1 100644 (file)
 #include <base/config.h>
 #include <grid/tria.h>
 
-#include <map>
 #include <vector>
 
 
-namespace internal
-{
-/**
- * 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
- * 3d case.
- *
- * @author Wolfgang Bangerth, 2000
- */
-  template <>
-  class GridReorderingInfo<3>
-  {
-    public:
-                                      /**
-                                       * Number of possible
-                                       * orientations of a cell in 3d,
-                                       * i.e. in how many ways can we
-                                       * arrange the eight vertices of
-                                       * a cube such that the numbering
-                                       * conventions of deal.II are
-                                       * satisfied. Since we can put
-                                       * the cube on each of its six
-                                       * faces and rotate it into four
-                                       * ways on each of these faces,
-                                       * the total number is 24.
-                                       */
-      static const unsigned int rotational_states_of_cells = 24;
-
-                                      /**
-                                       * Number of possible
-                                       * orientations of a face in
-                                       * 2d. It is the quadrilateral
-                                       * rotated in all four
-                                       * orientations, and also viewed
-                                       * from the back, i.e. vertices
-                                       * not numbered in
-                                       * counter-clockwise sense but in
-                                       * the other direction.
-                                       */
-      static const unsigned int rotational_states_of_faces = 8;
-  };
-}
-
-
-
-
-
-namespace internal
-{
-/**
- * Implement the algorithm described in the documentation of the
- * GridReordering<2> class.
- *
- * @author Michael Anderson, 2003
- */
-  namespace GridReordering2d
-  {
-
-                                     /**
-                                      * Check whether a given
-                                      * arrangement of cells is
-                                      * already consisten. If this is
-                                      * the case, then we skip the
-                                      * reordering pass.
-                                      *
-                                      * This function works by looping
-                                      * over all cells, checking
-                                      * whether one of its faces
-                                      * already exists in a list of
-                                      * edges, and if it already
-                                      * exists in reverse order, then
-                                      * return @p{false}. If it is not
-                                      * already in the list, or in the
-                                      * correct direction, then go on
-                                      * with the next faces or cell.
-                                      */
-    bool
-    is_consistent (const std::vector<CellData<2> > &cells);
-    
-    
-/**
- * Defines a variety of variables related to the connectivity of a
- * simple quad element. This includes the nodes on each edge, which
- * edges come into each node and what the default deal.II directions
- * are for the quad.
- *
- * @begin{verbatim}
- *       s2
- *       
- *     +-->--+       
- *     |3   2|     
- * s3  ^     ^ s1   
- *     |0   1|     
- *     +-->--+               
- *   
- *       s0           
- * @end{verbatim}
- *
- * @author Michael Anderson, 2003
- */ 
-    class ConnectGlobals
-    {
-      public:
-                                        /**    
-                                         * The nodes on each edge in
-                                         * anti-clockwise order
-                                         * { {0,1},{1,2},{2,3},{3,0} }
-                                         */
-       static const int EdgeToNode[4][2];
-
-                                        /**
-                                         * The edges comin into each
-                                         * node, in anti-clockwise
-                                         * order
-                                         * { {3,0},{0,1},{1,2},{2,3} }
-                                         */
-       static const int NodeToEdge[4][2];
-  
-                                        /**
-                                         * The nodes on each edge in
-                                         * "default direction order".
-                                         * {{0,1},{1,2},{3,2},{0,3}}
-                                         */
-       static const int DefaultOrientation[4][2];
-    };
-
-
-/**
- * An enriched quad with information about how the mesh fits together
- * so that we can move around the mesh efficiently.
- *
- * @author Michael Anderson, 2003
- */
-    class MQuad
-    {
-      public:
-                                        /**
-                                         * v0 - v3 are indexes of the vertices of the quad,
-                                         * s0 - s3   are indexes for the sides of the quad
-                                         */
-       MQuad (const unsigned int  v0,
-              const unsigned int  v1,
-              const unsigned int  v2,
-              const unsigned int  v3,
-              const unsigned int  s0,
-              const unsigned int  s1,
-              const unsigned int  s2,
-              const unsigned int  s3,
-              const CellData<2>  &cd);
-       
-                                        /**
-                                         * Stores the vertex numbers
-                                         */
-       unsigned int v[4];
-                                        /**
-                                         * Stores the side numbers
-                                         */
-       unsigned int side[4]; 
-
-                                        /**
-                                         * Copy of the @p{CellData} object
-                                         * from which we construct the
-                                         * data of this object.
-                                         */
-       CellData<2>  original_cell_data;
-       
-                                        /**
-                                         * Makes an MQuad from the
-                                         * given CellData and MSide
-                                         * list.  Is derived from
-                                         * binary_function to be
-                                         * usable with STL
-                                         * containers.
-                                         *
-                                         * Also assumes that the
-                                         * edges listed present in
-                                         * the CellData are already
-                                         * present in the elist
-                                         * vector.
-                                         */ 
-       struct MakeQuad;
-    };
-
-/**
- * The enriched side class containing connectivity information.
- * Orientation is from v0 to v1; Initially this should have v0<v1.
- * After global orientation could be either way.
- *
- * @author Michael Anderson, 2003
- */
-    struct MSide
-    {
-                                        /**
-                                         * Constructor.
-                                         */
-       MSide (const unsigned int initv0,
-              const unsigned int initv1);
-
-                                        /**
-                                         * Return whether the sides
-                                         * are equal, even if their
-                                         * ends are reversed.
-                                         */
-       bool operator==(const MSide& s2) const;
-
-                                        /**
-                                         * Return the opposite.
-                                         */
-       bool operator!=(const MSide& s2) const;
-       
-       unsigned int v0;
-       unsigned int v1;
-       unsigned int Q0;
-       unsigned int Q1;
-
-                                        /**
-                                         * Local side numbers on quads 0 and 1.
-                                         */
-       unsigned int lsn0, lsn1;
-       bool Oriented;
-
-                                        /**
-                                         * This class makes a MSide have v0<v1
-                                         */
-       struct SideRectify;
-
-                                        /**
-                                         * Provides a side ordering,
-                                         * s1<s2, without assuming
-                                         * v0<v1 in either of the
-                                         * sides.
-                                         */
-       struct SideSortLess;
-    };
-
-
-
-/**
- * Implement the 2d algorithm for grid reordering described in the
- * documentation of the @ref{GridReordering} class.
- *
- * @author Michael Anderson, 2003
- */   
-    class GridReordering
-    {
-      public:
-
-                                        /**
-                                         * Do the work intended by
-                                         * this class.
-                                         */
-       void reorient(std::vector<CellData<2> > &quads);
-      private:
-       
-                                        /** 
-                                         * Sets up the internal data
-                                         * structures so that the we can
-                                         * do side hopping and face
-                                         * switching efficiently. This
-                                         * means we need a whole bunch of
-                                         * connectivity information
-                                         */
-       void build_graph (const std::vector<CellData<2> > &inquads);
-
-                                        /** 
-                                         * Orient the internal data
-                                         * into deal.II format The
-                                         * orientation algorith is as
-                                         * follows
-                                         *
-                                         * 1) Find an unoriented quad (A)
-                                         *
-                                         * 2) Orient an un_oriented side (s) of (A)
-                                         *
-                                         * 3) side hop on (s) of (A) to get (B)
-                                         *
-                                         * 4) if opposite side to (s)
-                                         * of (B) is unoriented
-                                         * orient it
-                                         *
-                                         * 5) repeat 3) and 4) until
-                                         * side-hoppong fails (we've
-                                         * reached a boundary) or (s)
-                                         * has already been oriented
-                                         * (we've closed a loop or
-                                         * unoriented sides).
-                                         *
-                                         * 6) Repeat 2), 3) ,4) and
-                                         * 5) on other unoriented
-                                         * sides of (A)
-                                         *
-                                         * 7) Choose a new unoriented
-                                         * A.
-                                         */
-       void orient();
-  
-                                        /**
-                                         * Get the (now correctly
-                                         * oriented if we've called
-                                         * orient) quads.
-                                         */
-       void get_quads(std::vector<CellData<2> > &outquads) const;
-
-                                        /**
-                                         * Orient_side(qnum,lsn)
-                                         * orients the local side lsn
-                                         * of the quad qnum in the
-                                         * triangulation. If the side
-                                         * opposite lsn is oriented
-                                         * then lsn is oriented to
-                                         * match it. Otherwise it is
-                                         * oriented in the "default"
-                                         * direction for the quad.
-                                         */
-       void orient_side (const unsigned int quadnum,
-                         const unsigned int localsidenum);
-                                        /**
-                                         * Returns true if all sides
-                                         * of the quad quadnum are
-                                         * oriented.
-                                         */
-       bool is_fully_oriented_quad (const unsigned int quadnum) const;
-  
-                                        /**
-                                         * Returns true if the side lsn
-                                         * of the quad quadnum is
-                                         * oriented.
-                                         */
-       bool is_oriented_side (const unsigned int quadnum,
-                              const unsigned int lsn) const;
-
-                                        /**
-                                         * Returns true is the side is
-                                         * oriented in the "default"
-                                         * direction
-                                         */
-       bool is_side_default_oriented (const unsigned int qnum,
-                                      const unsigned int lsn) const;
-  
-                                        /**
-                                         * Increases UnOrQLoc from
-                                         * it's original value to the
-                                         * next quad with an
-                                         * unoriented side. Returns
-                                         * true if there was another
-                                         * unoriented quad.
-                                         */
-       bool get_unoriented_quad (unsigned int &UnOrQLoc) const;
-
-                                        /**
-                                         * Sets sidenum to the local
-                                         * sidenumber of an
-                                         * unoriented side of the
-                                         * quad quadnum. Returns true
-                                         * if such a side exists.
-                                         */
-       bool get_unoriented_side (const unsigned int quadnum,
-                                 unsigned int &sidenum) const;
-  
-                                        /**
-                                         * side_hop(&qnum, &lsn) has
-                                         * qnum being the quadnumber
-                                         * of a quad in the
-                                         * triangulation, and a local
-                                         * side number. side_hop then
-                                         * sets qnum to the
-                                         * quadnumber across the
-                                         * other side of the side,
-                                         * and sets lsn so that
-                                         * quads[qnum].sides[lsn] is
-                                         * the same before and after
-                                         * the call.  if there is no
-                                         * other quad on the other
-                                         * side of the current quad,
-                                         * then side_hop returns
-                                         * false.
-                                         */
-       bool side_hop (unsigned int &qnum,
-                      unsigned int &lsn) const;
-  
-                                        /**
-                                         * Sets lsn so that it points
-                                         * to the opposite side of
-                                         * the current quad (qnum)
-                                         * that it was originally
-                                         * pointing to.
-                                         */
-       bool switch_faces (unsigned int &qnum,
-                          unsigned int &lsn) const;
-
-                                        /**
-                                         * A list of enriched
-                                         * sides/edges of the mesh.
-                                         */
-       std::vector<MSide> sides;
-                                        /**
-                                         * A list of enriched quads
-                                         * in the mesh.
-                                         */
-       std::vector<MQuad> mquads;
-    };
-  }  // namespace GridReordering2d
-}  // namespace internal
-
-
 
 
 /**
@@ -1006,425 +586,6 @@ namespace internal
 template <int dim>
 class GridReordering
 {
-  public:
-                                    /**
-                                     * This is the main function,
-                                     * doing what is announced in the
-                                     * general documentation of this
-                                     * class.
-                                     */
-    static void reorder_cells (std::vector<CellData<dim> > &original_cells);
-    
-  private:
-
-                                    /**
-                                     * Forward declarations of local
-                                     * classes.
-                                     */
-    class Cell;
-    class Face;
-    class FaceData;
-
-                                    /**
-                                     * Typedef for a stack type that
-                                     * describes the rotational
-                                     * states of all cells that have
-                                     * already been fitted into the
-                                     * grid. It is mostly used like a
-                                     * stack, but sometimes we need
-                                     * random access into values
-                                     * below the top, so we can't use
-                                     * the @p{stack} adaptor from
-                                     * STL.
-                                     */
-    typedef std::vector<unsigned int> RotationStack;
-    
-                                    /**
-                                     * 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
-                                     */
-    class Cell : public CellData<dim>
-    {
-      public:
-                                        /**
-                                         * 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.
-                                         *
-                                         * Note that the order in
-                                         * which the faces for a
-                                         * specific rotational state
-                                         * appear is not specified,
-                                         * as this is not
-                                         * important. It is only
-                                         * important that each face
-                                         * in one rotational state or
-                                         * other appears once for
-                                         * each orientation of the
-                                         * cell.
-                                         */
-       typename std::map<Face,FaceData>::iterator
-       faces[internal::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);
-
-                                        /**
-                                         * Copy constructor to work
-                                         * around a bug in a gcc3.0
-                                         * snapshot.
-                                         */
-       Cell (const Cell &);
-       
-                                        /**
-                                         * 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 (std::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. The size of the array 
-                                         * is equal to
-                                         * @p{GeometryInfo<dim>::vertices_per_face}.
-                                         */
-       unsigned int vertices[1<<(dim-1)];
-
-                                        /**
-                                         * Comparison operator. Use
-                                         * the vertex indices as
-                                         * primary, secondary,
-                                         * ... criteria for
-                                         * comparison.
-                                         */
-       bool operator < (const Face &face) const;
-
-                                        /**
-                                         * Check for equality of
-                                         * vertex indices.
-                                         */
-       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.
-                                         *
-                                         * Note that the order in
-                                         * which the reverted faces
-                                         * appear here is not
-                                         * specified, as it is not
-                                         * important for the
-                                         * algorithm.
-                                         */
-       typename std::map<Face,FaceData >::const_iterator
-       reverse_faces[internal::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 (std::vector<Cell>  &cells,
-                           RotationStack               &rotation_states,
-                           const unsigned int          track_back_to_cell);
-
-    static bool try_rotate_single_neighbors (std::vector<Cell>  &cells,
-                                            RotationStack               &rotation_states);
-    
-                                    /**
-                                     * 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 (std::vector<Cell>           &cells,
-                                std::vector<CellData<dim> > &original_cells,
-                                const std::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
-    std::vector<unsigned int>
-    presort_cells (std::vector<Cell>       &cells,
-                  std::map<Face,FaceData> &faces);
-
-                                    /** 
-                                     * By the resolution of Defect
-                                     * Report 45 to the ISO C++ 1998
-                                     * standard, nested classes
-                                     * automatically have access to
-                                     * members of the enclosing
-                                     * class. Nevertheless, some
-                                     * compilers don't implement this
-                                     * resolution yet, so we have to
-                                     * make them @p{friend}, which
-                                     * doesn't hurt on the other
-                                     * compilers as well.
-                                     */
-    friend class Cell;
-    friend class Face;
-    friend class FaceData;
 };
 
 
@@ -1470,10 +631,28 @@ class GridReordering<2>
 
 
 
-/* -------------- declaration of explicit specializations ------------- */
-
+/**
+ * This specialization of the general template implements the
+ * 3d-algorithm described in the documentation of the general
+ * template.
+ *
+ * @author Michael Anderson, 2003
+ */
 template <>
-void GridReordering<3>::Cell::insert_faces (std::map<Face,FaceData> &global_faces);
+class GridReordering<3>
+{
+  public:
+                                    /**
+                                     *  This is the main function,
+                                     *  doing what is announced in
+                                     *  the general documentation of
+                                     *  this class.
+                                     */
+    static void reorder_cells (std::vector<CellData<3> > &original_cells,
+                              std::vector<Point<3> > * node_vec_ptr=NULL);
+};
+
+
 
 
 #endif
diff --git a/deal.II/deal.II/include/grid/grid_reordering_internal.h b/deal.II/deal.II/include/grid/grid_reordering_internal.h
new file mode 100644 (file)
index 0000000..bf4f190
--- /dev/null
@@ -0,0 +1,636 @@
+//----------------------------  grid_reordering_internal.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 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_internal.h  ---------------------------
+#ifndef __deal2__grid_reordering_internal_h
+#define __deal2__grid_reordering_internal_h
+
+
+#include <base/config.h>
+#include <grid/tria.h>
+
+#include <map>
+#include <vector>
+
+
+
+
+
+namespace internal
+{
+/**
+ * Implement the algorithm described in the documentation of the
+ * GridReordering<2> class.
+ *
+ * @author Michael Anderson, 2003
+ */
+  namespace GridReordering2d
+  {
+
+                                     /**
+                                      * Check whether a given
+                                      * arrangement of cells is
+                                      * already consisten. If this is
+                                      * the case, then we skip the
+                                      * reordering pass.
+                                      *
+                                      * This function works by looping
+                                      * over all cells, checking
+                                      * whether one of its faces
+                                      * already exists in a list of
+                                      * edges, and if it already
+                                      * exists in reverse order, then
+                                      * return @p{false}. If it is not
+                                      * already in the list, or in the
+                                      * correct direction, then go on
+                                      * with the next faces or cell.
+                                      */
+    bool
+    is_consistent (const std::vector<CellData<2> > &cells);
+    
+    
+/**
+ * Defines a variety of variables related to the connectivity of a
+ * simple quad element. This includes the nodes on each edge, which
+ * edges come into each node and what the default deal.II directions
+ * are for the quad.
+ *
+ * @begin{verbatim}
+ *       s2
+ *       
+ *     +-->--+       
+ *     |3   2|     
+ * s3  ^     ^ s1   
+ *     |0   1|     
+ *     +-->--+               
+ *   
+ *       s0           
+ * @end{verbatim}
+ *
+ * @author Michael Anderson, 2003
+ */ 
+    class ConnectGlobals
+    {
+      public:
+                                        /**    
+                                         * The nodes on each edge in
+                                         * anti-clockwise order
+                                         * { {0,1},{1,2},{2,3},{3,0} }
+                                         */
+       static const int EdgeToNode[4][2];
+
+                                        /**
+                                         * The edges comin into each
+                                         * node, in anti-clockwise
+                                         * order
+                                         * { {3,0},{0,1},{1,2},{2,3} }
+                                         */
+       static const int NodeToEdge[4][2];
+  
+                                        /**
+                                         * The nodes on each edge in
+                                         * "default direction order".
+                                         * {{0,1},{1,2},{3,2},{0,3}}
+                                         */
+       static const int DefaultOrientation[4][2];
+    };
+
+
+/**
+ * An enriched quad with information about how the mesh fits together
+ * so that we can move around the mesh efficiently.
+ *
+ * @author Michael Anderson, 2003
+ */
+    class MQuad
+    {
+      public:
+                                        /**
+                                         * v0 - v3 are indexes of the vertices of the quad,
+                                         * s0 - s3   are indexes for the sides of the quad
+                                         */
+       MQuad (const unsigned int  v0,
+              const unsigned int  v1,
+              const unsigned int  v2,
+              const unsigned int  v3,
+              const unsigned int  s0,
+              const unsigned int  s1,
+              const unsigned int  s2,
+              const unsigned int  s3,
+              const CellData<2>  &cd);
+       
+                                        /**
+                                         * Stores the vertex numbers
+                                         */
+       unsigned int v[4];
+                                        /**
+                                         * Stores the side numbers
+                                         */
+       unsigned int side[4]; 
+
+                                        /**
+                                         * Copy of the @p{CellData} object
+                                         * from which we construct the
+                                         * data of this object.
+                                         */
+       CellData<2>  original_cell_data;
+       
+                                        /**
+                                         * Makes an MQuad from the
+                                         * given CellData and MSide
+                                         * list.  Is derived from
+                                         * binary_function to be
+                                         * usable with STL
+                                         * containers.
+                                         *
+                                         * Also assumes that the
+                                         * edges listed present in
+                                         * the CellData are already
+                                         * present in the elist
+                                         * vector.
+                                         */ 
+       struct MakeQuad;
+    };
+
+/**
+ * The enriched side class containing connectivity information.
+ * Orientation is from v0 to v1; Initially this should have v0<v1.
+ * After global orientation could be either way.
+ *
+ * @author Michael Anderson, 2003
+ */
+    struct MSide
+    {
+                                        /**
+                                         * Constructor.
+                                         */
+       MSide (const unsigned int initv0,
+              const unsigned int initv1);
+
+                                        /**
+                                         * Return whether the sides
+                                         * are equal, even if their
+                                         * ends are reversed.
+                                         */
+       bool operator==(const MSide& s2) const;
+
+                                        /**
+                                         * Return the opposite.
+                                         */
+       bool operator!=(const MSide& s2) const;
+       
+       unsigned int v0;
+       unsigned int v1;
+       unsigned int Q0;
+       unsigned int Q1;
+
+                                        /**
+                                         * Local side numbers on quads 0 and 1.
+                                         */
+       unsigned int lsn0, lsn1;
+       bool Oriented;
+
+                                        /**
+                                         * This class makes a MSide have v0<v1
+                                         */
+       struct SideRectify;
+
+                                        /**
+                                         * Provides a side ordering,
+                                         * s1<s2, without assuming
+                                         * v0<v1 in either of the
+                                         * sides.
+                                         */
+       struct SideSortLess;
+    };
+
+
+
+/**
+ * Implement the 2d algorithm for grid reordering described in the
+ * documentation of the @ref{GridReordering} class.
+ *
+ * @author Michael Anderson, 2003
+ */   
+    class GridReordering
+    {
+      public:
+
+                                        /**
+                                         * Do the work intended by
+                                         * this class.
+                                         */
+       void reorient(std::vector<CellData<2> > &quads);
+      private:
+       
+                                        /** 
+                                         * Sets up the internal data
+                                         * structures so that the we can
+                                         * do side hopping and face
+                                         * switching efficiently. This
+                                         * means we need a whole bunch of
+                                         * connectivity information
+                                         */
+       void build_graph (const std::vector<CellData<2> > &inquads);
+
+                                        /** 
+                                         * Orient the internal data
+                                         * into deal.II format The
+                                         * orientation algorith is as
+                                         * follows
+                                         *
+                                         * 1) Find an unoriented quad (A)
+                                         *
+                                         * 2) Orient an un_oriented side (s) of (A)
+                                         *
+                                         * 3) side hop on (s) of (A) to get (B)
+                                         *
+                                         * 4) if opposite side to (s)
+                                         * of (B) is unoriented
+                                         * orient it
+                                         *
+                                         * 5) repeat 3) and 4) until
+                                         * side-hoppong fails (we've
+                                         * reached a boundary) or (s)
+                                         * has already been oriented
+                                         * (we've closed a loop or
+                                         * unoriented sides).
+                                         *
+                                         * 6) Repeat 2), 3) ,4) and
+                                         * 5) on other unoriented
+                                         * sides of (A)
+                                         *
+                                         * 7) Choose a new unoriented
+                                         * A.
+                                         */
+       void orient();
+  
+                                        /**
+                                         * Get the (now correctly
+                                         * oriented if we've called
+                                         * orient) quads.
+                                         */
+       void get_quads(std::vector<CellData<2> > &outquads) const;
+
+                                        /**
+                                         * Orient_side(qnum,lsn)
+                                         * orients the local side lsn
+                                         * of the quad qnum in the
+                                         * triangulation. If the side
+                                         * opposite lsn is oriented
+                                         * then lsn is oriented to
+                                         * match it. Otherwise it is
+                                         * oriented in the "default"
+                                         * direction for the quad.
+                                         */
+       void orient_side (const unsigned int quadnum,
+                         const unsigned int localsidenum);
+                                        /**
+                                         * Returns true if all sides
+                                         * of the quad quadnum are
+                                         * oriented.
+                                         */
+       bool is_fully_oriented_quad (const unsigned int quadnum) const;
+  
+                                        /**
+                                         * Returns true if the side lsn
+                                         * of the quad quadnum is
+                                         * oriented.
+                                         */
+       bool is_oriented_side (const unsigned int quadnum,
+                              const unsigned int lsn) const;
+
+                                        /**
+                                         * Returns true is the side is
+                                         * oriented in the "default"
+                                         * direction
+                                         */
+       bool is_side_default_oriented (const unsigned int qnum,
+                                      const unsigned int lsn) const;
+  
+                                        /**
+                                         * Increases UnOrQLoc from
+                                         * it's original value to the
+                                         * next quad with an
+                                         * unoriented side. Returns
+                                         * true if there was another
+                                         * unoriented quad.
+                                         */
+       bool get_unoriented_quad (unsigned int &UnOrQLoc) const;
+
+                                        /**
+                                         * Sets sidenum to the local
+                                         * sidenumber of an
+                                         * unoriented side of the
+                                         * quad quadnum. Returns true
+                                         * if such a side exists.
+                                         */
+       bool get_unoriented_side (const unsigned int quadnum,
+                                 unsigned int &sidenum) const;
+  
+                                        /**
+                                         * side_hop(&qnum, &lsn) has
+                                         * qnum being the quadnumber
+                                         * of a quad in the
+                                         * triangulation, and a local
+                                         * side number. side_hop then
+                                         * sets qnum to the
+                                         * quadnumber across the
+                                         * other side of the side,
+                                         * and sets lsn so that
+                                         * quads[qnum].sides[lsn] is
+                                         * the same before and after
+                                         * the call.  if there is no
+                                         * other quad on the other
+                                         * side of the current quad,
+                                         * then side_hop returns
+                                         * false.
+                                         */
+       bool side_hop (unsigned int &qnum,
+                      unsigned int &lsn) const;
+  
+                                        /**
+                                         * Sets lsn so that it points
+                                         * to the opposite side of
+                                         * the current quad (qnum)
+                                         * that it was originally
+                                         * pointing to.
+                                         */
+       bool switch_faces (unsigned int &qnum,
+                          unsigned int &lsn) const;
+
+                                        /**
+                                         * A list of enriched
+                                         * sides/edges of the mesh.
+                                         */
+       std::vector<MSide> sides;
+                                        /**
+                                         * A list of enriched quads
+                                         * in the mesh.
+                                         */
+       std::vector<MQuad> mquads;
+    };
+  }  // namespace GridReordering2d
+
+
+/**
+ * Implement the algorithm described in the documentation of the
+ * GridReordering<2> class.
+ *
+ * @author Michael Anderson, 2003
+ */
+  namespace GridReordering3d
+  {
+
+  /**
+   * During building the conectivity information we 
+   * dont need all the heavy duty information about 
+   * edges that we will need later. So we can save 
+   * memory and time by using these light-weight edges.
+   **/
+  class CheapEdge
+  {
+    public:
+      //! The first node
+      int node0;
+      //! The second node
+      int node1;
+      //! A simple constructor
+      CheapEdge(int n0, int n1);
+      //! Need a partial ordering for the STL
+      bool operator<(const CheapEdge & e2) const;
+  };
+
+  
+  class ElementInfo
+    {
+  public:
+    // The numbers of the edges coming into node i
+    // are given by edge_to_node[i][k] where k=0,1,2
+    int edge_to_node[8][3];
+
+    //the orientation of edge coming into node i
+    //is given by edge_to_node_orient[i][k] where k=0,1,2
+    //  1 means the given node is the start of the edge
+    // -1 means the end of the edge
+    int edge_to_node_orient[8][3];
+
+    // nodesonedge[i][0] is the start node for edge i
+    // nodesonedge[i][1] is the end node for edge i
+    int nodes_on_edge[12][2];
+    int nodes_on_face[6][4];
+};
+
+class DealElemInfo : public ElementInfo
+{
+  public:
+    DealElemInfo();
+};
+
+
+
+  //! A conectivity and orientation aware edge class
+  class Edge
+  {
+    public:
+      //! Simple constructor
+      Edge(int n0, int n1, int orient=0): orientation_flag(orient), group(0),
+                                          num_neighbouring_cubes(0), neighbouring_cubes(NULL)
+      {nodes[0]=n0; nodes[1]=n1;};
+      
+      //! Simple Destructor
+      ~Edge();
+      
+      //! The IDs for the end nodes
+      
+      int nodes[2];
+      /** 
+       * Whether the edge has been oriented (0), 
+       * points from node 0 to node 1 (1), 
+       * or the reverse (-1)
+       **/
+      int orientation_flag;
+
+      /** 
+       * Used to determine which "sheet" of parallel edges the edge falls in
+       * when oriented. 0 means not yet decided.  
+       **/
+      int group;
+
+      unsigned int num_neighbouring_cubes;
+      unsigned int * neighbouring_cubes;
+  };
+
+  //! A conectivity an orientation aware cell
+  /**
+   * The connectivity of the cell is not contained within
+   * (This was for flexability in using deal's ordering of 
+   * edges or the XDA format etc) For this information we 
+   * need the ElemInfo class. 
+   *
+   * One thing we do know is that the first four edges in 
+   * the edge class are parallel, as are the second four, 
+   * and the third four.
+   *
+   * \todo TODO: Need to move conectivity information out 
+   *             of cell and into edge.
+   **/
+  class Cell
+      {
+       public:
+         int edges[12]; //!< The IDs for each of the edges
+         int nodes[8];  //!< The IDs for each of the nodes
+         //! Which way do the edges point.
+         /**
+          * Whether node 0 of the edge is the base of 
+          * the edge in local element (1) or node 1 is 
+          * the base (-1)
+          **/
+         int local_orientation_flags[12]; 
+         /**
+          * An internal flag used to determine whether the cell is in the
+          * queue of cells to be orriented in the current sheet. 
+          **/
+         bool waiting_to_be_processed;  
+
+         //! Copy Constructor
+         Cell(const Cell& c)
+         {
+           for(int i=0;i<12;++i)
+           {
+             edges[i]=c.edges[i];
+             local_orientation_flags[i]=c.local_orientation_flags[i];
+           }
+           for(int i=0;i<8;++i)
+           {
+             nodes[i]=c.nodes[i];
+           }
+           waiting_to_be_processed=c.waiting_to_be_processed;
+         }
+
+         //! Default Constructor
+         Cell()
+         {
+           for(int i=0;i<12;++i)
+           {
+             edges[i]=-1;
+             local_orientation_flags[i]=1;
+           }
+           for(int i=0;i<8;++i)
+           {
+             nodes[i]=-1;
+           }
+           waiting_to_be_processed=false;
+         }
+      };
+
+
+  //! This holds all the pieces for orientation together..
+  /**
+   * Contains lists of nodes, edges and cells. 
+   * As well as the information about how they all connect together.
+   **/
+  class Mesh
+  {
+    public:
+      //! Information about how a cell is built up from nodes and edges.
+      const ElementInfo & info;
+      
+      //! The list of nodes
+      std::vector< Point<3> > node_list;
+      //! The list of edges 
+      std::vector<Edge> edge_list;
+      //! The list of cells
+      std::vector<Cell> cell_list;
+
+      //! Checks whether every cell in the mesh is sensible
+      /**
+       * By calling sanity_check(cell_num) on every cell.
+       **/
+      bool sanity_check() const;
+      //! Checks that every node matches with its edges
+      /**
+       * By calling sanity_check(cell_num,node_num) for each node
+       **/
+      bool sanity_check(int cell_num) const;
+      //! Checks that each edge going into a node is correctly setup
+      bool sanity_check_node(int cell_num, int i) const;
+
+      //! Default Constructor
+      Mesh(const ElementInfo & INFO): info(INFO) { }
+
+      //! Prints all information about the mesh
+      void dump() const;
+      //! Prints all information about the cell
+      void dump_cell(const Cell &c) const;
+      //! Writes edge information to a file.
+      void dump_edges(char const * const fname) const;
+
+    private:
+      //! Unimplemented private copy constructor to disable it.
+      Mesh(const Mesh&);
+      //! Unimplemented private assignemnet operator to disable it.
+      Mesh& operator=(const Mesh&);
+  };
+
+
+    class Orienter
+  {
+    public:
+
+      Orienter(){};
+
+      //The cube we're looking at now.
+      unsigned int cur_posn; 
+      //We have fully oriented all cubes before this one.
+      unsigned int marker_cube;
+
+      std::vector<int> SheetToProcess;
+
+      int cur_edge_group;
+
+      bool edge_orient_array[12];
+      
+      bool orient_edges(Mesh &m);
+      void orient_cubes(Mesh &m);
+      
+      bool GetNextUnorientedCube(Mesh &m);
+      bool is_oriented(const Mesh &m, int cell_num);
+
+      bool OrientEdgesInCurrentCube(Mesh &m);
+      bool OrientEdgeSetInCurrentCube(Mesh &m, int edge_set);
+      bool OrientNextUnorientedEdge(Mesh &m);
+      bool Consistant(Mesh &m, int cell_num);
+
+
+      void GetAdjacentCubes(Mesh &m);
+      bool GetNextActiveCube(Mesh &m);
+
+      bool CheckCellEdgeGroupConsistancy(const Mesh &m, const Cell &c) const;
+      bool CheckCellEdgeGroupConsistancy(const Mesh &m, const Cell & c, int egrp) const;
+
+  };
+
+
+  //! Creates the connectivity information for the mesh m.
+  void build_mesh(Mesh &m);
+
+
+    
+    
+  }  // namespace GridReordering3d
+}  // namespace internal
+
+
+
+#endif
index 67a45fad26440a58cf2fb6119e4b91e6d28e9b8a..ea7bcf0db38526026e024a06eacc2b5001f9e263 100644 (file)
 //
 //----------------------------  grid_reordering.cc  ---------------------------
 
-#include <base/thread_management.h>
 #include <grid/grid_reordering.h>
+#include <grid/grid_reordering_internal.h>
 
-#include <set>
 #include <algorithm>
+#include <set>
+#include <iostream>
+#include <fstream>
 
 
 
-namespace internal
-{
-// static variables
-#if deal_II_dimension == 3
-  const unsigned int GridReorderingInfo<3>::rotational_states_of_cells;
-  const unsigned int GridReorderingInfo<3>::rotational_states_of_faces;
-#endif
-
-}
-
-template <int dim>
-const unsigned int GridReordering<dim>::Cell::invalid_neighbor;
-
-template <int dim>
-const unsigned int GridReordering<dim>::FaceData::invalid_adjacent_cell;
-
 
 
 #if deal_II_dimension == 1
@@ -49,1143 +35,6 @@ void GridReordering<1>::reorder_cells (const std::vector<CellData<1> > &)
 
 
 
-#if deal_II_dimension == 3
-
-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>
-GridReordering<dim>::Cell::Cell (const Cell &c) :
-               CellData<dim> (c),
-                cell_no(c.cell_no),
-                track_back_to_cell(c.track_back_to_cell)
-{
-  for (unsigned int i=0; i<internal::GridReorderingInfo<dim>::rotational_states_of_cells; ++i)
-    for (unsigned int j=0; j<GeometryInfo<dim>::faces_per_cell; ++j)
-      faces[i][j]=c.faces[i][j];
-
-  for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
-    neighbors[i]=c.neighbors[i];
-}
-
-
-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 (std::map<Face,FaceData> &/*global_faces*/)
-{
-  Assert (false, ExcNotImplemented());
-}
-
-
-
-template <>
-void
-GridReordering<3>::Cell::insert_faces (std::map<Face,FaceData> &global_faces)
-{
-  const unsigned int dim = 3;
-
-                                  // first generate for each of the 6
-                                  // faces of a cell in 3d the four
-                                  // possible orientations and
-                                  // cross-link them among each other
-                                  //
-                                  // do this generation step by first
-                                  // only inserting each face in
-                                  // standard orientation and then
-                                  // fill in the other ones by
-                                  // rotation of these faces
-                                  //
-                                  // note that we have the indices
-                                  // reversed here compared to the
-                                  // Cell class, for simplicity
-  const Face new_faces_tmp[GeometryInfo<dim>::faces_per_cell]
-    = { { { this->vertices[0], this->vertices[1],
-            this->vertices[2], this->vertices[3] } },
-       { { this->vertices[4], this->vertices[5],
-            this->vertices[6], this->vertices[7] } },
-       { { this->vertices[0], this->vertices[1],
-            this->vertices[5], this->vertices[4] } },
-       { { this->vertices[1], this->vertices[5],
-            this->vertices[6], this->vertices[2] } },
-       { { this->vertices[3], this->vertices[2],
-            this->vertices[6], this->vertices[7] } },
-       { { this->vertices[0], this->vertices[4],
-            this->vertices[7], this->vertices[3] } } };
-  Face new_faces[GeometryInfo<dim>::faces_per_cell][internal::GridReorderingInfo<dim>::rotational_states_of_faces]
-    = { { new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0],
-          new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] },
-        { new_faces_tmp[1], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0],
-         new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] },
-       { new_faces_tmp[2], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0],
-          new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] },
-        { new_faces_tmp[3], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0],
-          new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] },
-        { new_faces_tmp[4], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0],
-          new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] },
-        { new_faces_tmp[5], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0],
-          new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] }};
-
-                                  // first do the faces in their
-                                  // usual direction
-  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-    for (unsigned int rot=1; rot<internal::GridReorderingInfo<dim>::rotational_states_of_faces/2; ++rot)
-      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
-       new_faces[face][rot].vertices[v]
-         = new_faces[face][0].vertices[(v+rot) % GeometryInfo<dim>::vertices_per_face];
-                                  // then do everything as viewed
-                                  // from the back. this is simple,
-                                  // as we only have to revert
-                                  // indices 1 and 3
-  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-    for (unsigned int rot=internal::GridReorderingInfo<dim>::rotational_states_of_faces/2;
-        rot<internal::GridReorderingInfo<dim>::rotational_states_of_faces;
-        ++rot)
-      {
-       for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
-         new_faces[face][rot].vertices[v]
-           = new_faces[face]
-           [rot-internal::GridReorderingInfo<dim>::rotational_states_of_faces/2].vertices[v];
-       std::swap (new_faces[face][rot].vertices[1],
-                  new_faces[face][rot].vertices[3]);
-      };
-  
-       
-                                  // now insert all the faces, by now
-                                  // without specific order with
-                                  // respect to the orientational
-                                  // states of the cell. note that we
-                                  // get the indices correct
-                                  // here. also remark that the face
-                                  // might already have been in the
-                                  // map, depending on whether a
-                                  // newighbor has already inserted
-                                  // it or not. we don't care about
-                                  // that here, though
-  std::map<Face,FaceData>::iterator
-    new_faces_ptr[internal::GridReorderingInfo<dim>::rotational_states_of_faces]
-    [GeometryInfo<dim>::faces_per_cell];
-  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-    for (unsigned int rot=0;
-        rot<internal::GridReorderingInfo<dim>::rotational_states_of_faces;
-        ++rot)
-      new_faces_ptr[rot][face]
-       = global_faces.insert (std::make_pair(new_faces[face][rot], FaceData())).first;
-  
-                                  // and crosslink them to each other
-  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-    for (unsigned int rot=0;
-        rot<internal::GridReorderingInfo<dim>::rotational_states_of_faces;
-        ++rot)
-      for (unsigned int other_rot=0;
-          other_rot<internal::GridReorderingInfo<dim>::rotational_states_of_faces;
-          ++other_rot)
-       {
-         if (other_rot < rot)
-           new_faces_ptr[rot][face]->second.reverse_faces[other_rot]
-             = new_faces_ptr[other_rot][face];
-         else
-           if (other_rot > rot)
-             new_faces_ptr[rot][face]->second.reverse_faces[other_rot-1]
-               = new_faces_ptr[other_rot][face];
-                                          // if rot==other_rot, then
-                                          // we need not link this
-                                          // cell to itself
-       };
-  
-
-                                  // for each of the faces (whether
-                                  // already inserted or not) note
-                                  // that the present cell is one of
-                                  // the neighbors
-  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-    {
-      if (new_faces_ptr[0][face]->second.adjacent_cells[0] ==
-         FaceData::invalid_adjacent_cell)
-       {
-                                          // no, faces had not been
-                                          // used before, so we are the
-                                          // first adjacent cell
-         for (unsigned int rot=0;
-              rot<internal::GridReorderingInfo<dim>::rotational_states_of_faces;
-              ++rot)
-           {
-             Assert (new_faces_ptr[rot][face]->second.adjacent_cells[0]
-                     == FaceData::invalid_adjacent_cell,
-                     ExcInternalError());
-             new_faces_ptr[rot][face]->second.adjacent_cells[0] = cell_no;
-           };
-       }
-      else
-       {       
-                                          // otherwise: cell had been
-                                          // entered before, so we are
-                                          // the second neighbor
-         const unsigned int
-           previous_neighbor = new_faces_ptr[0][face]->second.adjacent_cells[0];
-         for (unsigned int rot=0;
-              rot<internal::GridReorderingInfo<dim>::rotational_states_of_faces;
-              ++rot)
-           {
-             Assert (new_faces_ptr[rot][face]->second.adjacent_cells[0] ==
-                     previous_neighbor,
-                     ExcInternalError());
-             Assert (new_faces_ptr[rot][face]->second.adjacent_cells[1] ==
-                     FaceData::invalid_adjacent_cell,
-                     ExcInternalError());
-             new_faces_ptr[rot][face]->second.adjacent_cells[1] = cell_no;
-           };
-       };
-    };
-  
-  
-
-                                  // we still have to link cell in
-                                  // its 24 different orientations to
-                                  // the 6 faces in their
-                                  // orientations. now, there we
-                                  // could hardcode which faces in
-                                  // which rotation belong to the
-                                  // cell in each direction, but
-                                  // there a good reasons not to do
-                                  // so:
-                                  //
-                                  // first, this depends on that we
-                                  // know which orientation of the
-                                  // cell has which number, but this
-                                  // knowledge is hardcoded in the
-                                  // function CellData::rotate, so
-                                  // hardcoding it here again would
-                                  // mean redundancy, and would above
-                                  // that mean that we have to update
-                                  // two very different place if we
-                                  // chose to change one.
-                                  //
-                                  // second, finding out which face
-                                  // belongs to which cell is error
-                                  // prone, and one might get it
-                                  // wrong.
-                                  //
-                                  // the solution is: compute it once
-                                  // this function is first called
-                                  // based on the information from
-                                  // CellData::rotate and use that
-                                  // data in following calls to this
-                                  // function. the computed data has,
-                                  // of course, to be a static member
-                                  // function, and we store whether
-                                  // the data has been initialized
-                                  // already by checking the value of
-                                  // a special flag. furthermore, we
-                                  // guard the initialization by a
-                                  // thread mutex to make it
-                                  // thread-safe (in case someone
-                                  // wanted to read in two grids at
-                                  // the same time, for whatever
-                                  // reason).
-  static Threads::ThreadMutex initialization_lock;
-  initialization_lock.acquire ();
-
-  static bool already_initialized = false;
-  
-                                  // for each orientation of the
-                                  // cell, store in which orientation
-                                  // each of the six faces build the
-                                  // cell (store which face and which
-                                  // orientation):
-  static std::pair<unsigned int, unsigned int>
-    cell_orientation_faces[internal::GridReorderingInfo<dim>::rotational_states_of_cells][GeometryInfo<dim>::faces_per_cell];
-
-  if (already_initialized == false)
-    {
-      for (unsigned int rot=0; rot<internal::GridReorderingInfo<dim>::rotational_states_of_cells; ++rot)
-       {
-                                          // initialize a standard
-                                          // cell with the vertex
-                                          // numbers of the present
-                                          // cell we are working on
-         CellData<dim> standard_cell;
-         for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-           standard_cell.vertices[v] = this->vertices[v];
-
-                                          // then rotate it the given
-                                          // number of times
-         standard_cell.rotate (rot);
-
-                                          // then create the six
-                                          // faces of the thus
-                                          // rotated cell
-         const Face standard_faces[GeometryInfo<dim>::faces_per_cell]
-           = { { { standard_cell.vertices[0], standard_cell.vertices[1],
-                   standard_cell.vertices[2], standard_cell.vertices[3] } },
-               { { standard_cell.vertices[4], standard_cell.vertices[5],
-                   standard_cell.vertices[6], standard_cell.vertices[7] } },
-               { { standard_cell.vertices[0], standard_cell.vertices[1],
-                   standard_cell.vertices[5], standard_cell.vertices[4] } },
-               { { standard_cell.vertices[1], standard_cell.vertices[5],
-                   standard_cell.vertices[6], standard_cell.vertices[2] } },
-               { { standard_cell.vertices[3], standard_cell.vertices[2],
-                   standard_cell.vertices[6], standard_cell.vertices[7] } },
-               { { standard_cell.vertices[0], standard_cell.vertices[4],
-                   standard_cell.vertices[7], standard_cell.vertices[3] } } };
-
-                                          // then try to identify
-                                          // these faces in the ones
-                                          // we have already created
-         for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-           {
-             bool face_found = false;
-             for (unsigned int f=0;
-                  (!face_found) && (f<GeometryInfo<dim>::faces_per_cell); ++f)
-               for (unsigned int r=0;
-                    r<internal::GridReorderingInfo<dim>::rotational_states_of_faces;
-                    ++r)
-                 if (standard_faces[face] == new_faces[f][r])
-                   {
-                     cell_orientation_faces[rot][face] = std::make_pair(f,r);
-                     face_found = true;
-                     break;
-                   };
-
-                                              // make sure that we
-                                              // have found something
-                                              // indeed
-             Assert (face_found == true, ExcInternalError());
-           };
-
-                                          // more checks: make sure
-                                          // that each of the
-                                          // original faces appears
-                                          // in one rotation or other
-                                          // as face of the present
-                                          // cell in its orientation
-                                          // we currently check. as
-                                          // we don't call this part
-                                          // of the program too
-                                          // often, don't make
-                                          // differences between
-                                          // debug and optimized mode
-         std::vector<bool> face_used(GeometryInfo<dim>::faces_per_cell, false);
-         for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-           {
-                                              // ups, face already
-                                              // used? can't be!
-             Assert (face_used[face] == false, ExcInternalError());
-             face_used[face] = true;
-           };
-                                          // since we have checked
-                                          // that each face has not
-                                          // been used previously, we
-                                          // also know that all faces
-                                          // have been used exactly
-                                          // once, so no more checks
-                                          // necessary
-       };
-
-                                      // that's it: we now know which
-                                      // faces build up this cell in
-                                      // each of its possible
-                                      // orientations
-      already_initialized = true;
-    };
-                                  // initialization is done, so
-                                  // release the lock and let other
-                                  // threads run
-  initialization_lock.release ();
-
-                                  // now we can use the information:
-                                  // link the faces in their
-                                  // directions to the cell in each
-                                  // of its orientations
-  for (unsigned int rot=0;
-       rot<internal::GridReorderingInfo<dim>::rotational_states_of_cells;
-       ++rot)
-    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-      faces[rot][face] = new_faces_ptr
-                        [cell_orientation_faces[rot][face].second]
-                        [cell_orientation_faces[rot][face].first];
-}
-
-
-
-template <int dim>
-void GridReordering<dim>::Cell::fix_cell_neighbors ()
-{
-  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-    {
-                                      // 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
-  track_back_to_cell = cell_no-1;
-  if (cell_no == 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
-{
-                                  // 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<internal::GridReorderingInfo<dim>::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)
-{
-  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 <int dim>
-inline
-void GridReordering<dim>::Cell::mark_faces_unused (const unsigned int rot)
-{
-  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>
-bool GridReordering<dim>::Face::operator == (const Face &face) const
-{
-  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
-    if (vertices[v] != face.vertices[v])
-      return false;
-  return true;
-}
-
-
-
-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 (std::vector<Cell> &cells,
-                                     RotationStack     &rotation_states,
-                                     const unsigned     tbtc)
-{
-  unsigned int track_back_to_cell = tbtc;
-  
-  top_of_function:
-  
-  Assert (track_back_to_cell > 0, ExcInternalError());
-
-  unsigned int last_rotation_state = static_cast<unsigned int>(-1);
-  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.back();
-      
-                                      // 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_back();
-    };
-  Assert (last_rotation_state < internal::GridReorderingInfo<dim>::rotational_states_of_cells, ExcInternalError());
-  
-                                  // 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 < internal::GridReorderingInfo<dim>::rotational_states_of_cells-1)
-    {
-                                      // possible. push that state to
-                                      // the stack and leave
-      rotation_states.push_back (last_rotation_state+1);
-      return;
-    }
-  else
-    {
-                                      // last cell can't be rotated
-                                      // further. go on with
-                                      // backtracking
-      const typename std::vector<Cell>::iterator
-       try_cell = cells.begin() + rotation_states.size();
-      
-      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>
-bool GridReordering<dim>::try_rotate_single_neighbors (std::vector<Cell> &cells,
-                                                      RotationStack     &rotation_states)
-{
-                                  // the rotation state of the cell
-                                  // which we try to add by rotating
-                                  // neighbors has already been
-                                  // popped from the stack, so we get
-                                  // its number like this:
-  const unsigned int cell_no = rotation_states.size();
-
-                                  // now try each of the neighbors
-                                  // that have already been added to
-                                  // the grid. don't try the cell
-                                  // that we will track back to
-                                  // anyway if this operation should
-                                  // fail
-  for (unsigned int neighbor=0; neighbor<GeometryInfo<dim>::faces_per_cell; ++neighbor)
-    if (cells[cell_no].neighbors[neighbor] < cell_no)
-      if (cells[cell_no].neighbors[neighbor] != cells[cell_no].track_back_to_cell)
-       {
-         const unsigned int neighbor_no = cells[cell_no].neighbors[neighbor];
-         const unsigned int old_rotation_state = rotation_states[neighbor_no];
-         
-                                          // unlink faces used by the
-                                          // present rotation state
-         cells[neighbor_no].mark_faces_unused (old_rotation_state);
-
-                                          // then try all rotation
-                                          // states besides the ones
-                                          // that have already been
-                                          // tried:
-         for (unsigned int neighbor_rot=old_rotation_state+1;
-              neighbor_rot<internal::GridReorderingInfo<dim>::rotational_states_of_cells; ++neighbor_rot)
-           {
-                                              // first, if the
-                                              // neighbor itself does
-                                              // not fit in the grid,
-                                              // then there is
-                                              // nothing to do
-             if (! cells[neighbor_no].check_consistency (neighbor_rot))
-               continue;
-
-                                              // however, if the
-                                              // neighbor worked,
-                                              // then mark its faces
-                                              // as used
-                                              // preliminarily and
-                                              // try to fit in the
-                                              // present cell in some
-                                              // orientation
-             cells[neighbor_no].mark_faces_used (neighbor_rot);
-
-             for (unsigned int cell_rot=0; cell_rot<internal::GridReorderingInfo<dim>::rotational_states_of_cells; ++cell_rot)
-               if (cells[cell_no].check_consistency (cell_rot) == true)
-                 {
-                                                    // ah, see,
-                                                    // this
-                                                    // combination
-                                                    // of neighbor
-                                                    // rotation and
-                                                    // this cell
-                                                    // works. enter
-                                                    // the
-                                                    // respective
-                                                    // states into
-                                                    // the arrays
-                                                    // and leave
-                                                    // with success
-                   rotation_states[neighbor_no] = neighbor_rot;
-                   
-                   rotation_states.push_back (cell_rot);
-                   cells[cell_no].mark_faces_used (cell_rot);
-
-                   return true;
-                 };
-             
-                                              // no, there was no
-                                              // way to fit the
-                                              // present cell into
-                                              // the grid given
-                                              // this orientation
-                                              // of the
-                                              // neighbor. discard
-                                              // this attempt and
-                                              // try that neighbors
-                                              // next rotation
-             cells[neighbor_no].mark_faces_unused (neighbor_rot);
-           };
-         
-                                          // there was no way to
-                                          // rotate this neighbor so
-                                          // that the present cell
-                                          // fit into the
-                                          // grid. reinstantiate the
-                                          // old state and go on to
-                                          // the next neighbor
-         cells[neighbor_no].mark_faces_used (old_rotation_state);
-       };
-
-                                  // rotation of neighbors did not
-                                  // help this cell, there is no
-                                  // other way than to do a full
-                                  // backtracking
-  return false;
-}
-
-
-
-template <int dim>
-void GridReordering<dim>::find_reordering (std::vector<Cell>           &cells,
-                                          std::vector<CellData<dim> > &original_cells,
-                                          const std::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
-  RotationStack 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_back (0);
-  cells[0].mark_faces_used (rotation_states.back());
-  
-  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_back (0);
-
-                                      // check whether the present
-                                      // cell in the present
-                                      // orientation is valid
-      check_topmost_cell:
-      
-      const typename std::vector<Cell>::iterator
-       try_cell = cells.begin() + rotation_states.size()-1;
-      if (try_cell->check_consistency (rotation_states.back()))
-       {
-                                          // 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.back());
-         
-         continue;
-       }
-      else
-       {
-                                          // no, doesn't work. see if
-                                          // we can rotate the top
-                                          // cell so that it works
-         if (rotation_states.back()+1
-             < internal::GridReorderingInfo<dim>::rotational_states_of_cells)
-           {
-                                              // yes, can be
-                                              // done. then do so and
-                                              // check again
-             ++rotation_states.back();
-             goto check_topmost_cell;
-           }
-         else
-           {
-                                              // no, no more
-                                              // orientation of the
-                                              // top cell possible,
-                                              // we have to backtrack
-                                              // some way
-
-                                              // first pop rotational
-                                              // state of top cell,
-                                              // since for that no
-                                              // faces have been
-                                              // marked as used yet
-             rotation_states.pop_back();
-
-                                              // in general, if we
-                                              // fail to insert the
-                                              // present cell somehow
-                                              // into the existing
-                                              // part of the grid,
-                                              // then we track back
-                                              // to the neighbor of
-                                              // the failed cell with
-                                              // the highest cell
-                                              // index below the
-                                              // index of the present
-                                              // cell. however,
-                                              // before we do so, we
-                                              // try a simple
-                                              // heuristic: if
-                                              // rotating single
-                                              // neighbors a little
-                                              // helps the process
-                                              // somewhat:
-             const bool rotation_helps
-               = try_rotate_single_neighbors (cells, rotation_states);
-
-                                              // if rotation helped,
-                                              // then go on to the
-                                              // next cell. the
-                                              // called function has
-                                              // already marked the
-                                              // respective faces as
-                                              // used and has pushed
-                                              // the rotation state
-                                              // of the present cell
-                                              // to the stack
-             if (rotation_helps == true)
-               continue;
-
-                                              // if that failed to
-                                              // help, 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 = std::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.back());
-
-                                      // to check the correctness of
-                                      // the program up to here:
-                                      // unmark the cells' faces to
-                                      // check whether they have all
-                                      // correctly declared they
-                                      // use. checking this is done
-                                      // in the calling function, as
-                                      // only that has direct access
-                                      // to the map of faces (this
-                                      // function only accesses it
-                                      // through pointers stored in
-                                      // the cells)
-      cells[new_cell_number].mark_faces_unused (rotation_states.back());
-
-                                      // then delete this rotational
-                                      // state as we don't need it
-                                      // any more
-      rotation_states.pop_back ();
-    };
-}
-
-
-
-template <int dim>
-std::vector<unsigned int>
-GridReordering<dim>::presort_cells (std::vector<Cell>       &cells,
-                                   std::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;
-
-
-                                  // 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);
-  std::vector<unsigned int> new_cell_numbers (cells.size(), invalid_cell_number);
-
-  unsigned int next_free_new_number = 0;
-  
-                                   // loop over each connected part of
-                                   // the domain. since the domain may
-                                   // consist of different unconnected
-                                   // parts, we have to loop until
-                                   // there are no more unnumbered
-                                   // cells
-  while (next_free_new_number < cells.size())
-   {
-                                     // for initialization of
-                                     // min_neighbors, go to the
-                                     // first cell of this part of
-                                     // the domain that has an
-                                     // invalid cell number, i.e. has
-                                     // not yet been renumbered
-     for (unsigned int i=0; i<cells.size(); ++i)
-       if (new_cell_numbers[i]==invalid_cell_number)
-        {
-          min_neighbors           = cells[i].count_neighbors();
-          cell_with_min_neighbors = i;
-          break;
-        }
-     
-                                     // check, if we have an as yet
-                                     // unnumbered cell with less
-                                     // neighbors than the first
-                                     // found cell
-     for (unsigned int i=1; i<cells.size(); ++i)
-       if ((min_neighbors > cells[i].count_neighbors()) &&
-          (new_cell_numbers[i]==invalid_cell_number))
-        {  
-          min_neighbors = cells[i].count_neighbors();
-          cell_with_min_neighbors = i;
-          if (min_neighbors == 1)
-                                             // better is not possible
-            break;
-        };
-
-                                     // have an array of the next
-                                     // cells to be numbered (old numbers)
-     std::vector<unsigned int> next_round_cells (1, cell_with_min_neighbors);
-  
-                                     // 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
-        std::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::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]);
-        
-        
-                                         // 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
-        std::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));
-       };
-   };   // end of loop over subdomains
-  
-    
-  Assert (std::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:
-  std::vector<Cell> 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
-  std::swap (cells, new_cells);
-  
-                                  // now we still have to convert all
-                                  // old cell numbers to new cells
-                                  // numbers. non-existent neighbors
-                                  // (with index -1) are mapped to
-                                  // non-existent neighbors, so we
-                                  // need not touch these indices
-  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)
-       if (cells[c].neighbors[n] != Cell::invalid_neighbor)
-         {
-           Assert (cells[c].neighbors[n] < new_cell_numbers.size(),
-                   ExcIndexRange(cells[c].neighbors[n], 0, 
-                                 new_cell_numbers.size()));
-           cells[c].neighbors[n] = new_cell_numbers[cells[c].neighbors[n]];
-         };
-    };
-
-  for (typename std::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 (std::vector<CellData<dim> > &original_cells)
-{
-                                  // 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
-  std::vector<Cell> cells;
-  cells.reserve (original_cells.size());
-  for (unsigned int i=0; i<original_cells.size(); ++i)
-    cells.push_back (Cell(original_cells[i], i));
-  
-                                  // first generate all the faces
-                                  // possible, i.e. in each possible
-                                  // direction and rotational state
-  std::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)
-    {
-      Cell &cell = cells[cell_no];
-      cell.fix_cell_neighbors ();
-    };
-
-
-                                  // do a preordering step in order
-                                  // to make further backtracking
-                                  // more local
-  const std::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);
-
-
-  
-                                  // finally check the consistency of
-                                  // the program by ensuring that all
-                                  // faces have no use-marks any
-                                  // more. to this end, the
-                                  // find_reordering function has
-                                  // cleared all used marks it knows
-                                  // of
-  for (typename std::map<Face,FaceData>::iterator i=faces.begin(); i!=faces.end(); ++i)
-    Assert (i->second.use_count == 0, ExcInternalError());
-}
-
-#endif // deal_II_dimension == 3
-
 #if deal_II_dimension == 2
 
 namespace internal
@@ -1722,13 +571,821 @@ void GridReordering<2>::reorder_cells (std::vector<CellData<2> > &original_cells
 #endif
 
 
-// explicit instantiations. only require the main function, it should
-// then claim whatever templates it needs. note that in 1d, the
-// respective function is already specialized, and in 2d we have an
-// explicit specialization of the whole class
+
 #if deal_II_dimension == 3
-template
-void
-GridReordering<deal_II_dimension>::
-reorder_cells (std::vector<CellData<deal_II_dimension> > &);
-#endif
+
+namespace internal
+{
+  namespace GridReordering3d
+  {
+    DeclException1 (GridOrientError,
+                   std::string,
+                   <<  "Grid Orientation Error"<< arg1);
+
+
+    //Switched two ints so that v1<v2;
+    inline void sort2(int &v1, int &v2)
+    {
+      if(v1>v2)
+       {
+         int t=v1;
+         v1=v2;
+         v2=t;
+       }
+    }
+
+    CheapEdge::CheapEdge(int n0, int n1) : node0(n0), node1(n1)
+    {
+      //sort the entries so that node0<node1;
+      sort2(node0,node1);
+    }
+
+    bool CheapEdge::operator<(const CheapEdge & e2) const
+    {
+      if ((node0)<(e2.node0)) return true;
+      if ((node0)>(e2.node0)) return false;
+      if ((node1)<(e2.node1)) return true;
+      return false;
+    };
+
+  
+    // This is the guts of the matter...
+    void build_mesh(Mesh &m)
+    {
+      const ElementInfo &info = m.info; 
+      std::vector<Cell> & cell_list = m.cell_list;
+      std::vector<Edge> & edge_list = m.edge_list;
+
+      const unsigned int cell_list_length=cell_list.size();
+
+
+      unsigned int NumEdges=0;
+      // Corectly build the edge list
+      {
+       // edge_map stores the edge_number associated with a given CheapEdge
+       std::map<CheapEdge,int> edge_map;
+       unsigned int ctr=0;
+       for(unsigned int cur_cell_id=0; 
+           cur_cell_id<cell_list_length; 
+           ++cur_cell_id)
+         {
+           //Get the local node numbers on edge edge_num
+           Cell & cur_cell = cell_list[cur_cell_id];
+           //m.DumpCell(cur_cell);
+           for(unsigned short int edge_num=0; 
+               edge_num<12; 
+               ++edge_num)
+             {
+               unsigned int gl_edge_num=0;
+               int l_edge_orient=1;
+               //Construct the CheapEdge
+               int node0=cur_cell.nodes[info.nodes_on_edge[edge_num][0]];
+               int node1=cur_cell.nodes[info.nodes_on_edge[edge_num][1]];
+               CheapEdge cur_edge(node0,node1);
+               if(edge_map.count(cur_edge)==0) // Edge not in map
+                 {
+                   // put edge in hash map with ctr value;
+                   edge_map[cur_edge]=ctr;
+                   gl_edge_num=ctr;
+                   edge_list.push_back(Edge(node0,node1)); //put the edge into the global edge list.
+                   ctr++;
+                 }
+               else
+                 {
+                   gl_edge_num=edge_map[cur_edge]; // get edge_num from hash_map;
+                   if (edge_list[gl_edge_num].nodes[0]!=node0)
+                     { 
+                       l_edge_orient=-1;
+                     }
+                 }
+               cell_list[cur_cell_id].edges[edge_num]=gl_edge_num; // set edge number to edgenum;
+               cell_list[cur_cell_id].local_orientation_flags[edge_num]=l_edge_orient;
+             }
+         }
+       NumEdges=ctr;
+      }
+
+      //Count each of the edges.
+      {
+       std::vector<int> edge_count(NumEdges,0);
+      
+
+       //Count every time an edge occurs in a cube.
+       for(unsigned int cur_cell_id=0; 
+           cur_cell_id<cell_list_length; 
+           ++cur_cell_id)
+         {
+           Cell & cur_cell = cell_list[cur_cell_id];
+           for(unsigned short int edge_num=0; 
+               edge_num<12; 
+               ++edge_num)
+             {
+               edge_count[cur_cell.edges[edge_num]]++;
+             } 
+         }
+
+       // So we now know howmany cubes contain a given edge.
+       // Just need to store the list of cubes in the edge
+
+       //Alocate the space for the neighbour list
+       for(unsigned int cur_edge_id = 0;
+           cur_edge_id<NumEdges;
+           ++cur_edge_id)
+         {
+           Edge & cur_edge = edge_list[cur_edge_id];
+           unsigned int NN = edge_count[cur_edge_id];
+           cur_edge.num_neighbouring_cubes=NN;
+           cur_edge.neighbouring_cubes = new unsigned int [NN];
+         }
+      
+       //Stores the position of the current neighbour in the edge's neighbour list
+       std::vector<int> cur_cell_edge_list_posn(NumEdges,0);
+       for(unsigned int cur_cell_id =0;
+           cur_cell_id<cell_list_length;
+           ++cur_cell_id)
+         {
+           Cell & cur_cell = cell_list[cur_cell_id];
+           for(unsigned short int edge_num=0; 
+               edge_num<12; 
+               ++edge_num)
+             {
+               unsigned int gl_edge_id=cur_cell.edges[edge_num];
+               Edge & cur_edge = edge_list[gl_edge_id];
+               cur_edge.neighbouring_cubes[cur_cell_edge_list_posn[gl_edge_id]]=cur_cell_id;
+               cur_cell_edge_list_posn[gl_edge_id]++;
+             }
+         }
+      
+      }
+    }
+
+    // Deal Element Information
+    const int DealEdgeToNodeArray[8][3] = 
+      { 
+       {0,4,8},
+       {0,5,9},
+       {3,5,10},
+       {3,4,11},
+       {1,7,8},
+       {1,6,9},
+       {2,6,10},
+       {2,7,11}
+      };
+
+    const int DealEdgeOrientArray[8][3] =
+      {
+       { 1, 1, 1},
+       {-1, 1, 1},
+       {-1,-1, 1},
+       { 1,-1, 1},
+       { 1, 1,-1},
+       {-1, 1,-1},
+       {-1,-1,-1}, 
+       { 1,-1,-1}
+      };
+
+    const int DealNodesOnEdgeArray[12][2] =
+      {
+       {0,1},
+       {4,5},
+       {7,6},
+       {3,2},
+       {0,3},
+       {1,2},
+       {5,6},
+       {4,7},
+       {0,4},
+       {1,5},
+       {2,6},
+       {3,7}
+      };
+
+    //Starting at emination node (for edges) and chosing clockwise order
+    const int DealFaceNodeArray[6][4] = //TODO: HERE
+      {
+       {0,1,2,3},
+       {0,4,5,1},
+       {1,5,6,2},
+       {3,2,6,7},
+       {0,3,7,4},
+       {4,7,6,5}
+      };
+
+    DealElemInfo::DealElemInfo() : ElementInfo()
+    {
+      for(int node=0; node<8; ++node)
+       {
+         for(int i=0;i<3;++i)
+           {
+             edge_to_node[node][i]=DealEdgeToNodeArray[node][i];
+             edge_to_node_orient[node][i]=DealEdgeOrientArray[node][i];
+           }
+       }
+
+      for (int edge=0;edge<12; ++edge)
+       {
+         nodes_on_edge[edge][0]=DealNodesOnEdgeArray[edge][0];
+         nodes_on_edge[edge][1]=DealNodesOnEdgeArray[edge][1];
+       }
+
+      for(int facenum=0;facenum<6;++facenum)
+       for(int nodenum=0;nodenum<4;++nodenum)
+         nodes_on_face[facenum][nodenum]=DealFaceNodeArray[facenum][nodenum];
+    }
+
+
+
+
+
+
+    Edge::~Edge()
+    {
+      delete [] neighbouring_cubes;
+    }
+
+    bool Mesh::sanity_check() const
+    {
+      bool retval=true;
+      for(unsigned int i=0; i<cell_list.size(); ++i)
+       retval&=sanity_check(i);
+      return retval;
+    }
+
+    bool Mesh::sanity_check(int cellnum) const
+    {
+      //Should check that every edge coming into a node has the 
+      //same node value
+      bool retval=true;
+      for(int i=0;i<8;++i)
+       {
+         retval&=sanity_check_node(cellnum,i);
+       }
+      return retval;
+    }
+
+    bool Mesh::sanity_check_node(int cell_num, int local_node_num) const
+    {
+      // Get the Local Node Numbers of the incoming edges
+      int e0 = info.edge_to_node[local_node_num][0];
+      int e1 = info.edge_to_node[local_node_num][1]; 
+      int e2 = info.edge_to_node[local_node_num][2];
+
+      // Global Edge Numbers
+      const Cell& c = cell_list[cell_num];
+  
+      int ge0 = c.edges[e0];
+      int ge1 = c.edges[e1];
+      int ge2 = c.edges[e2];
+
+
+      //  std::cout<<"Coming Into Node "<<node_num<< " are the local edges : "
+      //    <<e0<<" "<<e1<<" "<<e2<<std::endl;
+
+      int or0 = info.edge_to_node_orient[local_node_num][0]*c.local_orientation_flags[e0];
+      int or1 = info.edge_to_node_orient[local_node_num][1]*c.local_orientation_flags[e1];
+      int or2 = info.edge_to_node_orient[local_node_num][2]*c.local_orientation_flags[e2];
+
+      //  std::cout<<"They have orientations : "<<or0<<" "<<or1<<" "<<or2<<std::endl;
+
+      //What each edge thinks the current node should be.
+
+      int curglobalnodenum0 = edge_list[ge0].nodes[or0==1 ? 0 : 1];
+      int curglobalnodenum1 = edge_list[ge1].nodes[or1==1 ? 0 : 1];
+      int curglobalnodenum2 = edge_list[ge2].nodes[or2==1 ? 0 : 1];
+
+      //  std::cout<<"This means the current node is "
+      //    <<curglobalnodenum0 <<" "
+      //    <<curglobalnodenum1 <<" "
+      //    <<curglobalnodenum2 <<std::endl; 
+
+      bool retval = ((curglobalnodenum0 == curglobalnodenum1)&&
+                    (curglobalnodenum1 == curglobalnodenum2) );
+
+
+      if (!retval)
+       {
+         std::cout<<"FAILED SANITY TEST";
+         dump_cell(c);
+       }
+      Assert (retval == true, ExcInternalError());
+  
+      return retval;
+    }
+
+    void Mesh::dump_cell(const Cell &c) const
+    {
+      std::cout<<std::endl
+              <<"===CELL NODES==="<<std::endl;
+      for(int i=0;i<8;++i)
+       std::cout<<"\t"<<c.nodes[i];
+      std::cout<<std::endl;
+      std::cout<<"===CELL EDGES==="<<std::endl;
+      for(int i=0;i<12;++i)
+       {
+         std::cout<<"\t"<<c.edges[i]<<" "<<c.local_orientation_flags[i];
+         if(c.edges[i]>=0)
+           {
+             std::cout<<":"<<edge_list[c.edges[i]].nodes[0]<<" ";
+             std::cout<<":"<<edge_list[c.edges[i]].nodes[1]<<" ";
+             std::cout<<":"<<edge_list[c.edges[i]].orientation_flag;
+           }
+         std::cout<<std::endl;
+       }
+    }
+
+
+    void Mesh::dump() const
+    {
+      std::cout<<std::endl
+              <<"===NODES==="<<std::endl;
+      const int nnodes=node_list.size();
+      for(int i=0;i<nnodes;++i)
+       {
+         std::cout<<i<<"\t"<<node_list[i](0)<<"\t"<<node_list[i](1)<<"\t"<<node_list[i](2)<<std::endl;
+       }
+
+      std::cout<<"===EDGES==="<<std::endl;
+      const unsigned int nedges=edge_list.size();
+      for(unsigned int i=0;i<nedges;++i)
+       {
+         const Edge &e = edge_list[i];
+         std::cout<<i<<"\t"<<e.orientation_flag<<"\t"<<e.nodes[0]<<"\t"<<e.nodes[1]<<std::endl;
+       }
+  
+      std::cout<<"===CELLS==="<<std::endl;
+      const int ncells=cell_list.size();
+      for(int i=0;i<ncells;++i)
+       {
+         const Cell & c = cell_list[i];
+         std::cout<<"cell "<<i<<std::endl<<"  nodes:\t";
+         std::cout<<c.nodes[0]<<"\t"<<c.nodes[1]<<"\t"<<c.nodes[2]<<"\t"<<c.nodes[3]<<std::endl;
+         std::cout<<"\t\t"<<c.nodes[4]<<"\t"<<c.nodes[5]<<"\t"<<c.nodes[6]<<"\t"<<c.nodes[7]<<std::endl;
+         std::cout<<"  edges:"<<std::endl;
+         for(int j=0;j<12;++j)
+           std::cout<<"\t\t"<<c.edges[j]<<" "<<c.local_orientation_flags[j]<<std::endl;
+       }
+    }
+       
+    void Mesh::dump_edges(char const * const fname) const
+    {
+      const int nedges = edge_list.size();
+      const int npoints = node_list.size();
+      //Only do this if we've have the extra information
+      if (npoints==0) 
+       return;
+  
+      std::ofstream outfile(fname);
+
+      outfile<<npoints<<" "<<nedges<<std::endl;
+      for(int i=0;i<npoints;++i)
+       {
+         const Point<3> & n = node_list[i];
+         outfile<<n(0)<<" "<<n(1)<<" "<<n(2)<<std::endl;
+       } 
+      for(int i=0;i<nedges;++i)
+       {
+         const Edge & e = edge_list[i];
+         outfile<<e.nodes[0]<<" "
+                <<e.nodes[1]<<" "
+                <<e.orientation_flag<<" "
+                <<e.group<<std::endl;
+       } 
+    }
+
+
+    /**
+     * This assignes an orientation to each edge so that 
+     * every cube is a rotated Deal.II cube.
+     */
+    bool Orienter::orient_edges(Mesh &m)
+    {
+  
+      // First check that the mesh is sensible
+      AssertThrow(m.sanity_check(),GridOrientError("Invalid Mesh Detected"));
+
+      // We start by looking only at the first cube.  
+      cur_posn=0; 
+      marker_cube=0;
+
+      // We mark each edge with a group number 
+      // (mostly for mesh debugging purposes)
+      cur_edge_group=1;
+      // While there are still cubes to orient
+      while(GetNextUnorientedCube(m))
+       {
+         // And there are edges in the cube to orient
+         while(OrientNextUnorientedEdge(m))
+           {
+             // Make all the sides in the current set match
+             OrientEdgesInCurrentCube(m);
+             // Add the adjacent cubes to the list for processing
+             GetAdjacentCubes(m);
+             // Start working on this list of cubes
+             while(GetNextActiveCube(m))
+               {
+                 // Make sure the Cube doesn't have a contradiction 
+                 if(!Consistant(m,cur_posn))
+                   {
+                     m.dump_edges("edgelist.dat");
+                   }
+                 AssertThrow(Consistant(m,cur_posn),GridOrientError("Mesh is Unorientable"));
+                 // If we needed to orient any edges in the current cube
+                 // then we may have to process the neighbour.
+                 if(OrientEdgesInCurrentCube(m)) GetAdjacentCubes(m);
+               }
+             cur_edge_group++;
+           }
+       }
+      return true;
+    }
+
+    bool Orienter::GetNextUnorientedCube(Mesh &m)
+    {
+      // The last cube in the list
+      unsigned int end_cube_num = m.cell_list.size();
+      // Keep shifting along the list until we find a cube which is not
+      // fully oriented or the end.
+      while( (marker_cube<end_cube_num)&&(is_oriented(m,marker_cube)) )
+       marker_cube++;
+      cur_posn=marker_cube;
+      //Return true if we now point at a valid cube.
+      return cur_posn<end_cube_num;
+    }
+
+    bool Orienter::is_oriented(const Mesh &m, int cell_num)
+    {
+      const Cell& c = m.cell_list[cell_num];
+      for(int i = 0; i<12; ++i)
+       {
+         int edgenum = c.edges[i];
+         if (m.edge_list[edgenum].orientation_flag==0) return false;
+       }
+      return true;
+    }
+
+    bool Orienter::Consistant(Mesh &m, int cell_num)
+    {
+
+      const Cell& c = m.cell_list[cell_num];
+  
+      //Checks that all oriented edges in the group are oriented consistantly.  
+      for(int group=0; group<3; ++group)
+       {
+         // When a nonzero orientation is first encoutered in the group it is 
+         // stored in this
+         int value=0;
+         //Loop over all parallel edges
+         for(int i=4*group;i<4*(group+1);++i)
+           {
+             //The local edge orientation within the cell
+             int LOR = c.local_orientation_flags[i] *
+               m.edge_list[c.edges[i]].orientation_flag; 
+             //If the edge has orientation
+             if (LOR!=0)
+               {
+                 //And we haven't seen an oriented edge before
+                 if (value==0) 
+                   {
+                     //Store it's value
+                     value=LOR;
+                   }
+                 else
+                   {
+                     //If we have seen a oriented edge in this group
+                     //we'd better have the same orientation.
+                     if (value!=LOR) return false;
+                   }
+               }
+           }
+       }
+      return true;
+    }
+
+    bool Orienter::OrientNextUnorientedEdge(Mesh &m)
+    {
+      cur_posn=marker_cube;
+      const Cell& c = m.cell_list[cur_posn];
+      int i=0;
+      // search for the unoriented side
+      while((i<12)&&(m.edge_list[c.edges[i]].orientation_flag!=0))
+       ++i;
+  
+      //if we found none then return false
+      if (i==12)
+       return false;
+  
+      //Which edge group we're in.
+      int egrp=i/4;
+
+      //A sanity check that none of the other edges in the group have been oriented yet
+      //Each of the edges in the group should be un-oriented
+      for(int j=egrp*4;j<egrp*4+4;++j)
+       {
+         Assert(m.edge_list[c.edges[j]].orientation_flag==0,
+                GridOrientError("Tried to orient edge when other edges in group already oriented!"));
+       }
+
+      //Make the edge alignment match that of the local cube.
+      m.edge_list[c.edges[i]].orientation_flag
+       = c.local_orientation_flags[i];
+      m.edge_list[c.edges[i]].group=cur_edge_group;
+
+      edge_orient_array[i]=true;
+
+      return true;
+    }
+
+    bool Orienter::OrientEdgesInCurrentCube(Mesh &m)
+    {
+      bool retval = false;
+      for(int i=0;i<3;++i)
+       retval=retval||OrientEdgeSetInCurrentCube(m,i);
+      return retval;
+    }
+
+    bool Orienter::OrientEdgeSetInCurrentCube(Mesh &m,int n)
+    {
+      const Cell& c = m.cell_list[cur_posn];
+  
+      //Check if Any edge is oriented
+      int num_oriented =0 ;
+      int glorient = 0;
+      unsigned int edge_flags = 0;
+      unsigned int cur_flag = 1;
+      for(int i = 4*n; i<4*(n+1); ++i)
+       {
+         int orient = m.edge_list[c.edges[i]].orientation_flag *
+           c.local_orientation_flags[i];
+         if (orient!=0)
+           {
+             num_oriented++;
+             if (glorient!=0)
+               {
+                 AssertThrow(orient==glorient,GridOrientError("Attempted to Orient Misaligned cube"));
+               }
+             else
+               {
+                 glorient=orient;
+               }
+           }
+         else
+           {
+             edge_flags|=cur_flag;
+           }
+         cur_flag*=2;
+       }
+
+      //were any of the sides oriented?
+      //were they all already oriented?
+      if ((glorient==0)||(num_oriented==4))
+       return false;
+
+      //If so orient all edges consistantly.
+      cur_flag = 1;
+      for(int i=4*n; i<4*(n+1); ++i)
+       {
+         //std::cout<<i<<" ORIENTING\n";
+         if ((edge_flags&cur_flag)!=0)
+           {
+             m.edge_list[c.edges[i]].orientation_flag 
+               = c.local_orientation_flags[i]*glorient;
+             m.edge_list[c.edges[i]].group=cur_edge_group;
+             edge_orient_array[i]=true;
+           }
+         cur_flag*=2;
+       } 
+      return true;
+    }
+
+    void Orienter::GetAdjacentCubes(Mesh &m)
+    {
+      const Cell& c = m.cell_list[cur_posn];
+      for(unsigned int e=0;e<12;++e)
+       {
+         if(edge_orient_array[e])
+           {
+             edge_orient_array[e]=false;
+             const unsigned int cur_local_edge_num = e;
+
+             Edge & the_edge=m.edge_list[c.edges[cur_local_edge_num]];
+             for(unsigned int local_cube_num = 0; 
+                 local_cube_num < the_edge.num_neighbouring_cubes; 
+                 ++local_cube_num)
+               {
+                 unsigned int global_cell_num = the_edge.neighbouring_cubes[local_cube_num];
+                 Cell& ncell=m.cell_list[global_cell_num];
+                 if (!ncell.waiting_to_be_processed)
+                   {
+                     SheetToProcess.push_back(global_cell_num);
+                     ncell.waiting_to_be_processed=true;
+                   }
+
+               }
+           }
+
+       }
+    }  
+
+    bool Orienter::GetNextActiveCube(Mesh &m)
+    {
+      //Mark the curent Cube as finnished with.
+      Cell &c = m.cell_list[cur_posn];
+      c.waiting_to_be_processed=false;
+      if(SheetToProcess.size()!=0)
+       {
+         cur_posn=SheetToProcess.back();
+         SheetToProcess.pop_back();
+         return true;
+       }
+      return false;
+    }
+
+    bool Orienter::CheckCellEdgeGroupConsistancy(const Mesh &m, const Cell & c, int egrp) const
+    {
+      int grp=0;
+      for(int i=4*egrp;i<4*egrp+4;++i)
+       {
+         int cgrp = m.edge_list[c.edges[i]].group;
+         if (cgrp!=0)
+           {
+             if(grp==0)
+               {
+                 grp=cgrp;
+               }
+             else if (grp!=cgrp)
+               {
+                 return false;
+               }
+           }
+       }
+      return true;
+    }
+      
+    
+    bool Orienter::CheckCellEdgeGroupConsistancy(const Mesh &m, const Cell &c) const
+    {
+      return ( 
+             CheckCellEdgeGroupConsistancy(m,c,0) &&
+             CheckCellEdgeGroupConsistancy(m,c,1) &&
+             CheckCellEdgeGroupConsistancy(m,c,2)
+             );
+    }
+
+    void Orienter::orient_cubes(Mesh & the_mesh)
+    {
+      // We assume that the mesh has all edges oriented already.
+  
+      const unsigned int numelems=the_mesh.cell_list.size();
+  
+      // This is a list of permutations that take node 0 to node i
+      // but only rotate the cube. 
+      // (This set is far from unique (there are 3 for each node - for our
+      // algorithm it doesn't matter which of the three we use)
+      const unsigned int CubePermutations[8][8] = {
+       {0,1,2,3,4,5,6,7},
+       {1,2,3,0,5,6,7,4},
+       {2,3,0,1,6,7,4,5},
+       {3,0,1,2,7,4,5,6},
+       {4,7,6,5,0,3,2,1},
+       {5,4,7,6,1,0,3,2},
+       {6,5,4,7,2,1,0,3},
+       {7,6,5,4,3,2,1,0}
+      };
+  
+      // So now we need to work out which node needs to be mapped to the 
+      // zero node. 
+      // The trick is that the node that should be the local zero node has
+      // three edges coming into it.
+      for(unsigned int i=0;i<numelems;++i)
+       {
+         Cell& the_cell = the_mesh.cell_list[i];
+    
+         // This stores whether the global oriented edge points in the same
+         // direction as it's local edge on the current cube. (for each edge on
+         // the curent cube)
+         int local_edge_orientation[12];
+         for(unsigned int j=0;j<12;++j)
+           {
+             // get the global edge
+             const Edge& the_edge = the_mesh.edge_list[the_cell.edges[j]];
+             // All edges should be oriented at this stage..
+             Assert(the_edge.orientation_flag!=0,GridOrientError("Unoriented edge encountered"));
+             // calculate whether it points the right way (1) or not (-1)
+             local_edge_orientation[j]=(the_cell.local_orientation_flags[j]*the_edge.orientation_flag);
+           }
+
+         // Here the number of incoming edges is tallied for each node.
+         int perm_num=-1;
+         for(int node_num=0;node_num<8;++node_num)
+           {
+             //The local edge numbers coming into the node
+             int iedg0 = the_mesh.info.edge_to_node[node_num][0];
+             int iedg1 = the_mesh.info.edge_to_node[node_num][1];
+             int iedg2 = the_mesh.info.edge_to_node[node_num][2];
+
+             //The local orientation of the edge coming into the node.
+             int isign0 = the_mesh.info.edge_to_node_orient[node_num][0];
+             int isign1 = the_mesh.info.edge_to_node_orient[node_num][1];
+             int isign2 = the_mesh.info.edge_to_node_orient[node_num][2];
+
+             //Add one to the total for each edge pointing in 
+             int Total  = ((local_edge_orientation[iedg0]*isign0==1)?1:0)
+               + ((local_edge_orientation[iedg1]*isign1==1)?1:0)
+               + ((local_edge_orientation[iedg2]*isign2==1)?1:0);
+      
+             //      std::cout<<"TOTAL : "<<Total<<std::endl;
+      
+             if (Total==3) 
+               {
+                 Assert(perm_num==-1, GridOrientError("More than one node with 3 incoming edges found in curent hex.")); 
+                 perm_num=node_num;
+               }
+           }
+         // We should now have a valid permutation number
+         Assert(perm_num!=-1, GridOrientError("No node having 3 incoming edges found in curent hex.")); 
+
+         // So use the apropriate rotation to get the new cube
+         int temp[8];
+         for(int i=0;i<8;++i)
+           {
+             temp[i]=the_cell.nodes[CubePermutations[perm_num][i]];
+           }
+         for(int i=0;i<8;++i)
+           {
+             the_cell.nodes[i]=temp[i];
+           }
+       }
+    };
+
+  }
+}
+
+
+
+void GridReordering<3>::reorder_cells(std::vector<CellData<3> >& incubes, std::vector<Point<3> > * node_vec_ptr)
+{
+
+  Assert(incubes.size()!=0,
+        ExcMessage("List of elements to orient was of zero length"));
+  
+  // This keeps track of all the local element conectivity information.
+  // e.g. what edges come into which nodes and with which orientation
+  internal::GridReordering3d::DealElemInfo deal_info; 
+  
+  //This does the real work
+  internal::GridReordering3d::Orienter orienter;
+
+  // This is the internal store for all global connectivity information
+  // it starts prety much empty.
+  internal::GridReordering3d::Mesh the_mesh(deal_info);
+
+  if(node_vec_ptr!=NULL)
+    {
+      the_mesh.node_list=*node_vec_ptr;
+    }
+  
+  //Copy the cells into our own internal data format.
+  const unsigned int numelems=incubes.size();
+  for(unsigned int i =0 ; i<numelems; ++i)
+    {
+      internal::GridReordering3d::Cell the_cell;
+      for(unsigned int j=0;j<8;j++)
+       {
+         the_cell.nodes[j]=incubes[i].vertices[j];
+       }
+      the_mesh.cell_list.push_back(the_cell);
+    }
+  
+  // Build the conectivity information
+  // This fills in the conectivity information in the internal structure
+  build_mesh(the_mesh);
+
+  //Orient the mesh
+  orienter.orient_edges(the_mesh);
+
+  // Now we have a bunch of oriented edges int the structure
+  // we only have to turn the cubes so thy match the edge orientation.
+
+  orienter.orient_cubes(the_mesh);
+
+  // Copy the elements from our internal structure back into 
+  // their original location.
+  
+  for(unsigned int i =0 ; i<numelems; ++i)
+    {
+      internal::GridReordering3d::Cell& the_cell = the_mesh.cell_list[i];
+      for(unsigned int j=0;j<8;j++)
+       {
+         incubes[i].vertices[j]=the_cell.nodes[j];
+       }
+    }
+}
+
+
+       
+
+#endif // deal_II_dimension == 3
+

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.