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

// ------------------------------------------------------------
//
// ggen.cc
//
// ------------------------------------------------------------


//  Grid_base::Grid_base

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

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

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

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

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

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

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


// for not active processes
void Grid_base::Dummy_grid_generation() {
 Initialize_hash0(hashtable0_lenght);
 Initialize_hash1(1);
 Initialize_hash2(1);
 Initialize_hash3(1);

  // I.2. Schritt: Vermeidung von typ d) Flaechen auf feinsten Zellen
  bool go_on;
  Start_for_face_correction_parallel();
  for(go_on=must_edges_be_send();go_on;go_on=must_edges_be_send()) {
  }
  End_for_face_correction_parallel();
}


void Grid_base::Grid_generation(int n_max) {
  int number_for_hashs;
  int l;

  // Teil O: Preparations
  // ---------------------

//  MPI_Barrier(comm);
// cout << "C: " << endl;

  // 0.1.Schritt: Maximal level
  if(n_max<=1)
    cout << " please choose maximal level at least 2 in constructor of Grid! "
         << endl;
  if(n_max<domain->Give_n_uniform())
    cout << " level is too small for the discretization of this domain! "
         << endl;
  
  // 0.2.Schritt: Construction of hash table 0:  
  Initialize_hash0(hashtable0_lenght);

  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "  I Basic grid generation by edges; " << endl;
  } 

#ifdef USE_PBE
  pbe_start(2,"I-Calc_type_of_edges_uniform");
#endif

  // Teil I: Basic grid generation by edges 
  //----------------------------------------
  // I.1. Schritt: Berechne Typ der Kanten
  // verbessern ( nur einmal ueber Kanten gehen!)
  Calc_type_of_edges_uniform();  // studies all edges of cells in finest level

#ifdef USE_PBE
  pbe_stop(2);
#endif

  // I.2. Schritt: Vermeidung von typ d) Flaechen auf feinsten Zellen
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   I.2 edge correction; " << endl;
  } 

#ifdef USE_PBE
  pbe_start(3,"I.2-edges_correction");
#endif

  bool go_on, send_yn;
  Start_for_face_correction_parallel();
  while(Correct_faces_uniform()) {};
  // -> parallel
  if(parallel_version) {
    for(go_on=must_edges_be_send();go_on;go_on=must_edges_be_send()) {
      send_yn = send_edges();
      if(send_yn) {
        while(Correct_faces_uniform()) {};
      }
    }
  }
  End_for_face_correction_parallel();

  // I.3. Schritt: Resize hashtable0
  Resize_hash0(Factor_hashtable0_lenght*hashtable0_occ+7);       

#ifdef USE_PBE
  pbe_stop(3);
#endif


  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "  II Calculation of points and cells; " << endl;
  } 
  // Teil II: Calculation of interior points, cell points and 
  //          points near the boundary
  //----------------------------------------
#ifdef USE_PBE
  pbe_start(4,"II-1_Construction_of_hash_table");
#endif

  // II.1.Schritt: Construction of hash table 1:  
  Initialize_hash1(Factor_hashtable1_lenght*hashtable0_occ+7);       

#ifdef USE_PBE
  pbe_stop(4);
#endif


  // II.2.Schritt: Uebertragung der Eckpunkte der Kanten Hash0 in Hash1
  //               typ: nearb
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   II.2 construct basic points; " << endl;
  } 

#ifdef USE_PBE
  pbe_start(5,"II-2_basic_points");
#endif

  Construct_points_hash1();

  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "        send them; " << endl;
  } 
#ifdef USE_PBE
  pbe_stop(5);
#endif


#ifdef USE_PBE
  pbe_start(6,"II-2_send_basic_points");
#endif

  Send_grid_points_direct_parallel(nearb);

#ifdef USE_PBE
  pbe_stop(6);
#endif


  // II.3.Schritt: add some exterior
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   II.3 fullfill B1B2; " << endl;
  } 

#ifdef USE_PBE
  pbe_start(7,"II-3_fullfull_B1B2");
#endif

  Fullfill_B1B2();                

#ifdef USE_PBE
  pbe_stop(7);
#endif


  // II.4.Schritt: construct cell points of Hash0; look at interior>= points
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   II-4_construct_cell_points; " << endl;
  } 

