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

// ------------------------------------------------------------
//
// print_avs.cc
//
// ------------------------------------------------------------


// Aendern fuer parallel:
//      if(point1->typ >= interior || point1->typ==parallel_p) {



// a)  for AVS
// b)  for Gnuplot


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

// a)  for AVS

void Grid_base::Print_ec_avs(IndexBo indexbo_A, IndexBo indexbo_B,
                                   ofstream *Datei, int number) {
  *Datei << number << "  1  line  ";
  *Datei << hashtable3_point(indexbo_A)->nummer  << "  ";
  *Datei << hashtable3_point(indexbo_B)->nummer  << "  ";
  *Datei << "\n"; 
}

// Ausgabe der Tetraeder
int Grid_base::avs_bo_cell(BoCell* bo,ofstream *datei,bool print_or_calc, 
                           int number) {
  Tetraeder_storage* tets;
  int num;
  Index3D I;
  Cell_type_points typ;

#ifdef PERIODIC
  if(bo->Give_Index().neighbour_non_periodic(WSDd).Is_non_periodic()) {
    return number;
  }
#endif

  if(print_or_calc==false) {
    for(tets = bo->Give_tets();tets!=NULL;tets = tets->Next()) {
      number++;
    }
  }
  else {
    I = bo->Give_Index();
    for(tets = bo->Give_tets();tets!=NULL;tets = tets->Next()) {
      ++number; 
      *datei << number << "  1  tet  ";
      // Punkt 0:
      num = tets->N0();
      typ = bo->edge_point(num);
      if(typ==edge_poi_typ) {
        *datei << Give_Nummer(I.neighbour(bo->corner(num)),
                              bo->edge_dir(num)) << "  ";
      }
      if(typ==corner_poi_typ) {
        *datei << Give_Nummer(I.neighbour(bo->corner(num)))
               << "  "; 
      }
      if(typ==cell_poi_typ) {
        *datei << Give_Nummer_cellpoi(I) << "  ";       
      }
      // Punkt 1:
      num = tets->N1();
      typ = bo->edge_point(num);
      if(typ==edge_poi_typ) {
        *datei << Give_Nummer(I.neighbour(bo->corner(num)),
                                    bo->edge_dir(num)) << "  ";
      }
      if(typ==corner_poi_typ) {
        *datei << Give_Nummer(I.neighbour(bo->corner(num)))
               << "  "; 
      }
      if(typ==cell_poi_typ) {
        *datei << Give_Nummer_cellpoi(I) << "  ";       
      }
      // Punkt 2:
      num = tets->N2();
      typ = bo->edge_point(num);
      if(typ==edge_poi_typ) {
        *datei << Give_Nummer(I.neighbour(bo->corner(num)),
                                    bo->edge_dir(num)) << "  ";
      }
      if(typ==corner_poi_typ) {
        *datei << Give_Nummer(I.neighbour(bo->corner(num)))
               << "  "; 
      }
      if(typ==cell_poi_typ) {
        *datei << Give_Nummer_cellpoi(I) << "  ";       
      }
      // Punkt 3:
      num = tets->N3();
      typ = bo->edge_point(num);
      if(typ==edge_poi_typ) {
        *datei << Give_Nummer(I.neighbour(bo->corner(num)),
                                    bo->edge_dir(num)) << "  ";
      }
      if(typ==corner_poi_typ) {
        *datei << Give_Nummer(I.neighbour(bo->corner(num)))
               << "  "; 
      }
      if(typ==cell_poi_typ) {
        *datei << Give_Nummer_cellpoi(I) << "  ";       
      }
      *datei << endl;
    }
  }
  return number;
}

int Grid_base::avs_bo_cell(BoCell* bo,ofstream *datei, 
                           int number,int number_cell_var) {
  Tetraeder_storage* tets;

  for(tets = bo->Give_tets();tets!=NULL;tets = tets->Next()) {
    ++number; 
    *datei << number << "  " << tets->Give_variable()[number_cell_var] << "\n";
  }
  return number;
}


