src/expde/grid/gridbase.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
//

// ------------------------------------------------------------
//
// gridbase.h
// bindet unten auch hash.h ein
//
// ------------------------------------------------------------

#ifndef GRIDBASE_H_
#define GRIDBASE_H_
         
#define Point_hashtable2 BoCell
#define Point_hashtable3 Bo2Point
#define Point_hashtable4 Data
class Point_hashtable0;
class Point_hashtable1;
class BoCell;
class Bo2Point;
class Data;
class Grid;
class P_boundary_tet;

#include <set>
#include <string>
#include <vector>
using namespace std;


class Grid_base : public Storage_manager, public Parallel_Info  {
 friend class Point_hashtable3;
 public:
  // Information ueber Hashtables
  void Information();


#ifdef WITH_BRICK
#define DIMENSION 3
  
  struct BInfoT {
    double localXMin[DIMENSION];
    double localXMax[DIMENSION];
    double eps;
    int rank;
  };

  struct eqBInfoT {
    bool operator()(const BInfoT b1, const BInfoT b2) const
    {
      return ( ABS(b1.localXMin[0] - b2.localXMin[0]) > b1.eps ||
               ABS(b1.localXMin[1] - b2.localXMin[1]) > b1.eps ||
               ABS(b1.localXMin[2] - b2.localXMin[2]) > b1.eps ||
               ABS(b1.localXMax[0] - b2.localXMax[0]) > b1.eps ||
               ABS(b1.localXMax[1] - b2.localXMax[1]) > b1.eps ||
               ABS(b1.localXMax[2] - b2.localXMax[2]) > b1.eps  );
    }
  };
  
  
  typedef set<BInfoT,eqBInfoT> Set_of_Processors;
  
  Set_of_Processors brickInfo_m;
  
  void getBrickInfo();   
  void boundbox(double localXMin[DIMENSION], double localXMax[DIMENSION]);
  void localbox(double localXMin[DIMENSION], double localXMax[DIMENSION], int node);


  inline int isInBrick(double x[DIMENSION]) {
      
    /*
      If x \in \Omega_k return rank of the processor
      which owns \Omega_k
      
      This must be optimized !!
    */
      Set_of_Processors::const_iterator i;
      for(i=brickInfo_m.begin();i!=brickInfo_m.end();i++) {
          if (((*i).localXMin[0] <= x[0]) &&
              ((*i).localXMax[0] >  x[0]) &&
              ((*i).localXMin[1] <= x[1]) &&
              ((*i).localXMax[1] >  x[1]) &&
              ((*i).localXMin[2] <= x[2]) &&
              ((*i).localXMax[2] >  x[2])
              )
              return (*i).rank;
      }
      return -1;  // coordinate does not belong to a processor
  }
#endif // WITH_BRICK


  // Test functions
  void Print_test_AVS(ofstream *Datei, int type);
  void Test(int type);
  void Print_Test_Cell(ofstream *Datei);
  void Print_Test_all_Cell(ofstream *Datei, int level);

  void Test_just_this();
  void Test_just_this2(int bla);


  // Visualization 
  // AVS:
  void Print_Variable_AVS(ofstream *Datei, int number_var);
  void Print_Variable_AVS_coarse(ofstream *Datei, int number_var, int level);
  void Print_Cell_Variable_AVS(ofstream *Datei, int number_var);
  void Print_Cell_Variable_AVS_parallel(ofstream *Datei, int number_var);
  void Print_Variable_AVS_parallel(ofstream *Datei, int number_var);
  void Print_Variable_AVS_parallel(ofstream *Datei, int number_var, int level);
  void Print_Variable_AVS_parallel(ofstream *Datei, int number_var,
                                   int number_varb, int number_varc);
  void Print_surface_Variable_AVS_parallel(ofstream *Datei, int number_var);

  void Print_Cell_Variable_AVS_moved(ofstream *Datei, int number_var,
                                     int var_a, int var_b, int var_c);
  void Print_Variable_AVS_moved(ofstream *Datei, int number_var,
                                int var_a, int var_b, int var_c);
  void Print_Variable_AVS_moved_parallel(ofstream *Datei, int number_var,
                                         int var_a, int var_b, int var_c);
  void Print_Variable_AVS(ofstream *Datei, int var_a, int var_b, int var_c);
  void Print_Domain_AVS(ofstream *Datei);
  void Print_processes_UCD(ofstream *Datei);
  void Print_region_processes_UCD(ofstream *Datei);
  void Print_surface_processes_UCD(ofstream *Datei);
 
  // Gnuplot
  void Print_Grid_Gnuplot(ofstream *Datei);
  void Print_Grid_Gnuplot_moved(ofstream *Datei,
                                int var_a, int var_b, int var_c);

