src/expde/grid/gridbase.cc

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
//

// ------------------------------------------------------------
//
// grid_base.cc
//
// ------------------------------------------------------------

//  A. Grid_base::Initialize_variable

// Include level 0:
#ifdef COMP_GNUOLD
 #include <iostream.h>
 #include <fstream.h>
 #include <time.h>
 #include <math.h>
#else
 #include <iostream>
 #include <fstream>
 #include <ctime>
 #include <cmath>
#endif
 
#ifdef USE_PBE
#include "pbedefs.h"
#endif

// Include level 0:
#include "../parser.h"

// Include level 1:
#include "../paramete.h"
#include "../abbrevi.h"
#include "../math_lib/math_lib.h"

// Include level 2:
#include "../basic/basic.h"

// Include level 3:
#include "../domain/domain.h"

// Include level 4:
#include "../formulas/boundy.h"
#include "../formulas/loc_sten.h"

// Include level 5:
#include "gpar.h"
#include "parallel.h"
#include "mgcoeff.h"
#include "sto_man.h"
#include "gridbase.h"
#include "grid.h"


//  Z. Grid_base::Grid_base

Grid_base::Grid_base(int n_max, All_Domains* dom,
                     Grid_gen_parameters& gpara, MPI_Comm comm) : 
  Parallel_Info(dom,comm) {

  if(dom->Is_periodic()) {
    cout << " please do not use this constructor! " << endl;
  }

  // for parallel:
  refine_level = gpara.Give_r_parallel();

  // calc offset
  offset_square = gpara.Give_offset_square();
  h_offset = offset_square * dom->GiveH() / Zweipotenz(n_max);

  // set bounding_box parameters
  A_bounding = dom->GiveA() + D3vector(-h_offset, -h_offset, -h_offset);
  H_bounding = dom->GiveH()+2.0*h_offset;
  H_bounding = gpara.Give_stretch_square() * H_bounding;

  max_level=n_max;

  // move A_bounding if necessary
  if(gpara.Is_there_grid_point()) {
    make_grid_with_grid_point(gpara.Give_grid_point(),
                              H_bounding / Zweipotenz(max_level));
  }

  // -> parallel
  if(parallel_version) {
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << " Generate_processes(); " << endl;
    }
    Generate_processes();
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << " Grid_generation(); " << endl;
    }
    if(I_am_active()) Grid_generation(n_max);
    else Dummy_grid_generation();
  }
  // -> seriell
  else {
    Generate_seriel_process();
    Grid_generation(n_max);
  }
}

Grid_base::Grid_base(double h_finest_level, All_Domains* dom,
                     Grid_gen_parameters& gpara, MPI_Comm comm) : 
  Parallel_Info(dom,comm) {
  int i;

  if(dom->Is_periodic()) {
    cout << " please do not use this constructor! " << endl;
  }

  // for parallel:
  refine_level = gpara.Give_r_parallel();

  // calc offset
  offset_square = gpara.Give_offset_square();
  h_offset = offset_square * h_finest_level;

  // set bounding_box parameters
  A_bounding = dom->GiveA() + D3vector(-h_offset, -h_offset, -h_offset);
  H_bounding = dom->GiveH()+2.0*h_offset;
  H_bounding = gpara.Give_stretch_square() * H_bounding;


  if(h_finest_level > H_bounding / 4.0) 
    cout << " Coose large meshsize h_finest_level!! " << endl;

  // move A_bounding if necessary
  if(gpara.Is_there_grid_point()) {
    make_grid_with_grid_point(gpara.Give_grid_point(),h_finest_level);
  }

  // h_finest_level = H_bounding / Zweipotenz(max_level)
  for(i=2;(i<Max_number_levels) && 
        (H_bounding > ((double)(Zweipotenz(i)) * h_finest_level));i++) { };
  max_level=i;
  H_bounding = Zweipotenz(max_level) * h_finest_level;


  // -> parallel
  if(parallel_version) {
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << " Generate_processes(); " << endl;
    }
    Generate_processes();
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << " Grid_generation(); " << endl;
    }
    if(I_am_active()) Grid_generation(max_level);
    else Dummy_grid_generation();
  }
  // -> seriell
  else {
    Generate_seriel_process();
    Grid_generation(max_level);
  }
}

Grid_base::Grid_base(int n_max, All_Domains* dom, MPI_Comm comm) : 
  Parallel_Info(dom,comm) {
  Grid_gen_parameters gpara;

  if(dom->Is_periodic()) {
    gpara.Set_offset_square(0.0);
    gpara.Set_stretch_square(1.0);
  }

  // for parallel:
  refine_level = gpara.Give_r_parallel();

  // calc offset
  offset_square = gpara.Give_offset_square();
  h_offset = offset_square * dom->GiveH() / Zweipotenz(n_max);


  // set bounding_box parameters
  A_bounding = dom->GiveA() + D3vector(-h_offset, -h_offset, -h_offset);
  H_bounding = dom->GiveH()+2.0*h_offset;
  H_bounding = gpara.Give_stretch_square() * H_bounding;

  max_level=n_max;

  // -> parallel
  if(parallel_version) {
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << " Generate_processes(); " << endl;
    }
    Generate_processes();
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << " Grid_generation(); " << endl;
    }
    if(I_am_active()) Grid_generation(max_level);
    else Dummy_grid_generation();
  }
  // -> seriell
  else {
    Generate_seriel_process();
    Grid_generation(max_level);
  }
}

Grid_base::Grid_base(double h_finest_level, All_Domains* dom, MPI_Comm comm) : 
  Parallel_Info(dom,comm) {
  int i;
  Grid_gen_parameters gpara;

  if(dom->Is_periodic()) {
    cout << " please do not use this constructor! " << endl;
  }

  // for parallel:
  refine_level = gpara.Give_r_parallel();

  // calc offset
  offset_square = gpara.Give_offset_square();
  h_offset = offset_square * h_finest_level;

  // set bounding_box parameters
  A_bounding = dom->GiveA() + D3vector(-h_offset, -h_offset, -h_offset);
  H_bounding = dom->GiveH()+2.0*h_offset;
  H_bounding = gpara.Give_stretch_square() * H_bounding;

  if(h_finest_level > H_bounding / 4.0) 
    cout << " Coose large meshsize h_finest_level!! " << endl;

  // es muss gelten: 
  // h_finest_level = H_bounding / Zweipotenz(max_level)
  for(i=2;(i<Max_number_levels) && 
        (H_bounding > ((double)(Zweipotenz(i)) * h_finest_level));i++) { };
  
  max_level=i;
  H_bounding = Zweipotenz(max_level) * h_finest_level;

  // -> parallel
  if(parallel_version) {
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << " Generate_processes(); " << endl;
    }
    Generate_processes();
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << " Grid_generation(); " << endl;
    }
    if(I_am_active()) Grid_generation(max_level);
    else Dummy_grid_generation();
  }
  // -> seriell
  else {
    Generate_seriel_process();
    Grid_generation(max_level);
  }
}


