sgrid.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#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/bigunsignedint.hh>
#include <dune/common/parallel/collectivecommunication.hh>
#include <dune/common/reservedvector.hh>
#include <dune/geometry/genericgeometry/topologytypes.hh>
#include <dune/geometry/axisalignedcubegeometry.hh>
#include <dune/grid/common/capabilities.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, class GridImp> class SEntitySeed;
  template<int codim, PartitionIteratorType, class GridImp> class SLevelIterator;
  template<int dim, int dimworld, class ctype> class SGrid;
  template<class GridImp> class SIntersection;
  template<class GridImp> class SIntersectionIterator;
  template<class GridImp> class SHierarchicIterator;
 
  namespace FacadeOptions
  {
 
    template<int dim, int dimworld, class ctype, int mydim, int cdim>
    struct StoreGeometryReference<mydim, cdim,
        SGrid<dim,dimworld,ctype>, SGeometry>
    {
      static const bool v = false;
    };
 
    template<int dim, int dimworld, class ctype, int mydim, int cdim>
    struct StoreGeometryReference<mydim, cdim,
        const SGrid<dim,dimworld,ctype>, SGeometry>
    {
      static const bool v = false;
    };
 
  }
 
  //************************************************************************
  template<int mydim, int cdim, class GridImp>
  class SGeometry
    : public AxisAlignedCubeGeometry<typename GridImp::ctype,mydim,cdim>
  {
  public:
    typedef typename GridImp::ctype ctype;
 
    void make (const FieldVector<ctype,cdim>& lower,
               const FieldMatrix<ctype,mydim,cdim>& A)
    {
      if (mydim==0) {
        // set up base class
        static_cast< AxisAlignedCubeGeometry<ctype,mydim,cdim> & >( *this ) = AxisAlignedCubeGeometry<ctype,mydim,cdim>(lower);
        return;
      }
 
      // construct the upper right corner of the cube geometry
      FieldVector<ctype, cdim> upper = lower;
      for (int i=0; i<mydim; i++)
        upper += A[i];
 
      // look for the directions where the cube is actually extended
      std::bitset<cdim> axes;
 
      for (size_t i=0; i<cdim; i++)
        if ((upper[i] - lower[i]) > 1e-10)
          axes[i] = true;
 
      // set up base class
      static_cast< AxisAlignedCubeGeometry<ctype,mydim,cdim> & >( *this ) = AxisAlignedCubeGeometry<ctype,mydim,cdim>(lower, upper, axes);
    }
 
    SGeometry ()
      : AxisAlignedCubeGeometry<ctype,mydim,cdim>(FieldVector<ctype,cdim>(0),FieldVector<ctype,cdim>(0))    // anything
    {}
  };
 
 
  //************************************************************************
  template<int codim, int dim, class GridImp, template<int,int,class> class EntityImp>
  class SEntityBase :
    public EntityDefaultImplementation<codim,dim,GridImp,EntityImp>
  {
    friend class SEntityPointer<codim,GridImp>;
    friend class SIntersectionIterator<GridImp>;
    enum { dimworld = GridImp::dimensionworld };
 
    typedef typename GridImp::Traits::template Codim< codim >::GeometryImpl GeometryImpl;
 
  public:
    typedef typename GridImp::ctype ctype;
    typedef typename GridImp::template Codim<codim>::Geometry Geometry;
    typedef typename GridImp::PersistentIndexType PersistentIndexType;
 
    int level () const
    {
      return l;
    }
 
    int globalIndex() const;
 
    SEntitySeed<codim, GridImp> seed () const {
      return SEntitySeed<codim, GridImp>(l, index);
    }
 
    GeometryType type () const
    {
      static const GeometryType cubeType(GeometryType::cube,dim-codim);
      return cubeType;
    }
 
    Geometry geometry () const
    {
      if (!builtgeometry) makegeometry();
 
      // return result
      return Geometry( geo );
    }
 
    PartitionType partitionType () const { return InteriorEntity; }
 
    SEntityBase (GridImp* _grid, int _l, int _index) :
      grid(_grid),
      l(_l),
      index(_index),
      z(grid->z(l,index,codim)),
      builtgeometry(false) {}
 
    SEntityBase () :
      builtgeometry(false) // mark geometry as not built
    {}
 
    SEntityBase ( const SEntityBase& other ) :
      grid(other.grid),
      l(other.l),
      index(other.index),
      z(other.z),
      geo(), // do not copy geometry
      builtgeometry(false) // mark geometry as not built
    {}
 
    void make (GridImp* _grid, int _l, int _id);
 
    void make (int _l, int _id);
 
    void makegeometry () const;
 
    PersistentIndexType persistentIndex () const
    {
      return grid->persistentIndex(l, codim, z);
    }
 
    int compressedIndex () const
    {
      return index;
    }
 
    int compressedLeafIndex () const
    {
      // codim != dim -> there are no copies of entities
      // maxlevel -> ids are fine
      if (codim<dim || l==grid->maxLevel())
        return compressedIndex();
 
      // 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);
    }
 
    int subCompressedIndex (int cd, int i) const
    {
      DUNE_THROW(NotImplemented,"subIndex for entities with codimension > 0 is not implemented");
      return -1;
    }
 
    int subCompressedLeafIndex (int cd, int i) const
    {
      DUNE_THROW(NotImplemented,"subIndex for entities with codimension > 0 is not implemented");
      return -1;
    }
 
  protected:
    // this is how we implement our elements
    GridImp* grid;       
    int l;               
    int index;           
    array<int,dim> z;    
    mutable GeometryImpl geo; 
    mutable bool builtgeometry; 
  };
 
 
  template<int codim, int dim, class GridImp>
  class SEntity : public SEntityBase<codim,dim,GridImp,SEntity>
  {
    typedef Dune::SEntityBase<codim,dim,GridImp,Dune::SEntity> SEntityBase;
    friend class SEntityPointer<codim,GridImp>;
    friend class SIntersectionIterator<GridImp>;
  public:
    SEntity (GridImp* _grid, int _l, int _id) :
      SEntityBase(_grid,_l,_id) {}
  };
 
  template<int dim, class GridImp>
  class SEntity<0,dim,GridImp> : public SEntityBase<0,dim,GridImp,SEntity>
  {
    enum { dimworld = GridImp::dimensionworld };
    typedef Dune::SEntityBase<0,dim,GridImp,Dune::SEntity> SEntityBase;
    using SEntityBase::grid;
    using SEntityBase::l;
    using SEntityBase::index;
    using SEntityBase::z;
 
    typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
    typedef typename GridImp::Traits::template Codim< 0 >::LocalGeometryImpl LocalGeometryImpl;
 
    friend class SEntityPointer<0,GridImp>;
    friend class SIntersectionIterator<GridImp>;
 
  public:
    typedef typename GridImp::ctype ctype;
    typedef typename GridImp::template Codim<0>::Geometry Geometry;
    typedef typename GridImp::template Codim<0>::LocalGeometry LocalGeometry;
    template <int cd>
    struct Codim
    {
      typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
    };
    typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
    typedef typename GridImp::LeafIntersectionIterator IntersectionIterator;
    typedef typename GridImp::HierarchicIterator HierarchicIterator;
    typedef typename GridImp::PersistentIndexType PersistentIndexType;
 
    friend class SHierarchicIterator<GridImp>;
 
    template<int cc> int count () const;
 
    template<int cc> typename Codim<cc>::EntityPointer subEntity (int i) const;
 
    int subCompressedIndex (int codim, int i) const
    {
      if (codim==0) return this->compressedIndex();
      // compute subIndex
      return (this->grid)->n(this->l, this->grid->subz(this->z,i,codim));
    }
 
    int subCompressedLeafIndex (int codim, int i) const
    {
      if (codim==0) return this->compressedLeafIndex();
 
      assert(this->l == this->grid->maxLevel());
      // compute subIndex
      return (this->grid)->n(this->l, this->grid->subz(this->z,i,codim));
    }
 
    PersistentIndexType subPersistentIndex (int codim, int i) const
    {
      if (codim==0) return this->persistentIndex();
      // compute subId
      return this->grid->persistentIndex(this->l, codim, this->grid->subz(this->z,i,codim));
    }
 
    IntersectionIterator ibegin () const;
    IntersectionIterator ileafbegin () const;
    IntersectionIterator ilevelbegin () const;
    IntersectionIterator iend () const;
    IntersectionIterator ileafend () const;
    IntersectionIterator ilevelend () const;
 
    EntityPointer father () const;
 
    bool hasFather () const
    {
      return (this->level()>0);
    }
 
    bool isLeaf () const
    {
      return ( this->grid->maxLevel() == this->level() );
    }
 
    LocalGeometry geometryInFather () const;
 
    HierarchicIterator hbegin (int maxLevel) const;
 
    HierarchicIterator hend (int maxLevel) const;
 
    // members specific to SEntity
    SEntity (GridImp* _grid, int _l, int _index) :
      SEntityBase(_grid,_l,_index),
      built_father(false)
    {}
 
    SEntity (const SEntity& other ) :
      SEntityBase(other.grid, other.l, other.index ),
      built_father(false)
    {}
 
    void make (GridImp* _grid, int _l, int _id)
    {
      SEntityBase::make(_grid,_l,_id);
      built_father = false;
    }
 
    void make (int _l, int _id)
    {
      SEntityBase::make(_l,_id);
      built_father = false;
    }
 
  private:
 
    SEntity();
 
    mutable bool built_father;
    mutable int father_index;
    mutable LocalGeometryImpl in_father_local;
    void make_father() const;
  };
 
 
  //************************************************************************
  struct SHierarchicStackElem {
    int l;
    int index;
    SHierarchicStackElem () : l(-1), index(-1) {}
    SHierarchicStackElem (int _l, int _index) {l=_l; index=_index;}
    bool operator== (const SHierarchicStackElem& s) const {return !operator!=(s);}
    bool operator!= (const SHierarchicStackElem& s) const {return l!=s.l || index!=s.index;}
  };
 
  template<class GridImp>
  class SHierarchicIterator :
    public Dune::SEntityPointer <0,GridImp>
  {
    friend class SHierarchicIterator<const GridImp>;
    enum { dim = GridImp::dimension };
    enum { dimworld = GridImp::dimensionworld };
    typedef Dune::SEntityPointer<0,GridImp> SEntityPointer;
    using SEntityPointer::realEntity;
    using SEntityPointer::grid;
    using SEntityPointer::l;
    using SEntityPointer::index;
  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) :
      SEntityPointer(_grid,_e.level(),_e.compressedIndex())
    {
      // 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_index = _grid->getRealImplementation(this->entity()).compressedIndex();
 
      // push original element on stack
      SHierarchicStackElem originalElement(orig_l, orig_index);
      stack.push(originalElement);
 
      // compute maxLevel
      maxLevel = std::min(_maxLevel,this->grid->maxLevel());
 
      // ok, push all the sons as well
      push_sons(orig_l,orig_index);
 
      // and pop the first son
      increment();
    }
 
  private:
    int maxLevel;              
    int orig_l, orig_index;       
 
    std::stack<SHierarchicStackElem, Dune::ReservedVector<SHierarchicStackElem,GridImp::MAXL> > stack;
 
    void push_sons (int level, int fatherid); 
  };
 
  //************************************************************************
  template<class GridImp>
  class SIntersectionIterator
  {
    enum { dim=GridImp::dimension };
    enum { dimworld=GridImp::dimensionworld };
 
    typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl;
    typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl;
 
    friend class SIntersection<GridImp>;
 
  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::SIntersection<GridImp> IntersectionImp;
    typedef Dune::Intersection< const GridImp, Dune::SIntersection< const GridImp > > 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 intersection;
    }
 
    EntityPointer inside() const;
 
    EntityPointer outside() const;
 
    bool boundary () const;
 
    int boundaryId () const {
      if (boundary()) return count + 1;
      return 0;
    }
 
    int boundarySegmentIndex () const {
      if (boundary())
        return grid->boundarySegmentIndex(self.level(), count, zred);
      return -1;
    }
 
    bool neighbor () const;
 
    LocalGeometry geometryInInside () const;
    LocalGeometry geometryInOutside () const;
    Geometry geometry () const;
 
    GeometryType type () const
    {
      static const GeometryType cubeType(GeometryType::cube,dim-1);
      return cubeType;
    }
 
    int indexInInside () const;
    int indexInOutside () const;
 
    SIntersectionIterator (GridImp* _grid, const SEntity<0,dim,GridImp >* _self, int _count) :
      self(*_self), ne(self), grid(_grid),
      partition(_grid->partition(grid->getRealImplementation(ne).l,_self->z)),
      zred(_grid->compress(grid->getRealImplementation(ne).l,_self->z)),
      intersection(IntersectionImp(*this))
    {
      // make neighbor
      make(_count);
    }
 
    SIntersectionIterator (const SIntersectionIterator & other) :
      self(other.self), ne(other.ne), grid(other.grid),
      partition(other.partition), zred(other.zred),
      count(other.count), valid_count(other.valid_count),
      valid_nb(other.valid_nb), is_on_boundary(other.is_on_boundary),
      built_intersections(false),
      intersection(IntersectionImp(*this))
    {}
 
    SIntersectionIterator& operator = (const SIntersectionIterator& other)
    {
      /* We can't assign the grid */
      assert(grid == other.grid);
 
      /* Assign data from other */
      self = other.self;
      ne = other.ne;
      partition = other.partition;
      zred = other.zred;
      count = other.count;
      valid_count = other.valid_count;
      valid_nb = other.valid_nb;
      is_on_boundary = other.is_on_boundary;
 
      /* mark cached data as invalid */
      built_intersections = false;
 
      return *this;
    }
 
  private:
    void make (int _count) const;         
    void makeintersections () const;      
    EntityPointer self;                   
    mutable EntityPointer ne;             
    const GridImp * grid;                 
    int partition;                        
    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 LocalGeometryImpl is_self_local;  
    mutable GeometryImpl is_global;           
    mutable LocalGeometryImpl is_nb_local;    
    Intersection intersection;
  };
 
  template<class GridImp>
  class SIntersection
  {
    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 Geometry::LocalCoordinate LocalCoordinate;
    typedef typename Geometry::GlobalCoordinate GlobalCoordinate;
    typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
    typedef Dune::Intersection< const GridImp, Dune::SIntersectionIterator< const GridImp > > Intersection;
    enum { dimension=dim };
    enum { dimensionworld=dimworld };
    typedef typename GridImp::ctype ctype;
 
    bool boundary () const
    {
      return is.boundary();
    }
 
    int boundaryId () const
    {
      return is.boundaryId();
    }
 
    size_t boundarySegmentIndex () const
    {
      return is.boundarySegmentIndex();
    }
 
    bool neighbor () const
    {
      return is.neighbor();
    }
 
    EntityPointer inside() const
    {
      return is.inside();
    }
 
    EntityPointer outside() const
    {
      return is.outside();
    }
 
    bool conforming () const
    {
      return true;
    }
 
    LocalGeometry geometryInInside () const
    {
      return is.geometryInInside();
    }
 
    LocalGeometry geometryInOutside () const
    {
      return is.geometryInOutside();
    }
 
    Geometry geometry () const
    {
      return is.geometry();
    }
 
    GeometryType type () const
    {
      return is.type();
    }
 
    int indexInInside () const
    {
      return is.indexInInside();
    }
 
    int indexInOutside () const
    {
      return is.indexInOutside();
    }
 
    GlobalCoordinate outerNormal (const LocalCoordinate& local) const
    {
      return centerUnitOuterNormal();
    }
 
    GlobalCoordinate integrationOuterNormal (const LocalCoordinate& local) const
    {
      FieldVector<ctype, dimworld> n = centerUnitOuterNormal();
      n *= is.geometry().integrationElement(local);
      return n;
    }
 
    GlobalCoordinate unitOuterNormal (const LocalCoordinate& local) const
    {
      return centerUnitOuterNormal();
    }
 
    GlobalCoordinate centerUnitOuterNormal () const
    {
      FieldVector<ctype, dimworld> normal(0.0);
      normal[is.count/2] =  (is.count%2) ? 1.0 : -1.0;
      return normal;
    }
 
    SIntersection (const SIntersectionIterator<GridImp> & is_) : is(is_) {}
 
  private:
#ifndef DOXYGEN // doxygen can't handle this recursive usage
    const SIntersectionIterator<GridImp> & is;
#endif
  };
 
  //************************************************************************
 
  template <class T>
  class AutoPtrStack : public std::stack<T*>
  {
  public:
    ~AutoPtrStack()
    {
      while(! this->empty())
      {
        T* e = this->top();
        delete e;
        this->pop();
      }
    }
  };
 
  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 _index) :
      grid(_grid), l(_l), index(_index),
      e(0)
    {}
 
    SEntityPointer (const SEntity<codim,dim,GridImp> & _e) :
      grid(_e.grid), l(_e.l), index(_e.index),
      e(0)
    {}
 
    SEntityPointer (const SEntityPointer<codim,GridImp>& other) :
      grid(other.grid), l(other.l), index(other.index),
      e( 0 )
    {}
 
    ~SEntityPointer()
    {
      if( e )
        enStack().push( e );
#ifndef NDEBUG
      index = -1;
#endif
    }
 
    SEntityPointer& operator = (const SEntityPointer& other)
    {
      grid = other.grid;
      l = other.l;
      index = other.index;
 
      // free current entity
      if( e )
        enStack().push( e );
      e = 0;
 
      return *this;
    }
 
  protected:
    inline SEntity<codim,dim,GridImp>& realEntity() const
    {
      return grid->getRealImplementation(entity());
    }
 
    inline Entity& entity() const
    {
      if( ! e )
      {
        e = getEntity( grid, l, index );
      }
      return *e;
    }
 
    typedef AutoPtrStack< Entity > EntityStackType;
    static inline EntityStackType& enStack()
    {
      static EntityStackType eStack;
      return eStack;
    }
 
    inline Entity* getEntity(GridImp* _grid, int _l, int _id ) const
    {
      // get stack reference
      EntityStackType& enSt = enStack();
 
      if( enSt.empty() )
      {
        return (new Entity(SEntity<codim,dim,GridImp>(_grid, _l, _id)));
      }
      else
      {
        Entity* e = enSt.top();
        enSt.pop();
        grid->getRealImplementation(*e).make(_grid, _l,_id);
        return e;
      }
    }
 
    GridImp* grid;               
    int l;                       
    mutable int index;           
    mutable Entity* e;           
  };
 
  template<int codim, class GridImp>
  class SEntitySeed
  {
    enum { dim = GridImp::dimension };
  public:
    enum { codimension = codim };
 
    SEntitySeed () :
      _l(-1), _index(0)
    {}
 
    SEntitySeed (int l, int index) :
      _l(l), _index(index)
    {}
 
    bool isValid() const
    {
      return _l != -1;
    }
 
    int level () const { return this->_l; }
    int index () const { return this->_index; }
 
  private:
    int _l;                       
    int _index;                   
  };
 
  //************************************************************************
 
 
  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 };
    typedef Dune::SEntityPointer<codim,GridImp> SEntityPointer;
    using SEntityPointer::realEntity;
    using SEntityPointer::l;
    using SEntityPointer::index;
  public:
    typedef typename GridImp::template Codim<codim>::Entity Entity;
 
    void increment();
 
    SLevelIterator (GridImp * _grid, int _l, int _id) :
      SEntityPointer(_grid,_l,_id) {}
  };
 
 
  //========================================================================
  //========================================================================
 
  template<class GridImp>
  class SGridLevelIndexSet : public IndexSet<GridImp,SGridLevelIndexSet<GridImp> >
  {
    typedef SGridLevelIndexSet< GridImp > This;
    typedef IndexSet< GridImp, This > Base;
 
    enum { dim = GridImp::dimension };
 
  public:
 
    SGridLevelIndexSet ( const GridImp &g, int l )
      : grid( g ),
        level( l )
    {
      // TODO move list of geometrytypes to grid, can be computed static (singleton)
      // 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.getRealImplementation(e).compressedIndex();
    }
 
    template< int cc >
    int subIndex ( const typename GridImp::Traits::template Codim< cc >::Entity &e,
                   int i, unsigned int codim ) const
    {
      if( cc == 0 )
        return grid.getRealImplementation(e).subCompressedIndex(codim, i);
      else
        DUNE_THROW( NotImplemented, "subIndex for higher codimension entity not implemented for SGrid." );
    }
 
    // return true if the given entity is contained in \f$E\f$.
    template< class EntityType >
    bool contains ( const EntityType &e ) const
    {
      return (e.level() == level);
    }
 
    int size (GeometryType type) const
    {
      return grid.size( level, type );
    }
 
    int size (int codim) const
    {
      return grid.size( level, codim );
    }
 
    const std::vector<GeometryType>& geomTypes (int codim) const
    {
      return mytypes[codim];
    }
 
  private:
    const GridImp& grid;
    int level;
    std::vector<GeometryType> mytypes[GridImp::dimension+1];
  };
 
 
 
  //========================================================================
  //========================================================================
 
  template<class GridImp>
  class SGridGlobalIdSet :
    public IdSet<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.
     */
  {
    typedef SGridGlobalIdSet< GridImp > This;
 
  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;
 
    /*
       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 GridImp::getRealImplementation(e).persistentIndex();
    }
 
    /*
       We use the remove_const to extract the Type from the mutable class,
       because the const class is not instantiated yet.
     */
    IdType subId ( const typename remove_const< GridImp >::type::Traits::template Codim< 0 >::Entity &e,
                   int i, unsigned int codim ) const
    {
      return GridImp::getRealImplementation(e).subPersistentIndex(codim, i);
    }
  };
 
 
  template<int dim, int dimworld, class ctype>
  struct SGridFamily
  {
    typedef GridTraits<dim,dimworld,Dune::SGrid<dim,dimworld,ctype>,
        SGeometry,SEntity,
        SEntityPointer,SLevelIterator,
        SIntersection,              // leaf intersection
        SIntersection,              // level intersection
        SIntersectionIterator,              // leaf  intersection iter
        SIntersectionIterator,              // level intersection iter
        SHierarchicIterator,
        SLevelIterator,
        SGridLevelIndexSet<const SGrid<dim,dimworld,ctype> >,
        SGridLevelIndexSet<const SGrid<dim,dimworld,ctype> >,
        SGridGlobalIdSet<const SGrid<dim,dimworld,ctype> >,
        bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits>,
        SGridGlobalIdSet<const SGrid<dim,dimworld,ctype> >,
        bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits>,
        CollectiveCommunication<Dune::SGrid<dim,dimworld,ctype> >,
        DefaultLevelGridViewTraits, DefaultLeafGridViewTraits,
        SEntitySeed>
    Traits;
  };
 
 
  //************************************************************************
  template<int dim, int dimworld, typename _ctype = sgrid_ctype>
  class SGrid : public GridDefaultImplementation <dim,dimworld,_ctype,SGridFamily<dim,dimworld,_ctype> >
  {
  public:
    typedef SGridFamily<dim,dimworld,_ctype> 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 SGridLevelIndexSet<SGrid<dim,dimworld> > LeafIndexSetType;
    typedef SGridGlobalIdSet<SGrid<dim,dimworld> > GlobalIdSetType;
 
    typedef typename SGridFamily<dim,dimworld,_ctype>::Traits Traits;
 
    enum { MAXL=32 };
 
    typedef _ctype ctype;
 
    // constructors
 
    SGrid (const int * const N_, const ctype * const H_);
 
    SGrid (const int * const N_, const ctype * const L_, const ctype * const H_);
 
    SGrid (FieldVector<int,dim> N_, FieldVector<ctype,dim> L_, FieldVector<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>
    typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lbegin (int level) const
    {
      return lbegin<cd,All_Partition>(level);
    }
 
    template<int cd>
    typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lend (int level) const
    {
      return lend<cd,All_Partition>(level);
    }
 
    template<int cd, PartitionIteratorType pitype>
    typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafbegin () const;
 
    template<int cd, PartitionIteratorType pitype>
    typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafend () const;
 
    template<int cd>
    typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafbegin () const
    {
      return leafbegin<cd,All_Partition>();
    }
 
    template<int cd>
    typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafend () const
    {
      return leafend<cd,All_Partition>();
    }
 
    // \brief obtain EntityPointer from EntitySeed. */
    template <typename Seed>
    typename Traits::template Codim<Seed::codimension>::EntityPointer
    entityPointer(const Seed& seed) const
    {
      enum { codim = Seed::codimension };
      return SEntityPointer<codim,const SGrid<dim,dimworld> >(this,
                                                              this->getRealImplementation(seed).level(),
                                                              this->getRealImplementation(seed).index());
    }
 
    template<class T, template<class> class P, int codim>
    void communicate (T& t, InterfaceType iftype, CommunicationDirection dir, int level)
    {
      // SGrid is sequential and has no periodic boundaries, so do nothing ...
      return;
    }
 
    int size (int level, int codim) const;
 
    int size (int codim) const
    {
      return size(maxLevel(),codim);
    }
 
    int size (int level, GeometryType type) const
    {
      return (type.isCube()) ? size(level,dim-type.dim()) : 0;
    }
 
    int size (GeometryType type) const
    {
      return size(maxLevel(),type);
    }
 
    size_t numBoundarySegments () const
    {
      return boundarysize;
    }
 
    int global_size (int codim) const;
 
    int overlapSize (int level, int codim)
    {
      return 0;
    }
 
    int ghostSize (int level, int codim)
    {
      return 0;
    }
 
    // these are all members specific to sgrid
 
    void globalRefine (int refCount);
 
    const array<int, dim>& dims(int level) const {
      return N[level];
    }
 
    const FieldVector<ctype, dimworld>& lowerLeft() const {
      return low;
    }
 
    FieldVector<ctype, dimworld> upperRight() const {
      return H;
    }
 
    bool adapt ()
    {
      globalRefine(1);
      return true;
    }
 
    // The new index sets from DDM 11.07.2005
    const typename Traits::GlobalIdSet& globalIdSet() const
    {
      return theglobalidset;
    }
 
    const typename Traits::LocalIdSet& localIdSet() const
    {
      return theglobalidset;
    }
 
    const typename Traits::LevelIndexSet& levelIndexSet(int level) const
    {
      assert(level>=0 && level<=maxLevel());
      return *(indexsets[level]);
    }
 
    const typename Traits::LeafIndexSet& leafIndexSet() const
    {
      return *indexsets.back();
    }
 
    template<class DataHandle>
    void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level) const
    {}
 
    template<class DataHandle>
    void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir) const
    {}
 
    const CollectiveCommunication<SGrid>& comm () const
    {
      return ccobj;
    }
 
    int overlapSize (int level, int codim) const
    {
      return 0;
    }
 
    int overlapSize (int codim) const
    {
      return 0;
    }
 
    int ghostSize (int level, int codim) const
    {
      return 0;
    }
 
    int ghostSize (int codim) const
    {
      return 0;
    }
 
    /*
       @}
     */
 
  private:
    /*
       Make associated classes friends to grant access to the real entity
     */
    friend class Dune::SGridLevelIndexSet<Dune::SGrid<dim,dimworld> >;
    friend class Dune::SGridGlobalIdSet<Dune::SGrid<dim,dimworld> >;
    friend class Dune::SIntersectionIterator<Dune::SGrid<dim,dimworld> >;
    friend class Dune::SHierarchicIterator<Dune::SGrid<dim,dimworld> >;
    friend class Dune::SEntity<0,dim,Dune::SGrid<dim,dimworld> >;
 
    friend class Dune::SGridLevelIndexSet<const Dune::SGrid<dim,dimworld> >;
    friend class Dune::SGridGlobalIdSet<const Dune::SGrid<dim,dimworld> >;
    friend class Dune::SIntersectionIterator<const Dune::SGrid<dim,dimworld> >;
    friend class Dune::SHierarchicIterator<const Dune::SGrid<dim,dimworld> >;
    friend class Dune::SEntity<0,dim,const Dune::SGrid<dim,dimworld> >;
 
    template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
    friend class Dune::SEntityBase;
 
    template<int codim_, class GridImp_>
    friend class Dune::SEntityPointer;
 
    template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
    friend class Entity;
 
    FieldVector<ctype, dimworld> pos (int level, array<int,dim>& z) const;
 
    int calc_codim (int level, const array<int,dim>& z) const;
 
    int n (int level, const array<int,dim>& z) const;
 
    array<int,dim> z (int level, int i, int codim) const;
 
    array<int,dim> subz (const array<int,dim> & z, int i, int codim) const;
 
    array<int,dim> compress (int level, const array<int,dim>& z) const;
 
    array<int,dim> expand (int level, const array<int,dim>& r, int b) const;
 
    int partition (int level, const array<int,dim>& z) const;
 
    bool exists (int level, const array<int,dim>& zred) const;
 
    // compute boundary segment index for a given zentity and a face
    int boundarySegmentIndex (int l, int face, const array<int,dim> & zentity) const
    {
      array<int,dim-1> zface;
      int dir = face/2;
      int side = face%2;
      // compute z inside the global face
      for (int i=0; i<dir; i++) zface[i] = zentity[i]/(1<<l);
      for (int i=dir+1; i<dim; i++) zface[i-1] = zentity[i]/(1<<l);
      zface = boundarymapper[dir].expand(zface, 0);
      // compute index in the face
      int index = boundarymapper[dir].n(zface);
      // compute offset
      for (int i=0; i<dir; i++)
        index += 2*boundarymapper[i].elements(0);
      index += side*boundarymapper[dir].elements(0);
      return index;
    }
 
    // compute persistent index for a given zentity
    PersistentIndexType persistentIndex (int l, int codim, const array<int,dim> & zentity) 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(zentity[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 (zentity[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(zentity[i]>>trailing);
        }
 
        return id;
      }
    }
 
    // disable copy and assign
    SGrid(const SGrid &) {}
    SGrid & operator = (const SGrid &) { return *this; }
    // generate SGrid
    void makeSGrid (const array<int,dim>& N_, const FieldVector<ctype, dim>& L_, const FieldVector<ctype, dim>& H_);
 
    /*
       internal data
     */
    CollectiveCommunication<SGrid> ccobj;
 
    ReservedVector<SGridLevelIndexSet<const SGrid<dim,dimworld> >*, MAXL> indexsets;
    SGridGlobalIdSet<const SGrid<dim,dimworld> > theglobalidset;
 
    int L;                        // number of levels in hierarchic mesh 0<=level<L
    FieldVector<ctype, dim> low;  // lower left corner of the grid
    FieldVector<ctype, dim> H;    // length of cube per direction
    std::vector<array<int,dim> > N;            // number of elements per direction for each level
    std::vector<FieldVector<ctype, dim> > h;   // mesh size per direction for each level
    mutable CubeMapper<dim> *mapper; // a mapper for each level
 
    // boundary segement index set
    array<CubeMapper<dim-1>, dim> boundarymapper; // a mapper for each coarse grid face
    int boundarysize;
  };
 
  namespace Capabilities
  {
 
    template<int dim, int dimw>
    struct hasSingleGeometryType< SGrid<dim,dimw> >
    {
      static const bool v = true;
      static const unsigned int topologyId = GenericGeometry :: CubeTopology< dim > :: type :: id ;
    };
 
    template<int dim, int dimw>
    struct isCartesian< SGrid<dim,dimw> >
    {
      static const bool v = true;
    };
 
    template<int dim, int dimw, int cdim>
    struct hasEntity< SGrid<dim,dimw>, cdim>
    {
      static const bool v = true;
    };
 
    template<int dim, int dimw>
    struct isLevelwiseConforming< SGrid<dim,dimw> >
    {
      static const bool v = true;
    };
 
    template<int dim, int dimw>
    struct isLeafwiseConforming< SGrid<dim,dimw> >
    {
      static const bool v = true;
    };
 
  } // end namespace Capabilities
 
} // end namespace Dune
 
#include "sgrid/sgrid.cc"
 
#endif