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

// ------------------------------------------------------------
// parallel.cc
//
// ------------------------------------------------------------

#define DEBUG false

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


// 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 5:
#include "gpar.h"
#include "parallel.h"

int MY_RANK;

Parallel_Info::~Parallel_Info() {
  for(int d=0;d<26;++d) {
    delete buffer_send[d];
    delete buffer_receive[d];
  }
  delete buffer_send_info;
  delete buffer_receive_info;

 //Loeschen des hashtables proc:
  int i_iter;
  Point_hashtable_proc* poi;
  for(i_iter=0;i_iter<hashtable_proc_leng;++i_iter) {
    poi=hashtable_proc_start[i_iter];
    delete poi;
  }
  delete hashtable_proc_start;

  Point_hashtable_proc_corner* poic;
  for(i_iter=0;i_iter<hashtable_proc_corner_leng;++i_iter) {
    poic=hashtable_proc_corner_start[i_iter];
    delete poic;
  }
  delete hashtable_proc_corner_start;
}

Parallel_Info::Parallel_Info(All_Domains* dom, MPI_Comm com) : 
  domain(dom), comm(com) {
  // -1.Step: refine_level is set in Grid_base::Grid_base
  // 0.Step: buffer initialization
  for(int d=0;d<26;++d) {
    length_send[d]    = 2;
    length_receive[d] = 2;
    buffer_send[d]    = new double[2];
    buffer_receive[d] = new double[2];
  }
  length_send_info = 2;
  length_receive_info = 2;
  buffer_send_info    = new int[2];
  buffer_receive_info = new int[2];

  // 1.step: MPI infos
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &p);

  // 2.step: initialize haschtable
  Initialize_hash_proc(2*p);
  Initialize_hash_proc_corner(2*p);

  MY_RANK=my_rank;
}

// constructor for testing how many processors are used
// no grid is generated
//------------
Parallel_Info::Parallel_Info(int n_max, All_Domains* dom,
                             int pp, Grid_gen_parameters& gpara) : 
  domain(dom) {
  //-1.Step: things usually done in Grid_base::Grid_base
    // for parallel:
  refine_level = gpara.Give_r_parallel();

  // calc offset
  double offset_square;
  double h_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));
  }

  // 0.Step: buffer initialization
  for(int d=0;d<26;++d) {
    length_send[d]    = 2;
    length_receive[d] = 2;
    buffer_send[d]    = new double[2];
    buffer_receive[d] = new double[2];
  }
  length_send_info = 2;
  length_receive_info = 2;
  buffer_send_info    = new int[2];
  buffer_receive_info = new int[2];

  // 1.step: MPI infos
  my_rank=0;
  p = pp;

  // 2.step: initialize haschtable
  Initialize_hash_proc(2*p);
  Initialize_hash_proc_corner(2*p);

  MY_RANK=my_rank;

  // 3.step: Generate_processes();
  Generate_processes();
}