double Grid_base::finest_mesh_size() const {
  return H_mesh() / Zweipotenz(Max_level());
}

//  A. Grid_base::Initialize_variable


void Grid_base::Initialize_variable() {
 if(print_status_of_grid_generator) if(my_rank==0) {
    cout << " Initialize_variable(); " << endl;
  }
  //  if(Is_storage_initialized() == false) {
  // Continue: Teil III:
  // III.3.1.Schritt: Randpunkte 2.Teil berechnen
  //                und Berechnung der FE-Zerlegung
  //                und Datenspeicher bereitstellen fuer Randpunkte
  //                and one list of P_boundary_tet
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   III.3.1 construct boundary points; " << endl;
  }
#ifdef USE_PBE
    pbe_start(19,"III-3.1_construct_boundary_points");
#endif
  auxiliary_P = NULL;
  auxiliary_number = 0;
  h_min_for_boundary_points = 
    finest_mesh_size() * finest_mesh_size() 
    / domain->GiveH() * factor_boundary_dis;
  Calc_boundary_2points(Give_max_num_var());
#ifdef USE_PBE
  pbe_stop(19);
#endif

  // III.3.2. Schritt: beseitige Kanten
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "         remove edges; " << endl;
  } 
#ifdef USE_PBE
    pbe_start(20,"III-3.2_remove_edges");
#endif
  Remove_edges();

  // III.3.3. Schritt: 
  Calc_multigrid_points_part2();
#ifdef USE_PBE
  pbe_stop(20);
#endif

  // III.4.1. Schritt: Datenspeicher bereitstellen fuer Punkte in Randzellen
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   III.4.1 storage for boundary cells; " << endl;
  } 
#ifdef USE_PBE
    pbe_start(21,"III-4.1_storage_for_bo_cells");
#endif
  iterate_hash2 {
    if(bocell->Exists_bocellpoint()) {
      if(Give_max_num_var()>0)  bocell->var = newdouble(Give_max_num_var());
      else bocell->var   = NULL;
    }
  }
#ifdef USE_PBE
  pbe_stop(21);
#endif
  // III.4.2. Schritt: Datenspeicher bereitstellen Zellen-Variablen 
  //                   und Stencils
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   III.4.2 storage for cells and stencils; " << endl;
  } 
#ifdef USE_PBE
    pbe_start(22,"III-4.2_storage_for_cells_and_stencils");
#endif
  iterate_hash0 {
    if(point0->Give_Index().Cell_index()) {
      if(point0->typ==bo_cell) {
        if(Give_max_num_bo_cells()>0) 
          point0->var = newdouble(Give_max_num_bo_cells());
        else point0->var = NULL;
      }
      if(point0->typ==int_cell) {
        if(point0->Give_Index().Tiefe()<=Max_level()) {
          if(Give_max_num_coarse_int_cells()>0) 
            point0->var = newdouble(Give_max_num_coarse_int_cells());
          else point0->var = NULL;
        }
      }
      // new:
      if(point0->typ==int_cell) {
        if(point0->Give_Index().Tiefe()==Max_level()+1) {
          if(Give_max_num_fine_all_cells()>0)
            point0->var = newdouble(Give_max_num_fine_all_cells());
          else point0->var = NULL;
        }
      }
    }
  }
#ifdef USE_PBE
  pbe_stop(22);
#endif
  // III.4.3. Schritt: Datenspeicher bereitstellen fuer Matrizen in Randzellen 
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   III.4.3 storage for boundary matrices; " << endl;
  }


  // new:
  // -> boundary tets: see Grid_base::Calc_boundary_2points
  /*
  Tetraeder_storage* tet;
  iterate_hash2 {
    for(tet=bocell->Give_tets();tet!=NULL;tet=tet->Next()) {
      tet->Set_var(newdouble(Give_max_num_fine_all_cells()));
    }
  }
  */
  //  }
  //  Storge_initialization_done();
}

// Schritt III.3.1
//-----------------

void Grid_base::Calc_boundary_2points(int number_var) {
  static int i, d;
  static IndexEdge edge;
  static Index3D I, Icorner;
  Bo_description  Bodes;
  Pointtype typpoi;
  static dir_3D oppdir;
  static double h;

  // for list P_boundary_tet:
  Tetraeder_storage* tet;
  double h_[24];
  D3vector V;
  //

  iterate_hash2 {
    I = bocell->Give_Index();
    Bodes.Put_zero();

    for(i=0;i<12;++i) {
      edge = I.edge((Edges_cell)i);
      if(Give_edge_typ(edge.I()) == edge_interior) {
        // Kante im Gebiet
        Bodes.Put_edge((Edges_cell)i);  
      }
      else {
        // Kante nicht im Gebiet
        for(d=0;d<2;++d) {
          Icorner = edge.corner((dir1D)d);       
          typpoi=Give_type(Icorner);
          // typpoi=exterior;
          if(typpoi!=exterior && typpoi!=multigrid &&  typpoi!=boundary) {
            oppdir = edge.opposite_direction((dir1D)d);
            Add_bo2point(Icorner,oppdir,I.Tiefe()-1,number_var);
            // Kantenpunkt im Gebiet
            h = hashtable3_point(Icorner,oppdir)->h;
            Bodes.Put_edge_point(Transform((Edges_cell)i,(dir1D)d),oppdir,h); 
          }
        }
      }
    }

    // Welche Eckpunkte sind im Gebiet?
    for(i=0;i<8;++i) {
      typpoi=Give_type(I.neighbour((dir_sons)i));
      if(typpoi!=exterior && typpoi!=multigrid) 
        Bodes.Put_corner_point((dir_sons)i); 
    }

    Bodes.Put_Meshsize(H_mesh()/Zweipotenz(I.Tiefe()-1));
    bocell->Calc_FE_elements(Bodes);    

    // for list P_boundary_tet:
    if(Give_max_num_fine_all_cells()>0) {
      bocell->Set_h(h_,Bodes);
      
      for(tet=bocell->Give_tets();tet!=NULL;tet=tet->Next()) {
        tet->Set_var(newdouble(Give_max_num_fine_all_cells()));

        V = (bocell->coord(tet->N0(),h_) +
             bocell->coord(tet->N1(),h_) +
             bocell->coord(tet->N2(),h_) +
             bocell->coord(tet->N3(),h_)   ) / 4.0 +
          transform_coord(I.neighbour(WSDd));
        
        auxiliary_P = 
          new P_boundary_tet(V,auxiliary_P,
                             tet->Give_variable(),tet,bocell);
        ++auxiliary_number;
      }
    }
  }
}



