src/expde/basic/basic.h

Go to the documentation of this file.

//    expde: expression templates for partial differential equations.
//    Copyright (C) 2001  Christoph Pflaum
//    This program is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation; either version 2 of the License, or
//    (at your option) any later version.
//
//    This program is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//                 SEE  Notice1.doc made by 
//                 LAWRENCE LIVERMORE NATIONAL LABORATORY
//

// ------------------------------------------------------------
//
// basic.h
//
// ------------------------------------------------------------


#ifndef BASIC_H_
#define BASIC_H_
             
class Index3D;
class Grid_base;
class Grid;
class IndexEdge;
class IndexSurface;
class Parallel_Info;


// 1D index
class Index1D {
 friend class Index3D;
 friend class Grid_base;
 friend class Parallel_Info;
 private:
  long_int index;
 public:
  Index1D(long_int ind) {
    /*
    if(developer_version) {
      if(ind==0) {
        cout << " error in Index1D::Index1D(long_int ind) " << endl;
      }
    }
    */
    index = ind; 
  };

#ifdef PERIODIC
  Index1D()                   { index = 1; };
#else
  Index1D()                   { index = 0; };
#endif
  inline long_int get() const { return index; };

  inline double coordinate() const;
  inline int    Tiefe() const;
  inline double h_loc() const;
  inline Index1D son(dir1D d) const;
  inline Index1D father();
  inline Index1D neighbour(dir1D d) const;
  inline Index1D neighbour_non_periodic(dir1D d) const;
  inline bool    Is_non_periodic() const;

  inline Index1D next(dir1D d, int l) const;

  inline long_int cart_index(int l) const;       // cartesian index of level l
                                            // is 0,1,...,2^l

  inline void Print() const { cout << index; }

  inline Index1D next(Ort1D ort, int l) const;

  inline int  operator!=(Index1D I) const;
  inline int  operator==(Index1D I) const;
  inline bool operator<=(Index1D I) const;  // point finer or equal as this
  inline bool operator< (Index1D I) const;  // I_this subset [I-h,I+h[
  inline bool operator<<(Index1D I) const;  // I_this subset [I-h,I+h]

  inline Index1D cut_level(int l) const;
};

// 3D index
class Index3D {
 friend class Grid_base;
 friend class Grid;
 friend class Parallel_Info;
 private:
  Index1D ind_x;
  Index1D ind_y;
  Index1D ind_z;
 public:
  Index3D()  { ind_x = Index1D(); ind_y = Index1D(); ind_z = Index1D(); };
  // Index3D(2,2,2) : middle  point: (0.5,0.5,0.5)
  // Index3D(3,3,3) : outside point: (1.5,1.5,1.5)
  Index3D(long_int intx, long_int inty, long_int intz);
  Index3D(Index1D indx,Index1D indy,Index1D indz);
  int     Tiefe() const;
  Index3D son(dir_sons) const;
  Index3D father();
  Index3D son_EW(dir1D);
  Index3D son_NS(dir1D);
  Index3D son_TD(dir1D);
  Index3D next(Ort1D ortx, Ort1D orty, Ort1D ortz, int l) const;
  Index3D next_NS(dir1D d, int l) const;
  Index3D next_EW(dir1D d, int l) const;
  Index3D next_TD(dir1D d, int l) const;

  Index3D next(dir_sons d, int l) const;
  Index3D next(dir_3D d, int l) const;
  Index3D next(Edges_cell ed, int l) const;

  Index3D next_N(int l) const;
  Index3D next_S(int l) const;
  Index3D next_W(int l) const;
  Index3D next_E(int l) const;
  Index3D next_T(int l) const;
  Index3D next_D(int l) const;

  Index3D next_EN(int l) const;
  Index3D next_WN(int l) const;
  Index3D next_ES(int l) const;
  Index3D next_WS(int l) const;

  Index3D next_ET(int l) const;
  Index3D next_WT(int l) const;
  Index3D next_ED(int l) const;
  Index3D next_WD(int l) const;