Parallel_Info::Parallel_Info(double h_finest_level,  All_Domains* dom,
                             int pp, Grid_gen_parameters& gpara) : 
  domain(dom) {
  double h_offset, offset_square;
  int i;

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

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

  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;

  if(gpara.Is_there_grid_point()) {
    make_grid_with_grid_point(gpara.Give_grid_point(),h_finest_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;

  //-1.Step: things usually done in Grid_base::Grid_base


  // 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));
  }

  // 0.Step: buffer initialization
  for(int d=0;d<26;++d) {
    length_send[d]    = 2;
    length_receive[d] = 2;
    buffer_send[d]    = new double[2];
    buffer_receive[d] = new double[2];
  }
  length_send_info = 2;
  length_receive_info = 2;
  buffer_send_info    = new int[2];
  buffer_receive_info = new int[2];

  // 1.step: MPI infos
  my_rank=0;
  p = pp;

  // 2.step: initialize haschtable
  Initialize_hash_proc(2*p);
  Initialize_hash_proc_corner(2*p);

  MY_RANK=my_rank;

  // 3.step: Generate_processes();
  Generate_processes();
}

void Parallel_Info::process_info() {
  cout << " Information about processes: " << endl;
  cout << " Number of processes: " << give_number_of_processes() << endl;
  cout << " Number of active processes: " 
       << give_number_of_active_processes() << endl;
#ifdef PERIODIC
  cout << " Periodic BC in z selected " << endl;
#endif
  if(parallel_version == false)  {
    cout << " This is the serial version! " << endl;
  }
}

Index3D Parallel_Info::Give_level_index_of(int rang, int level) {
  if(level > n_parallel)
    level = n_parallel;
  iterate_hash_proc {
    if(point_proc->num_proc == rang)
      if(point_proc->Give_Index().Tiefe()==level)
        return point_proc->Give_Index();
  }

  // neu periodisch
  return Index3D(3,3,3);
}

/*
int Parallel_Info::Give_coarsest_level_of(int rang, int level) {
  iterate_hash_proc {
    if(point_proc->num_proc == rang)
      if(point_proc->Give_Index().Tiefe()==level)
        return point_proc->Give_Index();
  }

  // neu periodisch
  return Index3D(3,3,3);
}
*/

void Parallel_Info::Generate_processes() {
  Partitioning* pointer_parti;
  Partitioning* pointer_parti_old;
  rough_domain *pointer_rough_domain;
  int opt_I,     opt_J,     opt_K,     n_nec,     n_interior;
  int opt_I_old, opt_J_old, opt_K_old, n_nec_old, n_interior_old;
  int i;

  // 3.step: find minimal value for p_level: p_level_min
  // result is: p_level=0,1,2,3,...
  for(p_level=0;Zweipotenz(p_level*3)<=p;++p_level) {};
  --p_level;



  // 4.step: Partitioning
  pointer_parti = NULL;
  pointer_parti_old = NULL;
  for(n_nec=p-1;n_nec<p && p_level<max_level;++p_level) {
    // a) construct rough domain
    pointer_rough_domain = new rough_domain(p_level,refine_level,domain,
                                            A_bounding,H_bounding);
    n_interior_old       = n_interior;
    n_interior           = pointer_rough_domain->Give_n_interior();


    // b) save old partitioning and make new one
    n_nec_old = n_nec;
    opt_I_old = opt_I; opt_J_old = opt_J; opt_K_old = opt_K;
    if(pointer_parti_old!=NULL) delete(pointer_parti_old);
    pointer_parti_old = pointer_parti;    
    pointer_parti = new Partitioning(p_level);
    // c) calc number of necessary processes
    n_nec =
      pointer_rough_domain->Calc_optimal_partitioning(p,opt_I,opt_J,opt_K,
                                                      pointer_parti);

    delete(pointer_rough_domain);

    if(DEBUG) {
      MPI_Barrier(MPI_COMM_WORLD);
      if(my_rank==0) {
        
        cout << " Iteration: " 
             << " p_level: " << p_level
             << " max_level: " << max_level
             << " n_nec: " << n_nec << endl;
        cout << " partitioning: " << endl;
        pointer_parti->Print();
        
      }
    MPI_Barrier(MPI_COMM_WORLD);
    }

  }
  // take actual partitioning
  if(n_nec==p) {
    if(pointer_parti_old!=NULL) delete(pointer_parti_old);
    p_level=p_level-1;
  }
  // take old partitioning
  else {
    delete(pointer_parti);
    pointer_parti = pointer_parti_old;
    n_nec = n_nec_old;
    opt_I = opt_I_old; opt_J = opt_J_old; opt_K = opt_K_old;
    n_interior = n_interior_old;
    p_level=p_level-2;
  }
  if(MY_RANK==0) if(p_level>=max_level-2)
    cout << "\n Choose large value for n_level! \n" << endl;


  if(n_nec > p)
    cout << " Not implemented: Error in Parallel_Info::Generate_processes() "
         << endl;

  // 5. step: Set bounding box and ...
  n_parallel=p_level+2;

   // periodic
  // work
  if(domain->Is_periodic()) {
  // if(false) {
    n_parallel=p_level+1;

  }
  else {
    n_parallel=p_level+2;

    max_level =  max_level+1;
    A_bounding = A_bounding -
      H_bounding/(double)n_interior 
      * D3vector((double)opt_I,(double)opt_J,(double)opt_K);
    H_bounding = H_bounding*2.0;
  }

  if(MY_RANK==0) if(n_parallel>max_level-2)
    cout << "\n Error : Choose a finer grid or less processors! \n" << endl;

if(DEBUG) {  
    cout << " opt_I " << opt_I << endl;
    cout << " opt_J " << opt_J << endl;
    cout << " opt_K " << opt_K << endl;

    if(MY_RANK==0) 
      cout << "   take p_level: " <<  p_level 
           << " max_level:    "   << max_level
           << " p_level:      "   << p_level
           << " n_parallel:    "  << n_parallel << endl;
}

  if(domain->Is_periodic()) {
    // periodic
    Ixl = Index1D(0);
    Iyl = Index1D(0);
    Izl = Index1D(0);

    Ixr = Index1D(1);
    Iyr = Index1D(1);
    Izr = Index1D(1);
  }
  else {
    Ixl = Index1D(0);
    Iyl = Index1D(0);
    Izl = Index1D(0);
    // Berchnung des tree-index (Index1) aus kartesischen index (opt_I)
    for(i=0;i<opt_I;++i) Ixl = Ixl.next(Rd,p_level+refine_level+1);
    for(i=0;i<opt_J;++i) Iyl = Iyl.next(Rd,p_level+refine_level+1);
    for(i=0;i<opt_K;++i) Izl = Izl.next(Rd,p_level+refine_level+1);
    
    
    Ixr = Ixr.next(Rd,1);
    Iyr = Iyr.next(Rd,1);
    Izr = Izr.next(Rd,1);
    for(i=0;i<opt_I;++i) Ixr = Ixr.next(Rd,p_level+refine_level+1);
    for(i=0;i<opt_J;++i) Iyr = Iyr.next(Rd,p_level+refine_level+1);
    for(i=0;i<opt_K;++i) Izr = Izr.next(Rd,p_level+refine_level+1);
  }

  pxl = Ixl.coordinate();   pyl = Iyl.coordinate();   pzl = Izl.coordinate();
  pxr = Ixr.coordinate();   pyr = Iyr.coordinate();   pzr = Izr.coordinate();

  // 6.step: search cells on level max_level_parallel
  Index3D ind(2,2,2);

  if(n_parallel==1) {
    Add_proc_cell(ind);
  }
  else {
    Recursion_study_cells(ind, pointer_parti);
  } 

if(DEBUG) {
  MPI_Barrier(MPI_COMM_WORLD);
  if(my_rank==0) {
    cout << " partitioning: " << endl;
    pointer_parti->Print();
  }
  MPI_Barrier(MPI_COMM_WORLD);
}

  // 7.step: stop if number of processes is too small
  if(p < hashtable_proc_occ) {
    cout << " Number of processes is too small! "
         << p << ", active: "
         << hashtable_proc_occ << endl;
    return;
  }

  // 8.step: my_index
  i=0;
  iterate_hash_proc {
    if(point_proc!=NULL) {
      point_proc->num_proc = i; 
      if(i==my_rank) my_index = point_proc->Give_Index();
      ++i;
    }
  }

  // 9.   Step: Calc coarse processors
  // 9.1. Step: set number of active processors
  number_of_active_procs = hashtable_proc_occ;
  // 9.2. Step: calc coarse processors
  //            and coarser index and level
  Calc_coarse_processors();

  // 10.Step:    Set rank3D and edges und faces of my_index
  // 10.1. Step: Set Rank for 6 main direction
  Set_rank_3D();

  // 10.2. Step: Set Rank for special direction
  if(I_am_active()) {
    for(i=0;i<6;++i) {
      coarse_faces[i]      = my_index.surface((dir_3D)i);
      info_coarse_faces[i] = give_next_rank_destination((dir_3D)i,
                                                        n_parallel-1) >= 0;
    }
    for(i=0;i<12;++i) {
      coarse_edges[i]      = my_index.edge((Edges_cell)i);
      info_coarse_edges[i] = give_next_rank_destination((Edges_cell)i) >= 0;
    }
    for(i=0;i<8;++i) {
      coarse_corners[i]      = my_index.neighbour((dir_sons)i);
      info_coarse_corners[i] = give_next_rank_destination((dir_sons)i) >= 0;
      info_coarse_corners[i] = true;
    }
  }

  /*
  cout << " my_rank: " << my_rank
       << " T: " << info_coarse_faces[Tdir]
       << " D: " << info_coarse_faces[Ddir]
       << " W: " << info_coarse_faces[Wdir]
       << " E: " << info_coarse_faces[Edir]
       << " N: " << info_coarse_faces[Ndir]
       << " S: " << info_coarse_faces[Sdir] << endl;
  */

if(DEBUG) {
 MPI_Barrier(MPI_COMM_WORLD);
 cout.precision(4);
  cout << " my_rank: " << my_rank
       << " my old bounding_box: " 
       << " xl: " << pxl << " yl: " << pyl << " zl: " << pzl
       << " xr: " << pxr << " yr: " << pyr << " zr: " << pzr
       << endl;
 cout.precision(8);
 MPI_Barrier(MPI_COMM_WORLD);
}

  // 11.Step: Delete
  delete(pointer_parti);
  //  cout << " \n My index: ";  my_index.coordinate().Print(); cout << endl;





if(DEBUG) {
  // Processor hashtable
  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);
  D3vector cooo;
  if(my_rank==0) {
    cout << " Physical Informations: " << endl;
    A_bounding.Print();
    cout << " \n " 
         << H_bounding
         << endl;

    cout << " procs von 0 physical: " << endl; 

    iterate_hash_proc {
      cooo = coord_in_domain(point_proc->Give_Index());
      cout << " processor "
           << " x: " << cooo.x
           << " y: " << cooo.y
           << " z: " << cooo.z
           << 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()
       << " \n " 
       << " x: " << my_index.I_x().get()
       << " y: " << my_index.I_y().get()
       << " z: " << my_index.I_z().get()
       << endl;
  MPI_Barrier(MPI_COMM_WORLD);
}




}


void Parallel_Info::Generate_seriel_process() {
  int i;

  // 1. step: Set bounding box and ...
  pxl = 0.0;   pyl = 0.0;   pzl = 0.0;
  pxr = 1.0;   pyr = 1.0;   pzr = 1.0;

  // 2.step: my_index, my_rank
  Index3D ind(2,2,2);
  my_index = ind;

  my_rank=0;
  n_parallel = 0;
  my_coarsest_level = 0;

  // 9.Step:  Set edges und faces of my_index
  for(i=0;i<6;++i) {
    coarse_faces[i]      = my_index.surface((dir_3D)i);
    info_coarse_faces[i] = give_next_rank_destination((dir_3D)i,
                                                      n_parallel-1) >= 0;
  }
  for(i=0;i<12;++i) {
    coarse_edges[i]      = my_index.edge((Edges_cell)i);
    info_coarse_edges[i] = give_next_rank_destination((Edges_cell)i) >= 0;
  }
  for(i=0;i<8;++i) {
    coarse_corners[i]      = my_index.neighbour((dir_sons)i);
    info_coarse_corners[i] = give_next_rank_destination((dir_sons)i) >= 0;
    info_coarse_corners[i] = true;
  }
}

void Parallel_Info::make_grid_with_grid_point(D3vector gp, double fms) {
  double new_q;
  D3vector move_v;
  if(A_bounding.x >= gp.x ||
     A_bounding.y >= gp.y ||
     A_bounding.z >= gp.z ||
     A_bounding.x+H_bounding <= gp.x ||
     A_bounding.y+H_bounding <= gp.y ||
     A_bounding.z+H_bounding <= gp.z)  {
    cout << " \n Error: Grid point is too far away from domain! \n" << endl;
  }
  else {
    for(new_q=gp.x;A_bounding.x<new_q;new_q=new_q-fms) {};
    move_v.x = A_bounding.x-new_q;

    for(new_q=gp.y;A_bounding.y<new_q;new_q=new_q-fms) {};
    move_v.y = A_bounding.y-new_q;

    for(new_q=gp.z;A_bounding.z<new_q;new_q=new_q-fms) {};
    move_v.z = A_bounding.z-new_q;

    H_bounding = H_bounding + L_infty(move_v);
    A_bounding = A_bounding - move_v;
  }
}

bool Parallel_Info::contained_in_old_bounding_box(Index3D I) const {
  if(I.I_x()==Ixl || I.I_x()==Ixr || 
     I.I_y()==Iyl || I.I_y()==Iyr || 
     I.I_z()==Izl || I.I_z()==Izr) return false;
  double s;
  s = I.I_x().coordinate();
  if(s<pxl || s>pxr) return false;
  s = I.I_y().coordinate();
  if(s<pyl || s>pyr) return false;
  s = I.I_z().coordinate();
  if(s<pzl || s>pzr) return false;
  return true;
}


// local and old bounding
bool Parallel_Info::contained_in_boxes(Index3D I) const {  
  int i;
  // weg falls periodic
#ifndef PERIODIC
   if(contained_in_old_bounding_box(I)) {
#endif
  
    if(my_index<=I) return true;
    for(i=0;i<6;++i) {
      if(coarse_faces[i]>>I) {
        if(info_coarse_faces[i]==false) return false;
        else if(coarse_faces[i]<=I)     return true;
      }
    }
    for(i=0;i<12;++i) {
      if(coarse_edges[i]>>I) {
        if(info_coarse_edges[i]==false) return false;
        else if(coarse_edges[i]<=I)     return true;
      }
    }
    for(i=0;i<8;++i) {
      if(coarse_corners[i]==I) {
        if(info_coarse_corners[i]==false) return false;
      }
    }


    return true;
#ifndef PERIODIC
   }
   return false;
#endif
}



int Parallel_Info::give_next_rank_source(Edges_cell ed) { // -1 if not exists
  int t;
  Index3D next, sp;
  
  ed = opposite3D(ed);

  t = my_index.Tiefe();
  sp = my_index.next(ed,t);

#ifndef PERIODIC
  if(sp.I_x()==Index1D(0) || sp.I_y()==Index1D(0) || sp.I_z()==Index1D(0) ||
     sp.I_x()==Index1D(1) || sp.I_y()==Index1D(1) || sp.I_z()==Index1D(1)) 
     return -1;
#endif

  next = sp.next(ed,t);

  iterate_hash_proc {
    if(point_proc!=NULL)
      if(point_proc->Give_Index() == next) {
        if(point_proc->num_proc == my_rank) return -1;
        return point_proc->num_proc; 
      }
  }
  return -1;
};

int Parallel_Info::give_next_rank_destination(Edges_cell ed) { 
                                              // -1 if not exists
  int t;
  Index3D next, sp;

  t = my_index.Tiefe();
  sp = my_index.next(ed,t);
#ifndef PERIODIC
  if(sp.I_x()==Index1D(0) || sp.I_y()==Index1D(0) || sp.I_z()==Index1D(0) ||
     sp.I_x()==Index1D(1) || sp.I_y()==Index1D(1) || sp.I_z()==Index1D(1)) 
     return -1;
#endif

  next = sp.next(ed,t);

  iterate_hash_proc {
    if(point_proc!=NULL)
      if(point_proc->Give_Index() == next) {
        if(point_proc->num_proc == my_rank) return -1;
        return point_proc->num_proc; 
      }
  }
  return -1;
};


int Parallel_Info::give_next_rank_destination(dir_sons d) { 
                                              // -1 if not exists
  int t;
  Index3D next, sp;

  t = my_index.Tiefe();
  sp = my_index.next(d,t);

#ifndef PERIODIC
  if(sp.I_x()==Index1D(0) || sp.I_y()==Index1D(0) || sp.I_z()==Index1D(0) ||
     sp.I_x()==Index1D(1) || sp.I_y()==Index1D(1) || sp.I_z()==Index1D(1)) 
    return -1;
#endif

  next = sp.next(d,t);

  iterate_hash_proc {
    if(point_proc!=NULL)
      if(point_proc->Give_Index() == next) {
        if(point_proc->num_proc == my_rank) return -1;
        return point_proc->num_proc; 
      }
  }
  return -1;
};

int Parallel_Info::give_next_rank_source(dir_sons d) { 
                                              // -1 if not exists
  int t;
  Index3D next, sp;
  dir_sons dd;

  dd = d;
  d = opposite3D(dd);

  t = my_index.Tiefe();
  sp = my_index.next(d,t);

#ifndef PERIODIC
    if(sp.I_x()==Index1D(0) || sp.I_y()==Index1D(0) || sp.I_z()==Index1D(0) ||
    sp.I_x()==Index1D(1) || sp.I_y()==Index1D(1) || sp.I_z()==Index1D(1)) 
    return -1;
#endif

  next = sp.next(d,t);

  iterate_hash_proc {
    if(point_proc!=NULL)
    if(point_proc->Give_Index() == next) {
        if(point_proc->num_proc == my_rank) return -1;
        return point_proc->num_proc; 
      }
  }
  return -1;
};

/*
void Parallel_Info::Set_coarser_index_and_level() {
  int lev;
  Index3D origin;
  // a) coarsest possible level
  origin = my_index.neighbour(WSDd);
  my_coarsest_level = origin.Tiefe();
  if(my_coarsest_level < 1)
    my_coarsest_level = 1;

  // b) my index on coarser levels
  my_level_index = new Index3D[n_parallel];
  // - set for too coarser levels
  for(lev=0;lev<=my_coarsest_level;++lev) {
    my_level_index[lev] = my_index;
  }
  // - calculate rest coarse levels
  for(lev=my_coarsest_level+1;lev<n_parallel;++lev) {
    my_level_index[lev] = origin.next_ENT(lev);
  }
}
*/

Index3D Parallel_Info::Give_next_on_level(int d,          // 0 <= d < 26
                                          Index3D me) {
  int t;
  Index3D nextS;

  t = me.Tiefe();
    
  if(d<6) {
    nextS = me.next((dir_3D)d,t);
    nextS = nextS.next((dir_3D)d,t);
  }
  else if(d<18) {
    nextS = me.next((Edges_cell)(d-6),t);
    nextS = nextS.next((Edges_cell)(d-6),t);
  }
  else {
    nextS = me.next((dir_sons)(d-18),t);
    nextS = nextS.next((dir_sons)(d-18),t);
  }
  return nextS;
}

void Parallel_Info::Calc_coarse_processors() {
  Index3D next, nnext;
  int i, l, rrr;
  bool found;
  Point_hashtable_proc* ppp;

  dir_sons priority[8];
  priority[0] = ENTd;  
  priority[1] = ENDd;   priority[2] = ESTd;   priority[3] = WNTd;
  priority[4] = ESDd;   priority[5] = WNDd;   priority[6] = WSTd;
  priority[7] = WSDd;

  min_level = 1;

  // new
  // 1. Calc all coarse interior processors
  Index3D ind(2,2,2);
  Recursion_calc_coarse_processors(ind);

  // 2. Calc corners of processors
  iterate_hash_proc {
    ind = point_proc->Give_Index();
    l = ind.Tiefe();
    for(i=0;i<8;++i)
      Add_proc_corner(ind.neighbour((dir_sons)i),l-1);
  }

  for(l=n_parallel;l>=min_level+1;--l) {
    // 3 set: used_on_coarser_grid=false
    iterate_hash_proc {
      point_proc->used_on_coarser_grid=false;
    }

    // 4. study proc corner, which exists also on coarser level
    /*

     idea of priority concept

     *---*---*
     | | | | |
     |---|---|
     | |2|1| |
     *---I---*
     | |4|3| |
     |---|---|
     | | | | |
     *---*---*

     */


    iterate_hash_proc_corner { 
      ind = point_proc_corner->ind;

      if(ind.Tiefe()<l-1 &&
         point_proc_corner->maximal_level>=l-1) {
        found=false;

        for(i=0;i<8 && found==false;++i) {
          next = ind.next(priority[i],l);
          if(Exists_proc(next)) {
            found = true;

            // work on fine proc
            ppp = hashtable_proc_point(next);
            ppp->used_on_coarser_grid=true;
            rrr = ppp->num_proc;

            // work on coarse proc
            next = ind.next(ENTd,l-1);
            if(Exists_proc(next)==false) {
              Add_proc_cell(next); // add  processor
              ppp = hashtable_proc_point(next);
              ppp->interior_processor = false; // make it exterior
            }
            else 
              ppp = hashtable_proc_point(next);
            ppp->num_proc = rrr;

            // work on proc corner
            point_proc_corner->num_proc = rrr;
          }
        }
      }
    }

    // 5. study proc corner, which exists only on this level
    iterate_hash_proc_corner { 
      ind = point_proc_corner->ind;
      if(ind.Tiefe()==l-1 &&
         point_proc_corner->maximal_level>=l-1) {
        found=false;
        for(i=0;i<8 && found==false;++i) {
          next = ind.next(priority[i],l);
          if(Exists_proc(next)) {
            found = true;

            // work on fine proc
            ppp = hashtable_proc_point(next);
            rrr = ppp->num_proc;

            // work on proc corner
            point_proc_corner->num_proc = rrr;
          }
        }
        if(developer_version) {
          if(found==false)
            cout << " error in Calc_coarse_processors(). 5. " << endl;
        }
      }
    }

    /*

      find out which of the fine proc 1,2,3,4
      are not used

     *---*
     |3 4|
     | I |
     |1 2|
     *---*

     */    


    // 6. study coarser procs
    iterate_hash_proc {
      ind = point_proc->Give_Index();
      if(ind.Tiefe()==l-1) {
        if(point_proc->num_proc <0) {
          // it is necessary to look for an unused processor 
          found=false;
          for(i=0;i<8 && found==false;++i) {
            next = ind.son((dir_sons)i); // zuerst WSDd
            if(Exists_proc(next)) {
              ppp = hashtable_proc_point(next);
              if(ppp->used_on_coarser_grid==false) {
                // found one
                found=true;
                point_proc->num_proc = ppp->num_proc;
              }
            }
          }
          if(found==false) {
            cout << " \n \n rise minimal level to: "
                 << l-1 << endl;
            if(min_level<l-1)
              min_level = l-1;
          }
        }
      }
    }
  }



  // 7. Calc my_level_index and my_coarsest_level
  if(I_am_active()) {
    int level;
    my_level_index = new Index3D[n_parallel];
    my_coarsest_level = n_parallel-1;
    for(level=n_parallel-1;level>min_level;level--) {
      my_level_index[level] = my_index.son(ENTd); // this means empty
      iterate_hash_proc {
        if(point_proc!=NULL)
          if(point_proc->Give_Index().Tiefe() == level &&
             point_proc->num_proc==my_rank) {
            my_level_index[level] = point_proc->Give_Index();
            my_coarsest_level = level-1;
          }
      }
    }
    if(my_coarsest_level<1)
      my_coarsest_level=1;
  }


  //  Print_processors();
}


void Parallel_Info::Recursion_calc_coarse_processors(Index3D I) {
  int level, i;
  bool no_sons;

  level = I.Tiefe();
  
  if(level<n_parallel-1) {
    // recursion
    for(i=0;i<8;++i) 
      Recursion_calc_coarse_processors(I.son((dir_sons)i));
  }
  // calc number of processor for this coarse cell
  no_sons = true;
  for(i=0;i<8;++i)
    if(exists_process(I.son((dir_sons)i))) no_sons = false;
  if(no_sons == false) {
    //  if(I.Tiefe()>1)
    Add_proc_cell(I); // add interior processor
  }
}

Index3D Parallel_Info::Give_my_level_index(int level) {
  if(level>=n_parallel) return my_index;
  return my_level_index[level];
}

void Parallel_Info::Print_processors() { 
  if(my_rank==0) {
    for(int l=n_parallel;l>=1;--l) {
      cout << " Processors on level: " << l << endl;
      iterate_hash_proc {
        if(point_proc->Give_Index().Tiefe() == l) {
          point_proc->Give_Index().coordinate().Print();
          cout << " rank: " << point_proc->num_proc << endl;
        }
      }
    }
  }
}

void Parallel_Info::Set_rank_3D() { 
  Index3D next, sp;
  bool stop;
  int d, num_pp;
  int lev;

  if(I_am_active()) {
    // 1. Part: faces
    //  **************
    // a) next processor for level >= n_parallel-1
    // rank on finest level
    for(d=0;d<6;++d) {
      stop=false;
      sp = my_index.next((dir_3D)d,my_index.Tiefe());

#ifndef PERIODIC
      if(sp.I_q(Main_direction((dir_3D)d))==OneDpart((dir_3D)d)) {
        next_rank_destination[d]        = -1;
        next_rank_source[opposite3D((dir_3D)d)] = -1;
        stop = true;
      }
#endif
      next = sp.next((dir_3D)d,my_index.Tiefe());
      
      if(stop==false) {
        iterate_hash_proc {
          if(point_proc!=NULL)
            if(point_proc->Give_Index() == next) {
              // processor found
              next_rank_destination[d]        = point_proc->num_proc;
              next_rank_source[opposite3D((dir_3D)d)] = point_proc->num_proc;
              stop = true;
            }
        }
      }
      if(stop==false) {
        // no processor found
        next_rank_destination[d]=-1;
        next_rank_source[opposite3D((dir_3D)d)]=-1;
      }
    }
    
    // b) next processor for level < n_parallel-1
    for(d=0;d<6;++d) {
      next_rank_source_coarse[d]      = new int[n_parallel-1];
      next_rank_destination_coarse[d] = new int[n_parallel-1];
    }
    // rank on coarser level
    // - set -1 for too coarser levels
    for(d=0;d<6;++d) for(lev=0;lev<my_coarsest_level;++lev) {
      next_rank_source_coarse[d][lev] = -1;
      next_rank_destination_coarse[d][lev] = -1;
    }
    // - calculate rest coarse levels
    for(d=0;d<6;++d) for(lev=my_coarsest_level;lev<n_parallel-1;++lev) {
      sp = my_level_index[lev+1].next((dir_3D)d,lev+1);
      // periodisch weg
#ifndef PERIODIC
      if(sp.I_q(Main_direction((dir_3D)d))==OneDpart((dir_3D)d)) {
        next_rank_destination_coarse[d][lev]                = -1;
        next_rank_source_coarse[opposite3D((dir_3D)d)][lev] = -1;
      }
      else
#endif
        {
          next = sp.next((dir_3D)d,lev+1);
          num_pp = Give_proc_number(next);
          next_rank_destination_coarse[d][lev]                = num_pp;
          next_rank_source_coarse[opposite3D((dir_3D)d)][lev] = num_pp;
        }
    }
    
    // 2. Part: edges
    // **************
    // a) next processor for level >= n_parallel-1
    // rank on finest level
    for(d=0;d<12;++d) {
      next_rank_source_edges[d] =
        give_next_rank_source((Edges_cell)d);
      next_rank_destination_edges[d] =
        give_next_rank_destination((Edges_cell)d);
    }

    // b) next processor for level < n_parallel-1
    for(d=0;d<12;++d) {
      next_rank_source_coarse_edges[d]      = new int[n_parallel-1];
      next_rank_destination_coarse_edges[d] = new int[n_parallel-1];
    }
    // rank on coarser level
    // - set -1 for too coarser levels
    for(d=0;d<12;++d) for(lev=0;lev<my_coarsest_level;++lev) {
      next_rank_source_coarse_edges[d][lev] = -1;
      next_rank_destination_coarse_edges[d][lev] = -1;
    }
    // - calculate rest coarse levels
    for(d=0;d<12;++d) for(lev=my_coarsest_level;lev<n_parallel-1;++lev) {
      sp = my_level_index[lev+1].next((Edges_cell)d,lev+1);
      next = sp.next((Edges_cell)d,lev+1);

      num_pp = Give_proc_number(next);
      next_rank_destination_coarse_edges[d][lev]                = num_pp;
      next_rank_source_coarse_edges[opposite3D((Edges_cell)d)][lev] = num_pp;
    }
    

    // 2. Part: corners
    // **************
    // a) next processor for level >= n_parallel-1
    // rank on finest level
    for(d=0;d<8;++d) {
      next_rank_source_corners[d] =
        give_next_rank_source((dir_sons)d);
      next_rank_destination_corners[d] =
        give_next_rank_destination((dir_sons)d);
    }

    // b) next processor for level < n_parallel-1
    for(d=0;d<8;++d) {
      next_rank_source_coarse_corners[d]      = new int[n_parallel-1];
      next_rank_destination_coarse_corners[d] = new int[n_parallel-1];
    }
    // rank on coarser level
    // - set -1 for too coarser levels
    for(d=0;d<8;++d) for(lev=0;lev<my_coarsest_level;++lev) {
      next_rank_source_coarse_corners[d][lev] = -1;
      next_rank_destination_coarse_corners[d][lev] = -1;
    }
    // - calculate rest coarse levels
    for(d=0;d<8;++d) for(lev=my_coarsest_level;lev<n_parallel-1;++lev) {
      sp = my_level_index[lev+1].next((dir_sons)d,lev+1);
      next = sp.next((dir_sons)d,lev+1);
      num_pp = Give_proc_number(next);
      next_rank_destination_coarse_corners[d][lev] 
        = num_pp;
      next_rank_source_coarse_corners[opposite3D((dir_sons)d)][lev]
        = num_pp;
    }




    /*


iterate_hash_proc {
  if(point_proc->Give_Index().Tiefe() == n_parallel) {
    cout << " rank: " << point_proc->num_proc << " , and: ";
    point_proc->Give_Index().coordinate().Print();
    cout << endl;
  }
}


  cout << " n_parallel: " << n_parallel
       << " my_coarsest_level: " << my_coarsest_level
       << " my_rank: " <<  my_rank
       << " \n source_corners[WST]: " 
       << next_rank_source_coarse_corners[WSTd][lev]
       << " des_corners[WST]: " 
       << next_rank_destination_coarse_corners[WSTd]
       << " source_corners[END]: "
       << next_rank_source_coarse_corners[ENDd]
       << " des_corners[END]: " 
       << next_rank_destination_coarse_corners[ENDd]
       << endl;
    */


    


  }
};

int* Parallel_Info::Give_buffer_send_info(int l) {
  if(length_send_info<l) {
    delete(buffer_send_info);
    buffer_send_info = new int[l];
  }
  return buffer_send_info;
}

int* Parallel_Info::Give_buffer_receive_info(int l) {
  if(length_receive_info<l) {
    delete(buffer_receive_info);
    buffer_receive_info = new int[l];
  }
  return buffer_receive_info;
}


double* Parallel_Info::Give_buffer_send(int d, int l) {
  if(length_send[d]<l) {
    delete(buffer_send[d]);
    buffer_send[d] = new double[l];
  }
  return buffer_send[d];
}

double* Parallel_Info::Give_buffer_receive(int d, int l) {
  if(length_receive[d]<l) {
    delete(buffer_receive[d]);
    buffer_receive[d] = new double[l];
  }
  return buffer_receive[d];
}

  

void Parallel_Info::Recursion_study_cells(Index3D I, Partitioning* parti) {
  int i;
  Index3D I_neigh;

  if(I.Tiefe() < n_parallel) {  // level too small-> recursion
    for(i=0;i<8;++i)
      Recursion_study_cells(I.son((dir_sons)i),parti);
  }
  else {                              // level we look for
    // look to neighbours
    if(parti->Is_this_a_proc(I)) {
      Add_proc_cell(I);
    }
  }
};

int Parallel_Info::Give_proc_number(Index3D I) {
  Point_hashtable_proc* ppp;
  ppp = hashtable_proc_point(I);
  if(ppp==NULL) return -1;
  return ppp->num_proc;
}

bool Parallel_Info::Exists_proc(Index3D I) {
  Point_hashtable_proc* ppp;
  ppp = hashtable_proc_point(I);
  if(ppp==NULL) return false;
  return true;
}

int Parallel_Info::Give_proc_corner_number(Index3D I) {
  Point_hashtable_proc_corner* ppp;
  ppp = hashtable_proc_corner_point(I);
  if(ppp==NULL) return -1;
  return ppp->num_proc;
}

bool Parallel_Info::Exists_proc_corner(Index3D I) {
  Point_hashtable_proc_corner* ppp;
  ppp = hashtable_proc_corner_point(I);
  if(ppp==NULL) return false;
  return true;
}

void Parallel_Info::Add_proc_cell(Index3D I) {
  static Point_hashtable_proc* point;
  point = hashtable_proc_start
    [hashtable_proc_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,
                         hashtable_proc_leng)];
  if(point==NULL) {
    hashtable_proc_occ++;
    point = new Point_hashtable_proc(I);
    hashtable_proc_start[hashtable_proc_function(I.ind_x.index,I.ind_y.index,
                                                 I.ind_z.index,
                                                 hashtable_proc_leng)] = point;
  }
  else {
    while(point->next!=NULL && point->ind!=I) point = point->next;
    if(point->ind!=I) {
      hashtable_proc_occ++;
      point->next = new Point_hashtable_proc(I);
    }
  }
}

void Parallel_Info::Add_proc_corner(Index3D I, int l) {
  Point_hashtable_proc_corner* point;
  Point_hashtable_proc_corner* pp;
  point = hashtable_proc_corner_start
    [hashtable_proc_corner_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,
                         hashtable_proc_corner_leng)];
  if(point==NULL) {
    hashtable_proc_corner_occ++;
    point = new Point_hashtable_proc_corner(I);
    hashtable_proc_corner_start[hashtable_proc_corner_function(I.ind_x.index,
                                                        I.ind_y.index,
                                                        I.ind_z.index,
                                      hashtable_proc_corner_leng)] = point;
    pp = point;
  }
  else {
    while(point->next!=NULL && point->ind!=I) point = point->next;
    if(point->ind!=I) {
      hashtable_proc_corner_occ++;
      point->next = new Point_hashtable_proc_corner(I);
      pp = point->next;
    }
    else 
      pp = point;
  }
  if(pp->maximal_level<l)
    pp->maximal_level=l;
}

