From b44a5e729cac287f81c0261be6e7425ef5442faa Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 1 Jul 2015 11:09:26 -0500 Subject: [PATCH] Reimplement the 2d implementation of edge orientation. The old implementation, written in 2003 by Mike Anderson, is quite concise and of high quality. But it is not extensively documented, and is unnecessarily hard to read because it uses different conventions for names etc than we use in the rest of the library. This reimplementation uses standard conventions, and it uses the language and symbols of a paper about to be submitted that concisely describes the algorithms and data structures used in this algorithm. As an additional benefit, it uses the standard vertex ordering of cells rather than the old-style ordering used in the existing algorithm. Following the algorithm without understanding the graph theoretical context of the problem may be difficult, but because it uses the same symbols as in the paper, should be easy enough if you have the paper. --- .../deal.II/grid/grid_reordering_internal.h | 267 ----- source/grid/grid_reordering.cc | 1044 +++++++++-------- 2 files changed, 567 insertions(+), 744 deletions(-) diff --git a/include/deal.II/grid/grid_reordering_internal.h b/include/deal.II/grid/grid_reordering_internal.h index e053e650e8..7c8d8c423a 100644 --- a/include/deal.II/grid/grid_reordering_internal.h +++ b/include/deal.II/grid/grid_reordering_internal.h @@ -27,273 +27,6 @@ DEAL_II_NAMESPACE_OPEN namespace internal { - /** - * Implement the algorithm described in the documentation of the - * GridReordering<2> class. - * - * @author Michael Anderson, 2003 - */ - namespace GridReordering2d - { - - /** - * Check whether a given arrangement of cells is already consistent. If - * this is the case, then we skip the reordering pass. - * - * This function works by looping over all cells, checking whether one of - * its faces already exists in a list of edges, and if it already exists - * in reverse order, then return @p false. If it is not already in the - * list, or in the correct direction, then go on with the next faces or - * cell. - */ - bool - is_consistent (const std::vector > &cells); - - - /** - * Defines a variety of variables related to the connectivity of a simple - * quad element. This includes the nodes on each edge, which edges come - * into each node and what the default deal.II directions are for the - * quad. - * - * @verbatim - * s2 - * - * +-->--+ - * |3 2| - * s3 ^ ^ s1 - * |0 1| - * +-->--+ - * - * s0 - * @endverbatim - * - * @author Michael Anderson, 2003 - */ - class ConnectGlobals - { - public: - /** - * The nodes on each edge in anti-clockwise order { - * {0,1},{1,2},{2,3},{3,0} } - */ - static const int EdgeToNode[4][2]; - - /** - * The edges common to each node, in anti-clockwise order { - * {3,0},{0,1},{1,2},{2,3} } - */ - static const int NodeToEdge[4][2]; - - /** - * The nodes on each edge in "default direction order". - * {{0,1},{1,2},{3,2},{0,3}} - */ - static const int DefaultOrientation[4][2]; - }; - - - /** - * An enriched quad with information about how the mesh fits together so - * that we can move around the mesh efficiently. - * - * @author Michael Anderson, 2003 - */ - class MQuad - { - public: - /** - * v0 - v3 are indexes of the vertices of the quad, s0 - s3 are indexes - * for the sides of the quad - */ - MQuad (const unsigned int v0, - const unsigned int v1, - const unsigned int v2, - const unsigned int v3, - const unsigned int s0, - const unsigned int s1, - const unsigned int s2, - const unsigned int s3, - const CellData<2> &cd); - - /** - * Stores the vertex numbers - */ - unsigned int v[4]; - /** - * Stores the side numbers - */ - unsigned int side[4]; - - /** - * Copy of the @p CellData object from which we construct the data of - * this object. - */ - CellData<2> original_cell_data; - }; - - /** - * The enriched side class containing connectivity information. - * Orientation is from v0 to v1; Initially this should have v0 > &quads); - private: - - /** - * Sets up the internal data structures so that the we can do side - * hopping and face switching efficiently. This means we need a whole - * bunch of connectivity information - */ - void build_graph (const std::vector > &inquads); - - /** - * Orient the internal data into deal.II format The orientation - * algorithm is as follows - * - * 1) Find an unoriented quad (A) - * - * 2) Orient an un_oriented side (s) of (A) - * - * 3) side hop on (s) of (A) to get (B) - * - * 4) if opposite side to (s) of (B) is unoriented orient it - * - * 5) repeat 3) and 4) until side-hopping fails (we've reached a - * boundary) or (s) has already been oriented (we've closed a loop or - * unoriented sides). - * - * 6) Repeat 2), 3), 4) and 5) on other unoriented sides of (A) - * - * 7) Choose a new unoriented A. - */ - void orient(); - - /** - * Get the (now correctly oriented if we've called orient) quads. - */ - void get_quads(std::vector > &outquads) const; - - /** - * Orient_side(qnum,lsn) orients the local side lsn of the quad qnum in - * the triangulation. If the side opposite lsn is oriented then lsn is - * oriented to match it. Otherwise it is oriented in the "default" - * direction for the quad. - */ - void orient_side (const unsigned int quadnum, - const unsigned int localsidenum); - - /** - * Return true if all sides of the quad quadnum are oriented. - */ - bool is_fully_oriented_quad (const unsigned int quadnum) const; - - /** - * Return true if the side lsn of the quad quadnum is oriented. - */ - bool is_oriented_side (const unsigned int quadnum, - const unsigned int lsn) const; - - /** - * Return true is the side is oriented in the "default" direction - */ - bool is_side_default_oriented (const unsigned int qnum, - const unsigned int lsn) const; - - /** - * Increases UnOrQLoc from it's original value to the next quad with an - * unoriented side. Returns true if there was another unoriented quad. - */ - bool get_unoriented_quad (unsigned int &UnOrQLoc) const; - - /** - * Sets sidenum to the local sidenumber of an unoriented side of the - * quad quadnum. Returns true if such a side exists. - */ - bool get_unoriented_side (const unsigned int quadnum, - unsigned int &sidenum) const; - - /** - * side_hop(&qnum, &lsn) has qnum being the quadnumber of a quad in the - * triangulation, and a local side number. side_hop then sets qnum to - * the quadnumber across the other side of the side, and sets lsn so - * that quads[qnum].sides[lsn] is the same before and after the call. If - * there is no other quad on the other side of the current quad, then - * side_hop returns false. - */ - bool side_hop (unsigned int &qnum, - unsigned int &lsn) const; - - /** - * A list of enriched sides/edges of the mesh. - */ - std::vector sides; - /** - * A list of enriched quads in the mesh. - */ - std::vector mquads; - }; - } // namespace GridReordering2d - - /** * Implement the algorithm described in the documentation of the * GridReordering<3> class. diff --git a/source/grid/grid_reordering.cc b/source/grid/grid_reordering.cc index 6bdfed9d34..b06963a894 100644 --- a/source/grid/grid_reordering.cc +++ b/source/grid/grid_reordering.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2000 - 2015 by the deal.II authors +// Copyright (C) 2000 - 2016 by the deal.II authors // // This file is part of the deal.II library. // @@ -29,585 +29,675 @@ DEAL_II_NAMESPACE_OPEN -template<> -void -GridReordering<1>::reorder_cells (std::vector > &, - const bool) -{ - // there should not be much to do - // in 1d... -} - - -template<> -void -GridReordering<1>::invert_all_cells_of_negative_grid(const std::vector > &, - std::vector > &) -{ - // nothing to be done in 1d -} - -template<> -void -GridReordering<1,2>::reorder_cells (std::vector > &, - const bool) -{ - // there should not be much to do - // in 1d... -} - - -template<> -void -GridReordering<1,2>::invert_all_cells_of_negative_grid(const std::vector > &, - std::vector > &) -{ - // nothing to be done in 1d -} - - -template<> -void -GridReordering<1,3>::reorder_cells (std::vector > &, - const bool) -{ - // there should not be much to do - // in 1d... -} - - -template<> -void -GridReordering<1,3>::invert_all_cells_of_negative_grid(const std::vector > &, - std::vector > &) -{ - // nothing to be done in 1d -} - - namespace internal { namespace GridReordering2d { -// -- Definition of connectivity information -- - const int ConnectGlobals::EdgeToNode[4][2] = - { {0,1},{1,2},{2,3},{3,0} }; - - const int ConnectGlobals::NodeToEdge[4][2] = - { {3,0},{0,1},{1,2},{2,3} }; - - const int ConnectGlobals::DefaultOrientation[4][2] = - {{0,1},{1,2},{3,2},{0,3}}; - - /** - * Simple data structure denoting - * an edge, i.e. the ordered pair - * of its vertices. This is only - * used in the is_consistent + * 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 Edge + struct CheapEdge { - Edge (const unsigned int v0, - const unsigned int v1) + /** + * Construct an edge from the global indices of its two vertices. + */ + CheapEdge (const unsigned int v0, + const unsigned int v1) : v0(v0), v1(v1) {} - const unsigned int v0, v1; - bool operator < (const Edge &e) const + /** + * 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 bool - is_consistent (const std::vector > &cells) + is_consistent (const std::vector > &cells) { - std::set edges; + std::set edges; - std::vector >::const_iterator c = cells.begin(); - for (; c != cells.end(); ++c) + for (typename std::vector >::const_iterator c = cells.begin(); + c != cells.end(); ++c) { - // construct the four edges - // in reverse order - const Edge reverse_edges[4] = { Edge (c->vertices[1], - c->vertices[0]), - Edge (c->vertices[2], - c->vertices[1]), - Edge (c->vertices[2], - c->vertices[3]), - Edge (c->vertices[3], - c->vertices[0]) - }; - // for each of them, check - // whether they are already - // in the set - if ((edges.find (reverse_edges[0]) != edges.end()) || - (edges.find (reverse_edges[1]) != edges.end()) || - (edges.find (reverse_edges[2]) != edges.end()) || - (edges.find (reverse_edges[3]) != edges.end())) - return false; - // ok, not. insert them - // in the order in which - // we want them - // (std::set eliminates - // duplicated by itself) - for (unsigned int i = 0; i<4; ++i) + // 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::lines_per_cell; ++l) { - const Edge e(reverse_edges[i].v1, reverse_edges[i].v0); - edges.insert (e); + const CheapEdge reverse_edge (c->vertices[GeometryInfo::line_to_cell_vertices(l, 1)], + c->vertices[GeometryInfo::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::line_to_cell_vertices(l, 0)], + c->vertices[GeometryInfo::line_to_cell_vertices(l, 1)]); + edges.insert (correct_edge); } - // then go on with next - // cell } - // no conflicts found, so - // return true + + // no conflicts found, so return true return true; } - /** - * Returns an MSide corresponding to the - * specified side of a deal.II CellData<2> object. + * 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. */ - MSide quadside(const CellData<2> &q, unsigned int i) + template + struct ParallelEdges { - Assert (i<4, ExcInternalError()); - return MSide(q.vertices[ConnectGlobals::EdgeToNode[i][0]], - q.vertices[ConnectGlobals::EdgeToNode[i][1]]); - } + static const unsigned int starter_edges[dim]; + static const unsigned int parallel_edges[GeometryInfo::lines_per_cell][(1<<(dim-1)) - 1]; + }; + 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 + }; + + + template + struct Edge; /** - * Wrapper class for the quadside() function + * 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. */ - struct QuadSide + template <> + struct Edge<2> { - typedef CellData<2> first_argument_type; - typedef int second_argument_type; - typedef MSide result_type; + static const unsigned int dim = 2; - MSide operator()(const CellData<2> &q, int i) const + /** + * Default constructor. Creates an invalid edge. + */ + Edge () + : + orientation_status (not_oriented) { - return quadside(q,i); + for (unsigned int i=0; i<2; ++i) + vertex_indices[i] = numbers::invalid_unsigned_int; + + for (unsigned int i=0; i<2; ++i) + { + static const AdjacentCell invalid_cell = { numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int + }; + adjacent_cells[i] = invalid_cell; + } } - }; + /** + * 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 &cell, + const unsigned int edge_number) + : + orientation_status (not_oriented) + { + Assert (edge_number < GeometryInfo::lines_per_cell, ExcInternalError()); + // copy vertices for this particular line + vertex_indices[0] = cell.vertices[GeometryInfo::line_to_cell_vertices(edge_number, 0)]; + vertex_indices[1] = cell.vertices[GeometryInfo::line_to_cell_vertices(edge_number, 1)]; - MQuad::MQuad (const unsigned int v0, - const unsigned int v1, - const unsigned int v2, - const unsigned int v3, - const unsigned int s0, - const unsigned int s1, - const unsigned int s2, - const unsigned int s3, - const CellData<2> &cd) - : - original_cell_data (cd) - { - v[0] = v0; - v[1] = v1; - v[2] = v2; - v[3] = v3; - side[0] = s0; - side[1] = s1; - side[2] = s2; - side[3] = s3; - } + // bring them into standard orientation + if (vertex_indices[0] > vertex_indices[1]) + std::swap (vertex_indices[0], vertex_indices[1]); + for (unsigned int i=0; i<2; ++i) + { + static const AdjacentCell invalid_cell = { numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int + }; + adjacent_cells[i] = invalid_cell; + } + } - MSide::MSide (const unsigned int initv0, - const unsigned int initv1) - : - v0(initv0), v1(initv1), - Q0(numbers::invalid_unsigned_int), - Q1(numbers::invalid_unsigned_int), - lsn0(numbers::invalid_unsigned_int), - lsn1(numbers::invalid_unsigned_int), - Oriented(false) - {} + /** + * Comparison operator for edges. It compares based on the + * lexicographic ordering of the two vertex indices. + */ + bool operator< (const Edge &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 &e) const + { + return ((vertex_indices[0] == e.vertex_indices[0]) + && + (vertex_indices[1] == e.vertex_indices[1])); + } + /** + * Store the given cell index as one of the cells as adjacent to + * the current edge. Also store that the current edge has index + * @p edge_within_cell within the given cell. + */ + void add_adjacent_cell (const unsigned int cell_index, + const unsigned int edge_within_cell) + { + const AdjacentCell adjacent_cell = { cell_index, edge_within_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; + } + } - bool - MSide::operator == (const MSide &s2) const - { - if ((v0 == s2.v0)&&(v1 == s2.v1)) - { - return true; - } - if ((v0 == s2.v1)&&(v1 == s2.v0)) - { - return true; - } - return false; - } + /** + * 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 + }; - bool - MSide::operator != (const MSide &s2) const - { - return !(*this == s2); - } + OrientationStatus orientation_status; + /** + * A structure that stores how this edge relates to the at most + * two cells that are next to it. + */ + struct AdjacentCell + { + unsigned int cell_index; + unsigned int edge_within_cell; + }; + + AdjacentCell adjacent_cells[2]; + }; - namespace + + + /** + * A data structure that represents a cell with all of its vertices + * and edges. + */ + template + struct Cell { - void side_rectify (MSide &s) + /** + * Default construct a cell. + */ + Cell () { - if (s.v0>s.v1) - std::swap (s.v0, s.v1); + for (unsigned int i=0; i::vertices_per_cell; ++i) + vertex_indices[i] = numbers::invalid_unsigned_int; + for (unsigned int i=0; i::lines_per_cell; ++i) + edge_indices[i] = numbers::invalid_unsigned_int; } - bool side_sort_less(const MSide &s1, const MSide &s2) + /** + * Construct a Cell object from a CellData object. Also take a (sorted) + * list of edges and to point into from the current object. + * @param c + * @param edge_list + */ + Cell (const CellData &c, + const std::vector > &edge_list) { - int s1vmin,s1vmax; - int s2vmin,s2vmax; - if (s1.v0::vertices_per_cell; ++i) + 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::lines_per_cell; ++l) { - s2vmin = s2.v1; - s2vmax = s2.v0; + const Edge 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()) } - - if (s1vmins2vmin) - return false; - return s1vmax &q, - const std::vector &elist) - { - // compute the indices of the four - // sides that bound this quad. note - // that the incoming list elist is - // sorted with regard to the - // side_sort_less criterion - unsigned int edges[4] = { numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int, - numbers::invalid_unsigned_int - }; + unsigned int vertex_indices[GeometryInfo::vertices_per_cell]; - for (unsigned int i=0; i<4; ++i) - edges[i] = (Utilities::lower_bound (elist.begin(), - elist.end(), - quadside(q,i), - side_sort_less) - - - elist.begin()); - - return MQuad(q.vertices[0],q.vertices[1], q.vertices[2], q.vertices[3], - edges[0], edges[1], edges[2], edges[3], - q); - } - } + /** + * 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::lines_per_cell]; + }; - void - GridReordering::reorient(std::vector > &quads) + + /** + * From a list of cells, build a sorted vector that contains all of the edges + * that exist in the mesh. + */ + template + std::vector > + build_edges (const std::vector > &cells) { - build_graph(quads); - orient(); - get_quads(quads); + // build the edge list for all cells. because each cell has + // GeometryInfo::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_list; + edge_list.reserve(cells.size()*GeometryInfo::lines_per_cell); + for (unsigned int i=0; i::lines_per_cell; ++l) + edge_list.push_back (Edge(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; } - - void - GridReordering::build_graph (const std::vector > &inquads) + /** + * Build the cell list. Update the edge array to let edges know + * which cells are adjacent to them. + */ + template + std::vector > + build_cells_and_connect_edges (const std::vector > &cells, + std::vector > &edges) { - //Reserve some space - sides.reserve(4*inquads.size()); - - //Insert all the sides into the side vector - for (int i = 0; i<4; ++i) + std::vector > cell_list; + cell_list.reserve(cells.size()); + for (unsigned int i=0; i(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::lines_per_cell; ++l) + edges[cell_list.back().edge_indices[l]].add_adjacent_cell (i, l); } + Assert (cell_list.size() == cells.size(), ExcInternalError()); - //Change each edge so that v0(sides).swap(sides); - - // Now assign the correct sides to - // each quads - mquads.reserve(inquads.size()); - std::transform(inquads.begin(), - inquads.end(), - std::back_inserter(mquads), - std_cxx11::bind(build_quad_from_vertices, - std_cxx11::_1, - std_cxx11::cref(sides)) ); - - // Assign the quads to their sides also. - int qctr = 0; - for (std::vector::iterator it = mquads.begin(); it != mquads.end(); ++it) - { - for (unsigned int i = 0; i<4; ++i) - { - MSide &ss = sides[(*it).side[i]]; - if (ss.Q0 == numbers::invalid_unsigned_int) - { - ss.Q0 = qctr; - ss.lsn0 = i; - } - else if (ss.Q1 == numbers::invalid_unsigned_int) - { - ss.Q1 = qctr; - ss.lsn1 = i; - } - else - AssertThrow (false, ExcInternalError()); - } - qctr++; - } + return cell_list; } - void GridReordering::orient() + /** + * Return the index within 'cells' of the first cell that has at least one + * edge that is not yet oriented. + */ + template + unsigned int + get_next_unoriented_quad(const std::vector > &cells, + const std::vector > &edges) { - // do what the comment in the - // class declaration says - unsigned int qnum = 0; - while (get_unoriented_quad(qnum)) + for (unsigned int c=0; c::lines_per_cell; ++l) + if (edges[cells[c].edge_indices[l]].orientation_status == Edge::not_oriented) + return c; + + return numbers::invalid_unsigned_int; + } + + + /** + * Given a set of cells and edges, orient all edges that are + * (globall) parallel to the one identified by the @p cell and + * within it the one with index @p local_edge. + */ + void + orient_one_set_of_parallel_edges (const std::vector > &cells, + std::vector > &edges, + const unsigned int cell, + const unsigned int local_edge) + { + const unsigned int dim = 2; + + // choose the direction of the first edge by chance + edges[cells[cell].edge_indices[local_edge]].orientation_status = Edge::forward; + + // 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. we indicate non-used elements + // of these sets by invalid unsigned ints; if the set has + // only one element, then we use the first + unsigned int Delta_k[2] = { numbers::invalid_unsigned_int, + numbers::invalid_unsigned_int + }; + unsigned int Delta_k_minus_1[2] = { cells[cell].edge_indices[local_edge], + numbers::invalid_unsigned_int + }; + while (Delta_k_minus_1[0] != numbers::invalid_unsigned_int) // while set is not empty { - unsigned int lsn = 0; - while (get_unoriented_side(qnum,lsn)) - { - orient_side(qnum,lsn); - unsigned int qqnum = qnum; - while (side_hop(qqnum,lsn)) - { - // switch this face - lsn = (lsn+2)%4; - if (!is_oriented_side(qqnum,lsn)) - orient_side(qqnum,lsn); - else - //We've found a - //cycle.. and - //oriented all - //quads in it. - break; - } - } + Delta_k[0] = Delta_k[1] = numbers::invalid_unsigned_int; + + for (unsigned int delta_element=0; delta_element<2; ++delta_element) + if (Delta_k_minus_1[delta_element] != numbers::invalid_unsigned_int) + { + // get the edge we are currently looking at. it must already + // have been oriented + const unsigned int delta = Delta_k_minus_1[delta_element]; + Assert (edges[delta].orientation_status != Edge::not_oriented, + ExcInternalError()); + + // now go through the cells adjacent to this edge + for (unsigned int K_element=0; K_element<2; ++K_element) + if (edges[delta].adjacent_cells[K_element].cell_index != numbers::invalid_unsigned_int) + { + const unsigned int K = edges[delta].adjacent_cells[K_element].cell_index; + const unsigned int delta_is_edge_in_K = edges[delta].adjacent_cells[K_element].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::forward + ? + edges[delta].vertex_indices[0] + : + edges[delta].vertex_indices[1]); + const unsigned int first_edge_vertex_in_K = cells[K].vertex_indices[GeometryInfo::face_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::face_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::face_to_cell_vertices( + ParallelEdges::parallel_edges[delta_is_edge_in_K][0], + (first_edge_vertex == first_edge_vertex_in_K + ? + 0 + : + 1))]; + + const Edge::OrientationStatus opposite_edge_orientation + = (edges[opposite_edge].vertex_indices[0] + == + first_opposite_edge_vertex + ? + Edge::forward + : + Edge::backward); + + // see if the opposite edge (there is only one in 2d) has already been + // oriented. + if (edges[opposite_edge].orientation_status == Edge::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; + if (Delta_k[0] == numbers::invalid_unsigned_int) + Delta_k[0] = opposite_edge; + else + { + Assert (Delta_k[1] == numbers::invalid_unsigned_int, ExcInternalError()); + Delta_k[1] = 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()); + } + } + } + + // finally copy the new set to the previous one + // (corresponding to increasing 'k' by one in the + // algorithm) + Delta_k_minus_1[0] = Delta_k[0]; + Delta_k_minus_1[1] = Delta_k[1]; } } + /** + * 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 raw_cells[cell_index]. + */ void - GridReordering::orient_side(const unsigned int quadnum, - const unsigned int localsidenum) + rotate_cell (const std::vector > &cell_list, + const std::vector > &edge_list, + const unsigned int cell_index, + std::vector > &raw_cells) { - MQuad &quad = mquads[quadnum]; - int op_side_l = (localsidenum+2)%4; - MSide &side = sides[mquads[quadnum].side[localsidenum]]; - const MSide &op_side = sides[mquads[quadnum].side[op_side_l]]; - - //is the opposite side oriented? - if (op_side.Oriented) + // 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) { - //YES - Make the orientations match - //Is op side in default orientation? - if (op_side.v0 == quad.v[ConnectGlobals::DefaultOrientation[op_side_l][0]]) - { - //YES - side.v0 = quad.v[ConnectGlobals::DefaultOrientation[localsidenum][0]]; - side.v1 = quad.v[ConnectGlobals::DefaultOrientation[localsidenum][1]]; - } + Assert (edge_list[cell_list[cell_index].edge_indices[e]].orientation_status + != Edge<2>::not_oriented, + ExcInternalError()); + if (edge_list[cell_list[cell_index].edge_indices[e]].orientation_status == Edge<2>::forward) + starting_vertex_of_edge[e] = edge_list[cell_list[cell_index].edge_indices[e]].vertex_indices[0]; else - { - //NO, its reversed - side.v0 = quad.v[ConnectGlobals::DefaultOrientation[localsidenum][1]]; - side.v1 = quad.v[ConnectGlobals::DefaultOrientation[localsidenum][0]]; - } + 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()); + + // 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) { - //NO - //Just use the default orientation - side.v0 = quad.v[ConnectGlobals::DefaultOrientation[localsidenum][0]]; - side.v1 = quad.v[ConnectGlobals::DefaultOrientation[localsidenum][1]]; + 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; } - side.Oriented = true; } - - bool - GridReordering::is_fully_oriented_quad(const unsigned int quadnum) const - { - return ( - (sides[mquads[quadnum].side[0]].Oriented)&& - (sides[mquads[quadnum].side[1]].Oriented)&& - (sides[mquads[quadnum].side[2]].Oriented)&& - (sides[mquads[quadnum].side[3]].Oriented) - ); - } - - - - bool - GridReordering::is_oriented_side(const unsigned int quadnum, - const unsigned int lsn) const + /** + * 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 + void reorient (std::vector > &cells) { - return (sides[mquads[quadnum].side[lsn]].Oriented); - } - - - + // first build the arrays that connect cells to edges and the other + // way around + std::vector > edge_list = build_edges(cells); + std::vector > 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; + while ((next_cell_with_unoriented_edge = get_next_unoriented_quad(cell_list, edge_list)) != + 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. there are only + // dim independent sets of edges + for (unsigned int l=0; l::starter_edges[l]]].orientation_status + == Edge::not_oriented) + orient_one_set_of_parallel_edges (cell_list, + edge_list, + next_cell_with_unoriented_edge, + ParallelEdges::starter_edges[l]); + + // ensure that we have really oriented all edges now, not just + // the starter edges + for (unsigned int l=0; l::lines_per_cell; ++l) + Assert (edge_list[cell_list[next_cell_with_unoriented_edge].edge_indices[l]].orientation_status + != Edge::not_oriented, + ExcInternalError()); + } - bool - GridReordering::get_unoriented_quad(unsigned int &UnOrQLoc) const - { - while ( (UnOrQLoc +void +GridReordering<1>::reorder_cells (std::vector > &, + const bool) +{ + // there should not be much to do + // in 1d... +} - bool - GridReordering::side_hop (unsigned int &qnum, unsigned int &lsn) const - { - const MQuad &mq = mquads[qnum]; - const MSide &s = sides[mq.side[lsn]]; - unsigned int opquad = 0; - if (s.Q0 == qnum) - { - opquad = s.Q1; - lsn = s.lsn1; - } - else - { - opquad = s.Q0; - lsn = s.lsn0; - } +template<> +void +GridReordering<1>::invert_all_cells_of_negative_grid(const std::vector > &, + std::vector > &) +{ + // nothing to be done in 1d +} - if (opquad != numbers::invalid_unsigned_int) - { - qnum = opquad; - return true; - } +template<> +void +GridReordering<1,2>::reorder_cells (std::vector > &, + const bool) +{ + // there should not be much to do + // in 1d... +} - return false; - } +template<> +void +GridReordering<1,2>::invert_all_cells_of_negative_grid(const std::vector > &, + std::vector > &) +{ + // nothing to be done in 1d +} - void - GridReordering::get_quads (std::vector > &outquads) const - { - outquads.clear(); - outquads.reserve(mquads.size()); - for (unsigned int qn = 0; qn q = mquads[qn].original_cell_data; - - // Are the sides oriented? - Assert (is_fully_oriented_quad(qn), ExcInternalError()); - bool s[4]; //whether side 1 ,2, 3, 4 are in the default orientation - for (int sn = 0; sn<4; sn++) - { - s[sn] = is_side_default_oriented(qn,sn); - } - // Are they oriented in the "deal way"? - Assert (s[0] == s[2], ExcInternalError()); - Assert (s[1] == s[3], ExcInternalError()); - // How much we rotate them by. - int rotn = 2*(s[0]?1:0)+ ((s[0]^s[1])?1:0); - for (int i = 0; i<4; ++i) - { - q.vertices[(i+rotn)%4] = mquads[qn].v[i]; - } - outquads.push_back(q); - } +template<> +void +GridReordering<1,3>::reorder_cells (std::vector > &, + const bool) +{ + // there should not be much to do + // in 1d... +} - } - bool - GridReordering::is_side_default_oriented (const unsigned int qnum, - const unsigned int lsn) const - { - return (sides[mquads[qnum].side[lsn]].v0 == - mquads[qnum].v[ConnectGlobals::DefaultOrientation[lsn][0]]); - } - } // namespace GridReordering2d -} // namespace internal +template<> +void +GridReordering<1,3>::invert_all_cells_of_negative_grid(const std::vector > &, + std::vector > &) +{ + // nothing to be done in 1d +} // anonymous namespace for internal helper functions @@ -676,21 +766,21 @@ void GridReordering<2>::reorder_cells (std::vector > &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); // check if grids are already // consistent. if so, do // nothing. if not, then do the // reordering if (!internal::GridReordering2d::is_consistent (cells)) - internal::GridReordering2d::GridReordering().reorient(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); } -- 2.39.5