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


// 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 "../domain/d_exam.h"

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

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

// Include level 5.1
#include "../evpar/evpar.h"


// calculation of computational time and MFlops
void Grid::Add_operations(long_int ops, int level)  {
  operations=operations+ops*number_interior_points[level];
}
void Grid::Start_calc_time()  {
  calc_yn = 1;
  operations=0;
  t0_all=clock();
  t_sum_interior = 0;
  t_sum_nearbb = 0;
}
void Grid::Stop_calc_time()  {
  calc_yn = 0;
  t1_all=clock();
}
void Grid::Start_calc_time_interior()  {
  if(calc_yn)
    t0_interior=clock();
}
void Grid::Stop_calc_time_interior()  {
  if(calc_yn) {
    t1_interior=clock();
    t_sum_interior = t_sum_interior + (t1_interior - t0_interior);
  }
}
void Grid::Start_calc_time_nearbb()  {
  if(calc_yn) {
    t0_nearbb=clock();
  }
}
void Grid::Stop_calc_time_nearbb()  {
  if(calc_yn) {
    t1_nearbb=clock();
    t_sum_nearbb = t_sum_nearbb + (t1_nearbb - t0_nearbb);
  }
}
void Grid::Report_calc_time( const char * what )  {
  /*
  double t_all = ((double)t1_all-t0_all)/CLOCKS_PER_SEC;
  double t_interior = ((double)t_sum_interior)/CLOCKS_PER_SEC;
  double t_nearbb = ((double)t_sum_nearbb)/CLOCKS_PER_SEC;
  calc_yn = 0;

  printf("\n Report for %s : ",what);
  printf("\n Time interior: %g seconds", t_interior);
  printf("\n Time near the boundary: %g seconds", t_nearbb);
  printf("\n Time all: %g seconds", t_all);
  printf("\n Number of operations at interior points: %f " , operations);
  printf("\n Interior effiency: %g Mflops\n", operations/t_interior/1000000 );
  */
}
          

// parallel constructors                                           
Grid::Grid(int n_max, All_Domains* dom, 
           Grid_gen_parameters& gpara, MPI_Comm comm) : 
  Grid_base(n_max, dom, gpara, comm) {
  non_bock_sending = true;
  Construct_storage();
}
Grid::Grid(int n_max, All_Domains* dom, MPI_Comm comm) : 
  Grid_base(n_max, dom, comm) {
  non_bock_sending = true;
  Construct_storage();
}
Grid::Grid(double h_finest_level, All_Domains* dom, 
           Grid_gen_parameters& gpara, MPI_Comm comm) : 
  Grid_base(h_finest_level, dom, gpara, comm) {
  non_bock_sending = true;
  Construct_storage();
}
Grid::Grid(double h_finest_level, All_Domains* dom, MPI_Comm comm) : 
  Grid_base(h_finest_level, dom, comm) {
  non_bock_sending = true;
  Construct_storage();
}


// seriell constructors                                           
Grid::Grid(int n_max, All_Domains* dom, 
           Grid_gen_parameters& gpara) : 
  Grid_base(n_max, dom, gpara, 0) {
  if(parallel_version)
    cout << " Do not use Grid::Grid(int n_max, All_Domains* dom, Grid_gen_parameters& gpara);  in parallel version!! " << endl;
  Construct_storage();
}
Grid::Grid(double h_finest_level, All_Domains* dom, 
           Grid_gen_parameters& gpara) : 
  Grid_base(h_finest_level, dom, gpara, 0) {
  if(parallel_version)
    cout << " Do not use Grid::Grid(int n_max, All_Domains* dom, Grid_gen_parameters& gpara);  in parallel version!! " << endl;
  Construct_storage();
}

Grid::Grid(int n_max, All_Domains* dom) : 
  Grid_base(n_max, dom, 0) {
  if(parallel_version)
    cout << " Do not use Grid::Grid(int n_max, All_Domains* dom);  in parallel version!! " << endl;
  Construct_storage();
}

Grid::Grid(double h_finest_level, All_Domains* dom) : 
  Grid_base(h_finest_level, dom, 0) {
  if(parallel_version)
    cout << " Do not use Grid::Grid(int n_max, All_Domains* dom);  in parallel version!! " << endl;
  Construct_storage();
}

void Grid::Do_non_bock_sending(bool yn) {
  non_bock_sending = yn;
}

void Grid::Construct_storage() {
  int i;
  nearb_ablage = new Nearb_Ablage;
  vector_ablage = new double[27];
  label_on_cell = new bool[28];
  eight_bools   = new bool[8];

  stencils_fine = new double*[8];
  for(i=0;i<8;++i)
    stencils_fine[i] = new double[64];
  weight_27 = new double[27];
  restriction_stencil = new Restriction_stencil_container;

  coarse_grid_functions = new int[216];
  for(i=0;i<8;++i)
    Set_coarse_function(&(coarse_grid_functions[i*27]),(dir_sons)i);
}

