UGGrid has the new intersection iterator functionality now.

leafNeighbor() and levelNeighbor() are the default implementations
of Robert based on the neighbor() method.

[[Imported from SVN: r4488]]

grid/uggrid/ugintersectionit.cc

@@ -3,25 +3,73 @@
 template<class GridImp>
 inline typename UGTypes<GridImp::dimension>::Element* UGGridIntersectionIterator< GridImp >::getNeighbor () const
 {
-  // if we have a neighbor on this level, then return it
-  if (UG_NS<dim>::NbElem(center_, neighborCount_)!=NULL)
-    return UG_NS<dim>::NbElem(center_, neighborCount_);
+  // if subcount is zero and there is a level neighbor return it
+  if (subCount_==0 && getLevelNeighbor()!=NULL)
+    return getLevelNeighbor();
+
+  // now we are either in subcount 1 or there is no level neighbor. If we are a leaf return leaf neighbor
+  if (UG_NS<dim>::isLeaf(center_))
+    return getLeafNeighbor();
+
+  // return nothing
+  return NULL;
+}
 
-  // now go down the stack of copies to find a lower level leaf neighbor
-  typename UGTypes<dim>::Element* father_ = UG_NS<GridImp::dimensionworld>::EFather(center_);
-  while (father_!=0)
+// returns a neighbor that is a leaf or nothing (neighbor might be on the same level)
+// works only on leaf elements!
+template<class GridImp>
+inline typename UGTypes<GridImp::dimension>::Element* UGGridIntersectionIterator< GridImp >::getLeafNeighbor () const
+{
+  // if the level neighbor exists and is a leaf then return it
+  typename UGTypes<dim>::Element* p = UG_NS<dim>::NbElem(center_, neighborCount_);
+  if (p!=NULL)
+    if (UG_NS<dim>::isLeaf(p))
+      return p;
+
+  // now I must be a leaf to proceed
+  if (!UG_NS<dim>::isLeaf(center_))
+    return NULL;
+
+  // up or down ?
+  if (p==NULL)
+  {
+    // I am a leaf and the neighbor does not exist: go down
+    typename UGTypes<dim>::Element* father_ = UG_NS<GridImp::dimensionworld>::EFather(center_);
+    while (father_!=0)
+    {
+      if (!UG_NS<dim>::hasCopy(father_)) break;             // father must be a copy
+      if (UG_NS<dim>::NbElem(father_, neighborCount_)!=NULL)             // check existence of neighbor
+        if (UG_NS<dim>::isLeaf(UG_NS<dim>::NbElem(father_, neighborCount_)))                 // check leafness
+          return UG_NS<dim>::NbElem(father_, neighborCount_);
+      father_ = UG_NS<dim>::EFather(father_);
+    }
+  }
+  else
   {
-    if (!UG_NS<dim>::hasCopy(father_)) break;       // father must be a copy
-    if (UG_NS<dim>::NbElem(father_, neighborCount_)!=NULL)       // check existence of neighbor
-      if (UG_NS<dim>::isLeaf(UG_NS<dim>::NbElem(father_, neighborCount_)))           // check leafness
-        return UG_NS<dim>::NbElem(father_, neighborCount_);
-    father_ = UG_NS<dim>::EFather(father_);
+    // I am a leaf and the neighbor exists and the neighbor is not a leaf: go up
+    while (p!=0)
+    {
+      if (!UG_NS<dim>::hasCopy(p)) break;             // element must be copy refined
+      typename UGTypes<dim>::Element *sons[32];
+      UG_NS<dim>::GetSons(p,sons);
+      p = sons[0];
+      if (UG_NS<dim>::isLeaf(p))
+        return p;
+    }
   }
 
   // nothing found, return 0 (might be a processor boundary)
   return NULL;
 }
 
+// return a neighbor that is on the same level or nothing (neighbor might be a leaf)
+template<class GridImp>
+inline typename UGTypes<GridImp::dimension>::Element* UGGridIntersectionIterator< GridImp >::getLevelNeighbor () const
+{
+  // return level neighbor or NULL
+  return UG_NS<dim>::NbElem(center_, neighborCount_);
+}
+
 template<class GridImp>
 inline bool UGGridIntersectionIterator< GridImp >::neighbor() const
 {

grid/uggrid/ugintersectionit.hh

@@ -46,7 +46,7 @@ namespace Dune {
         \todo Should be private
      */
     UGGridIntersectionIterator(typename TargetType<0,GridImp::dimensionworld>::T* center, int nb, int level)
-      : center_(center), level_(level), neighborCount_(nb)
+      : center_(center), level_(level), neighborCount_(nb), subCount_(0)
     {}
 
     //! The Destructor
@@ -58,12 +58,23 @@ namespace Dune {
 
     //! equality
     bool equals(const UGGridIntersectionIterator<GridImp>& i) const {
-      return center_==i.center_ && neighborCount_ == i.neighborCount_;
+      return center_==i.center_ && neighborCount_ == i.neighborCount_ && subCount_==i.subCount_;
     }
 
     //! prefix increment
     void increment() {
-      neighborCount_++;
+      typename UGTypes<dim>::Element* p = getLevelNeighbor();
+      bool secondnb=false;
+      if (subCount_==0 && p!=0)
+        if (UG_NS<dim>::isLeaf(p)==false && getLeafNeighbor()!=NULL)
+          secondnb=true;
+      if (secondnb)
+        subCount_++;
+      else
+      {
+        neighborCount_++;
+        subCount_=0;
+      }
       if (neighborCount_ >= UG_NS<GridImp::dimensionworld>::Sides_Of_Elem(center_))
         neighborCount_ = -1;
     }
@@ -131,6 +142,12 @@ namespace Dune {
     //! get neighbor on same or lower level or 0
     typename UGTypes<GridImp::dimension>::Element* getNeighbor () const;
 
+    //! returns a neighbor that is a leaf or nothing (neighbor might be on the same level)
+    typename UGTypes<GridImp::dimension>::Element* getLeafNeighbor () const;
+
+    //! return a neighbor that is on the same level or nothing (neighbor might be a leaf)
+    typename UGTypes<GridImp::dimension>::Element* getLevelNeighbor () const;
+
     //! vector storing the outer normal
     mutable FieldVector<UGCtype, dimworld> outerNormal_;
 
@@ -152,7 +169,8 @@ namespace Dune {
     //! count on which neighbor we are lookin' at. Note that this is interpreted in UG's ordering!
     int neighborCount_;
 
-
+    // count number of neighbors on current face that have been treated
+    int subCount_;
   };
 
 #include "ugintersectionit.cc"