#ifdef USE_PBE
  pbe_start(8,"II-4_cosntr_cell_points");
#endif

  Construct_cell_points_hash0();

#ifdef USE_PBE
  pbe_stop(8);
#endif


  // II.61.Schritt: Typ der Zellen berechnen auf feinsten Gitter (Part1) 
  //               int_cell, ex_cell, fine_bo_cell
  //               cell type depends on the type of corners
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   II.61 interior cell part 1; " << endl;
  }

#ifdef USE_PBE
  pbe_start(9,"II-61_interior_cell_part_1");
#endif

  Calc_interior_cells_Part1();

#ifdef USE_PBE
  pbe_stop(9);
#endif


  // II.62.Schritt: 
  // ^^^^^^^^ Parallel ^^^^^^^^
  // Send info about neighbour boundary cells and edges of them
  // -> parallel

  if(parallel_version) {
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << "   II.62 sends of boundary cells; " << endl;
    }

#ifdef USE_PBE
    pbe_start(10,"II-62_(P)_sends_of_boundary_cells");
#endif


#ifdef USE_PBE
    pbe_start(33,"II-62_(P)Send_boundary_cells_parallel");
#endif
    Send_boundary_cells_parallel();
#ifdef USE_PBE
  pbe_stop(33);
#endif

#ifdef USE_PBE
    pbe_start(34,"II-62_(P)Send_grid_points_direct_parallel(MG)");
#endif
    Send_grid_points_direct_parallel(multigrid);
#ifdef USE_PBE
  pbe_stop(34);
#endif

#ifdef USE_PBE
  pbe_stop(10);
#endif

  }

  // Test_just_this();

  // II.63.Schritt: Typ der Zellen berechnen bei restlichen Zellen (Part2) 
  //                int_cell, ex_cell, bo_cell fine_bo_cell
  //                vergroebert bis auf ganzes Gebiet bis n_parallel
  //                nicht eigene groebere Zellen werden exterior
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   II.63 interior cell part 2; " << endl;
  }

#ifdef USE_PBE
    pbe_start(11,"II-63_(P)_interior_cells_part_2");
#endif
  Calc_interior_cells_Part2();
#ifdef USE_PBE
  pbe_stop(11);
#endif

  // -> parallel
  if(parallel_version) {
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << "        send them; " << endl;
    }
#ifdef USE_PBE
    pbe_start(12,"II-62_(P)_send_them");
#endif
    // a) for finer cells
    // here l is the level of the cell point -1
    for(l=Max_level();l>=n_parallel-1;--l) 
      Send_cells_parallel(l);


    // b) for coarser cells
    // maximal coarsest level:
    //    cout << " my my_coarsest_level: " << my_coarsest_level 
    //   << " my rank: : "  << my_rank << endl;
    for(l=n_parallel-2;l>=my_coarsest_level;--l) {
    //    for(l=n_parallel-2;l>=1;--l) {
      Calc_interior_cells_Part3(l); 
      Send_cells_parallel(l);
    }
#ifdef USE_PBE
  pbe_stop(12);
#endif
  }

  // II.7.Schritt:
  // Einschub: beseitige Punkte die unnoetig sind
  // ---------
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   II.7 calc multigrid points; " << endl;
  }
  // Teil A. - Berechne Multigrid-Punkte
  //           Dies sind hoechstens alte exterior Punkte, 
  //           die Ecken einer nicht exterior Zelle sind
  //         - restlichen exterior Punkte werden in B. beseitigt
  //         - genauere Feststellung der Multigrid-Punkte
  //           erfolgt in mgcoeff.cc aus Interpolationsoperatoreigenschaft 
#ifdef USE_PBE
    pbe_start(13,"II-7_calc_mg_points");
#endif

  Calc_multigrid_points_part1();

#ifdef USE_PBE
  pbe_stop(13);
#endif

  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   II.7 remove exterior points; " << endl;
  }
#ifdef USE_PBE
    pbe_start(14,"II-7_(P)_remove_exterior_points");
#endif

  // Teil B. Beseitige exterior Punkte
  Remove_exterior_points();

#ifdef USE_PBE
  pbe_stop(14);