  // OpenDx
  void Print_Variable_OpenDx(ofstream *Datei, int number_var);
  void Print_Variable_OpenDx_parallel(ofstream *Datei, int number_var);
  void Print_Variable_OpenDx_moved(ofstream *Datei, int number_var,
                                   int var_a, int var_b, int var_c);
  // Special, not for everybody
  void Print_Special(ofstream *File, int number_var, 
                   int num_X_DISP, int num_Y_DISP, int num_Z_DISP);
  void Print_Special_streched(ofstream *File, int number_var,int num_z_correct,
                   int num_X_DISP, int num_Y_DISP, int num_Z_DISP);
  void Print_Special_streched(ofstream *File, int number_var, 
                   int num_x_coor, int num_y_coor, int num_z_coor, 
                   int num_X_DISP, int num_Y_DISP, int num_Z_DISP);
  void Print_Special_streched_half(ofstream *File, int number_var, 
                   int num_x_coor, int num_y_coor, int num_z_coor, 
                   int num_X_DISP, int num_Y_DISP, int num_Z_DISP);
  void Print_Special_streched_quad(ofstream *File, int number_var, 
                   int num_x_coor, int num_y_coor, int num_z_coor, 
                   int num_X_DISP, int num_Y_DISP, int num_Z_DISP);
  void Print_Special_streched_quad(ofstream *File, double hlen,
                   int number_var, 
                   int num_x_coor, int num_y_coor, int num_z_coor, 
                   int num_X_DISP, int num_Y_DISP, int num_Z_DISP,
                                   int num_heat);
  void Print_Special_streched_X_Y(ofstream *, 
                                  int, int, int, 
                                  int, int, int, int, int, int, int, int);


  void Pr_Var_AVS_mo_boundary(ofstream *Datei, int number_var, //Output of
                              int var_a, int var_b, int var_c);// boundary
  void Pr_Var_AVS_mo_surface(ofstream *Datei, int number_var, //Output of
                              int var_a, int var_b, int var_c);// surface
                                               // v Output of half zylinder
  void Pr_Var_AVS_mo_surface_half(ofstream *Datei, int number_var, //Output of
                              int var_a, int var_b, int var_c);// surface
  void Pr_Var_AVS_mo_surface_transv(ofstream *Datei, int number_var,//Output of
                      int var_a, int var_b, int var_c, double plane);// surface
  void Print_Variable_AVS_moved_half(ofstream *Datei, int number_var,
                                int var_a, int var_b, int var_c);

  // Meshsize and level
  //------------------------------
  int Max_level() const;                        // maximal level of grid
  int Min_level() const;                        // minimal level of grid
  double finest_mesh_size() const;              // finest_mesh_size calculated 
  double Give_finest_mesh_size() const          // finest_mesh_size stored
    { return finest_meshsize; } 

  inline double H_mesh() const;     // Maschenweite des umschreibenen Wuerfels
  inline D3vector Give_A() const;   // Anfangspunkt des umschreibenen Wuerfels
  double h_min_for_boundary_points; // used in 
                                    // Point_hashtable3::Point_hashtable3

  //Labels at the boundary:
  //-----------------------
  // a) at boundary points
 // true means here run at that points (Neumann boundary condition)
  bool Give_label_bo(Index3D I, dir_3D d,int num) const; 
  void Put_label_bo(bool (*Formula)(double x,double y,double z),int num);
  void Put_label_bo_complementary(int num_org, int num);
  void Put_label_bo(int num_var,int num);
  void Put_union_label(int num, int numa, int numb);
  // wird noch nicht verwendet:
  void Put_label_bo(bool lab, Index3D I, dir_3D d,int num) const; 
  // b) at multigrid points
  // true means here Dirichlet boundary condition
  bool Give_label_bo_D(Index3D I, int num, int level) const; // has point I
                                              // Dirichlet-condition on level
  void Restrict_label_bo(int num);
  void Set_label_bo_mg(Index3D Ind, int num, bool label, int level); // point I
                             // has Dirichlet-condition on level if label==true
  void Interpolation_label_bo(Index3D I, int level, int num, bool Dirichlet);

  // Funktionen um zu den Datenspeichern zu gelangen:
  //----------------------------------------------------
  // a) for points 
  // Gib Zeiger auf Datenspeicher im Inneren:
  double* Give_variable(Index3D I, int v_e) const;  
  double* Give_variable_slow(const Index3D I, int v_e) const; // Wie oben: 
                                             // Gib Zeiger auf Datenspeicher, 
                                             // aber NULL falls nicht existiert
                                             // daher etwas langsamer
  // Gib Zeiger auf Datenspeicher bei Randpunkten
  double* Give_variable(Index3D I, dir_3D d) const; 
  double  Give_h(Index3D I, dir_3D d) const;         // zugehoeriger Abstand h
  // Gib Zeiger auf  Datenspeicher bei Randzellfreiheitsgraden
  double* Give_variable_cellpoi(Index3D I) const;

