From: wolf Date: Fri, 21 Feb 2003 20:36:30 +0000 (+0000) Subject: Merge Michaels new reorientation scheme from the respective branch. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1c82d004a0647e9e7e6f2f474dccce10bc862b40;p=dealii-svn.git Merge Michaels new reorientation scheme from the respective branch. git-svn-id: https://svn.dealii.org/trunk@7219 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/grid_reordering.h b/deal.II/deal.II/include/grid/grid_reordering.h index 1201a228cb..df04870f54 100644 --- a/deal.II/deal.II/include/grid/grid_reordering.h +++ b/deal.II/deal.II/include/grid/grid_reordering.h @@ -35,39 +35,6 @@ namespace internal }; -/** - * 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 @@ -110,7 +77,341 @@ namespace internal } + + + + +namespace internal +{ +/** + * Implement the algorithm described in the documentation of the + * GridReordering<2> class. + * + * @author Michael Anderson, 2003 + */ + namespace GridReordering2d + { +/** + * Defines a variety of variables related to the connectivity of a + * simple quad element. This includes the nodes on each edge, which + * edges come into each node and what the default deal.II directions + * are for the quad. + * + * @begin{verbatim} + * s2 + * + * +-->--+ + * |3 2| + * s3 ^ ^ s1 + * |0 1| + * +-->--+ + * + * s0 + * @end{verbatim} + * + * @author Michael Anderson, 2003 + */ + class ConnectGlobals + { + public: + /** + * The nodes on each edge in + * anti-clockwise order + * { {0,1},{1,2},{2,3},{3,0} } + */ + static const int EdgeToNode[4][2]; + + /** + * The edges comin into each + * node, in anti-clockwise + * order + * { {3,0},{0,1},{1,2},{2,3} } + */ + static const int NodeToEdge[4][2]; + + /** + * The nodes on each edge in + * "default direction order". + * {{0,1},{1,2},{3,2},{0,3}} + */ + static const int DefaultOrientation[4][2]; + }; + + +/** + * An enriched quad with information about how the mesh fits together + * so that we can move around the mesh efficiently. + * + * @author Michael Anderson, 2003 + */ + class MQuad + { + public: + /** + * v0 - v3 are indexes of the vertices of the quad, + * s0 - s3 are indexes for the sides of the quad + */ + MQuad (const unsigned int v0, + const unsigned int v1, + const unsigned int v2, + const unsigned int v3, + const unsigned int s0, + const unsigned int s1, + const unsigned int s2, + const unsigned int s3, + const CellData<2> &cd); + + /** + * Stores the vertex numbers + */ + unsigned int v[4]; + /** + * Stores the side numbers + */ + unsigned int side[4]; + + /** + * Copy of the @p{CellData} object + * from which we construct the + * data of this object. + */ + CellData<2> original_cell_data; + + /** + * Makes an MQuad from the + * given CellData and MSide + * list. Is derived from + * binary_function to be + * usable with STL + * containers. + * + * Also assumes that the + * edges listed present in + * the CellData are already + * present in the elist + * vector. + */ + struct MakeQuad; + }; + +/** + * The enriched side class containing connectivity information. + * Orientation is from v0 to v1; Initially this should have v0 > &quads); + private: + + /** + * Sets up the internal data + * structures so that the we can + * do side hopping and face + * switching efficiently. This + * means we need a whole bunch of + * connectivity information + */ + void build_graph (const std::vector > &inquads); + + /** + * Orient the internal data + * into deal.II format The + * orientation algorith is as + * follows + * + * 1) Find an unoriented quad (A) + * + * 2) Orient an un_oriented side (s) of (A) + * + * 3) side hop on (s) of (A) to get (B) + * + * 4) if opposite side to (s) + * of (B) is unoriented + * orient it + * + * 5) repeat 3) and 4) until + * side-hoppong fails (we've + * reached a boundary) or (s) + * has already been oriented + * (we've closed a loop or + * unoriented sides). + * + * 6) Repeat 2), 3) ,4) and + * 5) on other unoriented + * sides of (A) + * + * 7) Choose a new unoriented + * A. + */ + void orient(); + + /** + * Get the (now correctly + * oriented if we've called + * orient) quads. + */ + void get_quads(std::vector > &outquads) const; + + /** + * Orient_side(qnum,lsn) + * orients the local side lsn + * of the quad qnum in the + * triangulation. If the side + * opposite lsn is oriented + * then lsn is oriented to + * match it. Otherwise it is + * oriented in the "default" + * direction for the quad. + */ + void orient_side (const unsigned int quadnum, + const unsigned int localsidenum); + + /** + * Returns true if all sides + * of the quad quadnum are + * oriented. + */ + bool is_fully_oriented_quad (const unsigned int quadnum) const; + + /** + * Returns true if the side lsn + * of the quad quadnum is + * oriented. + */ + bool is_oriented_side (const unsigned int quadnum, + const unsigned int lsn) const; + + /** + * Returns true is the side is + * oriented in the "default" + * direction + */ + bool is_side_default_oriented (const unsigned int qnum, + const unsigned int lsn) const; + + /** + * Increases UnOrQLoc from + * it's original value to the + * next quad with an + * unoriented side. Returns + * true if there was another + * unoriented quad. + */ + bool get_unoriented_quad (unsigned int &UnOrQLoc) const; + + /** + * Sets sidenum to the local + * sidenumber of an + * unoriented side of the + * quad quadnum. Returns true + * if such a side exists. + */ + bool get_unoriented_side (const unsigned int quadnum, + unsigned int &sidenum) const; + + /** + * side_hop(&qnum, &lsn) has + * qnum being the quadnumber + * of a quad in the + * triangulation, and a local + * side number. side_hop then + * sets qnum to the + * quadnumber across the + * other side of the side, + * and sets lsn so that + * quads[qnum].sides[lsn] is + * the same before and after + * the call. if there is no + * other quad on the other + * side of the current quad, + * then side_hop returns + * false. + */ + bool side_hop (unsigned int &qnum, + unsigned int &lsn) const; + /** + * Sets lsn so that it points + * to the opposite side of + * the current quad (qnum) + * that it was originally + * pointing to. + */ + bool switch_faces (unsigned int &qnum, + unsigned int &lsn) const; + + /** + * A list of enriched + * sides/edges of the mesh. + */ + std::vector sides; + /** + * A list of enriched quads + * in the mesh. + */ + std::vector mquads; + }; + } // namespace GridReordering2d +} // namespace internal + + /** @@ -171,7 +472,7 @@ namespace internal * | | | * o---o---o * @end{verbatim} - * (The reader is aked to try to find a conforming choice of line + * (The reader is asked to try to find a conforming choice of line * directions; it will soon be obvious that there can't exists such a * thing, even if we allow that there might be cells with clockwise * and counterclockwise orientation of the lines at the same time.) @@ -191,13 +492,29 @@ namespace internal * The purpose of this class is now to find an ordering for a given * set of cells such that the generated triangulation satisfies all * the requirements stated above. To this end, we will first show some - * examples why this is a difficult problem, and then develop an - * algorithm that finds such a reordering. Note that the algorithm + * examples why this is a difficult problem, and then develop + * algorithms that finds such a reordering. Note that the algorithm * operates on a set of @ref{CellData} objects that are used to * describe a mesh to the triangulation class. These objects are, for * example, generated by the @ref{GridIn} class, when reading in grids * from input files. * + * As a last question for this first section: is it guaranteed that + * such orientations of faces always exist for a given subdivision of + * a domain into cells? The linear complexity algorithm described + * below for 2d also proves that the answer is yes for 2d. For 3d, the + * answer is no (which also underlines that using such orientations + * might be an -- unfortunately uncurable -- misfeature of deal.II). A + * simple counter-example in 3d illustrates this: take a string of 3d + * cells and bend it together to a torus. Since opposing lines in a + * cell need to have the same direction, there is a simple ordering + * for them, for example all lines radially outward, tangentially + * clockwise, and axially upward. However, if before joining the two + * ends of the string of cells, the string is twisted by 180 degrees, + * then no such orientation is possible any more, as can easily be + * checked. In effect, some meshes cannot be used in deal.II, + * unfortunately. + * * * @sect3{Examples of problems} * @@ -305,8 +622,8 @@ namespace internal * requirements of deal.II triangulations are met. * * These two examples demonstrate that if we have added a certain - * number of cells in some oeirntation of faces and can't add the next - * one without introducingfaces that had already been added in another + * number of cells in some orientation of faces and can't add the next + * one without introducing faces that had already been added in another * direction, then it might not be sufficient to only rotate cells in * the neighborhood of the the cell that we failed to add. It might be * necessary to go back a long way and rotate cells that have been @@ -325,8 +642,8 @@ namespace internal * rotated cell 1, then we would have to rotate the cells 1 through * N-1 as well). * - * The only solution to this problem seems to be the following: if - * cell N can't be added, the try to rotate cell N-1. If we can't + * A brute force approach to this problem is the following: if + * cell N can't be added, then try to rotate cell N-1. If we can't * rotate cell N-1 any more, then try to rotate cell N-2 and try to * add cell N with all orientations of cell N-1. And so * on. Algorithmically, we can visualize this by a tree structure, @@ -361,7 +678,7 @@ namespace internal * that has already been added, then there are already only two * possible orientations left, so the total number of checks we have * to make until we find a valid way is significantly smaller than - * @p{4**N}. However, an algorithm is still exponential in time and + * @p{4**N}. However, the algorithm is still exponential in time and * linear in memory (we only have to store the information for the * present path in form of a stack of orientations of cells that have * already been added). @@ -371,26 +688,162 @@ namespace internal * very first cells there to find a way to add all cells in a * consistent fashion. * - * This discouraging situation is geatly improved by the fact that we - * can find an algorithm that in practice is usually only roughly - * linear in time and memory. We will describe this algorithm in the - * following. + * This discouraging situation is greatly improved by the fact that we + * have an alternative algorithm for 2d that is always linear in + * runtime (discovered and implemented by Michael Anderson of TICAM, + * University of Texas, in 2003), and that for 3d we can find an + * algorithm that in practice is usually only roughly linear in time + * and memory. We will describe these algorithms in the following. + * + * + * @sect3{The 2d linear complexity algorithm} + * + * The algorithm uses the fact that opposite faces of a cell need to + * have the same orientation. So you start with one arbitrary line, + * choose an orientation. Then the orientation of the opposite face is + * already fixed. Then go to the two cells across the two faces we + * have fixed: for them, one face is fixed, so we can also fix the + * opposite face. Go on with doing so. Eventually, we have done this + * for a string of cells. Then take one of the non-fixed faces of a + * cell which has already two fixed faces and do all this again. + * + * In more detail, the algorithm is best illustrated using an + * example. We consider the mesh below: + * @begin{verbatim} + * 9------10-------11 + * | | /| + * | | / | + * | | / | + * 6------7-----8 | + * | | | | + * | | | | + * | | | | + * 3------4-----5 | + * | | \ | + * | | \ | + * | | \| + * 0------1---------2 + * @end{verbatim} + * First a cell is chosen ( (0,1,4,3) in this case). A single side of the cell + * is oriented arbitrarily (3->4). This choice of orientation is then propogated + * through the mesh, across sides and elements. (0->1), (6->7) and (9->10). + * The involves edge-hopping and face hopping, giving a path through the mesh + * shown in dots. + * @begin{verbatim} + * 9-->--10-------11 + * | . | /| + * | . | / | + * | . | / | + * 6-->--7-----8 | + * | . | | | + * | . | | | + * | . | | | + * 3-->--4-----5 | + * | . | \ | + * | X | \ | + * | . | \| + * 0-->--1---------2 + * @end{verbatim} + * This is then repeated for the other sides of the chosen element, orienting + * more sides of the mesh. + * @begin{verbatim} + * 9-->--10-------11 + * | | /| + * v.....v.......V | + * | | /. | + * 6-->--7-----8 . | + * | | | . | + * | | | . | + * | | | . | + * 3-->--4-----5 . | + * | | \. | + * ^..X..^.......^ | + * | | \| + * 0-->--1---------2 + * @end{verbatim} + * Once an element has been completely oriented it need not be considered + * further. These elements are filled with o's in the diagrams. We then move + * to the next element. + * @begin{verbatim} + * 9-->--10->-----11 + * | ooo | . /| + * v ooo v . V | + * | ooo | . / | + * 6-->--7-->--8 | + * | | . | | + * | | . | | + * | | . | | + * 3-->--4-->--5 | + * | ooo | . \ | + * ^ ooo ^ X ^ | + * | ooo | . \| + * 0-->--1-->------2 + * @end{verbatim} + * Repeating this gives + * @begin{verbatim} + * 9-->--10->-----11 + * | ooo | oooooo /| + * v ooo v ooooo V | + * | ooo | oooo / | + * 6-->--7-->--8 | + * | | | | + * ^.....^..X..^...^ + * | | | | + * 3-->--4-->--5 | + * | ooo | oooo \ | + * ^ ooo ^ ooooo ^ | + * | ooo | oooooo \| + * 0-->--1-->------2 + * @end{verbatim} + * and the final oriented mesh is + * @begin{verbatim} + * 9-->--10->-----11 + * | | /| + * v v V | + * | | / | + * 6-->--7-->--8 | + * | | | | + * ^ ^ ^ ^ + * | | | | + * 3-->--4-->--5 | + * | | \ | + * ^ ^ ^ | + * | | \| + * 0-->--1-->-------2 + * @end{verbatim} + * It is obvious that this algorithm has linear run-time, since it + * only ever touches each face exactly once. + * + * The algorithm just described is implemented in a specialization of + * this class for the 2d case. Note that in principle, it should be + * possible to extend this algorithm to 3d as well, using sheets + * instead of strings of cells to work on. If a grid is reorientable, + * then such an algorithm should be able to do so in linear time; if + * it is not orientable, then it should abort in linear time as + * well. However, this has not yet been implemented. Rather, we use + * the backtracking and branch pruning algorithm for 3d described + * below; this algorithm predates the 2d linear complexity algorithm + * and was initially also used in 2d. + * + * + * @sect3{The 3d branch reshuffling and pruning algorithm} * * The first observation is that although there are counterexamples, - * problems are usually local. For example, in the second example, if - * we had numbered the cells in a way that neighboring cells have - * similar cell numbers, then the amount pf backtracking needed is - * greatly reduced. Therefore, in the implementation of the algorithm, - * the first step is to renumber the cells in a Cuthill-McKee fashion: - * start with the cell with the least number of neighbors and assign - * to it the cell number zero. Then find all neighbors of this cell - * and assign to them consecutive further numbers. Then find their - * neighbors that have not yet been numbered and assign to them - * numbers, and so on. Graphically, this represents finding zones of - * cells consecutively further away from the initial cells and number - * them in this front-marching way. This already greatly improves - * locality of problems and consequently reduced the necessary amount - * of backtracking. + * problems are usually local. For example, in the second example + * mentioned above, if we had numbered the cells in a way that + * neighboring cells have similar cell numbers, then the amount of + * backtracking needed is greatly reduced. Therefore, in the + * implementation of the algorithm, the first step is to renumber the + * cells in a Cuthill-McKee fashion: start with the cell with the + * least number of neighbors and assign to it the cell number + * zero. Then find all neighbors of this cell and assign to them + * consecutive further numbers. Then find their neighbors that have + * not yet been numbered and assign to them numbers, and so + * on. Graphically, this represents finding zones of cells + * consecutively further away from the initial cells and number them + * in this front-marching way. This already greatly improves locality + * of problems and consequently reduced the necessary amount of + * backtracking. * * The second point is that we can use some methods to prune the tree, * which usually lead to a valid orientation of all cells very @@ -408,13 +861,15 @@ namespace internal * neighbor of N with the largest cell index and which has already * been added. * - * Unfortunately, this method can fail to yield a valid path through the - * tree if not applied with care. Consider the following situation, - * initially extracted from a mesh of 950 cells generated + * Unfortunately, this method can fail to yield a valid path through + * the tree if not applied with care. Consider the following + * situation, initially extracted from a mesh of 950 cells generated * automatically by the program BAMG (this program usually generates * meshes that are quite badly balanced, often have many -- sometimes * 10 or more -- neighbors of one vertex, and exposed several problems - * in the initial algorithm): + * in the initial algorithm; note also that the example is in 2d where + * we now have the much better algorithm described above, but the same + * observations also apply to 3d): * @begin{verbatim} * 13----------14----15 * | \ | | @@ -518,7 +973,7 @@ namespace internal * the triangulation from this data using the * @ref{Triangulation}@p{::create_triangulation} function. * - * @author Wolfgang Bangerth, 2000 + * @author Wolfgang Bangerth, 2000, Michael Anderson 2003 */ template class GridReordering @@ -946,14 +1401,51 @@ class GridReordering -/* -------------- declaration of explicit specializations ------------- */ +/** + * This is the specialization of the general template for 1d. In 1d, + * there is actually nothing to be done. + * + * @author Wolfgang Bangerth, 2000 + */ +template <> +class GridReordering<1> +{ + public: + /** + * Do nothing, since in 1d no + * reordering is necessary. + */ + static void reorder_cells (const std::vector > &); +}; + + +/** + * This specialization of the general template implements the + * 2d-algorithm described in the documentation of the general + * template. + * + * @author Michael Anderson, 2003 + */ template <> -void GridReordering<2>::Cell::insert_faces (std::map &global_faces); +class GridReordering<2> +{ + public: + /** + * This is the main function, + * doing what is announced in + * the general documentation of + * this class. + */ + static void reorder_cells (std::vector > &original_cells); +}; + + + +/* -------------- declaration of explicit specializations ------------- */ + template <> void GridReordering<3>::Cell::insert_faces (std::map &global_faces); -template <> -void GridReordering<1>::reorder_cells (std::vector > &); #endif diff --git a/deal.II/deal.II/source/grid/grid_reordering.cc b/deal.II/deal.II/source/grid/grid_reordering.cc index bbb6ff1eaf..64713dfd69 100644 --- a/deal.II/deal.II/source/grid/grid_reordering.cc +++ b/deal.II/deal.II/source/grid/grid_reordering.cc @@ -11,7 +11,6 @@ // //---------------------------- grid_reordering.cc --------------------------- - #include #include @@ -22,11 +21,6 @@ namespace internal { // static variables -#if deal_II_dimension == 2 - const unsigned int GridReorderingInfo<2>::rotational_states_of_cells; - const unsigned int GridReorderingInfo<2>::rotational_states_of_faces; -#endif - #if deal_II_dimension == 3 const unsigned int GridReorderingInfo<3>::rotational_states_of_cells; const unsigned int GridReorderingInfo<3>::rotational_states_of_faces; @@ -42,6 +36,19 @@ const unsigned int GridReordering::FaceData::invalid_adjacent_cell; +#if deal_II_dimension == 1 + +void GridReordering<1>::reorder_cells (const std::vector > &) +{ + // there should not be much to do + // in 1d... +} + +#endif + + + +#if deal_II_dimension == 3 template GridReordering::Cell::Cell () : @@ -99,122 +106,6 @@ GridReordering::Cell::insert_faces (std::map &/*global_faces } -#if deal_II_dimension == 2 - -template <> -void -GridReordering<2>::Cell::insert_faces (std::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] - = { { { this->vertices[0], this->vertices[1] } }, - { { this->vertices[1], this->vertices[2] } }, - { { this->vertices[3], this->vertices[2] } }, - { { this->vertices[0], this->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 (std::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) - { - std::swap (new_faces[face].vertices[0], - new_faces[face].vertices[1]); - faces[2][face] = global_faces.insert (std::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; - }; - }; -} - -#endif - -#if deal_II_dimension == 3 template <> void @@ -543,7 +434,6 @@ GridReordering<3>::Cell::insert_faces (std::map &global_faces) [cell_orientation_faces[rot][face].first]; } -#endif template @@ -1293,25 +1183,466 @@ void GridReordering::reorder_cells (std::vector > &original_c Assert (i->second.use_count == 0, ExcInternalError()); } +#endif // deal_II_dimension == 3 +#if deal_II_dimension == 2 -#if deal_II_dimension == 1 +namespace internal +{ + namespace GridReordering2d + { +// -- Definition Of conectivity information -- + const int ConnectGlobals::EdgeToNode[4][2]= + { {0,1},{1,2},{2,3},{3,0} }; -template <> -void GridReordering<1>::reorder_cells (std::vector > &) + const int ConnectGlobals::NodeToEdge[4][2]= + { {3,0},{0,1},{1,2},{2,3} }; + + const int ConnectGlobals::DefaultOrientation[4][2]= + {{0,1},{1,2},{3,2},{0,3}}; + + + + + struct MSide::SideRectify : public std::unary_function + { + void operator() (MSide &s) const + { + if (s.v0>s.v1) + std::swap (s.v0, s.v1); + } + }; + + + struct MSide::SideSortLess : public std::binary_function + { + bool operator()(const MSide &s1, const MSide &s2) const + { + int s1vmin,s1vmax; + int s2vmin,s2vmax; + if (s1.v0s2vmin) + return false; + return s1vmax object. + */ + MSide quadside(const CellData<2> &q, unsigned int i) + { + Assert (i<4, ExcInternalError()); + return MSide(q.vertices[ConnectGlobals::EdgeToNode[i][0]], + q.vertices[ConnectGlobals::EdgeToNode[i][1]]); + } + + +/** + * Wrapper class for the quadside() function + */ + struct QuadSide: public std::binary_function,int,MSide> + { + MSide operator()(const CellData<2>& q, int i) const + { + return quadside(q,i); + } + }; + + + + MQuad::MQuad (const unsigned int v0, + const unsigned int v1, + const unsigned int v2, + const unsigned int v3, + const unsigned int s0, + const unsigned int s1, + const unsigned int s2, + const unsigned int s3, + const CellData<2> &cd) + : + original_cell_data (cd) + { + v[0]=v0; + v[1]=v1; + v[2]=v2; + v[3]=v3; + side[0]=s0; + side[1]=s1; + side[2]=s2; + side[3]=s3; + } + + + MSide::MSide (const unsigned int initv0, + const unsigned int initv1) + : + v0(initv0), v1(initv1), + Q0(static_cast(-1)),Q1(static_cast(-1)), + lsn0(static_cast(-1)),lsn1(static_cast(-1)), + Oriented(false) + {}; + + + + bool + MSide::operator== (const MSide& s2) const + { + if ((v0==s2.v0)&&(v1==s2.v1)) {return true;} + if ((v0==s2.v1)&&(v1==s2.v0)) {return true;} + return false; + } + + + struct MQuad::MakeQuad : public std::binary_function, + std::vector, + MQuad> + { + MQuad operator()(const CellData<2> &q, + const std::vector &elist) const + { + //Assumes that the sides + //are in the vector.. Bad + //things will happen if + //they are not! + return MQuad(q.vertices[0],q.vertices[1], q.vertices[2], q.vertices[3], + std::distance(elist.begin(), + std::lower_bound(elist.begin(), elist.end(), + quadside(q,0), + MSide::SideSortLess() )), + std::distance(elist.begin(), + std::lower_bound(elist.begin(), elist.end(), + quadside(q,1), + MSide::SideSortLess() )), + std::distance(elist.begin(), + std::lower_bound(elist.begin(), elist.end(), + quadside(q,2), + MSide::SideSortLess() )), + std::distance(elist.begin(), + std::lower_bound(elist.begin(), elist.end(), + quadside(q,3), + MSide::SideSortLess() )), + q); + } + + }; + + + + void + GridReordering::reorient(std::vector > &quads) + { + build_graph(quads); + orient(); + get_quads(quads); + } + + + void + GridReordering::build_graph (const std::vector > &inquads) + { + //Reserve some space + sides.reserve(4*inquads.size()); + mquads.reserve(inquads.size()); + + //Insert all the sides into the side vector + for (int i=0;i<4;++i) + { + std::transform(inquads.begin(),inquads.end(), + std::back_inserter(sides), std::bind2nd(QuadSide(),i)); + } + + //Change each edge so that v0(sides).swap(sides); + + //Assigns the correct sides to + //each quads + transform(inquads.begin(),inquads.end(), back_inserter(mquads), + std::bind2nd(MQuad::MakeQuad(),sides) ); + + // Assign the quads to their sides also. + int qctr=0; + for(std::vector::iterator it=mquads.begin(); it!=mquads.end(); ++it) + { + for(int i=0;i<4;++i) + { + MSide &ss =sides[(*it).side[i]]; + if(ss.Q0==static_cast(-1)) + { + ss.Q0=qctr; + ss.lsn0=i; + } + else if (ss.Q1==static_cast(-1)) + { + ss.Q1=qctr; + ss.lsn1=i; + } + else + { + exit(0); + } + } + qctr++; + } + } + + + void GridReordering::orient() + { + // do what the comment in the + // class declaration says + unsigned int qnum=0; + while(get_unoriented_quad(qnum)) + { + unsigned int lsn=0; + while(get_unoriented_side(qnum,lsn)) + { + orient_side(qnum,lsn); + unsigned int qqnum=qnum; + while(side_hop(qqnum,lsn)) + { + // switch this face + lsn = (lsn+2)%4; + if (!is_oriented_side(qqnum,lsn)) + orient_side(qqnum,lsn); + else + //We've found a + //cycle.. and + //oriented all + //quads in it. + break; + } + } + } + } + + + void + GridReordering::orient_side(const unsigned int quadnum, + const unsigned int localsidenum) + { + MQuad &quad = mquads[quadnum]; + int op_side_l = (localsidenum+2)%4; + MSide &side = sides[mquads[quadnum].side[localsidenum]]; + const MSide &op_side =sides[mquads[quadnum].side[op_side_l]]; + + //is the opposite side oriented? + if (op_side.Oriented) + { + //YES - Make the orientations match + //Is op side in default orientation? + if (op_side.v0==quad.v[ConnectGlobals::DefaultOrientation[op_side_l][0]]) + { + //YES + side.v0=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][0]]; + side.v1=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][1]]; + } + else + { + //NO, its reversed + side.v0=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][1]]; + side.v1=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][0]]; + } + } + else + { + //NO + //Just use the default orientation + side.v0=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][0]]; + side.v1=quad.v[ConnectGlobals::DefaultOrientation[localsidenum][1]]; + } + side.Oriented=true; + } + + + + bool + GridReordering::is_fully_oriented_quad(const unsigned int quadnum) const + { + return ( + (sides[mquads[quadnum].side[0]].Oriented)&& + (sides[mquads[quadnum].side[1]].Oriented)&& + (sides[mquads[quadnum].side[2]].Oriented)&& + (sides[mquads[quadnum].side[3]].Oriented) + ); + } + + + + bool + GridReordering::is_oriented_side(const unsigned int quadnum, + const unsigned int lsn) const + { + return (sides[mquads[quadnum].side[lsn]].Oriented); + } + + + + + bool + GridReordering::get_unoriented_quad(unsigned int &UnOrQLoc) const + { + while( (UnOrQLoc(-1)) + { + qnum = opquad; + return true; + } + + return false; + } + + + void + GridReordering::get_quads (std::vector > &outquads) const + { + outquads.clear(); + outquads.reserve(mquads.size()); + for(unsigned int qn=0;qn q = mquads[qn].original_cell_data; + + //Are the sides oriented? + assert(is_fully_oriented_quad(qn)); + bool s[4]; //whether side 1 ,2, 3, 4 are in the default orientation + for(int sn=0;sn<4;sn++) + { + s[sn]=is_side_default_oriented(qn,sn); + } + // Are they oriented in the "deal way"? + assert(s[0]==s[2]); + assert(s[1]==s[3]); + // How much we rotate them by. + int rotn = 2*(s[0]?1:0)+ ((s[0]^s[1])?1:0); + + for(int i=0;i<4;++i) + { + q.vertices[(i+rotn)%4]=mquads[qn].v[i]; + } + outquads.push_back(q); + } + + } + + bool + GridReordering::is_side_default_oriented (const unsigned int qnum, + const unsigned int lsn) const + { + return (sides[mquads[qnum].side[lsn]].v0 == + mquads[qnum].v[ConnectGlobals::DefaultOrientation[lsn][0]]); + } + } // namespace GridReordering2 +} // namespace internal + + +void GridReordering<2>::reorder_cells (std::vector > &original_cells) { - // there should not be much to do - // in 1d... + internal::GridReordering2d::GridReordering().reorient(original_cells); } #endif - // explicit instantiations. only require the main function, it should // then claim whatever templates it needs. note that in 1d, the -// respective function is already specialized -#if deal_II_dimension >= 2 +// respective function is already specialized, and in 2d we have an +// explicit specialization of the whole class +#if deal_II_dimension == 3 template void GridReordering::