  Index3D next_NT(int l) const;
  Index3D next_ST(int l) const;
  Index3D next_ND(int l) const;
  Index3D next_SD(int l) const;

  Index3D next_ENT(int l) const;
  Index3D next_WNT(int l) const;
  Index3D next_EST(int l) const;
  Index3D next_WST(int l) const;
  Index3D next_END(int l) const;
  Index3D next_WND(int l) const;
  Index3D next_ESD(int l) const;
  Index3D next_WSD(int l) const;

  void Print() const;

  Index3D corner_of_edge_index(dir1D d);

  Index3D neighbour(dir_sons) const;
  Index3D neighbour_non_periodic(dir_sons) const;
  bool    Is_non_periodic() const;

  inline Index3D neighbour_EW(dir1D);
  inline Index3D neighbour_NS(dir1D);
  inline Index3D neighbour_TD(dir1D);
  D3vector coordinate() const;
  double coordinate_x() const;
  double coordinate_y() const;
  double coordinate_z() const;

  bool Cell_index();  // is index a  cell index?
  bool Edge_index();  // is index an edge index?

  dir_3D direction(Index3D I);

  inline Index1D I_x() const { return ind_x; }
  inline Index1D I_y() const { return ind_y; }
  inline Index1D I_z() const { return ind_z; }

  inline Index1D I_q(main_dir_3D i) const { 
    if(i==WE_dir) return ind_x.index;
    if(i==SN_dir) return ind_y.index;
    return ind_z.index;
  }

  // functions which give back the edge or face of a cell
  // with center point of type Index3D
  IndexEdge edge(Edges_cell ec);
  Index3D I_edge(Edges_cell ec);
  IndexSurface surface(dir_3D d);

  dir_sons color(int level);

  bool Index_for_cell();
  int  operator!=(Index3D I) const;
  int  operator==(Index3D I) const;
  bool operator<=(Index3D I) const;  // point finer or equal as this
  bool operator< (Index3D I) const;  // I_this subset [I-h,I+h[ x ...
  bool operator<<(Index3D I) const;  // I_this subset [I-h,I+h] x ...
};

class IndexBo {
 public:
  IndexBo(Index3D Ind, dir_3D dir) : I(Ind), d(dir) {};
  IndexBo() : d(Wdir) {};
  Index3D I;
  dir_3D d;
};

class IndexEdge {
 friend class Index3D;
 public:
  IndexEdge();
  inline IndexEdge(Index3D in, dir_3D di);

  Index3D corner(dir1D di);
  dir_3D  opposite_direction(dir1D di);
  dir_3D           direction(dir1D di);
  Index3D I() { return ind; };
  bool operator>>(Index3D I) const;        // I subset Edge
  bool operator<=(Index3D I) const;  // point finer or equal edge
                                     // middle point
 private:
  Index3D ind;
  dir_3D d;
};

class IndexSurface {
 public:
  IndexSurface(Index3D I, dir_3D d) : ind(I), dir(d) {};
  IndexSurface() : dir(Wdir) {};
  Index3D corner(dir2D_sons cor);    // corner of face
  Index3D edge(dir_2D edg);          // edges of face
  dir_3D  direction(dir_2D d);       // 3D direction of 2D-face-direction
  Index3D I();                       // index of center point
  bool operator>>(Index3D I) const;  // I subset Surface
  bool operator<=(Index3D I) const;  // point finer or equal surface 
                                     // middle point
 private:
  Index3D ind;
  dir_3D dir;
};

// function for backtransfomation
// level ist Tiefe der Nachbarn der zu suchenden Zelle mit Index I
// in der x enthalten ist. Rueckgabe Index der Zelle.
inline Index1D back_coordinate(double x, int level); 
inline Index3D back_coordinate(D3vector v, int level); 


// inline implementation of some member functions

// coding of 1D index:
// 0 ist point 0
// ... + (2*c-1)*0.25 + (2*b-1)*0.5 + (2*a-1)*1.0   =  abc... = p
// where a=1 or p=0

//  Index1D:
//-----------