void Grid_base::Information() {
  cout << "Information about the grid: " << endl;
  Info_hashtable0();
  Info_hashtable1();
  Info_hashtable2();
  Info_hashtable3();
  Info_hashtable4();
}


int Grid_base::Give_Nummer(Index3D I) {
  static Point_hashtable1* point;

  if(developer_version) {
    point=hashtable1_point(I);
    if(point==NULL) {
      cout << " \n Punkt existiert nicht in Give_Nummer!" << endl;
      return 0;
    }
    else return point->nummer;
  }
  else return hashtable1_point(I)->nummer;
}

int Grid_base::Give_Nummer(Index3D I, dir_3D d) {
  /* old Vis :
  if(developer_version) {
    if(hashtable3_point(I,d)==NULL) {
      cout << " \n Grid_base::Fehler Give_Nummer!" << endl;
      return 0;
    }
    return hashtable3_point(I,d)->Nummer()+hashtable1_occ;
  }
  else return hashtable3_point(I,d)->Nummer()+hashtable1_occ;
  */
  if(developer_version) {
    if(hashtable3_point(I,d)==NULL) {
      cout << " \n Grid_base::Fehler Give_Nummer!" << endl;
      return 0;
    }
    return hashtable3_point(I,d)->Nummer();
  }
  else return hashtable3_point(I,d)->Nummer();
}

int Grid_base::Give_Nummer_cellpoi(Index3D I) {
  if(developer_version) {
    if(hashtable2_point(I)==NULL) {
      cout << " \n Fehler Give_Nummer_cellpoi!" << endl;
      return 0;
    }
    return hashtable2_point(I)->Number_avs();
  }
  else return hashtable2_point(I)->Number_avs();
}





Pointtype Grid_base::Give_type(Index3D I) const {
  Point_hashtable1* point;

  point=hashtable1_point(I);
  if(point==NULL) {
    return exterior;
  }
  else return point->typ;
}

Pointtype Grid_base::Give_global_type(Index3D I) const {
  Point_hashtable1* point;

  point=hashtable1_point(I);
  if(point==NULL) {
    return exterior;
  }
  else return point->global_typ;
  // else return point->typ;
}

bool Grid_base::Grid_point_on_finest_level(Index3D I) const {
  Point_hashtable1* point;

  point=hashtable1_point(I);
  if(point==NULL) {
    return false;
  }
  if(point->typ  >= interior) return true;
  if(point->typ  == parallel_p && // point->global_typ  != multigrid && 
     point->my_finest_parallel_level == max_level)
    return true;
  return false;
}


void Grid_base::Add_double(Index3D I, int varebene, int number_var) {
  // static Point_hashtable4* point;
  Put_finest_level_minimal(I,varebene);

  Point_hashtable4* point;

  point = hashtable4_start[hashtable4_function(
       I.ind_x.index,I.ind_y.index,I.ind_z.index,varebene,hashtable4_leng)];
  
  if(point==NULL) {
    hashtable4_occ++;
    point = new Point_hashtable4(I,varebene,number_var);
    hashtable4_start[hashtable4_function(
       I.ind_x.index,I.ind_y.index,I.ind_z.index,varebene,hashtable4_leng)] = 
       point;
  }
  else {
    while(point->next!=NULL && (point->ind!=I || point->level!=varebene))
      point = point->next;
    if(point->ind!=I || point->level!=varebene) {
      hashtable4_occ++;
      point->next = new Point_hashtable4(I,varebene,number_var);
    }
  }


  /* old_hash
  static Point_hashtable4* point;
  point = &(hashtable4_start
    [hashtable4_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,varebene,
                         hashtable4_leng)]);
  if(point->level==-2) {
    hashtable4_occ++;
    point->ind = I;
    point->level = varebene;
    point->var = newdouble(number_var);
    if(developer_version) {
      if(point->next!=NULL) cout << "\n Fehler 1 in Add_double!" << endl;
    }
  }
  else {
    while(point->next!=NULL && (point->ind!=I || point->level!=varebene)) 
      point = point->next;
    if(point->ind!=I || point->level!=varebene) {
      hashtable4_occ++;
      point->next = new Point_hashtable4(I,varebene,number_var);
      if(developer_version) {
        if(point->next->next!=NULL) cout << "\n Fehler 2 in Add_double!" << endl;
      }
    }
  }
  */
}


void Grid_base::Add_MG_double(int varebene, int number_var) {
  iterate_hash1 {
    if(point1->label4)
      Add_double(point1->ind,varebene,number_var);
  }
}

void Grid_base::False_Label4() {
  iterate_hash1 {
    point1->label4 = false;
  }
}


void Grid_base::Add_cell(Index3D I) {
  static Point_hashtable0* point;
  if(developer_version) {
    if(!I.Cell_index()) cout << "\n Error in Add_cell!" << endl;
  }
  point = hashtable0_start
    [hashtable0_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,
                         hashtable0_leng)];
  if(point==NULL) {
    hashtable0_occ++;
    point = new Point_hashtable0(I);
    hashtable0_start[hashtable0_function(I.ind_x.index,I.ind_y.index,
                                         I.ind_z.index,
                                         hashtable0_leng)] = point;
  }
  else {
    while(point->next!=NULL && point->ind!=I) point = point->next;
    if(point->ind!=I) {
      hashtable0_occ++;
      point->next = new Point_hashtable0(I);
    }
  }
}