#endif
  // II.8.Schritt: Berechnung der Ecken 2.Teil:
  //            von randnah auf randnah oder ausnahme oder hierarchical
  //            oder p_parallel ^^^^^^^^ Parallel ^^^^^^^^
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   II.8 calc interior points; " << endl;
  }
#ifdef USE_PBE
    pbe_start(15,"II-8_calc_interior_points");
#endif
  Calc_interior_points();
#ifdef USE_PBE
  pbe_stop(15);
#endif

  // II.9.Schritt: Send coarse grid points
  // -> parallel
  if(parallel_version) {
    if(print_status_of_grid_generator) if(my_rank==0) {
      cout << "   II.9 send coarse grid points; " << endl;
    }
#ifdef USE_PBE
    pbe_start(16,"II-9_(P)_sends_coarse_grid_points");
#endif

    Send_coarse_grid_points_parallel();
#ifdef USE_PBE
  pbe_stop(16);
#endif
  }
  
  
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "  III Calculation of boundary; " << endl;
  } 
  // Teil III: Berechnung der Randzellen und Randpunkte
  //----------------------------------------------------
#ifdef USE_PBE
    pbe_start(17,"III-1_Calc_of_boundary_(hash_init)");
#endif

  // 1.Schritt: Schritt Hash2, Hash3
  number_for_hashs = Calc_number_boundary_cells_for_hash2() + 7;
  Initialize_hash2(Factor_hashtable2_lenght * number_for_hashs + 7);
  Initialize_hash3(Factor_hashtable3_lenght * number_for_hashs + 7);

#ifdef USE_PBE
  pbe_stop(17);
#endif
  // 2.Schritt: Randzellen berechnen
  if(print_status_of_grid_generator) if(my_rank==0) {
    cout << "   III.2 boundary cells; " << endl;
  }

#ifdef USE_PBE
    pbe_start(18,"III-2_boundary_cells");
#endif

  Calc_boundary_cells();

#ifdef USE_PBE
  pbe_stop(18);
#endif
  // other Schritt's of Teil III: are contained in:
  //               Grid_base::Initialize_variable();
  // see also:     Grid::Initialize();

//Information();
//Info_hashtable0();
//Info_hashtable1();

}


// II.3.Schritt:
//--------------
void Grid_base::Fullfill_B1B2() {
  static int occ_old;
  // Schneller machen !!
  occ_old = 0;
  while(occ_old != hashtable1_occ) {
    occ_old = hashtable1_occ;
    iterate_hash1 {
      Add_B1B2(point1->ind);
    }
  }
}
void Grid_base::Add_B1B2(Index3D I) {
  static int t, i, j, k;
  static Index3D In;

  t = I.Tiefe();
  if(I.ind_x.Tiefe() > 0 && I.ind_z.Tiefe() > 0 && I.ind_y.Tiefe() > 0) {
    // B1 Bedingung: 
    if(I.ind_x.Tiefe() == t && I.ind_y.Tiefe() == t && I.ind_z.Tiefe() == t) {
      for(i=0;i<3;++i) for(j=0;j<3;++j) for(k=0;k<3;++k) {
        In = I.next((Ort1D)i,(Ort1D)j,(Ort1D)k,t);
        if(!Exists_Point(In)) Add_point(In);
      }
    }
    else {
      // B2 Bedingung: 
      if(I.ind_x.Tiefe() == t) {
        if(I.ind_y.Tiefe() == t) {
          // zwei Zellen in z_Richtung 
          for(i=0;i<2;++i) Add_points_of_cell(I.next_TD((dir1D)i,t));
        }
        else {
          if(I.ind_z.Tiefe() == t) {
            // zwei Zellen in y_Richtung 
            for(i=0;i<2;++i) Add_points_of_cell(I.next_NS((dir1D)i,t));
          }
          else {
            // vier Zellen in yz_Richtung 
            for(i=0;i<2;++i) for(j=0;j<2;++j) 
              Add_points_of_cell(I.next_NS((dir1D)i,t).next_TD((dir1D)j,t));
          }
        }
      }
      else {
        if(I.ind_y.Tiefe() == t) {
          if(I.ind_z.Tiefe() == t) {
            // zwei Zellen in x_Richtung 
            for(i=0;i<2;++i) Add_points_of_cell(I.next_EW((dir1D)i,t));
          }
          else {
            // vier Zellen in xz_Richtung 
            for(i=0;i<2;++i) for(j=0;j<2;++j) 
              Add_points_of_cell(I.next_EW((dir1D)i,t).next_TD((dir1D)j,t));
          }
        }
        else {
          // vier Zellen in xy_Richtung 
          for(i=0;i<2;++i) for(j=0;j<2;++j) 
            Add_points_of_cell(I.next_EW((dir1D)i,t).next_NS((dir1D)j,t));
        }
      }
    }
  }
}


