• 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

agrid.hh

#ifndef DUNE_ALBERTAGRID_IMP_HH
#define DUNE_ALBERTAGRID_IMP_HH

#if HAVE_ALBERTA

#include <iostream>
#include <fstream>
#include <dune/common/deprecated.hh>

#include <vector>
#include <assert.h>
#include <algorithm>

// Dune includes
#include <dune/common/misc.hh>
#include <dune/common/interfaces.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/common/stdstreams.hh>

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

#include <dune/common/exceptions.hh>

#include <dune/grid/common/grid.hh>
#include <dune/grid/common/defaultindexsets.hh>
#include <dune/grid/common/sizecache.hh>
#include <dune/grid/common/intersectioniteratorwrapper.hh>
#include <dune/grid/common/defaultgridview.hh>

// stack for index management 
#include <dune/grid/common/indexstack.hh>

//- Local includes 
// some cpp defines and include of alberta.h
#include "albertaheader.hh"

// grape data io 
#include <dune/grid/utility/grapedataioformattypes.hh>

// IndexManager defined in indexstack.hh 
// 10000 is the size of the finite stack used by IndexStack
typedef Dune::IndexStack<int,100000> IndexManagerType;

//#define CALC_COORD
// some extra functions for handling the Albert Mesh  
#include "albertaextra.hh"

#include <dune/grid/albertagrid/exceptions.hh>

// contains a simple memory management for some componds of this grid 
#include "agmemory.hh"

namespace Dune
{
  // i.e. double or float 
  typedef ALBERTA REAL albertCtype;
}

#include "referencetopo.hh"
#include "indexsets.hh"
#include "geometry.hh"

namespace Dune 
{
  template<int codim, int dim, class GridImp> class AlbertaGridEntity;
  template<int codim, PartitionIteratorType pitype, class GridImp> class AlbertaGridTreeIterator;
  template<int codim, PartitionIteratorType pitype, class GridImp> class AlbertaGridLeafIterator;
  template<int cd, class GridImp> class AlbertaGridEntityPointer;

  template<class GridImp>         class AlbertaGridHierarchicIterator;
  template<class GridImp>         class AlbertaGridIntersectionIterator;
  template<int dim, int dimworld> class AlbertaGrid;
  template<int dim, int dimworld> class AlbertaGridHierarchicIndexSet;

  template <int codim, int dim, class GridImp> 
  struct SelectEntityImp
  {
    typedef AlbertaGridEntity<codim,dim,GridImp> EntityImp;
    typedef Dune::Entity<codim, dim, const GridImp, AlbertaGridEntity> Entity;
    typedef MakeableInterfaceObject<Entity> EntityObject;
  };

  class AlbertaMarkerVector 
  {
    friend class AlbertaGrid<2,2>;
    friend class AlbertaGrid<2,3>;
    friend class AlbertaGrid<3,3>;

    enum { vxBufferSize_ = 10000 };
  public:
    AlbertaMarkerVector (bool meLevel=true) : up2Date_(false), meLevel_(meLevel) {} ;

    bool vertexNotOnElement(const int elIndex, const int vertex) const;

    bool edgeNotOnElement(const int elIndex, const int edge) const;

    bool faceNotOnElement(const int elIndex, const int face) const;

    template <class GridType>
    void markNewVertices(GridType &grid, int level);
   
    template <class GridType>
    void markNewLeafVertices(GridType &grid);
  
    bool up2Date () const { return up2Date_; }

    void unsetUp2Date () { up2Date_ = false; }
 
    void print() const;

  private:
    typedef std::vector<int> ArrayType;
    
    // built in array to mark on which element a vertex is reached
    ArrayType vec_;
    ArrayType edgevec_;
    ArrayType facevec_;

    // number of vertices 
    int numVertex_;

    // true is vertex marker is up to date 
    bool up2Date_;
    bool meLevel_; 
  };

  //**********************************************************************
  //
  // --AlbertaGridEntity
  // --Entity
  //
  template<int cd, int dim, class GridImp>
  class AlbertaGridEntity : 
    public EntityDefaultImplementation <cd,dim,GridImp,AlbertaGridEntity>
  {
    enum { dimworld = GridImp::dimensionworld };
    friend class AlbertaGrid < dim , dimworld >;
    friend class AlbertaGridEntity < 0, dim, GridImp>;
    friend class AlbertaGridTreeIterator < cd, All_Partition,GridImp>;
    friend class AlbertaGridEntityPointer<cd,GridImp>;

    typedef AlbertaGridGeometry<dim-cd,dimworld,GridImp> GeometryImp;
  public:
    template <int cc>
    struct Codim
    {
      typedef typename GridImp::template Codim<cc>::EntityPointer EntityPointer;
    };

    typedef typename GridImp::template Codim<cd>::Entity Entity;
    typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
    typedef typename GridImp::template Codim<cd>::Geometry Geometry;
    typedef typename GridImp::template Codim<cd>::LevelIterator LevelIterator;

    int level () const;

    int boundaryId () const ;
  
    AlbertaGridEntity(const GridImp &grid, int level, 
                      ALBERTA TRAVERSE_STACK * travStack);
 
    AlbertaGridEntity(const GridImp &grid, int level, bool);

    AlbertaGridEntity(const AlbertaGridEntity & org);

    PartitionType partitionType() const; 

    const Geometry & geometry () const;

    GeometryType type () const;

    //***********************************************
    //  End of Interface methods 
    //***********************************************
    ALBERTA EL_INFO *getElInfo () const;
    ALBERTA EL *getElement () const;

    void removeElInfo(); 

    int getFEVnum () const;

    bool equals ( const AlbertaGridEntity<cd,dim,GridImp> & i) const;
    
    // dummy function, only needed for codim 0 
    bool leafIt () const { return false; }
    
    // methods for setting the infos from the albert mesh
    void setTraverseStack (ALBERTA TRAVERSE_STACK *travStack);
    void setElInfo (ALBERTA EL_INFO *elInfo, int face,
                    int edge, int vertex );
  
    // same as setElInfo just with a entity given 
    void setEntity (const AlbertaGridEntity<cd,dim,GridImp> & org);

    // set level of entity 
    void setLevel ( int newLevel ); 
    void setNewLevel ( int newLevel , bool ) { setLevel(level); }
  
    const GridImp& grid() const { return grid_; }

  private:
    const GeometryImp &geoImp () const
    {
      return GridImp :: getRealImplementation( geo_ );
    }

    GeometryImp &geoImp ()
    {
      return GridImp :: getRealImplementation( geo_ );
    }

  private: 
    // the grid this entity belong to 
    const GridImp &grid_;

    // Alberta element info
    ALBERTA EL_INFO *elInfo_;

    // Alberta element 
    ALBERTA EL * element_;
  
    // current traverse stack this entity belongs too
    ALBERTA TRAVERSE_STACK * travStack_;

    int level_;

    // type of createable object, just derived from Geometry class
    typedef MakeableInterfaceObject<Geometry> GeometryObject;
      
    GeometryObject geo_;
  
    mutable bool builtgeometry_;  

    mutable FieldVector<albertCtype, dim> localFatherCoords_; 
    mutable bool localFCoordCalced_;

    int face_;

    int edge_;
  
    int vertex_;
  };

