• Home
• About DUNE
  • What is DUNE
  • News
  • Gallery
  • Publications
• Download
  • Core Modules (Stable)
    • Release Notes (v1.1)
    • Release Notes (v1.0)
  • Recent Changes
  • External Modules
  • Dune Logo
• Documentation
  • Installation
  • Using DUNE
  • Class Documentation
• Community
  • People
  • Mailinglists
  • FAQ
• Development
  • Developing DUNE
  • Unstable Modules
    • Download via SVN
    • Class Documentation
  • Bug Tracking System
  • Automated Test Environment

yaspgrid.hh

Go to the documentation of this file.

#ifndef DUNE_YASPGRID_HH
#define DUNE_YASPGRID_HH

#include<iostream>
#include<vector>
#include<map>
#include<algorithm>
#include<stack>

#include <dune/grid/common/grid.hh>     // the grid base classes
#include <dune/grid/yaspgrid/grids.hh>  // the yaspgrid base classes
#include <dune/grid/common/capabilities.hh> // the capabilities
#include <dune/common/misc.hh>
#include <dune/common/helpertemplates.hh>
#include <dune/common/bigunsignedint.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/collectivecommunication.hh>
#include <dune/grid/common/indexidset.hh>
#include <dune/grid/common/datahandleif.hh>


#if HAVE_MPI
#include <dune/common/mpicollectivecommunication.hh>
#endif

namespace Dune {

//************************************************************************
typedef double yaspgrid_ctype;
static const yaspgrid_ctype yasptolerance=1E-13; // tolerance in coordinate computations

/* some sizes for building global ids
 */
  const int yaspgrid_dim_bits = 24;   // bits for encoding each dimension
  const int yaspgrid_level_bits = 6;  // bits for encoding level number
  const int yaspgrid_codim_bits = 4;  // bits for encoding codimension


//************************************************************************
// forward declaration of templates

template<int dim, int dimworld>               class YaspGrid;
template<int mydim, int cdim, class GridImp>  class YaspGeometry;
template<int codim, int dim, class GridImp>   class YaspEntity;
template<int codim, class GridImp>            class YaspEntityPointer;
template<int codim, PartitionIteratorType pitype, class GridImp> class YaspLevelIterator;
template<class GridImp>            class YaspIntersectionIterator;
template<class GridImp>            class YaspHierarchicIterator;
template<class GridImp>            class YaspLevelIndexSet;
template<class GridImp>            class YaspLeafIndexSet;
template<class GridImp>            class YaspGlobalIdSet;

//========================================================================
//========================================================================

template<int mydim, int cdim, class GridImp>
class YaspSpecialGeometry : public Geometry<mydim, cdim, GridImp, YaspGeometry>
{
  typedef typename GridImp::ctype ctype;
public:
  YaspSpecialGeometry(const FieldVector<ctype, cdim>& p, const FieldVector<ctype, cdim>& h, int& m) :
    Geometry<mydim, cdim, GridImp, YaspGeometry>(YaspGeometry<mydim, cdim, GridImp>(p,h,m))
    {}
  YaspSpecialGeometry() :
    Geometry<mydim, cdim, GridImp, YaspGeometry>(YaspGeometry<mydim, cdim, GridImp>(false))
    {};
};

template<int mydim, class GridImp>
class YaspSpecialGeometry<mydim,mydim,GridImp> : public Geometry<mydim, mydim, GridImp, YaspGeometry>
{
  typedef typename GridImp::ctype ctype;
public:
  YaspSpecialGeometry(const FieldVector<ctype, mydim>& p, const FieldVector<ctype, mydim>& h) :
    Geometry<mydim, mydim, GridImp, YaspGeometry>(YaspGeometry<mydim, mydim, GridImp>(p,h))
    {}
  YaspSpecialGeometry() :
    Geometry<mydim, mydim, GridImp, YaspGeometry>(YaspGeometry<mydim, mydim, GridImp>(false))
    {};
};

template<int cdim, class GridImp>
class YaspSpecialGeometry<0,cdim,GridImp> : public Geometry<0, cdim, GridImp, YaspGeometry>
{
  typedef typename GridImp::ctype ctype;
public:
  YaspSpecialGeometry(const FieldVector<ctype, cdim>& p) :
    Geometry<0, cdim, GridImp, YaspGeometry>(YaspGeometry<0, cdim, GridImp>(p))
    {}
  YaspSpecialGeometry(const FieldVector<ctype, cdim>& p, const FieldVector<ctype, cdim>& h, int& m) :
    Geometry<0, cdim, GridImp, YaspGeometry>(YaspGeometry<0, cdim, GridImp>(p))
    {}
  YaspSpecialGeometry() :
    Geometry<0, cdim, GridImp, YaspGeometry>(YaspGeometry<0, cdim, GridImp>(false))
    {};
};

//========================================================================
// The transformation describing the refinement rule

template<int dim, class GridImp>
class YaspFatherRelativeLocalElement {
public:
  static FieldVector<yaspgrid_ctype, dim> midpoint;  // data neded for the refelem below
  static FieldVector<yaspgrid_ctype, dim> extension; // data needed for the refelem below
  static YaspSpecialGeometry<dim,dim,GridImp> geo;
  static YaspSpecialGeometry<dim,dim,GridImp>& getson (int i)
  {
    for (int k=0; k<dim; k++)
      if (i&(1<<k))
        midpoint[k] = 0.75;
      else
        midpoint[k] = 0.25;
    return geo;
  }
};

// initialize static variable with bool constructor (which makes reference elements)
template<int dim, class GridImp>
YaspSpecialGeometry<dim,dim,GridImp> 
YaspFatherRelativeLocalElement<dim,GridImp>::geo(YaspFatherRelativeLocalElement<dim,GridImp>::midpoint,
                                                 YaspFatherRelativeLocalElement<dim,GridImp>::extension);
template<int dim, class GridImp>
FieldVector<yaspgrid_ctype,dim> YaspFatherRelativeLocalElement<dim,GridImp>::midpoint(0.25);

template<int dim, class GridImp>
FieldVector<yaspgrid_ctype,dim> YaspFatherRelativeLocalElement<dim,GridImp>::extension(0.5);

template<int mydim,int cdim, class GridImp>
class YaspGeometry : public GeometryDefaultImplementation<mydim,cdim,GridImp,YaspGeometry>
{
public:
  typedef typename GridImp::ctype ctype;
  typedef Geometry<mydim,mydim,GridImp,Dune::YaspGeometry> ReferenceGeometry;
  
