// $Id: AlbertGrid.hh,v 1.12 2003/02/13 16:32:58 robertk Exp $

#ifndef DUNE_ALBERTGRID_HH
#define DUNE_ALBERTGRID_HH

//#include <array.hh>
#include <grid.hh>


namespace Albert
{
#define DIM MYDIM
#define DIM_OF_WORLD MYDIMOW
}

#include <albert.h>


using namespace std;

#define FILL_NOTHING               0X00L
#define FILL_COORDS                0X01L
#define FILL_BOUND                 0X02L
#define FILL_NEIGH                 0X04L
#define FILL_OPP_COORDS            0X08L
#define FILL_ORIENTATION           0X10L
#define CALL_LEAF_EL               0X0800L


namespace Albert
{

  typedef struct albert_leaf_data
  {
    /// Wegen Fehler in memory.cc, der noch behoben werden muss 
    S_CHAR reached[4]; // N_NEIGH ? Element durchlaufen ? 
    //S_CHAR reached[N_VERTICES]; // N_NEIGH ? Element durchlaufen ? 
  } AlbertLeafData;

  void setReached(const EL_INFO *elInfo)
  {
    AlbertLeafData * ldata = (AlbertLeafData *) elInfo->el->child[1];
    AlbertLeafData* ldataNeigh = NULL;
    for(int i=0; i<N_VERTICES; i++)
    {

      if(elInfo->neigh[i])
      {
        if(ldata->reached[i] != -1)
        {
          ldata->reached[i] = 1;
          ldataNeigh = (AlbertLeafData *) elInfo->neigh[i]->child[1];
          ldataNeigh->reached[elInfo->opp_vertex[i]] = -1;
        }
        else 
        {
          ldataNeigh = (AlbertLeafData *) elInfo->neigh[i]->child[1];
          ldataNeigh->reached[elInfo->opp_vertex[i]] = 1;
        }
      }
      else 
       ldata->reached[i] = 1;
    }  

  }
  void AlbertLeafRefine(EL *parent, EL *child[2])
  {
    AlbertLeafData *ldata;
    AlbertLeafData *ldataChi[2];

    ldata = (AlbertLeafData *) parent->child[1];
//    assert(ldata != NULL);
//    assert(ldata->volume != 0.0);
//    cout << "AlbertLeafRefine: Muss noch getestet werden \n";

    ldataChi[0] = (AlbertLeafData *) child[0]->child[1];
    // die neue Kante auf Child 0, siehe Albert Doc S. 11
#if DIM == 3
    ldataChi[0]->reached[0] = 1;
#else
    ldataChi[0]->reached[1] = 1;
#endif

    ldataChi[1] = (AlbertLeafData *) child[1]->child[1];
    // die neue Kante auf Child 1, siehe Albert Doc S. 11
    ldataChi[1]->reached[0] = -1;
  }

  void AlbertLeafCoarsen(EL *parent, EL *child[2])
  {
    // do notin' 
  }

  void initLeafData(LEAF_DATA_INFO *linfo)
  {
    cout << "AlbertLeafRefine: Muss noch getestet werden \n";
    linfo->leaf_data_size = sizeof(AlbertLeafData);
    linfo->refine_leaf_data = &AlbertLeafRefine;
    linfo->coarsen_leaf_data = &AlbertLeafCoarsen;

    return;
  }

const BOUNDARY *initBoundary(MESH *Spmesh, int bound)
{
  FUNCNAME("initBoundary");
  static const BOUNDARY Diet1 = {NULL,1};
  static const BOUNDARY PaulN1  = {NULL,-1};
  
  static const BOUNDARY Diet2 = {NULL,2};
  static const BOUNDARY PaulN2  = {NULL,-2};
 
  static const BOUNDARY Diet3 = {NULL,3};
  static const BOUNDARY PaulN3  = {NULL,-3};

  static const BOUNDARY Diet4 = {NULL,4};
  static const BOUNDARY PaulN4  = {NULL,-4};
 
  static const BOUNDARY Diet5 = {NULL,5};
  static const BOUNDARY PaulN5  = {NULL,-5};


  switch(bound)
  {
    case 1: return (&Diet1);
    case 2: return (&Diet2);
    case 3: return (&Diet3);
    case 4: return (&Diet4);
    case 5: return (&Diet5);

    case -1: return (&PaulN1);
    case -2: return (&PaulN2);
    case -3: return (&PaulN3);
    case -4: return (&PaulN4);
    case -5: return (&PaulN5);
    default: ERROR_EXIT("no Boundary for %d. Och! \n",bound);
  }
  return &Diet1;
}


}