  //***********************
  //  
  //  --AlbertaGridEntity
  //  --0Entity
  //
  //***********************
  template<int dim, class GridImp>
  class AlbertaGridEntity<0,dim,GridImp> : 
    public EntityDefaultImplementation <0,dim,GridImp,AlbertaGridEntity>
  {
    enum { dimworld = GridImp::dimensionworld };
    friend class AlbertaGrid < dim , GridImp::dimensionworld >;
    friend class AlbertaMarkerVector;
    friend class AlbertaGridIntersectionIterator <GridImp>;
    friend class AlbertaGridHierarchicIterator <GridImp>;
    friend class AlbertaGridTreeIterator <0,All_Partition,GridImp>;
    friend class AlbertaGridEntityPointer<0,GridImp>;
  public:
    template <int cd>
    struct Codim
    {
      typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
    };

    typedef typename GridImp::template Codim<0>::Entity Entity;
    typedef typename GridImp::template Codim<0>::Geometry Geometry;
    typedef AlbertaGridGeometry<dim,dimworld,GridImp> GeometryImp;

    typedef typename GridImp::template Codim<0>::LevelIterator LevelIterator;
    typedef typename GridImp::template Codim<0>::HierarchicIterator HierarchicIterator;
    typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;

    typedef LeafIntersectionIteratorWrapper<GridImp> AlbertaGridLeafIntersectionIteratorType;
    typedef AlbertaGridLeafIntersectionIteratorType  AlbertaGridIntersectionIteratorType;
    typedef AlbertaGridLeafIntersectionIteratorType  AlbertaGridLevelIntersectionIteratorType;

    enum { dimension = dim };
  
    ~AlbertaGridEntity() {};

    AlbertaGridEntity(const GridImp &grid, int level, bool leafIt );

    AlbertaGridEntity(const AlbertaGridEntity & org);
    
    int level () const;

    int boundaryId () const;

    const Geometry & geometry () const;

    GeometryType type () const;


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

    //template<int cc> void entity (AlbertaGridTreeIterator<cc,dim,dimworld>& it, int i);

    AlbertaGridLeafIntersectionIteratorType ileafbegin () const;
    AlbertaGridIntersectionIteratorType ibegin () const
    {
      return ileafbegin();
    }

    AlbertaGridLevelIntersectionIteratorType ilevelbegin () const
    {
      DUNE_THROW(NotImplemented,"method ilevelbegin not implemented!"); 
      return ileafbegin();
    }

    AlbertaGridIntersectionIteratorType ileafend () const;
    AlbertaGridIntersectionIteratorType iend () const{ return ileafend();}
    AlbertaGridLeafIntersectionIteratorType ilevelend () const
    {
      DUNE_THROW(NotImplemented,"method ilevelend not implemented!"); 
      return ibegin();
    }

    bool isLeaf () const ; 

    EntityPointer father () const;

    const Geometry & geometryInFather () const;
  
    AlbertaGridHierarchicIterator<GridImp> hbegin (int maxlevel) const;
  
    AlbertaGridHierarchicIterator<GridImp> hend (int maxlevel) const;

      bool wasRefined () const DUNE_DEPRECATED {return isNew();}

    bool mightBeCoarsened () const DUNE_DEPRECATED { return mightVanish(); }

    bool isNew () const ;

    bool mightVanish () const ;

    bool hasBoundaryIntersections () const ;
    
    PartitionType partitionType() const; 

    bool equals ( const AlbertaGridEntity<0,dim,GridImp> & i) const;

    //***************************************************************
    //  Interface for parallelisation 
    //***************************************************************
    // set leaf data with processor number  
    void setLeafData( int proc );

    // return true if this entity belong to master set of this grid   
    bool master() const; 

    // return 0 for elements 
    int getFEVnum () const { return 0; } 

    // needed for LevelIterator to compare 
    ALBERTA EL_INFO *getElInfo () const;

    // return element for equaltiy in EntityPointer 
    ALBERTA EL *getElement () const;

    // set elInfo and Element to nil
    void removeElInfo();
    
    // returns true if entity comes from LeafIterator 
    bool leafIt () const { return leafIt_; }
    
    // face, edge and vertex only for codim > 0, in this
    // case just to supply the same interface
    void setTraverseStack (ALBERTA TRAVERSE_STACK *travStack);
    void setElInfo (ALBERTA EL_INFO *elInfo,
                    int face = 0,
                    int edge = 0,
                    int vertex = 0 );

    // same as setElInfo just with a entity given 
    void setEntity (const AlbertaGridEntity<0,dim,GridImp> & org);

    const GridImp& grid() const { return grid_; }
  private: 
    int nChild () const;
    
    void makeDescription();

    const GridImp & grid_;

    int level_;

    ALBERTA TRAVERSE_STACK * travStack_;

    ALBERTA EL_INFO *elInfo_;
  
    ALBERTA EL *element_;

    // local coordinates within father 
    typedef MakeableInterfaceObject<Geometry> GeometryObject;

    mutable GeometryObject geoObj_;
    mutable GeometryImp & geo_;
    mutable bool builtgeometry_;  

    // is true if entity comes from leaf iterator 
    bool leafIt_; 
  
  }; // end of AlbertaGridEntity codim = 0