void Grid::Print_maximal_interior_angles() {
  double* h;
  int num;
  double max_edge_ang;
  double max_face_ang;
  double ang;

  if(!Is_storage_initialized()) Initialize();

  max_edge_ang = max_face_ang = 0.0;
  h = new double[30];

  Tetraeder_storage* tet;
  iterate_hash2 {
    for(num=0;num<bocell->Give_number_points();++num) {
      if(bocell->edge_point(num)==edge_poi_typ) {
        h[num] = Give_h(bocell->Give_Index().neighbour(bocell->corner(num)),
                        bocell->edge_dir(num));
      }
    }
    
    for(tet=bocell->Give_tets();tet!=NULL;tet=tet->Next()) {
      ang = calc_maximal_face_angle(bocell->coord(tet->N0(),h),
                                    bocell->coord(tet->N1(),h),
                                    bocell->coord(tet->N2(),h),
                                    bocell->coord(tet->N3(),h));
      if(ang>max_face_ang) max_face_ang = ang;
      if(Print_wrong_Tet) { // only for finding a bug
        if(ang > 155) {
          bocell->Print();
          cout << "Ecken des Tet: "
               << tet->N0() << ", "
               << tet->N1() << ", "
               << tet->N2() << ", "
               << tet->N3() << endl;
          cout << "Koordinaten des Tet: \n";
          bocell->coord(tet->N0(),h).Print(); cout << endl;
          bocell->coord(tet->N1(),h).Print(); cout << endl;
          bocell->coord(tet->N2(),h).Print(); cout << endl;
          bocell->coord(tet->N3(),h).Print(); cout << endl;
        }
      }
      ang = calc_maximal_edge_angle(bocell->coord(tet->N0(),h),
                                    bocell->coord(tet->N1(),h),
                                    bocell->coord(tet->N2(),h),
                                    bocell->coord(tet->N3(),h));
      if(Print_wrong_Tet) { // only for finding a bug
        if(ang > 136) {
          bocell->Print();
          cout << "Ecken des Tet: "
               << tet->N0() << ", "
               << tet->N1() << ", "
               << tet->N2() << ", "
               << tet->N3() << endl;
          cout << "Koordinaten des Tet: \n";
          bocell->coord(tet->N0(),h).Print(); cout << endl;
          bocell->coord(tet->N1(),h).Print(); cout << endl;
          bocell->coord(tet->N2(),h).Print(); cout << endl;
          bocell->coord(tet->N3(),h).Print(); cout << endl;
        }
      }
      if(ang>max_edge_ang) max_edge_ang = ang;
    }
  }
  
  if(max_edge_ang > 136 || max_face_ang > 155) {
    cout << "\n Error: Too large interior angle!" <<
      " Tell Christoph Pflaum this error! ";
  }
  
  cout << "\n Maximal interior face angle: " << max_face_ang;
  cout << "\n Maximal interior edge angle: " << max_edge_ang << endl;
  delete h;
}

double Grid::Total_number_interior_points(int lev) {
  double sum, total_sum;
  sum = (double)number_interior_points[lev];
  
  if(parallel_version)
    MPI_Allreduce(&sum,&total_sum,1,MPI_DOUBLE,MPI_SUM,comm);
  else total_sum = sum;
  return total_sum;
};

double Grid::Total_number_exteri_points(int lev) {
 double sum, total_sum;
  sum = (double)number_exteri_points[lev];
  
  if(parallel_version)
    MPI_Allreduce(&sum,&total_sum,1,MPI_DOUBLE,MPI_SUM,comm);
  else total_sum = sum;
  return total_sum;
};

double Grid::Total_number_nearb_points(int lev) {
 double sum, total_sum;
  sum = (double)number_nearb_points[lev];
  
  if(parallel_version)
    MPI_Allreduce(&sum,&total_sum,1,MPI_DOUBLE,MPI_SUM,comm);
  else total_sum = sum;
  return total_sum;
};

double Grid::Total_number_bo2p_points(int lev) {
 double sum, total_sum;
  sum = (double)number_bo2p_points[lev];
  
  if(parallel_version)
    MPI_Allreduce(&sum,&total_sum,1,MPI_DOUBLE,MPI_SUM,comm);
  else total_sum = sum;
  return total_sum;
};

double Grid::Total_number_cellpoi_points(int lev) {
 double sum, total_sum;
  sum = (double)number_cellpoi_points[lev];
  
  if(parallel_version)
    MPI_Allreduce(&sum,&total_sum,1,MPI_DOUBLE,MPI_SUM,comm);
  else total_sum = sum;
  return total_sum;
};