namespace Dune 
{

//  template <class BASE> class AlbertGridLevel;

////////////////////////////////////////////////////////////////////////
//
// --- AlbertGrid
//
////////////////////////////////////////////////////////////////////////
  template <class BASE, int CoDIM> class AlbertGridLevelAccess; 

  template <int DIMOFWORLD, int ADIM>
  class AlbertGrid : public Grid <AlbertGrid<DIMOFWORLD,ADIM> > 
  {
  public: 
    typedef AlbertGridLevelAccess<AlbertGrid<DIMOFWORLD, ADIM> , 0 > ELEMENTACCESS;
    typedef AlbertGridLevelAccess<AlbertGrid<DIMOFWORLD, ADIM> , 1 > FACEACCESS;
#if ADIM > 2    
    typedef AlbertGridLevelAccess<AlbertGrid<DIMOFWORLD, ADIM> , 2 >  EDGEACCESS;
#endif
    typedef AlbertGridLevelAccess<AlbertGrid<DIMOFWORLD, ADIM> , ADIM >  VERTEXACCESS;
    
    class GetAccess;
    
  public: 
    static const int BASEDIM = ADIM;
    static const int BASEDOW = DIMOFWORLD;
    Albert::MESH *mesh_;
  public:  
    class GetAccess : public Grid<AlbertGrid<DIMOFWORLD,ADIM> >::GetAccess 
    {
    private:
      AlbertGrid<DIMOFWORLD,ADIM> *grid_;
      
    public:
      GetAccess(AlbertGrid<DIMOFWORLD,ADIM> *grid) 
      {
        grid_ = grid;
      };

      GetAccess(const GetAccess &ga)
      {
        grid_ = ga.grid_;
      };
      
      ELEMENTACCESS getElementAccess()
      {
        ELEMENTACCESS access(*grid_);
        return access;
      }

      FACEACCESS getFaceAccess()
      {
        FACEACCESS access(*grid_);
        return access;
      }

      VERTEXACCESS getVertexAccess()
      {
        VERTEXACCESS access(*grid_);
        return access;
      }

#if ADIM > 2    
      EDGEACCESS getEdgeAccess()
      {
        EDGEACCESS access(*grid_);
        return access;
      }
#endif 

    }; // end GetAccess

    GetAccess *ga_;
  public: 
    AlbertGrid(char* MacroTriangFilename)
    {
      mesh_ = Albert::get_mesh("AlbertGrid",NULL,Albert::initLeafData);//initDofAdmin,initLeafData);
      Albert::read_macro(mesh_,MacroTriangFilename,Albert::initBoundary);

      Albert::mesh_traverse(mesh_,-1,CALL_LEAF_EL|FILL_NEIGH,
                            Albert::setReached);
      ga_ = new GetAccess (this);

    };

  
  public:
    /// Verfeinert das Gitter um refCount Levels
    void globalRefine(int refCount)
    {
      unsigned char flag;

      // trage auf jedem Element refCount ein 
      flag = Albert::global_refine(mesh_,refCount);
   
      // verfeinere 
      refine();
    };

    /// Verfeinert das Gitter um die auf jedem Element 
    /// festgelegte Zahl von Verfeinerungen 
    void refine()
    {
      unsigned char flag;
      flag = Albert::refine(mesh_);
    };
    
    void coarsen()
    {
      unsigned char flag;
      flag = Albert::coarsen(mesh_);
    };

    void writeGrid()
    {
      cout << "writeGrid () not implemented in this DIM \n";
      abort();
    };
  