inline long_int Index1D::cart_index(int l) const {
  long_int sum, ind;
  int level;
  sum  = 0;
  if(developer_version) {
    if(l<Tiefe()) cout << " Error in Index1D::cart_index!! " << endl;
  }  
  level= l-Tiefe();
  for(ind = index;ind>0;ind = ind >> 1) {
    sum = sum + Zweipotenz(level)*(2*(ind&1)-1);
    level=level+1;
  };
  return sum;
}



inline Index1D Index1D::cut_level(int l) const   {
  long_int ind, v;
  v=1;
  v = v << l;
  for(ind = index;ind>v;ind = ind >> 1) {};
  return Index1D(ind);
}


inline Index1D Index1D::neighbour(dir1D d) const   {
  long_int ind;

#ifdef PERIODIC
  for(ind = index;(ind&1)==d;ind = ind >> 1) {};
  if(developer_version) {
    if(ind==3) cout << " error 2 in Index1D::neighbour " << endl;
  }
  ind = ind >> 1;
  if(ind==0) return Index1D(1); 
  return Index1D(ind);
#else 
  if(d==Ld && index==0) return 3; // this means outside of [0,1]
  for(ind = index;(ind&1)==d;ind = ind >> 1) {};
  return Index1D(ind >> 1);
#endif
}


inline Index1D Index1D::neighbour_non_periodic(dir1D d) const   {
  long_int ind;
  if(d==Ld && index==0) return 3; // this means outside of [0,1]
  for(ind = index;(ind&1)==d;ind = ind >> 1) {};
  return Index1D(ind >> 1);
}

inline bool Index1D::Is_non_periodic() const {
// geaendert 22.5.03: index==0 statt index!=0
  return (index==0);
}

inline int Index1D::Tiefe() const {
  int T;
  long_int ind;

  ind = index;
  for(T=-1;ind!=0;++T) ind = ind >> 1;
  return T;
}

inline Index1D Index1D::next(dir1D d, int l) const {
#ifdef PERIODIC
    int i;
    long_int ind;
  if(developer_version) {
    if(index==0) cout << " error 1 in Index1D::next(dir1D d, int l) " << endl;
  }
  l = l - Tiefe();
  if(l==0) return neighbour(d);
  ind = (index << 1) + d ; 
  if(ind==3) ind=2;
  for(i=1;i<l;++i) ind = (ind << 1) + ((~d)&1) ; 
  return Index1D(ind);
#else
  int i;
  long_int ind;
  l = l - Tiefe();
  if(d==Ld && index==0) return 3; // this means outside of [0,1]
  if(l==0) return neighbour(d);
  ind = (index << 1) + d ; 
  for(i=1;i<l;++i) ind = (ind << 1) + ((~d)&1) ; 
  return Index1D(ind);
#endif
}

inline Index1D Index1D::next(Ort1D ort, int l) const {
#ifdef PERIODIC
  int i, dir_of_ort;
  long_int ind;

  if(developer_version) {
    if(index==0) cout << " error 1 in Index1D::next(dir1D d, int l) " << endl;
  }
  if(ort==MOrt) return Index1D(index);
  dir_of_ort = ort/2;

  l = l - Tiefe();
  if(l==0) return neighbour((dir1D)dir_of_ort);
  ind = (index << 1) + dir_of_ort ; 
  if(ind==3) ind=2;
  for(i=1;i<l;++i) ind = (ind << 1) + ((~dir_of_ort)&1) ; 
  return Index1D(ind);
#else
  int i, dir_of_ort;
  long_int ind;

  if(ort==1) return Index1D(index);
  dir_of_ort = ort/2;

  if(ort==LOrt && index==0) return 3; // this means outside of [0,1]
  l = l - Tiefe();
  if(l==0) return neighbour((dir1D)dir_of_ort);
  ind = (index << 1) + dir_of_ort ; 
  for(i=1;i<l;++i) ind = (ind << 1) + ((~dir_of_ort)&1) ; 
  return Index1D(ind);
#endif
}