  //**********************************************************************
  //
  // --AlbertaGridEntityPointer
  // --EntityPointer
  // --EnPointer
  template<int cd, class GridImp>
  class AlbertaGridEntityPointer
  {
    enum { dim       = GridImp::dimension };
    enum { dimworld  = GridImp::dimensionworld };

    friend class AlbertaGridEntity<cd,dim,GridImp>;
    friend class AlbertaGridEntity< 0,dim,GridImp>;
    friend class AlbertaGrid < dim , dimworld >;

  public:
    enum { codimension = cd };

    typedef typename GridImp::template Codim<cd>::Entity Entity;
    typedef typename SelectEntityImp<cd,dim,GridImp>::EntityImp EntityImp;
    typedef typename SelectEntityImp<cd,dim,GridImp>::EntityObject EntityObject;
    typedef AlbertaGridEntityPointer<cd,GridImp> EntityPointerImp;

    typedef AlbertaGridEntityPointer<cd,GridImp> AlbertaGridEntityPointerType;

    AlbertaGridEntityPointer(const GridImp & grid, 
                             int level,  ALBERTA EL_INFO *elInfo,int face,int edge,int vertex);

    AlbertaGridEntityPointer(const GridImp & grid, ALBERTA TRAVERSE_STACK * stack,
                             int level,  ALBERTA EL_INFO *elInfo,int face,int edge,int vertex);

    AlbertaGridEntityPointer(const GridImp & grid, int level , bool isLeaf, bool done);

    AlbertaGridEntityPointer(const EntityImp& entity);

    AlbertaGridEntityPointer(const AlbertaGridEntityPointerType & org);

    AlbertaGridEntityPointer(const GridImp & , const EntityImp & en);

    AlbertaGridEntityPointer& operator= (const AlbertaGridEntityPointer& org);

    ~AlbertaGridEntityPointer();

    bool equals (const AlbertaGridEntityPointerType& i) const;

    Entity & dereference () const ;

    int level () const ;

    void done ();

    void compactify() {}

  protected:
    bool leafIt () const { return isLeaf_; }

    EntityImp & entityImp (); 

    const EntityImp & entityImp () const; 
    
    // reference to grid 
    const GridImp & grid_;

    bool isLeaf_;

    // entity that this EntityPointer points to 
    EntityObject * entity_;
  };


  //**********************************************************************
  //
  // --AlbertaGridHierarchicIterator
  // --HierarchicIterator
  template<class GridImp>
  class AlbertaGridHierarchicIterator
  : public AlbertaGridEntityPointer<0,GridImp>
  {
  public:
    typedef typename GridImp::template Codim<0>::Entity Entity;
    typedef typename GridImp::ctype ctype;

    typedef typename SelectEntityImp<0,GridImp::dimension,GridImp>::EntityImp EntityImp;

    AlbertaGridHierarchicIterator(const GridImp &grid,
                                  ALBERTA TRAVERSE_STACK *travStack, 
                                  int actLevel, int maxLevel, bool leafIt );

    AlbertaGridHierarchicIterator(const GridImp &grid,
                                  int actLevel,int maxLevel);
  
    AlbertaGridHierarchicIterator(const AlbertaGridHierarchicIterator<GridImp> &org);
 
    AlbertaGridHierarchicIterator<GridImp> & operator = (const AlbertaGridHierarchicIterator<GridImp> &org);
 
    void increment();

  private:
    const int startLevel_;
    
    int level_;

    int maxlevel_;

    EntityImp & virtualEntity_;

    ALBERTA ManageTravStack manageStack_;

    bool end_;

    ALBERTA EL_INFO * recursiveTraverse(ALBERTA TRAVERSE_STACK * stack);

    ALBERTA EL_INFO * firstChild(ALBERTA TRAVERSE_STACK * stack);

    void makeIterator();
  };


  //**********************************************************************
  //
  // --AlbertaGridIntersectionIterator
  // --IntersectionIterator
  template<class GridImp>
  class AlbertaGridIntersectionIterator
  {
    enum { dim      = GridImp::dimension };
    enum { dimworld = GridImp::dimensionworld };

    friend class AlbertaGridEntity<0,dim,GridImp>;
    typedef AlbertaGridIntersectionIterator<GridImp> ThisType; 

  public:
    typedef Dune::Intersection< GridImp, Dune::AlbertaGridIntersectionIterator >
      Intersection;
    typedef ThisType ImplementationType;

    typedef AGMemoryProvider< ThisType > StorageType; 
    typedef typename GridImp::template Codim<0>::Entity Entity;
    typedef typename GridImp::template Codim<1>::Geometry Geometry;
    typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
    
    typedef typename SelectEntityImp<0,dim,GridImp>::EntityImp EntityImp;
    
    //typedef AlbertaGridMakeableGeometry<dim-1,dimworld,GridImp> LocalGeometryImp;
    typedef AlbertaGridGeometry<dim-1,dimworld,GridImp> LocalGeometryImp;
    typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;

    enum { dimension=dim };
    enum { dimensionworld=dimworld };
    typedef typename GridImp::ctype ctype;

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

    bool equals (const AlbertaGridIntersectionIterator<GridImp> & i) const;