// II.4.Schritt: construct cell points of Hash0
//----------------------------------------------
void Grid_base::Construct_cell_points_hash0() {
  int i;
  Index3D I, I_next;
  // part 1: add cells near interior points
  iterate_hash1 {
    I = point1->Give_Index();
    if(point1->typ>=interior && I<<my_index)
      // if(point1->typ>=interior)
      for(i=0;i<8;++i) {
        I_next = I.next((dir_sons)i,max_level+1);
        if(my_index<=I_next) {  // diese Bedingung war auch mal auskommentiert
          Add_cell(I_next);
          Add_points_of_cell(I_next);
          // part 2: add recursive coarse cells
          Recursion_Construct_cell_points_hash0(I_next);
        }
      }
  }
  // part 3: add basic cell
  Add_cell(my_index);
  Recursion_Construct_cell_points_hash0(my_index);
}

void Grid_base::Recursion_Construct_cell_points_hash0(Index3D I) {
  Index3D ind(2,2,2);

  if(I == Index3D(0,0,0)) {

    cout << "Hallo" ;

  }

  // Rekursiv weiter
  if(I!=ind) Recursion_Construct_cell_points_hash0(I.father());
  Add_cell(I);
  Add_points_of_cell(I);
}


// II.2.Schritt: Uebertragung der Eckpunkte der Kanten Hash0 in Hash1
//--------------
void Grid_base::Construct_points_hash1() {
  Index3D ind;
  int i;
  iterate_hash0 {
    ind=point0->Give_Index();
    if(developer_version) {
      if(point0->typ!=edge_interior ||
         !ind.Edge_index())
        cout << " Error in Grid_base::Construct_points_hash1() " << endl;
    }
    for(i=0;i<2;++i) {
      Add_point(ind.corner_of_edge_index((dir1D)i),nearb);
    }
  }
}


// I.1. Schritt: Berechne Typ der Kanten
//--------------
void Grid_base::Calc_type_of_edges_uniform() {
  Index3D ind;
  /* for seriel:
  ind = ind.son(ENTd);
  ind = ind.son(WSDd);
  */

  ind = my_index;

  // for pre_check_box_grid_generation
  A_domain_sp = domain->GiveA();
  B_domain_sp = A_domain_sp + domain->GiveVecH();
  Recursion_Calc_type_of_edges_uniform(ind);
}

void Grid_base::Recursion_Calc_type_of_edges_uniform(Index3D I) {
  int i, j;
  Index3D I_son;
  IndexEdge Iedge;
  Index3D indL, indR;
  D3vector VL, VR;
  Edgetype edge_now;

  if(I.Tiefe()<max_level+1) {
    // I ist nicht Index einer feinsten Zelle    
    for(i=0;i<8;++i) {
      I_son = I.son((dir_sons)i);
      // Rekursiv weiter
      if(pre_check_box_grid_generation) {
        VL = transform_coord(I_son.neighbour_non_periodic(WSDd));
        VR = transform_coord(I_son.neighbour_non_periodic(ENTd));

        if(B_domain_sp.x > VL.x && A_domain_sp.x < VR.x &&
           B_domain_sp.y > VL.y && A_domain_sp.y < VR.y &&
           B_domain_sp.z > VL.z && A_domain_sp.z < VR.z) {

          Recursion_Calc_type_of_edges_uniform(I_son);
        }
      }
      else
        Recursion_Calc_type_of_edges_uniform(I_son);
    }
  }
  else {
    // Untersuche Kanten
    for(j=0;j<12;++j) {
      Iedge = I.edge((Edges_cell)j);
      if(Give_edge_typ(Iedge.I())!=edge_interior) {
        indL = Iedge.corner(Ld);
        indR = Iedge.corner(Rd);
        // Keine Ecke der Kante ist am Rand von [0,1]^3 bzw. my_cell
        if(contained_in_boxes(indL) &&
           contained_in_boxes(indR)) {
          // Kante ist im Inneren von [0,1]^3
          VL = transform_coord(indL);
          VR = transform_coord(indR);
          if(domain->calc_edge(VL,VR)) {
            edge_now = (Edgetype)domain->edge(VL,VR);
            if(edge_now==edge_interior) {
              Add_edge(Iedge.I());
            }
          }
          else {
            cout <<
              "\n Dies ist in Recursion_Calc_type_of_edges noch nicht fertig!" 
                 << endl;
          }
        }
      }
    }
  }
}

