New CellAccessor<dim>::neighbor_of_coarser_neighbor function.

git-svn-id: https://svn.dealii.org/trunk@5230 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h
index 20271723070a5d6d32d625790c9206e234c349a3..e040cd209e4fd6ebfda54f2fd13a4368b5500f7a 100644 (file)
--- a/deal.II/deal.II/include/grid/tria_accessor.h
+++ b/deal.II/deal.II/include/grid/tria_accessor.h
@@ -18,6 +18,7 @@
 #include <base/exceptions.h>
 #include <grid/tria_iterator_base.h>
 
+
 template <int dim> class Point;
 
 template <int dim> class Triangulation;
@@ -27,6 +28,12 @@ template <int dim, typename Accessor> class TriaActiveIterator;
 class Line;
 class Quad;
 class Hexahedron;
+namespace std
+{
+  template<class T1, class T2>
+  struct pair;
+}
+
 
 
 // note: the file tria_accessor.templates.h is included at the end of
@@ -221,6 +228,10 @@ class TriaAccessor
                                      * Exception
                                      */
     DeclException0 (ExcNeighborIsCoarser);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNeighborIsNotCoarser);
                                     /*@}*/
        
   protected:
@@ -1857,15 +1868,25 @@ class CellAccessor :  public TriaObjectAccessor<dim,dim>
                                      * the present cell
                                      * (i.e. @p{cell->neighbor(neighbor)->level()}
                                      * needs to be equal to
-                                     * @p{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.
+                                     * @p{cell->level()}. Use the
+                                     * @p{neighbor_of_coarser_neighbor}
+                                     * function in that case.
                                      */
     unsigned int neighbor_of_neighbor (const unsigned int neighbor) const;
     
+                                    /**
+                                     * This function is a
+                                     * generalization of the
+                                     * @p{neighbor_of_neighbor}
+                                     * function for the case of a
+                                     * coarser neighbor. It returns a
+                                     * pair of numbers, face_no and
+                                     * subface_no, with the following
+                                     * property:
+                                     * @p{cell->neighbor(neighbor)->face(face_no)->child(subface_no)==cell}.
+                                     */
+    std::pair<unsigned int, unsigned int> neighbor_of_coarser_neighbor (const unsigned int neighbor) const;
+    
                                     /**
                                      *  Return whether the @p{i}th
                                      *  vertex or face (depending on

diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc
index 9701c1dd08e8ddfe14a111fdb5218681306504e1..11d45037db546b9b7a01db787d420643e95eba22 100644 (file)
--- a/deal.II/deal.II/source/grid/tria_accessor.cc
+++ b/deal.II/deal.II/source/grid/tria_accessor.cc
@@ -1806,6 +1806,71 @@ unsigned int CellAccessor<dim>::neighbor_of_neighbor (const unsigned int neighbo
 
 
 
+template <int dim>
+pair<unsigned int, unsigned int>
+CellAccessor<dim>::neighbor_of_coarser_neighbor (const unsigned int neighbor) const
+{
+                                  // make sure that the neighbor is
+                                  // on a coarser level
+  Assert (neighbor_level(neighbor) < present_level,
+         typename TriaAccessor<dim>::ExcNeighborIsNotCoarser());
+  Assert (neighbor < GeometryInfo<dim>::faces_per_cell,
+         typename TriaAccessor<dim>::ExcInvalidNeighbor(neighbor));
+
+  const TriaIterator<dim,TriaObjectAccessor<dim-1, dim> > this_face=face(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 face_no_guess
+    = GeometryInfo<dim>::opposite_face[neighbor];
+
+  const TriaIterator<dim,TriaObjectAccessor<dim-1, dim> > face_guess
+    =neighbor_cell->face(face_no_guess);
+  
+  if (face_guess->has_children())
+    for (unsigned int subface_no=0; subface_no<GeometryInfo<dim>::subfaces_per_face; ++subface_no)
+      if (face_guess->child(subface_no)==this_face)
+       return pair<unsigned int,unsigned int> (face_no_guess, subface_no);
+
+                                    // if the guess was false, then
+                                    // we need to loop over all faces
+                                    // and subfaces and find the
+                                    // number the hard way
+  for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+    {
+      if (face_no!=face_no_guess)
+       {
+         const TriaIterator<dim,TriaObjectAccessor<dim-1, dim> > face
+           =neighbor_cell->face(face_no);
+         if (face->has_children())
+           for (unsigned int subface_no=0; subface_no<GeometryInfo<dim>::subfaces_per_face; ++subface_no)
+             if (face->child(subface_no)==this_face)
+               return pair<unsigned int,unsigned int> (face_no, subface_no);
+       }
+    }
+  
+                                  // we should never get here,
+                                  // since then we did not find
+                                  // our way back...
+  Assert (false, ExcInternalError());
+  return pair<unsigned int,unsigned int> (static_cast<unsigned int>(-1),
+                                         static_cast<unsigned int>(-1));
+};
+
+
+
 template <int dim>
 bool CellAccessor<dim>::at_boundary (const unsigned int i) const
 {