• 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

sgrid.hh

Go to the documentation of this file.

#ifndef DUNE_SGRID_HH
#define DUNE_SGRID_HH

#include <limits>
#include <vector>
#include <stack>

#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/finitestack.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/common/bigunsignedint.hh>
#include <dune/common/collectivecommunication.hh>
#include <dune/grid/common/grid.hh>
#include <dune/grid/sgrid/numbering.hh>
#include <dune/grid/common/indexidset.hh>

namespace Dune {

//************************************************************************
  typedef double sgrid_ctype; 

  // globally define the persistent index type
  const int sgrid_dim_bits = 24;   // bits for encoding each dimension
  const int sgrid_level_bits = 6;  // bits for encoding level number
  const int sgrid_codim_bits = 4;  // bits for encoding codimension

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

template<int dim, int dimworld, class GridImp> class SGeometry;
template<int codim, int dim, class GridImp> class SEntity;
template<int codim, class GridImp> class SEntityPointer;
template<int codim, PartitionIteratorType, class GridImp> class SLevelIterator;
template<int dim, int dimworld> class SGrid;
template<class GridImp> class SIntersectionIterator;
template<class GridImp> class SHierarchicIterator;

//************************************************************************
template<int mydim, int cdim, class GridImp> 
class SGeometry : public GeometryDefaultImplementation<mydim,cdim,GridImp,SGeometry>
{
public:
        typedef sgrid_ctype ctype;

        GeometryType type () const;             

        int corners () const;                  

        const FieldVector<sgrid_ctype, cdim>& operator[] (int i) const;

        FieldVector<sgrid_ctype, cdim> global (const FieldVector<sgrid_ctype, mydim>& local) const;

        FieldVector<sgrid_ctype, mydim> local (const FieldVector<sgrid_ctype, cdim>& global) const;

        bool checkInside (const FieldVector<sgrid_ctype, mydim>& local) const;
  
        sgrid_ctype integrationElement (const FieldVector<sgrid_ctype, mydim>& local) const;

        const FieldMatrix<sgrid_ctype,cdim,mydim>& jacobianInverseTransposed (const FieldVector<sgrid_ctype, mydim>& local) const;

        void print (std::ostream& ss, int indent) const;

        void make (FieldMatrix<sgrid_ctype,mydim+1,cdim>& __As);

        SGeometry () : builtinverse(false) {};

private:
        FieldVector<sgrid_ctype, cdim> s;         
        FieldMatrix<sgrid_ctype,mydim,cdim> A;     
        array<FieldVector<sgrid_ctype, cdim>, 1<<mydim> c; 
        mutable FieldMatrix<sgrid_ctype,cdim,mydim> Jinv;       
        mutable bool builtinverse;
};

template<int cdim, class GridImp> 
class SGeometry<0,cdim,GridImp> : public GeometryDefaultImplementation<0,cdim,GridImp,SGeometry>
{
public:
        typedef sgrid_ctype ctype;

        GeometryType type () const;             

        int corners () const;

        const FieldVector<sgrid_ctype, cdim>& operator[] (int i) const;

        void print (std::ostream& ss, int indent) const;

        void make (FieldMatrix<sgrid_ctype,1,cdim>& __As);

        FieldVector<sgrid_ctype, cdim> global (const FieldVector<sgrid_ctype, 0>& local) const { return this->operator[] (0); }

        FieldVector<sgrid_ctype, 0> local (const FieldVector<sgrid_ctype, cdim>& global) const { return FieldVector<sgrid_ctype,0> (0.0); } 


        SGeometry () {};

        sgrid_ctype integrationElement(const FieldVector<sgrid_ctype, 0>& local) const {
            return 1.;
        }
        
        bool checkInside (const FieldVector<sgrid_ctype, 0>& local) const 
        {
          return true; 
        }
  
        const FieldMatrix<sgrid_ctype,cdim,0>& jacobianInverseTransposed (const FieldVector<sgrid_ctype, 0>& local) const
        {
          static FieldMatrix<sgrid_ctype,cdim,0> dummy( sgrid_ctype(0) );
          return dummy;
        }

protected:
        FieldVector<sgrid_ctype, cdim> s;         
};

template <int mydim, int cdim, class GridImp>
inline std::ostream& operator<< (std::ostream& s, SGeometry<mydim,cdim,GridImp>& e)
{
        e.print(s,0);
        return s;
}

template<int mydim, int cdim, class GridImp>
class SMakeableGeometry : public Geometry<mydim, cdim, GridImp, SGeometry>
{
public:
  SMakeableGeometry() :
    Geometry<mydim, cdim, GridImp, SGeometry>(SGeometry<mydim, cdim, GridImp>())
    {};