    GetAccess& getAccess()
    {
      return (*ga_);
    }
    
  }; // end Class AlbertGridGrid


//////////////////////////////////////////////////////////////////////////
//
//  AlbertGridLevelAccess 
//
/////////////////////////////////////////////////////////////////////////
    // BASE stellt hier das Gitter dar
    template <class BASE, int CoDIM>
    class AlbertGridLevelAccess : public 
    LevelAccess<AlbertGridLevelAccess<BASE,CoDIM> > 
    {  
    public:
      
      class RefElement;
      class Element;
      BASE *grid_;

      static const int BASEDIM = BASE::BASEDIM;
      static const int BASEDOW = BASE::BASEDOW;
      static const int numberOfVertices_ = BASEDIM+1;

      /// Nachbarn sind immer Elemente 
    public:  
      typedef int level;

      typedef typename LevelAccess<AlbertGridLevelAccess <BASE ,CoDIM > >::Element INTERelement;
    
      typedef typename LevelAccess<AlbertGridLevelAccess <BASE ,CoDIM > >::
                       Iterator INTERiterator;
      typedef typename AlbertGridLevelAccess<BASE , 0>::Element NEIGHS;
      typedef typename LevelAccess<AlbertGridLevelAccess <BASE ,0 > >
                       ::Element NEIGHtype;
    public:                   
      typedef double REAL;
      typedef REAL REAL_D[BASEDOW];
                       
      ///////////////////////////////////////////////////
    public:                   
      typedef NEIGHtype** NEIGHVEC;
      typedef REAL** VERTEXtype;

/// CoDIM 0 ==> Access Elements
/// CoDIM 1 ==> Access Faces
/// CoDIM 2 ==> Access Edges 
/// CoDIM 3 ==> Access Vertices
//
//
/////////////////////////////////////////////////////////////////////////////
//
//  Class Element 
//
////////////////////////////////////////////////////////////////////////////
    class Element : public LevelAccess<AlbertGridLevelAccess<BASE , CoDIM > >::Element 
    {
    private:
      /// always DIM+1-CoDIM Vertices on an Element
      /// Triangle ==> 3 Vertices (DIM =3, CoDIM =1)
      static const int numberOfVertices = BASEDIM+1-CoDIM; 
      
      /// number of neighbours diffrent for CoDIM s
      /// an Element has DIM+1 neighbours
      static const int numberOfNeighbours = BASEDIM+1-CoDIM*(BASEDOW-1);
      
      Albert::EL_INFO *elInfo_;// hier steht das Element des Gitters
      int elementNumber_;
      int face_;
      REAL** vertex_;
      NEIGHS** neighbours_;
      NEIGHS** realNeighbours_;

    public:
      ~Element() { cout << "Destroying Element \n"; };
      
      /// The Element is prescribed by the EL_INFO struct of Albert::MESH
      /// the pointer to this struct is set and get by setElInfo and
      /// getElInfo. 
      void makeElement()
      {
        //cout << "NumberOfVertices " << numberOfVertices << endl;
        //cout << "NumberOfNeighbours " << numberOfNeighbours << endl;
        //cout << "CoDIM " << CoDIM << endl;
        elementNumber_ = -1;
        face_ = 0;
        elInfo_ = NULL;
        neighbours_ = NULL;
        realNeighbours_ = NULL;
        
        // Jedes Element besitzt Vertices, also Speicher holen
        vertex_ = new REAL* [numberOfVertices];
        for(int i=0; i<numberOfVertices; i++)
          vertex_[i] = NULL;
      };

      void makeNeighbours()
      {
        neighbours_ = new NEIGHS* [numberOfNeighbours];
        //for(int i=0; i<numberOfNeighbours; i++)
        //  neighbours_[i] = new NEIGHS ();
          
        realNeighbours_ = new NEIGHS* [numberOfNeighbours];
        for(int i=0; i<numberOfNeighbours; i++)
          realNeighbours_[i] = new NEIGHS ();
      };

      Element (bool makeNeigh=false) 
      {  
        makeElement();
        if(makeNeigh)
          makeNeighbours();
      };  

#if 0
      Element (const Element& el)
      {
        cout << "Copy Element \n";
        makeElement();
        if(el.neighbours_)
          makeNeighbours();

        if(el.elInfo_)
          setElInfo(el.elInfo_,el.face_,el.elementNumber_);
      };
#endif

      Albert::EL_INFO* getElInfo()
      {
        return elInfo_;
      }

      void setElInfo(Albert::EL_INFO *elInfo,int f,int e)
      {
        face_ = f;
        elementNumber_ = -1;
        elInfo_ = elInfo;
        
        // Achtung spezielle Implementierung je nach DIM
        if(elInfo_)
        {
          elementNumber_ = elInfo_->el->index; // e;
          coordToVertex();
          setVerticesAndNeighs();
        }
      }

      void coordToVertex()
      { 
        // funktioniert fuer CoDIM 0 und CoDIM 1
        // fuer CoDIM 0 muss face_ = 0 sein 
        // Number of Vertices of EL_INFO is always BASEDIM+1 
        for(int i=0; i<numberOfVertices; i++)
          vertex_[i] =(REAL *) &elInfo_->coord[(face_+i)%numberOfVertices];
      }

      void setElementNeighbours()
      {
        for(int i=0; i<numberOfNeighbours; i++)
        {
          neighbours_[i] = NULL;
          realNeighbours_[i]->setupAsNeighbour(elInfo_,i); 
          if(elInfo_->neigh[i])
            neighbours_[i] = realNeighbours_[i];
        }
      }

      void setFaceNeighbours()
      {
        // numberOfNeighbours == 2
        realNeighbours_[0]->setupAsNeighbour(elInfo_,-1);
        neighbours_[0] = realNeighbours_[0];
        realNeighbours_[1]->setupAsNeighbour(elInfo_,face_);
        neighbours_[1] = realNeighbours_[1];
      }

      void setVerticesAndNeighs()
      {
        cout << "setVerticesAndNeighs() not implemented in this CoDIM \n";
      };
      
      void setupAsNeighbour(Albert::EL_INFO *elInfo,int neigh)
      {
        cout << "setupAsNeighbour(..)  not implemented in this CoDIM \n";
      };
     
      // nur fuer CoDIM 0 
      void setTheRealNeighbour(Albert::EL_INFO *elInfo,int neigh)
      {
        //cout << "setup Neighbour "<<neigh<<  "\n";
        if(elInfo)
        {
          face_ = 0;
          elInfo_ = elInfo;  
          if(neigh == -1)
          {
            elementNumber_ = elInfo->el->index;
            coordToVertex();  
          }
          else
          {
            vertex_[0] = (REAL *) &elInfo_->opp_coord[neigh];
            for(int i=1; i<numberOfVertices; i++)
              vertex_[i] = (REAL *) &elInfo_->coord[(neigh+i)%numberOfVertices];
            
            if(elInfo_->neigh[neigh])  
              elementNumber_ = elInfo->neigh[neigh]->index;
          }
        }
        //cout << "setup Neighbour done \n";
      };

      int getElNum()
      {
        return elementNumber_;
      };

      void printElInfo()
      {
        if(elInfo_)
        for(int i=0; i<numberOfVertices; i++)
        {
          printf("ElInfo %d: ",i);
          for(int j=0; j<BASEDOW; j++)
            printf(" %lf ",elInfo_->coord[i][j]);
          printf("\n");  
        }
      }

      void printVertex()
      {
        for(int i=0; i<numberOfVertices; i++)
        {
          if(vertex_[i])
          {
            printf("Vertex %d: ",i);
            for(int j=0; j<BASEDOW; j++)
              printf(" %lf ",vertex_[i][j]);
            printf("\n");  
          }
          else 
           printf("No Vertex %d \n",i);
        }
      }

      int getDof(int localNumber)
      {
        // 0 heisst hier Wert im Knoten, also lineare Elemente 
        return elInfo_->el->dof[localNumber][0];
      };

      VERTEXtype vertex()
      {
        return vertex_;
      };
      
      NEIGHtype*  getNeighbours (int neigh)
      {
        if(neighbours_[neigh])
          return (NEIGHtype *) neighbours_[neigh];
        else 
        {  
//          cout << "Kein Nachbar vorhanden \n";
          return NULL;
        }
      };

    };

    public:
      AlbertGridLevelAccess()
      {
        grid_ = NULL;
      };

      AlbertGridLevelAccess(BASE &grid)
      {
        grid_ = &grid;
      };


      //////////////////////////////////////////////////////////
      //
      //  AlbertGridLevelAccess::Iterator
      //
      //////////////////////////////////////////////////////////
      class Iterator : public
      LevelAccess<AlbertGridLevelAccess<BASE,CoDIM> >::Iterator
      {
      private:
        Element *el_;
        Albert::TRAVERSE_STACK *travStack_;
        int vertexCount_;
//        int elementNumber_;

      public:
        void makeIterator()
        {
          cout << "AlbertGridLevelAccess::Iterator \n";
          //elementNumber_ =0;
          vertexCount_ = 0;
          travStack_ = NULL;
          el_ = new Element ();
        };

        Iterator() 
        {  
          makeIterator();
        };

        Iterator(Albert::MESH *mesh,level travLevel) 
        {  
          vertexCount_ = 0;
          //elementNumber_ =0;
          if(mesh)
          {
            Albert::FLAGS travFlags = FILL_COORDS|FILL_NEIGH;
            if(travLevel >= 0)
              travFlags= travFlags|CALL_EL_LEVEL;
            else  
              travFlags= travFlags|CALL_LEAF_EL;
          
            // get traverse_stack
            travStack_ = Albert::get_traverse_stack();
            el_ = new Element (true);
            //cout << "Vertex " << el->vertex_ <<endl;
            //cout << "Iterator neues Element \n";
            el_->setElInfo( 
              Albert::traverse_first(travStack_,mesh,travLevel,travFlags),
              0,0);
          }
          else
            makeIterator();
          
        };

        Iterator (const Iterator& I)
        {
          //elementNumber_ = I.elementNumber_ ;
          vertexCount_ = I.vertexCount_;
          travStack_ = I.travStack_;
          el_ = I.el_; // Zeiger kopieren 
        };

        bool operator ==(const Iterator& I)
        {
          return (el_->getElInfo() == I.el_->getElInfo());
        };

        // gehe zum naechsten Element, wie auch immer
        Iterator& operator ++() 
        {
          //elementNumber_++;
          // die 0 ist wichtig, weil Face 0, heist hier jetzt Element
          el_->setElInfo(Albert::traverse_next(travStack_,el_->getElInfo()),0,0);//elementNumber_);
          return (*this);
        };
              
        Element& operator* ()
        {
          return (*el_);
        }

    private:
      Albert::EL_INFO* goNextFace(Albert::EL_INFO *elInfo)
      {
        // go next Element, Face 0
        vertexCount_++;
        if(vertexCount_ >= numberOfVertices_)  
        {
          elInfo = Albert::traverse_next(travStack_,elInfo);
          vertexCount_ = 0;
        }

        return elInfo;
      }

      Albert::EL_INFO* goNext(Albert::EL_INFO *elInfo)
      {
        Albert::AlbertLeafData * ldata = NULL;

        // go next Face // elInfo, Anz Vertices, DOW 
        elInfo = goNextFace(elInfo);
        if(!elInfo) return elInfo; // if no more Faces, return

        if(elInfo->el->child[1])
        {
          ldata = (Albert::AlbertLeafData *) elInfo->el->child[1];

          // if reached before, go next 
          if(ldata->reached[vertexCount_] != 1)
            elInfo = goNext(elInfo);
        }

        return elInfo;
      }

    };



    ////////////////////////////////////////////////////////
    //
    //  Hierarchischer Gitterdurchlauf
    //
    ////////////////////////////////////////////////////////
    class HierarchicIterator : public Iterator 
    {
    private:
      Element *el_;
      Albert::TRAVERSE_STACK *travStack_;
      int vertexCount_;
      int elementNumber_;
      
    public:
      HierarchicIterator ()   
      {
        makeIterator();
      };
         
      
      HierarchicIterator(Albert::MESH *mesh,level travLevel,
                 Element *macroElement) 
      {  
        vertexCount_ = 0;
        elementNumber_ =0;
        if(mesh)
        {
          Albert::FLAGS travFlags = FILL_COORDS|FILL_NEIGH;
          if(travLevel >= 0)
            travFlags= travFlags|CALL_EL_LEVEL;
          else  
            travFlags= travFlags|CALL_LEAF_EL;
          
          // get traverse_stack
          travStack_ = Albert::get_traverse_stack();
          travStack_->traverse_mesh = mesh;
          travStack_->traverse_level = travLevel;
          travStack_->traverse_fill_flag = travFlags;
          el_ = new Element (true);
          Albert::EL_INFO *elInfo = new Albert::EL_INFO (); 

          /// Hier kann ein beliebiges Element uebergeben werden,
          /// da jedes AlbertElement einen Zeiger auf das Macroelement
          /// enthaelt.
          const Albert::MACRO_EL *macroEl; 
          //=mesh->first_macro_el;  
          if(!macroEl)
          {
            cout << "Kein Macro Element \n";
            abort();
          }
          if(macroElement)
            macroEl = macroElement->getElInfo()->macro_el;
          else 
            macroEl = mesh->first_macro_el;

          Albert::fill_macro_info(macroEl,elInfo);
          el_->setElInfo(elInfo,0,0);
        }
        else
          makeIterator();
      };

      /// wird in Iterator implementiert 
      /// HierarchicIterator (const HierarchicIterator& I)
      
      /// wird in Iterator implementiert
      /// bool operator ==(const HierarchicIterator& I)
      
      HierarchicIterator& operator ++() 
      {
        elementNumber_++;
        // die 0 ist wichtig, weil Face 0, heist hier jetzt Element
        el_->setElInfo(recursiveTraverse(travStack_),0,elementNumber_);
                
        return (*this);
      };

    private:
      Albert::EL_INFO *recursiveTraverse(Albert::TRAVERSE_STACK *stack)
      {
        Albert::EL *el;
        int i;

        if (stack->stack_used == 0)   /* first call */
        {
          stack->traverse_mel = stack->traverse_mesh->first_macro_el;
          if (stack->traverse_mel == nil)  return(nil);

          stack->stack_used = 1;
          fill_macro_info(stack->traverse_mel,
              stack->elinfo_stack+stack->stack_used);
          stack->info_stack[stack->stack_used] = 0;

          el = stack->elinfo_stack[stack->stack_used].el;
          if ((el == nil) || (el->child[0] == nil)) 
          {
            return(stack->elinfo_stack+stack->stack_used);
          }
        }
        else
        {
          el = stack->elinfo_stack[stack->stack_used].el;

          /* go up in tree until we can go down again */
          while((stack->stack_used > 0) &&
            ((stack->info_stack[stack->stack_used] >= 2) || (el->child[0]==nil)))
          {
            stack->stack_used--;
            el = stack->elinfo_stack[stack->stack_used].el;
          }

          /* goto next macro element */
          if (stack->stack_used < 1) 
          {
            stack->traverse_mel = stack->traverse_mel->next;
            if (stack->traverse_mel == nil)  return NULL;

            stack->stack_used = 1;
            Albert::fill_macro_info(stack->traverse_mel,
            stack->elinfo_stack+stack->stack_used);
            stack->info_stack[stack->stack_used] = 0;

            el = stack->elinfo_stack[stack->stack_used].el;
            if ((el == NULL) || (el->child[0] == NULL)) 
              return(stack->elinfo_stack+stack->stack_used);
          }
        }

        /* go down next child */
        if(el->child[0])
        {
          if (stack->stack_used >= stack->stack_size-1)
            Albert::enlarge_traverse_stack(stack);

          i = stack->info_stack[stack->stack_used];
          el = el->child[i];
          stack->info_stack[stack->stack_used]++;
          Albert::fill_elinfo(i, stack->elinfo_stack+stack->stack_used,
              stack->elinfo_stack+stack->stack_used+1);
          stack->stack_used++;

          stack->info_stack[stack->stack_used] = 0;
        }

        return(stack->elinfo_stack+stack->stack_used);

      }; // recursive traverse
      
   }; // end Iterator 


    public:
      Iterator begin(int l=-1)
      {
        // first Macro Element
        cout << "getin begin Iterator \n";
        Iterator *it = 
            new Iterator (grid_->mesh_,l);
        return (*it);
      };

      Iterator end(int l=-1)
      {
        // behind last Macro Element, is NULL
        cout << "getin end Iterator \n";
        Iterator *it = new Iterator ();
        return (*it);
      };

      // hier wird ein MacroElement uebergeben
      HierarchicIterator firstChild(int l,INTERelement *el)
      {
        cout << "getin begin HierarchicIterator \n";
        HierarchicIterator *it = 
            new HierarchicIterator (grid_->mesh_,l,(Element *) el);
        return (*it);
      };
      
      HierarchicIterator lastChild(int l,INTERelement *el)
      {
        // behind last Macro Element, is NULL
        cout << "getin end HierarchicIterator \n";
        HierarchicIterator *it = new HierarchicIterator();
        return (*it);
      };

      
     
    }; //end LevelAccess 


//////////////////////////////////////////////////////////////
//
// Spezielle Spezifikationen fuer einzelne Dimensionen
//
//////////////////////////////////////////////////////////////
    void AlbertGridLevelAccess<AlbertGrid<3,3>, 0>::Element::
    setVerticesAndNeighs()
    {
      setElementNeighbours();
    }
    void AlbertGridLevelAccess<AlbertGrid<2,2>, 0>::Element::
    setVerticesAndNeighs()
    {
      setElementNeighbours();
    }