  GeometryType type () const
  {
      return GeometryType(GeometryType::cube,mydim);
  }

  int corners () const
  {
        return 1<<mydim;
  }

  const FieldVector<ctype, cdim>& operator[] (int i) const
  {
        assert( i >= 0 && i < (int) coord_.N() ); 
        FieldVector<ctype, cdim>& c = coord_[i];
        int bit=0;
        for (int k=0; k<cdim; k++) // run over all directions in world
          {
                if (k==missing)
                  {
                        c[k] = midpoint[k];
                        continue;
                  }
                //k is not the missing direction
                if (i&(1<<bit))  // check whether bit is set or not
                  c[k] = midpoint[k]+0.5*extension[k]; // bit is 1 in i
                else
                  c[k] = midpoint[k]-0.5*extension[k]; // bit is 0 in i
                bit++; // we have processed a direction
          }

        return c;
  }

  FieldVector<ctype, cdim> global (const FieldVector<ctype, mydim>& local) const
  {
      FieldVector<ctype, cdim> g;
      int bit=0;
      for (int k=0; k<cdim; k++)
          if (k==missing)
              g[k] = midpoint[k];
          else
          {
              g[k] = midpoint[k] + (local[bit]-0.5)*extension[k];
              bit++;
          }
      return g;
  }

  FieldVector<ctype, mydim> local (const FieldVector<ctype, cdim>& global) const
  {
      FieldVector<ctype, mydim> l; // result
      int bit=0;
      for (int k=0; k<cdim; k++)
          if (k!=missing)
          {
              l[bit] = (global[k]-midpoint[k])/extension[k] + 0.5;
              bit++;
          }
      return l;
  }

  ctype volume () const
  {
    ctype volume=1.0;
    for (int k=0; k<cdim; k++)
      if (k!=missing) volume *= extension[k];
    return volume;
  }

  ctype integrationElement (const FieldVector<ctype, mydim>& local) const
  {
    return volume();
  }

  FieldMatrix<ctype,cdim,mydim>& jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
  {
      Jinv = 0.0;
      int k=0;
      for (int i=0; i<cdim; ++i)
      {
          if (i != missing)
          {
              Jinv[i][k] = 1.0/extension[i]; // set diagonal element
              k++;
          }
      }
      return Jinv;
  }

  bool checkInside (const FieldVector<ctype, mydim>& local) const
  {
        for (int i=0; i<mydim; i++)
          if (local[i]<-yasptolerance || local[i]>1+yasptolerance) return false;
        return true;
  }

  YaspGeometry (const FieldVector<ctype, cdim>& p, const FieldVector<ctype, cdim>& h, int& m)
        : midpoint(p), extension(h), missing(m)
  {
        if (cdim!=mydim+1)
          DUNE_THROW(GridError, "general YaspGeometry assumes cdim=mydim+1");
  }

  YaspGeometry (const YaspGeometry& other) 
        : midpoint(other.midpoint), 
          extension(other.extension), 
          missing(other.missing)
  {
  }