inline Index1D Index1D::son(dir1D d) const {
#ifdef PERIODIC
  if(index==1) return Index1D(2); 
#endif
  return Index1D( (index << 1) + d ); 
};

inline Index1D Index1D::father() {
  return Index1D(index >> 1);
};

inline double Index1D::h_loc() const {
  static long_int ind;
  static double h;
  ind = index;
  h = 2.0;
  for(h = 2.0;ind!=0;h = h * 0.5) {
    ind = ind >> 1;
  }
  return h;
}

inline double Index1D::coordinate() const {
  int i;
  long_int ind;

  static double h, x;
  ind = index;
  h = h_loc();
  x = 0.0;
  for(i=Tiefe();i>=0;--i) {
    x = x + h * (2*(ind&1)-1);
    h = h * 2.0;
    ind = ind >> 1;
  }
  return x;
}

inline int Index1D::operator!=(Index1D I) const {
  return index!=I.index;
}
inline int Index1D::operator==(Index1D I) const {
  return index==I.index;
}


// I_this subset [I-h,I+h[
inline bool Index1D::operator<(Index1D I) const {
  long_int ind;
#ifdef PERIODIC
  if(I.index==1) return true;
#endif

  if(index==I.neighbour(Ld).index) return true;
  ind = index;
  while(I.index<ind) ind = ind>>1;
  return I.index==ind;
}

// I_this subset [I-h,I+h]
inline bool Index1D::operator<<(Index1D I) const {
  long_int ind;

#ifdef PERIODIC
  if(I.index==1) return true;
#else
  if(I.index==0 || I.index==1) return true;
#endif

  if(index==I.neighbour(Ld).index) return true;
  if(index==I.neighbour(Rd).index) return true;
  ind = index;
  while(I.index<ind) ind = ind>>1;
  return I.index==ind;
}

// point finer or equal as this
inline bool Index1D::operator<=(Index1D I) const {
  long_int ind;
  ind = I.index;
  while(index<ind) ind = ind>>1;
  return index==ind;
}

//  Index3D:
//----------
inline void Index3D::Print() const {
  ind_x.Print();
  cout << ", ";
  ind_y.Print();
  cout << ", ";
  ind_z.Print();
}

inline Index3D Index3D::corner_of_edge_index(dir1D d) {
  if(ind_x.Tiefe() > ind_y.Tiefe() &&
     ind_x.Tiefe() > ind_z.Tiefe()) return neighbour_EW(d);
  if(ind_y.Tiefe() > ind_x.Tiefe() &&
     ind_y.Tiefe() > ind_z.Tiefe()) return neighbour_NS(d);
  // if(ind_z.Tiefe() > ind_x.Tiefe() &&
  // ind_z.Tiefe() > ind_y.Tiefe()) 
  return neighbour_TD(d);
};

inline Index3D::Index3D(long_int intx,long_int inty,long_int intz) {
  ind_x = Index1D(intx); 
  ind_y = Index1D(inty); 
  ind_z = Index1D(intz); 
};

inline Index3D::Index3D(Index1D indx,Index1D indy,Index1D indz) {
  ind_x=indx; ind_y=indy; ind_z=indz; 
};


inline Index3D Index3D::son(dir_sons d) const {
  return Index3D(ind_x.son((dir1D)(d&1)),
                 ind_y.son((dir1D)((d >> 1)&1)),
                 ind_z.son((dir1D)((d >> 2)&1)));
};

inline Index3D Index3D::father() {
  return Index3D(ind_x.father(),
                 ind_y.father(),
                 ind_z.father());
};

inline Index3D Index3D::son_EW(dir1D d) {
  return Index3D(ind_x.son(d),ind_y,ind_z); 
};

inline Index3D Index3D::son_NS(dir1D d) {
  return Index3D(ind_x,ind_y.son(d),ind_z); 
};

inline Index3D Index3D::son_TD(dir1D d) {
  return Index3D(ind_x,ind_y,ind_z.son(d)); 
};