void Grid_base::Add_edge(Index3D I) {
  // cout << "\n Kante: "; I.coordinate().Print();

  static Point_hashtable0* point;
  if(developer_version) {
    if(!I.Edge_index()) {
      cout << "\n Error in Add_edge!" << endl;
      return;
    }
  }
  point = hashtable0_start
    [hashtable0_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,
                         hashtable0_leng)];
  if(point==NULL) {
    hashtable0_occ++;
    point = new Point_hashtable0(I);
    point->typ = edge_interior;
    hashtable0_start[hashtable0_function(I.ind_x.index,I.ind_y.index,
                                         I.ind_z.index,
                                         hashtable0_leng)] = point;
  }
  else {
    while(point->next!=NULL && point->ind!=I) point = point->next;
    if(point->ind!=I) {
      hashtable0_occ++;
      point->next = new Point_hashtable0(I);
      point->next->typ = edge_interior;
    }
  }
}

void Grid_base::Add_edge_parallel(Index3D I) {
  if((my_index<=I)==false) { 
    //    cout << " \n face correction not perfect implemented! " << endl;
    put_for_send_edge(I);
  }
  Add_edge(I);
}

Edgetype Grid_base::Give_edge_typ(Index3D I) {
  static Point_hashtable0* point;
  if(developer_version) {
    if(!I.Edge_index()) {
      cout << "\n Error 1 in Give_edge_typ!" << endl;
      cout << " Point: "; I.coordinate().Print(); cout << endl;
    }
  }
  point = hashtable0_point(I);
  if(point==NULL)
    return edge_not_interior;
  else  
    return (Edgetype)point->typ;
}

void Grid_base::Set_edge_typ(Index3D I, Edgetype type) {
  static Point_hashtable0* point;
  if(developer_version) {
    if(!I.Edge_index()) cout << "\n Error 1 in Set_edge_typ!" << endl;
  }
  if(developer_version) {
    if(type <0 || type>2) cout << "\n Error 3 in Set_edge_typ!" << endl;
  }
  if(type==edge_interior) {
    point = hashtable0_point(I);
    if(developer_version) {
      if(point==NULL) cout << "\n Error 2 in Set_edge_typ!" << endl;
    }
    point->typ = (Edgetype)type;
  }
}

void Grid_base::Set_cell_typ(Index3D I, Celltype type) {
  static Point_hashtable0* point;
  if(developer_version) {
    if(!I.Cell_index()) cout << "\n Error 1 in Set_cell_typ!" << endl;
  }
  if(developer_version) {
    if(type <0 || type>4) cout << "\n Error 3 in Set_cell_typ!" << endl;
  }
  point = hashtable0_point(I);
  if(developer_version) {
    if(point==NULL) cout << "\n Error 2 in Set_cell_typ!" << endl;
  }
  point->typ = (Celltype)type;
}


void Grid_base::Set_point_typ(Index3D I, Pointtype type) {
  static Point_hashtable1* point;
  if(developer_version) {
    if(type <0 || type>8) cout << "\n Error 1 in Set_point_typ!" << endl;
  }
  point = hashtable1_point(I);
  if(developer_version) {
    if(point==NULL) cout << "\n Error 2 in Set_point_typ!" << endl;
  }
  point->typ = type;
}

void Grid_base::Add_point(Index3D I, Pointtype type) {
  static Point_hashtable1* point;
  if(developer_version) {
    if(I<<Index3D(2,2,2) == false) {
      cout << " Error in Grid_base::Add_point(Index3D I, Pointtype type)"
           << endl;
    }
  }
  point = hashtable1_start
    [hashtable1_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,
                         hashtable1_leng)];
  if(point==NULL) {
    hashtable1_occ++;

    point = new Point_hashtable1(I);
    point->typ = type;
    hashtable1_start[hashtable1_function(I.ind_x.index,I.ind_y.index,
                                         I.ind_z.index,
                                         hashtable1_leng)] = point;
  }
  else {
    while(point->next!=NULL && point->ind!=I) {
      point = point->next;
    }
    if(point->ind!=I) {
      hashtable1_occ++;

      point->next = new Point_hashtable1(I);
      point->next->typ = type;
    }
  }
}



void Grid_base::Add_point(Index3D I) {
  static Point_hashtable1* point;
  point = hashtable1_start
    [hashtable1_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,
                         hashtable1_leng)];
  if(point==NULL) {
    hashtable1_occ++;
    point = new Point_hashtable1(I);
    hashtable1_start[hashtable1_function(I.ind_x.index,I.ind_y.index,
                                         I.ind_z.index,
                                         hashtable1_leng)] = point;
  }
  else {
    while(point->next!=NULL && point->ind!=I) point = point->next;
    if(point->ind!=I) {
      hashtable1_occ++;
      point->next = new Point_hashtable1(I);
    }
  }
}

void Grid_base::Set_point_typ(Index3D I, int finest_parallel_level,
                          Pointtype type) {
  static Point_hashtable1* point;
  if(developer_version) {
    if(type <0 || type>8) cout << "\n Error 3 in Set_point_typ!" << endl;
  }
  point = hashtable1_point(I);
  if(developer_version) {
    if(point==NULL) cout << "\n Error 4 in Set_point_typ!" << endl;
  }
  point->typ = type;
  point->my_finest_parallel_level = finest_parallel_level;
}

void Grid_base::Add_point_interior(Index3D I) {
  static Point_hashtable1* point;
  point = hashtable1_start
    [hashtable1_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,
                         hashtable1_leng)];
  if(point==NULL) {
    hashtable1_occ++;
    point = new Point_hashtable1(I,interior);
    hashtable1_start[hashtable1_function(I.ind_x.index,I.ind_y.index,
                                         I.ind_z.index,
                                         hashtable1_leng)] = point;
  }
  else {
    while(point->next!=NULL && point->ind!=I) point = point->next;
    if(point->ind!=I) {
      hashtable1_occ++;
      point->next = new Point_hashtable1(I,interior);
    }
  }
}

void Grid_base::Add_bocell(Index3D I) {
  static Point_hashtable2* bocell_sp;
  bocell_sp = hashtable2_start
    [hashtable2_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,
                         hashtable2_leng)];
  if(bocell_sp==NULL) {
    hashtable2_occ++;
    bocell_sp = new Point_hashtable2(I,this);
    hashtable2_start[hashtable2_function(I.ind_x.index,I.ind_y.index,
                                         I.ind_z.index,
                                         hashtable2_leng)] = bocell_sp;
  }
  else {
    while(bocell_sp->next!=NULL && bocell_sp->ind!=I)
      bocell_sp = bocell_sp->next;
    if(bocell_sp->ind!=I) {
      hashtable2_occ++;
      bocell_sp->next = new Point_hashtable2(I,this);
    }
  }
}


