]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Complete the implementation of the mesh orientation algorithm in 3d. 3263/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Thu, 13 Oct 2016 10:23:03 +0000 (04:23 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Wed, 19 Oct 2016 03:01:09 +0000 (21:01 -0600)
This replaces the existing implementation by one that is largely
dimension independent.

include/deal.II/grid/grid_reordering_internal.h [deleted file]
source/grid/grid_reordering.cc

diff --git a/include/deal.II/grid/grid_reordering_internal.h b/include/deal.II/grid/grid_reordering_internal.h
deleted file mode 100644 (file)
index 7c8d8c4..0000000
+++ /dev/null
@@ -1,367 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2003 - 2016 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-#ifndef dealii__grid_reordering_internal_h
-#define dealii__grid_reordering_internal_h
-
-
-#include <deal.II/base/config.h>
-#include <deal.II/grid/tria.h>
-
-#include <vector>
-
-DEAL_II_NAMESPACE_OPEN
-
-
-namespace internal
-{
-  /**
-   * Implement the algorithm described in the documentation of the
-   * GridReordering<3> class.
-   *
-   * @author Michael Anderson, 2003
-   */
-  namespace GridReordering3d
-  {
-    /**
-     * A structure indicating the direction of an edge. In the implementation
-     * file, we define three objects, <tt>unoriented_edge</tt>,
-     * <tt>forward_edge</tt>, and <tt>backward_edge</tt>, that denote whether
-     * an edge has already been oriented, whether it is in standard
-     * orientation, or whether it has reverse direction. The state that each
-     * of these objects encode is stored in the <tt>orientation</tt> member
-     * variable -- we would really need only three such values, which we pick
-     * in the implementation file, and make sure when we compare such objects
-     * that only these three special values are actually used.
-     *
-     * The reason for this way of implementing things is as follows. Usually,
-     * such a property would be implemented as an enum. However, in the
-     * previous implementation, a signed integer was used with unoriented=0,
-     * forward=+1, and backward=-1. A number of operations, such as equality
-     * of ordered edges were mapped to checking whether the product of two
-     * edge orientations equals +1. Such arithmetic isn't always portable and
-     * sometimes flagged when using -ftrapv with gcc. Using this class instead
-     * makes sure that there isn't going to be any arithmetic going on on edge
-     * orientations, just comparisons for equality or inequality.
-     *
-     * @author Wolfgang Bangerth, 2005
-     */
-    struct EdgeOrientation
-    {
-      /**
-       * A value indicating the orientation.
-       */
-      char orientation;
-
-      /**
-       * Comparison operator.
-       */
-      bool operator == (const EdgeOrientation &edge_orientation) const;
-
-      /**
-       * Comparison operator.
-       */
-      bool operator != (const EdgeOrientation &edge_orientation) const;
-    };
-
-    /**
-     * During building the connectivity information we don't need all the
-     * heavy duty information about edges that we will need later. So we can
-     * save memory and time by using a light-weight class for edges. It stores
-     * the two vertices, but no direction, so we make the optimization to
-     * store the vertex number in sorted order to allow for easier comparison
-     * of edge objects.
-     */
-    struct CheapEdge
-    {
-      /**
-       * The first node
-       */
-      const unsigned int node0;
-
-      /**
-       * The second node
-       */
-      const unsigned int node1;
-
-      /**
-       * Constructor. Take the vertex numbers and store them sorted.
-       */
-      CheapEdge (const unsigned int n0,
-                 const unsigned int n1);
-
-      /**
-       * Need a partial ordering for sorting algorithms
-       */
-      bool operator< (const CheapEdge &e2) const;
-    };
-
-
-
-    /**
-     * A connectivity and orientation aware edge class.
-     */
-    struct Edge
-    {
-      /**
-       * Simple constructor
-       */
-      Edge (const unsigned int n0,
-            const unsigned int n1);
-
-      /**
-       * The IDs for the end nodes
-       */
-      unsigned int nodes[2];
-
-      /**
-       * Whether the edge has not already been oriented, points from node 0 to
-       * node 1, or the reverse. The initial state of this flag is unoriented.
-       */
-      EdgeOrientation orientation_flag;
-
-      /**
-       * Used to determine which "sheet" or equivalence class of parallel
-       * edges the edge falls in when oriented. numbers::invalid_unsigned_int
-       * means not yet decided. This is also the default value after
-       * construction. Each edge will later be assigned an index greater than
-       * zero.
-       */
-      unsigned int group;
-
-      /**
-       * Indices of neighboring cubes.
-       */
-      std::vector<unsigned int> neighboring_cubes;
-    };
-
-    /**
-     * A connectivity and orientation aware cell.
-     *
-     * The connectivity of the cell is not contained within. (This was for
-     * flexibility in using deal.II'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: Need to move connectivity information out of cell and into edge.
-     */
-    struct Cell
-    {
-      /**
-       * Default Constructor
-       */
-      Cell ();
-
-      /**
-       * The IDs for each of the edges.
-       */
-      unsigned int edges[GeometryInfo<3>::lines_per_cell];
-
-      /**
-       * The IDs for each of the nodes.
-       */
-      unsigned int nodes[GeometryInfo<3>::vertices_per_cell];
-
-      /**
-       * 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).
-       */
-      EdgeOrientation local_orientation_flags[GeometryInfo<3>::lines_per_cell];
-
-      /**
-       * An internal flag used to determine whether the cell is in the queue
-       * of cells to be oriented in the current sheet.
-       */
-      bool waiting_to_be_processed;
-    };
-
-
-    /**
-     * 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:
-      /**
-       * Default Constructor
-       */
-      Mesh (const std::vector<CellData<3> > &incubes);
-
-      /**
-       * Export the data of this object to the deal.II format that the
-       * Triangulation class wants as input.
-       */
-      void
-      export_to_deal_format (std::vector<CellData<3> > &outcubes) const;
-
-    private:
-      /**
-       * The list of edges
-       */
-      std::vector<Edge> edge_list;
-
-      /**
-       * The list of cells
-       */
-      std::vector<Cell> cell_list;
-
-      /**
-       * Check whether every cell in the mesh is sensible.
-       */
-      void sanity_check() const;
-
-      /**
-       * Given the cell list, build the edge list and all the connectivity
-       * information and other stuff that we will need later.
-       */
-      void build_connectivity ();
-
-      /**
-       * Unimplemented private copy constructor to disable it.
-       */
-      Mesh (const Mesh &);
-
-      /**
-       * Unimplemented private assignment operator to disable it.
-       */
-      Mesh &operator=(const Mesh &);
-
-      /**
-       * Checks that each edge going into a node is correctly set up.
-       */
-      void sanity_check_node (const Cell        &cell,
-                              const unsigned int local_node_num) const;
-
-      /**
-       * Let the orienter access out private fields.
-       */
-      friend class Orienter;
-    };
-
-
-    /**
-     * The class that orients the edges of a triangulation in 3d. The member
-     * variables basically only store the present state of the algorithm.
-     */
-    class Orienter
-    {
-    public:
-      /**
-       * Orient the given mesh. Creates an object of the present type and lets
-       * that toil away at the task.
-       *
-       * This function is the single entry point to the functionality of this
-       * class.
-       *
-       * Return, whether a consistent orientation of lines was possible for
-       * the given mesh.
-       */
-      static
-      bool
-      orient_mesh (std::vector<CellData<3> > &incubes);
-
-    private:
-      /**
-       * Internal representation of the given list of cells, including
-       * connectivity information and the like.
-       */
-      Mesh mesh;
-
-      /**
-       * The cube we're looking at presently.
-       */
-      unsigned int cur_posn;
-
-      /**
-       * We have fully oriented all cubes before this one.
-       */
-      unsigned int marker_cube;
-
-      /**
-       * The index of the sheet or equivalence class we are presently
-       * processing.
-       */
-      unsigned int cur_edge_group;
-
-      /**
-       * Indices of the cells to be processed within the present sheet. If a
-       * cell is being processed presently, it is taken from this list.
-       */
-      std::vector<int> sheet_to_process;
-
-
-      /**
-       * Which edges of the current cell have been oriented during the current
-       * iteration. Is reset when moving on to the next cube.
-       */
-      bool edge_orient_array[12];
-
-      /**
-       * Constructor. Take a list of cells and set up the internal data
-       * structures of the mesh member variable.
-       *
-       * Since it is private, the only entry point of this class is the static
-       * function orient_mesh().
-       */
-      Orienter (const std::vector<CellData<3> > &incubes);
-
-      /**
-       * Orient all the edges of a mesh.
-       *
-       * Return, whether this action was carried out successfully.
-       */
-      bool orient_edges ();
-
-      /**
-       * Given oriented edges, rotate the cubes so that the edges are in
-       * standard direction.
-       */
-      void orient_cubes ();
-
-      bool get_next_unoriented_cube ();
-
-      /**
-       * Return whether the cell with cell number @p cell_num is fully
-       * oriented.
-       */
-      bool is_oriented (const unsigned int cell_num) const;
-
-      bool orient_edges_in_current_cube ();
-      bool orient_edge_set_in_current_cube (const unsigned int edge_set);
-      bool orient_next_unoriented_edge ();
-
-      /**
-       * Return whether the cell is consistently oriented at present (i.e.
-       * only considering those edges that are already oriented. This is a
-       * sanity check that should be called from inside an assert macro.
-       */
-      bool cell_is_consistent (const unsigned int cell_num) const;
-
-
-      void get_adjacent_cubes ();
-      bool get_next_active_cube ();
-    };
-  }  // namespace GridReordering3d
-}  // namespace internal
-
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
index e1e365f751dba5cc8b63fe7d0bc5ba380916e759..254f7e8a3a376b928c45261d8c6171f1f541ac4d 100644 (file)
 
 
 #include <deal.II/grid/grid_reordering.h>