    void make (FieldMatrix<sgrid_ctype,mydim+1,cdim>& __As) { this->realGeometry.make(__As); }
};

//************************************************************************
template<int codim, int dim, class GridImp>
class SEntityBase {
  friend class SEntityPointer<codim,GridImp>;
  friend class SIntersectionIterator<GridImp>;
  enum { dimworld = GridImp::dimensionworld };
public:
  typedef typename GridImp::template Codim<codim>::Geometry Geometry;
  typedef SMakeableGeometry<dim-codim, dimworld, const GridImp> MakeableGeometry;
  typedef typename GridImp::PersistentIndexType PersistentIndexType;

  int level () const
    {
      return l;
    }
  
  int index () const
    {
      return compressedIndex();
    }
  
  int globalIndex() const;
  
  const Geometry& geometry () const;
  
  SEntityBase (GridImp* _grid, int _l, int _id);
  SEntityBase ();

  void make (GridImp* _grid, int _l, int _id);

  void make (int _l, int _id);

  PersistentIndexType persistentIndex () const 
    { 
          if (codim!=dim)
                {
                  // encode codim, this would actually not be necessary
          // because z is unique in codim
                  PersistentIndexType id(codim);

                  // encode level
                  id = id << sgrid_level_bits;
                  id = id+PersistentIndexType(l);
        
                  // encode coordinates
                  for (int i=dim-1; i>=0; i--)
                        {
                          id = id << sgrid_dim_bits;
                          id = id+PersistentIndexType(z[i]);
                        }
                  
                  return id;
                }
          else
                {
                  // determine min number of trailing zeroes
                  // consider that z is on the doubled grid !
                  int trailing = 1000;
                  for (int i=0; i<dim; i++)
                        {
                          // count trailing zeros
                          int zeros = 0;
                          for (int j=0; j<l; j++)
                                if (z[i]&(1<<(j+1)))
                                  break;
                                else
                                  zeros++;
                          trailing = std::min(trailing,zeros);
                        }

                  // determine the level of this vertex
                  int level = l-trailing;

                  // encode codim
                  PersistentIndexType id(dim);

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

  int compressedIndex () const
    {
                return id;
    }

  int compressedLeafIndex () const
    {
          // codim != dim -> there are no copies of entities
          // maxlevel -> ids are fine
          if (codim<dim || l==grid->maxLevel())
                return id;
          
          // this is a vertex which is not on the finest level
          // move coordinates up to maxlevel (multiply by 2 for each level
          array<int,dim> coord;
          for (int k=0; k<dim; k++)
                coord[k] = z[k]*(1<<(grid->maxLevel()-l));

          // compute number with respect to maxLevel
          return grid->n(grid->maxLevel(),coord);
    }

protected:
  // this is how we implement our elements
  GridImp* grid;         
  int l;                 
  int id;                
  array<int,dim> z; 
  mutable MakeableGeometry geo; 
  mutable bool builtgeometry;   
};


template<int codim, int dim, class GridImp> 
class SEntity : public SEntityBase<codim,dim,GridImp>, 
                public EntityDefaultImplementation<codim,dim,GridImp,SEntity>
{
  enum { dimworld = GridImp::dimensionworld };
public:
  typedef typename GridImp::template Codim<codim>::Geometry Geometry;
  typedef typename GridImp::template Codim<codim>::LevelIterator LevelIterator;
  typedef typename GridImp::template Codim<0>::LeafIntersectionIterator IntersectionIterator;
  typedef typename GridImp::template Codim<0>::HierarchicIterator HierarchicIterator;
  
  // disambiguate member functions with the same name in both bases
  int level () const {return SEntityBase<codim,dim,GridImp>::level();}
  
  int index () const {return SEntityBase<codim,dim,GridImp>::index();}

  const Geometry& geometry () const { return SEntityBase<codim,dim,GridImp>::geometry(); }
  
  PartitionType partitionType () const { return InteriorEntity; }

  // specific to SEntity
  SEntity (GridImp* _grid, int _l, int _id) : SEntityBase<codim,dim,GridImp>::SEntityBase(_grid,_l,_id) {};
};

template<int dim, class GridImp>
class SEntity<0,dim,GridImp> : public SEntityBase<0,dim,GridImp>, 
                               public EntityDefaultImplementation<0,dim,GridImp,SEntity>
{
  enum { dimworld = GridImp::dimensionworld };
public:
  typedef typename GridImp::template Codim<0>::Geometry Geometry;
  typedef typename GridImp::template Codim<0>::LocalGeometry LocalGeometry;
  typedef SMakeableGeometry<dim, dimworld, const GridImp> MakeableGeometry;
  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>::LeafIntersectionIterator IntersectionIterator;
  typedef typename GridImp::template Codim<0>::HierarchicIterator HierarchicIterator;
  typedef typename GridImp::PersistentIndexType PersistentIndexType;

  friend class SHierarchicIterator<GridImp>;
  
  // disambiguate member functions with the same name in both bases
  int level () const {return SEntityBase<0,dim,GridImp>::level();}
  
  int index () const {return SEntityBase<0,dim,GridImp>::index();}
  
  PartitionType partitionType () const { return InteriorEntity; }

  const Geometry& geometry () const {return SEntityBase<0,dim,GridImp>::geometry();}

  template<int cc> int count () const; 

  template<int cc> typename Codim<cc>::EntityPointer entity (int i) const;

  template<int cc>
  int subCompressedIndex (int i) const
  {
        if (cc==0) return this->compressedIndex();
        return this->grid->template getRealEntity<cc>(*entity<cc>(i)).compressedIndex();
  }

  template<int cc>
  int subCompressedLeafIndex (int i) const
  {
        if (cc==0) return this->compressedLeafIndex();
        return this->grid->template getRealEntity<cc>(*entity<cc>(i)).compressedLeafIndex();
  }

  template<int cc>
  PersistentIndexType subPersistentIndex (int i) const
  {
        if (cc==0) return this->persistentIndex();
        return this->grid->template getRealEntity<cc>(*entity<cc>(i)).persistentIndex();
  }
  
  IntersectionIterator ibegin () const;
  IntersectionIterator ileafbegin () const;
  IntersectionIterator ilevelbegin () const;
  IntersectionIterator iend () const;
  IntersectionIterator ileafend () const;
  IntersectionIterator ilevelend () const;

  EntityPointer father () const;

  bool isLeaf () const
    {
      return ( this->grid->maxLevel() == level() );
    }

  const LocalGeometry& geometryInFather () const;

  HierarchicIterator hbegin (int maxLevel) const;

  HierarchicIterator hend (int maxLevel) const;

  // members specific to SEntity
  SEntity (GridImp* _grid, int _l, int _id) : 
    SEntityBase<0,dim,GridImp>::SEntityBase(_grid,_l,_id)
    {
      built_father = false;
    }

  SEntity ()
    {
      built_father = false;
    }

  void make (GridImp* _grid, int _l, int _id)
    {
      SEntityBase<0,dim,GridImp>::make(_grid,_l,_id);
      built_father = false;
    }

  void make (int _l, int _id)
    {
      SEntityBase<0,dim,GridImp>::make(_l,_id);
      built_father = false;
    }
  
private:

  mutable bool built_father;
  mutable int father_id;
    mutable SMakeableGeometry<dim,dim,const GridImp> in_father_local;
  void make_father() const;
};

template<int dim, class GridImp>
class SEntity<dim,dim,GridImp> : public SEntityBase<dim,dim,GridImp>, 
                                 public EntityDefaultImplementation <dim,dim,GridImp,SEntity>
{
  enum { dimworld = GridImp::dimensionworld };
public:
  typedef typename GridImp::template Codim<dim>::Geometry Geometry;
  typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;

  // disambiguate member functions with the same name in both bases
  int level () const {return SEntityBase<dim,dim,GridImp>::level();}

  int index () const {return SEntityBase<dim,dim,GridImp>::index();}
  
  PartitionType partitionType () const { return InteriorEntity; }

  const Geometry& geometry () const {return SEntityBase<dim,dim,GridImp>::geometry();}

        // members specific to SEntity
        SEntity (GridImp* _grid, int _l, int _id) : SEntityBase<dim,dim,GridImp>::SEntityBase(_grid,_l,_id)
        {
                built_father = false;
        }

        void make (GridImp* _grid, int _l, int _id)
        {
                SEntityBase<dim,dim,GridImp>::make(_grid, _l,_id);
                built_father = false;
        }

        void make (int _l, int _id)
        {
                SEntityBase<dim,dim,GridImp>::make(_l,_id);
                built_father = false;
        }

private:
  mutable bool built_father;
  mutable int father_id;
  mutable FieldVector<sgrid_ctype, dim> in_father_local;
  void make_father() const;
};

template<int codim, int dim, class GridImp>
class SMakeableEntity :
  public GridImp::template Codim<codim>::Entity
{
public:

  SMakeableEntity(GridImp* _grid, int _l, int _id) :
    GridImp::template Codim<codim>::Entity (SEntity<codim, dim, GridImp>(_grid,_l,_id))
    {};
  SMakeableEntity(const SEntity<codim, dim, GridImp>& e) :
    GridImp::template Codim<codim>::Entity (e)
    {};
  void make (int _l, int _id) { this->realEntity.make(_l, _id); }
  void make (GridImp* _grid, int _l, int _id) { this->realEntity.make(_grid, _l, _id); }
};

//************************************************************************
struct SHierarchicStackElem {
  int l;
  int id;
  SHierarchicStackElem () : l(-1), id(-1) {}
  SHierarchicStackElem (int _l, int _id) {l=_l; id=_id;}
  bool operator== (const SHierarchicStackElem& s) const {return !operator!=(s);}
  bool operator!= (const SHierarchicStackElem& s) const {return l!=s.l || id!=s.id;}
};

template<class GridImp>
class SHierarchicIterator :
  public Dune::SEntityPointer <0,GridImp>
{
  friend class SHierarchicIterator<const GridImp>;
  enum { dim = GridImp::dimension };
  enum { dimworld = GridImp::dimensionworld };
public:
  typedef typename GridImp::template Codim<0>::Entity Entity;
  typedef typename GridImp::ctype ctype;

  void increment();

  SHierarchicIterator (GridImp* _grid,
                       const Dune::SEntity<0,GridImp::dimension,GridImp>& _e,
                       int _maxLevel, bool makeend) :
    Dune::SEntityPointer<0,GridImp>(_grid,_e.level(),_e.index())
    {
      // without sons, we are done
      // (the end iterator is equal to the calling iterator)
      if (makeend) return;
      
      // remember element where begin has been called
      orig_l = this->entity().level();
      orig_id = _grid->template getRealEntity<0>(this->entity()).index();
      
      // push original element on stack
      SHierarchicStackElem originalElement(orig_l, orig_id);
      stack.push(originalElement);
      
      // compute maxLevel
      maxLevel = std::min(_maxLevel,this->grid->maxLevel());
      
      // ok, push all the sons as well
      push_sons(orig_l,orig_id);
      
      // and pop the first son
      increment();
    }

private:
  int maxLevel;                
  int orig_l, orig_id;         

  FiniteStack<SHierarchicStackElem,GridImp::MAXL> stack;      
  
  void push_sons (int level, int fatherid); 
};

//************************************************************************
template<class GridImp>
class SIntersectionIterator
{
  enum { dim=GridImp::dimension };
  enum { dimworld=GridImp::dimensionworld };
public:
  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 Dune::Intersection<const GridImp, Dune::SIntersectionIterator> Intersection;
  enum { dimension=dim };
  enum { dimensionworld=dimworld };
  typedef typename GridImp::ctype ctype;

  bool equals(const SIntersectionIterator<GridImp>& i) const;
  void increment();

  const Intersection & dereference() const
  {
    return reinterpret_cast<const Intersection&>(*this);
  }
  
  EntityPointer inside() const;

  EntityPointer outside() const;
  
  bool boundary () const;

  bool conforming () const;

  int boundaryId () const {
    if (boundary()) return count + 1;
    return 0;
  };
  bool neighbor () const;

  FieldVector<ctype, GridImp::dimensionworld> outerNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const
    {
      return unitOuterNormal(local);
    }
  FieldVector<ctype, GridImp::dimensionworld> unitOuterNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const
    {
      // while we are at it, compute normal direction
      FieldVector<sgrid_ctype, dimworld> normal(0.0);
      if (count%2)
        normal[count/2] =  1.0; // odd
      else
        normal[count/2] = -1.0; // even
      
      return normal; 
    }
  FieldVector<ctype, GridImp::dimensionworld> integrationOuterNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const
    {
        FieldVector<sgrid_ctype, dimworld> n = unitOuterNormal(local);
        n *= is_global.integrationElement(local);
        return n;
    }

  LocalGeometry& intersectionSelfLocal () const;
  LocalGeometry& intersectionNeighborLocal () const;
  Geometry& intersectionGlobal () const;
  int numberInSelf () const;
  int numberInNeighbor () const;

  SIntersectionIterator (GridImp* _grid, const SEntity<0,dim,GridImp >* _self, int _count);

  SIntersectionIterator& operator = (const SIntersectionIterator& it)
    {
      assert(grid == it.grid);
      
      /* Assert same Iterator Context */
      if (! self.equals(it.self))
        DUNE_THROW(GridError, "assignment of SIntersectionIterator "
                   << "with different inside Entity");
      if (partition != it.partition)
        DUNE_THROW(GridError, "assignment of SIntersectionIterator "
                   << "with different partition");

      /* Assign current position and update ne */
      ne = it.ne;
      count = it.count;
      make(count);

      return *this;
    }
  
private:
  void make (int _count) const;                 
  void makeintersections () const;              
  const EntityPointer self;                     
  mutable EntityPointer ne;                     
  const GridImp * grid;                         
  const int partition;                          
  const array<int,dim> zred;                 
  mutable int count;                              
  mutable bool valid_count;                       
  mutable bool valid_nb;                          
  mutable bool is_on_boundary;                    
  mutable bool built_intersections;               
  mutable SMakeableGeometry<dim-1,dim,const GridImp> is_self_local;      
  mutable SMakeableGeometry<dim-1,dimworld,const GridImp> is_global;     
  mutable SMakeableGeometry<dim-1,dim,const GridImp> is_nb_local;        
};

//************************************************************************

template<int codim, class GridImp>
class SEntityPointer
{
  enum { dim = GridImp::dimension };
  friend class SIntersectionIterator<GridImp>;
public:
  typedef SEntityPointer<codim,GridImp> EntityPointerImp;
  typedef typename GridImp::template Codim<codim>::Entity Entity;
  enum { codimension = codim };
  
  bool equals(const SEntityPointer<codim,GridImp>& i) const;
  Entity& dereference() const;
  int level () const;

  SEntityPointer (GridImp * _grid, int _l, int _id) :
    grid(_grid), l(_l), id(_id), 
    e(0) 
  {
  }

  SEntityPointer (const SEntity<codim,dim,GridImp> & _e) :
    grid(_e.grid), l(_e.l), id(_e.id),
    e(0)
  {
  }

  SEntityPointer (const SEntityPointer<codim,GridImp>& other) :
    grid(other.grid), l(other.l), id(other.id),
    e( 0 )
  {
  }

  ~SEntityPointer()
  {
    compactify();
#ifndef NDEBUG 
    id = -1;
#endif
  }

  SEntityPointer& operator = (const SEntityPointer& other)
  {
    grid = other.grid;
    l = other.l;
    id = other.id;

    compactify();
    return *this;
  }

  void compactify()
  {
    if( e )
    {
      enStack().push( e );
      e = 0;
    }
  }
  
protected:
  SMakeableEntity<codim,dim,GridImp>& entity() const
  {
    if( ! e )
    {
      e = getEntity( grid, l, id );
    }
    return *e;
  }

  typedef std::stack< SMakeableEntity<codim,dim,GridImp> * > EntityStackType;
  static inline EntityStackType& enStack()
  {
    static EntityStackType eStack;
    return eStack;
  }

  inline SMakeableEntity<codim,dim,GridImp>* getEntity(GridImp* _grid, int _l, int _id ) const
  {
    // get stack reference 
    EntityStackType& enSt = enStack();

    if( enSt.empty() )
    {
      return (new SMakeableEntity<codim,dim,GridImp>(_grid, _l, _id));
    }
    else
    {
      SMakeableEntity<codim,dim,GridImp>* e = enSt.top();
      enSt.pop();
      e->make(_grid, _l,_id);
      return e;
    }
  }

  GridImp* grid;                 
  int l;                         
  mutable int id;                
  mutable SMakeableEntity<codim,dim,GridImp>* e;  
};
//************************************************************************


template<int codim, PartitionIteratorType pitype, class GridImp>
class SLevelIterator :
  public Dune::SEntityPointer <codim,GridImp>
{
  friend class SLevelIterator<codim, pitype,const GridImp>;
  enum { dim = GridImp::dimension };
public:
  typedef typename GridImp::template Codim<codim>::Entity Entity;

  void increment();

  SLevelIterator (GridImp * _grid, int _l, int _id) :
    Dune::SEntityPointer<codim,GridImp>(_grid,_l,_id) {}
};


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

template <class GridImp> 
struct SGridLevelIndexSetTypes
{
  template<int cd>
  struct Codim
  {
        template<PartitionIteratorType pitype>
        struct Partition
        {
          typedef typename GridImp::Traits::template Codim<cd>::template Partition<pitype>::LevelIterator Iterator;
        };
  };
};

template<class GridImp>
class SGridLevelIndexSet : public IndexSetDefaultImplementation<GridImp,SGridLevelIndexSet<GridImp>,SGridLevelIndexSetTypes<GridImp> >
{
  typedef IndexSetDefaultImplementation<GridImp,SGridLevelIndexSet<GridImp>,SGridLevelIndexSetTypes<GridImp> > Base;
public:

  SGridLevelIndexSet (const GridImp& g, int l) : grid(g), level(l)
  {
    // contains a single element type;
        for (int codim=0; codim<=GridImp::dimension; codim++)
          mytypes[codim].push_back(GeometryType(GeometryType::cube,GridImp::dimension-codim));
  }

  template<int cd>
  int index (const typename GridImp::Traits::template Codim<cd>::Entity& e) const 
  {
        return grid.template getRealEntity<cd>(e).compressedIndex(); 
  }

  template<int cc>
  int subIndex (const typename GridImp::Traits::template Codim<0>::Entity& e, int i) const
  {
        return grid.template getRealEntity<0>(e).template subCompressedIndex<cc>(i);
  }

  int size (GeometryType type) const
  {
      return grid.size(level,GridImp::dimension-type.dim());
  }
  
  int size (int codim) const
  {
    return Base::size(codim);
  }

  const std::vector<GeometryType>& geomTypes (int codim) const
  {
        return mytypes[codim];
  }

  template<int cd, PartitionIteratorType pitype>
  typename Base::template Codim<cd>::template Partition<pitype>::Iterator begin () const
  {
        return grid.template lbegin<cd,pitype>(level);
  }
  
  template<int cd, PartitionIteratorType pitype>
  typename Base::template Codim<cd>::template Partition<pitype>::Iterator end () const
  {
        return grid.template lend<cd,pitype>(level);
  }

private:
  const GridImp& grid;
  int level;
  std::vector<GeometryType> mytypes[GridImp::dimension+1];
};

// Leaf Index Set

template <class GridImp> 
struct SGridLeafIndexSetTypes
{
  template<int cd>
  struct Codim
  {
        template<PartitionIteratorType pitype>
        struct Partition
        {
          typedef typename GridImp::Traits::template Codim<cd>::template Partition<pitype>::LeafIterator Iterator;
        };
  };
};

template<class GridImp>
class SGridLeafIndexSet : public IndexSetDefaultImplementation<GridImp,SGridLeafIndexSet<GridImp>,SGridLeafIndexSetTypes<GridImp> >
{
  typedef IndexSetDefaultImplementation<GridImp,SGridLeafIndexSet<GridImp>,SGridLeafIndexSetTypes<GridImp> > Base;
    enum {dim = remove_const<GridImp>::type::dimension};
public:

  SGridLeafIndexSet (const GridImp& g) : grid(g)
  {
    // contains a single element type;
        for (int codim=0; codim<=dim; codim++)
          mytypes[codim].push_back(GeometryType(GeometryType::cube,dim-codim));
  }

  /*
    We use the remove_const to extract the Type from the mutable class,
    because the const class is not instantiated yet.
  */
  template<int cd>
  int index (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const 
  {
        return grid.template getRealEntity<cd>(e).compressedLeafIndex(); 
  }

  /*
    We use the remove_const to extract the Type from the mutable class,
    because the const class is not instantiated yet.
  */
  template<int cc>
  int subIndex (const typename remove_const<GridImp>::type::Traits::template Codim<0>::Entity& e, int i) const
  {
        return grid.template getRealEntity<0>(e).template subCompressedLeafIndex<cc>(i);
  }

  int size (GeometryType type) const
  {
      return grid.size(grid.maxLevel(),GridImp::dimension-type.dim());
  }

  int size (int codim) const
  {
    return Base::size(codim);
  }

  const std::vector<GeometryType>& geomTypes (int codim) const
  {
        return mytypes[codim];
  }

  template<int cd, PartitionIteratorType pitype>
  typename Base::template Codim<cd>::template Partition<pitype>::Iterator begin () const
  {
        return grid.template leafbegin<cd,pitype>();
  }
  
  template<int cd, PartitionIteratorType pitype>
  typename Base::template Codim<cd>::template Partition<pitype>::Iterator end () const
  {
        return grid.template leafend<cd,pitype>();
  }

private:
  const GridImp& grid;
  std::vector<GeometryType> mytypes[dim+1];
};


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

template<class GridImp>
class SGridGlobalIdSet :
  public IdSetDefaultImplementation<GridImp,SGridGlobalIdSet<GridImp>, typename remove_const<GridImp>::type::PersistentIndexType>
  /*
    We used the remove_const to extract the Type from the mutable class,
    because the const class is not instantiated yet.
  */
{
public:
  /*
    We use the remove_const to extract the Type from the mutable class,
    because the const class is not instantiated yet.
  */
  typedef typename remove_const<GridImp>::type::PersistentIndexType IdType;

  SGridGlobalIdSet (const GridImp& g) : grid(g) {}

  /*
    We use the remove_const to extract the Type from the mutable class,
    because the const class is not instantiated yet.
  */
  template<int cd>
  IdType id (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const 
  {
          return grid.template getRealEntity<cd>(e).persistentIndex();
  }

  /*
    We use the remove_const to extract the Type from the mutable class,
    because the const class is not instantiated yet.
  */
  template<int cc>
  IdType subId (const typename remove_const<GridImp>::type::Traits::template Codim<0>::Entity& e, int i) const
  {
        return grid.template getRealEntity<0>(e).template subPersistentIndex<cc>(i);
  }

private:
  const GridImp& grid;
};


template<int dim, int dimworld>
struct SGridFamily
{
  typedef GridTraits<dim,dimworld,Dune::SGrid<dim,dimworld>,
                     SGeometry,SEntity,
                     SEntityPointer,SLevelIterator,
                     SIntersectionIterator, // leaf intersection
                     SIntersectionIterator, // level intersection
                     SIntersectionIterator, // leaf  intersection iter 
                     SIntersectionIterator, // level intersection iter
                     SHierarchicIterator,
                     SLevelIterator,
                                         SGridLevelIndexSet<const SGrid<dim,dimworld> >,
                                         SGridLevelIndexSetTypes<const SGrid<dim,dimworld> >,
                                         SGridLeafIndexSet<const SGrid<dim,dimworld> >,
                                         SGridLeafIndexSetTypes<const SGrid<dim,dimworld> >,
                                         SGridGlobalIdSet<const SGrid<dim,dimworld> >,
                                         bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits>,
                                         SGridGlobalIdSet<const SGrid<dim,dimworld> >,
                                         bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits>, 
                                         CollectiveCommunication<Dune::SGrid<dim,dimworld> > > 
  Traits;
};


//************************************************************************
template<int dim, int dimworld>
class SGrid : public GridDefaultImplementation <dim,dimworld,sgrid_ctype,SGridFamily<dim,dimworld> >
{
public:
  typedef SGridFamily<dim,dimworld> GridFamily;
  typedef bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits> PersistentIndexType;

  // need for friend declarations in entity
  typedef SGridLevelIndexSet<SGrid<dim,dimworld> > LevelIndexSetType;
  typedef SGridLeafIndexSet<SGrid<dim,dimworld> > LeafIndexSetType;
  typedef SGridGlobalIdSet<SGrid<dim,dimworld> > GlobalIdSetType;

  typedef typename SGridFamily<dim,dimworld>::Traits Traits;

  enum { MAXL=32 };

  typedef sgrid_ctype ctype;

  std::string name() const { return "SGrid"; };

  // constructors

  SGrid (const int* N_, const sgrid_ctype* H_);
  
  SGrid (const int* N_, const sgrid_ctype* L_, const sgrid_ctype* H_);

  SGrid (FieldVector<int,dim> N_, FieldVector<sgrid_ctype,dim> L_, FieldVector<sgrid_ctype,dim> H_);

  SGrid ();

  ~SGrid ();

  int maxLevel() const;

  template<int cd, PartitionIteratorType pitype>
  typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lbegin (int level) const;

  template<int cd, PartitionIteratorType pitype>
  typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lend (int level) const;

  template<int cd>