// I.2. Schritt: Vermeidung von typ d) Flaechen auf feinsten Zellen
//--------------

// new version
bool Grid_base::Correct_faces_uniform() {
  bool back;

  back = false;

  Index3D I;
  iterate_hash0 {
    I = point0->Give_Index();
    Correct_edge(I);
  }
  return back;
}

bool Grid_base::Correct_edge(Index3D I) {
  bool back;

  back = false;
  if(I << my_index) {
    if(I.I_x().Tiefe() == max_level+1) {
      back = back || Correct_face(IndexSurface(I.next_S(max_level+1),Ddir)); 
      back = back || Correct_face(IndexSurface(I.next_N(max_level+1),Ddir)); 
      back = back || Correct_face(IndexSurface(I.next_T(max_level+1),Sdir)); 
      back = back || Correct_face(IndexSurface(I.next_D(max_level+1),Sdir)); 
    }
    if(I.I_y().Tiefe() == max_level+1) {
      back = back || Correct_face(IndexSurface(I.next_E(max_level+1),Ddir)); 
      back = back || Correct_face(IndexSurface(I.next_W(max_level+1),Ddir)); 
      back = back || Correct_face(IndexSurface(I.next_T(max_level+1),Wdir)); 
      back = back || Correct_face(IndexSurface(I.next_D(max_level+1),Wdir)); 
    }
    if(I.I_z().Tiefe() == max_level+1) {
      back = back || Correct_face(IndexSurface(I.next_S(max_level+1),Wdir)); 
      back = back || Correct_face(IndexSurface(I.next_N(max_level+1),Wdir)); 
      back = back || Correct_face(IndexSurface(I.next_E(max_level+1),Sdir)); 
      back = back || Correct_face(IndexSurface(I.next_W(max_level+1),Sdir)); 
    }
  }
  return back;
}

bool Grid_base::Correct_face(IndexSurface Isur) {
  int i;
  int anz;
  // Anzahl der inneren Kanten zaehlen:
  anz=0;
  for(i=0;i<4;++i) {
    if(Give_edge_typ(Isur.edge((dir_2D)i))==edge_interior) ++anz;
  }
  if(anz==3) {
    // Typ d) Flache: also umformen zu Flaeche mit 4 inneren Kanten:
    for(i=0;i<4;++i) {
      if(Give_edge_typ(Isur.edge((dir_2D)i))!=edge_interior) {
        Add_edge_parallel(Isur.edge((dir_2D)i));
        Correct_edge(Isur.edge((dir_2D)i));
      }
    }
    return true;
  }
  else return false;
}