  // b) for cells
  // Gib Zeiger auf Datenspeicher fuer Stencils am Rand
  double* Give_bo_stencil(Index3D I, int num_stencil) const; 
  // Gib Zeiger auf Datenspeicher fuer allgemeine Stencils
  double* Give_stencil(Index3D I, int num_stencil) const;
  // Gib Zeiger auf Datenspeicher bei Variablen auf Zellen
  double* Give_cell_variable(Index3D I) const;
  // Gib Zeiger auf Datenspeicher fuer Stencils, Zellen, ...
  // Wird in commbos.cc verwendet
  double* Give_pointer_stencil(Index3D I) const; 

  // transform computational coordinate in physical coordinate
  //------------------------------------------------------------
  inline D3vector transform_coord(const D3vector V) const;
  inline D3vector back_transform_coord(const D3vector V) const;
  inline D3vector transform_coord(const Index3D  I) const;
  inline double   transform_coordX(const Index3D I) const;
  inline double   transform_coordY(const Index3D I) const;
  inline double   transform_coordZ(const Index3D I) const;
  inline D3vector transform_coord(const Index3D I, const dir_3D d) const;
  inline D3vector transform_coord_org_shift(const Index3D I, const dir_3D d,
                                            const double shift) const;
  inline D3vector transform_coord_cellpoi(const Index3D I) const;
  inline double   transform_coordX_cellpoi(const Index3D I) const;
  inline double   transform_coordY_cellpoi(const Index3D I) const;
  inline double   transform_coordZ_cellpoi(const Index3D I) const;

  // Funktionen die vermutlich nur zur Fehlersuche verwendet werden:
  // Ausdrucken was in den Hashtables steht:
  void Print_hashtable0();
  void Print_hashtable1();
  void Print_hashtable2();
  void Print_hashtable3();
  void Print_hashtable4();

  // Ausdrucken der Zellen:
  void Print_cell_info();

  // Informationen ueber die Punkte und Kanten und Zellen
  //-------------------------------------------------------
  bool Exists_Point(Index3D I) const;              // Existiert Punkt? 
  bool Point_in_domain(Index3D I) const;           // Punkt im Gebiet?
  bool Exists_Bo_Point(IndexBo I) const;           // Randpunkt im Gebiet?
  bool Exists_Bo_Point(Index3D I, dir_3D d) const; // Randpunkt im Gebiet?
  Pointtype Give_type(Index3D I) const;            // Typ des Punkt?
  Pointtype Give_global_type(Index3D I) const;     // globaler Typ des Punkt?

  bool Grid_point_on_finest_level(Index3D I) const;  // wird in mg.h verwendet
  Edgetype  Give_edge_typ(Index3D I);              // Typ der Kante?
  Celltype  Give_cell_typ(Index3D I) const;        // Typ der Zelle?
  bool Exists_Cell(Index3D I);                     // Gibt es Zelle?

  // Informationen ueber die Randzellen
  bool Is_Bo_cell(Index3D I);                   // Ist die Zelle Randzelle?
  BoCell *Give_Bo_cell(Index3D I) const;        // Gib Zeiger auf Randzelle
                                                // Falls nicht existiert: NUll

  // Alle Punkte werden durchnummeriert fuer AVS. Dies sind die Nummern:
  int Give_Nummer(Index3D I);           // Nummer eines regulaeren Punktes
  int Give_Nummer(Index3D I,dir_3D d);  // Nummer eines KantenPunktes
  int Give_Nummer_cellpoi(Index3D I);   // Nummer eines Randzellfreiheitsgrad

  // Ein einfacher Label
  void Label_true(Index3D I);                  // Label true an Index I
  bool Label_ask(Index3D I);                   // Wie ist Label an Index I ?
  void False_Label4();                         // Label4 in allen Punkten false

  // for iterate in void DResDiff_Bo<A>::Iterate_Calc_stencil(int r_bov);
  int Give_hashtable0_leng() { return hashtable0_leng; }
  Point_hashtable0** Give_hashtable0_start() {
    return hashtable0_start;
  }
  int Give_hashtable2_leng() { return hashtable2_leng; }
  Point_hashtable2** Give_hashtable2_start() {
    return hashtable2_start;
  }

  // protected:

 public:
  // auxiliary variable for later in grid.cc
  P_boundary_tet* auxiliary_P;
  int  auxiliary_number;


  // finest_meshsize and things for grid generation
  //  most things are in Parallel_Info
  double h_offset;       // offset for domain->GiveH()
  double finest_meshsize;

