From 1056bfbf4619c678f732873c39facecd1c3f4f0d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 27 Sep 2000 15:26:17 +0000 Subject: [PATCH] Add a class that implements reordering of grids whenever this is necessary to fulfill the assumptions of the triangulation class. git-svn-id: https://svn.dealii.org/trunk@3351 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/grid/grid_reordering.h | 448 ++++++++++ .../deal.II/source/grid/grid_reordering.cc | 805 ++++++++++++++++++ 2 files changed, 1253 insertions(+) create mode 100644 deal.II/deal.II/include/grid/grid_reordering.h create mode 100644 deal.II/deal.II/source/grid/grid_reordering.cc diff --git a/deal.II/deal.II/include/grid/grid_reordering.h b/deal.II/deal.II/include/grid/grid_reordering.h new file mode 100644 index 0000000000..67f7c8cd8a --- /dev/null +++ b/deal.II/deal.II/include/grid/grid_reordering.h @@ -0,0 +1,448 @@ +//---------------------------- grid_reordering.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- grid_reordering.h --------------------------- +#ifndef __deal2__grid_reordering_h +#define __deal2__grid_reordering_h + + +#include +#include +#include + + +/** + * Class declaring some dimension dependent numbers which are needed + * for the grid reordering class. + * + * @author Wolfgang Bangerth, 2000 + */ +template +class GridReorderingInfo +{}; + + +/** + * Class declaring some dimension dependent numbers which are needed + * for the grid reordering class. This is the specialization for the + * 2d case. + * + * @author Wolfgang Bangerth, 2000 + */ +template <> +class GridReorderingInfo<2> +{ + public: + /** + * Number of possible valid + * orientations of a cell. They + * are the state in which it was + * delivered and three possible + * rotations in counter-clockwise + * sense, thus a total of four. + */ + static const unsigned int rotational_states_of_cells = 4; + + /** + * Number of possible + * orientations of a face in + * 2d. It is the face and the + * face with vertices exchanged, + * thus two. + */ + static const unsigned int rotational_states_of_faces = 2; +}; + + + +/** + * Class declaring some dimension dependent numbers which are needed + * for the grid reordering class. This is the specialization for the + * 2d case. + * + * @author Wolfgang Bangerth, 2000 + */ +template <> +class GridReorderingInfo<3> +{ + public: + /** + * ??? + */ + static const unsigned int rotational_states_of_cells = static_cast(-1); + static const unsigned int rotational_states_of_faces = 4; +}; + + + + + +/** + * @author Wolfgang Bangerth, 2000 + */ +template +class GridReordering : private GridReorderingInfo +{ + public: + /** + * This is the main function, + * doing what is announced in the + * general documentation of this + * class. + */ + static void reorder_cells (vector > &original_cells); + + private: + + /** + * Forward declarations of local + * classes. + */ + class Cell; + class Face; + class FaceData; + + /** + * Class that describes the + * properties of cells beyond + * what is provided by the data + * that is available from the + * calling functions of this + * class. In particular, several + * fields are available that + * describe connections of cells + * to faces and to + * neighbors. These fields are + * filled in a first pass before + * the actual reoordering starts, + * as they are needed for the + * latter purpose. + * + * Since this class is derived + * from the @ref{CellData} class, + * it also contains all the + * information available + * beforehand. + * + * @author Wolfgang Bangerth, 2000 + */ + struct Cell : public CellData + { + /** + * Value to be used if a + * neighbor does not exist, + * i.e. if the cell is at the + * boundary of the domain + * with a certain face. + */ + static const unsigned int invalid_neighbor = static_cast(-1); + + /** + * Pointers to the faces of + * this cell and their + * rotations. If the first + * index is zero, then the + * faces denote the faces in + * their standard direction + * with respect to the + * ordering in this cell. If + * it is nonzero, then they + * denote the faces that + * would be needed if the + * cell were rotate so often. + */ + typename map::iterator + faces[GridReorderingInfo::rotational_states_of_cells][GeometryInfo::faces_per_cell]; + + /** + * Cell indices of the + * neighbors of this cell in + * the global array of cells. + */ + unsigned int neighbors[GeometryInfo::faces_per_cell]; + + /** + * The index of this cell in + * the global array of cells. + */ + unsigned int cell_no; + + /** + * If we fail to insert this + * cell, then we have to + * track back. We could track + * back right to the previous + * cell, but we can do better + * than that, by tracking + * back to the cell indicated + * by this field. Which value + * it has is described in the + * documentation of the + * @ref{GridReordering} + * class. + */ + unsigned int track_back_to_cell; + + /** + * Default + * constructor. Invalidate + * all data. + */ + Cell (); + + /** + * Constructor that copies + * the data of an object of + * the base class and + * requires to be given the + * index of this cell in the + * global array. + */ + Cell (const CellData &cd, + const unsigned int cell_no); + + /** + * Count the existing neighbors + * of this cell. + */ + unsigned int count_neighbors () const; + + /** + * Insert the faces of the + * present cell into the map + * of all faces. This + * function inserts them in + * all orientations possible + * if the given cell is + * rotated. The function also + * takes care to fill in the + * @p{adjacent_cells} field + * of the inserted faces. + */ + void insert_faces (map &global_faces); + + /** + * Find out the neighbors of the + * given cell by looking at the + * @p{adjacent_cells} field of + * the faces of this cell. Store + * the neighbor indices in the + * present object. + */ + void fix_cell_neighbors (); + + /** + * Compute back to which cell we + * have to backtrack in case we + * can't insert this cell in any + * orientation into the already + * existing part of the + * triangulation. The method of + * how to determine the point to + * which we have to backtrack is + * described in the documentation + * of the @ref{GridReordering} + * class. + */ + void find_backtracking_point (); + + /** + * Find out whether the cell + * could be inserted into the + * already existing part of + * the triangulation with + * orientation given by the + * parameter, by checking + * that no face would be + * inserted twice in + * different orientations. + */ + bool check_consistency (const unsigned int rot) const; + + /** + * Tag the faces of this cell in + * the given orientation as used + * by this cell. + */ + void mark_faces_used (const unsigned int rot); + + /** + * Remove the use tags on the + * faces of this cell in the + * given orientation by this + * cell. Tags may remain if + * there is another cell that + * uses a given face. + */ + void mark_faces_unused (const unsigned int rot); + }; + + + /** + * Structure describing a face of + * a cell. This class is used as + * key in a map storing all faces + * possible in a triangulation, + * i.e. all faces between cells + * in all possible orientations. + * + * @author Wolfgang Bangerth, 2000 + */ + struct Face + { + /** + * Indices of the vertices of + * this face. + */ + unsigned int vertices[GeometryInfo::vertices_per_face]; + + /** + * Comparison operator. Use + * the vertex indices as + * primary, secondary, + * ... criteria for + * comparison. + */ + bool operator < (const Face &face) const; + }; + + + /** + * Class describing some data to + * be stored on each + * face. Objects of this type are + * used as values in a map + * containing all faces. + * + * @author Wolfgang Bangerth, 2000 + */ + struct FaceData + { + /** + * Value denoting + * non-existing adjacent + * cells of this face, + * i.e. when the face is at + * the boundary of the domain + * and has only one adjacent + * cell. + */ + static const unsigned int invalid_adjacent_cell = static_cast(-1); + + /** + * Pointers to the same face + * but in all other + * orientations. Storing + * these pointers makes it + * much easier to find out + * whether a given faces has + * already been used in + * another direction, thus + * forbidding the present + * face to be used. + */ + typename map::const_iterator + reverse_faces[GridReorderingInfo::rotational_states_of_faces-1]; + + /** + * Indices of the one or two + * adjacent cells of this + * face in the global array + * of cells. + */ + unsigned int adjacent_cells[2]; + + /** + * Number of cells presently + * using this face in the + * orientation represented by + * this object. May be zero, + * one, or two. + */ + unsigned int use_count; + + /** + * Default constructor. + */ + FaceData (); + }; + + + /** + * If we couldn't insert a cell + * into the already existing part + * of the mesh, then we need to + * track back a while. This + * function does so, given the + * array of cells, the stack of + * rotation states, and the + * position to which backtracking + * shall take place. + * + * In some cases, it is possible + * that from the place where we + * backtracked to, there is no + * more possibility to orient a + * cell. Then we will have to + * backtrack further until we + * come to a place where further + * work is possible; this + * recursive backtracking is also + * done by this function, + * although it is not implemented + * as recursive calls but rather + * as eliminated tail-recursion. + */ + static void track_back (vector &cells, + stack > &rotation_states, + unsigned int track_back_to_cell); + + /** + * This is the main function that + * does the main work. It is + * called by the + * @p{reorder_cells} function + * after all the preparations + * have been completed and + * operates on the @p{cells} + * array. After a way to reorder + * the cells has been found, the + * @p{original_cells} are reorder + * accordingly, where the + * @p{new_cell_numbers} array is + * needed to find the connection + * between original cells and + * presorted cells. + */ + static void find_reordering (vector &cells, + vector > &original_cells, + const vector &new_cell_numbers); + + /** + * Preorder the incoming cells by + * some kind of Cuthill-McKee + * algorithm. The reason for the + * need to do so is described in + * the general documentation. + * + * Return a vector in which for + * each old cell the new index is + * stored. + */ + static + vector + presort_cells (vector &cells, + map &faces); +}; + + + +#endif diff --git a/deal.II/deal.II/source/grid/grid_reordering.cc b/deal.II/deal.II/source/grid/grid_reordering.cc new file mode 100644 index 0000000000..a5e9c9dfea --- /dev/null +++ b/deal.II/deal.II/source/grid/grid_reordering.cc @@ -0,0 +1,805 @@ +//---------------------------- grid_reordering.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- grid_reordering.cc --------------------------- + + +#include +#include + +#include + + + + +template +GridReordering::Cell::Cell () : + cell_no (invalid_neighbor) +{ + for (unsigned int i=0; i::faces_per_cell; ++i) + neighbors[i] = invalid_neighbor; +}; + + + +template +GridReordering::Cell::Cell (const CellData &cd, + const unsigned int cell_no) : + CellData (cd), cell_no(cell_no) +{ + for (unsigned int i=0; i::faces_per_cell; ++i) + neighbors[i] = invalid_neighbor; +}; + + + +template +inline +unsigned int GridReordering::Cell::count_neighbors () const +{ + unsigned int n = 0; + for (unsigned int i=0; i::faces_per_cell; ++i) + if (neighbors[i] != invalid_neighbor) + ++n; + return n; +}; + + + +template +void +GridReordering::Cell::insert_faces (map &/*global_faces*/) +{ + Assert (false, ExcNotImplemented()); +}; + + + +template <> +void +GridReordering<2>::Cell::insert_faces (map &global_faces) +{ + const unsigned int dim = 2; + + // first compute index numbers for + // the faces in usual order as + // defined by the order of vertices + // in the cell object + Face new_faces[GeometryInfo::faces_per_cell] + = { { {vertices[0], vertices[1]} }, + { {vertices[1], vertices[2]} }, + { {vertices[3], vertices[2]} }, + { {vertices[0], vertices[3]} } }; + + // then insert them into the global + // list and store iterators to + // them. note that if the face + // already exists, then the stored + // data is not touched. + for (unsigned int face=0; face::faces_per_cell; ++face) + faces[0][face] = global_faces.insert (make_pair(new_faces[face], + FaceData())).first; + + + // then for each of the faces also + // insert the reverse form and + // store pointers to them. note + // that the rotational state in + // which all faces are reverted is + // `2' + for (unsigned int face=0; face::faces_per_cell; ++face) + { + swap (new_faces[face].vertices[0], + new_faces[face].vertices[1]); + faces[2][face] = global_faces.insert (make_pair(new_faces[face], + FaceData())).first; + }; + + // then finally fill in rotational + // states 1 and 3 of the cell. the + // faces of these states can be + // obtained from states 0 and 2 + faces[1][0] = faces[2][0]; + faces[1][1] = faces[0][1]; + faces[1][2] = faces[2][2]; + faces[1][3] = faces[0][3]; + + faces[3][0] = faces[0][0]; + faces[3][1] = faces[2][1]; + faces[3][2] = faces[0][2]; + faces[3][3] = faces[2][3]; + + + // finally fill the crosslink and + // other fields of the new + // entries. note that since + // rotational states 0 and 2 of the + // cell are exactly reverted, we + // only have to operate on the face + // pointers of these two states to + // reach all possible faces and + // permutations thereof + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (faces[0][face]->second.adjacent_cells[0] == + FaceData::invalid_adjacent_cell) + { + // face had not been + // inserted by previous + // cells, since first + // adjacent cell is still + // untouched. provide + // xlinks to rotated faces + faces[0][face]->second.reverse_faces[0] = faces[2][face]; + faces[2][face]->second.reverse_faces[0] = faces[0][face]; + + // and insert this cell as + // adjacent_cell of the faces + faces[0][face]->second.adjacent_cells[0] = cell_no; + faces[2][face]->second.adjacent_cells[0] = cell_no; + } + else + { + // face had already been + // inserted. make sure that + // it was in the same way: + Assert (faces[0][face]->second.reverse_faces[0] == faces[2][face], + ExcInternalError()); + Assert (faces[2][face]->second.reverse_faces[0] == faces[0][face], + ExcInternalError()); + + // now insert ourselves as + // second + // adjacent_cell. the + // respective slots must + // necessarily be empty + // still + Assert (faces[0][face]->second.adjacent_cells[1] == + FaceData::invalid_adjacent_cell, + ExcInternalError()); + Assert (faces[2][face]->second.adjacent_cells[1] == + FaceData::invalid_adjacent_cell, + ExcInternalError()); + faces[0][face]->second.adjacent_cells[1] = cell_no; + faces[2][face]->second.adjacent_cells[1] = cell_no; + }; + }; +}; + + + +template +void GridReordering::Cell::fix_cell_neighbors () +{ + for (unsigned int face=0; face::faces_per_cell; ++face) + { + // first assert that the + // neighborship info of all + // versions of the same face is + // identical + for (unsigned int rot=1; rotsecond.adjacent_cells[adjacent_cell] + == + faces[0][face]->second.adjacent_cells[adjacent_cell], + ExcInternalError()); + + + // then insert the neighbor + // behind this face as neighbor + // of the present cell. note + // that it is not relevant to + // which permutation of a face + // we refer. note that it might + // well be that some of the + // neighbor indices are + // FaceData::invalid_adjacent_cell + if (faces[0][face]->second.adjacent_cells[0] == cell_no) + neighbors[face] = faces[0][face]->second.adjacent_cells[1]; + else + neighbors[face] = faces[0][face]->second.adjacent_cells[0]; + }; +}; + + + +template +void GridReordering::Cell::find_backtracking_point () +{ + // we know what neighbors we have, + // we can determine the neighbor + // with the maximal cell_no that is + // smaller than that of the present + // cell. we need this information + // in the backtracking process and + // don't want to compute it every + // time again + track_back_to_cell = FaceData::invalid_adjacent_cell; + for (unsigned int face=0; face::faces_per_cell; ++face) + if ((neighbors[face] != FaceData::invalid_adjacent_cell) + && + (neighbors[face] < cell_no) + && + ((neighbors[face] > track_back_to_cell) + || + (track_back_to_cell == FaceData::invalid_adjacent_cell))) + track_back_to_cell = neighbors[face]; + + // if this cell had no neighbors + // with lower cell numbers, we + // still need to know what cell to + // track back to in case some + // higher cell than the present one + // failed to coexist with the + // existing part of the mesh + // irrespective of the rotation + // state of this present cell. we + // then simply track back to the + // cell before this one, lacking a + // better alternative. this does, + // of course, not hold for cell 0, + // from which we should never be + // forced to track back + if (track_back_to_cell == 0) + track_back_to_cell = 0; + else + if (track_back_to_cell == FaceData::invalid_adjacent_cell) + track_back_to_cell = cell_no-1; +}; + + + +template +inline +bool GridReordering::Cell::check_consistency (const unsigned int /*rot*/) const +{ + // might be that we can use the + // specialization for dim==2 after + // some modification, but haven't + // thought about that yet + Assert (false, ExcNotImplemented()); + return false; +}; + + + +template <> +inline +bool GridReordering<2>::Cell::check_consistency (const unsigned int rot) const +{ + const unsigned int dim = 2; + // make sure that for each face of + // the cell the permuted faces are + // not already in use, as that + // would make the cell disallowed + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + const FaceData &face = faces[rot][face_no]->second; + + for (unsigned int face_rot=0; face_rotsecond; + if (reverse_face.use_count != 0) + return false; + }; + }; + + // no conflicts found + return true; +}; + + +template +inline +void GridReordering::Cell::mark_faces_used (const unsigned int /*rot*/) +{ + // might be that we can use the + // specialization for dim==2 after + // some modification, but haven't + // thought about that yet. in + // particular, I don't know whether + // we have to treat edges in 3d + // specially + Assert (false, ExcNotImplemented()); +}; + + + +template +inline +void GridReordering::Cell::mark_faces_unused (const unsigned int /*rot*/) +{ + // might be that we can use the + // specialization for dim==2 after + // some modification, but haven't + // thought about that yet. in + // particular, I don't know whether + // we have to treat edges in 3d + // specially + Assert (false, ExcNotImplemented()); +}; + + + +template <> +inline +void GridReordering<2>::Cell::mark_faces_used (const unsigned int rot) +{ + const unsigned int dim=2; + for (unsigned int face=0; face::faces_per_cell; ++face) + { + Assert (faces[rot][face]->second.use_count < 2, + ExcInternalError()); + ++faces[rot][face]->second.use_count; + }; +}; + + + +template <> +inline +void GridReordering<2>::Cell::mark_faces_unused (const unsigned int rot) +{ + const unsigned int dim=2; + for (unsigned int face=0; face::faces_per_cell; ++face) + { + Assert (faces[rot][face]->second.use_count > 0, + ExcInternalError()); + --faces[rot][face]->second.use_count; + }; +}; + + + +template +bool GridReordering::Face::operator < (const Face &face) const +{ + for (unsigned int v=0; v::vertices_per_face; ++v) + { + // if vertex index is smaller, + // then comparison is true + if (vertices[v] < face.vertices[v]) + return true; + else + // if vertex index is greater, + // then comparison is false + if (vertices[v] > face.vertices[v]) + return false; + // if indices are equal, then test + // next index + }; + + // if all indices are equal: + return false; +}; + + + +template +GridReordering::FaceData::FaceData () : + use_count (0) +{ + adjacent_cells[0] = adjacent_cells[1] = invalid_adjacent_cell; +}; + + + + + + +template +inline +void GridReordering::track_back (vector > &cells, + stack > &rotation_states, + unsigned int track_back_to_cell) +{ + top_of_function: + + Assert (track_back_to_cell > 0, ExcInternalError()); + + unsigned int last_rotation_state; + for (unsigned int cell_no=rotation_states.size()-1; cell_no>=track_back_to_cell; --cell_no) + { + // store rotation state of + // topmost cell, as we will + // have to advance that by one + last_rotation_state = rotation_states.top(); + + // first mark faces of that + // cell as no more used + cells[cell_no].mark_faces_unused (last_rotation_state); + + // then pop state from + // stack + rotation_states.pop(); + }; + + // now we will have to find out + // whether we can try the last cell + // we have popped from the stack in + // another rotation state, or will + // have to backtrack further: + if (last_rotation_state < rotational_states_of_cells-1) + { + // possible. push that state to + // the stack and leave + rotation_states.push (last_rotation_state+1); + return; + } + else + { + // last cell can't be rotated + // further. go on with + // backtracking + const typename vector::iterator + try_cell = cells.begin() + rotation_states.size(); +// cout << "Further backtracking from " << rotation_states.size(); +// cout << ". Neighbors are "; +// for (unsigned int i=0; i::faces_per_cell; ++i) +// cout << (int)try_cell->neighbors[i] << ' '; +// cout << "Will track back to cell " << try_cell->track_back_to_cell << endl; + + track_back_to_cell = try_cell->track_back_to_cell; + Assert (track_back_to_cell > 0, ExcInternalError()); + + // track further back. this + // could be done by recursive + // calls of this function, + // which in this case would + // represent a tail-recursion + // as there is nothing more to + // be done after calling the + // function recursively, but we + // prefer to write down the + // tail-recursion by hand using + // a goto, since the compiler + // seems to have problems to + // rewrite the tail recursion + // as a goto. + goto top_of_function; + }; +}; + + + +template +void GridReordering::find_reordering (vector > &cells, + vector > &original_cells, + const vector &new_cell_numbers) +{ + const unsigned int n_cells = cells.size(); + + // stack of value indicating that + // the nth cell needs to be rotated + // so-and-so often, where n is the + // position on the stack + stack > rotation_states; + + // for the first cell, the + // rotational state can never be + // important, since we can rotate + // all other cells + // accordingly. therefore preset + // the rotation state of the first + // cell + rotation_states.push (0); + cells[0].mark_faces_used (rotation_states.top()); + + while (true) + { + // if all cells have a coherent + // orientation, then we can + // exit the main loop + if (rotation_states.size() == n_cells) + break; + + // try to push back another + // cell in orientation zero + rotation_states.push (0); + + // check whether the present + // cell in the present + // orientation is valid + check_topmost_cell: + + static unsigned int max_size = 0; + if (rotation_states.size() > max_size) + { + if (max_size % 10 == 0) + cout << "New max size " << rotation_states.size() << endl; + max_size = rotation_states.size(); + }; + + const typename vector::iterator + try_cell = cells.begin() + rotation_states.size()-1; + if (try_cell->check_consistency (rotation_states.top())) + { + // yes, works, we found a + // way of how to add the + // present cell to the + // existing cells without + // violating any ordering + // constraints. now mark + // the respective faces as + // used and go on with the + // next cell + try_cell->mark_faces_used (rotation_states.top()); + +// cout << "Added cell " << try_cell->cell_no +// << " in rotation " << rotation_states.top() << endl; + // go on with next cell + continue; + } + else + { + // no, doesn't work. see if + // we can rotate the top + // cell so that it works + if (rotation_states.top()+1 < rotational_states_of_cells) + { + // yes, can be + // done. then do so and + // check again + ++rotation_states.top(); + goto check_topmost_cell; + } + else + { + // no, no more + // orientation of the + // top cell possible, + // we have to backtrack + // some way +// cout << "Failure with cell " << rotation_states.size()-1; +// cout << ". Neighbors are "; +// for (unsigned int i=0; i::faces_per_cell; ++i) +// cout << (int)try_cell->neighbors[i] << ' '; +// cout << "Will track back to cell " << try_cell->track_back_to_cell << endl; + + // first pop rotational + // state of top cell, + // since for that no + // faces have been + // marked as used yet + rotation_states.pop(); + // then track back + track_back (cells, rotation_states, try_cell->track_back_to_cell); + // and go on by + // checking the now + // topmost cell + goto check_topmost_cell; + }; + }; + }; + + + + // rotate the cells according to + // the results we have found. since + // we operate on a stack, we do the + // rotations from the back of the + // array to the front + while (rotation_states.size() != 0) + { + const unsigned int + new_cell_number = rotation_states.size()-1; + const unsigned int + old_cell_number = find (new_cell_numbers.begin(), + new_cell_numbers.end(), + new_cell_number) - new_cell_numbers.begin(); + Assert (old_cell_number < cells.size(), ExcInternalError()); + + original_cells[old_cell_number].rotate (rotation_states.top()); + rotation_states.pop (); + }; +}; + + + +template +vector +GridReordering::presort_cells (vector > &cells, + map &faces) +{ + // first find the cell with the + // least neighbors + unsigned int min_neighbors = cells[0].count_neighbors(); + unsigned int cell_with_min_neighbors = 0; + for (unsigned int i=1; i cells[i].count_neighbors()) + { + min_neighbors = cells[i].count_neighbors(); + cell_with_min_neighbors = i; + if (min_neighbors == 1) + // better is not possible + break; + }; + + // have an array into which we + // insert the new cells numbers of + // each cell + const unsigned int invalid_cell_number = static_cast(-1); + vector new_cell_numbers (cells.size(), invalid_cell_number); + + // and have an array of the next + // cells to be numbered (old numbers) + vector next_round_cells (1, cell_with_min_neighbors); + + unsigned int next_free_new_number = 0; + + // while there are still cells to + // be renumbered: + while (next_round_cells.size() != 0) + { + for (unsigned int i=0; i new_next_round_cells; + for (unsigned int i=0; i::faces_per_cell; ++n) + if (cells[next_round_cells[i]].neighbors[n] != Cell::invalid_neighbor) + if (new_cell_numbers[cells[next_round_cells[i]].neighbors[n]] + == invalid_cell_number) + new_next_round_cells.push_back (cells[next_round_cells[i]].neighbors[n]); + + + + // if no more cells have been + // found, then we must have + // renumbered all cells already + if (new_next_round_cells.size() == 0) + Assert (next_free_new_number == cells.size(), ExcInternalError()); + + // eliminate duplicates from + // the new_next_round_cells + // array. note that a cell + // which is entered into this + // array might have been + // entered more than once since + // it might be a neighbor of + // more than one cell of the + // present round + // + // in order to eliminate + // duplicates, we first sort + // tha array and then copy over + // only unique elements to the + // next_round_cells array, + // which is needed for the next + // loop iteration anyway + sort (new_next_round_cells.begin(), new_next_round_cells.end()); + next_round_cells.clear (); + unique_copy (new_next_round_cells.begin(), new_next_round_cells.end(), + back_inserter(next_round_cells)); + }; + Assert (find (new_cell_numbers.begin(), new_cell_numbers.end(), invalid_cell_number) + == + new_cell_numbers.end(), + ExcInternalError()); + + // now that we know in which order + // to sort the cells, do so: + vector > new_cells (cells.size()); + for (unsigned int i=0; i::faces_per_cell; ++n) + cells[c].neighbors[n] = new_cell_numbers[cells[c].neighbors[n]]; + }; + + for (typename map::iterator i=faces.begin(); i!=faces.end(); ++i) + for (unsigned int k=0; k<2; ++k) + if (i->second.adjacent_cells[k] != FaceData::invalid_adjacent_cell) + i->second.adjacent_cells[k] = new_cell_numbers[i->second.adjacent_cells[k]]; + + return new_cell_numbers; +}; + + + +template +void GridReordering::reorder_cells (vector > &original_cells) +{ + // the present function might + // actually work in 3d, but the + // ones it calls probably not. will + // have to think about what needs + // to be changed. however, rather + // than killing the program when + // calling functions that may not + // work, kill it here if the + // dimension is not appropriate + Assert (dim==2, ExcNotImplemented()); + + // we need more information than + // provided by the input parameter, + // in particular we need + // neighborship relations between + // cells. therefore copy over the + // old cells to another class that + // provides space to these + // informations + vector > cells; + cells.reserve (original_cells.size()); + for (unsigned int i=0; i(original_cells[i], i)); + + // first generate all the faces + // possible, i.e. in each possible + // direction and rotational state + map faces; + for (unsigned int cell_no=0; cell_no + new_cell_numbers = presort_cells (cells, faces); + + // finally do some preliminary work + // to make backtracking simpler + // later + for (unsigned int cell_no=0; cell_no +void GridReordering<1>::reorder_cells (vector > &) +{ + // there should not be much to do + // in 1d... +}; + + + + +// explicit instantiations. only require the main function, it should +// then claim whatever templates it needs +template +void +GridReordering:: +reorder_cells (vector > &); + -- 2.39.5