inline Index3D Index3D::next(Ort1D ortx, Ort1D orty, Ort1D ortz, int l) const {
  return Index3D(ind_x.next(ortx,l),ind_y.next(orty,l),ind_z.next(ortz,l));
}

inline Index3D Index3D::next(dir_sons d, int l) const {
  return Index3D(ind_x.next((dir1D)(d&1),l),
                 ind_y.next((dir1D)((d >> 1)&1),l),
                 ind_z.next((dir1D)((d >> 2)&1),l));
}
inline Index3D Index3D::next(dir_3D d, int l) const {
  if(d>Ndir) {
    return next_TD(OneDpart(d),l);
  }
  if(d<Sdir) {
    return next_EW(OneDpart(d),l);
  }
  return next_NS(OneDpart(d),l);
}

inline Index3D Index3D::next(Edges_cell ed, int t) const {
  if(ed==SDed) return next_SD(t);
  if(ed==NDed) return next_ND(t);
  if(ed==STed) return next_ST(t);
  if(ed==NTed) return next_NT(t);

  if(ed==EDed) return next_ED(t);
  if(ed==WDed) return next_WD(t);
  if(ed==ETed) return next_ET(t);
  if(ed==WTed) return next_WT(t);

  if(ed==NEed) return next_EN(t);
  if(ed==NWed) return next_WN(t);
  if(ed==SEed) return next_ES(t);
  // if(ed==SWed) 
  return next_WS(t);
}

inline Index3D Index3D::next_EW(dir1D d, int l) const  {
  return Index3D(ind_x.next(d,l),ind_y,ind_z);
}
inline Index3D Index3D::next_NS(dir1D d, int l) const {
  return Index3D(ind_x,ind_y.next(d,l),ind_z);
}
inline Index3D Index3D::next_TD(dir1D d, int l) const  {
  return Index3D(ind_x,ind_y,ind_z.next(d,l));
}

inline Index3D Index3D::next_N(int l) const {  return next_NS(Rd,l); };
inline Index3D Index3D::next_S(int l) const {  return next_NS(Ld,l); };
inline Index3D Index3D::next_E(int l) const {  return next_EW(Rd,l); };
inline Index3D Index3D::next_W(int l) const {  return next_EW(Ld,l); };
inline Index3D Index3D::next_T(int l) const {  return next_TD(Rd,l); };
inline Index3D Index3D::next_D(int l) const {  return next_TD(Ld,l); };

inline Index3D Index3D::next_EN(int l) const { 
  return next_EW(Rd,l).next_NS(Rd,l); };
inline Index3D Index3D::next_WN(int l) const {
  return next_EW(Ld,l).next_NS(Rd,l); };
inline Index3D Index3D::next_ES(int l) const {
  return next_EW(Rd,l).next_NS(Ld,l); };
inline Index3D Index3D::next_WS(int l) const {
  return next_EW(Ld,l).next_NS(Ld,l); };

inline Index3D Index3D::next_ET(int l) const {
  return next_EW(Rd,l).next_TD(Rd,l); };
inline Index3D Index3D::next_WT(int l) const {
  return next_EW(Ld,l).next_TD(Rd,l); };
inline Index3D Index3D::next_ED(int l) const {
  return next_EW(Rd,l).next_TD(Ld,l); };
inline Index3D Index3D::next_WD(int l) const {
  return next_EW(Ld,l).next_TD(Ld,l); };

inline Index3D Index3D::next_NT(int l) const {
  return next_NS(Rd,l).next_TD(Rd,l); };
inline Index3D Index3D::next_ST(int l) const {
  return next_NS(Ld,l).next_TD(Rd,l); };
inline Index3D Index3D::next_ND(int l) const {
  return next_NS(Rd,l).next_TD(Ld,l); };
inline Index3D Index3D::next_SD(int l) const {
  return next_NS(Ld,l).next_TD(Ld,l); };

inline Index3D Index3D::next_ENT(int l) const {
  return next_EW(Rd,l).next_NS(Rd,l).next_TD(Rd,l); };