void Grid::Print_poi_Info(int lev) {
  double sumr, sumin;
  sumr = 0;
  sumin = 0;
  cout << "\n on level: " << lev << ": " << endl;
  cout << "------------------------------------ \n";
  // Summe Innen
  sumin = sumin + (double)Total_number_interior_points(lev);
  cout << "number of interior points: " 
       << Total_number_interior_points(lev) << endl;
  // Summe Randnah
  sumr = sumr + number_nearb_points[lev];
  cout << "number of points near the boundary: " 
       << (double)number_nearb_points[lev] << endl;
  // Summe Randpunkte
  sumr = sumr + number_bo2p_points[lev];
  cout << "number of boundary points: " 
       << (double)number_bo2p_points[lev] << endl;
  // Summe Randzellen
  sumr = sumr + (double)number_cellpoi_points[lev];
  cout << "number of points in boundary cells: ";
  cout << (double)number_cellpoi_points[lev] << endl;
  cout << " \n ratio of points interior / near boundary: " 
       << (double)sumin / (double)sumr << endl;
  cout << "  sum of all points: " << sumin + sumr << endl;
}

void Grid::Print_poi_Sum_Info() {
  int i;
  double sumi;
  double sumr, sumin;
  cout << "\n sum over all levels: " << endl;
  cout << "-------------------------- \n";
  // Summe Innen
  sumi=0;
  for(i=0;i<Max_number_levels;++i) {
    sumi = Total_number_interior_points(i) + sumi;
  }
  cout << "number of interior points: " << sumi << endl;
  sumin = sumi;
  // Summe Randnah
  sumi=0;
  for(i=0;i<Max_number_levels;++i) {
    sumi = (double)number_nearb_points[i] + sumi;
  }
  cout << "number of points near the boundary: " << sumi << endl;
  sumr = sumi;
  // Summe Innen
  sumi=0;
  for(i=0;i<Max_number_levels;++i) {
    sumi = (double)number_bo2p_points[i] + sumi;
  }
  cout << "number of boundary points: " <<  sumi << endl;
  sumr = sumr + sumi;
  // Summe Innen
  sumi=0;
  for(i=0;i<Max_number_levels;++i) {
    sumi = (double)number_cellpoi_points[i] + sumi;
  }
  cout << "number of points in boundary cells: " <<  sumi << endl;
  sumr = sumr + sumi;
  cout << " \n ratio of points interior / near boundary: "
       << (double)sumin / (double)sumr << endl;
  cout << "  sum of all points: " << sumin + sumr << endl;
}


void Grid::Print_poi_Info_all(int lev) {
  int i;
  double sp;

  if(my_rank==0) {
    cout << "\n information until level: " << endl;
    cout << "------------------------------------ \n";
  }
  if(lev>=0 && lev<Max_number_levels) {
    // Innen
    for(i=1;i<=lev;++i) {
      sp = Total_number_interior_points(i);
      if(my_rank==0) {
        cout << "number of interior points on level " << i << " is: ";
        cout << sp << endl;
      }
    }
    // Randnah
    for(i=1;i<=lev;++i) {
      sp = Total_number_nearb_points(i);
      if(my_rank==0) {
        cout << "number of points near the boundary on level " << i << " is: ";
        cout << sp << endl;
      }
    }
    // Randpunkte
    for(i=1;i<=lev;++i) {
      sp = Total_number_bo2p_points(i);
      if(my_rank==0) {
        cout << "number of boundary points on level " << i << " is: ";
        cout << sp << endl;
      }
    }
    // Zellpunkte  
    for(i=1;i<=lev;++i) {
      sp = Total_number_cellpoi_points(i);
      if(my_rank==0) {
        cout << "number of points in boundary cells on level " << i << " is: ";
        cout << sp << endl;
      }
    }
    // Exterior mg points 
    for(i=1;i<=lev;++i) {
      sp = Total_number_exteri_points(i);
      if(my_rank==0) {
        cout << "number of exterior points for mg on level " << i << " is: ";
        cout << sp << endl;
      }
    }
  }
}

void Grid::Print_poi_Info() {
  double sumin, sumr, sp;

  if(my_rank==0) cout << "\n information about finest level: " << endl;
  if(my_rank==0) cout << "------------------------------------ \n";
  // Innen
  if(my_rank==0) cout << "number of interior points is:       ";
  sp = Total_number_interior_points(max_level);
  if(my_rank==0) cout << sp << endl;
  sumin = sp;
  // Randnah
  if(my_rank==0) cout << "number of points near the boundary: ";
  sp = Total_number_nearb_points(max_level);
  if(my_rank==0) cout << sp << endl;
  sumr = sp;
  // Randpunkte
  if(my_rank==0) cout << "number of boundary points:          ";
  sp = Total_number_bo2p_points(max_level);
  if(my_rank==0) cout << sp << endl;
  sumr += sp;
  // Zellpunkte  
  if(my_rank==0) cout << "number of points in boundary cells: ";
  sp = Total_number_cellpoi_points(max_level);
  if(my_rank==0) cout << sp << endl;
  sumr += sp;
  // Sum:
  if(my_rank==0) cout << "sum of all points:                  " 
                      << sumin + sumr << endl;
  if(my_rank==0) cout << "ratio of points interior / near boundary: "
                      << (double)sumin / (double)sumr << endl;

}

