From e473a5ab2de857f8abe6c9fecbdb9ad5efd8617d Mon Sep 17 00:00:00 2001 From: wolf Date: Thu, 5 Oct 2000 08:44:20 +0000 Subject: [PATCH] Implement the 3d part of grid_reordering. git-svn-id: https://svn.dealii.org/trunk@3375 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/grid/grid_reordering.h | 122 ++++++++---- deal.II/deal.II/include/grid/tria.h | 9 +- .../deal.II/source/grid/grid_reordering.cc | 174 +++++++++++++++++- .../source/grid/tria.all_dimensions.cc | 54 +++++- 4 files changed, 309 insertions(+), 50 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_reordering.h b/deal.II/deal.II/include/grid/grid_reordering.h index 415c8969f2..08a8121031 100644 --- a/deal.II/deal.II/include/grid/grid_reordering.h +++ b/deal.II/deal.II/include/grid/grid_reordering.h @@ -67,7 +67,7 @@ class GridReorderingInfo<2> /** * 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 */ @@ -76,9 +76,27 @@ 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. */ - static const unsigned int rotational_states_of_cells = static_cast(-1); static const unsigned int rotational_states_of_faces = 4; }; @@ -381,63 +399,91 @@ class GridReorderingInfo<3> * 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 @@ -455,7 +501,7 @@ class GridReorderingInfo<3> * * 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 @@ -716,6 +762,12 @@ class GridReordering : private GridReorderingInfo * comparison. */ bool operator < (const Face &face) const; + + /** + * Check for equality of + * vertex indices. + */ + bool operator == (const Face &face) const; }; diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index 1c8b78fa84..4ac2484c39 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -59,9 +59,9 @@ struct CellData * 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 @@ -75,7 +75,8 @@ struct CellData void rotate (const unsigned int times); /** - * Rotations are not possible in 1d. + * Rotations are not possible in + * 1d. */ DeclException0 (ExcNotPossible); /** diff --git a/deal.II/deal.II/source/grid/grid_reordering.cc b/deal.II/deal.II/source/grid/grid_reordering.cc index 1f39c2c6c5..de8572152d 100644 --- a/deal.II/deal.II/source/grid/grid_reordering.cc +++ b/deal.II/deal.II/source/grid/grid_reordering.cc @@ -13,6 +13,7 @@ #include +#include #include @@ -277,11 +278,165 @@ GridReordering<3>::Cell::insert_faces (map &global_faces) 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 + cell_orientation_faces[rotational_states_of_cells][GeometryInfo::faces_per_cell]; + + if (already_initialized == false) + { + for (unsigned int rot=0; rot standard_cell; + for (unsigned int v=0; v::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::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) + { + unsigned int f, r; + for (f=0; f::faces_per_cell; ++f) + for (r=0; r::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 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; 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]; }; #endif @@ -446,6 +601,17 @@ bool GridReordering::Face::operator < (const Face &face) const +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) diff --git a/deal.II/deal.II/source/grid/tria.all_dimensions.cc b/deal.II/deal.II/source/grid/tria.all_dimensions.cc index b5fcba1623..dfaf7bb6bd 100644 --- a/deal.II/deal.II/source/grid/tria.all_dimensions.cc +++ b/deal.II/deal.II/source/grid/tria.all_dimensions.cc @@ -32,7 +32,7 @@ template <> -void CellData<1>::rotate (unsigned int) +void CellData<1>::rotate (const unsigned int) { Assert (false, ExcNotPossible()); }; @@ -40,26 +40,66 @@ void CellData<1>::rotate (unsigned int) 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 -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::vertices_per_cell; ++i) + vertices[i] = old_vertices[rotations[times][i]]; }; -- 2.39.5