]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement the 3d part of grid_reordering.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 5 Oct 2000 08:44:20 +0000 (08:44 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 5 Oct 2000 08:44:20 +0000 (08:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@3375 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_reordering.h
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/source/grid/grid_reordering.cc
deal.II/deal.II/source/grid/tria.all_dimensions.cc

index 415c8969f21dbd2a56e5e58a5414f6e67b672125..08a8121031f8e5698af55034798b6daad38b50a8 100644 (file)
@@ -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<unsigned int>(-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<dim>
                                          * comparison.
                                          */
        bool operator < (const Face &face) const;
+
+                                        /**
+                                         * Check for equality of
+                                         * vertex indices.
+                                         */
+       bool operator == (const Face &face) const;
     };
 
 
index 1c8b78fa84065bcf4ae6f355b7506e9bca55a77c..4ac2484c39299b0cb13d81a56a9e3dae2103dd89 100644 (file)
@@ -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);
                                     /**
index 1f39c2c6c52c0da69cd7401806200b339e933669..de8572152d792997214042484c7e58ea5d1b8dd5 100644 (file)
@@ -13,6 +13,7 @@
 
 
 #include <algorithm>
+#include <base/thread_management.h>
 #include <grid/grid_reordering.h>
 
 
@@ -277,11 +278,165 @@ GridReordering<3>::Cell::insert_faces (map<Face,FaceData> &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<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
@@ -446,6 +601,17 @@ bool GridReordering<dim>::Face::operator < (const Face &face) const
 
 
 
+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)
index b5fcba1623316680cb0f9aec5d5330ad50bb7fc4..dfaf7bb6bdfd3cb49a302934c85bc34c7334942a 100644 (file)
@@ -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<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]];
 };
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.