// Elementares fuer hashtable_proc
// ----------------------------
void Parallel_Info::Initialize_hash_proc(int lenght) {
  int i;
  hashtable_proc_leng = Find_next_prime(lenght);
  hashtable_proc_start = new Point_hashtable_proc*[hashtable_proc_leng];
  for(i=0;i<hashtable_proc_leng;++i)
    hashtable_proc_start[i] = NULL;
  hashtable_proc_occ=0;
}

void Parallel_Info::Initialize_hash_proc_corner(int lenght) {
  int i;
  hashtable_proc_corner_leng = Find_next_prime(lenght);
  hashtable_proc_corner_start = 
    new Point_hashtable_proc_corner*[hashtable_proc_corner_leng];
  for(i=0;i<hashtable_proc_corner_leng;++i)
    hashtable_proc_corner_start[i] = NULL;
  hashtable_proc_corner_occ=0;
}


//  rough domain
// ---------------------------------
inline void rough_domain::set(int i, int j, int k) {
  start[i+l*j+ll*k] = true;
}

rough_domain::rough_domain(int P_level, int Refine_level, All_Domains* dom,
                           D3vector A_bounding, double H_bounding) :
  p_level(P_level), refine_level(Refine_level) {
  int i,j,k;
  D3vector A;
  double h;  
  bool nothing;

  // periodic
  is_periodic = dom->Is_periodic();

  // alloce storage
  n_interior = Zweipotenz(p_level+refine_level);
  n_side     = Zweipotenz(refine_level);
  l  = n_interior;
  ll = l*l;

  start = new bool[l*l*l];
  for(i=0;i<l*l*l;++i) start[i]=false;

  // set entries of rough_domain
  A = A_bounding;
  h = H_bounding / (double)n_interior;

  nothing=true;
  for(i=1;i<n_interior;++i) 
    for(j=1;j<n_interior;++j) 
      for(k=1;k<n_interior;++k) {
        if(dom->point_in_domain(D3vector(h*i,h*j,h*k)+A)) {
          nothing=false;
          set(i-1,j  ,k);         set(i,j  ,k);
          set(i-1,j-1,k);         set(i,j-1,k);

          set(i-1,j  ,k-1);       set(i,j  ,k-1);
          set(i-1,j-1,k-1);       set(i,j-1,k-1);
        }
      }

  if(nothing) {
    for(i=0;i<n_interior;++i)
      for(j=0;j<n_interior;++j)
        for(k=0;k<n_interior;++k)
          start[i+l*j+ll*k] =true;
  }

  /*
  cout << " \n Matrix: ";
  for(i=0;i<=n_interior;++i) {
    cout << " \n next:";
    for(j=0;j<=n_interior;++j) {
      cout << endl;
      for(k=0;k<=n_interior;++k) {
        cout << start[i+l*j+ll*k] << ", ";
      }
    }
  }
  cout << " \n : n_interior: " << n_interior
       << " \n : n_side: "     << n_side       << endl;
  */
};