void Grid::Refine(Index3D ind)  {
  Add_point(ind);
  
  cout << " \n Please do not use Grid::Refine! " << endl;

  // 4.Schritt: 
  Fullfill_B1B2();
}


void Grid::Initialize() {
  Pointtype typ_point;
  int i, j;
  Index3D I, Inext;
  int lev;

  if(Is_storage_initialized() == false) {
    if(I_am_active()) {
      nearb_ablage->Initialize_leer(Give_max_num_var());
      
      Initialize_variable();
    }
    // Teil IV: Berechnung Listen
    //-----------------------------
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << " IV. calculate lists; " << endl;
    } 
#ifdef USE_PBE
    pbe_start(23,"IV_calc_lists");
#endif
    
    // IV.1.Schritt: Initialisierung der Startpunkte
    start_P_interior = new P_interior**[anzahl_colors];
    for(j=0;j<anzahl_colors;++j) {
      start_P_interior[j] = new P_interior*[Max_number_levels];
      for(i=0;i < Max_number_levels;++i) start_P_interior[j][i]=NULL;
    }
    start_P_nearb = new P_nearb**[anzahl_colors];
    for(j=0;j<anzahl_colors;++j) {
      start_P_nearb[j] = new P_nearb*[Max_number_levels];
      for(i=0;i < Max_number_levels;++i) start_P_nearb[j][i]=NULL;
    }
    start_P_exteri = new P_nearb**[anzahl_colors];
    for(j=0;j<anzahl_colors;++j) {
      start_P_exteri[j] = new P_nearb*[Max_number_levels];
      for(i=0;i < Max_number_levels;++i) start_P_exteri[j][i]=NULL;
    }
    start_P_Bo2p = new P_Bo2p**[anzahl_colors];
    for(j=0;j<anzahl_colors;++j) {
      start_P_Bo2p[j] = new P_Bo2p*[Max_number_levels];
      for(i=0;i < Max_number_levels;++i) start_P_Bo2p[j][i]=NULL;
    }
    start_P_cellpoi = new P_cellpoi**[1];
    for(j=0;j<1;++j) {
      start_P_cellpoi[j] = new P_cellpoi*[Max_number_levels];
      for(i=0;i < Max_number_levels;++i) start_P_cellpoi[j][i]=NULL;
    }

    // for cell variables
    start_P_interior_cell = new P_interior_cell*[Max_number_levels];
    for(i=0;i < Max_number_levels;++i) start_P_interior_cell[i]=NULL;

    start_P_boundary_tet = new P_boundary_tet*[Max_number_levels];
    for(i=0;i < Max_number_levels;++i) start_P_boundary_tet[i]=NULL;

    // ^^^^^^^^ parallel ^^^^^^^^
    start_P_parallel = new P_parallel*[Max_number_levels];
    for(i=0;i < Max_number_levels;++i) start_P_parallel[i]=NULL;
    
    for(i=0;i < Max_number_levels;++i) number_interior_points[i] = 0;
    for(i=0;i < Max_number_levels;++i) number_nearb_points[i] = 0;
    for(i=0;i < Max_number_levels;++i) number_exteri_points[i] = 0;
    for(i=0;i < Max_number_levels;++i) number_cellpoi_points[i] = 0;
    for(i=0;i < Max_number_levels;++i) number_bo2p_points[i] = 0;

    for(i=0;i < Max_number_levels;++i) number_interior_cell_points[i] = 0;
    for(i=0;i < Max_number_levels;++i) number_boundary_cell_tet[i] = 0;
    // ^^^^^^^^ parallel ^^^^^^^^
    for(i=0;i < Max_number_levels;++i) number_parallel_points[i] = 0;

#ifdef USE_PBE
  pbe_stop(23);
#endif
    
    // IV.2.Schritt: innere Punkte und Randnahe Punkte 
    // und Datenspeicher bereitstellen, hashtable 4 !!
    // fuer hashtable 4:
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << "   IV. 2. interior points and points near boundary; " << endl;
    } 
#ifdef USE_PBE
    pbe_start(24,"IV-2._interior_and_near_boundary_points");