  // Feinster Level:
  // folgende zwei Funktionen sind eher fuer spaetere adaptive Version
  void Put_finest_level_minimal(Index3D I,    // Setze feinster moeglicher
                                int level);   // Level groesser gleich level
  int Give_finest_level(Index3D I);           // Gib feinsten moeglichen Level

  int Give_my_finest_parallel_level(Index3D I) const;  // Gib feinsten level
                                              // lokal fuer Prozessor

  void Construct_cell_points_hash0();         // Zellen von hash1 in hash0
  void Calc_type_of_edges_uniform();          // Typ der Kanten berechnen
                                              // abhaengig von domain
                                              // auf uniformen Gitter
  bool Correct_faces_uniform();              // Flaechen vom Typ d) vermeiden
  bool Correct_edge(Index3D I);              // Flae. an I vom Typ d) vermeiden
  bool Correct_face(IndexSurface I);          // Flaeche vom Typ d) vermeiden
  void Set_edge_typ(Index3D I, Edgetype typ); // Setze typ der Kante
  void Set_cell_typ(Index3D I, Celltype typ); // Setze typ der Zelle
  void Set_point_typ(Index3D I,
                     Pointtype typ);          // Setze typ des Punktes

  void Add_cell(Index3D I);                   // Fuege Zelle hinzu
  void Add_edge(Index3D I);                   // Fuege Kante hinzu
  void Add_edge_parallel(Index3D I);          // Fuege Kante hinzu
  void Add_point(Index3D I);                  // Fuege neuen Punkte hinzu
  void Set_point_typ(Index3D I,               // Setze typ des Punktes
                     int finest_level,        // and finest local level
                     Pointtype type);
  void Add_point(Index3D I, Pointtype type);  // Fuege neuen Punkte hinzu
  void Add_point_interior(Index3D I);         // Fuege inneren Punkte hinzu
  void Add_bocell(Index3D I);                 // Fuege Randzelle hinzu
  void Add_bo2point(Index3D I, dir_3D d,      // Fuege Randpunkt Klasse 2 hinzu
   int l, int number_var);                    // l ist Tiefe der zug. Zelle
  void Add_Uniform_Grid(int Max_tiefe);       // Fuege Punkte eines uniformen
                                              // Gitters hinzu
  // varebene=ve_no_multilevel: bedeutet ein Speicher ohne multilevel-Struktur
  // varebene > 0 dann ist varebene der Speicher auf dieser Ebene
  void Add_double(Index3D I, int varebene,
                  int number_var);            // Fuege Daten-Speicher hinzu
  void Add_MG_double(int varebene,            // Bei Punkten mit label4==true
  int number_var);                            // double Speicher mit varebene

  // for MG
  void Calc_multigrid_points_part1();
  void Calc_multigrid_points_part2();
  void Add_multigrid_point(Index3D I,int l);  // Fuege event. MG Punkte hinzu
  void Add_interpolation_point(Index3D I,     // Punkt muss interpoliert
                               int l);        // werden


  bool Is_Slave(Index3D I);                   // Ist dieser Punkt slave-point

  void Depth_hashtable0(int* max_depth, 
                        double* ave_depth);   // Tiefe des Hashtables
  void Depth_hashtable1(int* max_depth, 
                        double* ave_depth);   // Tiefe des Hashtables
  void Depth_hashtable2(int* max_depth, 
                        double* ave_depth);   // Tiefe des Hashtables
  void Depth_hashtable3(int* max_depth, 
                        double* ave_depth);   // Tiefe des Hashtables
  void Depth_hashtable4(int* max_depth, 
                        double* ave_depth);   // Tiefe des Hashtables


  D3vector local_coord_cellpoi(const Index3D I) const; // lokale Koordinate
                                             // eines Zellfreiheitsgrades
  Grid_base(int n_max, All_Domains* dom, 
            Grid_gen_parameters& gpara, MPI_Comm comm);
  Grid_base(double h_finest_level, All_Domains* dom, 
            Grid_gen_parameters& gpara, MPI_Comm comm);
  Grid_base(int n_max, All_Domains* dom, MPI_Comm comm);
  Grid_base(double h_finest_level, All_Domains* dom, MPI_Comm comm);
  void Grid_generation(int n_max);
  void Dummy_grid_generation();
  void Add_27Stencil(Index3D I);
  void Add_small_27Stencil(Index3D I, int level);
  void Fullfill_27Stencil();
  void Fullfill_B1B2();
  void Add_B1B2(Index3D I);
  void Add_points_of_cell(Index3D I);
  void Remove_edges();
  void Construct_points_hash1();
  void Remove_exterior_points();
  void Calc_interior_points();
  //  void Calc_slave_points();
  //  void Calc_nearb_points();
  void Calc_interior_cells_Part1();
  void Calc_interior_cells_Part2();
  void Calc_interior_cells_Part3(int level);
  void Calc_boundary_cells();
  void Calc_boundary_2points(int number_var);
  void Info_hashtable0();
  void Info_hashtable1();
  void Info_hashtable2();
  void Info_hashtable3();
  void Info_hashtable4();
  void Calc_Neumann_pro_type();
  void Remove_all_hashtables();