rough_domain::~rough_domain() {
  delete(start);
}

bool rough_domain::In_domain(int i, int j, int k) {
  if(developer_version) {
    if(i<1-n_side          || j<1-n_side          || k<1-n_side          ||
       i>n_interior+n_side || j>n_interior+n_side || k>n_interior+n_side) {
      cout << " \n error in rough_domain::In_domain!!"  
           << " i: " << i << " j: " << j << " k: " << k << endl;
    }
  }
  if(j<0          || i<0          || k<0   || 
     i>=n_interior || j>=n_interior || k>=n_interior) {
    return false;
  }
  return start[i+l*j+ll*k];
};


/*
int rough_domain::p_necessary(int I, int J, int K, int p) {
  int i,j,k,n_nec;
  int iS,jS,kS;
  int ii,jj,kk;
  bool intersec;

  n_nec=0;
  for(i=0;i<=n_interior/n_side;++i) 
    for(j=0;j<=n_interior/n_side;++j)  
      for(k=0;k<=n_interior/n_side;++k) {
        iS=-I+i*n_side;
        jS=-J+j*n_side;
        kS=-K+k*n_side;

        intersec=false;
        for(ii=0;ii<=n_side && intersec==false;++ii) 
          for(jj=0;jj<=n_side && intersec==false;++jj) 
            for(kk=0;kk<=n_side && intersec==false;++kk) {
              if(In_domain(iS+ii,jS+jj,kS+kk)) {
                intersec=true;
              }
            }

        if(intersec==true) {
          ++n_nec;
        }
        if(n_nec>p) return n_nec;
      }
  return n_nec;
};
*/