#endif

    Initialize_hash4(Factor_hashtable4_lenght*hashtable1_occ+7);

    // Beachte: wenn ein Prozessor fuer einen Punkt verantwortlich ist
    // (responsible)), dann ist er dies auch auf allen groeberen
    // Prozessoren, auf Grund der Konstruktion der groben Zuteilung
    // der Prozessoren

    if(I_am_active()) {    
      iterate_hash1 {
        typ_point = point1->typ;

        // cout << " point: " << typ_point << endl;

        if(typ_point==hierarchical || typ_point==ausnahme)
          Add_list_interior(point1->ind);
        if(typ_point==nearb)
          Add_list_nearb(point1->ind);
        // ^^^^^^^^ parallel ^^^^^^^^
        if(typ_point==parallel_p)
          Add_list_parallel(point1->ind);
        // ^^^^^^^^ exterior ^^^^^^^^
        if(typ_point==multigrid) {
            Add_list_exteri(point1->ind);
        }
      }
      // IV.3.1.Schritt: Randpunkte 2.Art; 
      //                nur Punkte des eigenen Prozesses
      iterate_hash3 {
        if(I_am_responsible_for(bo2point->Give_boindex())) { // ^^ parallel ^^
          Add_list_bo2p(bo2point);
        }
      }
      // IV.3.2.Schritt: Punkte der Randzellfreiheitsgrade;
      //                nur Punkte des eigenen Prozesses
      iterate_hash2 {
        if(I_am_responsible_for(bocell->Give_Index())) { // ^^ parallel ^^
          //  if(true) {
          if(bocell->Exists_bocellpoint()) {
            Add_list_cellpoi(bocell);
          }
        }
      }
      // IV.4.Schritt: Punkte der Zellen-Variablen und Initialisierung
      if(Give_max_num_fine_all_cells()>0) {
        iterate_hash0 {
          if(point0->Give_Index().Cell_index())
            // if(I_am_responsible_for(point0->Give_Index())) {
            if(point0->typ==int_cell || point0->typ==bo_cell) {
              Add_list_interior_cell(point0);
            }
        }
        // -> boundary tets: see Grid_base::Calc_boundary_2points
        start_P_boundary_tet[Max_level()]       = auxiliary_P;
        ++number_boundary_cell_tet[Max_level()] = auxiliary_number;
        /*
        Tetraeder_storage* tet;
        iterate_hash2 {
          for(tet=bocell->Give_tets();tet!=NULL;tet=tet->Next()) {
            Add_list_boundary_tet(bocell,tet);
          }
        }
        */
      }
#ifdef USE_PBE
  pbe_stop(24);
#endif
      // IV.5.Schritt: Zeiger auf double auf Randzellen uebertragen
      if(print_status_of_grid_generator) if(my_rank==0) {
        cout << "   IV. 5. pointer in boundary cells; " << endl;
      } 
#ifdef USE_PBE
    pbe_start(25,"IV-5_pointer_in_boundary_cells");
#endif

      Set_pointer_on_bocell();
      
      // IV.6.Schritt: Preparation for Prolongation
      if(print_status_of_grid_generator) if(my_rank==0) {
        cout << "   IV. 6. restriction weights; " << endl;
      } 
    }   
#ifdef USE_PBE
  pbe_stop(25);
#endif
    // IV.7.Schritt: Caculate coarsest level
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << "   IV. 7. coarsest level; " << endl;
    }

#ifdef USE_PBE
    pbe_start(26,"IV-7_coarsest_level");
#endif
    int min_level_local;
    min_level_local = my_coarsest_level;

    if(I_am_active()) {
      for(i=my_coarsest_level;i<max_level;++i) {
        if(number_nearb_points[i]+
           number_interior_points[i]+
           number_exteri_points[i] == 0)
          min_level_local = i+1;
      }
    }
    else {
      min_level_local=max_level;
    }
    MPI_Allreduce(&min_level_local,&min_level,1,MPI_INT,MPI_MIN,comm);
    

    /*
    if(min_level < n_parallel) {
      if(my_rank==0) cout << " min_level wir auf n_parallel gesetzt! \n " 
                          << " von: " << min_level 
                          << " auf: " << n_parallel
                          << endl;
      min_level = n_parallel;
    }
    */
#ifdef USE_PBE
  pbe_stop(26);