void Grid_base::Add_bo2point(Index3D I, dir_3D d, int l, int number_var) {
  static Point_hashtable3* bo2point_sp;
  bo2point_sp = hashtable3_start
    [hashtable3_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,d,
                         hashtable3_leng)];
  if(bo2point_sp==NULL) {
    hashtable3_occ++;
    bo2point_sp = new Point_hashtable3(I,d,l,this,number_var);
    hashtable3_start[hashtable3_function(I.ind_x.index,I.ind_y.index,
                                         I.ind_z.index,d,
                                         hashtable3_leng)] = bo2point_sp;
  }
  else {
    while(bo2point_sp->next!=NULL && (bo2point_sp->ind!=I 
          || bo2point_sp->direction!=d))
      bo2point_sp = bo2point_sp->next;
    if(bo2point_sp->ind!=I || bo2point_sp->direction!=d) {
      hashtable3_occ++;
      bo2point_sp->next = new Point_hashtable3(I,d,l,this,number_var);
    }
  }
}



// Informationen ueber die Punkte

void Grid_base::Label_true(Index3D I) {
  hashtable1_point(I)->label4=true;
}

bool Grid_base::Label_ask(Index3D I) {
  return hashtable1_point(I)->label4;;
}


// Informationen ueber die Randzellen
bool Grid_base::Is_Bo_cell(Index3D I) {
  return hashtable2_point(I)!=NULL;
}
BoCell *Grid_base::Give_Bo_cell(Index3D I) const {
  return hashtable2_point(I);
}


void Grid_base::Fullfill_27Stencil() {
  static int occ_old;
  // Schneller machen !!
  occ_old = 0;
  while(occ_old != hashtable1_occ) {
    occ_old = hashtable1_occ;
    iterate_hash1 {
      Add_27Stencil(point1->ind);
    }
  }
}

void Grid_base::Add_small_27Stencil(Index3D I, int level) {
  static int i, j, k;
  static Index3D Inext;

  for(i=0;i<3;++i)
    for(j=0;j<3;++j)
      for(k=0;k<3;++k) {
        Inext = I.next((Ort1D)i,(Ort1D)j,(Ort1D)k,level);
        if(!Exists_Point(Inext)) Add_point(Inext);
      }
}


void Grid_base::Add_27Stencil(Index3D I) {
  if(I.Tiefe()>0) Recursion_Add_27Stencil(I,I.Tiefe());
}

void Grid_base::Recursion_Add_27Stencil(Index3D I,int l) {
  if(!Exists_Point(I)) Add_point(I);
  if(I.ind_x.Tiefe() == l) {
    Recursion_Add_27Stencil(I.neighbour_EW(Ld),l);
    Recursion_Add_27Stencil(I.neighbour_EW(Rd),l);
  }
  if(I.ind_y.Tiefe() == l) {
    Recursion_Add_27Stencil(I.neighbour_NS(Ld),l);
    Recursion_Add_27Stencil(I.neighbour_NS(Rd),l);
  }
  if(I.ind_z.Tiefe() == l) {
    Recursion_Add_27Stencil(I.neighbour_TD(Ld),l);
    Recursion_Add_27Stencil(I.neighbour_TD(Rd),l);
  }
}



void Grid_base::Add_points_of_cell(Index3D I) {
  static int i;
  static Index3D In;

  for(i=0;i<8;++i) {
    In = I.neighbour((dir_sons)i);
    if(!Exists_Point(In)) Add_point(In);
  }
}



// hier fehlt was: Es reicht nicht wenn der 27-Stern im Gebiet ist;
// Man muss auch testen ob die Zellen im Gebiet sind.
/*
void Grid_base::Calc_interior_points_old() {
  static Pointtype point_typ;
  static Index3D I, Inext;
  static bool weiter, stop, no_27_S;
  static int t, i, j, k;
  iterate_hash1 {
    point_typ = point1->typ;
    if(point_typ != exterior && point_typ != multigrid &&
       point_typ != boundary) {
      // gut??
      point1->typ = interior;
      I = point1->ind;
      weiter = true;  // Annahme: es gibt einen Punkt in der Naehe
      for(t=I.Tiefe();weiter;++t) {
        no_27_S   = false;   // Annahme: es gibt 27-Stern
        weiter    = false;   // Annahme: nicht feinere Sterne untersuchen
        for(i=0;i<3;++i)     // fehlt: besser machen mit stop
          for(j=0;j<3;++j)
            for(k=0;k<3;++k) if(!((i==1) && (j==1) && (k==1)))  {
              Inext = I.next((Ort1D)i,(Ort1D)j,(Ort1D)k,t);
              if(Exists_Point(Inext)) {
                weiter    = true;          // Untersuche feinere Sterne
                point_typ = Give_type(Inext);
                // hier fehlt Test der Zellen!!!!!
                if(point_typ == exterior || point_typ == multigrid)
                  no_27_S=true;               // kein 27-Stern vorhanden
              }
              else {
                no_27_S   = true;             // stop kein 27-Stern vorhanden
              }
            }
        if(no_27_S==false) {
          weiter = false;                  // 27-Stern gefunden
          point1->level = t;
          if(t==I.Tiefe()) {
            point1->typ = hierarchical;
          }
          else {
            point1->typ = ausnahme;
          }
        }
      }
      if(no_27_S==true) {          // kein 27-Stern vorhanden
        point1->level = t-2;
        //  cout << "\n no: " << t; I.coordinate().Print();
      }
    }
  }
}
*/
/*
void Grid_base::Calc_slave_points() {
  iterate_hash1 {
    if(point1->typ == interior) {
      point1->typ = nearb;
    // vorerst sind alle Punkte hier nearb, keine slave points
    //  if(!Is_Slave(point1->ind)) {
    //  point1->typ = nearb;
    //      }
    }
  }
}
*/

/*
void Grid_base::Calc_nearb_points() {
  iterate_hash1 {
     if(point1->typ == nearb) {
       // Add_small_27Stencil(point1->ind,point1->ind.Tiefe());
       // Add_small_27Stencil(point1->ind,point1->level);
       //if(test<test_max) {
       //                Add_small_27Stencil(point1->ind,point1->ind.Tiefe());
       //        ++test;
       //    }
    }
  }
}
*/

bool Grid_base::Is_Slave(Index3D I) {
  bool slave_yn;
  Point_hashtable1* poi;

  if(I.Tiefe()>0) slave_yn = Calc_Is_Slave(I);
  else slave_yn = false;
  if(slave_yn==true) {
    //    cout << " slave ";
    poi = hashtable1_point(I);
    poi->typ = slave;
    poi->level = I.Tiefe();
  }
  return slave_yn;
}

