]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move GridReordering::reorder_cells(). 12082/head
authorDavid Wells <drwells@email.unc.edu>
Thu, 22 Apr 2021 15:48:15 +0000 (11:48 -0400)
committerDavid Wells <drwells@email.unc.edu>
Thu, 22 Apr 2021 20:13:09 +0000 (16:13 -0400)
doc/news/changes/incompatibilities/20210422DavidWells [new file with mode: 0644]
include/deal.II/grid/grid_reordering.h
include/deal.II/grid/grid_tools.h
source/grid/grid_generator.cc
source/grid/grid_in.cc
source/grid/grid_reordering.cc
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in

diff --git a/doc/news/changes/incompatibilities/20210422DavidWells b/doc/news/changes/incompatibilities/20210422DavidWells
new file mode 100644 (file)
index 0000000..d371ebd
--- /dev/null
@@ -0,0 +1,6 @@
+Deprecated: The GridReordering class as well as
+Triangulation::create_triangulation_compatibility have been deprecated.
+These functions use the old-style (before 5.2) numbering and have been
+unofficially deprecated since 2005.
+<br>
+(David Wells, 2021/04/22)
index 6d7d76ac78c6e8b1165aa562995a4ebf1feb8d07..430784db7aaf115df256b5b9202adc73466188d7 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
-/**
- * An exception that is thrown whenever the edges of a mesh are not
- * orientable.
- */
-DeclExceptionMsg(ExcMeshNotOrientable,
-                 "The edges of the mesh are not consistently orientable.");
-
-
 /**
  * A class implementing various grid reordering algorithms. For more information
  * see the @ref reordering "reordering module".
+ *
+ * @deprecated Use GridTools::invert_all_negative_measure_cells() or
+ * GridTools::consistently_order_cells() instead of the functions provided by
+ * this class. Usage of the old-style numbering is deprecated.
  */
 template <int dim, int spacedim = dim>
-class GridReordering
+class DEAL_II_DEPRECATED_EARLY GridReordering
 {
 public:
   /**
@@ -54,7 +50,10 @@ public:
    * vertices within a cell. If false (the default), then use the "old-style"
    * ordering of vertices within cells used by deal.II before version 5.2 and
    * as explained in the documentation of this class.
+   *
+   * @deprecated Use GridTools::consistently_order_cells() instead.
    */
+  DEAL_II_DEPRECATED_EARLY
   static void
   reorder_cells(std::vector<CellData<dim>> &original_cells,
                 const bool                  use_new_style_ordering = false);
index 723ed0f80235ce69b2e05b31440feb5ea0d19c44..f86d9d9882a34599621505c1e6187e6e93e9c4a0 100644 (file)
@@ -460,6 +460,19 @@ namespace GridTools
     const std::vector<Point<spacedim>> &all_vertices,
     std::vector<CellData<dim>> &        cells);
 
+  /**
+   * Given a vector of CellData objects describing a mesh, reorder their
+   * vertices so that all lines are consistently oriented.
+   *
+   * The expectations on orientation and a discussion of this function are
+   * available in the @ref reordering "reordering module".
+   *
+   * @param cells The array of CellData objects that describe the mesh's topology.
+   */
+  template <int dim>
+  void
+  consistently_order_cells(std::vector<CellData<dim>> &cells);
+
   /*@}*/
   /**
    * @name Rotating, stretching and otherwise transforming meshes
@@ -3311,12 +3324,24 @@ namespace GridTools
                  << "The given vertex with index " << arg1
                  << " is not used in the given triangulation.");
 
-
   /*@}*/
 
 } /*namespace GridTools*/
 
 
+/**
+ * An exception that is thrown whenever the edges of a mesh are not
+ * orientable.
+ *
+ * @note for backwards compatibility with the old GridReordering class this
+ * exception is not in the GridTools namespace.
+ *
+ * @ingroup Exceptions
+ */
+DeclExceptionMsg(ExcMeshNotOrientable,
+                 "The edges of the mesh are not consistently orientable.");
+
+
 
 /* ----------------- Template function --------------- */
 
index f7e99e478c6c17c8bb0868069e845f0e2c5d1741..1427244d7a7b892741fe571380f7b518efba3a49 100644 (file)
@@ -1997,10 +1997,7 @@ namespace GridGenerator
     cells[15].vertices[3] = 2;
     cells[15].material_id = 0;
 
-    // Must call this to be able to create a
-    // correct triangulation in dealii, read
-    // GridReordering<> doc
-    GridReordering<dim, spacedim>::reorder_cells(cells, true);
+    GridTools::consistently_order_cells(cells);
     tria.create_triangulation(vertices, cells, SubCellData());
 
     tria.set_all_manifold_ids(0);
@@ -2261,11 +2258,6 @@ namespace GridGenerator
   // Parallelepiped implementation in 1d, 2d, and 3d. @note The
   // implementation in 1d is similar to hyper_rectangle(), and in 2d is
   // similar to parallelogram().
-  //
-  // The GridReordering::reorder_grid is made use of towards the end of
-  // this function. Thus the triangulation is explicitly constructed for
-  // all dim here since it is slightly different in that respect
-  // (cf. hyper_rectangle(), parallelogram()).
   template <int dim, int spacedim>
   void
   subdivided_parallelepiped(Triangulation<dim, spacedim> &              tria,
@@ -2484,7 +2476,7 @@ namespace GridGenerator
     // Create triangulation
     // reorder the cells to ensure that they satisfy the convention for
     // edge and face directions
-    GridReordering<dim>::reorder_cells(cells, true);
+    GridTools::consistently_order_cells(cells);
     tria.create_triangulation(points, cells, SubCellData());
 
     // Finally assign boundary indicators according to hyper_rectangle
@@ -6438,7 +6430,7 @@ namespace GridGenerator
 
     // reorder the cells to ensure that they satisfy the convention for
     // edge and face directions
-    GridReordering<dim, spacedim>::reorder_cells(cells, true);
+    GridTools::consistently_order_cells(cells);
     result.clear();
     result.create_triangulation(vertices, cells, subcell_data);
   }
@@ -6671,7 +6663,7 @@ namespace GridGenerator
           1e-6 * input.begin_active()->diameter());
         // delete_duplicated_vertices also deletes any unused vertices
         // deal with any reordering issues created by delete_duplicated_vertices
-        GridReordering<dim>::reorder_cells(output_cell_data, true);
+        GridTools::consistently_order_cells(output_cell_data);
         // clean up the boundary ids of the boundary objects: note that we
         // have to do this after delete_duplicated_vertices so that boundary
         // objects are actually duplicated at this point
@@ -7069,12 +7061,12 @@ namespace GridGenerator
 
     // use all of this to finally create the extruded 3d
     // triangulation.  it is not necessary to call
-    // GridReordering<3,3>::reorder_cells because the cells we have
+    // GridTools::consistently_order_cells() because the cells we have
     // constructed above are automatically correctly oriented. this is
     // because the 2d base mesh is always correctly oriented, and
     // extruding it automatically yields a correctly oriented 3d mesh,
     // as discussed in the edge orientation paper mentioned in the