  void print (std::ostream& s) const
  {
        s << "YaspGeometry<"<<mydim<<","<<cdim<< "> ";
        s << "midpoint";
        for (int i=0; i<cdim; i++)
                  s << " " << midpoint[i];
                s << " extension";
        for (int i=0; i<cdim; i++)
                  s << " " << extension[i];
                s << " missing is " << missing;
  }
private:
  // the element is fully defined by its midpoint the extension
  // in each direction and the missing direction.
  // References are used because this information
  // is known outside the element in many cases.
  // Note cdim==cdim+1

  // IMPORTANT midpoint and extension are references,
  // YaspGeometry can't be copied
  const FieldVector<ctype, cdim> & midpoint;    // the midpoint
  const FieldVector<ctype, cdim> & extension;   // the extension
  int& missing;                         // the missing, i.e. constant direction

  // In addition we need memory in order to return references.
  // Possibly we should change this in the interface ...
  mutable FieldMatrix<ctype, cdim, mydim> Jinv;   // the jacobian inverse
  mutable FieldMatrix<ctype,  Power_m_p<2,mydim>::power, cdim> coord_;   // the coordinates 

  const YaspGeometry<mydim,cdim,GridImp>&
  operator = (const YaspGeometry<mydim,cdim,GridImp>& g);
  
};



template<int mydim, class GridImp>
class YaspGeometry<mydim,mydim,GridImp> : public GeometryDefaultImplementation<mydim,mydim,GridImp,YaspGeometry>
{
public:
  typedef typename GridImp::ctype ctype;
  typedef Geometry<mydim,mydim,GridImp,Dune::YaspGeometry> ReferenceGeometry;
  
  GeometryType type () const
  {
      return GeometryType(GeometryType::cube,mydim);
  }

  int corners () const
  {
        return 1<<mydim;
  }

  const FieldVector<ctype, mydim>& operator[] (int i) const
  {
        assert( i >= 0 && i < (int) coord_.N() ); 
        FieldVector<ctype, mydim>& c = coord_[i];
        for (int k=0; k<mydim; k++)
          if (i&(1<<k))
                c[k] = midpoint[k]+0.5*extension[k]; // kth bit is 1 in i
          else
                c[k] = midpoint[k]-0.5*extension[k]; // kth bit is 0 in i
        return c;
  }

  FieldVector<ctype, mydim> global (const FieldVector<ctype, mydim>& local) const
  {
      FieldVector<ctype,mydim> g;
        for (int k=0; k<mydim; k++)
          g[k] = midpoint[k] + (local[k]-0.5)*extension[k];
        return g;
  }

  FieldVector<ctype, mydim> local (const FieldVector<ctype,mydim>& global) const
  {
      FieldVector<ctype, mydim> l; // result
        for (int k=0; k<mydim; k++)
          l[k] = (global[k]-midpoint[k])/extension[k] + 0.5;
        return l;
  }

  ctype integrationElement (const FieldVector<ctype, mydim>& local) const
  {
    return volume();
  }

  ctype volume () const
  {
    ctype vol=1.0;
    for (int k=0; k<mydim; k++) vol *= extension[k];
    return vol;
  }

  FieldMatrix<ctype,mydim,mydim>& jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
  {
        for (int i=0; i<mydim; ++i)
          {
                Jinv[i] = 0.0;                // set column to zero
                Jinv[i][i] = 1.0/extension[i]; // set diagonal element
          }
        return Jinv;
  }

  bool checkInside (const FieldVector<ctype, mydim>& local) const
  {
        for (int i=0; i<mydim; i++)
          if (local[i]<-yasptolerance || local[i]>1+yasptolerance) return false;
        return true;
  }

  YaspGeometry (const FieldVector<ctype, mydim>& p, const FieldVector<ctype, mydim>& h)
        : midpoint(p), extension(h)
  {}

  YaspGeometry (const YaspGeometry& other) 
        : midpoint(other.midpoint), 
          extension(other.extension) 
  {
  }

  void print (std::ostream& s) const
  {
        s << "YaspGeometry<"<<mydim<<","<<mydim<< "> ";
        s << "midpoint";
        for (int i=0; i<mydim; i++)
                  s << " " << midpoint[i];
                s << " extension";
        for (int i=0; i<mydim; i++)
                  s << " " << extension[i];
  }

private:
  // the element is fully defined by midpoint and the extension
  // in each direction. References are used because this information
  // is known outside the element in many cases.
  // Note mydim==cdim

  // IMPORTANT midpoint and extension are references,
  // YaspGeometry can't be copied
  const FieldVector<ctype, mydim> & midpoint;   // the midpoint
  const FieldVector<ctype, mydim> & extension;  // the extension

  // In addition we need memory in order to return references.
  // Possibly we should change this in the interface ...
  mutable FieldMatrix<ctype, mydim, mydim> Jinv;   // the jacobian inverse
  mutable FieldMatrix<ctype, Power_m_p<2,mydim>::power, mydim> coord_;   // the coordinates 