/* old version
bool Grid_base::Correct_faces_uniform() {
  bool back;
  back = false;

  Index3D I;
  iterate_hash0 {
    I = point0->Give_Index();
    if(I.I_x().Tiefe() == max_level+1) {
      back = back || Correct_face(IndexSurface(I.next_S(max_level+1),Ddir)); 
      back = back || Correct_face(IndexSurface(I.next_N(max_level+1),Ddir)); 
      back = back || Correct_face(IndexSurface(I.next_T(max_level+1),Sdir)); 
      back = back || Correct_face(IndexSurface(I.next_D(max_level+1),Sdir)); 
    }
    if(I.I_y().Tiefe() == max_level+1) {
      back = back || Correct_face(IndexSurface(I.next_E(max_level+1),Ddir)); 
      back = back || Correct_face(IndexSurface(I.next_W(max_level+1),Ddir)); 
      back = back || Correct_face(IndexSurface(I.next_T(max_level+1),Wdir)); 
      back = back || Correct_face(IndexSurface(I.next_D(max_level+1),Wdir)); 
    }
    if(I.I_z().Tiefe() == max_level+1) {
      back = back || Correct_face(IndexSurface(I.next_S(max_level+1),Wdir)); 
      back = back || Correct_face(IndexSurface(I.next_N(max_level+1),Wdir)); 
      back = back || Correct_face(IndexSurface(I.next_E(max_level+1),Sdir)); 
      back = back || Correct_face(IndexSurface(I.next_W(max_level+1),Sdir)); 
    }
  }
  return back;
}

bool Grid_base::Correct_face(IndexSurface Isur) {
  int i;
  int anz;
  // Anzahl der inneren Kanten zaehlen:
  anz=0;
  for(i=0;i<4;++i) {
    if(Give_edge_typ(Isur.edge((dir_2D)i))==edge_interior) ++anz;
  }
  if(anz==3) {
    // Typ d) Flache: also umformen zu Flaeche mit 4 inneren Kanten:
    for(i=0;i<4;++i) {
      if(Give_edge_typ(Isur.edge((dir_2D)i))!=edge_interior)
        Add_edge_parallel(Isur.edge((dir_2D)i));
    }
    return true;
  }
  else return false;
}
*/


// Schritt II.61 nd 62
//---------------------


void Grid_base::Calc_interior_cells_Part1() {
  Index3D ind;
  Celltype ce_type;
  
  ind = my_index;

  ce_type = Recursion_calc_interior_cells_Part1(ind);
  Set_cell_typ(ind,ce_type);
}

void Grid_base::Calc_interior_cells_Part2() {
  Index3D ind(2,2,2);
  Celltype ce_type;

  // calc my cells
  ce_type = Recursion_calc_interior_cells_Part2(ind);

  Set_cell_typ(ind,ce_type);

  // Set not my cells exterior
  Index3D I;
  iterate_hash0 {
    I = point0->Give_Index();
    if(((I < my_index) == false) && (I.Tiefe()<= Max_level()))
      point0->typ=ex_cell;
  }
}

// I ist Index einer Zelle, die nicht im Gebiet ist
// Es wird der Zelltyp mit Mittelpunktindex I berechnet
Celltype Grid_base::Calc_cell_type(Index3D I) {
  static bool labelex;
  static bool labelin;
  static Pointtype poitype;
  static int k;

  //  1. Teil untersuche Ecken:
  labelin=false; //Annahme: Es gibt keine innere Ecke
  labelex=false; //Annahme: Es gibt keine aussere Ecke
  for(k=0;k<8;k++) {
    poitype = Give_type(I.neighbour((dir_sons)k));
    if(poitype!=exterior  && 
       poitype!=multigrid &&
       poitype!=boundary) labelin=true;
   if(poitype==exterior || poitype==multigrid) labelex=true;
  }
  if(labelin==false) return ex_cell;
  if(labelex==true)  return bo_cell;

  //  2. Teil Es ist notwendig zusaetzlich die Kanten zu untersuchen
  //          Zelle ist entweder innere Zelle oder Randzelle
  /*
  for(k=0;k<12;k++) {
    if(Give_edge_typ(I.I_edge((Edges_cell)k))==edge_not_interior)
      return bo_cell;
  }
  */
  return int_cell;
}




Celltype Grid_base::Recursion_calc_interior_cells_Part1(Index3D I) {
  Index3D I_son;
  Celltype ce_type;
  bool not_ex;
  bool not_in;
  int i;
  not_ex=false;                 //Annahme: Es gibt keinen Son nicht ex_cell
  not_in=false;                 //Annahme: Es gibt keinen Son nicht in_cell
  if(I.Tiefe()<max_level+1) {
    // I ist nicht Index einer feinsten Zelle    
    for(i=0;i<8;++i) {
      I_son = I.son((dir_sons)i);
      if(Exists_Cell(I_son)) { // Rekursion
        ce_type = Recursion_calc_interior_cells_Part1(I_son);
        // Setze Zelltyp
        Set_cell_typ(I_son,ce_type);      
        if(ce_type != int_cell) not_in=true;
        if(ce_type != ex_cell)  not_ex=true;
      }
      else not_in=true;
    }
    ce_type = bo_cell;
    if(not_in==false) ce_type = int_cell;
    if(not_ex==false) ce_type = ex_cell;
  }
  else {   // Berechne Typ der Zelle entsprechend der Ecken der Zelle
    ce_type = Calc_cell_type(I);
    if(ce_type==bo_cell) ce_type=fine_bo_cell;
    
    // Setze Zelltyp
    if(ce_type != ex_cell) {
      Set_cell_typ(I,ce_type);
    }
  }
  return ce_type;
}