quality rough_domain::calc_quality(int I, int J, int K, int p) {
  int i,j,k,n_nec;
  int iS,jS,kS;
  int ii,jj,kk;
  int num_inter_sec, minimal_inter_sec, maximal_inter_sec;
  bool intersec;

  minimal_inter_sec = (n_side+1)*(n_side+1)*(n_side+1);
  maximal_inter_sec = 0;
  n_nec=0;
  for(i=0;i<=n_interior/n_side;++i) 
    for(j=0;j<=n_interior/n_side;++j)  
      for(k=0;k<=n_interior/n_side;++k) {
        iS=-I+i*n_side;
        jS=-J+j*n_side;
        kS=-K+k*n_side;

        num_inter_sec = 0;
        intersec      = false;
        for(ii=0;ii<n_side;++ii) 
          for(jj=0;jj<n_side;++jj) 
            for(kk=0;kk<n_side;++kk) {
              if(In_domain(iS+ii,jS+jj,kS+kk)) {
                intersec=true;
                //      if(Full_in_domain(iS+ii,jS+jj,kS+kk))
                  ++num_inter_sec;
              }
            }
        if(intersec) {
          ++n_nec;
          if(num_inter_sec<minimal_inter_sec)
            minimal_inter_sec=num_inter_sec;
          if(num_inter_sec>maximal_inter_sec)
            maximal_inter_sec=num_inter_sec;
        }
        if(n_nec>p) return quality();
      }

  // cout << " intersec: " << minimal_inter_sec << endl;

  return quality(n_nec,minimal_inter_sec,maximal_inter_sec); // yes quality 
  // return quality(n_nec,0);  // no quality 
};