bool Grid_base::Decide_poi_interpolates(Index3D I) {
  Pointtype p_typ;

  p_typ = Give_type(I);
  if((p_typ == exterior) || (p_typ == multigrid) || (p_typ == nearb))
    // if((p_typ == exterior) || (p_typ == multigrid))
    return false;
  return true;
}


bool Grid_base::Calc_Is_Slave(Index3D I) {
  int t, i, j, k;
  Index3D Ix, Iy, Iz;

  t = I.Tiefe();
  if(I.ind_x.Tiefe() == t && I.ind_x.Tiefe() > 0) {
    if(I.ind_y.Tiefe() < t && I.ind_z.Tiefe() < t) return true;
    for(i=0;i<2;++i) {
      Ix = I.neighbour_EW((dir1D)i);
      if(I.ind_y.Tiefe() == t && I.ind_y.Tiefe() > 0) {
        for(j=0;j<2;++j) {
          Iy = Ix.neighbour_NS((dir1D)j);
          if(I.ind_z.Tiefe() == t && I.ind_z.Tiefe() > 0) {
            for(k=0;k<2;++k) {
              Iz = Iy.neighbour_TD((dir1D)k);
              if(!Decide_poi_interpolates(Iz)) return false;
            }
          }
          else {
            if(!Decide_poi_interpolates(Iy)) return false;
          }
        }
      }
      else {
        // Es kann nur dies eintreten:
        if(I.ind_z.Tiefe() == t && I.ind_z.Tiefe() > 0) {
          for(k=0;k<2;++k) {
            Iz = Ix.neighbour_TD((dir1D)k);
            if(!Decide_poi_interpolates(Iz)) return false;
          }
        }
        else {
          if(!Decide_poi_interpolates(Ix)) return false;
        }
      }
    }
  }
  else {
    if(I.ind_y.Tiefe() == t && I.ind_y.Tiefe() > 0) {
      if(I.ind_z.Tiefe() < t) return true;
        for(j=0;j<2;++j) {
          Iy = I.neighbour_NS((dir1D)j);
          // kann nur eintreten:
          if(I.ind_z.Tiefe() == t && I.ind_z.Tiefe() > 0) {
            for(k=0;k<2;++k) {
              Iz = Iy.neighbour_TD((dir1D)k);
              if(!Decide_poi_interpolates(Iz)) return false;
            }
          }
          else {
            if(!Decide_poi_interpolates(Iy)) return false;
          }
        }
      }
      else {
        return true;
      }
  }
  return true;
}



void Grid_base::Put_finest_level_minimal(Index3D I, int level) {
  Point_hashtable1* poi;
  poi = hashtable1_point(I);
  if(poi->finest_level<level) poi->finest_level = level;
}

int Grid_base::Give_finest_level(Index3D I) {
  return hashtable1_point(I)->finest_level;
}

/*
void Halte_mal() {
  cout << " halt mal" << endl;
}

Index1D Index1D::neighbour(dir1D d) const   {
  static int ind;
  if(developer_version) {
    if(d==Ld && index==0) cout << " error in Index1D::neighbour! "  << endl;
    Halte_mal();
  }
  for(ind = index;(ind&1)==d;ind = ind >> 1) {};
  return Index1D(ind >> 1);
}
*/

int Grid_base::Give_my_finest_parallel_level(Index3D I) const {
  //  return hashtable1_point(I)->my_finest_parallel_level;
  Point_hashtable1* poi;
  poi = hashtable1_point(I);
  if(poi==NULL) return -1;
  return poi->my_finest_parallel_level;
}


// Abstand zum Rand in Richting d
double Grid_base::distanceD(D3vector V, dir_3D d) {
  return domain->distance(transform_coord(V),d);
};






#ifdef WITH_BRICK
#define DIMENSION 3


void Grid_base::getBrickInfo() {   

    /* Used in /Particle/ParticleBase.h to setup
       the required datastructure
    */

    D3vector v_WSD;
    D3vector v_ENT;
    
    BInfoT     brick;

    
    brickInfo_m.clear();
      
    // set cells in brick
    if(parallel_version) {
      iterate_hash_proc {
        if(point_proc!=NULL) {
          if(point_proc->Give_Index().Tiefe()==n_parallel) {
            v_WSD = transform_coord(
                                    point_proc->Give_Index().neighbour(WSDd).coordinate());
            v_ENT = transform_coord(
                                    point_proc->Give_Index().neighbour(ENTd).coordinate());
            
            brick.localXMin[0]  = v_WSD.x - finest_meshsize * 0.5;
            brick.localXMin[1]  = v_WSD.y - finest_meshsize * 0.5;
            brick.localXMin[2]  = v_WSD.z - finest_meshsize * 0.5;
            
            brick.localXMax[0]  = v_ENT.x - finest_meshsize * 0.5;
            brick.localXMax[1]  = v_ENT.y - finest_meshsize * 0.5;
            brick.localXMax[2]  = v_ENT.z - finest_meshsize * 0.5;
            
            brick.eps           = finest_meshsize * 0.001;
            brick.rank = point_proc->Give_num_proc();
            brickInfo_m.insert(brick);
            
            //    Index3D d = point_proc->Give_Index().neighbour(WSDd);
            //for(int d=0;d<8;++d)
            //     int n = give_next_rank_destination((dir_sons)d);
            //cout << "next rank WSDd" << n << endl; 
            
          }
        }
      }
    }
    else {


      v_WSD = transform_coord(Index3D(2,2,2).neighbour(WSDd).coordinate());
      v_ENT = transform_coord(Index3D(2,2,2).neighbour(ENTd).coordinate());
      
      brick.localXMin[0]  = v_WSD.x - finest_meshsize * 0.5;
      brick.localXMin[1]  = v_WSD.y - finest_meshsize * 0.5;
      brick.localXMin[2]  = v_WSD.z - finest_meshsize * 0.5;
      
      brick.localXMax[0]  = v_ENT.x - finest_meshsize * 0.5;
      brick.localXMax[1]  = v_ENT.y - finest_meshsize * 0.5;
      brick.localXMax[2]  = v_ENT.z - finest_meshsize * 0.5;
      
      brick.eps           = finest_meshsize * 0.001;
      brick.rank = 0;
      brickInfo_m.insert(brick);
    }
  }
  