Celltype Grid_base::Recursion_calc_interior_cells_Part2(Index3D I) {
  Index3D I_son;
  Celltype ce_type;
  bool not_ex;
  bool not_in;
  int i;
  not_ex=false;                 //Annahme: Es gibt keinen Son nicht ex_cell
                                // d.h. Annahme: alle Son ex_cell
  not_in=false;                 //Annahme: Es gibt keinen Son nicht in_cell
                                // d.h. Annahme: alle Son in_cell
  if(I.Tiefe()<max_level+1) {
    // I ist nicht Index einer feinsten Zelle    
    for(i=0;i<8;++i) {
      I_son = I.son((dir_sons)i);
      if(Exists_Cell(I_son)) { // Rekursion
        ce_type = Recursion_calc_interior_cells_Part2(I_son);
        // Setze Zelltyp
        Set_cell_typ(I_son,ce_type);      
        if(ce_type != int_cell) not_in=true;
        if(ce_type != ex_cell)  not_ex=true;
      }
      else not_in=true;
    }
    ce_type = bo_cell;
    if(not_in==false) ce_type = int_cell;
    if(not_ex==false) ce_type = ex_cell;
  }
  else { 
    ce_type = Give_cell_typ(I);
    if(I<my_index == false && ce_type!=fine_bo_cell) { 
      // Berechne Typ der Zelle entsprechend der Ecken der Zelle
      ce_type = Calc_cell_type(I);
    }
  }
  // Setze Zelltyp
  if(ce_type != ex_cell) {
    Set_cell_typ(I,ce_type);
  }
  return ce_type;
}

void Grid_base::Calc_interior_cells_Part3(int level) {
  Index3D ind(2,2,2);
  Recursion_calc_interior_cells_Part3(ind,level);
}



void Grid_base::Recursion_calc_interior_cells_Part3(Index3D I, int level) {
  int i;
  Index3D I_son;
  Celltype ce_type;
  bool not_ex;
  bool not_in;

  if(I.Tiefe()<level+1) {
    // I not fine enough   
    for(i=0;i<8;++i) {
      I_son = I.son((dir_sons)i);
        Recursion_calc_interior_cells_Part3(I_son,level);
    }
  }
  else {
    // I is on the right level
    not_ex=false;                 //Annahme: Es gibt keinen Son nicht ex_cell
                                  // d.h. Annahme: alle Son ex_cell
    not_in=false;                 //Annahme: Es gibt keinen Son nicht in_cell
                                  // d.h. Annahme: alle Son in_cell
    
    for(i=0;i<8;++i) {
      I_son = I.son((dir_sons)i);
      if(Exists_Cell(I_son)) {
        ce_type = Give_cell_typ(I_son);
        if(ce_type != int_cell) not_in=true;
        if(ce_type != ex_cell)  not_ex=true;
      }
      else not_in=true;
    }
    ce_type = bo_cell;
    if(not_in==false) ce_type = int_cell;
    if(not_ex==false) ce_type = ex_cell;
    // Setze Zelltyp
    if(ce_type != ex_cell) {
      Add_cell(I);
      /*
      cout << " my rank: " << my_rank
           << " level: " << I.Tiefe() << endl;
      */
      Set_cell_typ(I,ce_type);
    }
  }
}