inline Index3D Index3D::next_WNT(int l) const {
  return next_EW(Ld,l).next_NS(Rd,l).next_TD(Rd,l); };
inline Index3D Index3D::next_EST(int l) const {
  return next_EW(Rd,l).next_NS(Ld,l).next_TD(Rd,l); };
inline Index3D Index3D::next_WST(int l) const {
  return next_EW(Ld,l).next_NS(Ld,l).next_TD(Rd,l); };
inline Index3D Index3D::next_END(int l) const {
  return next_EW(Rd,l).next_NS(Rd,l).next_TD(Ld,l); };
inline Index3D Index3D::next_WND(int l) const {
  return next_EW(Ld,l).next_NS(Rd,l).next_TD(Ld,l); };
inline Index3D Index3D::next_ESD(int l) const {
  return next_EW(Rd,l).next_NS(Ld,l).next_TD(Ld,l); };
inline Index3D Index3D::next_WSD(int l) const {
  return next_EW(Ld,l).next_NS(Ld,l).next_TD(Ld,l); };

inline dir_3D Index3D::direction(Index3D I) {
  if(ind_x != I.ind_x) {
    if(ind_x.coordinate() < I.ind_x.coordinate()) return Edir;
    else return Wdir;
  }
  if(ind_y != I.ind_y) {
    if(ind_y.coordinate() < I.ind_y.coordinate()) return Ndir;
    else return Sdir;
  }
  // if(ind_z != I.ind_z) {
  if(ind_z.coordinate() < I.ind_z.coordinate()) return Tdir;
  else return Ddir;
  // }
}

inline Index3D Index3D::neighbour(dir_sons d) const {
  return Index3D(ind_x.neighbour((dir1D)(d&1)),
                 ind_y.neighbour((dir1D)((d >> 1)&1)),
                 ind_z.neighbour((dir1D)((d >> 2)&1)));
};

inline bool Index3D::Is_non_periodic() const {
 // geaendert 22.5.03: || statt &&
   return (ind_x.Is_non_periodic() ||
          ind_y.Is_non_periodic() ||
          ind_z.Is_non_periodic());
}

inline Index3D Index3D::neighbour_non_periodic(dir_sons d) const {
  return Index3D(ind_x.neighbour_non_periodic((dir1D)(d&1)),
                 ind_y.neighbour_non_periodic((dir1D)((d >> 1)&1)),
                 ind_z.neighbour_non_periodic((dir1D)((d >> 2)&1)));
};

inline Index3D Index3D::neighbour_EW(dir1D d) {
  return Index3D(ind_x.neighbour(d),ind_y,ind_z); 
};

inline Index3D Index3D::neighbour_NS(dir1D d) {
  return Index3D(ind_x,ind_y.neighbour(d),ind_z); 
};

inline Index3D Index3D::neighbour_TD(dir1D d) {
  return Index3D(ind_x,ind_y,ind_z.neighbour(d)); 
};

inline int Index3D::Tiefe() const {
  return MAX(ind_x.Tiefe(),ind_y.Tiefe(),ind_z.Tiefe());
};

inline D3vector Index3D::coordinate() const {
  return D3vector(ind_x.coordinate(),ind_y.coordinate(),ind_z.coordinate());
};

inline double Index3D::coordinate_x() const {
  return ind_x.coordinate();
};

inline double Index3D::coordinate_y() const {
  return ind_y.coordinate();
};

inline double Index3D::coordinate_z() const {
  return ind_z.coordinate();
};

inline bool Index3D::Cell_index() {
  return (ind_x.Tiefe()==ind_y.Tiefe()) && (ind_y.Tiefe()==ind_z.Tiefe());
}

inline dir_sons Index3D::color(int level) {
  return (dir_sons)(((int)(ind_x.Tiefe()==level))
                    + ((int)(ind_y.Tiefe()==level)) * 2
                    + ((int)(ind_z.Tiefe()==level)) * 4);
}