    void increment();

    bool operator==(const AlbertaGridIntersectionIterator<GridImp>& i) const;

    EntityPointer outside() const;

    EntityPointer inside() const;

    AlbertaGridIntersectionIterator(const GridImp & grid,
                                    int level);

    AlbertaGridIntersectionIterator(const GridImp & grid,
                                    int level,
                                    ALBERTA EL_INFO *elInfo,
                                    bool leafIt );
    AlbertaGridIntersectionIterator(const AlbertaGridIntersectionIterator<GridImp> & org);

    void assign (const AlbertaGridIntersectionIterator<GridImp> & org);
  
    //~AlbertaGridIntersectionIterator();

    bool boundary () const;

    bool neighbor () const; 

    int boundaryId () const; 
      
    bool conforming () const;

    const LocalGeometry& intersectionSelfLocal () const;
    const LocalGeometry& intersectionNeighborLocal () const;
    const Geometry& intersectionGlobal () const;

    int numberInSelf () const;
    int numberInNeighbor () const;

    int twistInSelf() const;
    
    int twistInNeighbor() const;

    typedef FieldVector<albertCtype, GridImp::dimensionworld> NormalVecType;
    typedef FieldVector<albertCtype, GridImp::dimension-1> LocalCoordType;
  
    const NormalVecType & unitOuterNormal (const LocalCoordType & local) const;
  
    const NormalVecType & outerNormal (const LocalCoordType & local) const;

    const NormalVecType & integrationOuterNormal (const LocalCoordType & local) const;

    int level () const; 
    
    //**********************************************************
    //  private methods 
    //**********************************************************
    
    // reset IntersectionIterator  
    template <class EntityType> 
    void first(const EntityType & en, int level ); 

    // calls EntityPointer done and sets done_ to true
    void done (); 

  private:
    // returns true if actual neighbor has same level 
    bool neighborHasSameLevel () const;
    
    void makeBegin (const GridImp & grid,
                    int level, 
                    ALBERTA EL_INFO * elInfo ) const;
  
    void makeEnd (const GridImp & grid,int level ) const;

    // put objects on stack 
    void freeObjects () const;
  
    void setupVirtEn () const;

    void calcOuterNormal () const;

    // return whether the iterator was called from a LeafIterator entity or
    // LevelIterator entity
    bool leafIt () const { return leafIt_; }
    // private member variables 

    const GridImp& grid_;

    mutable int level_;
  
    mutable int neighborCount_;

    mutable bool builtNeigh_;

    bool leafIt_;

    mutable ALBERTA EL_INFO * elInfo_;

    typedef MakeableInterfaceObject<LocalGeometry> LocalGeometryObject;
   
    // the objects holding the real implementations 
    mutable LocalGeometryObject fakeNeighObj_;
    mutable LocalGeometryObject fakeSelfObj_;
    mutable LocalGeometryObject neighGlobObj_;
    
    mutable LocalGeometryImp & fakeNeigh_;
    mutable LocalGeometryImp & fakeSelf_;
    mutable LocalGeometryImp & neighGlob_;

    mutable ALBERTA EL_INFO neighElInfo_;

    mutable NormalVecType outNormal_;
    mutable NormalVecType unitNormal_;
  
    // tmp memory for normal calculation 
    mutable FieldVector<albertCtype, dimworld> tmpU_;
    mutable FieldVector<albertCtype, dimworld> tmpV_;

    // twist seen from the neighbor
    mutable int twist_;

    bool done_;
  };



  

  //**********************************************************************
  //
  // --AlbertaGridTreeIterator
  // --LevelIterator
  // --TreeIterator
  //
  
  namespace AlbertaTreeIteratorHelp {
    template <class IteratorImp, int dim, int codim>
    struct GoNextEntity;
  }
    
  template<int cd, PartitionIteratorType pitype, class GridImp>
  class AlbertaGridTreeIterator : 
    public AlbertaGridEntityPointer<cd,GridImp> 
  {
    enum { dim = GridImp::dimension };
    friend class AlbertaGridEntity<2,dim,GridImp>;
    friend class AlbertaGridEntity<1,dim,GridImp>;
    friend class AlbertaGridEntity<0,dim,GridImp>;
    friend class AlbertaGrid < dim , GridImp::dimensionworld >;


    typedef AlbertaGridTreeIterator<cd,pitype,GridImp>  AlbertaGridTreeIteratorType;
    typedef AlbertaGridTreeIteratorType ThisType;
    friend class AlbertaTreeIteratorHelp::GoNextEntity<ThisType,dim,cd>;
  public:
  
    typedef typename GridImp::template Codim<cd>::Entity Entity;
    typedef typename SelectEntityImp<cd,dim,GridImp>::EntityImp EntityImp;
  
    AlbertaGridTreeIterator(const AlbertaGridTreeIterator<cd,pitype,GridImp> & org ); 
  
    AlbertaGridTreeIterator<cd,pitype,GridImp> & operator = (const AlbertaGridTreeIterator<cd,pitype,GridImp> & org ); 
  
    AlbertaGridTreeIterator(const GridImp & grid, int
                            travLevel, int proc, bool leafIt=false );
 
    AlbertaGridTreeIterator(const GridImp & grid, 
                            const AlbertaMarkerVector * vec,
                            int travLevel, 
                            int proc,
                            bool leafIt=false);
  
    void increment();

  private:
    // private Methods
    void makeIterator();

    ALBERTA EL_INFO * goFirstElement(ALBERTA TRAVERSE_STACK *stack,
                                     ALBERTA MESH *mesh,
                                     int level, ALBERTA FLAGS fill_flag);
    ALBERTA EL_INFO * traverseElLevel(ALBERTA TRAVERSE_STACK * stack);
    ALBERTA EL_INFO * traverseElLevelInteriorBorder(ALBERTA TRAVERSE_STACK * stack);
    ALBERTA EL_INFO * traverseElLevelGhosts(ALBERTA TRAVERSE_STACK * stack);
 
    // the default is, go to next elInfo
    //template <int cc>  
    ALBERTA EL_INFO * goNextEntity(ALBERTA TRAVERSE_STACK *stack,ALBERTA EL_INFO *elinfo_old);
  
    // the real go next methods
    ALBERTA EL_INFO * goNextElInfo(ALBERTA TRAVERSE_STACK *stack,ALBERTA EL_INFO *elinfo_old);
    ALBERTA EL_INFO * goNextFace(ALBERTA TRAVERSE_STACK *stack,
                                 ALBERTA EL_INFO *elInfo);
    ALBERTA EL_INFO * goNextEdge(ALBERTA TRAVERSE_STACK *stack,
                                 ALBERTA EL_INFO *elInfo);
    ALBERTA EL_INFO * goNextVertex(ALBERTA TRAVERSE_STACK *stack,
                                   ALBERTA EL_INFO *elInfo);

    // search next macro el  
    ALBERTA MACRO_EL * nextGhostMacro(ALBERTA MACRO_EL *mel);

    int level_;
  
    int enLevel_;

    EntityImp & virtualEntity_; 
    
    // contains ALBERTA traversal stack 
    ALBERTA ManageTravStack manageStack_;

    int face_;
    int edge_;
    int vertex_;

    // knows on which element a point,edge,face is viewed
    const AlbertaMarkerVector * vertexMarker_;

    // variable for operator++
    bool okReturn_;

    // store processor number of elements
    // for ghost walktrough, i.e. walk over ghosts which belong
    // tp processor 2 
    const int proc_;
  };