    void AlbertGridLevelAccess<AlbertGrid<3,3>, 1>::Element::
    setVerticesAndNeighs()
    {
      setFaceNeighbours();
    }
    void AlbertGridLevelAccess<AlbertGrid<2,2>, 1>::Element::
    setVerticesAndNeighs()
    {
      setFaceNeighbours();
    }


    void AlbertGridLevelAccess<AlbertGrid<3,3>, 0>::Element::
    setupAsNeighbour(Albert::EL_INFO *elInfo,int neigh)
    {
      setTheRealNeighbour(elInfo,neigh);
    };
      
    void AlbertGridLevelAccess<AlbertGrid<2,2>, 0>::Element::
    setupAsNeighbour(Albert::EL_INFO *elInfo,int neigh)
    {
      setTheRealNeighbour(elInfo,neigh);
    };
      
    AlbertGridLevelAccess<AlbertGrid<3,3>, 1>::Iterator&
    AlbertGridLevelAccess<AlbertGrid<3,3>, 1>::Iterator::operator ++ ()
    {
      //elementNumber_++;
      el_->setElInfo(goNext(el_->getElInfo()),
                    vertexCount_,0);//elementNumber_);
      return (*this);
    }

    AlbertGridLevelAccess<AlbertGrid<2,2>, 1>::Iterator&
    AlbertGridLevelAccess<AlbertGrid<2,2>, 1>::Iterator::operator ++ ()
    {
      //elementNumber_++;
      el_->setElInfo(goNext(el_->getElInfo()),
                    vertexCount_,0);//elementNumber_);
      return (*this);
    }