#endif

    if(I_am_active()) {   
      
      // IV.8.Schritt: Initialize lists
      if(print_status_of_grid_generator) if(my_rank==0) {
        cout << "   IV. 8. initialize lists; " << endl;
      } 
#ifdef USE_PBE
    pbe_start(27,"IV-8_init_lists");
#endif
      int color;
      if(initialize_interior) {
        P_interior* iter_i;
        //      for(lev=Max_level()-1;lev <= Max_level();++lev)
        for(lev=0;lev <= Max_level();++lev)
          for(color=0;color < anzahl_colors;++color)
            for(iter_i=Start_P_interior(lev,color);
                iter_i!=NULL;iter_i=iter_i->Next()) {
              iter_i->initialize(this,lev);
            }
      }
    
      if(initialize_nearb) {
        P_nearb* iter_nb;
        // for(lev=Max_level();lev <= Max_level();++lev)
        for(lev=0;lev <= Max_level();++lev)
          for(color=0;color < anzahl_colors;++color)
            for(iter_nb=Start_P_nearb(lev,color);
                iter_nb!=NULL;iter_nb=iter_nb->Next()) {
              iter_nb->initialize(this,lev);
          }

        // for(lev=Max_level();lev <= Max_level();++lev)
        for(lev=0;lev <= Max_level();++lev)
          for(color=0;color < anzahl_colors;++color)
            for(iter_nb=Start_P_exteri(lev,color);
                iter_nb!=NULL;iter_nb=iter_nb->Next()) {
              iter_nb->initialize(this,lev);
          }
      }
      if(initialize_cellpoi) {
        P_cellpoi* iter_cf;
        for(lev=0;lev <= Max_level();++lev)
            for(iter_cf=Start_P_cellpoi(lev);
                iter_cf!=NULL;iter_cf=iter_cf->Next()) {
              iter_cf->initialize(this,lev);
            }
      }
      if(initialize_bo) {
        P_Bo2p* iter_bo;
        for(lev=0;lev <= Max_level();++lev)
          for(color=0;color < anzahl_colors;++color)
            for(iter_bo=Start_P_Bo2p(lev,color);
                iter_bo!=NULL;iter_bo=iter_bo->Next()) {
              iter_bo->initialize(this,lev);
            }
      }
      if(initialize_pa) {
        P_parallel* iter_pa;
        for(lev=Max_level();lev <= Max_level();++lev)
          //    for(lev=0;lev <= Max_level();++lev)
          for(iter_pa=Start_P_parallel(lev);
              iter_pa!=NULL;iter_pa=iter_pa->Next()) {
            iter_pa->initialize(this,lev);
          }
      }
#ifdef USE_PBE
  pbe_stop(27);
#endif
      // IV.9.Schritt: Initialize MG - storage
      if(print_status_of_grid_generator) if(my_rank==0) {
        cout << "   IV. 9. calc mg; " << endl;
      } 
#ifdef USE_PBE
    pbe_start(28,"IV-9_calc_mg");
#endif
      if(mg_coefficients) Calc_MG_Coefficients();
      
#ifdef USE_PBE
  pbe_stop(28);
#endif
      // V.Schritt: prepare evaluation  ^^^^^^^^ parallel ^^^^^^^^
      // if(my_rank==0) cout << " No Prepare_communication(); " << endl;
      if(print_status_of_grid_generator) if(my_rank==0) {
        cout << "   V.1 prepare evaluation; " << endl;
      } 
#ifdef USE_PBE
    pbe_start(29,"V1_prep_eval.");
#endif
      eval_par_object =
        new Evaluation_Parallelization_object(this,Give_max_num_var());
      // -> parallel
#ifdef USE_PBE
  pbe_stop(29);
#endif
      if(parallel_version) {
        if(print_status_of_grid_generator) if(my_rank==0) {
          cout << "   V.2 prepare parallel; " << endl;
          cout << "   a) ; " << endl;
        }

#ifdef USE_PBE
    pbe_start(30,"V2-1");
#endif

        Prepare_communication();
        if(print_status_of_grid_generator) if(my_rank==0) {
          cout << "   b) ; " << endl;
        }
#ifdef USE_PBE
        pbe_stop(30);
#endif

#ifdef USE_PBE
    pbe_start(31,"V2-2");
#endif

        Prepare_communication_boundary_stencils();
        if(print_status_of_grid_generator) if(my_rank==0) {
          cout << "   c) ; " << endl;
        }
#ifdef USE_PBE
        pbe_stop(31);
#endif

#ifdef USE_PBE
    pbe_start(32,"V2-3");
#endif

        Prepare_communication_interior_stencils();
        if(print_status_of_grid_generator) if(my_rank==0) {
          cout << "   end prepare ; " << endl;
        }
#ifdef USE_PBE
        pbe_stop(32);
#endif
      }
    }
  }
  
  // VI.  Schritt: store finest meshsize
  finest_meshsize = finest_mesh_size();

  // VII. Schritt:
  // cal min_interior_level
  for(i=max_level;i >= min_level;--i) {
    if(Total_number_interior_points(i) +
       Total_number_nearb_points(i)  > 1) {
      min_interior_level = i;
    }
  }

  Storge_initialization_done();
}