  // Parallel Functions
  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  void Send_boundary_cells_parallel();
  void Send_cells_parallel(int level);
  bool Send_boundary_cell_in_direction(Point_hashtable0* poi,
                                       Index3D next_index);
  bool Send_bo_with_cell_point_in_direction(dir_3D i, Point_hashtable0* poi);
  bool Send_cell_in_direction(Point_hashtable0* poi, 
                              int t, 
                              Index3D my_lev_index, Index3D next_index);


  // send points of type != exterior with I<<my_index
  void Send_grid_points_direct_parallel(Pointtype typ);
  bool Send_point_direct_in_direction(int i, Point_hashtable1* poi,
                                      Index3D next_index);

  // send points of type >=interior || parallel_p with I<<my_index
  // and gives them type: -typ        if   I<<my_index
  //                      -parallel_p else
  void Send_coarse_grid_points_parallel();
  void Send_coarse_grid_points_parallel(int level);
  bool Send_coarse_point_direct_in_direction(int i, Point_hashtable1* poi,
                                             int t, Index3D my_lev_index,
                                             Index3D next_index);

  // for labels
  bool Send_label_in_direction_A(int i, Point_hashtable1* poi,
                               int t, Index3D my_lev_index,
                               Index3D next_index);
  bool Send_label_in_direction_B(int i, Point_hashtable1* poi,
                                 int t, Index3D my_lev_index,
                               Index3D next_index);
  void Send_label_parallel(int level, int num);

  // for multi grid points
  void Send_multi_grid_points_parallel(int level);
  bool Send_multi_grid_point_in_direction(int d, Point_hashtable1* poi,
                                          int t, Index3D my_lev_index, 
                                          Index3D next_index);

  bool Send_prol_point_in_direction(int d, Point_hashtable1* poi, int t,
                                    Index3D my_lev_index, Index3D next_index);

  bool Send_I_point_in_direction(int d, Point_hashtable1* poi, int t,
                                 bool also_f_p, Index3D my_lev_index,
                                 Index3D next_index);
  bool Send_B_point_in_direction(int d, Bo2Point* poi, Index3D next_index);
  bool Send_Z_point_in_direction(int d, BoCell* poi, Index3D next_index);
  bool Send_boundary_stencils_in_direction(Point_hashtable0* poi, 
                                           int t, Index3D my_lev_index,
                                           Index3D next_index);

  bool Send_interior_stencils_in_direction(Point_hashtable0* poi, 
                                           int t, Index3D my_lev_index,
                                           Index3D next_index);

  void Prepare_communication();
  void Prepare_communication_coarser_grids();
  void Prepare_communication_all_grids();
  void Prepare_communication_boundary_stencils();
  void Prepare_communication_boundary_stencils(int level);
  void Prepare_communication_interior_stencils();
  void Prepare_communication_interior_stencils(int level);

  // for face correction in grid generation
  void Start_for_face_correction_parallel();
  void End_for_face_correction_parallel();
  
  int*        edges_to_be_sent[18];   // 12 edges and 6 faces
  int  lenght_edges_to_be_sent[18];
  int  number_edges_to_be_sent[18];
  void Test_array(int i);             // test if above arrays are large enough
  
  int direction_and_number[2];
  int *direction_and_number_process;
  int *receive_buffer;
  int length_receive_buffer;

  int* Give_receive_buffer(int length);

  void put_for_send_edge(Index3D I);
  bool must_edges_be_send();
  bool send_edges();


  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^


  // Information ueber das Gebiet: Eingabe in Rechenkoordinaten
  // Abstand zum Rand in Richting d, Ausgabe aber physikalischer Abstand:
  double distanceD(D3vector V, dir_3D d);
  double offset_square;

  void Initialize_variable();             // Initialisiere Speicher

  // fuer irgendeine  Iteration
  int i_iter;

  // Fuer hashtable0: Zellen, Kanten
  // --------------------------------
  // Funktionen zur Anwendung:
  Point_hashtable0* hashtable0_point(Index3D I) const;
  void Initialize_hash0(int lenght);

  // Speicher des hashtable0:
  Point_hashtable0* point0;             //fuer Iterationen
  int hashtable0_leng;                  // Laenge des hashtable
  int hashtable0_occ;                   // Besetzungszahl des hashtable
  Point_hashtable0** hashtable0_start;  //Zeiger auf Tabelle