-#include <deal.II/grid/grid_reordering_internal.h>
 #include <deal.II/grid/grid_tools.h>
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/std_cxx11/bind.h>
+#include <deal.II/base/timer.h>
 
 #include <algorithm>
 #include <set>
@@ -33,6 +33,9 @@ namespace internal
 {
   namespace GridReordering2d
   {
+    DeclExceptionMsg (ExcMeshNotOrientable,
+                      "The edges of the mesh are not consistently orientable.");
+
     /**
      * A simple data structure denoting an edge, i.e., the ordered pair
      * of its vertex indices. This is only used in the is_consistent()
@@ -117,8 +120,19 @@ namespace internal
     template <int dim>
     struct ParallelEdges
     {
+      /**
+       * An array that contains the indices of dim edges that can
+       * serve as (arbitrarily chosen) starting points for the
+       * dim sets of parallel edges within each cell.
+       */
       static const unsigned int starter_edges[dim];
-      static const unsigned int parallel_edges[GeometryInfo<dim>::lines_per_cell][(1<<(dim-1)) - 1];
+
+      /**
+       * Number and indices of all of those edges parallel to each of the
+       * edges in a cell.
+       */
+      static const unsigned int n_other_parallel_edges = (1<<(dim-1)) - 1;
+      static const unsigned int parallel_edges[GeometryInfo<dim>::lines_per_cell][n_other_parallel_edges];
     };
 
     template <>
@@ -189,6 +203,10 @@ namespace internal
     class AdjacentCells<2>
     {
     public:
+      /**
+       * An iterator that allows iterating over all cells adjacent
+       * to the edge represented by the current object.
+       */
       typedef const AdjacentCell *const_iterator;
 
       /**
@@ -251,6 +269,19 @@ namespace internal
     };
 
 
+
+    /**
+     * A class that represents all of the cells adjacent to a given edge.
+     * This class corresponds to the 3d case where each edge can have an
+     * arbitrary number of adjacent cells. We represent this as a
+     * std::vector<AdjacentCell>, from which class the current one is
+     * derived and from which it inherits all of its member functions.
+     */
+    template <>
+    class AdjacentCells<3> : public std::vector<AdjacentCell>
+    {};
+
+
     /**
      * A class that describes all of the relevant properties of an
      * edge. For the purpose of what we do here, that includes the
@@ -502,7 +533,7 @@ namespace internal
      * base class.
      */
     template <>
-    class EdgeDeltaSet<3> : std::set<unsigned int>
+    class EdgeDeltaSet<3> : public std::set<unsigned int>
     {};
 
 
@@ -590,14 +621,13 @@ namespace internal
      * (global) parallel to the one identified by the @p cell and
      * within it the one with index @p local_edge.
      */
+    template <int dim>
     void
-    orient_one_set_of_parallel_edges (const std::vector<Cell<2> > &cells,
-                                      std::vector<Edge<2> >       &edges,
-                                      const unsigned int           cell,
-                                      const unsigned int           local_edge)
+    orient_one_set_of_parallel_edges (const std::vector<Cell<dim> > &cells,
+                                      std::vector<Edge<dim> >       &edges,
+                                      const unsigned int             cell,
+                                      const unsigned int             local_edge)
     {
-      const unsigned int dim = 2;
-
       // choose the direction of the first edge. we have free choice
       // here and could simply choose "forward" if that's what pleases
       // us. however, for backward compatibility with the previous
@@ -614,12 +644,18 @@ namespace internal
       // same effect by modifying the rule above to choose the
       // direction of the starting edge of this parallel set
       // *opposite* to what it looks like in the current cell
+      //
+      // this bug only existed in the 2d implementation since there
+      // were different implementations for 2d and 3d. consequently,
+      // only replicate it for the 2d case and be "intuitive" in 3d.
       if (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0]
           ==
           cells[cell].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices (local_edge, 0)])
         // orient initial edge *opposite* to the way it is in the cell
         // (see above for the reason)
-        edges[cells[cell].edge_indices[local_edge]].orientation_status = Edge<dim>::backward;
+        edges[cells[cell].edge_indices[local_edge]].orientation_status = (dim == 2 ?
+            Edge<dim>::backward :
+            Edge<dim>::forward);
       else
         {
           Assert (edges[cells[cell].edge_indices[local_edge]].vertex_indices[0]
@@ -633,7 +669,9 @@ namespace internal
 
           // orient initial edge *opposite* to the way it is in the cell
           // (see above for the reason)
-          edges[cells[cell].edge_indices[local_edge]].orientation_status = Edge<dim>::forward;
+          edges[cells[cell].edge_indices[local_edge]].orientation_status = (dim == 2 ?
+              Edge<dim>::forward :
+              Edge<dim>::backward);
         }
 
       // walk outward from the given edge as described in
@@ -650,14 +688,14 @@ namespace internal
         {
           Delta_k.clear ();
 
-          for (EdgeDeltaSet<dim>::const_iterator delta = Delta_k_minus_1.begin();
+          for (typename EdgeDeltaSet<dim>::const_iterator delta = Delta_k_minus_1.begin();
                delta != Delta_k_minus_1.end(); ++delta)
             {
               Assert (edges[*delta].orientation_status != Edge<dim>::not_oriented,
                       ExcInternalError());
 
               // now go through the cells adjacent to this edge
-              for (AdjacentCells<dim>::const_iterator
+              for (typename AdjacentCells<dim>::const_iterator
                    adjacent_cell = edges[*delta].adjacent_cells.begin();
                    adjacent_cell != edges[*delta].adjacent_cells.end(); ++adjacent_cell)
                 {
@@ -672,48 +710,76 @@ namespace internal
                        edges[*delta].vertex_indices[0]
                        :
                        edges[*delta].vertex_indices[1]);
-                  const unsigned int first_edge_vertex_in_K = cells[K].vertex_indices[GeometryInfo<dim>::face_to_cell_vertices(delta_is_edge_in_K, 0)];
+                  const unsigned int first_edge_vertex_in_K
+                    = cells[K].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(delta_is_edge_in_K, 0)];
                   Assert (first_edge_vertex == first_edge_vertex_in_K
                           ||
-                          first_edge_vertex == cells[K].vertex_indices[GeometryInfo<dim>::face_to_cell_vertices(delta_is_edge_in_K, 1)],
+                          first_edge_vertex == cells[K].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
                           ExcInternalError());
 
-                  // now figure out which direction the opposite edge needs to be into.
-                  const unsigned int opposite_edge
-                    = cells[K].edge_indices[ParallelEdges<2>::parallel_edges[delta_is_edge_in_K][0]];
-                  const unsigned int first_opposite_edge_vertex
-                    =  cells[K].vertex_indices[GeometryInfo<dim>::face_to_cell_vertices(
-                                                 ParallelEdges<dim>::parallel_edges[delta_is_edge_in_K][0],
-                                                 (first_edge_vertex == first_edge_vertex_in_K
-                                                  ?
-                                                  0
-                                                  :
-                                                  1))];
-
-                  const Edge<dim>::OrientationStatus opposite_edge_orientation
-                    = (edges[opposite_edge].vertex_indices[0]
-                       ==
-                       first_opposite_edge_vertex
-                       ?
-                       Edge<dim>::forward
-                       :
-                       Edge<dim>::backward);
-
-                  // see if the opposite edge (there is only one in 2d) has already been
-                  // oriented.
-                  if (edges[opposite_edge].orientation_status == Edge<dim>::not_oriented)
+                  // now figure out which direction the each of the "opposite" edges
+                  // needs to be oriented into.
+                  for (unsigned int o_e=0; o_e<ParallelEdges<dim>::n_other_parallel_edges; ++o_e)
                     {
-                      // the opposite edge is not yet oriented. do orient it and add it to
-                      // Delta_k
-                      edges[opposite_edge].orientation_status = opposite_edge_orientation;
-                      Delta_k.insert (opposite_edge);
-                    }
-                  else
-                    {
-                      // the opposite edge has already been oriented. assert that it is
-                      // consistent with the current one
-                      Assert (edges[opposite_edge].orientation_status == opposite_edge_orientation,
-                              ExcInternalError());
+                      // get the index of the opposite edge and select which its first
+                      // vertex needs to be based on how the current edge is oriented
+                      // in the current cell
+                      const unsigned int opposite_edge
+                        = cells[K].edge_indices[ParallelEdges<dim>::parallel_edges[delta_is_edge_in_K][o_e]];
+                      const unsigned int first_opposite_edge_vertex
+                        =  cells[K].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(
+                                                     ParallelEdges<dim>::parallel_edges[delta_is_edge_in_K][o_e],
+                                                     (first_edge_vertex == first_edge_vertex_in_K
+                                                      ?
+                                                      0
+                                                      :
+                                                      1))];
+
+                      // then determine the orientation of the edge based on
+                      // whether the vertex we want to be the edge's first
+                      // vertex is already the first vertex of the edge, or
+                      // whether it points in the opposite direction
+                      const typename Edge<dim>::OrientationStatus opposite_edge_orientation
+                        = (edges[opposite_edge].vertex_indices[0]
+                           ==
+                           first_opposite_edge_vertex
+                           ?
+                           Edge<dim>::forward
+                           :
+                           Edge<dim>::backward);
+
+                      // see if the opposite edge (there is only one in 2d) has already been
+                      // oriented.
+                      if (edges[opposite_edge].orientation_status == Edge<dim>::not_oriented)
+                        {
+                          // the opposite edge is not yet oriented. do orient it and add it to
+                          // Delta_k
+                          edges[opposite_edge].orientation_status = opposite_edge_orientation;
+                          Delta_k.insert (opposite_edge);
+                        }
+                      else
+                        {
+                          // this opposite edge has already been oriented. it should be
+                          // consistent with the current one in 2d, while in 3d it may in fact
+                          // be mis-oriented, and in that case the mesh will not be
+                          // orientable. indicate this by throwing an exception that we can
+                          // catch further up; this has the advantage that we can propagate
+                          // through a couple of functions without having to do error
+                          // checking and without modifying the 'cells' array that the
+                          // user gave us
+                          if (dim == 2)
+                            {
+                              Assert (edges[opposite_edge].orientation_status == opposite_edge_orientation,
+                                      ExcMeshNotOrientable());
+                            }
+                          else if (dim == 3)
+                            {
+                              if (edges[opposite_edge].orientation_status != opposite_edge_orientation)
+                                throw ExcMeshNotOrientable ();
+                            }
+                          else
+                            Assert (false, ExcNotImplemented());
+                        }
                     }
                 }
             }
@@ -733,54 +799,135 @@ namespace internal
      * system matches the ones of the adjacent edges. Store the
      * rotated order of vertices in <code>raw_cells[cell_index]</code>.
      */
+    template <int dim>
     void
-    rotate_cell (const std::vector<Cell<2> > &cell_list,
-                 const std::vector<Edge<2> > &edge_list,
-                 const unsigned int           cell_index,
-                 std::vector<CellData<2> >   &raw_cells)
+    rotate_cell (const std::vector<Cell<dim> > &cell_list,
+                 const std::vector<Edge<dim> > &edge_list,
+                 const unsigned int             cell_index,
+                 std::vector<CellData<dim> >   &raw_cells)
     {
-      // find the first vertex of the cell. this is the
-      // vertex where two edges originate, so for
-      // each of the four edges record which the
-      // starting vertex is
-      unsigned int starting_vertex_of_edge[4];
-      for (unsigned int e=0; e<4; ++e)
+      // find the first vertex of the cell. this is the vertex where dim edges
+      // originate, so for each of the edges record which the starting vertex is
+      unsigned int starting_vertex_of_edge[GeometryInfo<dim>::lines_per_cell];
+      for (unsigned int e=0; e<GeometryInfo<dim>::lines_per_cell; ++e)
         {
           Assert (edge_list[cell_list[cell_index].edge_indices[e]].orientation_status
-                  != Edge<2>::not_oriented,
+                  != Edge<dim>::not_oriented,
                   ExcInternalError());
-          if (edge_list[cell_list[cell_index].edge_indices[e]].orientation_status == Edge<2>::forward)
+          if (edge_list[cell_list[cell_index].edge_indices[e]].orientation_status == Edge<dim>::forward)
             starting_vertex_of_edge[e] = edge_list[cell_list[cell_index].edge_indices[e]].vertex_indices[0];
           else
             starting_vertex_of_edge[e] = edge_list[cell_list[cell_index].edge_indices[e]].vertex_indices[1];
         }
 
-      // find the vertex number that appears twice. this must either be
-      // the first, second, or third vertex in the list. because edges
-      // zero and one don't share any vertices, and the same for edges
-      // two and three, the possibilities can easily be enumerated
-      unsigned int starting_vertex_of_cell = numbers::invalid_unsigned_int;
-      if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2])
-          ||
-          (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
-        starting_vertex_of_cell = starting_vertex_of_edge[0];
-      else if ((starting_vertex_of_edge[1] == starting_vertex_of_edge[2])
-               ||
-               (starting_vertex_of_edge[1] == starting_vertex_of_edge[3]))
-        starting_vertex_of_cell = starting_vertex_of_edge[1];
-      else
-        Assert (false, ExcInternalError());
+      // find the vertex number that appears dim times. this will then be
+      // the vertex at which we want to locate the origin of the cell's
+      // coordinate system (i.e., vertex 0)
+      unsigned int origin_vertex_of_cell = numbers::invalid_unsigned_int;
+      switch (dim)
+        {
+        case 2:
+        {
+          // in 2d, we can simply enumerate the possibilities where the
+          // origin may be located because edges zero and one don't share
+          // any vertices, and the same for edges two and three
+          if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2])
+              ||
+              (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
+            origin_vertex_of_cell = starting_vertex_of_edge[0];
+          else if ((starting_vertex_of_edge[1] == starting_vertex_of_edge[2])
+                   ||
+                   (starting_vertex_of_edge[1] == starting_vertex_of_edge[3]))
+            origin_vertex_of_cell = starting_vertex_of_edge[1];
+          else
+            Assert (false, ExcInternalError());
+
+          break;
+        }
+
+        case 3:
+        {
+          // one could probably do something similar in 3d, but that seems
+          // more complicated than one wants to write down. just go
+          // through the list of possible starting vertices and check
+          for (origin_vertex_of_cell=0;
+               origin_vertex_of_cell<GeometryInfo<dim>::vertices_per_cell;
+               ++origin_vertex_of_cell)
+            if (std::count (&starting_vertex_of_edge[0],
+                            &starting_vertex_of_edge[0]+GeometryInfo<dim>::lines_per_cell,
+                            cell_list[cell_index].vertex_indices[origin_vertex_of_cell])
+                == dim)
+              break;
+          Assert (origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell,
+                  ExcInternalError());
+
+          break;
+        }
 
-      // now rotate raw_cells[cell_index] until the starting indices match.
-      // take into account the ordering of vertices (not in clockwise
-      // or counter-clockwise sense)
-      while (raw_cells[cell_index].vertices[0] != starting_vertex_of_cell)
+        default:
+          Assert (false, ExcNotImplemented());
+        }
+
+      // now rotate raw_cells[cell_index] in such a way that its orientation
+      // matches that of cell_list[cell_index]
+      switch (dim)
+        {
+        case 2:
+        {
+          // in 2d, we can literally rotate the cell until its origin
+          // matches the one that we have determined above should be
+          // the origin vertex
+          //
+          // when doing a rotation, take into account the ordering of
+          // vertices (not in clockwise or counter-clockwise sense)
+          while (raw_cells[cell_index].vertices[0] != origin_vertex_of_cell)
+            {
+              const unsigned int tmp = raw_cells[cell_index].vertices[0];
+              raw_cells[cell_index].vertices[0] = raw_cells[cell_index].vertices[1];
+              raw_cells[cell_index].vertices[1] = raw_cells[cell_index].vertices[3];
+              raw_cells[cell_index].vertices[3] = raw_cells[cell_index].vertices[2];
+              raw_cells[cell_index].vertices[2] = tmp;
+            }
+          break;
+        }
+
+        case 3:
         {
-          const unsigned int tmp = raw_cells[cell_index].vertices[0];
-          raw_cells[cell_index].vertices[0] = raw_cells[cell_index].vertices[1];
-          raw_cells[cell_index].vertices[1] = raw_cells[cell_index].vertices[3];
-          raw_cells[cell_index].vertices[3] = raw_cells[cell_index].vertices[2];
-          raw_cells[cell_index].vertices[2] = tmp;
+          // in 3d, the situation is a bit more complicated. from above, we
+          // now know which vertex is at the origin (because 3 edges originate
+          // from it), but that still leaves 3 possible rotations of the cube.
+          // the important realization is that we can choose any of them:
+          // in all 3 rotations, all edges originate from the one vertex,
+          // and that fixes the directions of all 12 edges in the cube because
+          // these 3 cover all 3 equivalence classes! consequently, we can
+          // select an arbitrary one among the permutations -- for
+          // example the following ones:
+          static const unsigned int cube_permutations[8][8] =
+          {
+            {0,1,2,3,4,5,6,7},
+            {1,5,3,7,0,4,2,6},
+            {2,6,0,4,3,7,1,5},
+            {3,2,1,0,7,6,5,4},
+            {4,0,6,2,5,1,7,3},
+            {5,4,7,6,1,0,3,2},
+            {6,7,4,5,2,3,0,1},
+            {7,3,5,1,6,2,4,0}
+          };
+
+          unsigned int temp_vertex_indices[GeometryInfo<dim>::vertices_per_cell];
+          for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+            temp_vertex_indices[v]
+              = raw_cells[cell_index].vertices[cube_permutations[origin_vertex_of_cell][v]];
+          for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+            raw_cells[cell_index].vertices[v] = temp_vertex_indices[v];
+
+          break;
+        }
+
+        default:
+        {
+          Assert (false, ExcNotImplemented());
+        }
         }
     }
 