inline bool Index3D::Edge_index() {
  if((ind_x.Tiefe()>ind_y.Tiefe()) && (ind_x.Tiefe()>ind_z.Tiefe()))
    return true;
  if((ind_y.Tiefe()>ind_x.Tiefe()) && (ind_y.Tiefe()>ind_z.Tiefe()))
    return true;
  if((ind_z.Tiefe()>ind_x.Tiefe()) && (ind_z.Tiefe()>ind_y.Tiefe()))
    return true;
  return false;
}

inline int Index3D::operator!=(Index3D I) const {
  return ind_x!=I.ind_x || ind_y!=I.ind_y || ind_z!=I.ind_z;
}
inline int Index3D::operator==(Index3D I) const {
  return ind_x==I.ind_x && ind_y==I.ind_y && ind_z==I.ind_z;
}
inline bool Index3D::operator<(Index3D I) const {
  return ind_x<I.ind_x && ind_y<I.ind_y && ind_z<I.ind_z;
}
inline bool Index3D::operator<=(Index3D I) const {
  return ind_x<=I.ind_x && ind_y<=I.ind_y && ind_z<=I.ind_z;
}
inline bool Index3D::operator<<(Index3D I) const {
  return ind_x<<I.ind_x && ind_y<<I.ind_y && ind_z<<I.ind_z;
}

inline bool Index3D::Index_for_cell() {
  return (ind_x.Tiefe() == ind_y.Tiefe() &&
          ind_y.Tiefe() == ind_z.Tiefe());
}

// concept: ec has 4 bits. first two bits are direction of the edge
//          and the other two bits the location of the orthogonal
inline IndexEdge Index3D::edge(Edges_cell ec) {
  enum dir_3D dir;
  dir = (dir_3D) ((ec >> 2) << 1);
  if(dir<Sdir)
    return IndexEdge(Index3D(ind_x,
                             ind_y.neighbour((dir1D)(ec&1)),
                             ind_z.neighbour((dir1D)((ec>>1)&1))),Wdir);
  if(dir>Ndir)
    return IndexEdge(Index3D(ind_x.neighbour((dir1D)(ec&1)),
                             ind_y.neighbour((dir1D)((ec>>1)&1)),
                             ind_z),Ddir);

  return IndexEdge(Index3D(ind_x.neighbour((dir1D)(ec&1)),
                           ind_y,
                           ind_z.neighbour((dir1D)((ec>>1)&1))),Sdir);
};

inline Index3D Index3D::I_edge(Edges_cell ec) {
  dir_3D dir;
  dir = (dir_3D) ((ec >> 2) << 1);
  if(dir<Sdir)
    return Index3D(ind_x,
                   ind_y.neighbour((dir1D)(ec&1)),
                   ind_z.neighbour((dir1D)((ec>>1)&1)));
  if(dir>Ndir)
    return Index3D(ind_x.neighbour((dir1D)(ec&1)),
                   ind_y.neighbour((dir1D)((ec>>1)&1)),
                   ind_z);

  return Index3D(ind_x.neighbour((dir1D)(ec&1)),
                 ind_y,
                 ind_z.neighbour((dir1D)((ec>>1)&1)));
};

inline IndexSurface Index3D::surface(dir_3D d) {
  if(d<Sdir) 
    return IndexSurface(Index3D(ind_x.neighbour((dir1D)(d&1)),ind_y,ind_z),d);
  if(d>Ndir) 
    return IndexSurface(Index3D(ind_x,ind_y,ind_z.neighbour((dir1D)(d&1))),d);
  return IndexSurface(Index3D(ind_x,ind_y.neighbour((dir1D)(d&1)),ind_z),d);
};




//  IndexEdge:
//------------
inline IndexEdge::IndexEdge() : d(Wdir) {};
inline IndexEdge::IndexEdge(Index3D in, dir_3D di) : ind(in), d(di) {};

inline Index3D IndexEdge::corner(dir1D di) {
  if(d>Ndir) return ind.neighbour_TD(di);
  if(d<Sdir) return ind.neighbour_EW(di);
  return ind.neighbour_NS(di);
}

