From 56e00ba65552cfd5bcfb9cfbe983ce99cf694b5c Mon Sep 17 00:00:00 2001 From: mike Date: Tue, 23 Sep 2003 22:11:56 +0000 Subject: [PATCH] Check in first version of Mike's fabulous 3d reordering algorithm. git-svn-id: https://svn.dealii.org/trunk@8011 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/grid/grid_reordering.h | 863 +------ .../include/grid/grid_reordering_internal.h | 636 ++++++ .../deal.II/source/grid/grid_reordering.cc | 1985 +++++++---------- 3 files changed, 1478 insertions(+), 2006 deletions(-) create mode 100644 deal.II/deal.II/include/grid/grid_reordering_internal.h diff --git a/deal.II/deal.II/include/grid/grid_reordering.h b/deal.II/deal.II/include/grid/grid_reordering.h index e839f9b72e..7a71c0468b 100644 --- a/deal.II/deal.II/include/grid/grid_reordering.h +++ b/deal.II/deal.II/include/grid/grid_reordering.h @@ -17,429 +17,9 @@ #include #include -#include #include -namespace internal -{ -/** - * 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 - * 3d case. - * - * @author Wolfgang Bangerth, 2000 - */ - template <> - class GridReorderingInfo<3> - { - public: - /** - * Number of possible - * orientations of a cell in 3d, - * i.e. in how many ways can we - * arrange the eight vertices of - * a cube such that the numbering - * conventions of deal.II are - * satisfied. Since we can put - * the cube on each of its six - * faces and rotate it into four - * ways on each of these faces, - * the total number is 24. - */ - static const unsigned int rotational_states_of_cells = 24; - - /** - * Number of possible - * orientations of a face in - * 2d. It is the quadrilateral - * rotated in all four - * orientations, and also viewed - * from the back, i.e. vertices - * not numbered in - * counter-clockwise sense but in - * the other direction. - */ - static const unsigned int rotational_states_of_faces = 8; - }; - -} - - - - - -namespace internal -{ -/** - * Implement the algorithm described in the documentation of the - * GridReordering<2> class. - * - * @author Michael Anderson, 2003 - */ - namespace GridReordering2d - { - - /** - * Check whether a given - * arrangement of cells is - * already consisten. If this is - * the case, then we skip the - * reordering pass. - * - * This function works by looping - * over all cells, checking - * whether one of its faces - * already exists in a list of - * edges, and if it already - * exists in reverse order, then - * return @p{false}. If it is not - * already in the list, or in the - * correct direction, then go on - * with the next faces or cell. - */ - bool - is_consistent (const std::vector > &cells); - - -/** - * Defines a variety of variables related to the connectivity of a - * simple quad element. This includes the nodes on each edge, which - * edges come into each node and what the default deal.II directions - * are for the quad. - * - * @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 - - /** @@ -1006,425 +586,6 @@ namespace internal template class GridReordering { - 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); - - private: - - /** - * Forward declarations of local - * classes. - */ - class Cell; - class Face; - class FaceData; - - /** - * Typedef for a stack type that - * describes the rotational - * states of all cells that have - * already been fitted into the - * grid. It is mostly used like a - * stack, but sometimes we need - * random access into values - * below the top, so we can't use - * the @p{stack} adaptor from - * STL. - */ - typedef std::vector RotationStack; - - /** - * 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 - */ - class Cell : public CellData - { - public: - /** - * 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. - * - * Note that the order in - * which the faces for a - * specific rotational state - * appear is not specified, - * as this is not - * important. It is only - * important that each face - * in one rotational state or - * other appears once for - * each orientation of the - * cell. - */ - typename std::map::iterator - faces[internal::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); - - /** - * Copy constructor to work - * around a bug in a gcc3.0 - * snapshot. - */ - Cell (const Cell &); - - /** - * 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 (std::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. The size of the array - * is equal to - * @p{GeometryInfo::vertices_per_face}. - */ - unsigned int vertices[1<<(dim-1)]; - - /** - * Comparison operator. Use - * the vertex indices as - * primary, secondary, - * ... criteria for - * comparison. - */ - bool operator < (const Face &face) const; - - /** - * Check for equality of - * vertex indices. - */ - 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. - * - * Note that the order in - * which the reverted faces - * appear here is not - * specified, as it is not - * important for the - * algorithm. - */ - typename std::map::const_iterator - reverse_faces[internal::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 (std::vector &cells, - RotationStack &rotation_states, - const unsigned int track_back_to_cell); - - static bool try_rotate_single_neighbors (std::vector &cells, - RotationStack &rotation_states); - - /** - * 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 (std::vector &cells, - std::vector > &original_cells, - const std::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 - std::vector - presort_cells (std::vector &cells, - std::map &faces); - - /** - * By the resolution of Defect - * Report 45 to the ISO C++ 1998 - * standard, nested classes - * automatically have access to - * members of the enclosing - * class. Nevertheless, some - * compilers don't implement this - * resolution yet, so we have to - * make them @p{friend}, which - * doesn't hurt on the other - * compilers as well. - */ - friend class Cell; - friend class Face; - friend class FaceData; }; @@ -1470,10 +631,28 @@ class GridReordering<2> -/* -------------- declaration of explicit specializations ------------- */ - +/** + * This specialization of the general template implements the + * 3d-algorithm described in the documentation of the general + * template. + * + * @author Michael Anderson, 2003 + */ template <> -void GridReordering<3>::Cell::insert_faces (std::map &global_faces); +class GridReordering<3> +{ + 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, + std::vector > * node_vec_ptr=NULL); +}; + + #endif diff --git a/deal.II/deal.II/include/grid/grid_reordering_internal.h b/deal.II/deal.II/include/grid/grid_reordering_internal.h new file mode 100644 index 0000000000..bf4f190fcc --- /dev/null +++ b/deal.II/deal.II/include/grid/grid_reordering_internal.h @@ -0,0 +1,636 @@ +//---------------------------- grid_reordering_internal.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 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_internal.h --------------------------- +#ifndef __deal2__grid_reordering_internal_h +#define __deal2__grid_reordering_internal_h + + +#include +#include + +#include +#include + + + + + +namespace internal +{ +/** + * Implement the algorithm described in the documentation of the + * GridReordering<2> class. + * + * @author Michael Anderson, 2003 + */ + namespace GridReordering2d + { + + /** + * Check whether a given + * arrangement of cells is + * already consisten. If this is + * the case, then we skip the + * reordering pass. + * + * This function works by looping + * over all cells, checking + * whether one of its faces + * already exists in a list of + * edges, and if it already + * exists in reverse order, then + * return @p{false}. If it is not + * already in the list, or in the + * correct direction, then go on + * with the next faces or cell. + */ + bool + is_consistent (const std::vector > &cells); + + +/** + * Defines a variety of variables related to the connectivity of a + * simple quad element. This includes the nodes on each edge, which + * edges come into each node and what the default deal.II directions + * are for the quad. + * + * @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 + + +/** + * Implement the algorithm described in the documentation of the + * GridReordering<2> class. + * + * @author Michael Anderson, 2003 + */ + namespace GridReordering3d + { + + /** + * During building the conectivity information we + * dont need all the heavy duty information about + * edges that we will need later. So we can save + * memory and time by using these light-weight edges. + **/ + class CheapEdge + { + public: + //! The first node + int node0; + //! The second node + int node1; + //! A simple constructor + CheapEdge(int n0, int n1); + //! Need a partial ordering for the STL + bool operator<(const CheapEdge & e2) const; + }; + + + class ElementInfo + { + public: + // The numbers of the edges coming into node i + // are given by edge_to_node[i][k] where k=0,1,2 + int edge_to_node[8][3]; + + //the orientation of edge coming into node i + //is given by edge_to_node_orient[i][k] where k=0,1,2 + // 1 means the given node is the start of the edge + // -1 means the end of the edge + int edge_to_node_orient[8][3]; + + // nodesonedge[i][0] is the start node for edge i + // nodesonedge[i][1] is the end node for edge i + int nodes_on_edge[12][2]; + int nodes_on_face[6][4]; +}; + +class DealElemInfo : public ElementInfo +{ + public: + DealElemInfo(); +}; + + + + + //! A conectivity and orientation aware edge class + class Edge + { + public: + //! Simple constructor + Edge(int n0, int n1, int orient=0): orientation_flag(orient), group(0), + num_neighbouring_cubes(0), neighbouring_cubes(NULL) + {nodes[0]=n0; nodes[1]=n1;}; + + //! Simple Destructor + ~Edge(); + + //! The IDs for the end nodes + + int nodes[2]; + /** + * Whether the edge has been oriented (0), + * points from node 0 to node 1 (1), + * or the reverse (-1) + **/ + int orientation_flag; + + /** + * Used to determine which "sheet" of parallel edges the edge falls in + * when oriented. 0 means not yet decided. + **/ + int group; + + unsigned int num_neighbouring_cubes; + unsigned int * neighbouring_cubes; + }; + + //! A conectivity an orientation aware cell + /** + * The connectivity of the cell is not contained within + * (This was for flexability in using deal's ordering of + * edges or the XDA format etc) For this information we + * need the ElemInfo class. + * + * One thing we do know is that the first four edges in + * the edge class are parallel, as are the second four, + * and the third four. + * + * \todo TODO: Need to move conectivity information out + * of cell and into edge. + **/ + class Cell + { + public: + int edges[12]; //!< The IDs for each of the edges + int nodes[8]; //!< The IDs for each of the nodes + //! Which way do the edges point. + /** + * Whether node 0 of the edge is the base of + * the edge in local element (1) or node 1 is + * the base (-1) + **/ + int local_orientation_flags[12]; + /** + * An internal flag used to determine whether the cell is in the + * queue of cells to be orriented in the current sheet. + **/ + bool waiting_to_be_processed; + + //! Copy Constructor + Cell(const Cell& c) + { + for(int i=0;i<12;++i) + { + edges[i]=c.edges[i]; + local_orientation_flags[i]=c.local_orientation_flags[i]; + } + for(int i=0;i<8;++i) + { + nodes[i]=c.nodes[i]; + } + waiting_to_be_processed=c.waiting_to_be_processed; + } + + //! Default Constructor + Cell() + { + for(int i=0;i<12;++i) + { + edges[i]=-1; + local_orientation_flags[i]=1; + } + for(int i=0;i<8;++i) + { + nodes[i]=-1; + } + waiting_to_be_processed=false; + } + }; + + + //! This holds all the pieces for orientation together.. + /** + * Contains lists of nodes, edges and cells. + * As well as the information about how they all connect together. + **/ + class Mesh + { + public: + //! Information about how a cell is built up from nodes and edges. + const ElementInfo & info; + + //! The list of nodes + std::vector< Point<3> > node_list; + //! The list of edges + std::vector edge_list; + //! The list of cells + std::vector cell_list; + + //! Checks whether every cell in the mesh is sensible + /** + * By calling sanity_check(cell_num) on every cell. + **/ + bool sanity_check() const; + //! Checks that every node matches with its edges + /** + * By calling sanity_check(cell_num,node_num) for each node + **/ + bool sanity_check(int cell_num) const; + //! Checks that each edge going into a node is correctly setup + bool sanity_check_node(int cell_num, int i) const; + + //! Default Constructor + Mesh(const ElementInfo & INFO): info(INFO) { } + + //! Prints all information about the mesh + void dump() const; + //! Prints all information about the cell + void dump_cell(const Cell &c) const; + //! Writes edge information to a file. + void dump_edges(char const * const fname) const; + + private: + //! Unimplemented private copy constructor to disable it. + Mesh(const Mesh&); + //! Unimplemented private assignemnet operator to disable it. + Mesh& operator=(const Mesh&); + }; + + + class Orienter + { + public: + + Orienter(){}; + + //The cube we're looking at now. + unsigned int cur_posn; + //We have fully oriented all cubes before this one. + unsigned int marker_cube; + + std::vector SheetToProcess; + + int cur_edge_group; + + bool edge_orient_array[12]; + + bool orient_edges(Mesh &m); + void orient_cubes(Mesh &m); + + bool GetNextUnorientedCube(Mesh &m); + bool is_oriented(const Mesh &m, int cell_num); + + bool OrientEdgesInCurrentCube(Mesh &m); + bool OrientEdgeSetInCurrentCube(Mesh &m, int edge_set); + bool OrientNextUnorientedEdge(Mesh &m); + bool Consistant(Mesh &m, int cell_num); + + + void GetAdjacentCubes(Mesh &m); + bool GetNextActiveCube(Mesh &m); + + bool CheckCellEdgeGroupConsistancy(const Mesh &m, const Cell &c) const; + bool CheckCellEdgeGroupConsistancy(const Mesh &m, const Cell & c, int egrp) const; + + }; + + + //! Creates the connectivity information for the mesh m. + void build_mesh(Mesh &m); + + + + + } // namespace GridReordering3d +} // namespace internal + + + +#endif diff --git a/deal.II/deal.II/source/grid/grid_reordering.cc b/deal.II/deal.II/source/grid/grid_reordering.cc index 67a45fad26..ea7bcf0db3 100644 --- a/deal.II/deal.II/source/grid/grid_reordering.cc +++ b/deal.II/deal.II/source/grid/grid_reordering.cc @@ -11,30 +11,16 @@ // //---------------------------- grid_reordering.cc --------------------------- -#include #include +#include -#include #include +#include +#include +#include -namespace internal -{ -// static variables -#if deal_II_dimension == 3 - const unsigned int GridReorderingInfo<3>::rotational_states_of_cells; - const unsigned int GridReorderingInfo<3>::rotational_states_of_faces; -#endif - -} - -template -const unsigned int GridReordering::Cell::invalid_neighbor; - -template -const unsigned int GridReordering::FaceData::invalid_adjacent_cell; - #if deal_II_dimension == 1 @@ -49,1143 +35,6 @@ void GridReordering<1>::reorder_cells (const std::vector > &) -#if deal_II_dimension == 3 - -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 -GridReordering::Cell::Cell (const Cell &c) : - CellData (c), - cell_no(c.cell_no), - track_back_to_cell(c.track_back_to_cell) -{ - for (unsigned int i=0; i::rotational_states_of_cells; ++i) - for (unsigned int j=0; j::faces_per_cell; ++j) - faces[i][j]=c.faces[i][j]; - - for (unsigned int i=0; i::faces_per_cell; ++i) - neighbors[i]=c.neighbors[i]; -} - - -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 (std::map &/*global_faces*/) -{ - Assert (false, ExcNotImplemented()); -} - - - -template <> -void -GridReordering<3>::Cell::insert_faces (std::map &global_faces) -{ - const unsigned int dim = 3; - - // first generate for each of the 6 - // faces of a cell in 3d the four - // possible orientations and - // cross-link them among each other - // - // do this generation step by first - // only inserting each face in - // standard orientation and then - // fill in the other ones by - // rotation of these faces - // - // note that we have the indices - // reversed here compared to the - // Cell class, for simplicity - const Face new_faces_tmp[GeometryInfo::faces_per_cell] - = { { { this->vertices[0], this->vertices[1], - this->vertices[2], this->vertices[3] } }, - { { this->vertices[4], this->vertices[5], - this->vertices[6], this->vertices[7] } }, - { { this->vertices[0], this->vertices[1], - this->vertices[5], this->vertices[4] } }, - { { this->vertices[1], this->vertices[5], - this->vertices[6], this->vertices[2] } }, - { { this->vertices[3], this->vertices[2], - this->vertices[6], this->vertices[7] } }, - { { this->vertices[0], this->vertices[4], - this->vertices[7], this->vertices[3] } } }; - Face new_faces[GeometryInfo::faces_per_cell][internal::GridReorderingInfo::rotational_states_of_faces] - = { { new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], - new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] }, - { new_faces_tmp[1], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], - new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] }, - { new_faces_tmp[2], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], - new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] }, - { new_faces_tmp[3], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], - new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] }, - { new_faces_tmp[4], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], - new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] }, - { new_faces_tmp[5], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], - new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0], new_faces_tmp[0] }}; - - // first do the faces in their - // usual direction - for (unsigned int face=0; face::faces_per_cell; ++face) - for (unsigned int rot=1; rot::rotational_states_of_faces/2; ++rot) - for (unsigned int v=0; v::vertices_per_face; ++v) - new_faces[face][rot].vertices[v] - = new_faces[face][0].vertices[(v+rot) % GeometryInfo::vertices_per_face]; - // then do everything as viewed - // from the back. this is simple, - // as we only have to revert - // indices 1 and 3 - for (unsigned int face=0; face::faces_per_cell; ++face) - for (unsigned int rot=internal::GridReorderingInfo::rotational_states_of_faces/2; - rot::rotational_states_of_faces; - ++rot) - { - for (unsigned int v=0; v::vertices_per_face; ++v) - new_faces[face][rot].vertices[v] - = new_faces[face] - [rot-internal::GridReorderingInfo::rotational_states_of_faces/2].vertices[v]; - std::swap (new_faces[face][rot].vertices[1], - new_faces[face][rot].vertices[3]); - }; - - - // now insert all the faces, by now - // without specific order with - // respect to the orientational - // states of the cell. note that we - // get the indices correct - // here. also remark that the face - // might already have been in the - // map, depending on whether a - // newighbor has already inserted - // it or not. we don't care about - // that here, though - std::map::iterator - new_faces_ptr[internal::GridReorderingInfo::rotational_states_of_faces] - [GeometryInfo::faces_per_cell]; - for (unsigned int face=0; face::faces_per_cell; ++face) - for (unsigned int rot=0; - rot::rotational_states_of_faces; - ++rot) - new_faces_ptr[rot][face] - = global_faces.insert (std::make_pair(new_faces[face][rot], FaceData())).first; - - // and crosslink them to each other - for (unsigned int face=0; face::faces_per_cell; ++face) - for (unsigned int rot=0; - rot::rotational_states_of_faces; - ++rot) - for (unsigned int other_rot=0; - other_rot::rotational_states_of_faces; - ++other_rot) - { - if (other_rot < rot) - new_faces_ptr[rot][face]->second.reverse_faces[other_rot] - = new_faces_ptr[other_rot][face]; - else - if (other_rot > rot) - new_faces_ptr[rot][face]->second.reverse_faces[other_rot-1] - = new_faces_ptr[other_rot][face]; - // if rot==other_rot, then - // we need not link this - // cell to itself - }; - - - // for each of the faces (whether - // already inserted or not) note - // that the present cell is one of - // the neighbors - for (unsigned int face=0; face::faces_per_cell; ++face) - { - if (new_faces_ptr[0][face]->second.adjacent_cells[0] == - FaceData::invalid_adjacent_cell) - { - // no, faces had not been - // used before, so we are the - // first adjacent cell - for (unsigned int rot=0; - rot::rotational_states_of_faces; - ++rot) - { - Assert (new_faces_ptr[rot][face]->second.adjacent_cells[0] - == FaceData::invalid_adjacent_cell, - ExcInternalError()); - new_faces_ptr[rot][face]->second.adjacent_cells[0] = cell_no; - }; - } - else - { - // otherwise: cell had been - // entered before, so we are - // the second neighbor - const unsigned int - previous_neighbor = new_faces_ptr[0][face]->second.adjacent_cells[0]; - for (unsigned int rot=0; - rot::rotational_states_of_faces; - ++rot) - { - Assert (new_faces_ptr[rot][face]->second.adjacent_cells[0] == - previous_neighbor, - ExcInternalError()); - Assert (new_faces_ptr[rot][face]->second.adjacent_cells[1] == - FaceData::invalid_adjacent_cell, - ExcInternalError()); - new_faces_ptr[rot][face]->second.adjacent_cells[1] = cell_no; - }; - }; - }; - - - - // we still have to link cell in - // its 24 different orientations to - // the 6 faces in their - // orientations. now, there we - // could hardcode which faces in - // which rotation belong to the - // cell in each direction, but - // there a good reasons not to do - // so: - // - // first, this depends on that we - // know which orientation of the - // cell has which number, but this - // knowledge is hardcoded in the - // function CellData::rotate, so - // hardcoding it here again would - // mean redundancy, and would above - // that mean that we have to update - // two very different place if we - // chose to change one. - // - // second, finding out which face - // belongs to which cell is error - // prone, and one might get it - // wrong. - // - // the solution is: compute it once - // this function is first called - // based on the information from - // CellData::rotate and use that - // data in following calls to this - // function. the computed data has, - // of course, to be a static member - // function, and we store whether - // the data has been initialized - // already by checking the value of - // a special flag. furthermore, we - // guard the initialization by a - // thread mutex to make it - // thread-safe (in case someone - // wanted to read in two grids at - // the same time, for whatever - // reason). - static Threads::ThreadMutex initialization_lock; - initialization_lock.acquire (); - - static bool already_initialized = false; - - // for each orientation of the - // cell, store in which orientation - // each of the six faces build the - // cell (store which face and which - // orientation): - static std::pair - cell_orientation_faces[internal::GridReorderingInfo::rotational_states_of_cells][GeometryInfo::faces_per_cell]; - - if (already_initialized == false) - { - for (unsigned int rot=0; rot::rotational_states_of_cells; ++rot) - { - // initialize a standard - // cell with the vertex - // numbers of the present - // cell we are working on - CellData standard_cell; - for (unsigned int v=0; v::vertices_per_cell; ++v) - standard_cell.vertices[v] = this->vertices[v]; - - // then rotate it the given - // number of times - standard_cell.rotate (rot); - - // then create the six - // faces of the thus - // rotated cell - const Face standard_faces[GeometryInfo::faces_per_cell] - = { { { standard_cell.vertices[0], standard_cell.vertices[1], - standard_cell.vertices[2], standard_cell.vertices[3] } }, - { { standard_cell.vertices[4], standard_cell.vertices[5], - standard_cell.vertices[6], standard_cell.vertices[7] } }, - { { standard_cell.vertices[0], standard_cell.vertices[1], - standard_cell.vertices[5], standard_cell.vertices[4] } }, - { { standard_cell.vertices[1], standard_cell.vertices[5], - standard_cell.vertices[6], standard_cell.vertices[2] } }, - { { standard_cell.vertices[3], standard_cell.vertices[2], - standard_cell.vertices[6], standard_cell.vertices[7] } }, - { { standard_cell.vertices[0], standard_cell.vertices[4], - standard_cell.vertices[7], standard_cell.vertices[3] } } }; - - // then try to identify - // these faces in the ones - // we have already created - for (unsigned int face=0; face::faces_per_cell; ++face) - { - bool face_found = false; - for (unsigned int f=0; - (!face_found) && (f::faces_per_cell); ++f) - for (unsigned int r=0; - r::rotational_states_of_faces; - ++r) - if (standard_faces[face] == new_faces[f][r]) - { - cell_orientation_faces[rot][face] = std::make_pair(f,r); - face_found = true; - break; - }; - - // make sure that we - // have found something - // indeed - Assert (face_found == true, ExcInternalError()); - }; - - // more checks: make sure - // that each of the - // original faces appears - // in one rotation or other - // as face of the present - // cell in its orientation - // we currently check. as - // we don't call this part - // of the program too - // often, don't make - // differences between - // debug and optimized mode - std::vector face_used(GeometryInfo::faces_per_cell, false); - for (unsigned int face=0; face::faces_per_cell; ++face) - { - // ups, face already - // used? can't be! - Assert (face_used[face] == false, ExcInternalError()); - face_used[face] = true; - }; - // since we have checked - // that each face has not - // been used previously, we - // also know that all faces - // have been used exactly - // once, so no more checks - // necessary - }; - - // that's it: we now know which - // faces build up this cell in - // each of its possible - // orientations - already_initialized = true; - }; - // initialization is done, so - // release the lock and let other - // threads run - initialization_lock.release (); - - // now we can use the information: - // link the faces in their - // directions to the cell in each - // of its orientations - for (unsigned int rot=0; - rot::rotational_states_of_cells; - ++rot) - for (unsigned int face=0; face::faces_per_cell; ++face) - faces[rot][face] = new_faces_ptr - [cell_orientation_faces[rot][face].second] - [cell_orientation_faces[rot][face].first]; -} - - - -template -void GridReordering::Cell::fix_cell_neighbors () -{ - for (unsigned int face=0; face::faces_per_cell; ++face) - { - // 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 - track_back_to_cell = cell_no-1; - if (cell_no == 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 -{ - // 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_rot::rotational_states_of_faces-1; - ++face_rot) - { - const FaceData &reverse_face = face.reverse_faces[face_rot]->second; - 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) -{ - 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::Cell::mark_faces_unused (const unsigned int rot) -{ - 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 -bool GridReordering::Face::operator == (const Face &face) const -{ - for (unsigned int v=0; v::vertices_per_face; ++v) - if (vertices[v] != face.vertices[v]) - return false; - return true; -} - - - -template -GridReordering::FaceData::FaceData () : - use_count (0) -{ - adjacent_cells[0] = adjacent_cells[1] = invalid_adjacent_cell; -} - - - - - - -template -inline -void GridReordering::track_back (std::vector &cells, - RotationStack &rotation_states, - const unsigned tbtc) -{ - unsigned int track_back_to_cell = tbtc; - - top_of_function: - - Assert (track_back_to_cell > 0, ExcInternalError()); - - unsigned int last_rotation_state = static_cast(-1); - 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.back(); - - // 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_back(); - }; - Assert (last_rotation_state < internal::GridReorderingInfo::rotational_states_of_cells, ExcInternalError()); - - // 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 < internal::GridReorderingInfo::rotational_states_of_cells-1) - { - // possible. push that state to - // the stack and leave - rotation_states.push_back (last_rotation_state+1); - return; - } - else - { - // last cell can't be rotated - // further. go on with - // backtracking - const typename std::vector::iterator - try_cell = cells.begin() + rotation_states.size(); - - 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 -bool GridReordering::try_rotate_single_neighbors (std::vector &cells, - RotationStack &rotation_states) -{ - // the rotation state of the cell - // which we try to add by rotating - // neighbors has already been - // popped from the stack, so we get - // its number like this: - const unsigned int cell_no = rotation_states.size(); - - // now try each of the neighbors - // that have already been added to - // the grid. don't try the cell - // that we will track back to - // anyway if this operation should - // fail - for (unsigned int neighbor=0; neighbor::faces_per_cell; ++neighbor) - if (cells[cell_no].neighbors[neighbor] < cell_no) - if (cells[cell_no].neighbors[neighbor] != cells[cell_no].track_back_to_cell) - { - const unsigned int neighbor_no = cells[cell_no].neighbors[neighbor]; - const unsigned int old_rotation_state = rotation_states[neighbor_no]; - - // unlink faces used by the - // present rotation state - cells[neighbor_no].mark_faces_unused (old_rotation_state); - - // then try all rotation - // states besides the ones - // that have already been - // tried: - for (unsigned int neighbor_rot=old_rotation_state+1; - neighbor_rot::rotational_states_of_cells; ++neighbor_rot) - { - // first, if the - // neighbor itself does - // not fit in the grid, - // then there is - // nothing to do - if (! cells[neighbor_no].check_consistency (neighbor_rot)) - continue; - - // however, if the - // neighbor worked, - // then mark its faces - // as used - // preliminarily and - // try to fit in the - // present cell in some - // orientation - cells[neighbor_no].mark_faces_used (neighbor_rot); - - for (unsigned int cell_rot=0; cell_rot::rotational_states_of_cells; ++cell_rot) - if (cells[cell_no].check_consistency (cell_rot) == true) - { - // ah, see, - // this - // combination - // of neighbor - // rotation and - // this cell - // works. enter - // the - // respective - // states into - // the arrays - // and leave - // with success - rotation_states[neighbor_no] = neighbor_rot; - - rotation_states.push_back (cell_rot); - cells[cell_no].mark_faces_used (cell_rot); - - return true; - }; - - // no, there was no - // way to fit the - // present cell into - // the grid given - // this orientation - // of the - // neighbor. discard - // this attempt and - // try that neighbors - // next rotation - cells[neighbor_no].mark_faces_unused (neighbor_rot); - }; - - // there was no way to - // rotate this neighbor so - // that the present cell - // fit into the - // grid. reinstantiate the - // old state and go on to - // the next neighbor - cells[neighbor_no].mark_faces_used (old_rotation_state); - }; - - // rotation of neighbors did not - // help this cell, there is no - // other way than to do a full - // backtracking - return false; -} - - - -template -void GridReordering::find_reordering (std::vector &cells, - std::vector > &original_cells, - const std::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 - RotationStack 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_back (0); - cells[0].mark_faces_used (rotation_states.back()); - - 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_back (0); - - // check whether the present - // cell in the present - // orientation is valid - check_topmost_cell: - - const typename std::vector::iterator - try_cell = cells.begin() + rotation_states.size()-1; - if (try_cell->check_consistency (rotation_states.back())) - { - // 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.back()); - - continue; - } - else - { - // no, doesn't work. see if - // we can rotate the top - // cell so that it works - if (rotation_states.back()+1 - < internal::GridReorderingInfo::rotational_states_of_cells) - { - // yes, can be - // done. then do so and - // check again - ++rotation_states.back(); - goto check_topmost_cell; - } - else - { - // no, no more - // orientation of the - // top cell possible, - // we have to backtrack - // some way - - // first pop rotational - // state of top cell, - // since for that no - // faces have been - // marked as used yet - rotation_states.pop_back(); - - // in general, if we - // fail to insert the - // present cell somehow - // into the existing - // part of the grid, - // then we track back - // to the neighbor of - // the failed cell with - // the highest cell - // index below the - // index of the present - // cell. however, - // before we do so, we - // try a simple - // heuristic: if - // rotating single - // neighbors a little - // helps the process - // somewhat: - const bool rotation_helps - = try_rotate_single_neighbors (cells, rotation_states); - - // if rotation helped, - // then go on to the - // next cell. the - // called function has - // already marked the - // respective faces as - // used and has pushed - // the rotation state - // of the present cell - // to the stack - if (rotation_helps == true) - continue; - - // if that failed to - // help, 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 = std::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.back()); - - // to check the correctness of - // the program up to here: - // unmark the cells' faces to - // check whether they have all - // correctly declared they - // use. checking this is done - // in the calling function, as - // only that has direct access - // to the map of faces (this - // function only accesses it - // through pointers stored in - // the cells) - cells[new_cell_number].mark_faces_unused (rotation_states.back()); - - // then delete this rotational - // state as we don't need it - // any more - rotation_states.pop_back (); - }; -} - - - -template -std::vector -GridReordering::presort_cells (std::vector &cells, - std::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; - - - // have an array into which we - // insert the new cells numbers of - // each cell - const unsigned int invalid_cell_number = static_cast(-1); - std::vector new_cell_numbers (cells.size(), invalid_cell_number); - - unsigned int next_free_new_number = 0; - - // loop over each connected part of - // the domain. since the domain may - // consist of different unconnected - // parts, we have to loop until - // there are no more unnumbered - // cells - while (next_free_new_number < cells.size()) - { - // for initialization of - // min_neighbors, go to the - // first cell of this part of - // the domain that has an - // invalid cell number, i.e. has - // not yet been renumbered - for (unsigned int i=0; i cells[i].count_neighbors()) && - (new_cell_numbers[i]==invalid_cell_number)) - { - min_neighbors = cells[i].count_neighbors(); - cell_with_min_neighbors = i; - if (min_neighbors == 1) - // better is not possible - break; - }; - - // have an array of the next - // cells to be numbered (old numbers) - std::vector next_round_cells (1, cell_with_min_neighbors); - - // 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]); - - - // 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 - std::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)); - }; - }; // end of loop over subdomains - - - Assert (std::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: - std::vector new_cells (cells.size()); - for (unsigned int i=0; i::faces_per_cell; ++n) - if (cells[c].neighbors[n] != Cell::invalid_neighbor) - { - Assert (cells[c].neighbors[n] < new_cell_numbers.size(), - ExcIndexRange(cells[c].neighbors[n], 0, - new_cell_numbers.size())); - cells[c].neighbors[n] = new_cell_numbers[cells[c].neighbors[n]]; - }; - }; - - for (typename std::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 (std::vector > &original_cells) -{ - // 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 - std::vector cells; - cells.reserve (original_cells.size()); - for (unsigned int i=0; i 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::iterator i=faces.begin(); i!=faces.end(); ++i) - Assert (i->second.use_count == 0, ExcInternalError()); -} - -#endif // deal_II_dimension == 3 - #if deal_II_dimension == 2 namespace internal @@ -1722,13 +571,821 @@ void GridReordering<2>::reorder_cells (std::vector > &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, and in 2d we have an -// explicit specialization of the whole class + #if deal_II_dimension == 3 -template -void -GridReordering:: -reorder_cells (std::vector > &); -#endif + +namespace internal +{ + namespace GridReordering3d + { + DeclException1 (GridOrientError, + std::string, + << "Grid Orientation Error"<< arg1); + + + //Switched two ints so that v1v2) + { + int t=v1; + v1=v2; + v2=t; + } + } + + CheapEdge::CheapEdge(int n0, int n1) : node0(n0), node1(n1) + { + //sort the entries so that node0(e2.node0)) return false; + if ((node1)<(e2.node1)) return true; + return false; + }; + + + // This is the guts of the matter... + void build_mesh(Mesh &m) + { + const ElementInfo &info = m.info; + std::vector & cell_list = m.cell_list; + std::vector & edge_list = m.edge_list; + + const unsigned int cell_list_length=cell_list.size(); + + + unsigned int NumEdges=0; + // Corectly build the edge list + { + // edge_map stores the edge_number associated with a given CheapEdge + std::map edge_map; + unsigned int ctr=0; + for(unsigned int cur_cell_id=0; + cur_cell_id edge_count(NumEdges,0); + + + //Count every time an edge occurs in a cube. + for(unsigned int cur_cell_id=0; + cur_cell_id cur_cell_edge_list_posn(NumEdges,0); + for(unsigned int cur_cell_id =0; + cur_cell_id=0) + { + std::cout<<":"< & n = node_list[i]; + outfile<::reorder_cells(std::vector >& incubes, std::vector > * node_vec_ptr) +{ + + Assert(incubes.size()!=0, + ExcMessage("List of elements to orient was of zero length")); + + // This keeps track of all the local element conectivity information. + // e.g. what edges come into which nodes and with which orientation + internal::GridReordering3d::DealElemInfo deal_info; + + //This does the real work + internal::GridReordering3d::Orienter orienter; + + // This is the internal store for all global connectivity information + // it starts prety much empty. + internal::GridReordering3d::Mesh the_mesh(deal_info); + + if(node_vec_ptr!=NULL) + { + the_mesh.node_list=*node_vec_ptr; + } + + //Copy the cells into our own internal data format. + const unsigned int numelems=incubes.size(); + for(unsigned int i =0 ; i