void rough_domain::Set_partitioning(int I, int J, int K, Partitioning* parti) {
  int i,j,k;
  int iS,jS,kS;
  int ii,jj,kk;
  bool intersec;

  /*
  cout << " \n : n_interior: " << n_interior
       << " \n : n_side: "     << n_side       << endl;
  */

  for(i=0;i<=n_interior/n_side;++i) 
    for(j=0;j<=n_interior/n_side;++j)  
      for(k=0;k<=n_interior/n_side;++k) {
        iS=-I+i*n_side;
        jS=-J+j*n_side;
        kS=-K+k*n_side;

        intersec=false;
        for(ii=0;ii<n_side && intersec==false;++ii) 
          for(jj=0;jj<n_side && intersec==false;++jj) 
            for(kk=0;kk<n_side && intersec==false;++kk) {
              if(In_domain(iS+ii,jS+jj,kS+kk)) {
                intersec=true;
              }
            }

        if(intersec==true) {
          parti->Put_proc(i,j,k);
        }
      }
};


int rough_domain::Calc_optimal_partitioning(int p,
                                            int& opt_I,int& opt_J,int& opt_K,
                                            Partitioning* parti) {
  int I,J,K;

  quality q_new, q_opt;


  if(is_periodic) {
  // if(false) {
    I = 0; J = 0; K = 0;
    opt_I=0; opt_J=0; opt_K=0;

    q_new=calc_quality(opt_I,opt_J,opt_K,p);
    if(q_opt<q_new) {
      q_opt=q_new;
      opt_I=I; opt_J=J; opt_K=K;
    }
  }
  else {
    for(I=0;I<n_side;++I) for(J=0;J<n_side;++J) for(K=0;K<n_side;++K) {
      q_new=calc_quality(I,J,K,p);
      if(q_opt<q_new) {
        q_opt=q_new;
        opt_I=I; opt_J=J; opt_K=K;
      }
    }
  }

  Set_partitioning(opt_I,opt_J,opt_K,parti);
  if(q_opt.positive_quality()) return q_opt.Give_p_nec();
  else                         return p+1;
}