inline dir_3D IndexEdge::direction(dir1D di) {
  return (dir_3D)(d + di);
}

inline dir_3D IndexEdge::opposite_direction(dir1D di) {
  return (dir_3D)(d + ((~di)&1));
}

inline bool IndexEdge::operator>>(Index3D I) const {
  if(d>Ndir)
    return I.I_z()<<ind.I_z() && ind.I_y()==I.I_y() && ind.I_x()==I.I_x();
  if(d<Sdir)
    return ind.I_z()==I.I_z() && ind.I_y()==I.I_y() && I.I_x()<<ind.I_x();
  return   ind.I_z()==I.I_z() && I.I_y()<<ind.I_y() && ind.I_x()==I.I_x();
}

inline bool IndexEdge::operator<=(Index3D I) const {
  if(d>Ndir)
    return ind.I_z()<=I.I_z() && ind.I_y()==I.I_y() && ind.I_x()==I.I_x();
  if(d<Sdir)
    return ind.I_z()==I.I_z() && ind.I_y()==I.I_y() && ind.I_x()<=I.I_x();
  return   ind.I_z()==I.I_z() && ind.I_y()<=I.I_y() && ind.I_x()==I.I_x();
}


//  IndexSurface:
//---------------

inline bool IndexSurface::operator>>(Index3D I) const {
  return (I<<ind) && 
    (ind.I_q(Main_direction(dir))==I.I_q(Main_direction(dir)));
}

inline bool IndexSurface::operator<=(Index3D I) const {
  return (ind<=I) && 
    (ind.I_q(Main_direction(dir))==I.I_q(Main_direction(dir)));
}


inline Index3D IndexSurface::corner(dir2D_sons cor) {
  if(dir < Sdir)
    return ind.neighbour_NS((dir1D)(cor&1)).neighbour_TD((dir1D)(cor>>1));
  if(dir > Ndir)
    return ind.neighbour_EW((dir1D)(cor&1)).neighbour_NS((dir1D)(cor>>1));
  return ind.neighbour_EW((dir1D)(cor&1)).neighbour_TD((dir1D)(cor>>1));

};

inline Index3D IndexSurface::edge(dir_2D edg) {
  if(dir < Sdir) {
    if(edg < Sdir2D)
      return ind.neighbour_NS((dir1D)(edg&1));
    else
      return ind.neighbour_TD((dir1D)(edg&1));
  }
  if(dir > Ndir) {
    if(edg < Sdir2D)
      return ind.neighbour_EW((dir1D)(edg&1));
    else
      return ind.neighbour_NS((dir1D)(edg&1));
  }
  if(edg < Sdir2D)
    return ind.neighbour_EW((dir1D)(edg&1));
  else
    return ind.neighbour_TD((dir1D)(edg&1));
};

inline dir_3D IndexSurface::direction(dir_2D d) {
  if(dir < Sdir) return (dir_3D)(d+2);
  if(dir > Ndir) return (dir_3D)d;
  if(d < Sdir2D) return (dir_3D)d;
  else           return (dir_3D)(d+2);
}

inline Index3D IndexSurface::I() {
  return ind;
}



// back transformation function
/*
inline double Index1D::coordinate() const {
  int i;
  long_int ind;

  static double h, x;
  ind = index;
  h = h_loc();
  x = 0.0;
  for(i=Tiefe();i>=0;--i) {
    x = x + h * (2*(ind&1)-1);
    h = h * 2.0;
    ind = ind >> 1;
  }
  return x;
}
*/



// Kann man noch linear optimieren
// so ist es quadratisch
inline Index1D back_coordinate(double x, int level) {
  Index1D I(2);
  double xx;

  for(int t=1;t<=level;++t) {
    xx = I.coordinate();
    if(x<xx) 
      I = I.son(Ld);
    else
      I = I.son(Rd);
  }
  return I;
};


inline Index3D back_coordinate(D3vector v, int level) {
  return Index3D(back_coordinate(v.x,level),
                 back_coordinate(v.y,level),
                 back_coordinate(v.z,level));
} 









#endif
       

---