  // disable copy
  const YaspGeometry<mydim,mydim,GridImp>&
  operator = (const YaspGeometry<mydim,mydim,GridImp>& g);
};

template<int cdim, class GridImp>
class YaspGeometry<0,cdim,GridImp> : public GeometryDefaultImplementation<0,cdim,GridImp,YaspGeometry>
{
public:
  typedef typename GridImp::ctype ctype;
  
  GeometryType type () const
  {
      return GeometryType(GeometryType::cube,0);
  }

  int corners () const
  {
        return 1;
  }

  const FieldVector<ctype, cdim>& operator[] (int i) const
  {
        return position;
  }

  ctype integrationElement (const FieldVector<ctype, 0>& local) const
  {
        return 1.0;
  }

  FieldMatrix<ctype,cdim,0>& jacobianInverseTransposed (const FieldVector<ctype, 0>& local) const
  {
      static FieldMatrix<ctype,cdim,0> Jinv(0.0);
      return Jinv;
  }

  YaspGeometry (const FieldVector<ctype, cdim>& p) : position(p)
  {}

  void print (std::ostream& s) const
  {
        s << "YaspGeometry<"<<0<<","<<cdim<< "> ";
        s << "position " << position;
  }

private:
  // IMPORTANT position is a reference,
  // YaspGeometry can't be copied
  const FieldVector<ctype, cdim> & position; 