@@ -809,8 +956,14 @@ namespace internal
           // see which edge sets are still not oriented
           //
           // we do not need to look at each edge because if we orient edge
-          // 0, we will end up with edge 1 also oriented. there are only
-          // dim independent sets of edges
+          // 0, we will end up with edge 1 also oriented (in 2d; in 3d, there
+          // will be 3 other edges that are also oriented). there are only
+          // dim independent sets of edges, so loop over these.
+          //
+          // we need to check whether each one of these starter edges may
+          // already be oriented because the line (sheet) that connects
+          // globally parallel edges may be self-intersecting in the
+          // current cell
           for (unsigned int l=0; l<dim; ++l)
             if (edge_list[cell_list[next_cell_with_unoriented_edge].edge_indices[ParallelEdges<dim>::starter_edges[l]]].orientation_status
                 == Edge<dim>::not_oriented)
@@ -983,16 +1136,21 @@ void
 GridReordering<2,3>::reorder_cells (std::vector<CellData<2> > &cells,
                                     const bool use_new_style_ordering)
 {
-  // if necessary, convert to old-style format
-  if (use_new_style_ordering)
-    reorder_new_to_old_style(cells);
+  // if necessary, convert to old (compatibility) to new-style format
+  if (!use_new_style_ordering)
+    reorder_old_to_new_style(cells);
 
-  GridReordering<2>::reorder_cells(cells);
+  // check if grids are already
+  // consistent. if so, do
+  // nothing. if not, then do the
+  // reordering
+  if (!internal::GridReordering2d::is_consistent (cells))
+    internal::GridReordering2d::reorient(cells);
 
 
   // and convert back if necessary
-  if (use_new_style_ordering)
-    reorder_old_to_new_style(cells);
+  if (!use_new_style_ordering)
+    reorder_new_to_old_style(cells);
 }
 
 