  // Fuer hashtable1: Grid
  // ---------------------
  // Funktionen zur Anwendung:
  Point_hashtable1* hashtable1_point(Index3D I) const;
  void Initialize_hash1(int lenght);
  void Resize_hash0(int lenght);
  void Resize_hash1(int lenght);

  // Speicher des hashtable1:
  Point_hashtable1* point1;             //fuer Iterationen
  int hashtable1_leng;                  // Laenge des hashtable
  int hashtable1_occ;                   // Besetzungszahl des hashtable
  Point_hashtable1** hashtable1_start;  //Zeiger auf Tabelle

  // Fuer hashtable2: Randzellen
  // ----------------------------
  // Funktionen zur Anwendung:
  Point_hashtable2* hashtable2_point(Index3D I) const;
  void Initialize_hash2(int lenght);
  int Calc_number_boundary_cells_for_hash2();

  // Speicher des hashtable2:
  Point_hashtable2* bocell;             //fuer Iterationen
  int hashtable2_leng;                  // Laenge des hashtable
  int hashtable2_occ;                   // Besetzungszahl des hashtable
  Point_hashtable2** hashtable2_start;  //Zeiger auf Tabelle

  // Fuer hashtable3: Randpunkte 2.Teil
  // -----------------------------------
  // Funktionen zur Anwendung:
  Point_hashtable3* hashtable3_point(IndexBo IB) const;
  Point_hashtable3* hashtable3_point(Index3D Ind, dir_3D dir) const;
  void Initialize_hash3(int lenght);

  // Speicher des hashtable3:
  Point_hashtable3* bo2point;             //fuer Iterationen
  int hashtable3_leng;                  // Laenge des hashtable
  int hashtable3_occ;                   // Besetzungszahl des hashtable
  Point_hashtable3** hashtable3_start;  //Zeiger auf Tabelle

  // Fuer hashtable4: Datenspeicher
  // -------------------------------
  // Funktionen zur Anwendung:
  Point_hashtable4* hashtable4_point(Index3D, int level);
  Point_hashtable4* hashtable4_point(Index3D); // level=0
  void Initialize_hash4(int lenght);

  // Speicher des hashtable4:
  Point_hashtable4* varpoint;           // fuer Iterationen
  int hashtable4_leng;                  // Laenge des hashtable
  int hashtable4_occ;                   // Besetzungszahl des hashtable
  Point_hashtable4** hashtable4_start;   // Zeiger auf Tabelle
  

  // Sonstige Funktionen, die von obigen Funktionen aufgerufen werden: 
  bool Calc_Is_Slave(Index3D I);
  bool Decide_poi_interpolates(Index3D I);
  void Recursion_Add_27Stencil(Index3D I, int l);
  int  Recursion_Count_Cells(Index3D I);
  int  Recursion_Cells_parallel(Index3D I, int number, 
                                int number_cell_var,
                                double* buffer_send_double);
  int  Recursion_Cells_AVS(Index3D I, ofstream *Datei, int nummer);
  int  Recursion_Cells_AVS(Index3D I, ofstream *Datei, int nummer, 
                           int number_cell_variable);
  int  Recursion_Cells_AVS_parallel(Index3D I, int nummer, int* buffer);
  int      Write_Cells_AVS_parallel(Index3D I, int nummer, int* buffer);
  int  Recursion_Edges_AVS(Index3D I, ofstream *Datei, int nummer);
  int  Recursion_Cells_OPENDX(Index3D I, ofstream *Datei, int nummer);
  void Recursion_Print_cell_info(Index3D I);
  void Recursion_Construct_cell_points_hash0(Index3D I);  
  void Recursion_Calc_type_of_edges_uniform(Index3D I);
  // only for above function
  D3vector A_domain_sp;
  D3vector B_domain_sp;

  Celltype Recursion_calc_interior_cells_Part1(Index3D I);
  Celltype Recursion_calc_interior_cells_Part2(Index3D I);
  void     Recursion_calc_interior_cells_Part3(Index3D I, int level);
  Celltype Calc_cell_type(Index3D I);
  void Recursion_add_Uniform_Grid(Index3D I,int Max_tiefe);

  // fuer AVS:
  void Print_Cell_avs(Index3D I, ofstream *Datei);
  void Print_ec_avs(IndexBo indexbo_A, IndexBo indexbo_B, 
                              ofstream *Datei, int number);
  int avs_bo_cell(BoCell* bo,ofstream *Datei,bool print_or_calc,int number);
  int avs_bo_cell(BoCell* bo,ofstream *Datei,int number,
                  int number_cell_var);
  int avs_bo_cell_parallel_surface(BoCell* bo,int number, 
                                   int* buffer, bool print_or_calc);
  int avs_bo_cell_parallel(BoCell* bo,int number, int* buffer);
  int avs_bo_cell_parallel(BoCell* bo,int number, int number_cell_var,
                           double* buffer_send_double);

