]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Merge Michaels new reorientation scheme from the respective branch.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 21 Feb 2003 20:36:30 +0000 (20:36 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 21 Feb 2003 20:36:30 +0000 (20:36 +0000)
git-svn-id: https://svn.dealii.org/trunk@7219 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_reordering.h
deal.II/deal.II/source/grid/grid_reordering.cc

index 1201a228cbc33fdecf1e053fc52f2a0f43bd3a3f..df04870f5432d5f83654a90ee32ab212aba9093c 100644 (file)
@@ -35,39 +35,6 @@ namespace internal
   };
 
 
-/**
- * Class declaring some dimension dependent numbers which are needed
- * for the grid reordering class. This is the specialization for the
- * 2d case.
- *
- * @author Wolfgang Bangerth, 2000
- */
-  template <>
-  class GridReorderingInfo<2>
-  {
-    public:
-                                      /**
-                                       * Number of possible valid
-                                       * orientations of a cell. They
-                                       * are the state in which it was
-                                       * delivered and three possible
-                                       * rotations in counter-clockwise
-                                       * sense, thus a total of four.
-                                       */
-      static const unsigned int rotational_states_of_cells = 4;
-
-                                      /**
-                                       * Number of possible
-                                       * orientations of a face in
-                                       * 2d. It is the face and the
-                                       * face with vertices exchanged,
-                                       * thus two.
-                                       */
-      static const unsigned int rotational_states_of_faces = 2;
-  };
-
-
-
 /**
  * Class declaring some dimension dependent numbers which are needed
  * for the grid reordering class. This is the specialization for the
@@ -110,7 +77,341 @@ namespace internal
  
 }
 
+
+
+
+
+namespace internal
+{
+/**
+ * Implement the algorithm described in the documentation of the
+ * GridReordering<2> class.
+ *
+ * @author Michael Anderson, 2003
+ */
+  namespace GridReordering2d
+  {
+/**
+ * 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;
+  
+       unsigned int v0;
+       unsigned int v1;
+       unsigned int Q0;
+       unsigned int Q1;
+
+                                        /**
+                                         * Local side numbers on quads 0 and 1.
+                                         */
+       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
+
+
 
 
 /**
@@ -171,7 +472,7 @@ namespace internal
  *   |   |   |
  *   o---o---o
  * @end{verbatim}
- * (The reader is aked to try to find a conforming choice of line
+ * (The reader is asked to try to find a conforming choice of line
  * directions; it will soon be obvious that there can't exists such a
  * thing, even if we allow that there might be cells with clockwise
  * and counterclockwise orientation of the lines at the same time.)
@@ -191,13 +492,29 @@ namespace internal
  * The purpose of this class is now to find an ordering for a given
  * set of cells such that the generated triangulation satisfies all
  * the requirements stated above. To this end, we will first show some
- * examples why this is a difficult problem, and then develop an
- * algorithm that finds such a reordering. Note that the algorithm
+ * examples why this is a difficult problem, and then develop
+ * algorithms that finds such a reordering. Note that the algorithm
  * operates on a set of @ref{CellData} objects that are used to
  * describe a mesh to the triangulation class. These objects are, for
  * example, generated by the @ref{GridIn} class, when reading in grids
  * from input files.
  *
+ * As a last question for this first section: is it guaranteed that
+ * such orientations of faces always exist for a given subdivision of
+ * a domain into cells? The linear complexity algorithm described
+ * below for 2d also proves that the answer is yes for 2d. For 3d, the
+ * answer is no (which also underlines that using such orientations
+ * might be an -- unfortunately uncurable -- misfeature of deal.II). A
+ * simple counter-example in 3d illustrates this: take a string of 3d
+ * cells and bend it together to a torus. Since opposing lines in a
+ * cell need to have the same direction, there is a simple ordering
+ * for them, for example all lines radially outward, tangentially
+ * clockwise, and axially upward. However, if before joining the two
+ * ends of the string of cells, the string is twisted by 180 degrees,
+ * then no such orientation is possible any more, as can easily be
+ * checked. In effect, some meshes cannot be used in deal.II,
+ * unfortunately.
+ *
  *
  * @sect3{Examples of problems}
  *
@@ -305,8 +622,8 @@ namespace internal
  * requirements of deal.II triangulations are met.
  *
  * These two examples demonstrate that if we have added a certain
- * number of cells in some oeirntation of faces and can't add the next
- * one without introducingfaces that had already been added in another
+ * number of cells in some orientation of faces and can't add the next
+ * one without introducing faces that had already been added in another
  * direction, then it might not be sufficient to only rotate cells in
  * the neighborhood of the the cell that we failed to add. It might be
  * necessary to go back a long way and rotate cells that have been
@@ -325,8 +642,8 @@ namespace internal
  * rotated cell 1, then we would have to rotate the cells 1 through
  * N-1 as well).
  *
- * The only solution to this problem seems to be the following: if
- * cell N can't be added, the try to rotate cell N-1. If we can't
+ * A brute force approach to this problem is the following: if
+ * cell N can't be added, then try to rotate cell N-1. If we can't
  * rotate cell N-1 any more, then try to rotate cell N-2 and try to
  * add cell N with all orientations of cell N-1. And so
  * on. Algorithmically, we can visualize this by a tree structure,
@@ -361,7 +678,7 @@ namespace internal
  * that has already been added, then there are already only two
  * possible orientations left, so the total number of checks we have
  * to make until we find a valid way is significantly smaller than
- * @p{4**N}. However, an algorithm is still exponential in time and
+ * @p{4**N}. However, the algorithm is still exponential in time and
  * linear in memory (we only have to store the information for the
  * present path in form of a stack of orientations of cells that have
  * already been added).
@@ -371,26 +688,162 @@ namespace internal
  * very first cells there to find a way to add all cells in a
  * consistent fashion.
  *
- * This discouraging situation is geatly improved by the fact that we
- * can find an algorithm that in practice is usually only roughly
- * linear in time and memory. We will describe this algorithm in the
- * following.
+ * This discouraging situation is greatly improved by the fact that we
+ * have an alternative algorithm for 2d that is always linear in
+ * runtime (discovered and implemented by Michael Anderson of TICAM,
+ * University of Texas, in 2003), and that for 3d we can find an
+ * algorithm that in practice is usually only roughly linear in time
+ * and memory. We will describe these algorithms in the following.
+ *
+ *
+ * @sect3{The 2d linear complexity algorithm}
+ *
+ * The algorithm uses the fact that opposite faces of a cell need to
+ * have the same orientation. So you start with one arbitrary line,
+ * choose an orientation. Then the orientation of the opposite face is
+ * already fixed. Then go to the two cells across the two faces we
+ * have fixed: for them, one face is fixed, so we can also fix the
+ * opposite face. Go on with doing so. Eventually, we have done this
+ * for a string of cells. Then take one of the non-fixed faces of a
+ * cell which has already two fixed faces and do all this again.
+ *
+ * In more detail, the algorithm is best illustrated using an
+ * example. We consider the mesh below:
+ * @begin{verbatim}
+ *   9------10-------11
+ *   |      |        /|
+ *   |      |       / |
+ *   |      |      /  |
+ *   6------7-----8   |
+ *   |      |     |   |
+ *   |      |     |   |
+ *   |      |     |   |
+ *   3------4-----5   |
+ *   |      |      \  |
+ *   |      |       \ |
+ *   |      |        \|
+ *   0------1---------2
+ * @end{verbatim}
+ * First a cell is chosen ( (0,1,4,3) in this case). A single side of the cell
+ * is oriented arbitrarily (3->4). This choice of orientation is then propogated
+ * through the mesh, across sides and elements. (0->1), (6->7) and (9->10).
+ * The involves edge-hopping and face hopping, giving a path through the mesh
+ * shown in dots.
+ * @begin{verbatim}
+ *   9-->--10-------11
+ *   |  .  |        /|
+ *   |  .  |       / |
+ *   |  .  |      /  |
+ *   6-->--7-----8   |
+ *   |  .  |     |   |
+ *   |  .  |     |   |
+ *   |  .  |     |   |
+ *   3-->--4-----5   |
+ *   |  .  |      \  |
+ *   |  X  |       \ |
+ *   |  .  |        \|
+ *   0-->--1---------2
+ * @end{verbatim}
+ * This is then repeated for the other sides of the chosen element, orienting
+ * more sides of the mesh.
+ * @begin{verbatim}
+ *   9-->--10-------11
+ *   |     |        /|
+ *   v.....v.......V |
+ *   |     |      /. |
+ *   6-->--7-----8 . |
+ *   |     |     | . |
+ *   |     |     | . |
+ *   |     |     | . |
+ *   3-->--4-----5 . |
+ *   |     |      \. |
+ *   ^..X..^.......^ |
+ *   |     |        \|
+ *   0-->--1---------2
+ * @end{verbatim}
+ * Once an element has been completely oriented it need not be considered
+ * further. These elements are filled with o's in the diagrams. We then move
+ * to the next element.
+ * @begin{verbatim}
+ *   9-->--10->-----11
+ *   | ooo |  .     /|
+ *   v ooo v  .    V |
+ *   | ooo |  .   /  |
+ *   6-->--7-->--8   |
+ *   |     |  .  |   |
+ *   |     |  .  |   |
+ *   |     |  .  |   |
+ *   3-->--4-->--5   |
+ *   | ooo |  .   \  |
+ *   ^ ooo ^  X    ^ |
+ *   | ooo |  .     \|
+ *   0-->--1-->------2
+ * @end{verbatim}
+ * Repeating this gives
+ * @begin{verbatim}
+ *   9-->--10->-----11
+ *   | ooo | oooooo /|
+ *   v ooo v ooooo V |
+ *   | ooo | oooo /  |
+ *   6-->--7-->--8   |
+ *   |     |     |   |
+ *   ^.....^..X..^...^
+ *   |     |     |   |
+ *   3-->--4-->--5   |
+ *   | ooo | oooo \  |
+ *   ^ ooo ^ ooooo ^ |
+ *   | ooo | oooooo \|
+ *   0-->--1-->------2
+ * @end{verbatim}
+ * and the final oriented mesh is
+ * @begin{verbatim}
+ *   9-->--10->-----11
+ *   |     |        /|
+ *   v     v       V |
+ *   |     |      /  |
+ *   6-->--7-->--8   |
+ *   |     |     |   |
+ *   ^     ^     ^   ^
+ *   |     |     |   |
+ *   3-->--4-->--5   |
+ *   |     |      \  |
+ *   ^     ^       ^ |
+ *   |     |        \|
+ *   0-->--1-->-------2
+ * @end{verbatim}
+ * It is obvious that this algorithm has linear run-time, since it
+ * only ever touches each face exactly once.
+ *
+ * The algorithm just described is implemented in a specialization of
+ * this class for the 2d case. Note that in principle, it should be
+ * possible to extend this algorithm to 3d as well, using sheets
+ * instead of strings of cells to work on. If a grid is reorientable,
+ * then such an algorithm should be able to do so in linear time; if
+ * it is not orientable, then it should abort in linear time as
+ * well. However, this has not yet been implemented. Rather, we use
+ * the backtracking and branch pruning algorithm for 3d described
+ * below; this algorithm predates the 2d linear complexity algorithm
+ * and was initially also used in 2d.
+ *
+ *
+ * @sect3{The 3d branch reshuffling and pruning algorithm}
  *
  * The first observation is that although there are counterexamples,
- * problems are usually local. For example, in the second example, if
- * we had numbered the cells in a way that neighboring cells have
- * similar cell numbers, then the amount pf backtracking needed is
- * greatly reduced. Therefore, in the implementation of the algorithm,
- * the first step is to renumber the cells in a Cuthill-McKee fashion:
- * start with the cell with the least number of neighbors and assign
- * to it the cell number zero. Then find all neighbors of this cell
- * and assign to them consecutive further numbers. Then find their
- * neighbors that have not yet been numbered and assign to them
- * numbers, and so on. Graphically, this represents finding zones of
- * cells consecutively further away from the initial cells and number
- * them in this front-marching way. This already greatly improves
- * locality of problems and consequently reduced the necessary amount
- * of backtracking.
+ * problems are usually local. For example, in the second example
+ * mentioned above, if we had numbered the cells in a way that
+ * neighboring cells have similar cell numbers, then the amount of
+ * backtracking needed is greatly reduced. Therefore, in the
+ * implementation of the algorithm, the first step is to renumber the
+ * cells in a Cuthill-McKee fashion: start with the cell with the
+ * least number of neighbors and assign to it the cell number
+ * zero. Then find all neighbors of this cell and assign to them
+ * consecutive further numbers. Then find their neighbors that have
+ * not yet been numbered and assign to them numbers, and so
+ * on. Graphically, this represents finding zones of cells
+ * consecutively further away from the initial cells and number them
+ * in this front-marching way. This already greatly improves locality
+ * of problems and consequently reduced the necessary amount of
+ * backtracking.
  *
  * The second point is that we can use some methods to prune the tree,
  * which usually lead to a valid orientation of all cells very
@@ -408,13 +861,15 @@ namespace internal
  * neighbor of N with the largest cell index and which has already
  * been added.
  *
- * Unfortunately, this method can fail to yield a valid path through the
- * tree if not applied with care. Consider the following situation,
- * initially extracted from a mesh of 950 cells generated
+ * Unfortunately, this method can fail to yield a valid path through
+ * the tree if not applied with care. Consider the following
+ * situation, initially extracted from a mesh of 950 cells generated
  * automatically by the program BAMG (this program usually generates
  * meshes that are quite badly balanced, often have many -- sometimes
  * 10 or more -- neighbors of one vertex, and exposed several problems
- * in the initial algorithm):
+ * in the initial algorithm; note also that the example is in 2d where
+ * we now have the much better algorithm described above, but the same
+ * observations also apply to 3d):
  * @begin{verbatim}
  * 13----------14----15
  * | \         |     |
@@ -518,7 +973,7 @@ namespace internal
  * the triangulation from this data using the
  * @ref{Triangulation}@p{<dim>::create_triangulation} function.
  *
- * @author Wolfgang Bangerth, 2000
+ * @author Wolfgang Bangerth, 2000, Michael Anderson 2003
  */
 template <int dim>
 class GridReordering
@@ -946,14 +1401,51 @@ class GridReordering
 
 
 
-/* -------------- declaration of explicit specializations ------------- */
+/**
+ * This is the specialization of the general template for 1d. In 1d,
+ * there is actually nothing to be done.
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+template <>
+class GridReordering<1>
+{
+  public:
+                                    /**
+                                     * Do nothing, since in 1d no
+                                     * reordering is necessary.
+                                     */
+    static void reorder_cells (const std::vector<CellData<1> > &);
+};
+
+
 