    //template <int DIMOFWORLD, int ADIM>
    void AlbertGrid<3,3>::writeGrid()
    {
      cout << "\nWriting 3d Grid \n";

      AlbertGridLevelAccess<AlbertGrid<3,3> ,0 >::Iterator endit;
      AlbertGridLevelAccess<AlbertGrid<3,3> ,0 >::Iterator it;
      AlbertGridLevelAccess<AlbertGrid<3,3>,0 >
        access ( (*this));
             
      double **coord = new double* [mesh_->n_vertices]; 
      for(int i=0; i<mesh_->n_vertices; i++)
        coord[i] = new double [BASEDOW];

      //LEVit endit = access.end(-1);
      endit = access.end(-1);
      /// die eigentlichen ElementNummer auf die Zahlen 0 bis n_elements-1
      /// abbilden 
      int *lookup    = new int [mesh_->n_hier_elements];
      int count = 0;
      for(it = access.begin(-1); it != endit; ++it)
      {
        lookup[(*it).getElNum()] = count;
        count++;
      }

      int **vertex   = new int* [mesh_->n_elements];
      for(int i=0; i<mesh_->n_elements; i++)
        vertex[i] = new int [BASEDIM+1];

      /// setup the Wesenber 3d Grid 
      for(it = access.begin(-1); it != endit; ++it)
      {
        int elNum = lookup[(*it).getElNum()];
        
        for(int i=0; i<BASEDIM+1; i++)
        {
          int k = (*it).getDof(i); 

          vertex[elNum][i] = k;

          double ** elCoords = (*it).vertex();
          for(int j = 0; j<BASEDOW; j++)
            coord[k][j] = elCoords[i][j];
        }           
      }

      /// write the USPM Mesh
      FILE *file = fopen("test3dOUT.0","w");
      if(!file)
      {
        cout << "Couldnt open test3dOUT.0 \n";
        abort();
      }

      // die Zeit 
      fprintf(file,"0.0 \n");
      fprintf(file,"%d \n",mesh_->n_vertices);

      for(int i=0; i<mesh_->n_vertices; i++)
      {
        for(int j =0 ; j<BASEDOW; j++)
          fprintf(file,"%le ",coord[i][j]);
        fprintf(file,"\n");  
      }

      fprintf(file,"%d \n",mesh_->n_elements);
     
      for(int i=0; i<mesh_->n_elements; i++)
      {
        for(int j =0 ; j<BASEDIM+1; j++)
          fprintf(file,"%d ",vertex[i][j]);
        double a = (double) i;  
        fprintf(file,"%f \n",a); 
      }

      fclose(file);
      system("gzip test3dOUT.0");
      
    };
      
