From: David Wells Date: Thu, 22 Apr 2021 15:48:15 +0000 (-0400) Subject: Move GridReordering::reorder_cells(). X-Git-Tag: v9.3.0-rc1~189^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12082%2Fhead;p=dealii.git Move GridReordering::reorder_cells(). --- diff --git a/doc/news/changes/incompatibilities/20210422DavidWells b/doc/news/changes/incompatibilities/20210422DavidWells new file mode 100644 index 0000000000..d371ebdc05 --- /dev/null +++ b/doc/news/changes/incompatibilities/20210422DavidWells @@ -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. +
+(David Wells, 2021/04/22) diff --git a/include/deal.II/grid/grid_reordering.h b/include/deal.II/grid/grid_reordering.h index 6d7d76ac78..430784db7a 100644 --- a/include/deal.II/grid/grid_reordering.h +++ b/include/deal.II/grid/grid_reordering.h @@ -25,20 +25,16 @@ 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 -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> &original_cells, const bool use_new_style_ordering = false); diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 723ed0f802..f86d9d9882 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -460,6 +460,19 @@ namespace GridTools const std::vector> &all_vertices, std::vector> & 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 + void + consistently_order_cells(std::vector> &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 --------------- */ diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index f7e99e478c..1427244d7a 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -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::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 void subdivided_parallelepiped(Triangulation & 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::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::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::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(); diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc index 0ee5e9b82e..5fbf7c3aec 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -561,7 +561,7 @@ GridIn::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::read_vtk(std::istream &in) if (dim == spacedim) GridTools::invert_all_negative_measure_cells(vertices, cells); - GridReordering::reorder_cells(cells, true); + GridTools::consistently_order_cells(cells); tria->create_triangulation(vertices, cells, subcelldata); } else @@ -862,7 +862,7 @@ GridIn::read_unv(std::istream &in) if (dim == spacedim) GridTools::invert_all_negative_measure_cells(vertices, cells); - GridReordering::reorder_cells(cells, true); + GridTools::consistently_order_cells(cells); tria->create_triangulation(vertices, cells, subcelldata); } @@ -1095,7 +1095,7 @@ GridIn::read_ucd(std::istream &in, // ... and cells if (dim == spacedim) GridTools::invert_all_negative_measure_cells(vertices, cells); - GridReordering::reorder_cells(cells, true); + GridTools::consistently_order_cells(cells); tria->create_triangulation(vertices, cells, subcelldata); } @@ -1348,7 +1348,7 @@ GridIn::read_dbmesh(std::istream &in) GridTools::delete_unused_vertices(vertices, cells, subcelldata); // ...and cells GridTools::invert_all_negative_measure_cells(vertices, cells); - GridReordering::reorder_cells(cells, true); + GridTools::consistently_order_cells(cells); tria->create_triangulation(vertices, cells, subcelldata); } @@ -1416,7 +1416,7 @@ GridIn::read_xda(std::istream &in) GridTools::delete_unused_vertices(vertices, cells, subcelldata); // ... and cells GridTools::invert_all_negative_measure_cells(vertices, cells); - GridReordering::reorder_cells(cells, true); + GridTools::consistently_order_cells(cells); tria->create_triangulation(vertices, cells, subcelldata); } @@ -2096,7 +2096,7 @@ GridIn::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::read_msh(std::istream &in) // ... and cells if (dim == spacedim) GridTools::invert_all_negative_measure_cells(vertices, cells); - GridReordering::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::reorder_cells(cells, true); + GridTools::consistently_order_cells(cells); tria->create_triangulation(vertices, cells, subcelldata); } @@ -2954,7 +2954,7 @@ GridIn::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::reorder_cells(cells, true); + GridTools::consistently_order_cells(cells); tria->create_triangulation(vertices, cells, subcelldata); #else diff --git a/source/grid/grid_reordering.cc b/source/grid/grid_reordering.cc index 76f532e1d5..5dcc6c0157 100644 --- a/source/grid/grid_reordering.cc +++ b/source/grid/grid_reordering.cc @@ -29,984 +29,6 @@ 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 - bool - is_consistent(const std::vector> &cells) - { - std::set edges; - - for (typename std::vector>::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::lines_per_cell; ++l) - { - 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); - } - } - - // 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 - 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::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 - 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, 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 - * 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 - 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 &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)]; - - // 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 &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])); - } - - /** - * 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 adjacent_cells; - }; - - - - /** - * A data structure that represents a cell with all of its vertices - * and edges. - */ - template - 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 &c, const std::vector> &edge_list) - { - for (const unsigned int i : GeometryInfo::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::lines_per_cell; ++l) - { - 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()) - } - } - - /** - * A list of global indices for the vertices that bound this cell. - */ - unsigned int vertex_indices[GeometryInfo::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::lines_per_cell]; - }; - - - - template - 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 - {}; - - - - /** - * 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 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 < cells.size(); ++i) - for (unsigned int l = 0; l < GeometryInfo::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 - std::vector> - build_cells_and_connect_edges(const std::vector> &cells, - std::vector> & edges) - { - std::vector> 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::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 - unsigned int - get_next_unoriented_cell(const std::vector> &cells, - const std::vector> &edges, - const unsigned int current_cell) - { - for (unsigned int c = current_cell; c < cells.size(); ++c) - for (unsigned int l = 0; l < GeometryInfo::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 - * (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) - { - // 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::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::backward : Edge::forward); - else - { - Assert( - edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] == - cells[cell].vertex_indices[GeometryInfo::line_to_cell_vertices( - local_edge, 1)], - ExcInternalError()); - Assert( - edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] == - cells[cell].vertex_indices[GeometryInfo::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::forward : Edge::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 Delta_k; - EdgeDeltaSet 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::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 (typename AdjacentCells::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::forward ? - edges[*delta].vertex_indices[0] : - edges[*delta].vertex_indices[1]); - 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< - 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::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::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()); - } - } - } - } - - // 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 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) - { - // 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 < GeometryInfo::lines_per_cell; ++e) - { - Assert(edge_list[cell_list[cell_index].edge_indices[e]] - .orientation_status != Edge::not_oriented, - ExcInternalError()); - 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 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::vertices_per_cell; - ++origin_vertex_of_cell) - if (std::count(starting_vertex_of_edge, - starting_vertex_of_edge + - 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; - } - - 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::vertices_per_cell]; - for (const unsigned int v : GeometryInfo::vertex_indices()) - temp_vertex_indices[v] = - raw_cells[cell_index] - .vertices[cube_permutations[origin_vertex_of_cell][v]]; - for (const unsigned int v : GeometryInfo::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 - void - reorient(std::vector> &cells) - { - // 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 = 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::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 < GeometryInfo::lines_per_cell; ++l) - Assert( - edge_list[cell_list[next_cell_with_unoriented_edge].edge_indices[l]] - .orientation_status != Edge::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> &) - {} -} // namespace - - // anonymous namespace for internal helper functions namespace { @@ -1089,21 +111,7 @@ GridReordering::reorder_cells(std::vector> &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) diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 8b125a912f..f3561ae4b7 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -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 + bool + is_consistent(const std::vector> &cells) + { + std::set edges; + + for (typename std::vector>::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::lines_per_cell; ++l) + { + 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); + } + } + + // 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 + 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::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 + 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, 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 + * 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 + 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 &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)]; + + // 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 &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])); + } + + /** + * 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 adjacent_cells; + }; + + + + /** + * A data structure that represents a cell with all of its vertices + * and edges. + */ + template + 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 &c, const std::vector> &edge_list) + { + for (const unsigned int i : GeometryInfo::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::lines_per_cell; ++l) + { + 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()) + } + } + + /** + * A list of global indices for the vertices that bound this cell. + */ + unsigned int vertex_indices[GeometryInfo::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::lines_per_cell]; + }; + + + + template + 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 + {}; + + + + /** + * 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 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 < cells.size(); ++i) + for (unsigned int l = 0; l < GeometryInfo::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 + std::vector> + build_cells_and_connect_edges(const std::vector> &cells, + std::vector> & edges) + { + std::vector> 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::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 + unsigned int + get_next_unoriented_cell(const std::vector> &cells, + const std::vector> &edges, + const unsigned int current_cell) + { + for (unsigned int c = current_cell; c < cells.size(); ++c) + for (unsigned int l = 0; l < GeometryInfo::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 + * (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) + { + // 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::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::backward : Edge::forward); + else + { + Assert( + edges[cells[cell].edge_indices[local_edge]].vertex_indices[0] == + cells[cell].vertex_indices + [GeometryInfo::line_to_cell_vertices(local_edge, 1)], + ExcInternalError()); + Assert( + edges[cells[cell].edge_indices[local_edge]].vertex_indices[1] == + cells[cell].vertex_indices + [GeometryInfo::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::forward : Edge::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 Delta_k; + EdgeDeltaSet 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::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 (typename AdjacentCells::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::forward ? + edges[*delta].vertex_indices[0] : + edges[*delta].vertex_indices[1]); + 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< + 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::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::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::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()); + } + } + } + } + + // 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 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) + { + // 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 < GeometryInfo::lines_per_cell; ++e) + { + Assert(edge_list[cell_list[cell_index].edge_indices[e]] + .orientation_status != Edge::not_oriented, + ExcInternalError()); + 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 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::vertices_per_cell; + ++origin_vertex_of_cell) + if (std::count(starting_vertex_of_edge, + starting_vertex_of_edge + + 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; + } + + 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::vertices_per_cell]; + for (const unsigned int v : GeometryInfo::vertex_indices()) + temp_vertex_indices[v] = + raw_cells[cell_index] + .vertices[cube_permutations[origin_vertex_of_cell][v]]; + for (const unsigned int v : GeometryInfo::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 + void + reorient(std::vector> &cells) + { + // 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 = 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::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 < GeometryInfo::lines_per_cell; ++l) + Assert(edge_list[cell_list[next_cell_with_unoriented_edge] + .edge_indices[l]] + .orientation_status != Edge::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> &) + {} + } // namespace + + template + void + consistently_order_cells(std::vector> &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::reorder_cells(cells_to_add, true); + GridTools::consistently_order_cells(cells_to_add); // Save manifolds auto manifold_ids = tria.get_manifold_ids(); diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 145b5922a9..86091f40d6 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -250,6 +250,11 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const std::vector> &, std::vector> &); +# if deal_II_dimension == deal_II_space_dimension + template void + consistently_order_cells(std::vector> &); +# endif + template void shift( const Tensor<1, deal_II_space_dimension> &,