+/**
+ * This specialization of the general template implements the
+ * 2d-algorithm described in the documentation of the general
+ * template.
+ *
+ * @author Michael Anderson, 2003
+ */
 template <>
-void GridReordering<2>::Cell::insert_faces (std::map<Face,FaceData> &global_faces);
+class GridReordering<2>
+{
+  public:
+                                    /**
+                                     *  This is the main function,
+                                     *  doing what is announced in
+                                     *  the general documentation of
+                                     *  this class.
+                                     */
+    static void reorder_cells (std::vector<CellData<2> > &original_cells);
+};
+
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
 template <>
 void GridReordering<3>::Cell::insert_faces (std::map<Face,FaceData> &global_faces);
-template <>
-void GridReordering<1>::reorder_cells (std::vector<CellData<1> > &);
 
 
 #endif
index bbb6ff1eafac645858303c0a799d267bc04b789c..64713dfd69baf709cfbde220b68cffb21ad79fa3 100644 (file)
@@ -11,7 +11,6 @@
 //
 //----------------------------  grid_reordering.cc  ---------------------------
 
-
 #include <base/thread_management.h>
 #include <grid/grid_reordering.h>
 
 namespace internal
 {
 // static variables
-#if deal_II_dimension == 2
-  const unsigned int GridReorderingInfo<2>::rotational_states_of_cells;
-  const unsigned int GridReorderingInfo<2>::rotational_states_of_faces;
-#endif
-
 #if deal_II_dimension == 3
   const unsigned int GridReorderingInfo<3>::rotational_states_of_cells;
   const unsigned int GridReorderingInfo<3>::rotational_states_of_faces;
@@ -42,6 +36,19 @@ const unsigned int GridReordering<dim>::FaceData::invalid_adjacent_cell;
 
 
 
+#if deal_II_dimension == 1
+
+void GridReordering<1>::reorder_cells (const std::vector<CellData<1> > &)
+{
+                                  // there should not be much to do
+                                  // in 1d...
+}
+
+#endif
+
+
+
+#if deal_II_dimension == 3
 
 template <int dim>
 GridReordering<dim>::Cell::Cell () :
@@ -99,122 +106,6 @@ GridReordering<dim>::Cell::insert_faces (std::map<Face,FaceData> &/*global_faces
 }
 
 
-#if deal_II_dimension == 2
-
-template <>
-void
-GridReordering<2>::Cell::insert_faces (std::map<Face,FaceData> &global_faces)
-{
-  const unsigned int dim = 2;
-
-                                  // first compute index numbers for
-                                  // the faces in usual order as
-                                  // defined by the order of vertices
-                                  // in the cell object
-  Face new_faces[GeometryInfo<dim>::faces_per_cell]
-    = { { { this->vertices[0], this->vertices[1] } },
-        { { this->vertices[1], this->vertices[2] } },
-       { { this->vertices[3], this->vertices[2] } },
-       { { this->vertices[0], this->vertices[3] } } };
-
-                                  // then insert them into the global
-                                  // list and store iterators to
-                                  // them. note that if the face
-                                  // already exists, then the stored
-                                  // data is not touched.
-  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-    faces[0][face] = global_faces.insert (std::make_pair(new_faces[face],
-                                                        FaceData())).first;
-
-
-                                  // then for each of the faces also
-                                  // insert the reverse form and
-                                  // store pointers to them. note
-                                  // that the rotational state in
-                                  // which all faces are reverted is
-                                  // `2'
-  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-    {
-      std::swap (new_faces[face].vertices[0],
-                new_faces[face].vertices[1]);
-      faces[2][face] = global_faces.insert (std::make_pair(new_faces[face],
-                                                          FaceData())).first;
-    };
-
-                                  // then finally fill in rotational
-                                  // states 1 and 3 of the cell. the
-                                  // faces of these states can be
-                                  // obtained from states 0 and 2
-  faces[1][0] = faces[2][0];
-  faces[1][1] = faces[0][1];
-  faces[1][2] = faces[2][2];
-  faces[1][3] = faces[0][3];
-  
-  faces[3][0] = faces[0][0];
-  faces[3][1] = faces[2][1];
-  faces[3][2] = faces[0][2];
-  faces[3][3] = faces[2][3];
-  
-
-                                  // finally fill the crosslink and
-                                  // other fields of the new
-                                  // entries. note that since
-                                  // rotational states 0 and 2 of the
-                                  // cell are exactly reverted, we
-                                  // only have to operate on the face
-                                  // pointers of these two states to
-                                  // reach all possible faces and
-                                  // permutations thereof
-  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-    {
-      if (faces[0][face]->second.adjacent_cells[0] ==
-         FaceData::invalid_adjacent_cell)
-       {
-                                          // face had not been
-                                          // inserted by previous
-                                          // cells, since first
-                                          // adjacent cell is still
-                                          // untouched. provide
-                                          // xlinks to rotated faces
-         faces[0][face]->second.reverse_faces[0] = faces[2][face];
-         faces[2][face]->second.reverse_faces[0] = faces[0][face];
-
-                                          // and insert this cell as
-                                          // adjacent_cell of the faces
-         faces[0][face]->second.adjacent_cells[0] = cell_no;
-         faces[2][face]->second.adjacent_cells[0] = cell_no;
-       }
-      else
-       {
-                                          // face had already been
-                                          // inserted. make sure that
-                                          // it was in the same way:
-         Assert (faces[0][face]->second.reverse_faces[0] == faces[2][face],
-                 ExcInternalError());    
-         Assert (faces[2][face]->second.reverse_faces[0] == faces[0][face],
-                 ExcInternalError());
-
-                                          // now insert ourselves as
-                                          // second
-                                          // adjacent_cell. the
-                                          // respective slots must
-                                          // necessarily be empty
-                                          // still
-         Assert (faces[0][face]->second.adjacent_cells[1] ==
-                 FaceData::invalid_adjacent_cell,
-                 ExcInternalError());
-         Assert (faces[2][face]->second.adjacent_cells[1] ==
-                 FaceData::invalid_adjacent_cell,
-                 ExcInternalError());
-         faces[0][face]->second.adjacent_cells[1] = cell_no;
-         faces[2][face]->second.adjacent_cells[1] = cell_no;
-       };
-    };
-}
-
-#endif
-
-#if deal_II_dimension == 3
 
 template <>
 void
@@ -543,7 +434,6 @@ GridReordering<3>::Cell::insert_faces (std::map<Face,FaceData> &global_faces)
                         [cell_orientation_faces[rot][face].first];
 }
 
-#endif
 
 
 template <int dim>
@@ -1293,25 +1183,466 @@ void GridReordering<dim>::reorder_cells (std::vector<CellData<dim> > &original_c
     Assert (i->second.use_count == 0, ExcInternalError());
 }
 
+#endif // deal_II_dimension == 3
 
+#if deal_II_dimension == 2
 
-#if deal_II_dimension == 1
+namespace internal
+{
+  namespace GridReordering2d
+  {
+// -- Definition Of conectivity information --
+    const int ConnectGlobals::EdgeToNode[4][2]=
+    { {0,1},{1,2},{2,3},{3,0} };
 
-template <>
-void GridReordering<1>::reorder_cells (std::vector<CellData<1> > &)
+    const int ConnectGlobals::NodeToEdge[4][2]=
+    { {3,0},{0,1},{1,2},{2,3} };
+
+    const int ConnectGlobals::DefaultOrientation[4][2]=
+    {{0,1},{1,2},{3,2},{0,3}};
+
+
+
+
+    struct MSide::SideRectify : public std::unary_function<MSide,void>
+    {
+       void operator() (MSide &s) const
+         {
+           if (s.v0>s.v1)
+             std::swap (s.v0, s.v1);
+         }         
+    };
+
+
+    struct MSide::SideSortLess : public std::binary_function<MSide,MSide,bool>
+    {
+       bool operator()(const MSide &s1, const MSide &s2) const
+         {
+           int s1vmin,s1vmax;
+           int s2vmin,s2vmax;
+           if (s1.v0<s1.v1)
+             {
+               s1vmin=s1.v0;
+               s1vmax=s1.v1;
+             }
+           else
+             {
+               s1vmin=s1.v1;
+               s1vmax=s1.v0;
+             }
+           if (s2.v0<s2.v1)
+             {
+               s2vmin=s2.v0;
+               s2vmax=s2.v1;
+             }
+           else
+             {
+               s2vmin=s2.v1;
+               s2vmax=s2.v0;
+             }
+
+           if(s1vmin<s2vmin)
+             return true;
+           if(s1vmin>s2vmin)
+             return false;
+           return s1vmax<s2vmax;
+         }
+    };
+    
+
+/**
+ * Returns an MSide corresponding to the
+ * specified side of a deal.II CellData<2> object.
+ */
+    MSide quadside(const CellData<2> &q, unsigned int i)
+    {
+      Assert (i<4, ExcInternalError());
+      return MSide(q.vertices[ConnectGlobals::EdgeToNode[i][0]],
+                  q.vertices[ConnectGlobals::EdgeToNode[i][1]]);
+    }
+
+
+/**
+ * Wrapper class for the quadside() function
+ */
+    struct QuadSide: public std::binary_function<CellData<2>,int,MSide>
+    {
+       MSide operator()(const CellData<2>& q, int i) const
+         {
+           return quadside(q,i);
+         }
+    };
+    
+
+    MQuad::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)
+                   :
+                   original_cell_data (cd)
+    {
+      v[0]=v0;
+      v[1]=v1;
+      v[2]=v2;
+      v[3]=v3;
+      side[0]=s0;
+      side[1]=s1;
+      side[2]=s2;
+      side[3]=s3;
+    }
+
+
+    MSide::MSide (const unsigned int initv0,
+                 const unsigned int initv1)
+                   :
+                   v0(initv0), v1(initv1),
+                   Q0(static_cast<unsigned int>(-1)),Q1(static_cast<unsigned int>(-1)),
+                   lsn0(static_cast<unsigned int>(-1)),lsn1(static_cast<unsigned int>(-1)),
+                   Oriented(false)
+    {};
+
+    
+    
+    bool
+    MSide::operator== (const MSide& s2) const
+    {
+      if ((v0==s2.v0)&&(v1==s2.v1)) {return true;}
+      if ((v0==s2.v1)&&(v1==s2.v0)) {return true;}
+      return false;
+    }
+
+    
+    struct MQuad::MakeQuad : public std::binary_function<CellData<2>,
+                                                        std::vector<MSide>,
+                                                        MQuad>
+    {
+       MQuad operator()(const CellData<2> &q,
+                        const std::vector<MSide> &elist) const
+         {
+                                            //Assumes that the sides
+                                            //are in the vector.. Bad
+                                            //things will happen if
+                                            //they are not!
+           return MQuad(q.vertices[0],q.vertices[1], q.vertices[2], q.vertices[3],
+                        std::distance(elist.begin(),
+                                      std::lower_bound(elist.begin(), elist.end(),
+                                                       quadside(q,0),
+                                                       MSide::SideSortLess() )),
+                        std::distance(elist.begin(),
+                                      std::lower_bound(elist.begin(), elist.end(),
+                                                       quadside(q,1),
+                                                       MSide::SideSortLess() )),
+                        std::distance(elist.begin(),
+                                      std::lower_bound(elist.begin(), elist.end(),
+                                                       quadside(q,2),
+                                                       MSide::SideSortLess() )),
+                        std::distance(elist.begin(),
+                                      std::lower_bound(elist.begin(), elist.end(),
+                                                       quadside(q,3),
+                                                       MSide::SideSortLess() )),
+                        q);
+         }
+           
+    };
+
+
+    
+    void
+    GridReordering::reorient(std::vector<CellData<2> > &quads)
+    {
+      build_graph(quads);
+      orient();
+      get_quads(quads);
+    }
+
+
+    void
+    GridReordering::build_graph (const std::vector<CellData<2> > &inquads)
+    {
+                                      //Reserve some space 
+      sides.reserve(4*inquads.size());
+      mquads.reserve(inquads.size());
+  
+                                      //Insert all the sides into the side vector
+      for (int i=0;i<4;++i)
+       {
+         std::transform(inquads.begin(),inquads.end(),
+                        std::back_inserter(sides), std::bind2nd(QuadSide(),i));
+       }
+  
+                                      //Change each edge so that v0<v1
+      std::for_each(sides.begin(),sides.end(),
+                   MSide::SideRectify() );
+  
+                                      //Sort them by Sidevertices.
+      std::sort(sides.begin(),sides.end(),
+               MSide::SideSortLess());
+  
+                                      //Remove duplicates 
+      sides.erase(std::unique(sides.begin(),sides.end()),
+                 sides.end());
+
+                                      // Swap trick to shrink the
+                                      // side vector
+      std::vector<MSide>(sides).swap(sides);
+  
+                                      //Assigns the correct sides to
+                                      //each quads
+      transform(inquads.begin(),inquads.end(), back_inserter(mquads),
+               std::bind2nd(MQuad::MakeQuad(),sides) );
+  
+                                      // Assign the quads to their sides also.
+      int qctr=0;
+      for(std::vector<MQuad>::iterator it=mquads.begin(); it!=mquads.end(); ++it)
+       {
+         for(int i=0;i<4;++i)
+           {
+             MSide &ss =sides[(*it).side[i]];
+             if(ss.Q0==static_cast<unsigned int>(-1))
+               {
+                 ss.Q0=qctr;
+                 ss.lsn0=i;
+               }
+             else if (ss.Q1==static_cast<unsigned int>(-1))
+               {
+                 ss.Q1=qctr;
+                 ss.lsn1=i;
+               }
+             else
+               {
+                 exit(0);
+               }
+           }
+         qctr++;
+       }
+    }
+
+
+    void GridReordering::orient()
+    {
+                                      // do what the comment in the
+                                      // class declaration says
+      unsigned int qnum=0;
+      while(get_unoriented_quad(qnum))
+       {
+         unsigned int lsn=0;
+         while(get_unoriented_side(qnum,lsn))
+           {
+             orient_side(qnum,lsn);
+             unsigned int qqnum=qnum;
+             while(side_hop(qqnum,lsn))
+               {
+                                                  // switch this face
+                 lsn = (lsn+2)%4;
+                 if (!is_oriented_side(qqnum,lsn))
+                   orient_side(qqnum,lsn);
+                 else
+                                                    //We've found a
+                                                    //cycle.. and
+                                                    //oriented all
+                                                    //quads in it.
+                   break;
+               }
+           }
+       }
+    }
+
+
+    void
+    GridReordering::orient_side(const unsigned int quadnum,
+                               const unsigned int localsidenum)
+    {
+      MQuad &quad = mquads[quadnum];
+      int op_side_l = (localsidenum+2)%4;
+      MSide &side = sides[mquads[quadnum].side[localsidenum]];
+      const MSide &op_side =sides[mquads[quadnum].side[op_side_l]]; 
+  
+                                      //is the opposite side oriented?    
+      if (op_side.Oriented)
+       {
+                                          //YES - Make the orientations match
+                                          //Is op side in default orientation?
+         if (op_side.v0==quad.v[ConnectGlobals::DefaultOrientation[op_side_l][0]])
+           {
+                                              //YES
+             side.v0=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][0]];
+             side.v1=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][1]];
+           }
+         else
+           {
+                                              //NO, its reversed
+             side.v0=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][1]];
+             side.v1=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][0]];
+           }
+       }
+      else
+       {
+                                          //NO
+                                          //Just use the default orientation      
+         side.v0=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][0]];
+         side.v1=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][1]];
+       }
+      side.Oriented=true;  
+    }
+
+
+
+    bool
+    GridReordering::is_fully_oriented_quad(const unsigned int quadnum) const
+    {
+      return (
+       (sides[mquads[quadnum].side[0]].Oriented)&&
+       (sides[mquads[quadnum].side[1]].Oriented)&&
+       (sides[mquads[quadnum].side[2]].Oriented)&&
+       (sides[mquads[quadnum].side[3]].Oriented) 
+      );
+    }
+
+
+
+    bool
+    GridReordering::is_oriented_side(const unsigned int quadnum,
+                                    const unsigned int lsn) const
+    {
+      return (sides[mquads[quadnum].side[lsn]].Oriented);
+    }
+
+
+
+
+    bool
+    GridReordering::get_unoriented_quad(unsigned int &UnOrQLoc) const
+    {
+      while( (UnOrQLoc<mquads.size()) &&
+            is_fully_oriented_quad(UnOrQLoc) )
+       UnOrQLoc++;
+      return (UnOrQLoc!=mquads.size());
+    }
+
+
+
+    bool
+    GridReordering::get_unoriented_side (const unsigned int quadnum,
+                                        unsigned int &lsn) const
+    {
+      const MQuad &mq = mquads[quadnum];
+      if(!sides[mq.side[0]].Oriented)
+       {
+         lsn=0;
+         return true;
+       }
+      if(!sides[mq.side[1]].Oriented)
+       {
+         lsn=1;
+         return true;
+       }
+      if(!sides[mq.side[2]].Oriented)
+       {
+         lsn=2;
+         return true;
+       }
+      if(!sides[mq.side[3]].Oriented)
+       {
+         lsn=3;
+         return true;
+       }
+      return false;
+    }
+
+
+    bool
+    GridReordering::side_hop (unsigned int &qnum, unsigned int &lsn) const
+    {
+      const MQuad &mq=mquads[qnum];
+      const MSide &s = sides[mq.side[lsn]];
+      unsigned int opquad=0;
+      if (s.Q0==qnum)
+       {
+         opquad=s.Q1;
+         lsn =s.lsn1;
+       }
+      else
+       {
+         opquad=s.Q0;
+         lsn=s.lsn0;
+       }
+  
+      if (opquad!=static_cast<unsigned int>(-1))
+       {
+         qnum = opquad;
+         return true;
+       }
+  
+      return false;
+    }
+
+
+    void
+    GridReordering::get_quads (std::vector<CellData<2> > &outquads) const
+    {
+      outquads.clear();
+      outquads.reserve(mquads.size());
+      for(unsigned int qn=0;qn<mquads.size();++qn)
+       {
+                                          // initialize CellData object with
+                                          // previous contents, and the
+                                          // overwrite all the fields that
+                                          // might have changed in the
+                                          // process of rotating things
+         CellData<2> q = mquads[qn].original_cell_data;
+         
+                                          //Are the sides oriented? 
+         assert(is_fully_oriented_quad(qn));
+         bool s[4]; //whether side 1 ,2, 3, 4 are in the default orientation
+         for(int sn=0;sn<4;sn++)
+           {
+             s[sn]=is_side_default_oriented(qn,sn);
+           }
+                                          // Are they oriented in the "deal way"?
+         assert(s[0]==s[2]);
+         assert(s[1]==s[3]);
+                                          // How much we rotate them by.
+         int rotn = 2*(s[0]?1:0)+ ((s[0]^s[1])?1:0);
+
+         for(int i=0;i<4;++i)
+           {
+             q.vertices[(i+rotn)%4]=mquads[qn].v[i];
+           }
+         outquads.push_back(q);
+       }
+
+    }
+
+    bool
+    GridReordering::is_side_default_oriented (const unsigned int qnum,
+                                             const unsigned int lsn) const
+    {
+      return (sides[mquads[qnum].side[lsn]].v0 ==
+             mquads[qnum].v[ConnectGlobals::DefaultOrientation[lsn][0]]);
+    }
+  } // namespace GridReordering2
+} // namespace internal
+
+
+void GridReordering<2>::reorder_cells (std::vector<CellData<2> > &original_cells)
 {
-                                  // there should not be much to do
-                                  // in 1d...
+  internal::GridReordering2d::GridReordering().reorient(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
-#if deal_II_dimension >= 2
+// 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>::

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.