From 881b5b822437b6076c350e814d985a907e249020 Mon Sep 17 00:00:00 2001
From: hartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 21 Nov 2001 10:55:21 +0000
Subject: [PATCH] New CellAccessor<dim>::neighbor_of_coarser_neighbor function.

git-svn-id: https://svn.dealii.org/trunk@5230 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/deal.II/include/grid/tria_accessor.h | 33 ++++++++--
 deal.II/deal.II/source/grid/tria_accessor.cc | 65 ++++++++++++++++++++
 2 files changed, 92 insertions(+), 6 deletions(-)

diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h
index 2027172307..e040cd209e 100644
--- 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 9701c1dd08..11d45037db 100644
--- 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
 {
-- 
2.39.5