-    // introduction to the GridReordering class.
+    // introduction to the @ref reordering "reordering module".
     result.create_triangulation(points, cells, subcell_data);
 
     for (auto manifold_id_it = priorities.rbegin();
index 0ee5e9b82e7c5c55d8eec420f7b98995d9704fdd..5fbf7c3aeccea326c46fcb03fe6ac5eaec4c6d75 100644 (file)
@@ -561,7 +561,7 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
 
       // TODO: the functions below (GridTools::delete_unused_vertices(),
       // GridTools::invert_all_negative_measure_cells(),
-      // GridReordering::reorder_cells()) need to be
+      // GridTools::consistently_order_cells()) need to be
       // revisited for simplex/mixed meshes
 
       if (dim == 1 || (is_quad_or_hex_mesh && !is_tria_or_tet_mesh))
@@ -571,7 +571,7 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
           if (dim == spacedim)
             GridTools::invert_all_negative_measure_cells(vertices, cells);
 
-          GridReordering<dim, spacedim>::reorder_cells(cells, true);
+          GridTools::consistently_order_cells(cells);
           tria->create_triangulation(vertices, cells, subcelldata);
         }
       else
@@ -862,7 +862,7 @@ GridIn<dim, spacedim>::read_unv(std::istream &in)
   if (dim == spacedim)
     GridTools::invert_all_negative_measure_cells(vertices, cells);
 
-  GridReordering<dim, spacedim>::reorder_cells(cells, true);
+  GridTools::consistently_order_cells(cells);
 
   tria->create_triangulation(vertices, cells, subcelldata);
 }
@@ -1095,7 +1095,7 @@ GridIn<dim, spacedim>::read_ucd(std::istream &in,
   // ... and cells
   if (dim == spacedim)
     GridTools::invert_all_negative_measure_cells(vertices, cells);
-  GridReordering<dim, spacedim>::reorder_cells(cells, true);
+  GridTools::consistently_order_cells(cells);
   tria->create_triangulation(vertices, cells, subcelldata);
 }
 
@@ -1348,7 +1348,7 @@ GridIn<dim, spacedim>::read_dbmesh(std::istream &in)
   GridTools::delete_unused_vertices(vertices, cells, subcelldata);
   // ...and cells
   GridTools::invert_all_negative_measure_cells(vertices, cells);
-  GridReordering<dim, spacedim>::reorder_cells(cells, true);
+  GridTools::consistently_order_cells(cells);
   tria->create_triangulation(vertices, cells, subcelldata);
 }
 
@@ -1416,7 +1416,7 @@ GridIn<dim, spacedim>::read_xda(std::istream &in)
   GridTools::delete_unused_vertices(vertices, cells, subcelldata);
   // ... and cells
   GridTools::invert_all_negative_measure_cells(vertices, cells);
-  GridReordering<dim>::reorder_cells(cells, true);
+  GridTools::consistently_order_cells(cells);
   tria->create_triangulation(vertices, cells, subcelldata);
 }
 
@@ -2096,7 +2096,7 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
 
   // TODO: the functions below (GridTools::delete_unused_vertices(),
   // GridTools::invert_all_negative_measure_cells(),
-  // GridReordering::reorder_cells()) need to be revisited
+  // GridTools::consistently_order_cells()) need to be revisited
   // for simplex/mixed meshes
 
   if (dim == 1 || (is_quad_or_hex_mesh && !is_tria_or_tet_mesh))
@@ -2106,7 +2106,7 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
       // ... and cells
       if (dim == spacedim)
         GridTools::invert_all_negative_measure_cells(vertices, cells);
-      GridReordering<dim, spacedim>::reorder_cells(cells, true);
+      GridTools::consistently_order_cells(cells);
     }
   tria->create_triangulation(vertices, cells, subcelldata);
 
@@ -2799,7 +2799,7 @@ GridIn<2>::read_tecplot(std::istream &in)
 
   // do some cleanup on cells
   GridTools::invert_all_negative_measure_cells(vertices, cells);
-  GridReordering<dim, spacedim>::reorder_cells(cells, true);
+  GridTools::consistently_order_cells(cells);
   tria->create_triangulation(vertices, cells, subcelldata);
 }
 