//  Partitioning
// ---------------------------------
Partitioning::Partitioning(int level) {
  int i,j,k;
  
  size    = Zweipotenz(level+1);

  start = new bool[size*size*size];

  // set false of Partitioning
  for(i=0;i<size;++i) 
    for(j=0;j<size;++j) 
      for(k=0;k<size;++k) {
        start[i+size*j+size*size*k] = false;
      }
};

Partitioning::~Partitioning() {
  delete(start);
}

void Partitioning::Put_proc(int i, int j, int k) {
  if(developer_version) {
    if(i>=size || j>=size || k>=size) {
      cout << " Error in Partitioning::Put_proc at: "
           << i << ", "  << j << ", "  << k << ", "  << size << endl;
    }
  }
  start[i+size*j+size*size*k]=true;
}

bool Partitioning::Is_this_a_proc(Index3D Ind) {
  Index3D I;
  int t;
  t = Ind.Tiefe()-1;

  /*  
  if(developer_version) {
    if((size-1)*2!=Zweipotenz(t))
      cout << " Error in Partitioning::Is_this_a_proc " << endl;
  }
  */

  // periodic
  I = Ind.neighbour_non_periodic(WSDd);
  return
    start[I.I_x().cart_index(t)+
         size*I.I_y().cart_index(t)+
         size*size*I.I_z().cart_index(t)];
}

void Partitioning::Print() {
  int i,j,k;
 
  cout << " Partitioning: " << endl;
  for(i=0;i<size;++i) {
    cout << " \n new z-level: " << endl;
    for(j=0;j<size;++j) {
      cout << endl;
      for(k=0;k<size;++k) {
        cout << start[i+size*j+size*size*k] << ", ";
      }
    }
  }
  cout << endl;
}

---