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