@@ -2954,7 +2954,7 @@ GridIn<dim, spacedim>::read_assimp(const std::string &filename,
   GridTools::delete_unused_vertices(vertices, cells, subcelldata);
   if (dim == spacedim)
     GridTools::invert_all_negative_measure_cells(vertices, cells);
-  GridReordering<dim, spacedim>::reorder_cells(cells, true);
+  GridTools::consistently_order_cells(cells);
   tria->create_triangulation(vertices, cells, subcelldata);
 
 #else
index 76f532e1d53ea625f879559135aaee7c5f53d17b..5dcc6c0157c18f2d81aaabf89d4849f83feeae83 100644 (file)
 DEAL_II_NAMESPACE_OPEN
 
 
-namespace
-{
-  /**
-   * A simple data structure denoting an edge, i.e., the ordered pair
-   * of its vertex indices. This is only used in the is_consistent()
-   * function.
-   */
-  struct CheapEdge
-  {
-    /**
-     * Construct an edge from the global indices of its two vertices.
-     */
-    CheapEdge(const unsigned int v0, const unsigned int v1)
-      : v0(v0)
-      , v1(v1)
-    {}
-
-    /**
-     * Comparison operator for edges. It compares based on the
-     * lexicographic ordering of the two vertex indices.
-     */
-    bool
-    operator<(const CheapEdge &e) const
-    {
-      return ((v0 < e.v0) || ((v0 == e.v0) && (v1 < e.v1)));
-    }
-
-  private:
-    /**
-     * The global indices of the vertices that define the edge.
-     */
-    const unsigned int v0, v1;
-  };
-
-
-  /**
-   * A function that determines whether the edges in a mesh are
-   * already consistently oriented. It does so by adding all edges
-   * of all cells into a set (which automatically eliminates
-   * duplicates) but before that checks whether the reverse edge is
-   * already in the set -- which would imply that a neighboring cell
-   * is inconsistently oriented.
-   */
-  template <int dim>
-  bool
-  is_consistent(const std::vector<CellData<dim>> &cells)
-  {
-    std::set<CheapEdge> edges;
-
-    for (typename std::vector<CellData<dim>>::const_iterator c = cells.begin();
-         c != cells.end();
-         ++c)
-      {
-        // construct the edges in reverse order. for each of them,
-        // ensure that the reverse edge is not yet in the list of
-        // edges (return false if the reverse edge already *is* in
-        // the list) and then add the actual edge to it; std::set
-        // eliminates duplicates automatically
-        for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
-          {
-            const CheapEdge reverse_edge(
-              c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 1)],
-              c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 0)]);
-            if (edges.find(reverse_edge) != edges.end())
-              return false;
-
-
-            // ok, not. insert edge in correct order
-            const CheapEdge correct_edge(
-              c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 0)],
-              c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 1)]);
-            edges.insert(correct_edge);
-          }
-      }
-
-    // no conflicts found, so return true
-    return true;
-  }
-
-
-  /**
-   * A structure that describes some properties of parallel edges
-   * such as what starter edges are (i.e., representative elements
-   * of the sets of parallel edges within a cell) and what the set
-   * of parallel edges to each edge is.
-   */
-  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];
-
-    /**
-     * 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 <>
-  const unsigned int ParallelEdges<2>::starter_edges[2] = {0, 2};
-
-  template <>
-  const unsigned int ParallelEdges<2>::parallel_edges[4][1] = {{1},
-                                                               {0},
-                                                               {3},
-                                                               {2}};
-
-  template <>
-  const unsigned int ParallelEdges<3>::starter_edges[3] = {0, 2, 8};
-
-  template <>
-  const unsigned int ParallelEdges<3>::parallel_edges[12][3] = {
-    {1, 4, 5},   // line 0
-    {0, 4, 5},   // line 1
-    {3, 6, 7},   // line 2
-    {2, 6, 7},   // line 3
-    {0, 1, 5},   // line 4
-    {0, 1, 4},   // line 5
-    {2, 3, 7},   // line 6
-    {2, 3, 6},   // line 7
-    {9, 10, 11}, // line 8
-    {8, 10, 11}, // line 9
-    {8, 9, 11},  // line 10
-    {8, 9, 10}   // line 11
-  };
-
-
-  /**
-   * A structure that store the index of a cell and, crucially, how a
-   * given edge relates to this cell.
-   */
-  struct AdjacentCell
-  {
-    /**
-     * Default constructor. Initialize the fields with invalid values.
-     */
-    AdjacentCell()
-      : cell_index(numbers::invalid_unsigned_int)
-      , edge_within_cell(numbers::invalid_unsigned_int)
-    {}
-
-    /**
-     * Constructor. Initialize the fields with the given values.
-     */
-    AdjacentCell(const unsigned int cell_index,
-                 const unsigned int edge_within_cell)
-      : cell_index(cell_index)
-      , edge_within_cell(edge_within_cell)
-    {}
-
-
-    unsigned int cell_index;
-    unsigned int edge_within_cell;
-  };
-
-
-
-  template <int dim>
-  class AdjacentCells;
-
-  /**
-   * A class that represents all of the cells adjacent to a given edge.
-   * This class corresponds to the 2d case where each edge has at most
-   * two adjacent cells.
-   */
-  template <>
-  class AdjacentCells<2>
-  {
-  public:
-    /**
-     * An iterator that allows iterating over all cells adjacent
-     * to the edge represented by the current object.
-     */
-    using const_iterator = const AdjacentCell *;
-
-    /**
-     * Add the given cell to the collection of cells adjacent to
-     * the edge this object corresponds to. Since we are covering
-     * the 2d case, the set of adjacent cells currently
-     * represented by this object must have either zero or
-     * one element already, since we can not add more than two
-     * adjacent cells for each edge.
-     */
-    void
-    push_back(const AdjacentCell &adjacent_cell)
-    {
-      if (adjacent_cells[0].cell_index == numbers::invalid_unsigned_int)
-        adjacent_cells[0] = adjacent_cell;
-      else
-        {
-          Assert(adjacent_cells[1].cell_index == numbers::invalid_unsigned_int,
-                 ExcInternalError());
-          adjacent_cells[1] = adjacent_cell;
-        }
-    }
-
-
-    /**
-     * Return an iterator to the first valid cell stored as adjacent to the
-     * edge represented by the current object.
-     */
-    const_iterator
-    begin() const
-    {
-      return adjacent_cells;
-    }
-
-
-    /**
-     * Return an iterator to the element past the last valid cell stored
-     * as adjacent to the edge represented by the current object.
-     * @return
-     */
-    const_iterator
-    end() const
-    {
-      // check whether the current object stores zero, one, or two
-      // adjacent cells, and use this to point to the element past the
-      // last valid one
-      if (adjacent_cells[0].cell_index == numbers::invalid_unsigned_int)
-        return adjacent_cells;
-      else if (adjacent_cells[1].cell_index == numbers::invalid_unsigned_int)
-        return adjacent_cells + 1;
-      else
-        return adjacent_cells + 2;
-    }
-
-  private:
-    /**
-     * References to the (at most) two cells that are adjacent to
-     * the edge this object corresponds to. Unused elements are
-     * default-initialized and have invalid values; in particular,
-     * their cell_index field equals numbers::invalid_unsigned_int.
-     */
-    AdjacentCell adjacent_cells[2];
-  };
-
-
-
-  /**
-   * 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
-   * indices of the two vertices, and the indices of the adjacent
-   * cells (together with a description *where* in each of the
-   * adjacent cells the edge is located). It also includes the
-   * (global) direction of the edge: either from the first vertex to
-   * the second, the other way around, or so far undetermined.
-   */
-  template <int dim>
-  class Edge
-  {
-  public:
-    /**
-     * Constructor. Create the edge based on the information given
-     * in @p cell, and selecting the edge with number @p edge_number
-     * within this cell. Initialize the edge as unoriented.
-     */
-    Edge(const CellData<dim> &cell, const unsigned int edge_number)
-      : orientation_status(not_oriented)
-    {
-      Assert(edge_number < GeometryInfo<dim>::lines_per_cell,
-             ExcInternalError());
-
-      // copy vertices for this particular line
-      vertex_indices[0] =
-        cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(edge_number, 0)];
-      vertex_indices[1] =
-        cell.vertices[GeometryInfo<dim>::line_to_cell_vertices(edge_number, 1)];
-
-      // bring them into standard orientation
-      if (vertex_indices[0] > vertex_indices[1])
-        std::swap(vertex_indices[0], vertex_indices[1]);
-    }
-
-    /**
-     * Comparison operator for edges. It compares based on the
-     * lexicographic ordering of the two vertex indices.
-     */
-    bool
-    operator<(const Edge<dim> &e) const
-    {
-      return ((vertex_indices[0] < e.vertex_indices[0]) ||
-              ((vertex_indices[0] == e.vertex_indices[0]) &&
-               (vertex_indices[1] < e.vertex_indices[1])));
-    }
-
-    /**
-     * Compare two edges for equality based on their vertex indices.
-     */
-    bool
-    operator==(const Edge<dim> &e) const
-    {
-      return ((vertex_indices[0] == e.vertex_indices[0]) &&
-              (vertex_indices[1] == e.vertex_indices[1]));
-    }
-
-    /**
-     * The global indices of the two vertices that bound this edge. These
-     * will be ordered so that the first index is less than the second.
-     */
-    unsigned int vertex_indices[2];
-
-    /**
-     * An enum that indicates the direction of this edge with
-     * regard to the two vertices that bound it.
-     */
-    enum OrientationStatus
-    {
-      not_oriented,
-      forward,
-      backward
-    };
-
-    OrientationStatus orientation_status;
-
-    /**
-     * Store the set of cells adjacent to this edge (these cells then
-     * also store *where* in the cell the edge is located).
-     */
-    AdjacentCells<dim> adjacent_cells;
-  };
-
-
-
-  /**
-   * A data structure that represents a cell with all of its vertices
-   * and edges.
-   */
-  template <int dim>
-  struct Cell
-  {
-    /**
-     * Construct a Cell object from a CellData object. Also take a
-     * (sorted) list of edges and to point the edges of the current
-     * object into this list of edges.
-     */
-    Cell(const CellData<dim> &c, const std::vector<Edge<dim>> &edge_list)
-    {
-      for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
-        vertex_indices[i] = c.vertices[i];
-
-      // now for each of the edges of this cell, find the location inside the
-      // given edge_list array and store than index
-      for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
-        {
-          const Edge<dim> e(c, l);
-          edge_indices[l] =
-            (std::lower_bound(edge_list.begin(), edge_list.end(), e) -
-             edge_list.begin());
-          Assert(edge_indices[l] < edge_list.size(), ExcInternalError());
-          Assert(edge_list[edge_indices[l]] == e, ExcInternalError())
-        }
-    }
-
-    /**
-     * A list of global indices for the vertices that bound this cell.
-     */
-    unsigned int vertex_indices[GeometryInfo<dim>::vertices_per_cell];
-
-    /**
-     * A list of indices into the 'edge_list' array passed to the constructor
-     * for the edges of the current cell.
-     */
-    unsigned int edge_indices[GeometryInfo<dim>::lines_per_cell];
-  };
-
-
-
-  template <int dim>
-  class EdgeDeltaSet;
-
-  /**
-   * A class that represents by how much the set of parallel edges
-   * grows in each step. In the graph orientation paper, this set is
-   * called $\Delta_k$, thus the name.
-   *
-   * In 2d, this set can only include zero, one, or two elements.
-   * Consequently, the appropriate data structure is one in which
-   * we store at most 2 elements in a fixed sized data structure.
-   */
-  template <>
-  class EdgeDeltaSet<2>
-  {
-  public:
-    /**
-     * Iterator type for the elements of the set.
-     */
-    using const_iterator = const unsigned int *;
-
-    /**
-     * Default constructor. Initialize both slots as unused, corresponding
-     * to an empty set.
-     */
-    EdgeDeltaSet()
-    {
-      edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int;
-    }
-
-
-    /**
-     * Delete the elements of the set by marking both slots as unused.
-     */
-    void
-    clear()
-    {
-      edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int;
-    }
-
-    /**
-     * Insert one element into the set. This will fail if the set already
-     * has two elements.
-     */
-    void
-    insert(const unsigned int edge_index)
-    {
-      if (edge_indices[0] == numbers::invalid_unsigned_int)
-        edge_indices[0] = edge_index;
-      else
-        {
-          Assert(edge_indices[1] == numbers::invalid_unsigned_int,
-                 ExcInternalError());
-          edge_indices[1] = edge_index;
-        }
-    }
-
-
-    /**
-     * Return an iterator pointing to the first element of the set.
-     */
-    const_iterator
-    begin() const
-    {
-      return edge_indices;
-    }
-
-
-    /**
-     * Return an iterator pointing to the element past the last used one.
-     */
-    const_iterator
-    end() const
-    {
-      // check whether the current object stores zero, one, or two
-      // indices, and use this to point to the element past the
-      // last valid one
-      if (edge_indices[0] == numbers::invalid_unsigned_int)
-        return edge_indices;
-      else if (edge_indices[1] == numbers::invalid_unsigned_int)
-        return edge_indices + 1;
-      else
-        return edge_indices + 2;
-    }
-
-  private:
-    /**
-     * Storage space to store the indices of at most two edges.
-     */
-    unsigned int edge_indices[2];
-  };
-
-
-
-  /**
-   * A class that represents by how much the set of parallel edges
-   * grows in each step. In the graph orientation paper, this set is
-   * called $\Delta_k$, thus the name.
-   *
-   * In 3d, this set can have arbitrarily many elements, unlike the
-   * 2d case specialized above. Consequently, we simply represent
-   * the data structure with a std::set. Class derivation ensures
-   * that we simply inherit all of the member functions of the
-   * base class.
-   */
-  template <>
-  class EdgeDeltaSet<3> : public std::set<unsigned int>
-  {};
-
-
-
-  /**
-   * From a list of cells, build a sorted vector that contains all of the edges
-   * that exist in the mesh.
-   */
-  template <int dim>
-  std::vector<Edge<dim>>
-  build_edges(const std::vector<CellData<dim>> &cells)
-  {
-    // build the edge list for all cells. because each cell has
-    // GeometryInfo<dim>::lines_per_cell edges, the total number
-    // of edges is this many times the number of cells. of course
-    // some of them will be duplicates, and we throw them out below
-    std::vector<Edge<dim>> edge_list;
-    edge_list.reserve(cells.size() * GeometryInfo<dim>::lines_per_cell);
-    for (unsigned int i = 0; i < cells.size(); ++i)
-      for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
-        edge_list.emplace_back(cells[i], l);
-
-    // next sort the edge list and then remove duplicates
-    std::sort(edge_list.begin(), edge_list.end());
-    edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
-                    edge_list.end());
-
-    return edge_list;
-  }
-
-
-
-  /**
-   * Build the cell list. Update the edge array to let edges know
-   * which cells are adjacent to them.
-   */
-  template <int dim>
-  std::vector<Cell<dim>>
-  build_cells_and_connect_edges(const std::vector<CellData<dim>> &cells,
-                                std::vector<Edge<dim>> &          edges)
-  {
-    std::vector<Cell<dim>> cell_list;
-    cell_list.reserve(cells.size());
-    for (unsigned int i = 0; i < cells.size(); ++i)
-      {
-        // create our own data structure for the cells and let it
-        // connect to the edges array
-        cell_list.emplace_back(cells[i], edges);
-
-        // then also inform the edges that they are adjacent
-        // to the current cell, and where within this cell
-        for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
-          edges[cell_list.back().edge_indices[l]].adjacent_cells.push_back(
-            AdjacentCell(i, l));
-      }
-    Assert(cell_list.size() == cells.size(), ExcInternalError());
-
-    return cell_list;
-  }
-
-
-
-  /**
-   * Return the index within 'cells' of the first cell that has at least one
-   * edge that is not yet oriented.
-   */
-  template <int dim>
-  unsigned int
-  get_next_unoriented_cell(const std::vector<Cell<dim>> &cells,
-                           const std::vector<Edge<dim>> &edges,
-                           const unsigned int            current_cell)
-  {
-    for (unsigned int c = current_cell; c < cells.size(); ++c)
-      for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
-        if (edges[cells[c].edge_indices[l]].orientation_status ==
-            Edge<dim>::not_oriented)
-          return c;
-
-    return numbers::invalid_unsigned_int;
-  }
-
-
-
-  /**
-   * Given a set of cells and edges, orient all edges that are
-   * (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<dim>> &cells,
-                                   std::vector<Edge<dim>> &      edges,
-                                   const unsigned int            cell,
-                                   const unsigned int            local_edge)
-  {
-    // 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
-    // implementation used till 2016, let us just choose the
-    // direction so that it matches what we have in the given cell.
-    //
-    // in fact, in what can only be assumed to be a bug in the
-    // original implementation, after orienting all edges, the code
-    // that rotates the cells so that they match edge orientations
-    // (see the rotate_cell() function below) rotated the cell two
-    // more times by 90 degrees. this is ok -- it simply flips all
-    // edge orientations, which leaves them valid. rather than do
-    // the same in the current implementation, we can achieve the
-    // 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 =
-        (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
-    else
-      {
-        Assert(
-          edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
-            cells[cell].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(
-              local_edge, 1)],
-          ExcInternalError());
-        Assert(
-          edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] ==
-            cells[cell].vertex_indices[GeometryInfo<dim>::line_to_cell_vertices(
-              local_edge, 0)],
-          ExcInternalError());
-
-        // 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 =
-          (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
-      }
-
-    // walk outward from the given edge as described in
-    // the algorithm in the paper that documents all of
-    // this
-    //
-    // note that in 2d, each of the Deltas can at most
-    // contain two elements, whereas in 3d it can be arbitrarily many
-    EdgeDeltaSet<dim> Delta_k;
-    EdgeDeltaSet<dim> Delta_k_minus_1;
-    Delta_k_minus_1.insert(cells[cell].edge_indices[local_edge]);
-
-    while (Delta_k_minus_1.begin() !=
-           Delta_k_minus_1.end()) // while set is not empty
-      {
-        Delta_k.clear();
-
-        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 (typename AdjacentCells<dim>::const_iterator adjacent_cell =
-                   edges[*delta].adjacent_cells.begin();
-                 adjacent_cell != edges[*delta].adjacent_cells.end();
-                 ++adjacent_cell)
-              {
-                const unsigned int K = adjacent_cell->cell_index;
-                const unsigned int delta_is_edge_in_K =
-                  adjacent_cell->edge_within_cell;
-
-                // figure out the direction of delta with respect to the cell K
-                // (in the orientation in which the user has given it to us)
-                const unsigned int first_edge_vertex =
-                  (edges[*delta].orientation_status == Edge<dim>::forward ?
-                     edges[*delta].vertex_indices[0] :
-                     edges[*delta].vertex_indices[1]);
-                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>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
-                  ExcInternalError());
-
-                // 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)
-                  {
-                    // 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());
-                      }
-                  }
-              }
-          }
-
-        // finally copy the new set to the previous one
-        // (corresponding to increasing 'k' by one in the
-        // algorithm)
-        Delta_k_minus_1 = Delta_k;
-      }
-  }
-
-
-  /**
-   * Given data structures @p cell_list and @p edge_list, where
-   * all edges are already oriented, rotate the cell with
-   * index @p cell_index in such a way that its local coordinate
-   * 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<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 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<dim>::not_oriented,
-               ExcInternalError());
-        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 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,
-                             starting_vertex_of_edge +
-                               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;
-          }
-
-        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:
-          {
-            // 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 (const unsigned int v : GeometryInfo<dim>::vertex_indices())
-              temp_vertex_indices[v] =
-                raw_cells[cell_index]
-                  .vertices[cube_permutations[origin_vertex_of_cell][v]];
-            for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
-              raw_cells[cell_index].vertices[v] = temp_vertex_indices[v];
-
-            break;
-          }
-
-        default:
-          {
-            Assert(false, ExcNotImplemented());
-          }
-      }
-  }
-
-
-  /**
-   * Given a set of cells, find globally unique edge orientations
-   * and then rotate cells so that the coordinate system of the cell
-   * coincides with the coordinate systems of the adjacent edges.
-   */
-  template <int dim>
-  void
-  reorient(std::vector<CellData<dim>> &cells)
-  {
-    // first build the arrays that connect cells to edges and the other
-    // way around
-    std::vector<Edge<dim>> edge_list = build_edges(cells);
-    std::vector<Cell<dim>> cell_list =
-      build_cells_and_connect_edges(cells, edge_list);
-
-    // then loop over all cells and start orienting parallel edge sets
-    // of cells that still have non-oriented edges
-    unsigned int next_cell_with_unoriented_edge = 0;
-    while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
-              cell_list, edge_list, next_cell_with_unoriented_edge)) !=
-           numbers::invalid_unsigned_int)
-      {
-        // 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 (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)
-            orient_one_set_of_parallel_edges(
-              cell_list,
-              edge_list,
-              next_cell_with_unoriented_edge,
-              ParallelEdges<dim>::starter_edges[l]);
-
-        // ensure that we have really oriented all edges now, not just
-        // the starter edges
-        for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
-          Assert(
-            edge_list[cell_list[next_cell_with_unoriented_edge].edge_indices[l]]
-                .orientation_status != Edge<dim>::not_oriented,
-            ExcInternalError());
-      }
-
-    // now that we have oriented all edges, we need to rotate cells
-    // so that the edges point in the right direction with the now
-    // rotated coordinate system
-    for (unsigned int c = 0; c < cells.size(); ++c)
-      rotate_cell(cell_list, edge_list, c, cells);
-  }
-
-
-  // overload of the function above for 1d -- there is nothing
-  // to orient in that case
-  void reorient(std::vector<CellData<1>> &)
-  {}
-} // namespace
-
-
 // anonymous namespace for internal helper functions
 namespace
 {
@@ -1089,21 +111,7 @@ GridReordering<dim, spacedim>::reorder_cells(std::vector<CellData<dim>> &cells,
   if (use_new_style_ordering == false)
     reorder_old_to_new_style(cells);
 
-  // check if grids are already consistent. if so, do
-  // nothing. if not, then do the reordering
-  if (!is_consistent(cells))
-    try
-      {
-        reorient(cells);
-      }
-    catch (const ExcMeshNotOrientable &)
-      {
-        // the mesh is not orientable. this is acceptable if we are in 3d,
-        // as class Triangulation knows how to handle this, but it is
-        // not in 2d; in that case, re-throw the exception
-        if (dim < 3)
-          throw;
-      }
+  GridTools::consistently_order_cells(cells);
 
   // and convert back if necessary
   if (use_new_style_ordering == false)
index 8b125a912f692eab35aaa2de99848942722a59bb..f3561ae4b74f32f777f4cdc16fc97c5b3e80183f 100644 (file)
@@ -929,6 +929,1024 @@ namespace GridTools
 
 
 
+  // Functions and classes for consistently_order_cells
+  namespace
+  {
+    /**
+     * A simple data structure denoting an edge, i.e., the ordered pair
+     * of its vertex indices. This is only used in the is_consistent()
+     * function.
+     */
+    struct CheapEdge
+    {
+      /**
+       * Construct an edge from the global indices of its two vertices.
+       */
+      CheapEdge(const unsigned int v0, const unsigned int v1)
+        : v0(v0)
+        , v1(v1)
+      {}
+
+      /**
+       * Comparison operator for edges. It compares based on the
+       * lexicographic ordering of the two vertex indices.
+       */
+      bool
+      operator<(const CheapEdge &e) const
+      {
+        return ((v0 < e.v0) || ((v0 == e.v0) && (v1 < e.v1)));
+      }
+
+    private:
+      /**
+       * The global indices of the vertices that define the edge.
+       */
+      const unsigned int v0, v1;
+    };
+
+
+    /**
+     * A function that determines whether the edges in a mesh are
+     * already consistently oriented. It does so by adding all edges
+     * of all cells into a set (which automatically eliminates
+     * duplicates) but before that checks whether the reverse edge is
+     * already in the set -- which would imply that a neighboring cell
+     * is inconsistently oriented.
+     */
+    template <int dim>
+    bool
+    is_consistent(const std::vector<CellData<dim>> &cells)
+    {
+      std::set<CheapEdge> edges;
+
+      for (typename std::vector<CellData<dim>>::const_iterator c =
+             cells.begin();
+           c != cells.end();
+           ++c)
+        {
+          // construct the edges in reverse order. for each of them,
+          // ensure that the reverse edge is not yet in the list of
+          // edges (return false if the reverse edge already *is* in
+          // the list) and then add the actual edge to it; std::set
+          // eliminates duplicates automatically
+          for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+            {
+              const CheapEdge reverse_edge(
+                c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 1)],
+                c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 0)]);
+              if (edges.find(reverse_edge) != edges.end())
+                return false;
+
+
+              // ok, not. insert edge in correct order
+              const CheapEdge correct_edge(
+                c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 0)],
+                c->vertices[GeometryInfo<dim>::line_to_cell_vertices(l, 1)]);
+              edges.insert(correct_edge);
+            }
+        }
+
+      // no conflicts found, so return true
+      return true;
+    }
+
+
+    /**
+     * A structure that describes some properties of parallel edges
+     * such as what starter edges are (i.e., representative elements
+     * of the sets of parallel edges within a cell) and what the set
+     * of parallel edges to each edge is.
+     */
+    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];
+
+      /**
+       * 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 <>
+    const unsigned int ParallelEdges<2>::starter_edges[2] = {0, 2};
+
+    template <>
+    const unsigned int ParallelEdges<2>::parallel_edges[4][1] = {{1},
+                                                                 {0},
+                                                                 {3},
+                                                                 {2}};
+
+    template <>
+    const unsigned int ParallelEdges<3>::starter_edges[3] = {0, 2, 8};
+
+    template <>
+    const unsigned int ParallelEdges<3>::parallel_edges[12][3] = {
+      {1, 4, 5},   // line 0
+      {0, 4, 5},   // line 1
+      {3, 6, 7},   // line 2
+      {2, 6, 7},   // line 3
+      {0, 1, 5},   // line 4
+      {0, 1, 4},   // line 5
+      {2, 3, 7},   // line 6
+      {2, 3, 6},   // line 7
+      {9, 10, 11}, // line 8
+      {8, 10, 11}, // line 9
+      {8, 9, 11},  // line 10
+      {8, 9, 10}   // line 11
+    };
+
+
+    /**
+     * A structure that store the index of a cell and, crucially, how a
+     * given edge relates to this cell.
+     */
+    struct AdjacentCell
+    {
+      /**
+       * Default constructor. Initialize the fields with invalid values.
+       */
+      AdjacentCell()
+        : cell_index(numbers::invalid_unsigned_int)
+        , edge_within_cell(numbers::invalid_unsigned_int)
+      {}
+
+      /**
+       * Constructor. Initialize the fields with the given values.
+       */
+      AdjacentCell(const unsigned int cell_index,
+                   const unsigned int edge_within_cell)
+        : cell_index(cell_index)
+        , edge_within_cell(edge_within_cell)
+      {}
+
+
+      unsigned int cell_index;
+      unsigned int edge_within_cell;
+    };
+
+
+
+    template <int dim>
+    class AdjacentCells;
+
+    /**
+     * A class that represents all of the cells adjacent to a given edge.
+     * This class corresponds to the 2d case where each edge has at most
+     * two adjacent cells.
+     */
+    template <>
+    class AdjacentCells<2>
+    {
+    public:
+      /**
+       * An iterator that allows iterating over all cells adjacent
+       * to the edge represented by the current object.
+       */
+      using const_iterator = const AdjacentCell *;
+
+      /**
+       * Add the given cell to the collection of cells adjacent to
+       * the edge this object corresponds to. Since we are covering
+       * the 2d case, the set of adjacent cells currently
+       * represented by this object must have either zero or
+       * one element already, since we can not add more than two
+       * adjacent cells for each edge.
+       */
+      void
+      push_back(const AdjacentCell &adjacent_cell)
+      {
+        if (adjacent_cells[0].cell_index == numbers::invalid_unsigned_int)
+          adjacent_cells[0] = adjacent_cell;
+        else
+          {
+            Assert(adjacent_cells[1].cell_index ==
+                     numbers::invalid_unsigned_int,
+                   ExcInternalError());
+            adjacent_cells[1] = adjacent_cell;
+          }
+      }
+
+
+      /**
+       * Return an iterator to the first valid cell stored as adjacent to the
+       * edge represented by the current object.
+       */
+      const_iterator
+      begin() const
+      {
+        return adjacent_cells;
+      }
+
+
+      /**
+       * Return an iterator to the element past the last valid cell stored
+       * as adjacent to the edge represented by the current object.
+       * @return
+       */
+      const_iterator
+      end() const
+      {
+        // check whether the current object stores zero, one, or two
+        // adjacent cells, and use this to point to the element past the
+        // last valid one
+        if (adjacent_cells[0].cell_index == numbers::invalid_unsigned_int)
+          return adjacent_cells;
+        else if (adjacent_cells[1].cell_index == numbers::invalid_unsigned_int)
+          return adjacent_cells + 1;
+        else
+          return adjacent_cells + 2;
+      }
+
+    private:
+      /**
+       * References to the (at most) two cells that are adjacent to
+       * the edge this object corresponds to. Unused elements are
+       * default-initialized and have invalid values; in particular,
+       * their cell_index field equals numbers::invalid_unsigned_int.
+       */
+      AdjacentCell adjacent_cells[2];
+    };
+
+
+
+    /**
+     * 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
+     * indices of the two vertices, and the indices of the adjacent
+     * cells (together with a description *where* in each of the
+     * adjacent cells the edge is located). It also includes the
+     * (global) direction of the edge: either from the first vertex to
+     * the second, the other way around, or so far undetermined.
+     */
+    template <int dim>
+    class Edge
+    {
+    public:
+      /**
+       * Constructor. Create the edge based on the information given
+       * in @p cell, and selecting the edge with number @p edge_number
+       * within this cell. Initialize the edge as unoriented.
+       */
+      Edge(const CellData<dim> &cell, const unsigned int edge_number)
+        : orientation_status(not_oriented)
+      {
+        Assert(edge_number < GeometryInfo<dim>::lines_per_cell,
+               ExcInternalError());
+
+        // copy vertices for this particular line
+        vertex_indices[0] =
+          cell
+            .vertices[GeometryInfo<dim>::line_to_cell_vertices(edge_number, 0)];
+        vertex_indices[1] =
+          cell
+            .vertices[GeometryInfo<dim>::line_to_cell_vertices(edge_number, 1)];
+
+        // bring them into standard orientation
+        if (vertex_indices[0] > vertex_indices[1])
+          std::swap(vertex_indices[0], vertex_indices[1]);
+      }
+
+      /**
+       * Comparison operator for edges. It compares based on the
+       * lexicographic ordering of the two vertex indices.
+       */
+      bool
+      operator<(const Edge<dim> &e) const
+      {
+        return ((vertex_indices[0] < e.vertex_indices[0]) ||
+                ((vertex_indices[0] == e.vertex_indices[0]) &&
+                 (vertex_indices[1] < e.vertex_indices[1])));
+      }
+
+      /**
+       * Compare two edges for equality based on their vertex indices.
+       */
+      bool
+      operator==(const Edge<dim> &e) const
+      {
+        return ((vertex_indices[0] == e.vertex_indices[0]) &&
+                (vertex_indices[1] == e.vertex_indices[1]));
+      }
+
+      /**
+       * The global indices of the two vertices that bound this edge. These
+       * will be ordered so that the first index is less than the second.
+       */
+      unsigned int vertex_indices[2];
+
+      /**
+       * An enum that indicates the direction of this edge with
+       * regard to the two vertices that bound it.
+       */
+      enum OrientationStatus
+      {
+        not_oriented,
+        forward,
+        backward
+      };
+
+      OrientationStatus orientation_status;
+
+      /**
+       * Store the set of cells adjacent to this edge (these cells then
+       * also store *where* in the cell the edge is located).
+       */
+      AdjacentCells<dim> adjacent_cells;
+    };
+
+
+
+    /**
+     * A data structure that represents a cell with all of its vertices
+     * and edges.
+     */
+    template <int dim>
+    struct Cell
+    {
+      /**
+       * Construct a Cell object from a CellData object. Also take a
+       * (sorted) list of edges and to point the edges of the current
+       * object into this list of edges.
+       */
+      Cell(const CellData<dim> &c, const std::vector<Edge<dim>> &edge_list)
+      {
+        for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
+          vertex_indices[i] = c.vertices[i];
+
+        // now for each of the edges of this cell, find the location inside the
+        // given edge_list array and store than index
+        for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+          {
+            const Edge<dim> e(c, l);
+            edge_indices[l] =
+              (std::lower_bound(edge_list.begin(), edge_list.end(), e) -
+               edge_list.begin());
+            Assert(edge_indices[l] < edge_list.size(), ExcInternalError());
+            Assert(edge_list[edge_indices[l]] == e, ExcInternalError())
+          }
+      }
+
+      /**
+       * A list of global indices for the vertices that bound this cell.
+       */
+      unsigned int vertex_indices[GeometryInfo<dim>::vertices_per_cell];
+
+      /**
+       * A list of indices into the 'edge_list' array passed to the constructor
+       * for the edges of the current cell.
+       */
+      unsigned int edge_indices[GeometryInfo<dim>::lines_per_cell];
+    };
+
+
+
+    template <int dim>
+    class EdgeDeltaSet;
+
+    /**
+     * A class that represents by how much the set of parallel edges
+     * grows in each step. In the graph orientation paper, this set is
+     * called $\Delta_k$, thus the name.
+     *
+     * In 2d, this set can only include zero, one, or two elements.
+     * Consequently, the appropriate data structure is one in which
+     * we store at most 2 elements in a fixed sized data structure.
+     */
+    template <>
+    class EdgeDeltaSet<2>
+    {
+    public:
+      /**
+       * Iterator type for the elements of the set.
+       */
+      using const_iterator = const unsigned int *;
+
+      /**
+       * Default constructor. Initialize both slots as unused, corresponding
+       * to an empty set.
+       */
+      EdgeDeltaSet()
+      {
+        edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int;
+      }
+
+
+      /**
+       * Delete the elements of the set by marking both slots as unused.
+       */
+      void
+      clear()
+      {
+        edge_indices[0] = edge_indices[1] = numbers::invalid_unsigned_int;
+      }
+
+      /**
+       * Insert one element into the set. This will fail if the set already
+       * has two elements.
+       */
+      void
+      insert(const unsigned int edge_index)
+      {
+        if (edge_indices[0] == numbers::invalid_unsigned_int)
+          edge_indices[0] = edge_index;
+        else
+          {
+            Assert(edge_indices[1] == numbers::invalid_unsigned_int,
+                   ExcInternalError());
+            edge_indices[1] = edge_index;
+          }
+      }
+
+
+      /**
+       * Return an iterator pointing to the first element of the set.
+       */
+      const_iterator
+      begin() const
+      {
+        return edge_indices;
+      }
+
+
+      /**
+       * Return an iterator pointing to the element past the last used one.
+       */
+      const_iterator
+      end() const
+      {
+        // check whether the current object stores zero, one, or two
+        // indices, and use this to point to the element past the
+        // last valid one
+        if (edge_indices[0] == numbers::invalid_unsigned_int)
+          return edge_indices;
+        else if (edge_indices[1] == numbers::invalid_unsigned_int)
+          return edge_indices + 1;
+        else
+          return edge_indices + 2;
+      }
+
+    private:
+      /**
+       * Storage space to store the indices of at most two edges.
+       */
+      unsigned int edge_indices[2];
+    };
+
+
+
+    /**
+     * A class that represents by how much the set of parallel edges
+     * grows in each step. In the graph orientation paper, this set is
+     * called $\Delta_k$, thus the name.
+     *
+     * In 3d, this set can have arbitrarily many elements, unlike the
+     * 2d case specialized above. Consequently, we simply represent
+     * the data structure with a std::set. Class derivation ensures
+     * that we simply inherit all of the member functions of the
+     * base class.
+     */
+    template <>
+    class EdgeDeltaSet<3> : public std::set<unsigned int>
+    {};
+
+
+
+    /**
+     * From a list of cells, build a sorted vector that contains all of the
+     * edges that exist in the mesh.
+     */
+    template <int dim>
+    std::vector<Edge<dim>>
+    build_edges(const std::vector<CellData<dim>> &cells)
+    {
+      // build the edge list for all cells. because each cell has
+      // GeometryInfo<dim>::lines_per_cell edges, the total number
+      // of edges is this many times the number of cells. of course
+      // some of them will be duplicates, and we throw them out below
+      std::vector<Edge<dim>> edge_list;
+      edge_list.reserve(cells.size() * GeometryInfo<dim>::lines_per_cell);
+      for (unsigned int i = 0; i < cells.size(); ++i)
+        for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+          edge_list.emplace_back(cells[i], l);
+
+      // next sort the edge list and then remove duplicates
+      std::sort(edge_list.begin(), edge_list.end());
+      edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
+                      edge_list.end());
+
+      return edge_list;
+    }
+
+
+
+    /**
+     * Build the cell list. Update the edge array to let edges know
+     * which cells are adjacent to them.
+     */
+    template <int dim>
+    std::vector<Cell<dim>>
+    build_cells_and_connect_edges(const std::vector<CellData<dim>> &cells,
+                                  std::vector<Edge<dim>> &          edges)
+    {
+      std::vector<Cell<dim>> cell_list;
+      cell_list.reserve(cells.size());
+      for (unsigned int i = 0; i < cells.size(); ++i)
+        {
+          // create our own data structure for the cells and let it
+          // connect to the edges array
+          cell_list.emplace_back(cells[i], edges);
+
+          // then also inform the edges that they are adjacent
+          // to the current cell, and where within this cell
+          for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+            edges[cell_list.back().edge_indices[l]].adjacent_cells.push_back(
+              AdjacentCell(i, l));
+        }
+      Assert(cell_list.size() == cells.size(), ExcInternalError());
+
+      return cell_list;
+    }
+
+
+
+    /**
+     * Return the index within 'cells' of the first cell that has at least one
+     * edge that is not yet oriented.
+     */
+    template <int dim>
+    unsigned int
+    get_next_unoriented_cell(const std::vector<Cell<dim>> &cells,
+                             const std::vector<Edge<dim>> &edges,
+                             const unsigned int            current_cell)
+    {
+      for (unsigned int c = current_cell; c < cells.size(); ++c)
+        for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+          if (edges[cells[c].edge_indices[l]].orientation_status ==
+              Edge<dim>::not_oriented)
+            return c;
+
+      return numbers::invalid_unsigned_int;
+    }
+
+
+
+    /**
+     * Given a set of cells and edges, orient all edges that are
+     * (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<dim>> &cells,
+                                     std::vector<Edge<dim>> &      edges,
+                                     const unsigned int            cell,
+                                     const unsigned int            local_edge)
+    {
+      // 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
+      // implementation used till 2016, let us just choose the
+      // direction so that it matches what we have in the given cell.
+      //
+      // in fact, in what can only be assumed to be a bug in the
+      // original implementation, after orienting all edges, the code
+      // that rotates the cells so that they match edge orientations
+      // (see the rotate_cell() function below) rotated the cell two
+      // more times by 90 degrees. this is ok -- it simply flips all
+      // edge orientations, which leaves them valid. rather than do
+      // the same in the current implementation, we can achieve the
+      // 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 =
+          (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
+      else
+        {
+          Assert(
+            edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] ==
+              cells[cell].vertex_indices
+                [GeometryInfo<dim>::line_to_cell_vertices(local_edge, 1)],
+            ExcInternalError());
+          Assert(
+            edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] ==
+              cells[cell].vertex_indices
+                [GeometryInfo<dim>::line_to_cell_vertices(local_edge, 0)],
+            ExcInternalError());
+
+          // 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 =
+            (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
+        }
+
+      // walk outward from the given edge as described in
+      // the algorithm in the paper that documents all of
+      // this
+      //
+      // note that in 2d, each of the Deltas can at most
+      // contain two elements, whereas in 3d it can be arbitrarily many
+      EdgeDeltaSet<dim> Delta_k;
+      EdgeDeltaSet<dim> Delta_k_minus_1;
+      Delta_k_minus_1.insert(cells[cell].edge_indices[local_edge]);
+
+      while (Delta_k_minus_1.begin() !=
+             Delta_k_minus_1.end()) // while set is not empty
+        {
+          Delta_k.clear();
+
+          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 (typename AdjacentCells<dim>::const_iterator adjacent_cell =
+                     edges[*delta].adjacent_cells.begin();
+                   adjacent_cell != edges[*delta].adjacent_cells.end();
+                   ++adjacent_cell)
+                {
+                  const unsigned int K = adjacent_cell->cell_index;
+                  const unsigned int delta_is_edge_in_K =
+                    adjacent_cell->edge_within_cell;
+
+                  // figure out the direction of delta with respect to the cell
+                  // K (in the orientation in which the user has given it to us)
+                  const unsigned int first_edge_vertex =
+                    (edges[*delta].orientation_status == Edge<dim>::forward ?
+                       edges[*delta].vertex_indices[0] :
+                       edges[*delta].vertex_indices[1]);
+                  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>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
+                    ExcInternalError());
+
+                  // 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)
+                    {
+                      // 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());
+                        }
+                    }
+                }
+            }
+
+          // finally copy the new set to the previous one
+          // (corresponding to increasing 'k' by one in the
+          // algorithm)
+          Delta_k_minus_1 = Delta_k;
+        }
+    }
+
+
+    /**
+     * Given data structures @p cell_list and @p edge_list, where
+     * all edges are already oriented, rotate the cell with
+     * index @p cell_index in such a way that its local coordinate
+     * 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<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 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<dim>::not_oriented,
+                 ExcInternalError());
+          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 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,
+                               starting_vertex_of_edge +
+                                 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;
+            }
+
+          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:
+            {
+              // 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 (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+                temp_vertex_indices[v] =
+                  raw_cells[cell_index]
+                    .vertices[cube_permutations[origin_vertex_of_cell][v]];
+              for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+                raw_cells[cell_index].vertices[v] = temp_vertex_indices[v];
+
+              break;
+            }
+
+          default:
+            {
+              Assert(false, ExcNotImplemented());
+            }
+        }
+    }
+
+
+    /**
+     * Given a set of cells, find globally unique edge orientations
+     * and then rotate cells so that the coordinate system of the cell
+     * coincides with the coordinate systems of the adjacent edges.
+     */
+    template <int dim>
+    void
+    reorient(std::vector<CellData<dim>> &cells)
+    {
+      // first build the arrays that connect cells to edges and the other
+      // way around
+      std::vector<Edge<dim>> edge_list = build_edges(cells);
+      std::vector<Cell<dim>> cell_list =
+        build_cells_and_connect_edges(cells, edge_list);
+
+      // then loop over all cells and start orienting parallel edge sets
+      // of cells that still have non-oriented edges
+      unsigned int next_cell_with_unoriented_edge = 0;
+      while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
+                cell_list, edge_list, next_cell_with_unoriented_edge)) !=
+             numbers::invalid_unsigned_int)
+        {
+          // 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 (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)
+              orient_one_set_of_parallel_edges(
+                cell_list,
+                edge_list,
+                next_cell_with_unoriented_edge,
+                ParallelEdges<dim>::starter_edges[l]);
+
+          // ensure that we have really oriented all edges now, not just
+          // the starter edges
+          for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+            Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
+                               .edge_indices[l]]
+                       .orientation_status != Edge<dim>::not_oriented,
+                   ExcInternalError());
+        }
+
+      // now that we have oriented all edges, we need to rotate cells
+      // so that the edges point in the right direction with the now
+      // rotated coordinate system
+      for (unsigned int c = 0; c < cells.size(); ++c)
+        rotate_cell(cell_list, edge_list, c, cells);
+    }
+
+
+    // overload of the function above for 1d -- there is nothing
+    // to orient in that case
+    void reorient(std::vector<CellData<1>> &)
+    {}
+  } // namespace
+
+  template <int dim>
+  void
+  consistently_order_cells(std::vector<CellData<dim>> &cells)
+  {
+    Assert(cells.size() != 0,
+           ExcMessage(
+             "List of elements to orient must have at least one cell"));
+
+    // there is nothing for us to do in 1d
+    if (dim == 1)
+      return;
+
+    // check if grids are already consistent. if so, do
+    // nothing. if not, then do the reordering
+    if (!is_consistent(cells))
+      try
+        {
+          reorient(cells);
+        }
+      catch (const ExcMeshNotOrientable &)
+        {
+          // the mesh is not orientable. this is acceptable if we are in 3d,
+          // as class Triangulation knows how to handle this, but it is
+          // not in 2d; in that case, re-throw the exception
+          if (dim < 3)
+            throw;
+        }
+  }
+
+
   // define some transformations
   namespace internal
   {
@@ -4373,7 +5391,7 @@ namespace GridTools
     GridTools::delete_unused_vertices(vertices,
                                       cells_to_add,
                                       subcelldata_to_add);
-    GridReordering<dim, spacedim>::reorder_cells(cells_to_add, true);
+    GridTools::consistently_order_cells(cells_to_add);
 
     // Save manifolds
     auto manifold_ids = tria.get_manifold_ids();
index 145b5922a9f8c32162521240529e025a48c905c9..86091f40d644449445b636136358e4060dc934a3 100644 (file)
@@ -250,6 +250,11 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
         const std::vector<Point<deal_II_space_dimension>> &,
         std::vector<CellData<deal_II_dimension>> &);
 
+#  if deal_II_dimension == deal_II_space_dimension
+      template void
+      consistently_order_cells(std::vector<CellData<deal_II_dimension>> &);
+#  endif
+
       template void
       shift<deal_II_dimension>(
         const Tensor<1, deal_II_space_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.