  const YaspGeometry<0,cdim,GridImp>&
  operator = (const YaspGeometry<0,cdim,GridImp>& g);
};

// operator<< for all YaspGeometrys
template <int mydim, int cdim, class GridImp>
inline
std::ostream& operator<< (std::ostream& s, YaspGeometry<mydim,cdim,GridImp>& e)
{
  e.print(s);
  return s;
}

//========================================================================
//========================================================================

template<int codim, int dim, class GridImp>
class YaspSpecialEntity :
  public GridImp::template Codim<codim>::Entity
{
public:
  typedef typename GridImp::ctype ctype;
  
  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;

  YaspSpecialEntity(const YGLI& g, const TSI& it) :
    GridImp::template Codim<codim>::Entity (YaspEntity<codim, dim, GridImp>(g,it))
    {};
  YaspSpecialEntity(const YaspEntity<codim, dim, GridImp>& e) :
    GridImp::template Codim<codim>::Entity (e)
    {};
  const TSI& transformingsubiterator () const
  {
        return this->realEntity.transformingsubiterator();
  }
  const YGLI& gridlevel () const
  {
        return this->realEntity.gridlevel();
  }
};

template<int codim, int dim, class GridImp>
class YaspEntity 
:  public EntityDefaultImplementation <codim,dim,GridImp,YaspEntity>
{
public:
  typedef typename GridImp::ctype ctype;
  
  typedef typename GridImp::template Codim<codim>::Geometry Geometry;
  int level () const
  {
        DUNE_THROW(GridError, "YaspEntity not implemented");
  }

  int index () const
  {
        DUNE_THROW(GridError, "YaspEntity not implemented");
  }

  const Geometry& geometry () const
  {
        DUNE_THROW(GridError, "YaspEntity not implemented");
  }

  PartitionType partitionType () const
  {
        DUNE_THROW(GridError, "YaspEntity not implemented");
  }

  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
  YaspEntity (const YGLI& g, const TSI& it)
  {
        DUNE_THROW(GridError, "YaspEntity not implemented");
  }

  // IndexSets needs access to the private index methods
  friend class Dune::YaspLevelIndexSet<GridImp>;
  friend class Dune::YaspLeafIndexSet<GridImp>;
  friend class Dune::YaspGlobalIdSet<GridImp>;
  typedef typename GridImp::PersistentIndexType PersistentIndexType;

  PersistentIndexType persistentIndex () const 
  { 
    DUNE_THROW(GridError, "YaspEntity not implemented");
  }

  int compressedIndex () const 
  {
    DUNE_THROW(GridError, "YaspEntity not implemented");
  }

  int compressedLeafIndex () const 
  {
    DUNE_THROW(GridError, "YaspEntity not implemented");
  }
};


// specialization for codim=0
template<int dim, class GridImp>
class YaspEntity<0,dim,GridImp> 
: public EntityDefaultImplementation <0,dim,GridImp,YaspEntity>
{
  enum { dimworld = GridImp::dimensionworld };
public:
  typedef typename GridImp::ctype ctype;
  
  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;

  typedef YaspSpecialGeometry<dim-0,dim,GridImp> SpecialGeometry;
  
  typedef typename GridImp::template Codim<0>::Geometry Geometry;
  template <int cd>
  struct Codim
  {
    typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
  };
  typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
  typedef typename GridImp::template Codim<0>::LevelIntersectionIterator IntersectionIterator;
  typedef typename GridImp::template Codim<0>::LevelIntersectionIterator LevelIntersectionIterator;
  typedef typename GridImp::template Codim<0>::LeafIntersectionIterator  LeafIntersectionIterator;
  typedef typename GridImp::template Codim<0>::HierarchicIterator HierarchicIterator;

  typedef typename GridImp::PersistentIndexType PersistentIndexType;

  typedef typename YGrid<dim,ctype>::iTupel iTupel;

  // constructor
  YaspEntity (const YGLI& g, const TSI& it)
    : _it(it), _g(g), _geometry(it.position(),it.meshsize())
  {
  }

  int level () const {return _g.level();}

  int index () const { return _it.superindex();} // superindex works also for iteration over subgrids

  int globalIndex () const {
    return _g.cell_global().index(_it.coord());
  }

  PartitionType partitionType () const
  {
        if (_g.cell_interior().inside(_it.coord())) return InteriorEntity;
        if (_g.cell_overlap().inside(_it.coord())) return OverlapEntity;
        return GhostEntity;
  }

  const Geometry& geometry () const { return _geometry; }

  template<int cc> int count () const
  {
    if (cc==dim) return 1<<dim;
    if (cc==1) return 2*dim;
    if (cc==dim-1) return dim*(1<<(dim-1));
    if (cc==0) return 1;
    DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
  }

  template<int cc>
  typename Codim<cc>::EntityPointer entity (int i) const
  {
      dune_static_assert( cc == dim || cc == 0 ,
          "YaspGrid only supports Entities with codim=dim and codim=0");
      // coordinates of the cell == coordinates of lower left corner
      if (cc==dim)
          {
                iTupel coord = _it.coord();

                // get corner from there
                for (int k=0; k<dim; k++)
                  if (i&(1<<k)) (coord[k])++;

                return YaspEntityPointer<cc,GridImp>(_g,_g.vertex_overlapfront().tsubbegin(coord));
          }
      if (cc==0)
          {
                return YaspEntityPointer<cc,GridImp>(_g,_it);
          }
      DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
  }

  EntityPointer father () const
  {
        // check if coarse level exists
        if (_g.level()<=0)
          DUNE_THROW(GridError, "tried to call father on level 0");

        // yes, get iterator to it
        YGLI cg = _g.coarser();

        // coordinates of the cell
        iTupel coord = _it.coord();

        // get coordinates on next coarser level
        for (int k=0; k<dim; k++) coord[k] = coord[k]/2;

        return YaspEntityPointer<0,GridImp>(cg,cg.cell_overlap().tsubbegin(coord));
  }

  const Geometry& geometryInFather () const
  {
    // determine which son we are
    int son = 0;
    for (int k=0; k<dim; k++)
      if (_it.coord(k)%2)
        son += (1<<k);
    
    // configure one of the 2^dim transformations
    return YaspFatherRelativeLocalElement<dim,GridImp>::getson(son);
  }

  const TSI& transformingsubiterator () const
  {
        return _it;
  }

  const YGLI& gridlevel () const
  {
        return _g;
  }

  bool isLeaf() const
  {
    return (_g.level() == _g.mg()->maxlevel());
  }
  
  IntersectionIterator ibegin () const
  {
    return YaspIntersectionIterator<GridImp>(*this,false);
  }

  LeafIntersectionIterator ileafbegin () const
  {
    // only if entity is leaf this iterator delivers intersections
    return YaspIntersectionIterator<GridImp>(*this, ! isLeaf() );
  }

  LevelIntersectionIterator ilevelbegin () const
  {
    return ibegin(); 
  }

  IntersectionIterator iend () const
  {
    return YaspIntersectionIterator<GridImp>(*this,true);
  }

  LeafIntersectionIterator ileafend () const
  {
    return iend();
  }

  LevelIntersectionIterator ilevelend () const
  {
    return iend();
  }

  HierarchicIterator hbegin (int maxlevel) const
  {
        return YaspHierarchicIterator<GridImp>(_g,_it,maxlevel);
  }

  HierarchicIterator hend (int maxlevel) const
  {
        return YaspHierarchicIterator<GridImp>(_g,_it,_g.level());
  }

private:
  // IndexSets needs access to the private index methods
  friend class Dune::YaspLevelIndexSet<GridImp>;
  friend class Dune::YaspLeafIndexSet<GridImp>;
  friend class Dune::YaspGlobalIdSet<GridImp>;

  PersistentIndexType persistentIndex () const 
  { 
    // get size of global grid
    const iTupel& size =  _g.cell_global().size();

    // get coordinate correction for periodic boundaries
    int coord[dim];
    for (int i=0; i<dim; i++)
      {
        coord[i] = _it.coord(i);
        if (coord[i]<0) coord[i] += size[i];
        if (coord[i]>=size[i]) coord[i] -= size[i];
      }

    // encode codim
    PersistentIndexType id(0);

    // encode level
    id = id << yaspgrid_level_bits;
    id = id+PersistentIndexType(_g.level());
    

    // encode coordinates
    for (int i=dim-1; i>=0; i--)
      {
        id = id << yaspgrid_dim_bits;
        id = id+PersistentIndexType(coord[i]);
      }
    
    return id;
  }

  int compressedIndex () const 
  {
    return _it.superindex();
  }

  int compressedLeafIndex () const 
  {
    return _it.superindex();
  }

  template<int cc>
  PersistentIndexType subPersistentIndex (int i) const
  {
    if (cc==0)
      return persistentIndex();

    // get position of cell, note that global origin is zero
    // adjust for periodic boundaries
    int coord[dim];
    for (int k=0; k<dim; k++)
      {
        coord[k] = _it.coord(k);
        if (coord[k]<0) coord[k] += _g.cell_global().size(k);
        if (coord[k]>=_g.cell_global().size(k)) coord[k] -= _g.cell_global().size(k);
      }
    
    if (cc==dim)
      {
        // transform to vertex coordinates
        for (int k=0; k<dim; k++)
          if (i&(1<<k)) (coord[k])++;

        // determine min number of trailing zeroes
        int trailing = 1000;
        for (int i=0; i<dim; i++)
          {
            // count trailing zeros
            int zeros = 0;
            for (int j=0; j<_g.level(); j++)
              if (coord[i]&(1<<j))
                break;
              else
                zeros++;
            trailing = std::min(trailing,zeros);
          } 

        // determine the level of this vertex
        int level = _g.level()-trailing;

        // encode codim
        PersistentIndexType id(dim);
        
        // encode level
        id = id << yaspgrid_level_bits;
        id = id+PersistentIndexType(level);
        
        // encode coordinates
        for (int i=dim-1; i>=0; i--)
          {
            id = id << yaspgrid_dim_bits;
            id = id+PersistentIndexType(coord[i]>>trailing);
          }
        
        return id;
      }

    if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
      {
        // Idea: Use the doubled grid to assign coordinates to faces

        // ivar is the direction that varies
        int ivar=i/2; 

        // compute position from cell position
        for (int k=0; k<dim; k++)
          coord[k] = coord[k]*2 + 1; // the doubled grid
        if (i%2) 
          coord[ivar] += 1;
        else
          coord[ivar] -= 1;

        // encode codim
        PersistentIndexType id(1);
        
        // encode level
        id = id << yaspgrid_level_bits;
        id = id+PersistentIndexType(_g.level());
        
        // encode coordinates
        for (int i=dim-1; i>=0; i--)
          {
            id = id << yaspgrid_dim_bits;
            id = id+PersistentIndexType(coord[i]);
          }
        
        return id;
      }

    if (cc==dim-1) // edges, exist only for dim>2
      {
        // Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction

        // number of entities per direction
        int m=1<<(dim-1);

        // ifix is the direction that is fixed
        int ifix=(dim-1)-(i/m); 

        // compute position from cell position
        int bit=1;
        for (int k=0; k<dim; k++)
          {
            coord[k] = coord[k]*2+1; // cell position in doubled grid
            if (k==ifix) continue;
            if ((i%m)&bit) coord[k] += 1; else coord[k] -= 1;
            bit *= 2;
          } 

        // encode codim
        PersistentIndexType id(dim-1);
        
        // encode level
        id = id << yaspgrid_level_bits;
        id = id+PersistentIndexType(_g.level());
        
        // encode coordinates
        for (int i=dim-1; i>=0; i--)
          {
            id = id << yaspgrid_dim_bits;
            id = id+PersistentIndexType(coord[i]);
          }
        
        return id;
      }

    DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
  }

  template<int cc>
  int subCompressedIndex (int i) const
  {
    if (cc==0)
      return compressedIndex();
    
    // get cell position relative to origin of local cell grid
    iTupel coord;
    for (int k=0; k<dim; ++k) 
      coord[k] = _it.coord(k)-_g.cell_overlap().origin(k);

    if (cc==dim) // vertices
      {
        // transform cell coordinate to corner coordinate
        for (int k=0; k<dim; k++)
          if (i&(1<<k)) (coord[k])++;

        // do lexicographic numbering
        int index = coord[dim-1];
        for (int k=dim-2; k>=0; --k)
          index = (index*(_g.cell_overlap().size(k)+1))+coord[k];
        return index;
      }

    if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
      {
        // Idea: direction ivar varies, all others are fixed, i.e. 2 possibilities per direction

        // ivar is the direction that varies
        int ivar=i/2; 

        // compute position from cell position
        if (i%2) coord[ivar] += 1; 

        // do lexicographic numbering
        int index = coord[dim-1];
        for (int k=dim-2; k>=0; --k)
          if (k==ivar)
            index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
          else
            index = (index*(_g.cell_overlap().size(k)))+coord[k];

        // add size of all subsets for smaller directions
        for (int j=0; j<ivar; j++)
          {
            int n=_g.cell_overlap().size(j)+1;
            for (int l=0; l<dim; l++)
              if (l!=j) n *= _g.cell_overlap().size(l);
            index += n;
          }

        return index;
      }

    if (cc==dim-1) // edges, exist only for dim>2
      {
        // Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction

        // number of entities per direction
        int m=1<<(dim-1);

        // ifix is the direction that is fixed
        int ifix=(dim-1)-(i/m); 

        // compute position from cell position
        int bit=1;
        for (int k=0; k<dim; k++)
          {
            if (k==ifix) continue;
            if ((i%m)&bit) coord[k] += 1;
            bit *= 2;
          } 

        // do lexicographic numbering
        int index = coord[dim-1];
        for (int k=dim-2; k>=0; --k)
          if (k!=ifix)
            index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
          else
            index = (index*(_g.cell_overlap().size(k)))+coord[k];

        // add size of all subsets for smaller directions
        for (int j=dim-1; j>ifix; j--)
          {
            int n=_g.cell_overlap().size(j);
            for (int l=0; l<dim; l++)
              if (l!=j) n *= _g.cell_overlap().size(l)+1;
            index += n;
          }

        return index;
      }

    DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
  }

  template<int cc>
  int subCompressedLeafIndex (int i) const
  {
    if (cc==0)
      return compressedIndex();
    
    // get cell position relative to origin of local cell grid
    iTupel coord;
    for (int k=0; k<dim; ++k) 
      coord[k] = _it.coord(k)-_g.cell_overlap().origin(k);

    if (cc==dim) // vertices
      {
        // transform cell coordinate to corner coordinate
        for (int k=0; k<dim; k++)
          if (i&(1<<k)) (coord[k])++;

        // move coordinates up to maxlevel
        for (int k=0; k<dim; k++)
          coord[k] = coord[k]<<(_g.mg()->maxlevel()-_g.level());

        // do lexicographic numbering
        int index = coord[dim-1];
        for (int k=dim-2; k>=0; --k)
          index = (index*(_g.mg()->rbegin().cell_overlap().size(k)+1))+coord[k];
        return index;
      }

    if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
      {
        // Idea: direction ivar varies, all others are fixed, i.e. 2 possibilities per direction

        // ivar is the direction that varies
        int ivar=i/2; 

        // compute position from cell position
        if (i%2) coord[ivar] += 1; 

        // do lexicographic numbering
        int index = coord[dim-1];
        for (int k=dim-2; k>=0; --k)
          if (k==ivar)
            index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
          else
            index = (index*(_g.cell_overlap().size(k)))+coord[k];

        // add size of all subsets for smaller directions
        for (int j=0; j<ivar; j++)
          {
            int n=_g.cell_overlap().size(j)+1;
            for (int l=0; l<dim; l++)
              if (l!=j) n *= _g.cell_overlap().size(l);
            index += n;
          }

        return index;
      }

    if (cc==dim-1) // edges, exist only for dim>2
      {
        // Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction

        // number of entities per direction
        int m=1<<(dim-1);

        // ifix is the direction that is fixed
        int ifix=(dim-1)-(i/m); 

        // compute position from cell position
        int bit=1;
        for (int k=0; k<dim; k++)
          {
            if (k==ifix) continue;
            if ((i%m)&bit) coord[k] += 1;
            bit *= 2;
          } 

        // do lexicographic numbering
        int index = coord[dim-1];
        for (int k=dim-2; k>=0; --k)
          if (k!=ifix)
            index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
          else
            index = (index*(_g.cell_overlap().size(k)))+coord[k];

        // add size of all subsets for smaller directions
        for (int j=dim-1; j>ifix; j--)
          {
            int n=_g.cell_overlap().size(j);
            for (int l=0; l<dim; l++)
              if (l!=j) n *= _g.cell_overlap().size(l)+1;
            index += n;
          }

        return index;
      }

    DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
  }
  
  const TSI& _it;           // position in the grid level
  const YGLI& _g;           // access to grid level
  SpecialGeometry _geometry; // the element geometry
};


// specialization for codim=dim (vertex)
template<int dim, class GridImp>
class YaspEntity<dim,dim,GridImp> 
: public EntityDefaultImplementation <dim,dim,GridImp,YaspEntity>
{
  enum { dimworld = GridImp::dimensionworld };
public:
  typedef typename GridImp::ctype ctype;
  
  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;

  typedef YaspSpecialGeometry<dim-dim,dim,GridImp> SpecialGeometry;
  
  typedef typename GridImp::template Codim<dim>::Geometry Geometry;
  template <int cd>
  struct Codim
  {
    typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
  };
  typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;

  typedef typename GridImp::PersistentIndexType PersistentIndexType;

  typedef typename YGrid<dim,ctype>::iTupel iTupel;

  // constructor
  YaspEntity (const YGLI& g, const TSI& it)
        : _it(it), _g(g), _geometry(it.position())
  {  }

  int level () const {return _g.level();}

  int index () const {return _it.superindex();}

  int globalIndex () const { return _g.cell_global().index(_it.coord()); }


  const Geometry& geometry () const { return _geometry; }

  PartitionType partitionType () const
  {
        if (_g.vertex_interior().inside(_it.coord())) return InteriorEntity;
        if (_g.vertex_interiorborder().inside(_it.coord())) return BorderEntity;
        if (_g.vertex_overlap().inside(_it.coord())) return OverlapEntity;
        if (_g.vertex_overlapfront().inside(_it.coord())) return FrontEntity;
        return GhostEntity;
  }

private:
  // IndexSets needs access to the private index methods
  friend class Dune::YaspLevelIndexSet<GridImp>;
  friend class Dune::YaspLeafIndexSet<GridImp>;
  friend class Dune::YaspGlobalIdSet<GridImp>;

  PersistentIndexType persistentIndex () const 
  { 
    // get coordinate and size of global grid
    const iTupel& size =  _g.vertex_global().size();
    int coord[dim];

    // correction for periodic boundaries
    for (int i=0; i<dim; i++)
      {
        coord[i] = _it.coord(i);
        if (coord[i]<0) coord[i] += size[i];
        if (coord[i]>=size[i]) coord[i] -= size[i];
      }

    // determine min number of trailing zeroes
    int trailing = 1000;
    for (int i=0; i<dim; i++)
      {
        // count trailing zeros
        int zeros = 0;
        for (int j=0; j<_g.level(); j++)
          if (coord[i]&(1<<j))
            break;
          else
            zeros++;
        trailing = std::min(trailing,zeros);
      } 

    // determine the level of this vertex
    int level = _g.level()-trailing;

    // encode codim
    PersistentIndexType id(dim);

    // encode level
    id = id << yaspgrid_level_bits;
    id = id+PersistentIndexType(level);
    
    // encode coordinates
    for (int i=dim-1; i>=0; i--)
      {
        id = id << yaspgrid_dim_bits;
        id = id+PersistentIndexType(coord[i]>>trailing);
      }
    
    return id;
  }

  int compressedIndex () const { return _it.superindex();}

  int compressedLeafIndex () const 
  { 
    if (_g.level()==_g.mg()->maxlevel())
      return _it.superindex();

    // not on leaf level, interpolate to finest grid
    int coord[dim];
    for (int i=0; i<dim; i++) coord[i] = _it.coord(i)-(_g).vertex_overlap().origin(i);

    // move coordinates up to maxlevel (multiply by 2 for each level
    for (int k=0; k<dim; k++)
      coord[k] = coord[k]*(1<<(_g.mg()->maxlevel()-_g.level()));

    // do lexicographic numbering
    int index = coord[dim-1];
    for (int k=dim-2; k>=0; --k)
      index = (index*(_g.mg()->rbegin().cell_overlap().size(k)+1))+coord[k];
    return index;
  }

public:  
  const TSI& transformingsubiterator() const { return _it; }
  const YGLI& gridlevel() const { return _g; }
protected:
  const TSI& _it;                // position in the grid level
  const YGLI& _g;                // access to grid level
  SpecialGeometry  _geometry;     // the element geometry
  // temporary object
  mutable FieldVector<ctype, dim> loc;   // always computed before being returned
};


//========================================================================
//========================================================================

template<class GridImp>
class YaspIntersectionIterator
{
  enum { dim=GridImp::dimension };
  enum { dimworld=GridImp::dimensionworld };  
  typedef typename GridImp::ctype ctype;
  YaspIntersectionIterator();
public:
  // types used from grids
  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
  typedef typename GridImp::template Codim<0>::Entity Entity;
  typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
  typedef typename GridImp::template Codim<1>::Geometry Geometry;
  typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
  typedef YaspSpecialEntity<0,dim,GridImp> SpecialEntity;
  typedef YaspSpecialGeometry<dim-1,dimworld,GridImp> SpecialGeometry;
  typedef YaspSpecialGeometry<dim-1,dim,GridImp> SpecialLocalGeometry;
  typedef Dune::Intersection<const GridImp, Dune::YaspIntersectionIterator> Intersection;
  
  void increment()
  {
        // check end
        if (_count==2*dim) return; // end iterator reached, we are done
        // update count, check end
        _count++;

        // update intersec