  // for test AVS
  int avs_bo_cell_typ(BoCell* bo,ofstream *datei, int number,double value);
  int Recursion_Cell_typ_AVS(Index3D I, ofstream *Datei, int nummer);
  int  Recursion_Count_all_Cells(Index3D I);

  // not for everybody:
  bool is_negative(BoCell* bo, Tetraeder_storage* tets);
  bool is_negative_surf(BoCell* bo, Tetraeder_storage* tets);
  int avs_bo_cell_half(BoCell* bo,ofstream *Datei,bool print_or_calc,int number);
  int count_edgepoint(Tetraeder_storage* tet, BoCell* bo);
  int count_nearb_surf_half(BoCell* bo, bool posi);
  int count_nearb_surf_transv(BoCell* bo, bool posi);
  int avs_bo_cell_surf(BoCell* bo,ofstream *datei,bool print_or_calc, 
                       int number); 
  int avs_bo_cell_surf_half(BoCell* bo,ofstream *datei,bool print_or_calc, 
                       int number); 
  int avs_bo_cell_surf_transv(BoCell* bo,ofstream *datei,bool print_or_calc, 
                       int number); 

  int  Recursion_Count_Cells_half(Index3D I);
  int  Recursion_Count_Cells_surf_half(Index3D I);
  int  Recursion_Count_Cells_surf_transv(Index3D I);

  int  Recursion_Cells_AVS_half(Index3D I, ofstream *Datei, int nummer);
  int  Recursion_Cells_AVS_surf_half(Index3D I, ofstream *Datei, int nummer);
  int  Recursion_Cells_AVS_surf_transv(Index3D I, ofstream *Datei, int nummer);


  // fuer OpenDx:
  void Print_Cell_opendx(Index3D I, ofstream *Datei);
  void Print_ec_opendx(IndexBo indexbo_A, IndexBo indexbo_B, 
                              ofstream *Datei, int number);
  int opendx_bo_cell(BoCell* bo,ofstream *Datei,bool print_or_calc,int number);
};


// Hashtable einbinden
#include "hash.h"

// inline Funktionen fuer schnellen Zugriff

inline bool Grid_base::Point_in_domain(Index3D I) const {
  Point_hashtable1* point;
  point =  hashtable1_point(I);
  if(point==NULL) return false;
  // old
  //  if(point->typ < interior) return false;
  // new 
  if(point->typ < interior && point->typ!=parallel_p) return false;
  return true;
}

inline bool Grid_base:: Exists_Cell(Index3D I) {
  if(developer_version) {
    if(!I.Cell_index()) cout << "\n Error in Exists_Cell!" << endl;
  }
  return hashtable0_point(I)!=NULL;
}

inline bool Grid_base::Exists_Point(Index3D I) const {
  return hashtable1_point(I)!=NULL;
}

inline bool Grid_base::Exists_Bo_Point(IndexBo I) const {
  return hashtable3_point(I)!=NULL;
}

inline bool Grid_base::Exists_Bo_Point(Index3D I, dir_3D d) const {
  return hashtable3_point(I,d)!=NULL;
}


inline int  Grid_base::Max_level() const { return max_level; };
inline int  Grid_base::Min_level() const { return min_level; };

inline double* Grid_base::Give_cell_variable(Index3D I) const {
  return hashtable0_point(I)->var;
}

inline double* Grid_base::Give_pointer_stencil(Index3D I) const {
  return hashtable0_point(I)->var;
}

inline double* Grid_base::Give_bo_stencil(Index3D I,int num_stencil) const {
  return &(hashtable0_point(I)->var[cell_number_of_bo_stencil(num_stencil)]);
}

inline double* Grid_base::Give_stencil(Index3D I,int num_stencil) const {
  return &(hashtable0_point(I)->var[cell_number_of_stencil(num_stencil)]);
}

inline double* Grid_base::Give_variable(Index3D Ind, dir_3D dir) const {
  static Point_hashtable3* point;
  point = hashtable3_start
    [hashtable3_function(Ind.ind_x.index,Ind.ind_y.index,
                         Ind.ind_z.index, dir, hashtable3_leng)];
  while(point!=NULL) {
    if(point->ind==Ind && point->direction==dir) return point->var;
    point = point->next;
  }
  // old
  //  return point->var;
  // new
  return NULL;
}