  template<int cd, PartitionIteratorType pitype, class GridImp>
  class AlbertaGridLevelIterator
  : public AlbertaGridTreeIterator<cd,pitype,GridImp>
  {
  public:  
    typedef typename GridImp::template Codim<cd>::Entity Entity;

    AlbertaGridLevelIterator(const GridImp & grid, int level, int proc) :
      AlbertaGridTreeIterator<cd,pitype,GridImp> (grid,level,proc) 
    {
    }

    AlbertaGridLevelIterator(const GridImp & grid, 
                             const AlbertaMarkerVector * vec, int level, int proc) : 
      AlbertaGridTreeIterator<cd,pitype,GridImp> (grid,vec,level,proc) 
    {
    }

    void increment () 
    {
      AlbertaGridTreeIterator<cd,pitype,GridImp>::increment();
    }
  };

  //**********************************************************************
  //
  //  AlbertaGridLeafIterator 
  //  --LeafIterator 
  //
  //**********************************************************************
  template<int codim, PartitionIteratorType pitype, class GridImp>
  class AlbertaGridLeafIterator
  : public AlbertaGridTreeIterator<codim, pitype, GridImp>
  {
  public:  
    typedef typename GridImp::template Codim<codim>::Entity Entity;

    AlbertaGridLeafIterator(const GridImp & grid, int level, int proc) :
      AlbertaGridTreeIterator<codim,pitype,GridImp> (grid,level,proc,true) 
    {
    }

    AlbertaGridLeafIterator(const GridImp & grid, 
                            const AlbertaMarkerVector * vec, int level, int proc) : 
      AlbertaGridTreeIterator<codim, pitype, GridImp> (grid,vec,level,proc,true) 
    {
    }

    void increment () 
    {
      AlbertaGridTreeIterator<codim, pitype, GridImp>::increment();
    }
  };

  //**********************************************************************
  //
  // --AlbertaGrid
  // --Grid
  //
  //**********************************************************************
  
  template <int dim, int dimworld> 
  struct AlbertaGridFamily 
  {
    typedef AlbertaGrid<dim,dimworld> GridImp;

    typedef DefaultLevelIndexSet< AlbertaGrid<dim,dimworld> > LevelIndexSetImp;
    typedef DefaultLeafIndexSet< AlbertaGrid<dim,dimworld> > LeafIndexSetImp;

    typedef AlbertaGridIdSet<dim,dimworld> IdSetImp; 
    typedef int IdType;

    struct Traits
    {
      typedef GridImp Grid;

      typedef Dune :: Intersection< const GridImp, LeafIntersectionIteratorWrapper > LeafIntersection;
      typedef Dune :: Intersection< const GridImp, LeafIntersectionIteratorWrapper > LevelIntersection;
      typedef Dune::IntersectionIterator<const GridImp, LeafIntersectionIteratorWrapper, LeafIntersectionIteratorWrapper > LeafIntersectionIterator;
      typedef Dune::IntersectionIterator<const GridImp, LeafIntersectionIteratorWrapper, LeafIntersectionIteratorWrapper > LevelIntersectionIterator;

      typedef Dune::HierarchicIterator<const GridImp, AlbertaGridHierarchicIterator> HierarchicIterator;

      typedef IdType GlobalIdType;
      typedef IdType LocalIdType;

      template <int cd>
      struct Codim
      {
        // IMPORTANT: Codim<codim>::Geometry == Geometry<dim-codim,dimw>
        typedef Dune::Geometry<dim-cd, dimworld, const GridImp, AlbertaGridGeometry> Geometry;
        typedef Dune::Geometry<dim-cd, dim, const GridImp, AlbertaGridGeometry> LocalGeometry;
        // we could - if needed - introduce an other struct for dimglobal of Geometry

        typedef typename SelectEntityImp<cd,dim,GridImp>::Entity Entity; 

        typedef Dune::LevelIterator<cd,All_Partition,const GridImp,AlbertaGridLevelIterator> LevelIterator;

        typedef Dune::LeafIterator<cd,All_Partition,const GridImp,AlbertaGridLeafIterator> LeafIterator;

        typedef AlbertaGridEntityPointer< cd, const GridImp > EntityPointerImpl;
        typedef Dune::EntityPointer< const GridImp, EntityPointerImpl > EntityPointer;