int Grid_base::Recursion_Count_Cells(Index3D I) {
  Index3D Ison;

  int summe, i;
  summe = 0;
  for(i=0;i<8;++i) {
    if(Exists_Point(I.son((dir_sons)i)))
      summe = summe + Recursion_Count_Cells(I.son((dir_sons)i));
    else {
      Ison = I.son((dir_sons)i);
      // periodic
#ifdef PERIODIC
      if(Ison.neighbour_non_periodic(WSDd).Is_non_periodic()
         && Give_cell_typ(Ison) == int_cell) {
#else
      if(Give_cell_typ(Ison) == int_cell) {
#endif
        summe = summe + 6;
      }
    }
  }
  return summe;
}

// just now not used function
void Grid_base::Print_Cell_avs(Index3D I, ofstream *Datei) {
  *Datei << Give_Nummer(I.neighbour_EW(Ld).neighbour_NS(Ld).neighbour_TD(Rd)) << "  ";
  *Datei << Give_Nummer(I.neighbour_EW(Rd).neighbour_NS(Ld).neighbour_TD(Rd)) << "  ";
  *Datei << Give_Nummer(I.neighbour_EW(Rd).neighbour_NS(Rd).neighbour_TD(Rd)) << "  ";
  *Datei << Give_Nummer(I.neighbour_EW(Ld).neighbour_NS(Rd).neighbour_TD(Rd)) << "  ";

  *Datei << Give_Nummer(I.neighbour_EW(Ld).neighbour_NS(Ld).neighbour_TD(Ld)) << "  ";
  *Datei << Give_Nummer(I.neighbour_EW(Rd).neighbour_NS(Ld).neighbour_TD(Ld)) << "  ";
  *Datei << Give_Nummer(I.neighbour_EW(Rd).neighbour_NS(Rd).neighbour_TD(Ld)) << "  ";
  *Datei << Give_Nummer(I.neighbour_EW(Ld).neighbour_NS(Rd).neighbour_TD(Ld)) << "  ";
}


int Grid_base::Recursion_Cells_AVS(Index3D I, ofstream *Datei, int nummer) {
  int i;
  Index3D Ison;

  for(i=0;i<8;++i) {
    if(Exists_Point(I.son((dir_sons)i)))
      nummer = Recursion_Cells_AVS(I.son((dir_sons)i),Datei,nummer);
    else {
      Ison = I.son((dir_sons)i);
      // periodic
#ifdef PERIODIC
      if(Ison.neighbour_non_periodic(WSDd).Is_non_periodic()
         && Give_cell_typ(Ison) == int_cell) {
#else
      if(Give_cell_typ(Ison) == int_cell) {
#endif
        nummer++;

        *Datei << nummer << "  1  tet  ";
        // 0. Tet: WND, WNT, WST, EST 
        *Datei << Give_Nummer(Ison.neighbour(WNTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(WNDd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(WSTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(ESTd))  << "  ";
        *Datei << "\n";
        
        nummer++;
        *Datei << nummer << "  1  tet  ";
        // 1. Tet: EST, WND, WST, ESD  
        *Datei << Give_Nummer(Ison.neighbour(ESTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(WNDd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(WSTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(ESDd))  << "  ";
        *Datei << "\n";
        
        nummer++;
        *Datei << nummer << "  1  tet  ";
        // 2. Tet: WND, WSD, WST, ESD        
        *Datei << Give_Nummer(Ison.neighbour(WNDd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(WSDd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(WSTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(ESDd))  << "  ";
        *Datei << "\n";
        
        nummer++;
        *Datei << nummer << "  1  tet  ";
        // 3. Tet: EST, WND, ESD, END
        *Datei << Give_Nummer(Ison.neighbour(ESTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(WNDd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(ESDd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(ENDd))  << "  ";
        *Datei << "\n";
        
        nummer++;
        *Datei << nummer << "  1  tet  ";
        // 4. Tet: ENT, WNT, EST, END
        *Datei << Give_Nummer(Ison.neighbour(ENTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(WNTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(ESTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(ENDd))  << "  ";
        *Datei << "\n";
        
        nummer++;
        *Datei << nummer << "  1  tet  ";
        // 5. Tet: WNT, WND, EST, END
        *Datei << Give_Nummer(Ison.neighbour(WNTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(WNDd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(ESTd))  << "  ";
        *Datei << Give_Nummer(Ison.neighbour(ENDd))  << "  ";
        *Datei << "\n";
      }
    }
  }
  return nummer;
}


int Grid_base::Recursion_Cells_AVS(Index3D I, ofstream *Datei, int nummer,
                                   int number_cell_var) {
  int i;
  Index3D Ison;
  double value;

  for(i=0;i<8;++i) {
    if(Exists_Point(I.son((dir_sons)i)))
      nummer = Recursion_Cells_AVS(I.son((dir_sons)i),Datei,nummer,
                                   number_cell_var);
    else {
      if(Give_cell_typ(I.son((dir_sons)i)) == int_cell) {
        Ison = I.son((dir_sons)i);
        value = Give_cell_variable(Ison)[number_cell_var];

        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 0. Tet: WND, WNT, WST, EST 
        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 1. Tet: EST, WND, WST, ESD  
        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 2. Tet: WND, WSD, WST, ESD        
        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 3. Tet: EST, WND, ESD, END
        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 4. Tet: ENT, WNT, EST, END
        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 5. Tet: WNT, WND, EST, END
      }
    }
  }
  return nummer;
}

void Grid_base::Print_Domain_AVS(ofstream *Datei) {
  static int number_cells;
  static Index3D ind;

  if(Is_storage_initialized()==false) { 
    cout << " \n Fehler in Print_Domain_AVS! Initialisierung fehlt! "
         << endl;
    return;
  }
  // Teil 1:
  ind = ind.son(ENTd);
  ind = ind.son(WSDd);
  number_cells = Recursion_Count_Cells(ind);
  *Datei << hashtable1_occ << " " << number_cells << " 1 0 0" << endl;

  // Teil 2:
  iterate_hash1 {
    *Datei << point1->nummer << " " ;
    transform_coord(point1->ind.coordinate()).Print(Datei);
    *Datei  << "\n";
  }

  // Teil 3:
  ind = Index3D(2,2,2);
  Recursion_Cells_AVS(ind,Datei,0);

  // Teil 4:
  *Datei << "1 1" << "\n";

  // Teil 5:
  *Datei << "Konstante " << "\n";

  // Teil 6:
  iterate_hash1 {
    *Datei << point1->nummer << "  1.0 \n";
  }
}

// There are the main print functions
void Grid_base::Print_Variable_AVS(ofstream *Datei, int number_var) {
  int number_cells, number, ebene;
  int number_points;
  Pointtype typ_point;
  Index3D ind(2,2,2);

  if(Is_storage_initialized()==false) { 
    cout << " \n Fehler in Print_Variable_AVS! Initialisierung fehlt! "
         << endl;
    return;
  }
  if(Give_max_num_var() <= number_var) {
    cout << " \n Fehler in Print_Variable_AVS! number_var zu gross! " << endl;
  }
  else {
    // Teil 0: Write information
    *Datei << "# File created by ExPDE " << endl;
    *Datei << "# UCD file format for AVS " << endl;
    *Datei << "# scalar value " << endl;

    // Teil 1: Anzahl der Zellen zaehlen
    // Teil 1.1: Innerer Zellen
    ind = Index3D(2,2,2);

    number_cells = Recursion_Count_Cells(ind);
    // Teil 1.2: Randzellen
    iterate_hash2 {
      number_cells = avs_bo_cell(bocell,Datei,false,number_cells);
    }
    // Teil 1.3: Anzahl der Punkte
    // a) innen und randnah
    number_points=0;
    iterate_hash1 {
      if(point1->typ >= interior || point1->typ==parallel_p) {
        point1->nummer = number_points;
        number_points++;
      }
      else point1->nummer = -1;
    }
    // b) Randpunkte
    iterate_hash3 {
      bo2point->Set_number_avs(number_points);
      ++number_points;
    }
    // c) Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        bocell->Set_number_avs(number_points);
        number_points++;
      }
    }
    // Teil 1.4: Kopfzeile schreiben
    *Datei << number_points << " "
           << number_cells << " 1 0 0" << endl;

    // Teil 2: Koordinaten der Punkte eingeben
    // Teil 2.1 alle inneren Punkte
    iterate_hash1 {
      if(point1->nummer != -1) {
        *Datei << point1->nummer << " " ;
        transform_coord(point1->ind.coordinate()).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 2.2 alle Randpunkte (Kantenpunkte)
    iterate_hash3 {
      *Datei << bo2point->nummer << " " ;
      (bo2point->transform_coord(Give_A(),
                                 H_mesh())).Print(Datei);
      *Datei  << "\n";
    }
    // Teil 2.3 alle Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        *Datei << bocell->Number_avs() << " " ;
        (transform_coord(bocell->ind.coordinate()) + 
         bocell->Local_coord_bocellpoint()).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 3: Zellen ausgeben
    // Teil 3.1: Innere Zellen ausgeben
    ind = Index3D(2,2,2);
    number=Recursion_Cells_AVS(ind,Datei,0);
    // Teil 3.2: Randzellen
    iterate_hash2 {
      number = avs_bo_cell(bocell,Datei,true,number);
    }

    // Teil 4:
    *Datei << "1 1" << "\n";
    
    // Teil 5:
    *Datei << "variable, unit " << "\n";
    
    // Teil 6: Ausgabe der Werte der Variablen
    // Teil 6.1: Werte im Inneren
    iterate_hash1 {
      if(point1->nummer != -1) {
        *Datei << point1->nummer << "  ";
        typ_point = point1->typ;
        if(typ_point==exterior)
          cout << " Fehler in Print_Variable_AVS! " << endl;
        if(typ_point!=multigrid) {
          ebene = Give_finest_level(point1->ind);
          *Datei << Give_variable(point1->ind,ebene)[number_var];
          //     *Datei << 2.0;
          *Datei << "\n";
        }
        else *Datei << " 1.0 \n";
      }
    }
    // Teil 6.2: Werte am Rand
    iterate_hash3 {
      *Datei << bo2point->nummer << "  ";
      *Datei << Give_variable(bo2point->ind,bo2point->direction)[number_var];
      *Datei << " \n";
    }
    // Teil 6.3 Werte bei Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        *Datei << bocell->Number_avs() << " " ;
        *Datei << bocell->Give_var()[number_var];
        *Datei  << "\n";
      }
    }
  }
}

void Grid_base::Print_Variable_AVS_coarse(ofstream *Datei, int number_var,
                                          int level) {
  int number_cells, number;
  int number_points;
  Index3D ind;
  Pointtype typ_point;

  if(Is_storage_initialized()==false) { 
    cout << " \n Fehler in Print_Variable_AVS_coarse! Initialisierung fehlt! "
         << endl;
    return;
  }
  if(Give_max_num_var() <= number_var) {
    cout << " \n Fehler in Print_Variable_AVS_coarse! number_var zu gross! "
         << endl;
  }
  else {
    // Teil 0: Write information
    *Datei << "# File created by ExPDE " << endl;
    *Datei << "# UCD file format for AVS " << endl;
    *Datei << "# scalar value " << endl;

    // Teil 1: Anzahl der Zellen zaehlen
    // Teil 1.1: Innerer Zellen
    number_cells = 0;
    iterate_hash0 {
      ind = point0->Give_Index();
      if(ind.Cell_index()) {
        if(ind.Tiefe()==level+1 && 
           (point0->typ == int_cell || point0->typ == bo_cell)) {
          number_cells = number_cells + 6;
        }
      }
    }
    // Teil 1.3: Anzahl der Punkte
    // a) innen und randnah
    number_points=0;
    iterate_hash1 {
      if(point1->typ != exterior 
         && point1->ind.Tiefe()<=level     // for excluding too fine points   
         && point1->finest_level >= level) // for excluding too coarse
                                           //     multigrid points
        {
        point1->nummer = number_points;
        number_points++;
      }
      else point1->nummer = -1;
    }
    // Teil 1.4: Kopfzeile schreiben
    *Datei << number_points << " "
           << number_cells << " 1 0 0" << endl;
    
    // Teil 2: Koordinaten der Punkte eingeben
    // Teil 2.1 alle inneren Punkte
    iterate_hash1 {
      if(point1->nummer != -1) {
        *Datei << point1->nummer << " " ;
        transform_coord(point1->ind.coordinate()).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 3: Zellen ausgeben
    // Teil 3.1: Innere Zellen ausgeben
    number = 0;
    iterate_hash0 {
      ind = point0->Give_Index();
      if(ind.Cell_index()) {
        if(ind.Tiefe()==level+1 && 
           (point0->typ == int_cell || point0->typ == bo_cell)) {
          number++;
          
          *Datei << number << "  1  tet  ";
          // 0. Tet: WND, WNT, WST, EST 
          *Datei << Give_Nummer(ind.neighbour(WNTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(WNDd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(WSTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(ESTd))  << "  ";
          *Datei << "\n";
          
          number++;
          *Datei << number << "  1  tet  ";
          // 1. Tet: EST, WND, WST, ESD  
          *Datei << Give_Nummer(ind.neighbour(ESTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(WNDd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(WSTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(ESDd))  << "  ";
          *Datei << "\n";
        
          number++;
          *Datei << number << "  1  tet  ";
          // 2. Tet: WND, WSD, WST, ESD        
          *Datei << Give_Nummer(ind.neighbour(WNDd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(WSDd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(WSTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(ESDd))  << "  ";
          *Datei << "\n";
          
          number++;
          *Datei << number << "  1  tet  ";
          // 3. Tet: EST, WND, ESD, END
          *Datei << Give_Nummer(ind.neighbour(ESTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(WNDd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(ESDd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(ENDd))  << "  ";
          *Datei << "\n";
          
          number++;
          *Datei << number << "  1  tet  ";
          // 4. Tet: ENT, WNT, EST, END
          *Datei << Give_Nummer(ind.neighbour(ENTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(WNTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(ESTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(ENDd))  << "  ";
          *Datei << "\n";
          
          number++;
          *Datei << number << "  1  tet  ";
          // 5. Tet: WNT, WND, EST, END
          *Datei << Give_Nummer(ind.neighbour(WNTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(WNDd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(ESTd))  << "  ";
          *Datei << Give_Nummer(ind.neighbour(ENDd))  << "  ";
          *Datei << "\n";
        }
      }
    }
    
    // Teil 4:
    *Datei << "1 1" << "\n";
    
    // Teil 5:
    *Datei << "variable, unit " << "\n";
    
    // Teil 6: Ausgabe der Werte der Variablen
    // Teil 6.1: Werte im Inneren
    double* vars;
    iterate_hash1 {
      if(point1->nummer != -1) {
        *Datei << point1->nummer << "  ";
        typ_point = point1->typ;
        if(typ_point==exterior)
          cout << " Fehler in Print_Variable_AVS_coarse! " << endl;
        vars = Give_variable_slow(point1->ind,level);
        if(vars!=NULL)
          *Datei << vars[number_var];
        else
          *Datei << 0.0;
        *Datei << "\n";
      }
    }
  }
}



void Grid_base::Print_Variable_AVS(ofstream *Datei, 
                                   int number_var_a,
                                   int number_var_b,
                                   int number_var_c) {
  int number_cells, number, ebene;
  int number_points;
  Index3D ind;
  Pointtype typ_point;

  if(Is_storage_initialized()==false) { 
    cout << " \n Fehler in Print_Variable_AVS! Initialisierung fehlt! "
         << endl;
    return;
  }
  if(Give_max_num_var() <= number_var_a || 
     Give_max_num_var() <= number_var_b || 
     Give_max_num_var() <= number_var_c ) {
    cout << " \n Fehler in Print_Variable_AVS! number_var zu gross! " << endl;
  }
  else {
    // Teil 0: Write information
    *Datei << "# File created by ExPDE " << endl;
    *Datei << "# UCD file format for AVS " << endl;
    *Datei << "# vector field " << endl;

    // Teil 1: Anzahl der Zellen zaehlen
    // Teil 1.1: Innerer Zellen
    ind = Index3D(2,2,2);
    number_cells = Recursion_Count_Cells(ind);
    // Teil 1.2: Randzellen
    iterate_hash2 {
      number_cells = avs_bo_cell(bocell,Datei,false,number_cells);
    }
    // Teil 1.3: Anzahl der Punkte der Randzellfreiheitsgrad
    // a) innen und randnah
    number_points=0;
    iterate_hash1 {
      if(point1->typ >= interior || point1->typ==parallel_p) {
        point1->nummer = number_points;
        number_points++;
      }
      else point1->nummer = -1;
    }
    // b) Randpunkte
    iterate_hash3 {
      bo2point->Set_number_avs(number_points);
      ++number_points;
    }
    // c) Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        bocell->Set_number_avs(number_points);
        number_points++;
      }
    }
    // Teil 1.4: Kopfzeile schreiben
    *Datei << number_points << " "
           << number_cells << " 3 0 0" << endl;

    // Teil 2: Koordinaten der Punkte eingeben
    // Teil 2.1 alle inneren Punkte
    iterate_hash1 {
      if(point1->nummer != -1) {
        *Datei << point1->nummer << " " ;
        transform_coord(point1->ind.coordinate()).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 2.2 alle Randpunkte (Kantenpunkte)
    iterate_hash3 {
      *Datei << bo2point->nummer << " " ;
      (bo2point->transform_coord(Give_A(),
                                 H_mesh())).Print(Datei);
      *Datei  << "\n";
    }
    // Teil 2.3 alle Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        *Datei << bocell->Number_avs() << " " ;
        /*
        transform_coord(bocell->ind.coordinate() + 
                        bocell->Local_coord_bocellpoint()).Print(Datei);
        */
        (transform_coord(bocell->ind.coordinate()) + 
         bocell->Local_coord_bocellpoint()).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 3: Zellen ausgeben
    // Teil 3.1: Innere Zellen ausgeben
    ind = Index3D(2,2,2);
    number=Recursion_Cells_AVS(ind,Datei,0);
    // Teil 3.2: Randzellen
    iterate_hash2 {
      number = avs_bo_cell(bocell,Datei,true,number);
    }

    // Teil 4:
    *Datei << "3 1 1 1" << "\n";
    
    // Teil 5:
    *Datei << "x-coord " << "\n";
    *Datei << "y-coord " << "\n";
    *Datei << "z-coord " << "\n";
    
    // Teil 6: Ausgabe der Werte der Variablen
    // Teil 6.1: Werte im Inneren
    iterate_hash1 {
      if(point1->nummer != -1) {
        *Datei << point1->nummer << "  ";
        typ_point = point1->typ;
        if(typ_point==exterior)
          cout << " Fehler in Print_Variable_AVS! " << endl;
        if(typ_point!=multigrid) {
          // fehlt ebene richtig angeben
          ebene = Give_finest_level(point1->ind);
          *Datei << Give_variable(point1->ind,ebene)[number_var_a] << "  ";
          *Datei << Give_variable(point1->ind,ebene)[number_var_b] << "  ";
          *Datei << Give_variable(point1->ind,ebene)[number_var_c];
          //     *Datei << 2.0;
          *Datei << "\n";
        }
        else *Datei << " 1.0 1.0 1.0 \n";
      }
    }
    // Teil 6.2: Werte am Rand
    iterate_hash3 {
      *Datei << bo2point->nummer << "  ";
      *Datei << Give_variable(bo2point->ind,bo2point->direction)[number_var_a]
             << "   ";
      *Datei << Give_variable(bo2point->ind,bo2point->direction)[number_var_b]
             << "   ";
      *Datei << Give_variable(bo2point->ind,bo2point->direction)[number_var_c];
      *Datei << " \n";
    }
    // Teil 6.3 Werte bei Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        *Datei << bocell->Number_avs() << " " ;
        *Datei << bocell->Give_var()[number_var_a] << "  ";
        *Datei << bocell->Give_var()[number_var_b] << "  ";
        *Datei << bocell->Give_var()[number_var_c];
        *Datei  << "\n";
      }
    }
  }
}

void Grid_base::Print_Variable_AVS_moved(ofstream *Datei, int number_var,
                                         int number_var_a,
                                         int number_var_b, 
                                         int number_var_c) {
  int number_cells, number, ebene;
  int number_points;
  Index3D ind;
  Pointtype typ_point;

  if(Is_storage_initialized()==false) { 
    cout << " \n Fehler in Print_Variable_AVS! Initialisierung fehlt! "
         << endl;
    return;
  }
  if(Give_max_num_var() <= number_var) {
    cout << " \n Fehler in Print_Variable_AVS! number_var zu gross! " << endl;
  }
  else {
    // Teil 0: Write information
    *Datei << "# File created by ExPDE " << endl;
    *Datei << "# UCD file format for AVS " << endl;
    *Datei << "# scalar value moved by a vector" << endl;

    // Teil 1: Anzahl der Zellen zaehlen
    // Teil 1.1: Innerer Zellen
    ind = Index3D(2,2,2);
    number_cells = Recursion_Count_Cells(ind);
    // Teil 1.2: Randzellen
    iterate_hash2 {
      number_cells = avs_bo_cell(bocell,Datei,false,number_cells);
    }
    // Teil 1.3: Anzahl der Punkte
    // a) innen und randnah
    number_points=0;
    iterate_hash1 {
      if(point1->typ >= interior || point1->typ==parallel_p) {
        point1->nummer = number_points;
        number_points++;
      }
      else point1->nummer = -1;
    }
    // b) Randpunkte
    iterate_hash3 {
      bo2point->Set_number_avs(number_points);
      ++number_points;
    }
    // c) Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        bocell->Set_number_avs(number_points);
        number_points++;
      }
    }
    // Teil 1.4: Kopfzeile schreiben
    *Datei << number_points << " "
           << number_cells << " 1 0 0" << endl;

    // Teil 2: Koordinaten der Punkte eingeben
    // Teil 2.1 alle inneren Punkte
    iterate_hash1 {
      if(point1->nummer != -1) {
        *Datei << point1->nummer << " " ;
        (transform_coord(point1->ind.coordinate()) +
         D3vector(Give_variable(point1->ind,Max_level())[number_var_a],
                  Give_variable(point1->ind,Max_level())[number_var_b],
                  Give_variable(point1->ind,Max_level())[number_var_c])
         ).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 2.2 alle Randpunkte (Kantenpunkte)
    iterate_hash3 {
      *Datei << bo2point->nummer << " " ;
      ((bo2point->transform_coord(Give_A(),
                                 H_mesh())) +
       D3vector(Give_variable(bo2point->ind,bo2point->direction)[number_var_a],
                Give_variable(bo2point->ind,bo2point->direction)[number_var_b],
                Give_variable(bo2point->ind,bo2point->direction)[number_var_c])
       ).Print(Datei);
      *Datei  << "\n";
    }
    // Teil 2.3 alle Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        *Datei << bocell->Number_avs() << " " ;
        (transform_coord(bocell->ind.coordinate()) + 
         bocell->Local_coord_bocellpoint() +
         D3vector(bocell->Give_var()[number_var_a],
                  bocell->Give_var()[number_var_b],
                  bocell->Give_var()[number_var_c])
         ).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 3: Zellen ausgeben
    // Teil 3.1: Innere Zellen ausgeben
    ind = Index3D(2,2,2);
    number=Recursion_Cells_AVS(ind,Datei,0);
    // Teil 3.2: Randzellen
    iterate_hash2 {
      number = avs_bo_cell(bocell,Datei,true,number);
    }

    // Teil 4:
    *Datei << "1 1" << "\n";
    
    // Teil 5:
    *Datei << "variable, unit " << "\n";
    
    // Teil 6: Ausgabe der Werte der Variablen
    // Teil 6.1: Werte im Inneren
    iterate_hash1 {
      if(point1->nummer != -1) {
        *Datei << point1->nummer << "  ";
        typ_point = point1->typ;
        if(typ_point==exterior)
          cout << " Fehler in Print_Variable_AVS! " << endl;
        if(typ_point!=multigrid) {
          ebene = Give_finest_level(point1->ind);
          *Datei << Give_variable(point1->ind,ebene)[number_var];
          //     *Datei << 2.0;
          *Datei << "\n";
        }
        else *Datei << " 1.0 \n";
      }
    }
    // Teil 6.2: Werte am Rand
    iterate_hash3 {
      *Datei << bo2point->nummer << "  ";
      *Datei << Give_variable(bo2point->ind,bo2point->direction)[number_var];
      *Datei << " \n";
    }
    // Teil 6.3 Werte bei Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        *Datei << bocell->Number_avs() << " " ;
        *Datei << bocell->Give_var()[number_var];
        *Datei  << "\n";
      }
    }
  }
}
// for cell variable

void Grid_base::Print_Cell_Variable_AVS(ofstream *Datei, int number_cell_var) {
  int number_cells, number;
  int number_points;
  Index3D ind;

  if(Is_storage_initialized()==false) { 
    cout << " \n Fehler in Print_Cell_Variable_AVS! Initialisierung fehlt! "
         << endl;
    return;
  }
  if(Give_max_num_cell_var() <= number_cell_var) {
    cout << " \n Fehler in Print_Cell_Variable_AVS!  " << endl;
  }
  else {
    // Teil 0: Write information
    *Datei << "# File created by ExPDE " << endl;
    *Datei << "# UCD file format for AVS " << endl;
    *Datei << "# scalar value on cells" << endl;

    // Teil 1: Anzahl der Zellen zaehlen
    // Teil 1.1: Innerer Zellen
    ind = Index3D(2,2,2);
    number_cells = Recursion_Count_Cells(ind);
    // Teil 1.2: Randzellen
    iterate_hash2 {
      number_cells = avs_bo_cell(bocell,Datei,false,number_cells);
    }
    // Teil 1.3: Anzahl der Punkte
    // a) innen und randnah
    number_points=0;
    iterate_hash1 {
      if(point1->typ >= interior || point1->typ==parallel_p) {
        point1->nummer = number_points;
        number_points++;
      }
      else point1->nummer = -1;
    }
    // b) Randpunkte
    iterate_hash3 {
      bo2point->Set_number_avs(number_points);
      ++number_points;
    }
    // c) Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        bocell->Set_number_avs(number_points);
        number_points++;
      }
    }
    // Teil 1.4: Kopfzeile schreiben
    *Datei << number_points << " "
           << number_cells << " 0 1 0" << endl;

    // Teil 2: Koordinaten der Punkte eingeben
    // Teil 2.1 alle inneren Punkte
    iterate_hash1 {
      if(point1->nummer != -1) {
        *Datei << point1->nummer << " " ;
        transform_coord(point1->ind.coordinate()).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 2.2 alle Randpunkte (Kantenpunkte)
    iterate_hash3 {
      *Datei << bo2point->nummer << " " ;
      (bo2point->transform_coord(Give_A(),
                                 H_mesh())).Print(Datei);
      *Datei  << "\n";
    }
    // Teil 2.3 alle Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        *Datei << bocell->Number_avs() << " " ;
        (transform_coord(bocell->ind.coordinate()) + 
         bocell->Local_coord_bocellpoint()).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 3: Zellen ausgeben
    // Teil 3.1: Innere Zellen ausgeben
    ind = Index3D(2,2,2);
    number=Recursion_Cells_AVS(ind,Datei,0);
    // Teil 3.2: Randzellen
    iterate_hash2 {
      number = avs_bo_cell(bocell,Datei,true,number);
    }

    // Teil 4:
    *Datei << "1 1" << "\n";
    
    // Teil 5:
    *Datei << "variable, unit " << "\n";
    
    // Teil 6: Ausgabe der Werte der Zell-Variablen
    // Teil 3.1: Innere Zellen ausgeben
    ind = Index3D(2,2,2);
    number=Recursion_Cells_AVS(ind,Datei,0,number_cell_var);
    // Teil 3.2: Randzellen
    iterate_hash2 {
      number = avs_bo_cell(bocell,Datei,number,number_cell_var);
    }
  }
}

void Grid_base::Print_Cell_Variable_AVS_moved(ofstream *Datei, 
                                              int number_cell_var,
                                              int number_var_a,
                                              int number_var_b, 
                                              int number_var_c) {
  int number_cells, number;
  int number_points;
  Index3D ind;

  if(Is_storage_initialized()==false) { 
    cout << " \n Fehler in Print_Cell_Variable_AVS! Initialisierung fehlt! "
         << endl;
    return;
  }
  if(Give_max_num_cell_var() <= number_cell_var) {
    cout << " \n Fehler in Print_Cell_Variable_AVS!  " << endl;
  }
  else {
    // Teil 0: Write information
    *Datei << "# File created by ExPDE " << endl;
    *Datei << "# UCD file format for AVS " << endl;
    *Datei << "# scalar value on cells" << endl;

    // Teil 1: Anzahl der Zellen zaehlen
    // Teil 1.1: Innerer Zellen
    ind = Index3D(2,2,2);
    number_cells = Recursion_Count_Cells(ind);
    // Teil 1.2: Randzellen
    iterate_hash2 {
      number_cells = avs_bo_cell(bocell,Datei,false,number_cells);
    }
    // Teil 1.3: Anzahl der Punkte
    // a) innen und randnah
    number_points=0;
    iterate_hash1 {
      if(point1->typ >= interior || point1->typ==parallel_p) {
        point1->nummer = number_points;
        number_points++;
      }
      else point1->nummer = -1;
    }
    // b) Randpunkte
    iterate_hash3 {
      bo2point->Set_number_avs(number_points);
      ++number_points;
    }
    // c) Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        bocell->Set_number_avs(number_points);
        number_points++;
      }
    }
    // Teil 1.4: Kopfzeile schreiben
    *Datei << number_points << " "
           << number_cells << " 0 1 0" << endl;

    // Teil 2: Koordinaten der Punkte eingeben
    // Teil 2.1 alle inneren Punkte
    iterate_hash1 {
      if(point1->nummer != -1) {
        *Datei << point1->nummer << " " ;
        (transform_coord(point1->ind.coordinate()) +
         D3vector(Give_variable(point1->ind,Max_level())[number_var_a],
                  Give_variable(point1->ind,Max_level())[number_var_b],
                  Give_variable(point1->ind,Max_level())[number_var_c])
         ).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 2.2 alle Randpunkte (Kantenpunkte)
    iterate_hash3 {
      *Datei << bo2point->nummer << " " ;
      ((bo2point->transform_coord(Give_A(),
                                  H_mesh())) +
       D3vector(Give_variable(bo2point->ind,bo2point->direction)[number_var_a],
                Give_variable(bo2point->ind,bo2point->direction)[number_var_b],
                Give_variable(bo2point->ind,bo2point->direction)[number_var_c])
       ).Print(Datei);
      *Datei  << "\n";
    }
    // Teil 2.3 alle Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Exists_bocellpoint()) {
        *Datei << bocell->Number_avs() << " " ;
        (transform_coord(bocell->ind.coordinate()) + 
         bocell->Local_coord_bocellpoint() +
         D3vector(bocell->Give_var()[number_var_a],
                  bocell->Give_var()[number_var_b],
                  bocell->Give_var()[number_var_c])
         ).Print(Datei);
        *Datei  << "\n";
      }
    }
    // Teil 3: Zellen ausgeben
    // Teil 3.1: Innere Zellen ausgeben
    ind = Index3D(2,2,2);
    number=Recursion_Cells_AVS(ind,Datei,0);
    // Teil 3.2: Randzellen
    iterate_hash2 {
      number = avs_bo_cell(bocell,Datei,true,number);
    }

    // Teil 4:
    *Datei << "1 1" << "\n";
    
    // Teil 5:
    *Datei << "variable, unit " << "\n";
    
    // Teil 6: Ausgabe der Werte der Zell-Variablen
    // Teil 3.1: Innere Zellen ausgeben
    ind = Index3D(2,2,2);
    number=Recursion_Cells_AVS(ind,Datei,0,number_cell_var);
    // Teil 3.2: Randzellen
    iterate_hash2 {
      number = avs_bo_cell(bocell,Datei,number,number_cell_var);
    }
  }
}






// b)  for Gnuplot

// Print simple gnuplot-file
void Grid_base::Print_Grid_Gnuplot(ofstream *Datei) {
  *Datei << "# File created by ExPDE " << endl;
  *Datei << "# file format for Gnuplot " << endl;
  *Datei << "# visualization of the surface ";
  
  Index3D I;

  // Print minimal and maximal points
  *Datei << "\n";
  Give_A().Print(Datei);
  *Datei << "\n \n" ;
  (Give_A()+D3vector(H_mesh(),H_mesh(),H_mesh())).Print(Datei);
  *Datei << "\n" ;

  //  Print surface Tetrahedral
  Tetraeder_storage* tet;
  iterate_hash2 {
    for(tet=bocell->Give_boundary_tets();tet!=NULL;tet=tet->Next()) {
      I = bocell->Give_Index();
      *Datei << "\n" ;
      *Datei << "\n" ;
      transform_coord(I.neighbour(bocell->corner(tet->N0())),
                      bocell->edge_dir(tet->N0())).Print(Datei);
      *Datei << endl;
      transform_coord(I.neighbour(bocell->corner(tet->N1())),
                      bocell->edge_dir(tet->N1())).Print(Datei);
      *Datei << endl;
      transform_coord(I.neighbour(bocell->corner(tet->N2())),
                      bocell->edge_dir(tet->N2())).Print(Datei);
    }
  }
}


inline D3vector MAX_of_vectors(D3vector& A, D3vector& B) {
  return D3vector(MAX(A.x,B.x),MAX(A.y,B.y),MAX(A.z,B.z));
}
inline D3vector MIN_of_vectors(D3vector& A, D3vector& B) {
  return D3vector(MIN(A.x,B.x),MIN(A.y,B.y),MIN(A.z,B.z));
}

void Grid_base::Print_Grid_Gnuplot_moved(ofstream *Datei,
                                         int number_var_a,
                                         int number_var_b, 
                                         int number_var_c) {
  *Datei << "# File created by ExPDE " << endl;
  *Datei << "# file format for Gnuplot " << endl;
  *Datei << "# visualization of the surface ";
  
  Index3D I;
  D3vector loc0, loc1, loc2;
  D3vector locmax, locmin;

  Tetraeder_storage* tet;

  locmax = D3vector(0.0,0.0,0.0);
  locmin = D3vector(0.0,0.0,0.0);

  //  Calc corner of moved maximal square
  iterate_hash2 {
    for(tet=bocell->Give_boundary_tets();tet!=NULL;tet=tet->Next()) {
      I = bocell->Give_Index();
      loc0 = transform_coord(I.neighbour(bocell->corner(tet->N0())),
                             bocell->edge_dir(tet->N0())) +
        D3vector(bocell->vars[tet->N0()][number_var_a],
                 bocell->vars[tet->N0()][number_var_b],
                 bocell->vars[tet->N0()][number_var_c]);
      loc1 = transform_coord(I.neighbour(bocell->corner(tet->N1())),
                             bocell->edge_dir(tet->N1())) +
        D3vector(bocell->vars[tet->N1()][number_var_a],
                 bocell->vars[tet->N1()][number_var_b],
                 bocell->vars[tet->N1()][number_var_c]);
      loc2 = transform_coord(I.neighbour(bocell->corner(tet->N2())),
                             bocell->edge_dir(tet->N2())) +
        D3vector(bocell->vars[tet->N2()][number_var_a],
                 bocell->vars[tet->N2()][number_var_b],
                 bocell->vars[tet->N2()][number_var_c]);

      locmax=MAX_of_vectors(locmax,loc0);
      locmax=MAX_of_vectors(locmax,loc1);
      locmax=MAX_of_vectors(locmax,loc2);

      locmin=MIN_of_vectors(locmin,loc0);
      locmin=MIN_of_vectors(locmin,loc1);
      locmin=MIN_of_vectors(locmin,loc2);
    }
  }

  // Print minimal and maximal points
  /*
  *Datei << "\n";
  Give_A().Print(Datei);
  *Datei << "\n \n" ;
  (Give_A()+D3vector(H_mesh(),H_mesh(),H_mesh())).Print(Datei);
  *Datei << "\n";
  */
  double h_maxx;
  h_maxx = MAX(locmax.x-locmin.x,locmax.y-locmin.y,locmax.z-locmin.z);
  *Datei << "\n";
  locmin.Print(Datei);
  *Datei << "\n \n" ;
  (locmin+D3vector(h_maxx,h_maxx,h_maxx)).Print(Datei);
  *Datei << "\n";

  //  Print surface Tetrahedral
  iterate_hash2 {
    for(tet=bocell->Give_boundary_tets();tet!=NULL;tet=tet->Next()) {
      I = bocell->Give_Index();
      *Datei << "\n" ;
      *Datei << "\n" ;

      loc0 = transform_coord(I.neighbour(bocell->corner(tet->N0())),
                             bocell->edge_dir(tet->N0())) +
        D3vector(bocell->vars[tet->N0()][number_var_a],
                 bocell->vars[tet->N0()][number_var_b],
                 bocell->vars[tet->N0()][number_var_c]);
      loc1 = transform_coord(I.neighbour(bocell->corner(tet->N1())),
                             bocell->edge_dir(tet->N1())) +
        D3vector(bocell->vars[tet->N1()][number_var_a],
                 bocell->vars[tet->N1()][number_var_b],
                 bocell->vars[tet->N1()][number_var_c]);
      loc2 = transform_coord(I.neighbour(bocell->corner(tet->N2())),
                             bocell->edge_dir(tet->N2())) +
        D3vector(bocell->vars[tet->N2()][number_var_a],
                 bocell->vars[tet->N2()][number_var_b],
                 bocell->vars[tet->N2()][number_var_c]);
      
      loc0.Print(Datei);  *Datei << endl;
      loc1.Print(Datei);  *Datei << endl;
      loc2.Print(Datei);
    }
  }
}

//  Test

void Grid_base::Print_test_AVS(ofstream *Datei, int type) {
  int number_cells, number;
  int number_points;
  Index3D ind;
  Pointtype typ_point;

  // Teil 0: Write information
  *Datei << "# File created by ExPDE " << endl;
  *Datei << "# UCD file format for AVS " << endl;
  *Datei << "# scalar value " << endl;
  
  // Teil 1: Anzahl der Zellen zaehlen
  // Teil 1.1: Innerer Zellen
  ind = Index3D(2,2,2);
  number_cells = Recursion_Count_Cells(ind);
  // Teil 1.2: Randzellen
  iterate_hash2 {
    number_cells = avs_bo_cell(bocell,Datei,false,number_cells);
  }
  // Teil 1.3: Anzahl der Punkte
  // a) innen und randnah
  number_points=0;
  iterate_hash1 {
    //    if(point1->typ >= interior || point1->typ==parallel_p) {
    if(point1->typ != exterior) {
      point1->nummer = number_points;
      number_points++;
    }
    else point1->nummer = -1;
  }
  // b) Randpunkte
  iterate_hash3 {
    bo2point->Set_number_avs(number_points);
    ++number_points;
  }
  // c) Randzellfreiheitsgrad
  iterate_hash2 {
    if(bocell->Exists_bocellpoint()) {
      bocell->Set_number_avs(number_points);
      number_points++;
    }
  }
  // Teil 1.4: Kopfzeile schreiben
  *Datei << number_points << " "
         << number_cells << " 1 0 0" << endl;
  
  // Teil 2: Koordinaten der Punkte eingeben
  // Teil 2.1 alle inneren Punkte
  iterate_hash1 {
    if(point1->nummer != -1) {
      *Datei << point1->nummer << " " ;
      transform_coord(point1->ind.coordinate()).Print(Datei);
      *Datei  << "\n";
    }
  }
  // Teil 2.2 alle Randpunkte (Kantenpunkte)
  iterate_hash3 {
    *Datei << bo2point->nummer << " " ;
    (bo2point->transform_coord(Give_A(),
                               H_mesh())).Print(Datei);
    *Datei  << "\n";
  }
  // Teil 2.3 alle Randzellfreiheitsgrad
  iterate_hash2 {
    if(bocell->Exists_bocellpoint()) {
      *Datei << bocell->Number_avs() << " " ;
      (transform_coord(bocell->ind.coordinate()) + 
       bocell->Local_coord_bocellpoint()).Print(Datei);
      *Datei  << "\n";
    }
  }
  // Teil 3: Zellen ausgeben
  // Teil 3.1: Innere Zellen ausgeben
  ind = Index3D(2,2,2);
  number=Recursion_Cells_AVS(ind,Datei,0);
  // Teil 3.2: Randzellen
  iterate_hash2 {
    number = avs_bo_cell(bocell,Datei,true,number);
  }
  
  // Teil 4:
  *Datei << "1 1" << "\n";
  
  // Teil 5:
  *Datei << "variable, unit " << "\n";
  
  // Teil 6: Ausgabe der Werte der Variablen
  // Teil 6.1: Werte im Inneren
  iterate_hash1 {
    if(point1->nummer != -1) {
      *Datei << point1->nummer << "  ";
      typ_point = point1->typ;
      if(typ_point==exterior)
        cout << " Fehler in Print_Variable_AVS! " << endl;
      //      if(typ_point!=multigrid) {
        if(type==0) *Datei << 2.0;
        if(type==1) *Datei << point1->typ;
        *Datei << "\n";
        //      }
        //      else *Datei << " 1.0 \n";
    }
  }
  // Teil 6.2: Werte am Rand
  iterate_hash3 {
    *Datei << bo2point->nummer << "  ";
    if(type==0) *Datei << 2.0;
    if(type==1) *Datei << -2.0;
    *Datei << " \n";
  }
  // Teil 6.3 Werte bei Randzellfreiheitsgrad
  iterate_hash2 {
    if(bocell->Exists_bocellpoint()) {
      *Datei << bocell->Number_avs() << " " ;
      if(type==0) *Datei << 2.0;
      if(type==1) *Datei << -1.0;
      *Datei  << "\n";
    }
  }
}


void Grid_base::Test(int type) {
  if(type==0) {
    cout << " \n Print boundary cells: " << endl;
    iterate_hash2 {
      cout << " center point: "; bocell->ind.coordinate().Print(); 
      cout << endl;
    }
  }
  if(type==1) {
    cout << " \n Print cells: " << endl;
    iterate_hash0 {
      if(point0->Give_Index().Cell_index()
         && point0->typ==int_cell
         && point0->Give_Index().Tiefe()==4) {
        cout << " cell point: " << point0->typ << " where: ";
        transform_coord(point0->Give_Index().coordinate()).Print();
        cout << endl;
      }
    }
  }
  if(type==7) {
    cout << " \n Print cells boundary 0: " << endl;
    int Anzahl;
    Anzahl=0;
    iterate_hash0 {
      if(point0->Give_Index().Cell_index())
        if(point0->typ==fine_bo_cell) {
          Anzahl++;
          cout << " cell point: " << point0->typ << " where: ";
          point0->Give_Index().coordinate().Print();
          cout << endl;
      }
    }
    cout << " Anzahl: " << Anzahl << endl;
  }
  if(type==2) {
    cout << " \n Test Corners of boundary cells: " << endl;
    bool error;
    int k;
    iterate_hash2 {
      error=false;
      for(k=0;k<8;++k) {
        if(Exists_Point(bocell->ind.neighbour((dir_sons)k))==false)
          error=true;
      }
      if(error) {
        cout << " error at center point: "; bocell->ind.coordinate().Print(); 
        cout << endl;
      }
    }
  }
  if(type==3) {
    cout << " \n Test Corners of cells of type fine_bo_cell: " << endl;
    bool error;
    int k;
    iterate_hash0 {
      if(point0->Give_Index().Cell_index()) {
        if(point0->typ==fine_bo_cell) {
          error=false;
          for(k=0;k<8;++k) {
            if(Exists_Point(point0->ind.neighbour((dir_sons)k))==false)
              error=true;
          }
          if(error) {
            cout << " error at center point: "; point0->ind.coordinate().Print(); 
            point0->ind.Print();
            cout << endl;
          }
        }
      }
    }
  }
  if(type==4) {
    if(my_rank==0) cout << " \n boundary points: " << endl;
    for(int p=0;p<give_number_of_processes();++p) {
      if(my_rank==p) {
        cout << " process: " << my_rank << endl;
        iterate_hash3 {
          if(I_am_responsible_for(bo2point->Give_boindex())) {
            cout << " point: ";
            transform_coord(bo2point->ind.coordinate()).Print();
            cout << ", " << bo2point->direction << endl;
          }
        }
      }
      MPI_Barrier(comm);
    }
  }
  if(type==5) {
    if(my_rank==0) cout << " \n interior points: " << endl;
    for(int p=0;p<give_number_of_processes();++p) {
      if(my_rank==p) {
        cout << " process: " << my_rank << endl;
        iterate_hash1 {
          if(I_am_responsible_for(point1->Give_Index()) && 
             point1->typ>=interior) {
            cout << " point: ";
            transform_coord(point1->ind.coordinate()).Print();
            cout <<  " typ: " << point1->typ << endl;
          }
        }
      }
      MPI_Barrier(comm);
    }
  }
 if(type==6) {
    if(my_rank==0) cout << " \n interior points: " << endl;
    for(int p=0;p<give_number_of_processes();++p) {
      if(my_rank==p) {
        cout << " process: " << my_rank << endl;
        iterate_hash1 {
          //      if(I_am_responsible_for(point1->Give_Index()) && 
          //         point1->typ>=interior) {
          //      if(point1->typ>=interior || point1->typ==parallel_p) {
          if(true) {
            cout << " point: ";
            transform_coord(point1->ind.coordinate()).Print();
            cout <<  " typ: " << point1->typ << endl;
          }
        }
      }
      MPI_Barrier(comm);
    }
  }
}


// Test for cells

void Grid_base::Print_Test_Cell(ofstream *Datei) {
  int number_cells, number;
  int number_points;
  Index3D ind;

  // Teil 0: Write information
  *Datei << "# File created by ExPDE " << endl;
  *Datei << "# UCD file format for AVS " << endl;
  *Datei << "# test cell" << endl;

  // Teil 1: Anzahl der Zellen zaehlen
  // Teil 1.1: Innerer Zellen
  ind = Index3D(2,2,2);
  number_cells = Recursion_Count_Cells(ind);
  // Teil 1.2: Randzellen
  iterate_hash2 {
    number_cells = avs_bo_cell(bocell,Datei,false,number_cells);
  }
  // Teil 1.3: Anzahl der Punkte
  // a) innen und randnah
  number_points=0;
  iterate_hash1 {
    if(point1->typ >= interior || point1->typ==parallel_p) {
      point1->nummer = number_points;
      number_points++;
    }
    else point1->nummer = -1;
  }
  // b) Randpunkte
  iterate_hash3 {
    bo2point->Set_number_avs(number_points);
    ++number_points;
  }
  // c) Randzellfreiheitsgrad
  iterate_hash2 {
    if(bocell->Exists_bocellpoint()) {
      bocell->Set_number_avs(number_points);
      number_points++;
    }
  }
  // Teil 1.4: Kopfzeile schreiben
  *Datei << number_points << " "
         << number_cells << " 0 1 0" << endl;
  
  // Teil 2: Koordinaten der Punkte eingeben
  // Teil 2.1 alle inneren Punkte
  iterate_hash1 {
    if(point1->nummer != -1) {
      *Datei << point1->nummer << " " ;
      transform_coord(point1->ind.coordinate()).Print(Datei);
      *Datei  << "\n";
    }
  }
  // Teil 2.2 alle Randpunkte (Kantenpunkte)
  iterate_hash3 {
    *Datei << bo2point->nummer << " " ;
    (bo2point->transform_coord(Give_A(),
                               H_mesh())).Print(Datei);
    *Datei  << "\n";
  }
  // Teil 2.3 alle Randzellfreiheitsgrad
  iterate_hash2 {
    if(bocell->Exists_bocellpoint()) {
      *Datei << bocell->Number_avs() << " " ;
      (transform_coord(bocell->ind.coordinate()) + 
       bocell->Local_coord_bocellpoint()).Print(Datei);
      *Datei  << "\n";
    }
  }
  // Teil 3: Zellen ausgeben
  // Teil 3.1: Innere Zellen ausgeben
  ind = Index3D(2,2,2);
  number=Recursion_Cells_AVS(ind,Datei,0);
  // Teil 3.2: Randzellen
  iterate_hash2 {
    number = avs_bo_cell(bocell,Datei,true,number);
  }
  
  // Teil 4:
  *Datei << "1 1" << "\n";
  
  // Teil 5:
  *Datei << "variable, unit " << "\n";
  
  // Teil 6: Ausgabe der Werte der Zell-Variablen
  // Teil 3.1: Innere Zellen ausgeben
  ind = Index3D(2,2,2);
  number=Recursion_Cell_typ_AVS(ind,Datei,0);
  // Teil 3.2: Randzellen
  iterate_hash2 {
    number = avs_bo_cell_typ(bocell,Datei,number,-1.0);
  }
}

int Grid_base::avs_bo_cell_typ(BoCell* bo,ofstream *datei, 
                           int number,double value) {
  Tetraeder_storage* tets;

  for(tets = bo->Give_tets();tets!=NULL;tets = tets->Next()) {
    ++number; 
    *datei << number << "  " << value << "\n";
  }
  return number;
}


int Grid_base::Recursion_Cell_typ_AVS(Index3D I, ofstream *Datei, int nummer) {
  int i;
  Index3D Ison;
  double value;

  for(i=0;i<8;++i) {
    if(Exists_Point(I.son((dir_sons)i)))
      nummer = Recursion_Cell_typ_AVS(I.son((dir_sons)i),Datei,nummer);
    else {
      if(Give_cell_typ(I.son((dir_sons)i)) == int_cell) {
        Ison = I.son((dir_sons)i);
        value = Give_cell_typ(Ison);

        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 0. Tet: WND, WNT, WST, EST 
        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 1. Tet: EST, WND, WST, ESD  
        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 2. Tet: WND, WSD, WST, ESD        
        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 3. Tet: EST, WND, ESD, END
        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 4. Tet: ENT, WNT, EST, END
        nummer++;
        *Datei << nummer << "  " << value << "\n";
        // 5. Tet: WNT, WND, EST, END
      }
    }
  }
  return nummer;
}

void Grid_base::Print_Test_all_Cell(ofstream *Datei, int level) {
  int number_cells, i;
  int number_points;
  Index3D I;
  D3vector v;
    
  // Teil 0: Write information
  *Datei << "# File created by ExPDE " << endl;
  *Datei << "# UCD file format for AVS " << endl;
  *Datei << "# test cell" << endl;

  // Teil 1: Anzahl der Zellen zaehlen
  // Teil 1.1: Innerer Zellen
  number_cells = 0;
  iterate_hash0 {
    I = point0->Give_Index();
    if(I.Cell_index() && I.Tiefe() == level) {
      number_cells = number_cells + 1;
    }
  }
  // Teil 1.3: Anzahl der Punkte
  number_points = number_cells * 8;

  // Teil 1.4: Kopfzeile schreiben
  *Datei << number_points << " "
         << number_cells << " 0 1 0" << endl;
  
  // Teil 2: Koordinaten der Punkte eingeben
  // Teil 2.1 alle inneren Punkte
  number_cells = 0;
  iterate_hash0 {
    I = point0->Give_Index();
    if(I.Cell_index() && I.Tiefe() == level) {
      for(i=0;i<8;++i) {
        v = transform_coord(I.neighbour((dir_sons)i).coordinate());
        *Datei << number_cells*8+i
               << " " << v.x  
               << " " << v.y  
               << " " << v.z
               << "\n"; 
      }
      number_cells = number_cells + 1;
    }
  }

  // Teil 3: Zellen ausgeben
  // Teil 3.1: Innere Zellen ausgeben
  number_cells = 0;
  iterate_hash0 {
    I = point0->Give_Index();
    if(I.Cell_index() && I.Tiefe() == level) {
      *Datei << number_cells << "  1  hex  "
             << number_cells*8+WNTd << " "
             << number_cells*8+WSTd << " "
             << number_cells*8+ESTd << " "
             << number_cells*8+ENTd << " "
             << number_cells*8+WNDd << " "
             << number_cells*8+WSDd << " "
             << number_cells*8+ESDd << " "
             << number_cells*8+ENDd << " "
             << "\n";
      number_cells = number_cells + 1;
    }
  }
  
  // Teil 4:
  *Datei << "1 1" << "\n";
  
  // Teil 5:
  *Datei << "variable, unit " << "\n";
  
  // Teil 6: Ausgabe der Werte der Zell-Variablen
  // Teil 3.1: Innere Zellen ausgeben
  number_cells = 0;
  iterate_hash0 {
    I = point0->Give_Index();
    if(I.Cell_index() && I.Tiefe() == level) {
      
      *Datei << number_cells << "  " << Give_cell_typ(I) << "\n";
      number_cells = number_cells + 1;
    }
  }
}

---