void Grid_base::boundbox(double localXMin[DIMENSION], double localXMax[DIMENSION]) {
    Set_of_Processors::const_iterator i;
    for(unsigned int d=0; d<DIMENSION; d++) {
      localXMin[d] = 100000000.0; /* FIXME: VERY ugly */
      localXMax[d] =-100000000.0; /* FIXME: VERY ugly */
    }
    for(i=brickInfo_m.begin();i!=brickInfo_m.end();i++) {
      for(unsigned int d=0; d<DIMENSION; d++) {
        localXMin[d] = MIN(localXMin[d],(*i).localXMin[d]);
        localXMax[d] = MAX(localXMax[d],(*i).localXMax[d]);
      }
    }
  }
  
void Grid_base::localbox(double localXMin[DIMENSION], double localXMax[DIMENSION], int node) {
    
    Set_of_Processors::const_iterator i;
    for(i=brickInfo_m.begin();i!=brickInfo_m.end();i++) {
      if ((*i).rank==node) {
        for(unsigned int d=0; d<DIMENSION; d++) {
          localXMin[d] = (*i).localXMin[d];
          localXMax[d] = (*i).localXMax[d];
        }
      }
    }
  }


#endif // WITH_BRICK



































void Grid_base::Test_just_this2(int bla) {
  D3vector v;
  Index3D I;

  MPI_Barrier(MPI_COMM_WORLD);

  for(int i=0;i<8;++i) {
    if(my_rank==i) {
      cout << " \n rank: " << i << endl;
      iterate_hash1 {
        I = point1->Give_Index();
        if(Give_type(I)==multigrid 
           && I.Tiefe()<=bla
           && point1->finest_level>=bla 
           && ((point1->Give_label(5)>=bla) == true)
           ) {
        
//      if(Give_type(I)==multigrid && I==Index3D(0,4,0)) {

          //    if(Give_type(I)==multigrid) {
          //  v = transform_coord(I);
          v = I.coordinate();
          cout << " point: " << point1->Give_label(5) 
            //         << " my_rank: " << my_rank 
               << " resp: " << (I<my_index)
               << " vector: x: " << v.x << " y: " << v.y << " z: " << v.z 
               << " level: " << I.Tiefe()
               << " x: " << I.I_x().get()
               << " y: " << I.I_y().get()
               << " z: " << I.I_z().get()
               << " lab: " << point1->Give_label(5) 
               <<  endl;
        }
      }
    }
    MPI_Barrier(MPI_COMM_WORLD);
  }
}