        template <PartitionIteratorType pitype>
        struct Partition
        {
          typedef Dune::LevelIterator<cd,pitype,const GridImp,AlbertaGridLevelIterator> LevelIterator;
          typedef Dune::LeafIterator<cd,pitype,const GridImp,AlbertaGridLeafIterator> LeafIterator;
        };

      };

      template <PartitionIteratorType pitype>
      struct Partition
      {
        typedef Dune::GridView<DefaultLevelGridViewTraits<const GridImp,pitype> > 
        LevelGridView;
        typedef Dune::GridView<DefaultLeafGridViewTraits<const GridImp,pitype> > 
        LeafGridView;
      };

      typedef IndexSet<GridImp,LevelIndexSetImp,DefaultLevelIteratorTypes<GridImp> > LevelIndexSet;
      typedef IndexSet<GridImp,LeafIndexSetImp,DefaultLeafIteratorTypes<GridImp> > LeafIndexSet;
      typedef IdSet<GridImp,IdSetImp,IdType> GlobalIdSet;
      typedef IdSet<GridImp,IdSetImp,IdType> LocalIdSet;

//#if HAVE_MPI 
      // use collective communciation with MPI
//      typedef CollectiveCommunication<MPI_Comm> CollectiveCommunication;
//#else 
      // use dummy collective communication 
      typedef Dune :: CollectiveCommunication< GridImp >
        CollectiveCommunication;
//#endif
    };
  };

