/**
* Class declaring some dimension dependent numbers which are needed
* for the grid reordering class. This is the specialization for the
- * 2d case.
+ * 3d case.
*
* @author Wolfgang Bangerth, 2000
*/
{
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.
*/
- static const unsigned int rotational_states_of_cells = static_cast<unsigned int>(-1);
static const unsigned int rotational_states_of_faces = 4;
};
* 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
+ * 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 -- somtimes
+ * 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):
* @begin{verbatim}
- * 12----13----14----15
+ * 13----------14----15
+ * | \ | |
+ * | \ 4 | 5 |
+ * | \ | |
+ * | 12-----10----11
* | | | |
- * | 5 | 7 | 6 |
+ * | | | 7 |
* | | | |
- * 8-----9-----10----11
+ * | 3 | 8-----9
* | | | |
- * | 3 | | 4 |
+ * | | | 6 |
* | | | |
* 4-----5-----6-----7
* | | | |
- * | 0 | 1 | 2 |
+ * | 2 | 1 | 0 |
* | | | |
* 0-----1-----2-----3
* @end{verbatim}
- */
-/*
* Note that there is a hole in the middle. Assume now that the user
- * described the first cell 0 by the vertex numbers @p{0 1 5 4}, cell
- * 5 by @p{12 8 9 13}, and cell 6 by @p{10 11 15 14}. All other cells
- * are numbered in the usual way, i.e. starting at the bottom left and
- * counting counterclockwise. Cell 5 therefore is the only one that
- * does not follow this order; however, note that the bottom line of
- * cell 5 given by this order of cell 5 does match with the top line
- * of cell 4 in that orientation. Given this description of cells, the
- * algorithm will start with cell zero and add one cell after the
- * other, up until the sixth one. Then the situation will be the
- * following:
+ * described the first cell 0 by the vertex numbers @p{2 3 7 6}, and
+ * cell 5 by @p{15 14 10 11}, and assume that cells 1, 2, 3, and 4 are
+ * numbered such that 5 can be added in initial rotation. All other
+ * cells are numbered in the usual way, i.e. starting at the bottom
+ * left and counting counterclockwise. Given this description of
+ * cells, the algorithm will start with cell zero and add one cell
+ * after the other, up until the sixth one. Then the situation will be
+ * the following:
* @begin{verbatim}
- * 12->--13----14->--15
+ * 13----->---14--<--15
+ * | \ | |
+ * | > 4 v 5 v
+ * | \ | |
+ * | 12->--10--<--11
* | | | |
- * v 5 v 7 ^ 6 ^
+ * ^ | | 7 |
* | | | |
- * 8-->--9-----10->--11
+ * | 3 ^ 8-->--9
* | | | |
- * ^ 3 ^ ^ 4 ^
+ * | | ^ 6 ^
* | | | |
* 4-->--5-->--6-->--7
* | | | |
- * ^ 0 ^ 1 ^ 2 ^
+ * ^ 2 ^ 1 ^ 0 ^
* | | | |
* 0-->--1-->--2-->--3
* @end{verbatim}
- *
- * Coming now to cell 7, we see that the two
- * opposite lines to its left and right have different directions; we
- * will therefore find no orientation of cell 7 in which it can be
- * added without violation of the consistency of the
- * triangulation. According to the rule stated above, we track back to
- * the neighbor with greatest index, which is cell 6......
+ * Coming now to cell 7, we see that the two opposite lines at its top
+ * and bottom have different directions; we will therefore find no
+ * orientation of cell 7 in which it can be added without violation of
+ * the consistency of the triangulation. According to the rule stated
+ * above, we track back to the neighbor with greatest index, which is
+ * cell 6, but since its bottom line is to the right, its top line
+ * must be to the right as well, so we won't be able to find an
+ * orientation of cell 6 such that 7 will fit into the
+ * triangulation. Then, if we have finished all possible orientations
+ * of cell 6, we track back to the neighbor of 6 with the largest
+ * index and which has been added already. This would be cell
+ * 0. However, we know that the orientation of cell 0 can't be
+ * important, so we conclude that there is no possible way to orient
+ * all the lines of the given cells such that they satisfy the
+ * requirements if deal.II triangulations. We know that this can't be,
+ * so it results in an exception be thrown.
+ *
+ * The bottom line of this example is that when we looked at all
+ * possible orientations of cell 6, we couldn't find one such that
+ * cell 7 could be added, and then decided to track back to cell 0. We
+ * did not even attempt to turn cell 5, after which it would be simple
+ * to add cell 7. Thus, the algorithm described above has to be
+ * modified: we are only allowed to track back to that neighbor that
+ * has already been added, with the largest cell index, if we fail to
+ * add a cell in any orientation. If we track back further because we
+ * have exhausted all possible orientations but could add the cell
+ * (i.e. we track back since another cell, further down the road
+ * couldn't be added, irrespective of the orientation of the cell
+ * which we are presently considering), then we are not allowed to
+ * track back to one of its neighbors, but have to track back only one
+ * cell index.
*
* The second method to prune the tree is that usually we cannot add a
* new cell since the orientation of one of its neighbors that have
*
* These two methods have proven extremely efficient. We have been
* able to read very large grids (several ten thousands of cells)
- * without the need to backtrack much. In particular, the time to find
+ * without the need to track back much. In particular, the time to find
* an ordering of the cells was found to be mostly linear in the
* number of cells, and the time to reorder them is usually much
* smaller (for example by one order of magnitude) than the time
* comparison.
*/
bool operator < (const Face &face) const;
+
+ /**
+ * Check for equality of
+ * vertex indices.
+ */
+ bool operator == (const Face &face) const;
};
* in 1d, and that there are four
* possible orientations of a
* cell in 2d, which are numbered
- * counter-clockwise. The
- * function is presently not
- * implemented for 3d.
+ * counter-clockwise. In 3d, a
+ * cell can be rotated in 24
+ * ways.
*
* Note that in general the
* result will be different if
void rotate (const unsigned int times);
/**
- * Rotations are not possible in 1d.
+ * Rotations are not possible in
+ * 1d.
*/
DeclException0 (ExcNotPossible);
/**
#include <algorithm>
+#include <base/thread_management.h>
#include <grid/grid_reordering.h>
new_faces_ptr[rot][face]->second.adjacent_cells[1] = cell_no;
};
- // still have to link cell in
- // different orientations to the
- // faces
+ // 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).
+ Threads::ThreadMutex initialization_lock;
+ initialization_lock.acquire ();
+
+ static bool already_initialized = false;
- Assert (false, ExcNotImplemented());
+ // 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 pair<unsigned int, unsigned int>
+ cell_orientation_faces[rotational_states_of_cells][GeometryInfo<dim>::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<dim> standard_cell;
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+ standard_cell.vertices[v] = 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<dim>::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<GeometryInfo<dim>::faces_per_cell; ++face)
+ {
+ unsigned int f, r;
+ for (f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ for (r=0; r<rotational_states_of_cells; ++r)
+ if (standard_faces[face] == new_faces[f][r])
+ {
+ cell_orientation_faces[rot][face] = make_pair(f,r);
+ break;
+ };
+
+ // make sure that we
+ // have found something
+ // indeed
+ Assert ( ! ((f == GeometryInfo<dim>::faces_per_cell) &&
+ (r == rotational_states_of_cells)),
+ 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
+ vector<bool> face_used(GeometryInfo<dim>::faces_per_cell, false);
+ for (unsigned int face=0; face<GeometryInfo<dim>::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; rotational_states_of_cells; ++rot)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ faces[rot][face] = new_faces_ptr
+ [cell_orientation_faces[rot][face].second]
+ [cell_orientation_faces[rot][face].first];
};
#endif
+template <int dim>
+bool GridReordering<dim>::Face::operator == (const Face &face) const
+{
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+ if (vertices[v] != face.vertices[v])
+ return false;
+ return true;
+};
+
+
+
template <int dim>
GridReordering<dim>::FaceData::FaceData () :
use_count (0)
template <>
-void CellData<1>::rotate (unsigned int)
+void CellData<1>::rotate (const unsigned int)
{
Assert (false, ExcNotPossible());
};
template <>
-void CellData<2>::rotate (unsigned int times)
+void CellData<2>::rotate (const unsigned int times)
{
- while (times != 0)
+ Assert (times < 4, ExcInvalidRotation(times));
+
+ for (unsigned int i=0; i<times; ++i)
{
const unsigned int x = vertices[0];
vertices[0] = vertices[1];
vertices[1] = vertices[2];
vertices[2] = vertices[3];
vertices[3] = x;
-
- --times;
};
};
template <>
-void CellData<3>::rotate (unsigned int)
+void CellData<3>::rotate (const unsigned int times)
{
- Assert (false, ExcNotImplemented());
+ Assert (times < 24, ExcInvalidRotation(times));
+
+ // list the 24 ways to rotate a
+ // cell, by the ways by which the
+ // vertices are then permuted
+ static unsigned int rotations[24][GeometryInfo<3>::vertices_per_cell]
+ = { {0,1,2,3,4,5,6,7},
+ {1,5,6,2,0,4,7,3},
+ {5,4,7,6,1,0,3,2},
+ {4,0,3,7,5,1,2,6},
+ {2,3,0,1,6,7,4,5},
+ {3,7,4,0,2,6,5,1},
+ {7,6,5,4,3,2,1,0},
+ {6,2,1,5,7,3,0,4},
+ {0,3,7,4,1,2,6,5},
+ {3,2,6,7,0,1,5,4},
+ {2,1,6,5,3,0,7,4},
+ {1,0,4,5,2,3,7,6},
+ {1,2,3,0,5,6,7,4},
+ {2,6,7,3,1,5,4,0},
+ {6,5,4,7,2,1,0,3},
+ {5,1,0,4,6,2,3,7},
+ {5,6,2,1,4,7,3,0},
+ {6,7,3,2,5,4,0,1},
+ {7,4,0,3,6,5,1,2},
+ {4,5,1,0,7,6,2,3},
+ {3,0,1,2,7,4,5,6},
+ {0,4,5,1,3,7,6,2},
+ {4,7,6,5,0,3,2,1},
+ {7,3,2,6,4,0,1,5} };
+
+ // first copy over the old vertex
+ // numbers
+ const unsigned int old_vertices[GeometryInfo<3>::vertices_per_cell]
+ = { vertices[0], vertices[1], vertices[2], vertices[3],
+ vertices[4], vertices[5], vertices[6], vertices[7] };
+
+ // then copy them back in the new
+ // order
+ for (unsigned int i=0; i<GeometryInfo<3>::vertices_per_cell; ++i)
+ vertices[i] = old_vertices[rotations[times][i]];
};
projects / dealii-svn.git / commitdiff