]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement and use CellAccessor::neighbor_of_neighbor
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 2 Feb 2000 12:02:30 +0000 (12:02 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 2 Feb 2000 12:02:30 +0000 (12:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@2325 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/geometry_info.h
deal.II/deal.II/include/grid/tria_accessor.h
deal.II/deal.II/source/grid/geometry_info.cc
deal.II/deal.II/source/grid/tria.cc
deal.II/deal.II/source/grid/tria_accessor.cc
deal.II/deal.II/source/numerics/error_estimator.cc

index eec80241045eb09c0259ad85cfbe504969f127b9..f79fc7993e0d6670f4abf4affcbc167860f87716 100644 (file)
@@ -127,6 +127,15 @@ struct GeometryInfo
     static const unsigned int hexes_per_cell = (2*GeometryInfo<dim-1>::hexes_per_cell
                                                + GeometryInfo<dim-1>::quads_per_cell);
 
+                                    /**
+                                     * List of numbers which is
+                                     * denote which face is opposite
+                                     * to a given face. In 1d, this
+                                     * list is #{1,0}#, in 2d #{2, 3, 0, 1}#,
+                                     * in 3d #{1, 0, 4, 5, 2, 3}#.
+                                     */
+    static const unsigned int opposite_face[faces_per_cell];
+    
                                     /**
                                      * This field store which child cells
                                      * are adjacent to a certain face of
index 6d87bb6143b7b9d06411a68a1da51376cf5c38b5..29a92ffbebe2520217a072262912cb7ab6b0894f 100644 (file)
@@ -207,6 +207,10 @@ class TriaAccessor
                                      * Exception
                                      */
     DeclException0 (ExcCantCompareIterators);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNeighborIsCoarser);
                                     /*@}*/
        
   protected:
@@ -1711,10 +1715,11 @@ class TriaObjectAccessor<3, dim> :  public TriaAccessor<dim>
  * the possibility to check whether they are at the boundary etc. This class
  * offers access to all this data.
  *
- * @author Wolfgang Bangerth, 1998
+ * @author Wolfgang Bangerth, 1998, 1999, 2000
  */
 template <int dim>
-class CellAccessor :  public TriaObjectAccessor<dim,dim> {
+class CellAccessor :  public TriaObjectAccessor<dim,dim>
+{
   public:
                                     /**
                                      * Propagate the AccessorData type
@@ -1761,6 +1766,33 @@ class CellAccessor :  public TriaObjectAccessor<dim,dim> {
     void set_neighbor (const unsigned int i,
                       const TriaIterator<dim,CellAccessor<dim> > &pointer) const;
 
+                                    /**
+                                     * Return the how-many'th
+                                     * neighbor this cell is of
+                                     * #cell->neighbor(neighbor)#,
+                                     * i.e. return the number #n#
+                                     * such that
+                                     * #cell->neighbor(neighbor)->neighbor(n)==cell#. This
+                                     * function is the right one if
+                                     * you want to know how to get
+                                     * back from a neighbor to the
+                                     * present cell.
+                                     *
+                                     * Note that this operation is
+                                     * only useful if the neighbor is
+                                     * not on a coarser level than
+                                     * the present cell
+                                     * (i.e. #cell->neighbor(neighbor)->level()#
+                                     * needs to be equal to
+                                     * #cell->level()#, since
+                                     * otherwise the neighbors of the
+                                     * neighbor cell are on a coarser
+                                     * level than the present one and
+                                     * you can't get back from there
+                                     * to this cell.
+                                     */
+    unsigned int neighbor_of_neighbor (const unsigned int neighbor) const;
+    
                                     /**
                                      *  Return whether the #i#th vertex or
                                      *  face (depending on the dimension) is
index b68cc09f4f56b17cd64901a3733acd39ce71180a..a5cfea3080d7a4a13a2e68a737a38cfca1f3a917 100644 (file)
 //const unsigned int GeometryInfo<deal_II_dimension>::children_per_cell;
 
 
+template <>
+const unsigned int GeometryInfo<1>::opposite_face[GeometryInfo<1>::faces_per_cell]
+= { 0, 1 };
+
+
+template <>
+const unsigned int GeometryInfo<2>::opposite_face[GeometryInfo<2>::faces_per_cell]
+= { 2, 3, 0, 1 };
+
+
+template <>
+const unsigned int GeometryInfo<3>::opposite_face[GeometryInfo<3>::faces_per_cell]
+= { 1, 0, 4, 5, 2, 3 };
 
index 4834e6bbcd2d6c839be9f89b95c838b6496a47ef..dae91a404f5a7afdc31ce5fe0ec0815cc68c5f56 100644 (file)
@@ -4242,43 +4242,41 @@ void Triangulation<2>::execute_refinement () {
                                                     // two children which
                                                     // we can use.
                    cell_iterator neighbor = cell->neighbor(nb);
-                   for (unsigned int nb_nb=0; nb_nb<4; ++nb_nb)
-                     if (neighbor->neighbor(nb_nb)==cell)
                                                         // this cell is the nb_nb-th
                                                         // neighbor or neighbor(nb)
-                       {
-                         neighbors_neighbor[2*nb] = neighbors_neighbor[2*nb+1] = nb_nb;
-                                                          // vertex 1 of child 0
-                                                          // is always the interior
-                                                          // one
-                         new_vertices[2*nb+1] = neighbor->line(nb_nb)
-                                                ->child(0)->vertex_index(1);
-
-                         if (nb < 2) 
-                           {
-                             new_lines[2*nb]  = neighbor->line(nb_nb)->child(0);
-                             new_lines[2*nb+1]= neighbor->line(nb_nb)->child(1);
-                           } else {
-                                                              // lines 2 and 3 have
-                                                              // opposite sense
-                             new_lines[2*nb]  = neighbor->line(nb_nb)->child(1);
-                             new_lines[2*nb+1]= neighbor->line(nb_nb)->child(0);
-                           };
-
-                                                          // finally find out which
-                                                          // are the two neighbor
-                                                          // subcells, adjacent to
-                                                          // the two sublines
-                         const int child_mapping[4][2] = {{0,1},{1,2},{3,2},{0,3}};     
-                         if (nb < 2) 
-                           {
-                             neighbors[2*nb]  = neighbor->child(child_mapping[nb_nb][0]);
-                             neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][1]);
-                           } else {
-                             neighbors[2*nb]  = neighbor->child(child_mapping[nb_nb][1]);
-                             neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][0]);
-                           };
-                       };
+                   const unsigned int nb_nb = cell->neighbor_of_neighbor (nb);
+
+                   neighbors_neighbor[2*nb] = neighbors_neighbor[2*nb+1] = nb_nb;
+                                                    // vertex 1 of child 0
+                                                    // is always the interior
+                                                    // one
+                   new_vertices[2*nb+1] = neighbor->line(nb_nb)
+                                          ->child(0)->vertex_index(1);
+
+                   if (nb < 2) 
+                     {
+                       new_lines[2*nb]  = neighbor->line(nb_nb)->child(0);
+                       new_lines[2*nb+1]= neighbor->line(nb_nb)->child(1);
+                     } else {
+                                                        // lines 2 and 3 have
+                                                        // opposite sense
+                       new_lines[2*nb]  = neighbor->line(nb_nb)->child(1);
+                       new_lines[2*nb+1]= neighbor->line(nb_nb)->child(0);
+                     };
+                   
+                                                    // finally find out which
+                                                    // are the two neighbor
+                                                    // subcells, adjacent to
+                                                    // the two sublines
+                   static const unsigned int child_mapping[4][2] = {{0,1},{1,2},{3,2},{0,3}};     
+                   if (nb < 2) 
+                     {
+                       neighbors[2*nb]  = neighbor->child(child_mapping[nb_nb][0]);
+                       neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][1]);
+                     } else {
+                       neighbors[2*nb]  = neighbor->child(child_mapping[nb_nb][1]);
+                       neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][0]);
+                     };
                  }
                else 
            
@@ -5560,7 +5558,7 @@ void Triangulation<3>::execute_refinement () {
            for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell;
                 ++face)
              {
-               cell_iterator neighbor = hex->neighbor(face);
+               const cell_iterator neighbor = hex->neighbor(face);
 
                                                 // if no neighbor
                if (neighbor.state() != valid)
@@ -5604,11 +5602,7 @@ void Triangulation<3>::execute_refinement () {
                                                         // of the neighbor
                                                         // adjacent to which
                                                         // the present cell is
-                       unsigned int nb_nb;
-                       for (nb_nb=0; nb_nb<GeometryInfo<dim>::faces_per_cell;
-                            ++nb_nb)
-                         if (neighbor->neighbor(nb_nb) == hex)
-                           break;
+                       const unsigned int nb_nb = hex->neighbor_of_neighbor(face);
                        Assert (nb_nb<GeometryInfo<dim>::faces_per_cell,
                                ExcInternalError());
 
index 301e7093b46fd23431b4ca451d433e664c76f3f3..ae058d1562e3f7a112d81fbbb4d6c1485c75e8ed 100644 (file)
@@ -1599,6 +1599,58 @@ void CellAccessor<dim>::set_neighbor (const unsigned int i,
 
 
 
+template <int dim>
+unsigned int CellAccessor<dim>::neighbor_of_neighbor (const unsigned int neighbor) const
+{
+                                  // make sure that the neighbor is
+                                  // not on a coarser level
+  Assert (neighbor_level(neighbor) == present_level,
+         ExcNeighborIsCoarser());
+  Assert (neighbor < GeometryInfo<dim>::faces_per_cell,
+         ExcInvalidNeighbor(neighbor));
+
+  const TriaIterator<dim,CellAccessor<dim> > neighbor_cell = this->neighbor(neighbor);
+  
+                                  // usually, on regular patches of
+                                  // the grid, this cell is just on
+                                  // the opposite side of the
+                                  // neighbor that the neighbor is of
+                                  // this cell. for example in 2d, if
+                                  // we want to know the
+                                  // neighbor_of_neighbor if
+                                  // neighbor==1 (the right
+                                  // neighbor), then we will get 3
+                                  // (the left neighbor) in most
+                                  // cases. look up this relationship
+                                  // in the table provided by
+                                  // GeometryInfo and try it
+  const unsigned int neighbor_guess
+    = GeometryInfo<dim>::opposite_face[neighbor];
+  
+  if ((neighbor_cell->neighbor_index (neighbor_guess) == present_index) &&
+      (neighbor_cell->neighbor_level (neighbor_guess) == present_level))
+    return neighbor_guess;
+  else
+                                    // if the guess was false, then
+                                    // we need to loop over all
+                                    // neighbors and find the number
+                                    // the hard way
+    {
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+       if ((neighbor_cell->neighbor_index (face) == present_index) &&
+           (neighbor_cell->neighbor_level (face) == present_level))
+         return face;
+
+                                      // we should never get here,
+                                      // since then we did not find
+                                      // our way back...
+      Assert (false, ExcInternalError());
+      return static_cast<unsigned int>(-1);
+    };
+};
+
+
+
 template <int dim>
 bool CellAccessor<dim>::at_boundary (const unsigned int i) const {
   Assert (used(), ExcCellNotUsed());
index f4ca76d3636b1b57ce4659171e825653942fcd85..7f08d1532fcecba98195e3a4b59ec7fba6ac3f9f 100644 (file)
@@ -578,18 +578,14 @@ integrate_over_regular_face (Data &data,
                                     // of gradient across this face
     {
       Assert (cell->neighbor(face_no).state() == valid,
-             ExcInternalError());
-      unsigned int neighbor_neighbor;
-      DoFHandler<dim>::active_cell_iterator neighbor = cell->neighbor(face_no);
+             ExcInternalError());      
+      
+      const DoFHandler<dim>::active_cell_iterator neighbor = cell->neighbor(face_no);
       
                                       // find which number the current
                                       // face has relative to the neighboring
                                       // cell
-      for (neighbor_neighbor=0; neighbor_neighbor<GeometryInfo<dim>::faces_per_cell;
-          ++neighbor_neighbor)
-       if (neighbor->neighbor(neighbor_neighbor) == cell)
-         break;
-      
+      const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
       Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell, ExcInternalError());
       
                                       // get restriction of finite element
@@ -720,6 +716,7 @@ integrate_over_irregular_face (Data &data,
                               FESubfaceValues<dim>       &fe_subface_values)
 {
   const DoFHandler<dim>::cell_iterator neighbor = cell->neighbor(face_no);
+
   Assert (neighbor.state() == valid, ExcInternalError());
   Assert (neighbor->has_children(), ExcInternalError());
                                   // set up a vector of the gradients
@@ -736,11 +733,7 @@ integrate_over_irregular_face (Data &data,
   
                                   // store which number #cell# has in the
                                   // list of neighbors of #neighbor#
-  unsigned int neighbor_neighbor;
-  for (neighbor_neighbor=0; neighbor_neighbor<GeometryInfo<dim>::faces_per_cell;
-       ++neighbor_neighbor)
-    if (neighbor->neighbor(neighbor_neighbor) == cell)
-      break;
+  const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
   Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell, ExcInternalError());
   
                                   // loop over all subfaces

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.