// IV.5.Schritt: 
void Grid::Set_pointer_on_bocell() {
  int num_points, num;
  Index3D I;
  iterate_hash2 {
    num_points = bocell->Give_number_points();
    I = bocell->ind;
    for(num=0;num<num_points;++num) {
      if(bocell->edge_point(num)==edge_poi_typ) {
        bocell->Set_variable_pointer(num,
                    Give_variable(I.neighbour(bocell->corner(num)),
                                  bocell->edge_dir(num)));
      }
      if(bocell->edge_point(num)==corner_poi_typ) {
        bocell->Set_variable_pointer(num,
                    Give_variable_slow(I.neighbour(bocell->corner(num)),
                                       I.Tiefe()-1));
      }
    }
    if(bocell->Exists_bocellpoint()) {
      bocell->Set_variable_pointer(num_points,bocell->var);
    }
  }
}


void Grid::Test_init() {
  if(Is_storage_initialized() == false) Initialize();
}



// IV.2.Schritt:
void Grid::Add_list_interior(Index3D I) {
  static Index3D Inext;
  static bool weiter;
  static int t, i;

  t = I.Tiefe();
  if(t<my_coarsest_level)
    t = my_coarsest_level;

  for(;t<=Max_level();++t) {
    weiter = true;     // Suppose: no neighbour boundary cell

    for(i=0;i<8 && weiter;++i) {
      Inext = I.next((dir_sons)i,t+1);
      if(Give_cell_typ(Inext)==bo_cell ||
         Give_cell_typ(Inext)==fine_bo_cell) weiter=false;
    }
    if(weiter == true) {
      Add_list_interior(I,t);    // interior point with respect to this level
    }
    else {
      Add_list_nearb(I,t);       // nearb point with respect to this level
    }
  }
}


void Grid::Add_list_interior(Index3D I, int level) {  
  dir_sons color;
  //cout << " \n interior " ; I.coordinate().Print(); cout << endl;

  if(developer_version)
    if(level >= Max_number_levels) {
      cout << " \n level zu gross bei Add_list_interior " << endl;
    }

  color = I.color(level);
  color = reordering_formula(color);

  start_P_interior[color][level] = 
    new P_interior(I,start_P_interior[color][level]);
  ++number_interior_points[level];

  Add_double(I,level,Give_max_num_var());
}


void Grid::Add_list_nearb(Index3D I) {
  static Index3D Inext;
  static int t;

  t = I.Tiefe();
  if(t<my_coarsest_level)
    t = my_coarsest_level;
  
  for(;t<=Max_level();++t) {
    //    if(I_am_responsible_for(I,t))
    Add_list_nearb(I,t);        // nearb point with respect to this level
  }
}


void Grid::Add_list_nearb(Index3D I, int level) {
  dir_sons color;
  color = I.color(level);
  color = reordering_formula(color);

  start_P_nearb[color][level] = new P_nearb(I,start_P_nearb[color][level]);
  ++number_nearb_points[level];

  Add_double(I,level,Give_max_num_var());
}

void Grid::Add_list_exteri(Index3D I) {
  Index3D Inext;
  int t, t_max;

  t = I.Tiefe();
  if(t<my_coarsest_level)
    t = my_coarsest_level;

  t_max = Give_finest_level(I);

  //  for(;t<Max_level();++t) {
  for(;t<=t_max;++t) {
    Add_list_exteri(I,t);        // exteri point with respect to this level
  }
}


void Grid::Add_list_exteri(Index3D I, int level) {
  dir_sons color;
  color = I.color(level);
  color = reordering_formula(color);

  start_P_exteri[color][level] = new P_nearb(I,start_P_exteri[color][level]);
  ++number_exteri_points[level];

  Add_double(I,level,Give_max_num_var());
}


void Grid::Add_list_bo2p(Bo2Point* bopoint) {
  dir_sons color;
  /* old ordering
  color = bopoint->ind.color(bopoint->level);
  color = reordering_formula(color);
  */
  color = reordering_boundary_formula(bopoint->ind,
                                      bopoint->direction,
                                      my_index);

  start_P_Bo2p[color][bopoint->level] = 
    new P_Bo2p(bopoint->ind,bopoint->direction,
               start_P_Bo2p[color][bopoint->level]);
  ++number_bo2p_points[bopoint->level];
}





void Grid::Add_list_cellpoi(BoCell* bocellf) {
  int lev;
  lev = bocell->ind.Tiefe() - 1;
  start_P_cellpoi[0][lev] = 
    new P_cellpoi(bocell->ind,start_P_cellpoi[0][lev]);
  ++number_cellpoi_points[lev];
}

void Grid::Add_list_interior_cell(Point_hashtable0* cellp) {
  int lev;
  lev = cellp->ind.Tiefe() - 1;

  start_P_interior_cell[lev] = 
    new P_interior_cell(cellp->ind,start_P_interior_cell[lev],
                        cellp->Give_var());
  ++number_interior_cell_points[lev];
}

