Fix CellAccessor<dim>::neighbor_of_neighbor() function in case of cells neighboring...

git-svn-id: https://svn.dealii.org/trunk@14964 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 89f3505212a2c02cc9a65b0c5471aa41f335744c..7680f26837641a36b7d9c032f0d7a864d82a086a 100644 (file)
--- a/deal.II/deal.II/include/grid/tria_accessor.h
+++ b/deal.II/deal.II/include/grid/tria_accessor.h
@@ -3348,6 +3348,18 @@ class CellAccessor :  public TriaObjectAccessor<dim,dim>
     TriaIterator<dim,TriaObjectAccessor<dim-1, dim> >
     face (const unsigned int i) const;
 
+                                    /**
+                                     * Return the (global) index of the
+                                     * @p ith face of this cell.
+                                     *
+                                     * This function is not
+                                     * implemented in 1D, and maps to
+                                     * line_index in 2D and quad_index
+                                     * in 3D.
+                                     */
+    unsigned int
+    face_index (const unsigned int i) const;
+
                                      /**
                                       * Return an iterator to that
                                       * cell that neighbors the

diff --git a/deal.II/deal.II/include/grid/tria_accessor.templates.h b/deal.II/deal.II/include/grid/tria_accessor.templates.h
index 2e184ee1d7a2b00735e79f424d0444ac6f41588d..da7eb4ab1fcdeeb8cf7101e2341dca86ff91b09d 100644 (file)
--- a/deal.II/deal.II/include/grid/tria_accessor.templates.h
+++ b/deal.II/deal.II/include/grid/tria_accessor.templates.h
@@ -1301,6 +1301,37 @@ CellAccessor<3>::face (const unsigned int i) const
 
 
 
+template <>
+inline
+unsigned int
+CellAccessor<1>::face_index (const unsigned int) const 
+{
+  Assert (false, ExcImpossibleInDim(1));
+  return deal_II_numbers::invalid_unsigned_int;
+}
+
+
+
+template <>
+inline
+unsigned int
+CellAccessor<2>::face_index (const unsigned int i) const 
+{
+  return this->line_index(i);
+}
+
+
+
+template <>
+inline
+unsigned int
+CellAccessor<3>::face_index (const unsigned int i) const 
+{
+  return this->quad_index(i);
+}
+
+
+
 template <int dim>
 inline
 int

diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc
index 26efef1a272a6289099dff732f50ea514dfeda59..49a1f469827ec562623f4b93d033497ec0ce21f6 100644 (file)
--- a/deal.II/deal.II/source/grid/tria_accessor.cc
+++ b/deal.II/deal.II/source/grid/tria_accessor.cc
@@ -2156,6 +2156,9 @@ unsigned int CellAccessor<dim>::neighbor_of_neighbor (const unsigned int neighbo
   Assert (neighbor < GeometryInfo<dim>::faces_per_cell,
          TriaAccessorExceptions::ExcInvalidNeighbor(neighbor));
 
+  if (dim==1)
+    return GeometryInfo<dim>::opposite_face[neighbor];
+  
   const TriaIterator<dim,CellAccessor<dim> > neighbor_cell = this->neighbor(neighbor);
   
                                   // usually, on regular patches of
@@ -2171,12 +2174,14 @@ unsigned int CellAccessor<dim>::neighbor_of_neighbor (const unsigned int neighbo
                                   // cases. look up this relationship
                                   // in the table provided by
                                   // GeometryInfo and try it
+
+  const unsigned int this_face_index=face_index(neighbor);
+  
   const unsigned int neighbor_guess
     = GeometryInfo<dim>::opposite_face[neighbor];
   
-  if ((neighbor_cell->neighbor_index (neighbor_guess) == this->present_index) &&
-      (neighbor_cell->neighbor_level (neighbor_guess) == this->present_level))
-    return neighbor_guess;
+  if (neighbor_cell->face_index (neighbor_guess) == this_face_index)
+      return neighbor_guess;
   else
                                     // if the guess was false, then
                                     // we need to loop over all
@@ -2184,8 +2189,7 @@ unsigned int CellAccessor<dim>::neighbor_of_neighbor (const unsigned int neighbo
                                     // the hard way
     {
       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-       if ((neighbor_cell->neighbor_index (face) == this->present_index) &&
-           (neighbor_cell->neighbor_level (face) == this->present_level))
+       if (neighbor_cell->face_index (face) == this_face_index)
          return face;
 
                                       // we should never get here,
@@ -2193,7 +2197,7 @@ unsigned int CellAccessor<dim>::neighbor_of_neighbor (const unsigned int neighbo
                                       // our way back...
       Assert (false, ExcInternalError());
       return deal_II_numbers::invalid_unsigned_int;
-    };
+    }
 }
 
 
@@ -2208,8 +2212,9 @@ CellAccessor<dim>::neighbor_of_coarser_neighbor (const unsigned int neighbor) co
          TriaAccessorExceptions::ExcNeighborIsNotCoarser());
   Assert (neighbor < GeometryInfo<dim>::faces_per_cell,
          TriaAccessorExceptions::ExcInvalidNeighbor(neighbor));
-
-  const TriaIterator<dim,TriaObjectAccessor<dim-1, dim> > this_face=face(neighbor);
+  Assert (dim>1, ExcImpossibleInDim(dim));
+  
+  const int this_face_index=face_index(neighbor);
   const TriaIterator<dim,CellAccessor<dim> > neighbor_cell = this->neighbor(neighbor);
   
                                   // usually, on regular patches of
@@ -2233,7 +2238,7 @@ CellAccessor<dim>::neighbor_of_coarser_neighbor (const unsigned int neighbor) co
   
   if (face_guess->has_children())
     for (unsigned int subface_no=0; subface_no<face_guess->n_children(); ++subface_no)
-      if (face_guess->child(subface_no)==this_face)
+      if (face_guess->child_index(subface_no)==this_face_index)
        return std::make_pair (face_no_guess, subface_no);
 
                                     // if the guess was false, then
@@ -2248,7 +2253,7 @@ CellAccessor<dim>::neighbor_of_coarser_neighbor (const unsigned int neighbor) co
            =neighbor_cell->face(face_no);
          if (face->has_children())
            for (unsigned int subface_no=0; subface_no<face->n_children(); ++subface_no)
-             if (face->child(subface_no)==this_face)
+             if (face->child_index(subface_no)==this_face_index)
                return std::make_pair (face_no, subface_no);
        }
     }