void Grid_base::Test_just_this() {

  Index3D I;

  MPI_Barrier(MPI_COMM_WORLD);
  if(my_rank==0) cout << " Test_just_this: " << endl;
  MPI_Barrier(MPI_COMM_WORLD);



  MPI_Barrier(MPI_COMM_WORLD);
  if(my_rank==0) {
    cout << " procs von 0: " << endl; 
    iterate_hash_proc {
      cout << " processor "
           << " x: " << point_proc->Give_Index().I_x().coordinate()
           << " y: " << point_proc->Give_Index().I_y().coordinate()
           << " z: " << point_proc->Give_Index().I_z().coordinate()
           << endl;
    }
  }

  MPI_Barrier(MPI_COMM_WORLD);
  cout << " my_rank:  " << my_rank
           << " x: " << my_index.I_x().coordinate()
           << " y: " << my_index.I_y().coordinate()
           << " z: " << my_index.I_z().coordinate()
           << endl;
  MPI_Barrier(MPI_COMM_WORLD);

  /*
  // I = Index3D(0,4,8);
  // I = Index3D(0,2,0);
  // I = Index3D(9,8,16);
  //  I = Index3D(8,35,35);
  I = Index3D(37,17,33);
  // I = Index3D(37,35,8);
  if(my_rank==0) {
    v = I.coordinate();
    cout << " vector: x: " 
         << v.x << " y: " << v.y << " z: " << v.z 
         << endl;
  }

  if(I_am_active()) {
    if(Exists_Point(I)) {
      cout << " Point: type: " << Give_type(I)
           << " my rank: " << my_rank 
           << " in my finest index: " 
           << (I << Give_my_level_index(n_parallel))
           << " my_coarsest_level: " 
           << my_coarsest_level
           << " finest level: " 
           << Give_finest_level(I)
           << " Tiefe: " << I.Tiefe()
           << endl;
    }
  }
  MPI_Barrier(MPI_COMM_WORLD);
  */

  /*
  I = Index3D(9,8,18);
  if(Exists_Bo_Point(I,Edir)) {
    cout << " Bo Point: type: " << Give_type(I)
         << " my rank: " << my_rank << endl;
  }
  MPI_Barrier(MPI_COMM_WORLD);
  */  

  /*
  MPI_Barrier(MPI_COMM_WORLD);
  if(my_rank==0) {
    I = Index3D(37,38,36);
    //    v = transform_coord(I);
    v = I.coordinate();
    cout        << " Punkt D: " 
                << v.x << " y: " << v.y << " z: " << v.z 
                << endl;
  }
  */
  /*
  MPI_Barrier(MPI_COMM_WORLD);
  if(my_rank==0) {
    cout << " procs von 0: " << endl; 
    iterate_hash_proc {
      cout << " processor "
           << " x: " << point_proc->Give_Index().I_x().get()
           << " y: " << point_proc->Give_Index().I_y().get()
           << " z: " << point_proc->Give_Index().I_z().get()
           << endl;
    }
  }

  MPI_Barrier(MPI_COMM_WORLD);
  if(my_rank==7) {
    cout << " procs von 7: " << endl; 
    iterate_hash_proc {
      cout << " processor "
           << " x: " << point_proc->Give_Index().I_x().get()
           << " y: " << point_proc->Give_Index().I_y().get()
           << " z: " << point_proc->Give_Index().I_z().get()
           << endl;
    }
  }

  */
  
  /*
  MPI_Barrier(MPI_COMM_WORLD);
  //  v = transform_coord(my_index);
  v = my_index.coordinate();
  cout  << " my_rank: " << my_rank 
        << " vector: x: " 
        << v.x << " y: " << v.y << " z: " << v.z 
        << " x: " << my_index.I_x().get()
        << " y: " << my_index.I_y().get()
        << " z: " << my_index.I_z().get()
        << "Give_my_coarsest_level(): " << Give_my_coarsest_level() 
        << endl;
  MPI_Barrier(MPI_COMM_WORLD);

  int lll;
  lll = n_parallel-1;
  if(my_rank==0)  cout << " Grob: cell level " << lll << endl;

  MPI_Barrier(MPI_COMM_WORLD);
  //  v = transform_coord(my_index);
  if(Give_my_coarsest_level() <= lll-1) {
    v = Give_my_level_index(lll).coordinate();
    cout        << " my_rank: " << my_rank 
                << " vector: x: " 
                << v.x << " y: " << v.y << " z: " << v.z 
                << " x: " << Give_my_level_index(lll).I_x().get()
                << " y: " << Give_my_level_index(lll).I_y().get()
                << " z: " << Give_my_level_index(lll).I_z().get()
                << endl;
  }
  MPI_Barrier(MPI_COMM_WORLD);
  

  
  lll = n_parallel-2;
  if(my_rank==0)  cout << " Grob: cell level " << lll << endl;

  MPI_Barrier(MPI_COMM_WORLD);
  //  v = transform_coord(my_index);
  if(Give_my_coarsest_level() <= lll-1) {
    v = Give_my_level_index(lll).coordinate();
    cout        << " my_rank: " << my_rank 
                << " vector: x: " 
                << v.x << " y: " << v.y << " z: " << v.z 
                << " x: " << Give_my_level_index(lll).I_x().get()
                << " y: " << Give_my_level_index(lll).I_y().get()
                << " z: " << Give_my_level_index(lll).I_z().get()
                << endl;
  }
  MPI_Barrier(MPI_COMM_WORLD);
  */


  /*
  MPI_Barrier(MPI_COMM_WORLD);
  
  I = Index3D(4,2,19).next(ESTd,max_level+1).next(WSTd,max_level+1);
  v = transform_coord(I);
  if(Exists_Point(I)) {
    cout << " point: " << Give_type(I) 
         << " my_rank: " << my_rank 
         << " resposible: " << (I<my_index)
         << " vector: x: " << v.x << " y: " << v.y << " z: " << v.z 
         << " x: " << I.I_x().get()
         << " y: " << I.I_y().get()
         << " z: " << I.I_z().get()
         <<  endl;
  }
  */  
  
  
  /*
  //0.995 0.995 0.995
  iterate_hash1 {
    I = point1->ind;
    if(point1->typ==multigrid && I.Tiefe()==3) {
      v = transform_coord(point1->ind);
      
      cout << " point: " << Give_type(I) 
           << " my_rank: " << my_rank 
           << " resposible: " << (I << my_index)
           << " vector: x: " << v.x << " y: " << v.y << " z: " << v.z 
           << " x: " << I.I_x().get()
           << " y: " << I.I_y().get()
           << " z: " << I.I_z().get()
           <<  endl;
    }
  }
  */

  /*
  ofstream DATEI;
  // Print FE-solution in file with UCD_file format

  if(my_rank==0) {
    DATEI.open("../../point0.dat",ios :: out);
    Print_test_AVS(&DATEI,1);
    DATEI.close();
  }
  if(my_rank==1) {
    DATEI.open("../../point1.dat",ios :: out);
    Print_test_AVS(&DATEI,1);
    DATEI.close();
  }
  if(my_rank==2) {
    DATEI.open("../../point2.dat",ios :: out);
    Print_test_AVS(&DATEI,1);
    DATEI.close();
  }
  if(my_rank==3) {
    DATEI.open("../../point3.dat",ios :: out);
    Print_test_AVS(&DATEI,1);
    DATEI.close();
  }
  */

  /*
MPI_Barrier(MPI_COMM_WORLD);
 if(my_rank==0) {
   cout << " Problem at: " << endl;
   transform_coord(Index3D(17,17,17)).Print();
   cout << endl;
 }

MPI_Barrier(MPI_COMM_WORLD);

  if(Exists_Cell(Index3D(17,17,17))) {
     cout << " cell typ: " 
           << Give_cell_typ(Index3D(17,17,17))
           << " my_rank: " << my_rank << endl;
  }


  for(int r=0;r<20;++r) {
    MPI_Barrier(MPI_COMM_WORLD);

    if(my_rank==r) { 
      cout << endl;
      for(int i=0;i<8;++i) {
        I = Index3D(17,17,17).neighbour((dir_sons)i);
        v = transform_coord(I);
        if(Exists_Point(I)) {
          cout << " typ: " 
               << Give_type(I) 
               << " my_rank: " << my_rank 
               << " resposible: " << (I<my_index)
               << " vector: x: " << v.x << " y: " << v.y << " z: " << v.z 
               << " x: " << I.I_x().get()
               << " y: " << I.I_y().get()
               << " z: " << I.I_z().get()
               <<  endl;
        }
      }
    }
  }
  */

  /*
  int number_sp;
  number_sp=0;
  if(my_rank==3) 
    iterate_hash1 {
    if(point1->typ==hierarchical) {
      ++number_sp;
      v = transform_coord(point1->ind);
      if(my_rank==3) {
        if(v.y<0.45 && v.z<0.45)
          cout << " EEEErrror! " << endl;

        cout << " hiera.: x: " << v.x << " y: " << v.y << " z: " << v.z 
           << " my rank: : "  << my_rank << endl; 

        
      }
      
    }
  }
  
  cout << " my rank:   "  << my_rank
       << " number_sp: "  << number_sp << endl; 
  */


  /*
  if(my_rank==1) {
    if(Exists_Cell(Index3D(34,35,35))) {
      cout << " \n \n typ: " << Give_cell_typ(Index3D(34,35,35)) << endl;
    }
  }
  */

  /*
  for(int r=0;r<12;++r) {
    MPI_Barrier(comm);
    if(my_rank==r) {
      v = transform_coord(my_index);
      cout << " my index: x: " << v.x << " y: " << v.y << " z: " << v.z 
           << " my rank: : "  << my_rank << endl; 

      v = transform_coord(Give_my_level_index(n_parallel-1));
      cout << " unten: x: " << v.x << " y: " << v.y << " z: " << v.z 
           << " my rank: : "  << my_rank << endl; 


    }
  }
  */


  /*
  d = Wdir;

  for(int r=0;r<16;++r) {
    MPI_Barrier(comm);
    if(my_rank==r) {
      v = transform_coord(my_index);
      cout << " my index: x: " << v.x << " y: " << v.y << " z: " << v.z 
           << " my rank: : "  << my_rank << endl; 
      for(int l=0;l<=max_level;++l) {
        cout << " \n level: " << l << endl;
        cout << " next rank destination: " 
             << give_next_rank_destination(d,l) << endl;
        cout << " next rank source: " 
             << give_next_rank_source(d,l) << endl;
      }
    }
  }
  */
}

---