/*
void Grid::Add_list_boundary_tet(BoCell* cellp, Tetraeder_storage* tet) {
  int lev;
  lev = cellp->ind.Tiefe() - 1;

  start_P_boundary_tet[lev] = 
    new P_boundary_tet(cellp->ind,start_P_boundary_tet[lev],
                        tet->Give_variable());
  ++number_boundary_cell_tet[lev];
}
*/

int Grid::Min_interior_level() {
  // return max_level-1;

  return min_interior_level;
  
  // return max_level;
  // return max_level-1;
  
  //  return min_level;
};  

double Grid::Give_number_interior_grid_points(int lev) {
  if(developer_version) {
     if(lev<0 || lev >max_level)
       cout << "error in Grid::Give_number_interior_gid_points!! " << endl;
  }
  return Total_number_interior_points(lev) +
    Total_number_nearb_points(lev);
};  


void Grid::Set_type_of_update(int number_variable, int level, 
                              type_of_update type) const {
  eval_par_object->Set_type_of_update(number_variable,level,type); 
}

// delete storage

void Grid::Delete_storage() {
  // a) hashtables
  Remove_all_hashtables();

  // b) delete things which were constructed in  Grid::Construct_storage()
  delete nearb_ablage;
  delete vector_ablage;
  delete label_on_cell;
  delete eight_bools;
  
  for(int i=0;i<8;++i)
    delete stencils_fine[i];
  delete stencils_fine;

  delete weight_27 ;
  delete restriction_stencil;
  delete coarse_grid_functions;

  // c) lists
  int color, lev;
  P_interior *iter_i, *iter_i_old;
  for(lev=0;lev <= Max_level();++lev)
    for(color=0;color < anzahl_colors;++color)
      for(iter_i=Start_P_interior(lev,color);
          iter_i!=NULL;) {
        iter_i_old = iter_i;
        iter_i=iter_i->Next();
        delete iter_i_old;
      }

  P_nearb *iter_nb, *iter_nb_old;
  for(lev=0;lev <= Max_level();++lev)
    for(color=0;color < anzahl_colors;++color)
      for(iter_nb=Start_P_nearb(lev,color);
          iter_nb!=NULL;) {
        iter_nb_old = iter_nb;
        iter_nb=iter_nb->Next();
        delete iter_nb_old;
      }

  P_cellpoi *iter_cf, *iter_cf_old;
  for(lev=0;lev <= Max_level();++lev)
    for(iter_cf=Start_P_cellpoi(lev);
        iter_cf!=NULL;) {
      iter_cf_old = iter_cf;
      iter_cf=iter_cf->Next();
        delete iter_cf_old;
    }

  P_Bo2p *iter_bo, *iter_bo_old;
  for(lev=0;lev <= Max_level();++lev)
    for(color=0;color < anzahl_colors;++color)
      for(iter_bo=Start_P_Bo2p(lev,color);
          iter_bo!=NULL;iter_bo=iter_bo->Next()) {
        iter_bo_old = iter_bo;
        iter_bo=iter_bo->Next();
        delete iter_bo_old;
      }

  P_parallel *iter_pa, *iter_pa_old;
  for(lev=Max_level();lev <= Max_level();++lev)
    for(iter_pa=Start_P_parallel(lev);
        iter_pa!=NULL;iter_pa=iter_pa->Next()) {
      iter_pa_old = iter_pa;
      iter_pa=iter_pa->Next();
      delete iter_pa_old;
    }

  // for cell variables
  int i;
  
  for(i=0;i < Max_number_levels;++i) delete start_P_interior_cell[i];
  for(i=0;i < Max_number_levels;++i) delete start_P_boundary_tet[i];
  delete start_P_interior_cell;
  delete start_P_boundary_tet;

  for(i=0;i<anzahl_colors;++i)
      delete start_P_interior[i];
  delete start_P_interior;

  for(i=0;i<anzahl_colors;++i)
      delete start_P_nearb[i];
  delete start_P_nearb;

  for(i=0;i<1;++i) 
      delete start_P_cellpoi[i];
  delete start_P_cellpoi;

 for(i=0;i<anzahl_colors;++i) 
      delete start_P_Bo2p[i];
  delete start_P_Bo2p;

  delete start_P_parallel;

  for(lev=0;lev <= Max_level();++lev)
    for(color=0;color < anzahl_colors;++color)
      for(iter_nb=Start_P_exteri(lev,color);
          iter_nb!=NULL;) {
        iter_nb_old = iter_nb;
        iter_nb=iter_nb->Next();
        delete iter_nb_old;
      }
  for(i=0;i<anzahl_colors;++i)
      delete start_P_exteri[i];
  delete start_P_exteri;
}

---