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

// ------------------------------------------------------------
//
// gparallel.cc
//
// parallel grid generation Part 1.
//
// ------------------------------------------------------------


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

// 1.  Send neighbour fine boundary cells
// 2.  Send neighbour cells all levels
// 3.  Send neighbour points direct
// 4.  Send neighbour points all levels
// 5.  Send multi grid points
// 6.  member functions for parallel edge correction
// 7.  Lists


// 1. Send neighbour boundary cells


// II.9.Schritt of grid generation:

bool Grid_base::Send_boundary_cell_in_direction(Point_hashtable0* poi,
                                                Index3D next_index) {
  Index3D I, Icell;
  int co;
  
  I = poi->ind;
  if(poi->isCell()) {
    if(poi->typ==fine_bo_cell) {
      // first part: for send on levels
      Icell = I;
      for(co=0;co<8;++co) {
        // study corners of cell
        if(Icell.neighbour((dir_sons)co) << next_index)
          return true;
      }
      // second part: for restriction
      Icell = I.father();
      for(co=0;co<8;++co) {
        // study corners of cell
        if(Icell.neighbour((dir_sons)co) << next_index) {
          return true;
        }
      }
    }
  }
  return false;
}


void Grid_base::Send_boundary_cells_parallel() {
  int k,l,d;
  Index3D next, I, I_edge, Icorner;
  Index3D next_index;

  int number_bo_send, number_bo_rec, num_message, num, sp;
  IndexEdge edge;

  int* send_info;
  int* receive_info;

  MPI_Request req[2];
  MPI_Status status[2];

  int rank_source, rank_destination;

  for(d=0;d<26;++d) {
    // 0. index of receiving next processor 
    next_index = Give_next_on_level(d,my_index);

    // 1. part: basic infos about send and recieve
    if(d<6) {
      rank_source      = give_next_rank_source((dir_3D)d,n_parallel-1);
      rank_destination = give_next_rank_destination((dir_3D)d,n_parallel-1);  
    }
    else if(d<18) {
      rank_source      = give_next_rank_source((Edges_cell)(d-6),n_parallel-1);
      rank_destination = give_next_rank_destination((Edges_cell)(d-6),
                                                    n_parallel-1);  
    }
    else {
      rank_source      = give_next_rank_source((dir_sons)(d-18),n_parallel-1);
      rank_destination = give_next_rank_destination((dir_sons)(d-18),
                                                    n_parallel-1);
    }

    // 2. part: count number of boundary cells to be send
    number_bo_send = 0;
    iterate_hash0 {
      if(Send_boundary_cell_in_direction(point0,next_index)) {
        ++number_bo_send;
      }
    }

    // 3. part: send informations about number_bo_send
    number_bo_rec = 0;
    num_message = 0;
    if(rank_source != -1) {
      MPI_Irecv(&number_bo_rec,1,MPI_INT,rank_source,1,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1) {
      MPI_Isend(&number_bo_send,1,MPI_INT,rank_destination,1,comm,
                &req[num_message]);
      ++num_message;
    }
    MPI_Waitall(num_message,req,status);
    
    // 4. part: build informations about boundary cells + edges
    // How informations are stored:
    // Index of cell          edges (Edges_cell)  corners  >
    // ind_x, ind_y, ind_z,    ...              ,   ...    >   5 int
    send_info = Give_buffer_send_info(5*number_bo_send);
    receive_info = Give_buffer_receive_info(5*number_bo_rec);
    
    num = 0;
    iterate_hash0 {
      if(Send_boundary_cell_in_direction(point0,next_index)) {
        // cell index
        I = point0->ind;
        send_info[num*5]   = I.I_x().get();
        send_info[num*5+1] = I.I_y().get();
        send_info[num*5+2] = I.I_z().get();
        // info edges
        sp = 0;
        for(k=0;k<12;++k) {
          edge = I.edge((Edges_cell)k);
          if(Give_edge_typ(edge.I()) == edge_interior) {        
            // Kante im Gebiet
            sp = (sp << 1) + 1;
          }
          else {
            // Kante nicht im Gebiet
            sp = (sp << 1);
          }
        }
        send_info[num*5+3] = sp;
        // info corners
        sp = 0;
        for(k=0;k<8;++k) {
          Icorner = I.neighbour((dir_sons)k);
          if(Give_type(Icorner) == nearb) {     
            // Ecke im Gebiet
            sp = (sp << 1) + 1;
          }
          else {
            // Ecke nicht im Gebiet
            sp = (sp << 1);
          }
        }
        send_info[num*5+4] = sp;
        num++;
      }
    }
    
    // 5. part: send boundary cells + edge info
    num_message = 0;
    if(rank_source != -1 && number_bo_rec>0) {
      MPI_Irecv(receive_info,5*number_bo_rec,MPI_INT,rank_source,2,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1 && number_bo_send>0) {
      MPI_Isend(send_info,5*number_bo_send,MPI_INT,rank_destination,2,comm,
                &req[num_message]);
      ++num_message;
    }
    // MPI_Barrier(comm);
    MPI_Waitall(num_message,req,status);
    
    
    // 6. part: set boundary cells + edge info
    for(num=0;num<number_bo_rec;++num) {
      I = Index3D(receive_info[num*5],
                  receive_info[num*5+1],
                  receive_info[num*5+2]);
      if(developer_version)
        if(I.Cell_index()==false) 
          cout << "\n Mistake in Send_boundary_cells_parallel() " << endl;
      // Add fine bo cell
      Add_cell(I);
      Set_cell_typ(I,fine_bo_cell);
      // Add neighbours of fine bo cell
      for(k=0;k<8;++k) {
        Add_point(I.neighbour((dir_sons)k));
      }
      // Add edges of fine bo cell (and set type of their corners)
      sp = receive_info[num*5+3];
      for(k=11;k>=0;--k) {
        if((sp&1)==1) { 
          // Kante im Gebiet
          I_edge = I.I_edge((Edges_cell)k);
          Add_edge(I_edge);
          Set_edge_typ(I_edge, edge_interior);
          
          // important for boundary cells in second layer
          for(l=0;l<2;++l) {
            next = I_edge.corner_of_edge_index((dir1D)l);

            /*
if(next == Index3D(9,8,18)) {
  D3vector v;
  v = transform_coord(next);
  cout << " point: " << Give_type(next) 
       << " f level: " << Give_my_finest_parallel_level(next)
       << " my_rank: " << Give_my_rank()
       << " max_level: " << Max_level() 
       << " resposible: " << (next << Give_my_index())
    //  << " vector: x: " << v.x << " y: " << v.y << " z: " << v.z 
       << " x: " << next.I_x().get()
       << " y: " << next.I_y().get()
       << " z: " << next.I_z().get()
       << " x: " << I.I_x().get()
       << " y: " << I.I_y().get()
       << " z: " << I.I_z().get()
       <<  endl;
}
            */


            if(Give_type(next)<interior)
              Set_point_typ(next, parallel_p);
          }
        }
        sp = (sp >> 1);
      }
      // Add corners of fine bo cell
      sp = receive_info[num*5+4];
      for(k=7;k>=0;--k) {
        if((sp&1)==1) { 
          // Ecke im Gebiet
          Icorner = I.neighbour((dir_sons)k);
          Add_point(Icorner);
          if(Give_type(Icorner)<interior)
            Set_point_typ(Icorner, parallel_p);
        }
        sp = (sp >> 1);
      }
    }
  }
}





// 2. Send neighbour cells all levels

// II.9.Schritt of grid generation:


inline bool Grid_base::Send_cell_in_direction(Point_hashtable0* poi,
                                       int t,
                                       Index3D my_lev_index,
                                       Index3D next_index) {
  Index3D I;
  I = poi->ind;
  
  /* old
  if(poi->isCell()) {
    if(poi->Give_Tiefe() == t+1) {
      I = poi->ind;
      if(I < my_lev_index) {
        //  && I < my_lev_index ist wohl unnoetig
        if(poi->typ==int_cell || poi->typ==bo_cell) {
          // send one layer in right and left direction
          */
  int co;
  for(co=0;co<8;++co) {
    // study corners of cell
    if(I.neighbour((dir_sons)co) << next_index) {
      poi->will_be_sent(true);
      return true;
    }
  }
  /*
    }
    }
    }
    }
  */

  poi->will_be_sent(false);
  return false;
}



void Grid_base::Send_cells_parallel(int level) {
  // here level is the level of the cell point -1
  Index3D I, next_index;
  Index3D my_lev_index;

  int number_bo_send, number_bo_rec, num_message, num;

  int d;

  int* send_info;
  int* receive_info;
  
  MPI_Request req[2];
  MPI_Status status[2];
  
  int rank_source, rank_destination;


  // -1. set my_lev_index 
  my_lev_index = Give_my_level_index(level+1);

  // construct extra array for invariant points which are needed.
  int lenght_invar_array = 0;
  iterate_hash0 {
    if(point0->isCell()) {
      if(point0->Give_Tiefe() == level+1) {
        if(point0->ind < my_lev_index) {
          if(point0->typ==int_cell || point0->typ==bo_cell) {
            lenght_invar_array++;
          }
        }
      }
    }
  }
  Point_hashtable0** invariant_array = new Point_hashtable0*[lenght_invar_array];
  int array_i = 0;
  iterate_hash0 {
    if(point0->isCell()) {
      if(point0->Give_Tiefe() == level+1) {
        if(point0->ind < my_lev_index) {
          if(point0->typ==int_cell || point0->typ==bo_cell) {
            invariant_array[array_i] = point0;
            array_i++;
          }
        }
      }
    }
  }


  for(d=0;d<26;++d) { 
    // 0. index of receiving next processor 
    next_index = Give_next_on_level(d,my_lev_index);

    // 1. part: basic infos about send and recieve
    if(d<6) {
      rank_source      = give_next_rank_source((dir_3D)d,level);
      rank_destination = give_next_rank_destination((dir_3D)d,level);  
    }
    else if(d<18) {
      rank_source      = give_next_rank_source((Edges_cell)(d-6),level);
      rank_destination = give_next_rank_destination((Edges_cell)(d-6),level);  
    }
    else {
      rank_source      = give_next_rank_source((dir_sons)(d-18),level);
      rank_destination = give_next_rank_destination((dir_sons)(d-18),level);
    }

    // 2. part: count number of boundary cells to be send
    number_bo_send = 0;


    for(int i = 0;i < lenght_invar_array;i++) {
      if(Send_cell_in_direction(invariant_array[i],level,my_lev_index,next_index))
        ++number_bo_send;
    }


    /* old
    iterate_hash0 {
      if(Send_cell_in_direction(point0,level,my_lev_index,next_index))
        ++number_bo_send;
    }
    */

    // 3. part: send informations about number_bo_send
    number_bo_rec = 0;
    num_message = 0;
    if(rank_source != -1) {
      MPI_Irecv(&number_bo_rec,1,MPI_INT,rank_source,3,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1) {
      MPI_Isend(&number_bo_send,1,MPI_INT,rank_destination,3,comm,
                &req[num_message]);
      ++num_message;
    }
    MPI_Waitall(num_message,req,status);

    // 4. part: build informations about boundary cells + edges
    // How informations are stored:
    // Index of cell          type  >
    // ind_x, ind_y, ind_z,    ...  >   4 int
    send_info = Give_buffer_send_info(4*number_bo_send);
    receive_info = Give_buffer_receive_info(4*number_bo_rec);
    

    num = 0;
    for(int i = 0;i < lenght_invar_array;i++) {
      point0 = invariant_array[i];
      if(point0->doSend()) {
        I = point0->ind;
        // cell index
        send_info[num*4]   = I.I_x().get();
        send_info[num*4+1] = I.I_y().get();
        send_info[num*4+2] = I.I_z().get();
        // type
        send_info[num*4+3] = (int)point0->typ;
        
        ++num;
      }
    }

    /* old
    num = 0;
    iterate_hash0 {
  // if(Send_cell_in_direction(point0,level,my_lev_index,next_index)) {
      if(point0->doSend()) {
        I = point0->ind;
        // cell index
        send_info[num*4]   = I.I_x().get();
        send_info[num*4+1] = I.I_y().get();
        send_info[num*4+2] = I.I_z().get();
        // type
        send_info[num*4+3] = (int)point0->typ;
        
        ++num;
      }
    }
    */
    
    // 5. part: send cells + type
    num_message = 0;
    if(rank_source != -1 && number_bo_rec>0) {
      MPI_Irecv(receive_info,4*number_bo_rec,MPI_INT,rank_source,4,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1 && number_bo_send>0) {
      MPI_Isend(send_info,4*number_bo_send,MPI_INT,rank_destination,4,comm,
                &req[num_message]);
      ++num_message;
    }
    // MPI_Barrier(comm);
    MPI_Waitall(num_message,req,status);
    
    // 6. part: set cells + type
    for(num=0;num<number_bo_rec;++num) {
      I = Index3D(receive_info[num*4],
                  receive_info[num*4+1],
                  receive_info[num*4+2]);
      if(developer_version)
        if(I.Cell_index()==false) 
          cout << "\n Mistake in Send_cells_parallel " << endl;
      // Add cell
      Add_cell(I);
      Set_cell_typ(I,(Celltype)receive_info[num*4+3]);
    }
  }
  delete invariant_array;
}


// 3. Send neighbour points direct


// Teste ob Punkt poi mit Index poi->ind fuer Prozessor mit Index next_index
// ein Ghostpoint ist. Dazu verschiebe poi->ind um "h" in Richtung d
// in der der naechste Prozessor liegt
bool Grid_base::Send_point_direct_in_direction(int d, 
                                               Point_hashtable1* poi,
                                               Index3D next_index) {
 Index3D I;
 Index3D nextS;

  if(poi->typ!=exterior) {

    /*
if(t<Max_level()) {
  if(Index3D(4,4,9)==poi->ind) 
    cout << " my rank: " << my_rank 
         << " typ: " <<poi->typ <<  endl;
}
    */

    I = poi->ind;
    /*
     * Berkeley Feb 2004 and Aenderung 4.6.2003:
     *      if(I < my_index) {
     */
    if(I << my_index) {
      if(d<6) {
        nextS = I.next((dir_3D)d,max_level);
      }
      else if(d<18) {
        nextS = I.next((Edges_cell)(d-6),max_level);
      }
      else
        nextS = I.next((dir_sons)(d-18),max_level);
      
      if(nextS < next_index)
        return true;
    }
  }
  return false;
}

// used only on the finest level
void Grid_base::Send_grid_points_direct_parallel(Pointtype typ) {
  int d;
  Pointtype type_gotten;
  Index3D next_index, I;
  int number_bo_send, number_bo_rec, num_message, num;

  int* send_info;
  int* receive_info;

  MPI_Request req[2];
  MPI_Status status[2];

  int rank_source, rank_destination;
 
  for(d=0;d<26;++d) {
    // 0. index of receiving next processor 
    next_index = Give_next_on_level(d,my_index);

    // 1. part: basic infos about send and recieve
    if(d<6) {
      rank_source      = give_next_rank_source((dir_3D)d,n_parallel-1);
      rank_destination = give_next_rank_destination((dir_3D)d,n_parallel-1);  
    }
    else if(d<18) {
      rank_source      = give_next_rank_source((Edges_cell)(d-6),n_parallel-1);
      rank_destination = give_next_rank_destination((Edges_cell)(d-6),
                                                    n_parallel-1);  
    }
    else {
      rank_source      = give_next_rank_source((dir_sons)(d-18),n_parallel-1);
      rank_destination = give_next_rank_destination((dir_sons)(d-18),
                                                    n_parallel-1);
    }

/*
  cout << " my_rank: " << my_rank
  << " d: " << d
  << " rank_source: " << rank_source
  << " rank_destination: " << rank_destination << endl;
  MPI_Barrier(comm);
*/

    // 2. part: count number of points to be send
    number_bo_send = 0;
    if(rank_destination != -1) {
      iterate_hash1 {
        if(Send_point_direct_in_direction(d, point1,next_index))
          ++number_bo_send;
      }
    }
    
    // 3. part: send informations about number_bo_send
    number_bo_rec = 0;
    num_message = 0;
    if(rank_source != -1) {
      MPI_Irecv(&number_bo_rec,1,MPI_INT,rank_source,160+d,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1) {
      MPI_Isend(&number_bo_send,1,MPI_INT,rank_destination,160+d,comm,
                &req[num_message]);
      ++num_message;
    }
    MPI_Waitall(num_message,req,status);
  
    // 4. part: build informations about points to be send
    // How informations are stored:
    // Index of cell       
    // ind_x, ind_y, ind_z, typ     >   4 int
    send_info = Give_buffer_send_info(4*number_bo_send);
    receive_info = Give_buffer_receive_info(4*number_bo_rec);
    
    num = 0; 
    if(rank_destination != -1) {
      iterate_hash1 {
        if(Send_point_direct_in_direction(d, point1, next_index)) {
          I = point1->ind;
          
          send_info[num*4]   = I.I_x().get();
          send_info[num*4+1] = I.I_y().get();
          send_info[num*4+2] = I.I_z().get();
          send_info[num*4+3] = point1->typ;
          num++;
        }
      }
    }
  
    // 4. part: send points
    num_message = 0;
    if(rank_source != -1 && number_bo_rec>0) {
      MPI_Irecv(receive_info,4*number_bo_rec,MPI_INT,rank_source,180+d,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1 && number_bo_send>0) {
      MPI_Isend(send_info,4*number_bo_send,MPI_INT,rank_destination,180+d,comm,
                &req[num_message]);
      ++num_message;
    }
    // MPI_Barrier(comm);
    MPI_Waitall(num_message,req,status);
    
    // 5. part: set points
    for(num=0;num<number_bo_rec;++num) {
      I = Index3D(receive_info[num*4],
                  receive_info[num*4+1],
                  receive_info[num*4+2]);

      type_gotten = (Pointtype)receive_info[num*4+3];
      // Add point
      Add_point(I); // if a new point was added it gets type exterior
      if(Give_type(I)==exterior) {
        if(type_gotten==typ)
          Set_point_typ(I, typ);
        else
          Set_point_typ(I, parallel_p);
      }
    }
  }
}


// 4.  Send neighbour points all levels


bool Grid_base::Send_coarse_point_direct_in_direction(int d, 
                                                      Point_hashtable1* poi,
                                                      int t, 
                                                      Index3D my_lev_index,
                                                      Index3D next_index) {
  Index3D I;
  Index3D nextS;
  
  if(poi->typ!=exterior) {
    I = poi->ind;
    // here << for multigrid points
    if(poi->Give_Tiefe() <= t) {
      if(I << my_lev_index) {
        if(I << next_index)
          return true;
        
        if(d<6) {
          nextS = I.next((dir_3D)d,t);
        }
        else if(d<18) {
          nextS = I.next((Edges_cell)(d-6),t);
        }
        else
          nextS = I.next((dir_sons)(d-18),t);
        
        // here << for multigrid points
        if(nextS << next_index)
          return true;
      }
    }
  }
  return false;
}


void Grid_base::Send_coarse_grid_points_parallel() {
  int l; 

  // set my_finest_parallel_level for finest points: Max_level()
  iterate_hash1 {
    point1->my_finest_parallel_level = Max_level();
  }

  // here set my_finest_parallel_level for new points: l 
  for(l=Max_level()-1;l>=my_coarsest_level;--l) {
    Send_coarse_grid_points_parallel(l);
  }
}

void Grid_base::Send_coarse_grid_points_parallel(int level) {
  int d;
  Pointtype type_gotten;
  Pointtype type_old;

  Index3D I, next_index;
  int number_bo_send, number_bo_rec, num_message, num;
  Index3D my_lev_index;

  int* send_info;
  int* receive_info;

  MPI_Request req[2];
  MPI_Status status[2];

  int rank_source, rank_destination;

  // -1. set my_lev_index 
  my_lev_index = Give_my_level_index(level+1);

  for(d=0;d<26;++d) { 
    // 0. index of receiving next processor 
    next_index = Give_next_on_level(d,my_lev_index);

    // 1. part: basic infos about send and recieve
    if(d<6) {
      rank_source      = give_next_rank_source((dir_3D)d,level);
      rank_destination = give_next_rank_destination((dir_3D)d,level);  
    }
    else if(d<18) {
      rank_source      = give_next_rank_source((Edges_cell)(d-6),level);
      rank_destination = give_next_rank_destination((Edges_cell)(d-6),level);  
    }
    else {
      rank_source      = give_next_rank_source((dir_sons)(d-18),level);
      rank_destination = give_next_rank_destination((dir_sons)(d-18),level);
    }
    // 2. part: count number of points to be send
    number_bo_send = 0;
    if(rank_destination != -1) {
      iterate_hash1 {
        if(Send_coarse_point_direct_in_direction(d, point1,level,
                                                 my_lev_index,next_index))
          ++number_bo_send;
      }
    }
    
    // 3. part: send informations about number_bo_send
    number_bo_rec = 0;
    num_message = 0;
    if(rank_source != -1) {
      MPI_Irecv(&number_bo_rec,1,MPI_INT,rank_source,100+d,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1) {
      MPI_Isend(&number_bo_send,1,MPI_INT,rank_destination,100+d,comm,
                &req[num_message]);
      ++num_message;
    }
    MPI_Waitall(num_message,req,status);
  
    // 4. part: build informations about points to be send
    // How informations are stored:
    // Index of cell       
    // ind_x, ind_y, ind_z, typ     >   4 int
    send_info = Give_buffer_send_info(4*number_bo_send);
    receive_info = Give_buffer_receive_info(4*number_bo_rec);
    
    num = 0; 
    if(rank_destination != -1) {
      iterate_hash1 {
        if(Send_coarse_point_direct_in_direction(d, point1,level,
                                                 my_lev_index,next_index)) {
          I = point1->ind;
          
          send_info[num*4]   = I.I_x().get();
          send_info[num*4+1] = I.I_y().get();
          send_info[num*4+2] = I.I_z().get();
          send_info[num*4+3] = point1->typ;
          num++;
        }
      }
    }
  
    // 4. part: send points
    num_message = 0;
    if(rank_source != -1 && number_bo_rec>0) {
      MPI_Irecv(receive_info,4*number_bo_rec,MPI_INT,rank_source,120+d,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1 && number_bo_send>0) {
      MPI_Isend(send_info,4*number_bo_send,MPI_INT,rank_destination,120+d,comm,
                &req[num_message]);
      ++num_message;
    }
    // MPI_Barrier(comm);
    MPI_Waitall(num_message,req,status);
    
    // 5. part: set points
    for(num=0;num<number_bo_rec;++num) {
      I = Index3D(receive_info[num*4],
                  receive_info[num*4+1],
                  receive_info[num*4+2]);

      type_gotten = (Pointtype)receive_info[num*4+3];
      // Add point
      Add_point(I); // if a new point was added it gets type exterior
      type_old = Give_type(I);
      if(type_old==exterior || type_old==multigrid) {
        if(type_gotten==multigrid)
          Set_point_typ(I,level, multigrid);
        else
          Set_point_typ(I,level, parallel_p);
        /*
        if(developer_version) {
          if((I << my_lev_index)) {
            cout << " error in Send_coarse_grid_points_parallel!" << endl;
          }
        }
        */
      }
    }
  }
}




// 5. Multigrid points


bool Grid_base::Send_multi_grid_point_in_direction(int d, 
                                                   Point_hashtable1* poi,
                                                   int t, 
                                                   Index3D my_lev_index, 
                                                   Index3D next_index) {
  Index3D I;
  Index3D nextS;
  
  /*
if(t<Max_level()) {
  if(Index3D(4,4,9)==poi->ind) 
    cout << " my rank: " << my_rank 
         << " typ: " <<poi->typ
         << " fl: " <<poi->finest_level <<  endl;
}
  */

  if(poi->typ==multigrid) {
    I = poi->ind;
    if(poi->Give_Tiefe() <= t) {
      if(I << my_lev_index) {
        if(I << next_index)
        return true;
        
        if(d<6) {
          nextS = I.next((dir_3D)d,t);
        }
        else if(d<18) {
          nextS = I.next((Edges_cell)(d-6),t);
        }
        else
          nextS = I.next((dir_sons)(d-18),t);
        
        //      if((nextS < my_lev_index)==false)
        if(nextS << next_index)
          return true;
      }
    }
  }

  return false;
}


void Grid_base::Send_multi_grid_points_parallel(int level) {
  int d;
  int finest_level;
  Index3D I, next_index;
  int number_bo_send, number_bo_rec, num_message, num;
  Index3D my_lev_index;

  int* send_info;
  int* receive_info;

  MPI_Request req[2];
  MPI_Status status[2];

  int rank_source, rank_destination;

  // 0. set my_lev_index 
  my_lev_index = Give_my_level_index(level+2);
  //  my_lev_index = Give_my_level_index(level+1);

  for(d=0;d<26;++d) {
    // 0. part: basic infos about send and recieve
    if(d<6) {
      rank_source      = give_next_rank_source((dir_3D)d,level+1);
      rank_destination = give_next_rank_destination((dir_3D)d,level+1);  
    }
    else if(d<18) {
      rank_source      = give_next_rank_source((Edges_cell)(d-6),level+1);
      rank_destination = give_next_rank_destination((Edges_cell)(d-6),level+1);  
    }
    else {
      rank_source      = give_next_rank_source((dir_sons)(d-18),level+1);
      rank_destination = give_next_rank_destination((dir_sons)(d-18),level+1);
    }
    
    /*
    if(d<6) {
      rank_source      = give_next_rank_source((dir_3D)d,level);
      rank_destination = give_next_rank_destination((dir_3D)d,level);  
    }
    else if(d<18) {
      rank_source      = give_next_rank_source((Edges_cell)(d-6),level);
      rank_destination = give_next_rank_destination((Edges_cell)(d-6),level);  
    }
    else {
      rank_source      = give_next_rank_source((dir_sons)(d-18),level);
      rank_destination = give_next_rank_destination((dir_sons)(d-18),level);
    }
    */
    
    // 1. index of receiving next processor and correct rank's
    //next_index = Give_next_on_level(d,my_lev_index);
    next_index = Give_level_index_of(rank_destination,level+1);

    // periodisch neu:
    if(next_index == Index3D(3,3,3))
      rank_destination = -1;

    if(Give_my_coarsest_level()>=level+1)
      rank_source      = -1;


    // 2. part: count number of points to be send
    number_bo_send = 0;
    if(rank_destination != -1) {
      iterate_hash1 {
        if(Send_multi_grid_point_in_direction(d, point1,level,
                                              my_lev_index, next_index))
          ++number_bo_send;
      }
    }
    
    // 3. part: send informations about number_bo_send
    number_bo_rec = 0;
    num_message = 0;
    if(rank_source != -1) {
      MPI_Irecv(&number_bo_rec,1,MPI_INT,rank_source,140+d,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1) {
      MPI_Isend(&number_bo_send,1,MPI_INT,rank_destination,140+d,comm,
                &req[num_message]);
      ++num_message;
    }
    MPI_Waitall(num_message,req,status);
  
    // 4. part: build informations about points to be send
    // How informations are stored:
    // Index of cell       
    // ind_x, ind_y, ind_z, typ     >   4 int
    send_info = Give_buffer_send_info(4*number_bo_send);
    receive_info = Give_buffer_receive_info(4*number_bo_rec);
    
    num = 0; 
    if(rank_destination != -1) {
      iterate_hash1 {
        if(Send_multi_grid_point_in_direction(d, point1,level,
                                              my_lev_index, next_index)) {
          I = point1->ind;
          
          send_info[num*4]   = I.I_x().get();
          send_info[num*4+1] = I.I_y().get();
          send_info[num*4+2] = I.I_z().get();
          send_info[num*4+3] = point1->finest_level;
          num++;
        }
      }
    }
  
    // 4. part: send points
    num_message = 0;
    if(rank_source != -1 && number_bo_rec>0) {
      MPI_Irecv(receive_info,4*number_bo_rec,MPI_INT,rank_source,160+d,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1 && number_bo_send>0) {
      MPI_Isend(send_info,4*number_bo_send,MPI_INT,rank_destination,160+d,comm,
                &req[num_message]);
      ++num_message;
    }
    // MPI_Barrier(comm);
    MPI_Waitall(num_message,req,status);
    
    // 5. part: set points
    for(num=0;num<number_bo_rec;++num) {
      I = Index3D(receive_info[num*4],
                  receive_info[num*4+1],
                  receive_info[num*4+2]);

      finest_level = (Pointtype)receive_info[num*4+3];
      // Add point
      /*
      if(developer_version) {
        if(Exists_Point(I)==false) {
          cout << " Error in Send_multi_grid_points_parallel " 
               << endl;
        }
      }
      */
      Add_point(I); // if a new point was added it gets type exterior
      if(Give_type(I)==exterior)
        Set_point_typ(I,multigrid); 
      Add_multigrid_point(I, finest_level);
    }
  }
}









// 6. member functions for parallel edge correction

void Grid_base::Start_for_face_correction_parallel() {
  int i;
  for(i=0;i<18;++i) {
    edges_to_be_sent[i] = new int[30];
    lenght_edges_to_be_sent[i] = 10;
    number_edges_to_be_sent[i] = 0;
  }
  receive_buffer = new int[30];
  length_receive_buffer = 30;
  direction_and_number_process = new int[2*give_number_of_processes()];
}

void Grid_base::End_for_face_correction_parallel() {
  int i;
  for(i=0;i<18;++i) {
    delete(edges_to_be_sent[i]);
  }
  delete(direction_and_number_process);
  delete(receive_buffer);
}

void Grid_base::Test_array(int i) {
  int k;
  int *sp;
  if(number_edges_to_be_sent[i] >= lenght_edges_to_be_sent[i]) {
    lenght_edges_to_be_sent[i] = lenght_edges_to_be_sent[i] * 2;
    sp =  new int[3*lenght_edges_to_be_sent[i]];
    for(k=0;k<lenght_edges_to_be_sent[i]/2*3;++k) {
      sp[k] = edges_to_be_sent[i][k];
    }
    delete(edges_to_be_sent[i]);
    edges_to_be_sent[i] = sp;
  }
}



void Grid_base::put_for_send_edge(Index3D I) {
  int i;
  if(developer_version) {
    if(!I.Edge_index()) {
      cout << "\n Error in Grid_base::put_for_send_edge!" << endl;
      return;
    }
  }
  for(i=0;i<6;++i) {
    if(my_index.surface((dir_3D)i)>>I) {
      Test_array(i);
      edges_to_be_sent[i][3*number_edges_to_be_sent[i]]   = I.I_x().get();
      edges_to_be_sent[i][3*number_edges_to_be_sent[i]+1] = I.I_y().get();
      edges_to_be_sent[i][3*number_edges_to_be_sent[i]+2] = I.I_z().get();
      ++number_edges_to_be_sent[i];
    }
  }
  for(i=6;i<18;++i) {
    if(my_index.edge((Edges_cell)(i-6))>>I) {
      Test_array(i);
      edges_to_be_sent[i][3*number_edges_to_be_sent[i]]   = I.I_x().get();
      edges_to_be_sent[i][3*number_edges_to_be_sent[i]+1] = I.I_y().get();
      edges_to_be_sent[i][3*number_edges_to_be_sent[i]+2] = I.I_z().get();
      ++number_edges_to_be_sent[i];
    }
  }
}
int *Grid_base::Give_receive_buffer(int l) {
  if(length_receive_buffer<l) {
    delete(receive_buffer);
    receive_buffer = new int[l];
  }
  return receive_buffer;
}

bool Grid_base::must_edges_be_send() {
  int i, I_must_send,global_must_send;

  // 1.step: must I send?
  I_must_send = false;
  //  for(i=0;i<18;++i) {
//  Achtung!!!
  for(i=0;i<6;++i) {
    if(number_edges_to_be_sent[i]>0) I_must_send=true;
  }
 
  // 2.step: send-receive
  MPI_Allreduce(&I_must_send,&global_must_send,1,MPI_INT,MPI_LOR,comm);

  return global_must_send;
}

// used only on the finest level
bool Grid_base::send_edges() {
  int d;
  bool edges_where_added;
  Index3D I;
  int number_bo_send, number_bo_rec, num_message, num;

  int* send_info;
  int* receive_info;

  MPI_Request req[2];
  MPI_Status status[2];

  int rank_source, rank_destination;
 
  edges_where_added = false;
  for(d=0;d<18;++d) {
    // 1. part: basic infos about send and recieve
    if(d<6) {
      rank_source      = give_next_rank_source((dir_3D)d,n_parallel-1);
      rank_destination = give_next_rank_destination((dir_3D)d,n_parallel-1);  
    }
    else {
      rank_source      = give_next_rank_source((Edges_cell)(d-6));
      rank_destination = give_next_rank_destination((Edges_cell)(d-6));  
    }
    // 2.1 part: count number of points to be send
    number_bo_send = number_edges_to_be_sent[d];
  
    // 2.2 part: send informations about number_bo_send
    number_bo_rec = 0;
    num_message = 0;
    if(rank_source != -1) {
      MPI_Irecv(&number_bo_rec,1,MPI_INT,rank_source,11,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1) {
      MPI_Isend(&number_bo_send,1,MPI_INT,rank_destination,11,comm,
                &req[num_message]);
      ++num_message;
    }
    MPI_Waitall(num_message,req,status);

    /*
    cout << " myrank: " << my_rank 
         << " number_bo_rec "  << number_bo_rec 
         << " number_bo_send " << number_bo_send
         << endl;
    */

    if(number_bo_rec>0) edges_where_added = true;  

    // 3. part: build informations about points to be send
    // How informations are stored:
    // Index of edges      
    // ind_x, ind_y, ind_z     >   3 int
    send_info = edges_to_be_sent[d];
    receive_info = Give_receive_buffer(3*number_bo_rec);

    // 4. part: send grid points
    num_message = 0;
    if(rank_source != -1 && number_bo_rec>0) {
      MPI_Irecv(receive_info,3*number_bo_rec,MPI_INT,rank_source,12,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1 && number_bo_send>0) {
      MPI_Isend(send_info,3*number_bo_send,MPI_INT,rank_destination,12,comm,
                &req[num_message]);
      ++num_message;
    }
    MPI_Waitall(num_message,req,status);
  
    // 5. part: set points
    for(num=0;num<number_bo_rec;++num) {
      I = Index3D(receive_info[num*3],
                  receive_info[num*3+1],
                  receive_info[num*3+2]);
      // Add point
      Add_edge(I);
    }

    // 6. part: Set zero
    number_edges_to_be_sent[d]=0;
  }

  //  if(edges_where_added) cout << " Edges send " << endl;

  return edges_where_added;
}

// 7. Lists
void Grid::Add_list_parallel(Index3D I) {
  Index3D Inext;
  int t;
  int my_finest_parallel_level;

  t = I.Tiefe();
  if(t<my_coarsest_level-1)
    t = my_coarsest_level-1;
  // above my_coarsest_level-1 is necessary for prolongation
  my_finest_parallel_level = Give_my_finest_parallel_level(I);
  //  my_finest_parallel_level =Max_level();

  //  for(;t<=Max_level();++t) {
  for(;t<=my_finest_parallel_level;++t) {
    if(developer_version)
      if(t >= Max_number_levels) {
        cout << " \n level zu gross bei Add_list_interior " << endl;
      }
    start_P_parallel[t] = new P_parallel(I,start_P_parallel[t]);
    ++number_parallel_points[t];
    
    Add_double(I,t,Give_max_num_var());
  }
}

---