// Schritt II.8
//--------------
void Grid_base::Calc_interior_points() {
  Pointtype point_typ;
  Index3D I, Inext;
  int i;
  bool weiter;

  iterate_hash1 {
    point_typ = point1->typ;
    if(point_typ != exterior && point_typ != multigrid) {
      I = point1->ind;
      if(I<my_index) {//  ^^^^^^^^ parallel ^^^^^^^^^^^^^
        weiter = true;     // Suppose: no neighbour boundary cell
        for(i=0;i<8 && weiter;++i) {
          Inext = I.next((dir_sons)i,Max_level()+1);
          if(my_index<=Inext) { //  ^^^^^^^^ parallel ^^^^^^^^^^^^^
            if(Give_cell_typ(Inext)!=int_cell) weiter=false;
          }
        }
        if(weiter == true)
          point1->typ = hierarchical;
      }
      else 
        point1->typ = parallel_p;
    }
    /* adaptive version, later!!
    point_typ = point1->typ;
    if(point_typ != exterior && point_typ != multigrid) {
      // gut??
      //      point1->typ = nearb;
      I = point1->ind;
      weiter = true;  // Annahme: es gibt einen Punkt in der Naehe
      for(t=I.Tiefe();weiter;++t) {
        no_27_S   = false;   // Annahme: es gibt 27-Stern
        weiter    = false;   // Annahme: nicht feinere Sterne untersuchen
        for(i=0;i<3;++i)     // fehlt: besser machen mit stop
          for(j=0;j<3;++j)
            for(k=0;k<3;++k) if(!((i==1) && (j==1) && (k==1)))  {
              Inext = I.next((Ort1D)i,(Ort1D)j,(Ort1D)k,t);
              if(Exists_Point(Inext)) {
                weiter    = true;          // Untersuche feinere Sterne
                point_typ = Give_type(Inext);
                // hier fehlt Test der Zellen!!!!!
                if(point_typ == exterior || point_typ == multigrid)
                  no_27_S=true;               // kein 27-Stern vorhanden
              }
              else {
                no_27_S   = true;             // stop kein 27-Stern vorhanden
              }
            }
        if(no_27_S==false) {
          weiter = false;                  // 27-Stern gefunden

          point1->level = t;
          if(t==I.Tiefe()) {
            point1->typ = hierarchical;
          }
          else {
            point1->typ = ausnahme;
          }
          //      point1->typ = interior;
        }
      }
      if(no_27_S==true) {          // kein 27-Stern vorhanden
        point1->level = t-2;
        //  cout << "\n no: " << t; I.coordinate().Print();
      }
    }
    */
  }
}


// Schritt III.2
//---------------
void Grid_base::Calc_boundary_cells() {
  iterate_hash0 {
    if(point0->isCell()) {
      if(point0->typ==fine_bo_cell)
        Add_bocell(point0->ind);
    }
  }
}


// Teil II.1: 
//----------
int Grid_base::Calc_number_boundary_cells_for_hash2() {
  int number;
  number=0;
  iterate_hash0 {
    if(point0->isCell())
      if(point0->typ==fine_bo_cell) ++number;
  }
  return number;
}



// Schritt II.7.B
//---------------
void Grid_base::Calc_multigrid_points_part1() {
  int i;
  Index3D I, Icell;
  Pointtype typI;

  iterate_hash0 {
    Icell = point0->Give_Index();
    if(point0->isCell()) {
      if((point0->typ > ex_cell) && ((Icell.Tiefe()-1) >= my_coarsest_level)) {
        for(i=0;i<8;++i) {
          // study corners of cell
          I = Icell.neighbour((dir_sons)i);
          Add_point(I);
          typI = Give_type(I);
          if(typI==exterior) {
            Set_point_typ(I,multigrid);
            /*
            t = point0->ind.Tiefe()-1; 
            if(t<max_level) // not points on finest level
              Put_finest_level_minimal(I,t);
            */
          }
          else if(typI==multigrid) {
            /*
            t = point0->ind.Tiefe()-1;
            if(t<max_level) // not points on finest level
              Put_finest_level_minimal(I,t);
            */
          }
        }
      }
    }
  }
}



//-----------------------------------------------------
// some old member functions

void Grid_base::Add_Uniform_Grid(int Max_tiefe) {
  Index3D ind(2,2,2);

  if(Max_tiefe < 1) Max_tiefe = 1;
  Recursion_add_Uniform_Grid(ind,Max_tiefe);
}

void Grid_base::Recursion_add_Uniform_Grid(Index3D I,int Max_tiefe) {
  int i;

  // lassen wir zur Zeit weg:
  //  V = transform_coord(I);
  //  if(domain->point_in_domain(V)) {
  Add_point(I);
  //  }

  if(I.Tiefe() <  Max_tiefe) 
    for(i=0;i<8;++i)
      Recursion_add_Uniform_Grid(I.son((dir_sons)i),Max_tiefe);
}

---