diff --git a/grid/uggrid/ugintersectionit.cc b/grid/uggrid/ugintersectionit.cc index 554d114d787e40838b4c9a29c405f0c387566159..a517ae700ecb0e444c8331b65c82ec9a253163ea 100644 --- a/grid/uggrid/ugintersectionit.cc +++ b/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 { diff --git a/grid/uggrid/ugintersectionit.hh b/grid/uggrid/ugintersectionit.hh index 434b2857c08631a77c6a3712529185901b4ab85b..d3ddc4220c29ce5cf04581d6494d7bb07aba9381 100644 --- a/grid/uggrid/ugintersectionit.hh +++ b/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"