  template <int dim, int dimworld>
  class AlbertaGrid : 
    public GridDefaultImplementation <dim,dimworld,albertCtype, AlbertaGridFamily<dim,dimworld> >,
    public HasObjectStream , 
    public HasHierarchicIndexSet 
  {
    friend class AlbertaGridEntity <0,dim,const AlbertaGrid<dim,dimworld> >;
    friend class AlbertaGridEntity <1,dim,const AlbertaGrid<dim,dimworld> >;
    friend class AlbertaGridEntity <2,dim,const AlbertaGrid<dim,dimworld> >;
    friend class AlbertaGridEntity <dim,dim,const AlbertaGrid<dim,dimworld> >;

    friend class AlbertaGridEntityPointer <0,const AlbertaGrid<dim,dimworld> >;
    friend class AlbertaGridEntityPointer <1,const AlbertaGrid<dim,dimworld> >;
    friend class AlbertaGridEntityPointer <2,const AlbertaGrid<dim,dimworld> >;
    friend class AlbertaGridEntityPointer <3,const AlbertaGrid<dim,dimworld> >;

    // friends because of fillElInfo
    friend class AlbertaGridTreeIterator<0,All_Partition,AlbertaGrid<dim,dimworld> >;
    friend class AlbertaGridTreeIterator<1,All_Partition,AlbertaGrid<dim,dimworld> >;
    friend class AlbertaGridTreeIterator<2,All_Partition,AlbertaGrid<dim,dimworld> >;
    friend class AlbertaGridTreeIterator<3,All_Partition,AlbertaGrid<dim,dimworld> >;

    friend class AlbertaGridHierarchicIterator<AlbertaGrid<dim,dimworld> >;
  
    friend class AlbertaGridIntersectionIterator<AlbertaGrid<dim,dimworld> >;
    friend class AlbertaGridIntersectionIterator<const AlbertaGrid<dim,dimworld> >;

    //CompileTimeChecker<dimworld != 1>   Do_not_use_AlbertaGrid_for_1d_Grids;   
 
    typedef AlbertaGrid<dim,dimworld> MyType; 

    friend class AlbertaMarkerVector;
    friend class AlbertaGridHierarchicIndexSet<dim,dimworld>;

    // minimum number of elements assumed to be created during adaption 
    enum { defaultElementChunk_ = 100 };

    //**********************************************************
    // The Interface Methods
    //**********************************************************
  public: 
    typedef AlbertaGridFamily<dim,dimworld> GridFamily;  
    typedef GridDefaultImplementation <dim,dimworld,albertCtype,
            AlbertaGridFamily<dim,dimworld> > BaseType;

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

    typedef AlbertaGridHierarchicIndexSet<dim,dimworld> HierarchicIndexSet;

  private:  
    typedef typename Traits:: CollectiveCommunication  CollectiveCommunicationType;

    typedef typename Traits::template Codim<0>::LeafIterator LeafIterator;

    typedef typename GridFamily:: LevelIndexSetImp LevelIndexSetImp;
    typedef typename GridFamily:: LeafIndexSetImp LeafIndexSetImp;
    
    typedef typename Traits :: LeafIndexSet LeafIndexSet;

    typedef AlbertaGridIdSet<dim,dimworld> IdSetImp; 
    typedef typename Traits :: GlobalIdSet GlobalIdSet; 
    typedef typename Traits :: LocalIdSet  LocalIdSet; 

  public:  
    typedef typename ALBERTA AlbertHelp::AlbertLeafData<dimworld,dim+1> LeafDataType;

    class ObjectStream
    {
      public:
        class EOFException {} ;
        template <class T>
        void readObject (T &) {}
        void readObject (int) {}
        void readObject (double) {}
        template <class T>
        void writeObject (T &) {}
        void writeObject (int) {}
        void writeObject (double) {}

        template <class T>
        void read (T &) const {}
        template <class T>
        void write (const T &) {}
    };

    typedef ObjectStream ObjectStreamType;

    enum { 
      numCodim = dim+1 };

    enum { 
      MAXL = 64 };

    /*
       levInd = true means that a consecutive level index is generated 
       if levInd == true the the element number of first macro element is
       set to 1 so hasLevelIndex_ can be identified we grid is read from
       file */
    AlbertaGrid(const std::string macroTriangFilename);
  
    /* (for internal use only) 
       Constructor which reads an ALBERTA macro triangulation file 
       or given GridFile , proc is the number of domain , 
       levInd = true means that a consecutive level index is generated 
       if levInd == true the the element number of first macro element is
       set to 1 so hasLevelIndex_ can be identified we grid is read from
       file */
    //AlbertaGrid(AlbertaGrid<dim,dimworld> & oldGrid, int proc);
  
    AlbertaGrid();
 
    ~AlbertaGrid();
   
    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;

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

    template <int codim>
    typename Traits::template Codim<codim>::LeafIterator 
    leafbegin () const;
  
    template <int codim, PartitionIteratorType pitype>
    typename Traits::template Codim<codim>::template Partition<pitype>::LeafIterator 
    leafend   () const;

    template <int codim>
    typename Traits::template Codim<codim>::LeafIterator 
    leafend   () const;

  private:  
    template <int codim, PartitionIteratorType pitype>
    typename Traits::template Codim<codim>::template Partition<pitype>::LeafIterator 
    leafbegin ( int maxlevel, int proc = -1 ) const;

    template <int codim>
    typename Traits::template Codim<codim>::LeafIterator 
    leafbegin ( int maxlevel, int proc = -1 ) const;
  
    template <int codim, PartitionIteratorType pitype>
    typename Traits::template Codim<codim>::template Partition<pitype>::LeafIterator 
    leafend   ( int maxlevel, int proc = -1 ) const;

    template <int codim>
    typename Traits::template Codim<codim>::LeafIterator 
    leafend   ( int maxlevel, int proc = -1 ) const;

    LeafIterator leafbegin ( int maxlevel, int proc = -1 ) const;
  
    LeafIterator leafend   ( int maxlevel, int proc = -1 ) const;

    LeafIterator leafbegin () const;
  
    LeafIterator leafend   () const;
    
  public:  
    int size (int level, int codim) const; 

    int size (int level, GeometryType type) const;

    int size (int codim) const; 

    int size (GeometryType type) const;

  public:
    //***************************************************************
    //  Interface for Adaptation
    //***************************************************************
    bool mark( int refCount , const typename Traits::template Codim<0>::EntityPointer & en ) const DUNE_DEPRECATED; 

    int getMark( const typename Traits::template Codim<0>::EntityPointer & ) const DUNE_DEPRECATED; 
    
    int getMark( const typename Traits::template Codim<0>::Entity & ) const; 

    bool mark( int refCount , const typename Traits::template Codim<0>::Entity & en ) const; 

public:
    bool globalRefine(int refCount);
  
    bool adapt ( );

    template <class DofManagerType, class RestrictProlongOperatorType>
    bool adapt (DofManagerType &, RestrictProlongOperatorType &, bool verbose=false );
  
    bool preAdapt ();

    bool postAdapt();

    const CollectiveCommunicationType & comm () const
    {
      return comm_;
    }
    
    std::string name () const { return "AlbertaGrid"; };
  
    //**********************************************************
    // End of Interface Methods
    //**********************************************************
    template <GrapeIOFileFormatType ftype>
    bool writeGrid( const std::basic_string<char> filename, albertCtype time ) const;
  
    template <GrapeIOFileFormatType ftype>
    bool readGrid( const std::basic_string<char> filename, albertCtype & time );

    /* returns size of mesh include all levels 
       max Index of grid entities with given codim 
       for outside the min index is 0, the shift has to done inside 
       the grid which is of minor cost 
    */
    int global_size (int codim) const;

    // return number of my processor 
    int myRank () const { return myRank_; }; 

    void setNewCoords(const FieldVector<albertCtype, dimworld> & trans, const albertCtype scalar);

    // return hierarchic index set 
    const HierarchicIndexSet & hierarchicIndexSet () const { return hIndexSet_; }

    const typename Traits :: LevelIndexSet & levelIndexSet (int level) const;

    const typename Traits :: LeafIndexSet & leafIndexSet () const;
  
    const GlobalIdSet & globalIdSet () const { return globalIdSet_; }
 
    const LocalIdSet & localIdSet () const { return globalIdSet_; }

    // access to mesh pointer, needed by some methods
    ALBERTA MESH* getMesh () const { return mesh_; }; 

    // return real entity implementation 
    template <int cd>
    AlbertaGridEntity<cd,dim,const AlbertaGrid<dim,dimworld> >& 
    getRealEntity(typename Traits::template Codim<cd>::Entity& entity) 
    {
      return this->getRealImplementation(entity);
    }

  private:
    template <int cd>
    const AlbertaGridEntity<cd,dim,const AlbertaGrid<dim,dimworld> >& 
    getRealEntity(const typename Traits::template Codim<cd>::Entity& entity) const 
    {
      return this->getRealImplementation(entity);
    }

  public:
    const std::vector < GeometryType > & geomTypes (int codim) const 
    { 
      assert( codim >= 0 );
      assert( codim < dim+1 );
      return geomTypes_[codim]; 
    }

  private:
    friend class Conversion<AlbertaGrid<dim, dimworld>, HasObjectStream>;  
    friend class Conversion<const AlbertaGrid<dim, dimworld>, HasObjectStream>;

    friend class Conversion<AlbertaGrid<dim, dimworld>, HasHierarchicIndexSet >;  
    friend class Conversion<const AlbertaGrid<dim, dimworld>, HasHierarchicIndexSet>;

    // do not use copy constructor 
    AlbertaGrid(const MyType& other);
    // do not use assigment 
    MyType& operator=(const MyType& other);

  private:
    typedef std::vector<int> ArrayType;
    
    ArrayType ghostFlag_; // store ghost information 
  
    // initialize of some members 
    void initGrid(int proc);
  
    // make the calculation of indexOnLevel and so on.
    // extra method because of Reihenfolge
    void calcExtras(); 
 
    // write ALBERTA mesh file  
    bool writeGridXdr  ( const std::basic_string<char> filename, albertCtype time ) const;

    bool readGridXdr   ( const std::basic_string<char> filename, albertCtype & time );

    bool readGridAscii ( const std::basic_string<char> filename, albertCtype & time );
    
    // delete mesh and all vectors 
    void removeMesh();

    // pointer to an Albert Mesh, which contains the data
    ALBERTA MESH *mesh_;

    // object of collective communication 
    CollectiveCommunicationType comm_;

    // number of maxlevel of the mesh
    int maxlevel_;

    // true if grid was refined or coarsend
    bool wasChanged_; 

    // help vector for setNewCoords 
    mutable ArrayType macroVertices_;

  public:  
    // this method is new fill_elinfo from ALBERTA but here the neighbor
    // relations are calced diffrent, on ervery level there are neighbor
    // realtions ( in ALBERTA only on leaf level ), so we needed a new
    // fill_elinfo. 
    void fillElInfo(int ichild, int actLevel ,const ALBERTA EL_INFO *elinfo_old, 
                    ALBERTA EL_INFO *elinfo, bool hierachical, bool leaf=false ) const;
  
    // calc the neigh[0] 
    void firstNeigh(const int ichild,const ALBERTA EL_INFO *elinfo_old, 
                    ALBERTA EL_INFO *elinfo, const bool leafLevel) const;
  
    // calc the neigh[1] 
    void secondNeigh(const int ichild, const ALBERTA EL_INFO *elinfo_old, 
                     ALBERTA EL_INFO *elinfo, const bool leafLevel) const;
  
    // calc the neigh[2] 
    void thirdNeigh(const int ichild, const ALBERTA EL_INFO *elinfo_old, 
                    ALBERTA EL_INFO *elinfo, const bool leafLevel) const;
 
  private:  
    // needed for VertexIterator, mark on which element a vertex is treated 
    mutable AlbertaMarkerVector vertexMarkerLeaf_; 

    // needed for VertexIterator, mark on which element a vertex is treated 
    mutable AlbertaMarkerVector vertexMarkerLevel_[MAXL]; 
    
    //***********************************************************************
    //  MemoryManagement for Entitys and Geometrys 
    //**********************************************************************
    typedef typename SelectEntityImp<0,dim,const MyType>::EntityObject EntityObject;
  
  public:
    typedef AGMemoryProvider< EntityObject > EntityProvider;

    typedef AlbertaGridIntersectionIterator< const MyType > IntersectionIteratorImp;
    typedef IntersectionIteratorImp LeafIntersectionIteratorImp;
    typedef AGMemoryProvider< LeafIntersectionIteratorImp > LeafIntersectionIteratorProviderType;
    friend class LeafIntersectionIteratorWrapper< const MyType > ;

    typedef LeafIntersectionIteratorWrapper<const MyType >
      AlbertaGridIntersectionIteratorType; 

    LeafIntersectionIteratorProviderType & leafIntersetionIteratorProvider() const { return leafInterItProvider_; }

  private:  
    mutable EntityProvider               entityProvider_;
    mutable LeafIntersectionIteratorProviderType leafInterItProvider_;
   
  public:  
    template< class IntersectionInterfaceType >
    const typename BaseType
      :: template ReturnImplementationType< IntersectionInterfaceType >
      :: ImplementationType & DUNE_DEPRECATED
    getRealIntersectionIterator ( const IntersectionInterfaceType &iterator ) const
    {
      return this->getRealImplementation( iterator );
    }
    
    template< class IntersectionType >
    const typename BaseType
      :: template ReturnImplementationType< IntersectionType >
      :: ImplementationType &
    getRealIntersection ( const IntersectionType &intersection ) const
    {
      return this->getRealImplementation( intersection );
    }
   
    // (for internal use only) return obj pointer to EntityImp 
    template <int codim> 
    typename SelectEntityImp<codim,dim,const MyType>::EntityObject * 
    getNewEntity (int level, bool leafIt ) const;

    // (for internal use only) free obj pointer of EntityImp 
    template <int codim>
    void freeEntity (typename SelectEntityImp<codim,dim,const MyType>::EntityObject * en) const;

  private:
    //*********************************************************************
    // organisation of the global index 
    //*********************************************************************
    // provides the indices for the elements 
    IndexManagerType indexStack_[AlbertHelp::numOfElNumVec];
  
    // the DOF_INT_VECs we need 
    // * change to mutable here
    mutable ALBERTA AlbertHelp::DOFVEC_STACK dofvecs_;
        
    const ALBERTA DOF_ADMIN * elAdmin_;

    const ALBERTA REAL_D * coordsVec_;
    
    // pointer to vec of elNumbers_
    const int * elNewVec_;

    // for access in the elNewVec and ownerVec 
    const int nv_;
    const int dof_; 
 
  public:  
    // make some shortcuts 
    void arrangeDofVec(); 

    // return true if el is new 
    bool checkElNew ( const ALBERTA EL * el ) const;
  
    // read global element number from elNumbers_  
    const ALBERTA REAL_D & getCoord ( const ALBERTA EL_INFO * elInfo, int vx ) const
    {
      assert( vx>= 0);
      assert( vx < dim+1 );
#ifndef CALC_COORD 
      assert(coordsVec_);
      return coordsVec_[elInfo->el->dof[vx][0]];
#else 
      return elInfo->coord[vx];
#endif
    } 

    // read level from elNewCehck vector   
    int getLevelOfElement ( const ALBERTA EL * el ) const; 

    // read global element number from elNumbers_  
    int getElementNumber ( const ALBERTA EL * el ) const; 

    // read global element number from elNumbers_  
    int getEdgeNumber ( const ALBERTA EL * el, int edge ) const; 

    // read global element number from elNumbers_  
    int getFaceNumber ( const ALBERTA EL * el, int face ) const; 

    // read global element number from elNumbers_  
    int getVertexNumber ( const ALBERTA EL * el, int vx ) const; 

  private:
    // rank of my thread, i.e. number of my processor 
    const