@@ -1063,854 +1221,6 @@ GridReordering<2,3>::invert_all_cells_of_negative_grid(const std::vector<Point<3
 
 
 
-namespace internal
-{
-  namespace GridReordering3d
-  {
-    DeclException1 (ExcGridOrientError,
-                    char *,
-                    <<  "Grid Orientation Error: " << arg1);
-
-    const EdgeOrientation unoriented_edge = {'u'};
-    const EdgeOrientation forward_edge    = {'f'};
-    const EdgeOrientation backward_edge   = {'b'};
-
-
-    inline
-    bool
-    EdgeOrientation::
-    operator == (const EdgeOrientation &edge_orientation) const
-    {
-      Assert ((orientation == 'u') || (orientation == 'f') || (orientation == 'b'),
-              ExcInternalError());
-      return orientation == edge_orientation.orientation;
-    }
-
-
-
-    inline
-    bool
-    EdgeOrientation::
-    operator != (const EdgeOrientation &edge_orientation) const
-    {
-      return ! (*this == edge_orientation);
-    }
-
-
-
-    namespace ElementInfo
-    {
-      /**
-       * The numbers of the edges
-       * coming into node i are
-       * given by
-       * edge_to_node[i][k] where
-       * k=0,1,2.
-       */
-      static const unsigned int edge_to_node[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}
-      };
-
-
-      /**
-       * 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.
-       */
-      static const EdgeOrientation edge_to_node_orient[8][3] =
-      {
-        {forward_edge,  forward_edge,  forward_edge},
-        {backward_edge, forward_edge,  forward_edge},
-        {backward_edge, backward_edge, forward_edge},
-        {forward_edge,  backward_edge, forward_edge},
-        {forward_edge,  forward_edge,  backward_edge},
-        {backward_edge, forward_edge,  backward_edge},
-        {backward_edge, backward_edge, backward_edge},
-        {forward_edge,  backward_edge, backward_edge}
-      };
-
-      /**
-       * nodesonedge[i][0] is the
-       * start node for edge i.
-       * nodesonedge[i][1] is the
-       * end node for edge i.
-       */
-      static const unsigned int nodes_on_edge[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}
-      };
-    }
-
-
-    CheapEdge::CheapEdge (const unsigned int n0,
-                          const unsigned int n1)
-      :
-      // sort the
-      // entries so
-      // that
-      // node0<node1
-      node0(std::min (n0, n1)),
-      node1(std::max (n0, n1))
-    {}
-
-
-
-    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;
-    }
-
-
-    Edge::Edge (const unsigned int n0,
-                const unsigned int n1)
-      :
-      orientation_flag (unoriented_edge),
-      group (numbers::invalid_unsigned_int)
-    {
-      nodes[0] = n0;
-      nodes[1] = n1;
-    }
-
-
-
-    Cell::Cell ()
-    {
-      for (unsigned int i=0; i<GeometryInfo<3>::lines_per_cell; ++i)
-        {
-          edges[i] = numbers::invalid_unsigned_int;
-          local_orientation_flags[i] = forward_edge;
-        }
-
-      for (unsigned int i=0; i<GeometryInfo<3>::vertices_per_cell; ++i)
-        nodes[i] = numbers::invalid_unsigned_int;
-
-      waiting_to_be_processed = false;
-    }
-
-
-
-    Mesh::Mesh (const std::vector<CellData<3> > &incubes)
-    {
-      // copy the cells into our own
-      // internal data format.
-      const unsigned int numelems = incubes.size();
-      for (unsigned int i=0; i<numelems; ++i)
-        {
-          Cell the_cell;
-          std::copy (&incubes[i].vertices[0],
-                     &incubes[i].vertices[GeometryInfo<3>::vertices_per_cell],
-                     &the_cell.nodes[0]);
-
-          cell_list.push_back(the_cell);
-        }
-
-      // then build edges and
-      // connectivity
-      build_connectivity ();
-    }
-
-
-
-    void
-    Mesh::sanity_check () const
-    {
-      for (unsigned int i=0; i<cell_list.size(); ++i)
-        for (unsigned int j=0; j<8; ++j)
-          sanity_check_node (cell_list[i], j);
-    }
-
-
-
-    void
-    Mesh::sanity_check_node (const Cell         &c,
-                             const unsigned int local_node_num) const
-    {
-#ifdef DEBUG
-      // check that every edge
-      // coming into a node has the
-      // same node value
-
-      // Get the Local Node Numbers
-      // of the incoming edges
-      const unsigned int e0 = ElementInfo::edge_to_node[local_node_num][0];
-      const unsigned int e1 = ElementInfo::edge_to_node[local_node_num][1];
-      const unsigned int e2 = ElementInfo::edge_to_node[local_node_num][2];
-
-      // Global Edge Numbers
-      const unsigned int ge0 = c.edges[e0];
-      const unsigned int ge1 = c.edges[e1];
-      const unsigned int ge2 = c.edges[e2];
-
-      const EdgeOrientation or0 = ElementInfo::edge_to_node_orient[local_node_num][0] ==
-                                  c.local_orientation_flags[e0] ?
-                                  forward_edge : backward_edge;
-      const EdgeOrientation or1 = ElementInfo::edge_to_node_orient[local_node_num][1] ==
-                                  c.local_orientation_flags[e1] ?
-                                  forward_edge : backward_edge;
-      const EdgeOrientation or2 = ElementInfo::edge_to_node_orient[local_node_num][2] ==
-                                  c.local_orientation_flags[e2] ?
-                                  forward_edge : backward_edge;
-
-      // Make sure that edges agree
-      // what the current node should
-      // be.
-      Assert ((edge_list[ge0].nodes[or0 == forward_edge ? 0 : 1] ==
-               edge_list[ge1].nodes[or1 == forward_edge ? 0 : 1])
-              &&
-              (edge_list[ge1].nodes[or1 == forward_edge ? 0 : 1] ==
-               edge_list[ge2].nodes[or2 == forward_edge ? 0 : 1]),
-              ExcMessage ("This message does not satisfy the internal "
-                          "consistency check"));
-#else
-      (void)c;
-      (void)local_node_num;
-#endif
-    }
-
-
-
-    // This is the guts of the matter...
-    void Mesh::build_connectivity ()
-    {
-      const unsigned int n_cells = cell_list.size();
-
-      unsigned int n_edges = 0;
-      // Correctly build the edge
-      // list
-      {
-        // edge_map stores the
-        // edge_number associated
-        // with a given CheapEdge
-        std::map<CheapEdge,unsigned int> edge_map;
-        unsigned int ctr = 0;
-        for (unsigned int cur_cell_id = 0;
-             cur_cell_id<n_cells;
-             ++cur_cell_id)
-          {
-            // Get the local node
-            // numbers on edge
-            // edge_num
-            const Cell &cur_cell = cell_list[cur_cell_id];
-
-            for (unsigned short int edge_num = 0;
-                 edge_num<12;
-                 ++edge_num)
-              {
-                unsigned int gl_edge_num = 0;
-                EdgeOrientation l_edge_orient = forward_edge;
-
-                // Construct the
-                // CheapEdge
-                const unsigned int
-                node0 = cur_cell.nodes[ElementInfo::nodes_on_edge[edge_num][0]],
-                node1 = cur_cell.nodes[ElementInfo::nodes_on_edge[edge_num][1]];
-                const 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;
-
-                    // put the edge
-                    // into the
-                    // global edge
-                    // list
-                    edge_list.push_back(Edge(node0,node1));
-                    ctr++;
-                  }
-                else
-                  {
-                    // get edge_num
-                    // from hash_map
-                    gl_edge_num = edge_map[cur_edge];
-                    if (edge_list[gl_edge_num].nodes[0] != node0)
-                      l_edge_orient = backward_edge;
-                  }
-                // set edge number to
-                // edgenum
-                cell_list[cur_cell_id].edges[edge_num] = gl_edge_num;
-                cell_list[cur_cell_id].local_orientation_flags[edge_num]
-                  = l_edge_orient;
-              }
-          }
-        n_edges = ctr;
-      }
-
-      // Count each of the edges.
-      {
-        std::vector<int> edge_count(n_edges,0);
-
-
-        // Count every time an edge
-        // occurs in a cube.
-        for (unsigned int cur_cell_id=0; cur_cell_id<n_cells; ++cur_cell_id)
-          for (unsigned short int edge_num = 0; edge_num<12; ++edge_num)
-            ++edge_count[cell_list[cur_cell_id].edges[edge_num]];
-
-        // So we now know how many
-        // cubes contain a given
-        // edge. Just need to store
-        // the list of cubes in the
-        // edge
-
-        // Allocate the space for the
-        // neighbor list
-        for (unsigned int cur_edge_id=0; cur_edge_id<n_edges; ++cur_edge_id)
-          edge_list[cur_edge_id].neighboring_cubes
-          .resize (edge_count[cur_edge_id]);
-
-        // Store the position of the
-        // current neighbor in the
-        // edge's neighbor list
-        std::vector<int> cur_cell_edge_list_posn(n_edges,0);
-        for (unsigned int cur_cell_id=0; cur_cell_id<n_cells; ++cur_cell_id)
-          for (unsigned short int edge_num=0; edge_num<12; ++edge_num)
-            {
-              const unsigned int
-              gl_edge_id = cell_list[cur_cell_id].edges[edge_num];
-              Edge &cur_edge = edge_list[gl_edge_id];
-              cur_edge.neighboring_cubes[cur_cell_edge_list_posn[gl_edge_id]]
-                = cur_cell_id;
-              cur_cell_edge_list_posn[gl_edge_id]++;
-            }
-      }
-    }
-
-
-
-    void
-    Mesh::export_to_deal_format (std::vector<CellData<3> > &outcubes) const
-    {
-      Assert (outcubes.size() == cell_list.size(),
-              ExcInternalError());
-
-      // simply overwrite the output
-      // array with the new
-      // information
-      for (unsigned int i=0; i<cell_list.size(); ++i)
-        std::copy (&cell_list[i].nodes[0],
-                   &cell_list[i].nodes[GeometryInfo<3>::vertices_per_cell],
-                   &outcubes[i].vertices[0]);
-    }
-
-
-
-    Orienter::Orienter (const std::vector<CellData<3> > &incubes)
-      :
-      mesh (incubes),
-      cur_posn (0),
-      marker_cube (0),
-      cur_edge_group  (0)
-    {
-      for (unsigned int i = 0; i<12; ++i)
-        edge_orient_array[i] = false;
-    }
-
-
-
-    bool Orienter::orient_mesh (std::vector<CellData<3> > &incubes)
-    {
-      Orienter orienter (incubes);
-
-      // First check that the mesh is
-      // sensible
-      orienter.mesh.sanity_check ();
-
-      // Orient the mesh
-
-      // if not successful, break here, else go
-      // on
-      if (!orienter.orient_edges ())
-        return false;
-
-      // Now we have a bunch of oriented
-      // edges int the structure we only
-      // have to turn the cubes so they
-      // match the edge orientation.
-      orienter.orient_cubes ();
-
-      // Copy the elements from our
-      // internal structure back into
-      // their original location.
-      orienter.mesh.export_to_deal_format (incubes);
-      // reordering was successful
-      return true;
-    }
-
-    /**
-     * This assigns an orientation
-     * to each edge so that every
-     * cube is a rotated Deal.II
-     * cube.
-     */
-    bool Orienter::orient_edges ()
-    {
-      // While there are still cubes
-      // to orient
-      while (get_next_unoriented_cube())
-        // And there are edges in
-        // the cube to orient
-        while (orient_next_unoriented_edge())
-          {
-            // Make all the sides
-            // in the current set
-            // match
-            orient_edges_in_current_cube();
-
-            // Add the adjacent
-            // cubes to the list
-            // for processing
-            get_adjacent_cubes();
-            // Start working on
-            // this list of cubes
-            while (get_next_active_cube())
-              {
-                // Make sure the
-                // Cube doesn't
-                // have a
-                // contradiction
-                if (!cell_is_consistent(cur_posn))
-                  return false;
-
-                // If we needed to
-                // orient any edges
-                // in the current
-                // cube then we may
-                // have to process
-                // the neighbor.
-                if (orient_edges_in_current_cube())
-                  get_adjacent_cubes();
-              }
-
-            // start the next sheet
-            // (equivalence class
-            // of edges)
-            ++cur_edge_group;
-          }
-      return true;
-    }
-
-
-
-    bool Orienter::get_next_unoriented_cube ()
-    {
-      // The last cube in the list
-      const unsigned int n_cubes = mesh.cell_list.size();
-      // Keep shifting along the list
-      // until we find a cube which
-      // is not fully oriented or the
-      // end.
-      while ( (marker_cube<n_cubes) &&
-              (is_oriented(marker_cube)) )
-        ++marker_cube;
-      cur_posn = marker_cube;
-      // Return true if we now point
-      // at a valid cube.
-      return (cur_posn < n_cubes);
-    }
-
-
-
-    bool Orienter::is_oriented (const unsigned int cell_num) const
-    {
-      for (unsigned int i=0; i<12; ++i)
-        if (mesh.edge_list[mesh.cell_list[cell_num].edges[i]].orientation_flag
-            == unoriented_edge)
-          return false;
-      return true;
-    }
-
-
-
-    bool
-    Orienter::cell_is_consistent(const unsigned int cell_num) const
-    {
-
-      const Cell &c = mesh.cell_list[cell_num];
-
-      // Checks that all oriented
-      // edges in the group are
-      // oriented consistently.
-      for (unsigned int group=0; group<3; ++group)
-        {
-          // When a nonzero
-          // orientation is first
-          // encountered in the group
-          // it is stored in this
-          EdgeOrientation value = unoriented_edge;
-          // Loop over all parallel
-          // edges
-          for (unsigned int i=4*group; i<4*(group+1); ++i)
-            {
-              // If the edge has
-              // orientation
-              if ((c.local_orientation_flags[i] !=
-                   unoriented_edge)
-                  &&
-                  (mesh.edge_list[c.edges[i]].orientation_flag !=
-                   unoriented_edge))
-                {
-                  const EdgeOrientation this_edge_direction
-                    = (c.local_orientation_flags[i]
-                       == mesh.edge_list[c.edges[i]].orientation_flag  ?
-                       forward_edge : backward_edge);
-
-                  // If we haven't
-                  // seen an oriented
-                  // edge before,
-                  // then store its
-                  // value:
-                  if (value == unoriented_edge)
-                    value = this_edge_direction;
-                  else
-                    // If we have
-                    // seen an
-                    // oriented edge
-                    // in this group
-                    // we'd better
-                    // have the same
-                    // orientation.
-                    if (value != this_edge_direction)
-                      return false;
-                }
-            }
-        }
-      return true;
-    }
-
-
-
-    bool Orienter::orient_next_unoriented_edge ()
-    {
-      cur_posn = marker_cube;
-      const Cell &c = mesh.cell_list[cur_posn];
-      unsigned int edge = 0;
-
-      // search for the unoriented
-      // side
-      while ((edge<12) &&
-             (mesh.edge_list[c.edges[edge]].orientation_flag !=
-              unoriented_edge))
-        ++edge;
-
-      // if we found none then return
-      // false
-      if (edge == 12)
-        return false;
-
-      // Which edge group we're in.
-      const unsigned int edge_group = edge/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 (unsigned int j = edge_group*4; j<edge_group*4+4; ++j)
-        Assert (mesh.edge_list[c.edges[j]].orientation_flag ==
-                unoriented_edge,
-                ExcGridOrientError("Tried to orient edge when other edges "
-                                   "in group are already oriented!"));
-
-      // Make the edge alignment
-      // match that of the local
-      // cube.
-      mesh.edge_list[c.edges[edge]].orientation_flag
-        = c.local_orientation_flags[edge];
-      mesh.edge_list[c.edges[edge]].group = cur_edge_group;
-
-      // Remember that we have oriented
-      // this edge in the current cell.
-      edge_orient_array[edge] = true;
-
-      return true;
-    }
-
-
-
-    bool Orienter::orient_edges_in_current_cube ()
-    {
-      for (unsigned int edge_group=0; edge_group<3; ++edge_group)
-        if (orient_edge_set_in_current_cube(edge_group) == true)
-          return true;
-
-      return false;
-    }
-
-
-
-    bool
-    Orienter::orient_edge_set_in_current_cube (const unsigned int n)
-    {
-      const Cell &c = mesh.cell_list[cur_posn];
-
-      // Check if any edge is
-      // oriented
-      unsigned int n_oriented = 0;
-      EdgeOrientation glorient   = unoriented_edge;
-      unsigned int edge_flags = 0;
-      unsigned int cur_flag   = 1;
-      for (unsigned int i = 4*n; i<4*(n+1); ++i, cur_flag<<=1)
-        {
-          if ((mesh.edge_list[c.edges[i]].orientation_flag !=
-               unoriented_edge)
-              &&
-              (c.local_orientation_flags[i] !=
-               unoriented_edge))
-            {
-              ++n_oriented;
-
-              const EdgeOrientation orient
-                = (mesh.edge_list[c.edges[i]].orientation_flag ==
-                   c.local_orientation_flags[i] ?
-                   forward_edge : backward_edge);
-
-              if (glorient == unoriented_edge)
-                glorient = orient;
-              else
-                AssertThrow(orient == glorient,
-                            ExcGridOrientError("Attempted to Orient Misaligned cube"));
-            }
-          else
-            edge_flags |= cur_flag;
-        }
-
-      // were any of the sides
-      // oriented?  were they all
-      // already oriented?
-      if ((glorient == unoriented_edge) || (n_oriented == 4))
-        return false;
-
-      // If so orient all edges
-      // consistently.
-      cur_flag = 1;
-      for (unsigned int i=4*n; i<4*(n+1); ++i, cur_flag<<=1)
-        if ((edge_flags & cur_flag) != 0)
-          {
-            mesh.edge_list[c.edges[i]].orientation_flag
-              = (c.local_orientation_flags[i] == glorient ?
-                 forward_edge : backward_edge);
-
-            mesh.edge_list[c.edges[i]].group = cur_edge_group;
-            // Remember that we have oriented
-            // this edge in the current cell.
-            edge_orient_array[i] = true;
-          }
-
-      return true;
-    }
-
-
-
-    void Orienter::get_adjacent_cubes ()
-    {
-      const Cell &c = mesh.cell_list[cur_posn];
-      for (unsigned int e=0; e<12; ++e)
-        // Only need to add the adjacent
-        // cubes for edges we recently
-        // oriented
-        if (edge_orient_array[e] == true)
-          {
-            const Edge &the_edge = mesh.edge_list[c.edges[e]];
-            for (unsigned int local_cube_num = 0;
-                 local_cube_num < the_edge.neighboring_cubes.size();
-                 ++local_cube_num)
-              {
-                const unsigned int
-                global_cell_num = the_edge.neighboring_cubes[local_cube_num];
-                Cell &ncell = mesh.cell_list[global_cell_num];
-
-                // If the cell is waiting to be
-                // processed we dont want to add
-                // it to the list a second time.
-                if (!ncell.waiting_to_be_processed)
-                  {
-                    sheet_to_process.push_back(global_cell_num);
-                    ncell.waiting_to_be_processed = true;
-                  }
-              }
-          }
-      // we're done with this cube so
-      // clear its processing flags.
-      for (unsigned int e=0; e<12; ++e)
-        edge_orient_array[e] = false;
-
-    }
-
-
-
-    bool Orienter::get_next_active_cube ()
-    {
-      // Mark the curent Cube as
-      // finished with.
-      Cell &c = mesh.cell_list[cur_posn];
-      c.waiting_to_be_processed = false;
-      if (sheet_to_process.empty() == false)
-        {
-          cur_posn = sheet_to_process.back();
-          sheet_to_process.pop_back();
-          return true;
-        }
-      return false;
-    }
-
-
-    void Orienter::orient_cubes ()
-    {
-      // We assume that the mesh has
-      // all edges oriented already.
-
-      // 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)
-      static 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<mesh.cell_list.size(); ++i)
-        {
-          Cell &the_cell = 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)
-          EdgeOrientation local_edge_orientation[12];
-          for (unsigned int j = 0; j<12; ++j)
-            {
-              // get the global edge
-              const Edge &the_edge = mesh.edge_list[the_cell.edges[j]];
-              // All edges should be
-              // oriented at this
-              // stage..
-              Assert (the_edge.orientation_flag != unoriented_edge,
-                      ExcGridOrientError ("Unoriented edge encountered"));
-              // calculate whether it
-              // points the right way
-              // or not
-              local_edge_orientation[j] = (the_cell.local_orientation_flags[j] ==
-                                           the_edge.orientation_flag ?
-                                           forward_edge : backward_edge);
-            }
-
-          // Here the number of
-          // incoming edges is
-          // tallied for each node.
-          unsigned int perm_num = numbers::invalid_unsigned_int;
-          for (unsigned int node_num=0; node_num<8; ++node_num)
-            {
-              // The local edge
-              // numbers coming into
-              // the node
-              const unsigned int e0 = ElementInfo::edge_to_node[node_num][0];
-              const unsigned int e1 = ElementInfo::edge_to_node[node_num][1];
-              const unsigned int e2 = ElementInfo::edge_to_node[node_num][2];
-
-              // The local
-              // orientation of the
-              // edge coming into the
-              // node.
-              const EdgeOrientation sign0 = ElementInfo::edge_to_node_orient[node_num][0];
-              const EdgeOrientation sign1 = ElementInfo::edge_to_node_orient[node_num][1];
-              const EdgeOrientation sign2 = ElementInfo::edge_to_node_orient[node_num][2];
-
-              // Add one to the total
-              // for each edge
-              // pointing in
-              Assert (local_edge_orientation[e0] != unoriented_edge,
-                      ExcInternalError());
-              Assert (local_edge_orientation[e1] != unoriented_edge,
-                      ExcInternalError());
-              Assert (local_edge_orientation[e2] != unoriented_edge,
-                      ExcInternalError());
-
-              const unsigned int
-              total  = (((local_edge_orientation[e0] == sign0) ? 1 : 0)
-                        +((local_edge_orientation[e1] == sign1) ? 1 : 0)
-                        +((local_edge_orientation[e2] == sign2) ? 1 : 0));
-
-              if (total == 3)
-                {
-                  Assert (perm_num == numbers::invalid_unsigned_int,
-                          ExcGridOrientError("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 != numbers::invalid_unsigned_int,
-                  ExcGridOrientError("No node having 3 incoming edges found in curent hex."));
-
-          // So use the appropriate
-          // rotation to get the new
-          // cube
-          unsigned int temp[8];
-          for (unsigned int v=0; v<8; ++v)
-            temp[v] = the_cell.nodes[CubePermutations[perm_num][v]];
-          for (unsigned int v=0; v<8; ++v)
-            the_cell.nodes[v] = temp[v];
-        }
-    }
-  } // namespace GridReordering3d
-} // namespace internal
-
-
-
 template<>
 void
 GridReordering<3>::reorder_cells (std::vector<CellData<3> > &cells,
@@ -1919,28 +1229,25 @@ GridReordering<3>::reorder_cells (std::vector<CellData<3> > &cells,
   Assert (cells.size() != 0,
           ExcMessage("List of elements to orient must have at least one cell"));
 
-  // if necessary, convert to old-style format
-  if (use_new_style_ordering)
-    reorder_new_to_old_style(cells);
-
-  // create a backup to use if GridReordering
-  // was not successful
-  std::vector<CellData<3> > backup=cells;
-
-  // This does the real work
-  const bool success=
-    internal::GridReordering3d::Orienter::orient_mesh (cells);
+  // if necessary, convert to new-style format
+  if (use_new_style_ordering == false)
+    reorder_old_to_new_style(cells);
 
-  // if reordering was not successful use
-  // original connectivity, otherwise do
-  // nothing (i.e. use the reordered
-  // connectivity)
-  if (!success)
-    cells=backup;
+  // check if grids are already consistent. if so, do
+  // nothing. if not, then do the reordering
+  if (!internal::GridReordering2d::is_consistent (cells))
+    try
+      {
+        internal::GridReordering2d::reorient(cells);
+      }
+    catch (const internal::GridReordering2d::ExcMeshNotOrientable &)
+      {
+        // the mesh is not orientable
+      }
 
   // and convert back if necessary
-  if (use_new_style_ordering)
-    reorder_old_to_new_style(cells);
+  if (use_new_style_ordering == false)
+    reorder_new_to_old_style(cells);
 }
 
 

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.