};
-/**
- * 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
}
+
+
+
+
+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<v1.
+ * After global orientation could be either way.
+ *
+ * @author Michael Anderson, 2003
+ */
+ struct MSide
+ {
+ /**
+ * Constructor.
+ */
+ MSide (const unsigned int initv0,
+ const unsigned int initv1);
+
+ /**
+ * Return whether the sides
+ * are equal, even if their
+ * ends are reversed.
+ */
+ bool operator==(const MSide& s2) const;
+
+ unsigned int v0;
+ unsigned int v1;
+ unsigned int Q0;
+ unsigned int Q1;
+
+ /**
+ * Local side numbers on quads 0 and 1.
+ */
+ int lsn0, lsn1;
+ bool Oriented;
+
+ /**
+ * This class makes a MSide have v0<v1
+ */
+ struct SideRectify;
+
+ /**
+ * Provides a side ordering,
+ * s1<s2, without assuming
+ * v0<v1 in either of the
+ * sides.
+ */
+ struct SideSortLess;
+ };
+
+
+
+/**
+ * Implement the 2d algorithm for grid reordering described in the
+ * documentation of the @ref{GridReordering} class.
+ *
+ * @author Michael Anderson, 2003
+ */
+ class GridReordering
+ {
+ public:
+
+ /**
+ * Do the work intended by
+ * this class.
+ */
+ void reorient(std::vector<CellData<2> > &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<CellData<2> > &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<CellData<2> > &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<MSide> sides;
+ /**
+ * A list of enriched quads
+ * in the mesh.
+ */
+ std::vector<MQuad> mquads;
+ };
+ } // namespace GridReordering2d
+} // 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.)
* 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}
*
* 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
* 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,
* 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).
* 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
* 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
* | \ | |
* the triangulation from this data using the
* @ref{Triangulation}@p{<dim>::create_triangulation} function.
*
- * @author Wolfgang Bangerth, 2000
+ * @author Wolfgang Bangerth, 2000, Michael Anderson 2003
*/
template <int dim>
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<CellData<1> > &);
+};
+
+
+/**
+ * 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<Face,FaceData> &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<CellData<2> > &original_cells);
+};
+
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
template <>
void GridReordering<3>::Cell::insert_faces (std::map<Face,FaceData> &global_faces);
-template <>
-void GridReordering<1>::reorder_cells (std::vector<CellData<1> > &);
#endif
//
//---------------------------- grid_reordering.cc ---------------------------
-
#include <base/thread_management.h>
#include <grid/grid_reordering.h>
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;
+#if deal_II_dimension == 1
+
+void GridReordering<1>::reorder_cells (const std::vector<CellData<1> > &)
+{
+ // there should not be much to do
+ // in 1d...
+}
+
+#endif
+
+
+
+#if deal_II_dimension == 3
template <int dim>
GridReordering<dim>::Cell::Cell () :
}
-#if deal_II_dimension == 2
-
-template <>
-void
-GridReordering<2>::Cell::insert_faces (std::map<Face,FaceData> &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<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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
[cell_orientation_faces[rot][face].first];
}
-#endif
template <int dim>
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<CellData<1> > &)
+ 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<MSide,void>
+ {
+ void operator() (MSide &s) const
+ {
+ if (s.v0>s.v1)
+ std::swap (s.v0, s.v1);
+ }
+ };
+
+
+ struct MSide::SideSortLess : public std::binary_function<MSide,MSide,bool>
+ {
+ bool operator()(const MSide &s1, const MSide &s2) const
+ {
+ int s1vmin,s1vmax;
+ int s2vmin,s2vmax;
+ if (s1.v0<s1.v1)
+ {
+ s1vmin=s1.v0;
+ s1vmax=s1.v1;
+ }
+ else
+ {
+ s1vmin=s1.v1;
+ s1vmax=s1.v0;
+ }
+ if (s2.v0<s2.v1)
+ {
+ s2vmin=s2.v0;
+ s2vmax=s2.v1;
+ }
+ else
+ {
+ s2vmin=s2.v1;
+ s2vmax=s2.v0;
+ }
+
+ if(s1vmin<s2vmin)
+ return true;
+ if(s1vmin>s2vmin)
+ return false;
+ return s1vmax<s2vmax;
+ }
+ };
+
+
+/**
+ * Returns an MSide corresponding to the
+ * specified side of a deal.II CellData<2> 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<CellData<2>,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<unsigned int>(-1)),Q1(static_cast<unsigned int>(-1)),
+ lsn0(static_cast<unsigned int>(-1)),lsn1(static_cast<unsigned int>(-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<CellData<2>,
+ std::vector<MSide>,
+ MQuad>
+ {
+ MQuad operator()(const CellData<2> &q,
+ const std::vector<MSide> &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<CellData<2> > &quads)
+ {
+ build_graph(quads);
+ orient();
+ get_quads(quads);
+ }
+
+
+ void
+ GridReordering::build_graph (const std::vector<CellData<2> > &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<v1
+ std::for_each(sides.begin(),sides.end(),
+ MSide::SideRectify() );
+
+ //Sort them by Sidevertices.
+ std::sort(sides.begin(),sides.end(),
+ MSide::SideSortLess());
+
+ //Remove duplicates
+ sides.erase(std::unique(sides.begin(),sides.end()),
+ sides.end());
+
+ // Swap trick to shrink the
+ // side vector
+ std::vector<MSide>(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<MQuad>::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<unsigned int>(-1))
+ {
+ ss.Q0=qctr;
+ ss.lsn0=i;
+ }
+ else if (ss.Q1==static_cast<unsigned int>(-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<mquads.size()) &&
+ is_fully_oriented_quad(UnOrQLoc) )
+ UnOrQLoc++;
+ return (UnOrQLoc!=mquads.size());
+ }
+
+
+
+ bool
+ GridReordering::get_unoriented_side (const unsigned int quadnum,
+ unsigned int &lsn) const
+ {
+ const MQuad &mq = mquads[quadnum];
+ if(!sides[mq.side[0]].Oriented)
+ {
+ lsn=0;
+ return true;
+ }
+ if(!sides[mq.side[1]].Oriented)
+ {
+ lsn=1;
+ return true;
+ }
+ if(!sides[mq.side[2]].Oriented)
+ {
+ lsn=2;
+ return true;
+ }
+ if(!sides[mq.side[3]].Oriented)
+ {
+ lsn=3;
+ return true;
+ }
+ return false;
+ }
+
+
+ bool
+ GridReordering::side_hop (unsigned int &qnum, unsigned int &lsn) const
+ {
+ const MQuad &mq=mquads[qnum];
+ const MSide &s = sides[mq.side[lsn]];
+ unsigned int opquad=0;
+ if (s.Q0==qnum)
+ {
+ opquad=s.Q1;
+ lsn =s.lsn1;
+ }
+ else
+ {
+ opquad=s.Q0;
+ lsn=s.lsn0;
+ }
+
+ if (opquad!=static_cast<unsigned int>(-1))
+ {
+ qnum = opquad;
+ return true;
+ }
+
+ return false;
+ }
+
+
+ void
+ GridReordering::get_quads (std::vector<CellData<2> > &outquads) const
+ {
+ outquads.clear();
+ outquads.reserve(mquads.size());
+ for(unsigned int qn=0;qn<mquads.size();++qn)
+ {
+ // initialize CellData object with
+ // previous contents, and the
+ // overwrite all the fields that
+ // might have changed in the
+ // process of rotating things
+ CellData<2> 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<CellData<2> > &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<deal_II_dimension>::
projects / dealii-svn.git / commitdiff