    void AlbertGrid<2,2>::writeGrid()
    {
      cout << "\nWriting USPM grid \n\n";
      
      /// USPM 2d 
      AlbertGridLevelAccess<AlbertGrid<2,2>,0 >::Iterator endit;
      AlbertGridLevelAccess<AlbertGrid<2,2>,0 >::Iterator it;
      AlbertGridLevelAccess<AlbertGrid<2,2>,0 > 
         access ( (*this));

      double **coord = new double* [mesh_->n_vertices]; 
      for(int i=0; i<mesh_->n_vertices; i++)
        coord[i] = new double [BASEDOW];
      cout << "Number of Vertices " << mesh_->n_vertices << endl;

      endit = access.end(-1);


      /// die eigentlichen ElementNummer auf die Zahlen 0 bis n_elements-1
      /// abbilden 
      int *lookup    = new int [mesh_->n_hier_elements];
      int count = 0;
      for(it = access.begin(-1); it != endit; ++it)
      {
        lookup[(*it).getElNum()] = count;
        count++;
      }

      int **nb       = new int* [mesh_->n_elements];
      for(int i=0; i<mesh_->n_elements; i++)
        nb[i] = new int [BASEDIM+1];
       

      
      int **vertex   = new int* [mesh_->n_elements];
      for(int i=0; i<mesh_->n_elements; i++)
        vertex[i] = new int [BASEDIM+1];

      /// setup the USPM Mesh 
      for(it = access.begin(-1); it != endit; ++it)
      {
        int elNum = lookup[(*it).getElNum()];
        
        for(int i=0; i<BASEDIM+1; i++)
        {
          int k = (*it).getDof(i); 
          //printf("Vertex Nr %d Dof Nr. %d\n",i,k);
          //printf("El Num %d \n",elNum);

          vertex[elNum][i] = k;

          if((*it).getNeighbours(i))
            nb[elNum][i] = lookup[((*it).getNeighbours(i))->getElNum()];
          else nb[elNum][i] = -1;

          double ** elCoords = (*it).vertex();
          for(int j = 0; j<BASEDOW; j++)
            coord[k][j] = elCoords[i][j];
        }           
        //cout << "done \n";
      }

      /// write the USPM Mesh
      FILE *file = fopen("out.uspm","w");
      if(!file)
      {
        cout << "Couldnt open out.uspm \n";
        abort();
      }
      fprintf(file,"USPM 2\n");
      fprintf(file,"%d %d \n",N_VERTICES,N_VERTICES*mesh_->n_elements);
      fprintf(file,"%d %d 0\n",mesh_->n_elements,mesh_->n_vertices);

      for(int i=0; i<mesh_->n_vertices; i++)
      {
        fprintf(file,"%d ",i);  
        for(int j =0 ; j<BASEDOW; j++)
          fprintf(file,"%le ",coord[i][j]);
        fprintf(file,"\n");  
      }
     
      for(int i=0; i<mesh_->n_elements; i++)
      {
        fprintf(file,"%d ",i);  
        for(int j =0 ; j<BASEDIM+1; j++)
          fprintf(file,"%d ",vertex[i][j]);
        for(int j =0 ; j<BASEDIM+1; j++)
          fprintf(file,"%d ",nb[i][j]);
        fprintf(file,"\n"); 
      }

      fclose(file);
    };
    



}; // namespace Dune

#endif