inline double* Grid_base::Give_variable_cellpoi(Index3D Ind) const {
  static Point_hashtable2* point;
  point = hashtable2_start
    [hashtable2_function(Ind.ind_x.index,Ind.ind_y.index,
                         Ind.ind_z.index, hashtable2_leng)];
  while(point!=NULL) {
    if(point->ind==Ind) return point->var;
    point = point->next;
  }
  if(developer_version) {
    if(point==NULL) {
      cout << " Error 1 in Give_variable_cellpoi! my rank"
           << my_rank 
           << " x: " << Ind.I_x().get()
           << " y: " << Ind.I_y().get()
           << " z: " << Ind.I_z().get()
           << endl;
      return NULL;
    }
    else 
      if(point->var==NULL)   
        cout << " Error 2 in Give_variable_cellpoi! " << endl;
  }
  return point->var;
}


inline double Grid_base::Give_h(Index3D Ind, dir_3D dir) const {
  static Point_hashtable3* point;
  point = hashtable3_start
    [hashtable3_function(Ind.ind_x.index,Ind.ind_y.index,
                         Ind.ind_z.index, dir, hashtable3_leng)];
  while(point!=NULL) {
    if(point->ind==Ind && point->direction==dir) return point->h;
    point = point->next;
  }
  return point->h;
}

// Maschenweite
inline double Grid_base::H_mesh() const {
  return H_bounding;  
};

inline D3vector Grid_base::Give_A() const {
  return A_bounding;  
}

// transform computational coordinate in physical coordinate
inline D3vector Grid_base::transform_coord(const D3vector V) const {
  return D3vector(Give_A().x + V.x*H_mesh(),
                  Give_A().y + V.y*H_mesh(),
                  Give_A().z + V.z*H_mesh());
};

// transform physical coordinate in computational coordinate
inline D3vector Grid_base::back_transform_coord(const D3vector V) const {
  return D3vector((-Give_A().x + V.x)/H_mesh(),
                  (-Give_A().y + V.y)/H_mesh(),
                  (-Give_A().z + V.z)/H_mesh());
};

// transform Index3D in physical coordinate
inline D3vector Grid_base::transform_coord(const Index3D I) const {
  return transform_coord(I.coordinate());
};
inline double Grid_base::transform_coordX(const Index3D I) const {
  return Give_A().x + H_mesh() * I.coordinate_x();
};
inline double Grid_base::transform_coordY(const Index3D I) const {
  return Give_A().y + H_mesh() * I.coordinate_y();
};
inline double Grid_base::transform_coordZ(const Index3D I) const {
  return Give_A().z + H_mesh() * I.coordinate_z();
};

// transform computational coordinate in physical coordinate
inline D3vector Grid_base::transform_coord(const Index3D I, const dir_3D d) 
     const {
  return hashtable3_point(I,d)->transform_coord(Give_A(),H_mesh());
};

inline D3vector Grid_base::transform_coord_org_shift(const Index3D I,
                                                     const dir_3D d,
                                                     const double shift) 
     const {
  return hashtable3_point(I,d)->transform_coord_org_shift(Give_A(),
                          H_mesh()/Zweipotenz(Max_level()),
                          H_mesh(),shift);
};


// transform Index3D of cellpoint in physical coordinate
inline D3vector Grid_base::local_coord_cellpoi(const Index3D Ind) const {
  static Point_hashtable2* point;
  point = hashtable2_start
    [hashtable2_function(Ind.ind_x.index,Ind.ind_y.index,
                         Ind.ind_z.index, hashtable2_leng)];
  while(point!=NULL) {
    if(point->ind==Ind) return point->Local_coord_bocellpoint();
    point = point->next;
  }
  return point->Local_coord_bocellpoint();
};

inline D3vector Grid_base::transform_coord_cellpoi(const Index3D I) const {
  return transform_coord(I.coordinate()) + local_coord_cellpoi(I);
};

inline double Grid_base::transform_coordX_cellpoi(const Index3D I) const {
  return Give_A().x + H_mesh() * I.coordinate_x()
    + local_coord_cellpoi(I).x;
};

inline double Grid_base::transform_coordY_cellpoi(const Index3D I) const {
  return Give_A().y + H_mesh() * I.coordinate_y()
    + local_coord_cellpoi(I).y;
};

inline double Grid_base::transform_coordZ_cellpoi(const Index3D I) const {
  return Give_A().z + H_mesh() * I.coordinate_z()
    + local_coord_cellpoi(I).z;
};


inline Celltype Grid_base::Give_cell_typ(Index3D I) const {
  static Point_hashtable0* point;
  if(developer_version) {
    if(!I.Cell_index()) {
      cout << "\n Error 1 in Give_cell_typ!" << endl;
      cout << " Point: "; I.coordinate().Print(); cout << endl;
    }
  }
  point = hashtable0_point(I);
  if(point==NULL) {
    return no_cell_exists;
  }
  return (Celltype)point->typ;
}


#endif

      

---