src/expde/extemp/variable.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
//

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


// 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 "../grid/gpar.h"
#include "../grid/parallel.h"
#include "../grid/mgcoeff.h"
#include "../grid/sto_man.h"
#include "../grid/gridbase.h"
#include "../grid/grid.h"
#include "../grid/input.h"

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

// Include level 6:
#include "variable.h"


               // 1.: Maximum
               // 2.: Scalarproduct
               // 3.: print
               // 4.: calculate normal vector
               // 5.: Member functions


// 1.: Maximum


/* Binary operators */
/* product */


double L_infty(Variable& v) {
  bool run;
  int run_bo;
  int lev, color;
  double max, total_max;
  P_interior *iter_i;
  P_nearb    *iter_n;
  P_Bo2p     *iter_b;
  P_cellpoi  *iter_cf;
                         
  max=0;

  v.Test_init();  // Grid storage initialization
  lev = v.Level();
  if(lev < 0)          { cout << "\n Level too coarse! "   << endl; }
  else if (lev > v.Max_Level())  { cout << "\n Level too fine! " << endl; }
  else { // Start
    run = v.run_interior();
      if(run) {
        // number of operations
        v.Give_grid()->Add_operations(1,lev);
        // innere Punkte
        v.Give_grid()->Start_calc_time_interior();
        for(color=0;color<8;++color)
          for(iter_i=v.Give_grid()->Start_P_interior(lev,color);
              iter_i!=NULL;iter_i=iter_i->Next()) {
            if(max<ABS(v.Give_interior_simple(iter_i,v.Give_grid(),lev)))
              max = ABS(v.Give_interior_simple(iter_i,v.Give_grid(),lev));
          }
        v.Give_grid()->Stop_calc_time_interior();
      }
      run = v.run_nearb();
      v.Give_grid()->Start_calc_time_nearbb();
      if(run) {
        // Randzellfreiheitsgrade
        for(iter_cf=v.Give_grid()->Start_P_cellpoi(lev);
            iter_cf!=NULL;iter_cf=iter_cf->Next()) {
          if(max<ABS(v.Give_cellpoi(iter_cf,v.Give_grid())))
            max = ABS(v.Give_cellpoi(iter_cf,v.Give_grid()));
        }
        // randnahe Punkte
        v.Give_grid()->Start_calc_time_nearbb();
        for(color=0;color<8;++color)
          for(iter_n=v.Give_grid()->Start_P_nearb(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(max<ABS(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev)))
              max = ABS(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev));
          }
      }
      run_bo = v.run_boundary();
      if(run_bo==true) {
        // Randpunkte
        for(color=0;color<8;++color)
          for(iter_b=v.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            if(max<ABS(v.Give_Bo2p(iter_b,v.Give_grid(),lev)))
              max = ABS(v.Give_Bo2p(iter_b,v.Give_grid(),lev));
          }
        // exterior grobe Randpunkte
        for(color=0;color<8;++color)
          for(iter_n=v.Give_grid()->Start_P_exteri(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(max<ABS(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev)))
              max = ABS(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev));
          }
      }
      else if(run_bo<false) {
        // Randpunkte
        for(color=0;color<8;++color)
          for(iter_b=v.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            if(iter_b->Give_Label(-run_bo,v.Give_grid())) {
              if(max<ABS(v.Give_Bo2p(iter_b,v.Give_grid(),lev)))
                max = ABS(v.Give_Bo2p(iter_b,v.Give_grid(),lev));
            }
          }
        // exterior grobe Randpunkte
        for(color=0;color<8;++color)
          for(iter_n=v.Give_grid()->Start_P_exteri(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(iter_n->Give_Label(-run_bo,lev,v.Give_grid())) {
              if(max<ABS(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev)))
                max = ABS(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev));
            }
          }
      }
      v.Give_grid()->Stop_calc_time_nearbb();
  }
  MPI_Allreduce(&max,&total_max,1,MPI_DOUBLE,MPI_MAX,
                v.Give_grid()->Give_comm());

  return total_max;
}


double Find_one_value_of_an_active_point(Variable& v, double sign) {
  bool run;
  int run_bo;
  int lev, color;
  P_interior *iter_i;
  P_nearb    *iter_n;
  P_Bo2p     *iter_b;


  v.Test_init();  // Grid storage initialization
  lev = v.Level();
  if(lev < 0)          { cout << "\n Level too coarse! "   << endl; }
  else if (lev > v.Max_Level())  { cout << "\n Level too fine! " << endl; }
  else { // Start
    run = v.run_interior();
    if(run) {
      iter_i=v.Give_grid()->Start_P_interior(lev,0);
      if(iter_i!=NULL)
        return v.Give_interior_simple(iter_i,v.Give_grid(),lev);
    }
    run = v.run_nearb();
    if(run) {
      // randnahe Punkte
      iter_n=v.Give_grid()->Start_P_nearb(lev,0);
      if(iter_n!=NULL)
        return v.Give_nearb(iter_n,v.Give_grid(),0.0,lev);
    }
    run_bo = v.run_boundary();
    if(run_bo==true) {
      // Randpunkte
      iter_b=v.Give_grid()->Start_P_Bo2p(lev,0);
      if(iter_b!=NULL)
        return v.Give_Bo2p(iter_b,v.Give_grid(),lev);
      // exterior grobe Randpunkte
      iter_n=v.Give_grid()->Start_P_exteri(lev,0);
      if(iter_n!=NULL)
        return v.Give_nearb(iter_n,v.Give_grid(),0.0,lev);
    }
    else if(run_bo<false) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=v.Give_grid()->Start_P_Bo2p(lev,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          if(iter_b->Give_Label(-run_bo,v.Give_grid())) {
            return v.Give_Bo2p(iter_b,v.Give_grid(),lev);
          }
        }
      // exterior grobe Randpunkte
      for(color=0;color<8;++color)
        for(iter_n=v.Give_grid()->Start_P_exteri(lev,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          if(iter_n->Give_Label(-run_bo,lev,v.Give_grid())) {
            return v.Give_nearb(iter_n,v.Give_grid(),0.0,lev);
          }
        }
    }
  }
  return 1.0e99 * sign;
}


double Maximum(Variable& v) {
  bool run;
  int color;
  int run_bo;
  int lev;
  double max, total_max;
  P_interior *iter_i;
  P_nearb    *iter_n;
  P_Bo2p     *iter_b;
  P_cellpoi  *iter_cf;
                         
  max = Find_one_value_of_an_active_point(v,-1.0);

  v.Test_init();  // Grid storage initialization
  lev = v.Level();
  if(lev < 0)          { cout << "\n Level too coarse! "   << endl; }
  else if (lev > v.Max_Level())  { cout << "\n Level too fine! " << endl; }
  else { // Start
    run = v.run_interior();
      if(run) {
        // number of operations
        v.Give_grid()->Add_operations(1,lev);
        // innere Punkte
        v.Give_grid()->Start_calc_time_interior();
        for(color=0;color<8;++color)
          for(iter_i=v.Give_grid()->Start_P_interior(lev,color);
              iter_i!=NULL;iter_i=iter_i->Next()) {
            if(max<(v.Give_interior_simple(iter_i,v.Give_grid(),lev)))
              max = (v.Give_interior_simple(iter_i,v.Give_grid(),lev));
          }
        v.Give_grid()->Stop_calc_time_interior();
      }
      run = v.run_nearb();
      v.Give_grid()->Start_calc_time_nearbb();
      if(run) {
        // Randzellfreiheitsgrade
        for(iter_cf=v.Give_grid()->Start_P_cellpoi(lev);
            iter_cf!=NULL;iter_cf=iter_cf->Next()) {
          if(max<(v.Give_cellpoi(iter_cf,v.Give_grid())))
            max = (v.Give_cellpoi(iter_cf,v.Give_grid()));
        }
        // randnahe Punkte
        for(color=0;color<8;++color)
          for(iter_n=v.Give_grid()->Start_P_nearb(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(max<(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev)))
              max = (v.Give_nearb(iter_n,v.Give_grid(),0.0,lev));
          }
      }
      run_bo = v.run_boundary();
      if(run_bo==true) {
        // Randpunkte
        for(color=0;color<8;++color)
          for(iter_b=v.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            if(max<(v.Give_Bo2p(iter_b,v.Give_grid(),lev)))
              max = v.Give_Bo2p(iter_b,v.Give_grid(),lev);
          }
        // exterior grobe Randpunkte
        for(color=0;color<8;++color)
          for(iter_n=v.Give_grid()->Start_P_exteri(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(max<(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev)))
              max = (v.Give_nearb(iter_n,v.Give_grid(),0.0,lev));
          }
      }
      else if(run_bo<false) {
        // Randpunkte
        for(color=0;color<8;++color)
          for(iter_b=v.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            if(iter_b->Give_Label(-run_bo,v.Give_grid())) {
              if(max<(v.Give_Bo2p(iter_b,v.Give_grid(),lev)))
                max = v.Give_Bo2p(iter_b,v.Give_grid(),lev);
            }
          }
        // exterior grobe Randpunkte
        for(color=0;color<8;++color)
          for(iter_n=v.Give_grid()->Start_P_exteri(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(iter_n->Give_Label(-run_bo,lev,v.Give_grid())) {
              if(max<(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev)))
                max = (v.Give_nearb(iter_n,v.Give_grid(),0.0,lev));
            }
          }
      }
      v.Give_grid()->Stop_calc_time_nearbb();
  }
  MPI_Allreduce(&max,&total_max,1,MPI_DOUBLE,MPI_MAX,
                v.Give_grid()->Give_comm());

  return total_max;
}


double Minimum(Variable& v) {
  bool run;
  int color;
  int run_bo;
  int lev;
  double min, total_min;
  P_interior *iter_i;
  P_nearb    *iter_n;
  P_Bo2p     *iter_b;
  P_cellpoi  *iter_cf;
                         
  min = Find_one_value_of_an_active_point(v,1.0);

  v.Test_init();  // Grid storage initialization
  lev = v.Level();
  if(lev < 0)          { cout << "\n Level too coarse! "   << endl; }
  else if (lev > v.Max_Level())  { cout << "\n Level too fine! " << endl; }
  else { // Start
    run = v.run_interior();
      if(run) {
        // number of operations
        v.Give_grid()->Add_operations(1,lev);
        // innere Punkte
        v.Give_grid()->Start_calc_time_interior();
        for(color=0;color<8;++color)
          for(iter_i=v.Give_grid()->Start_P_interior(lev,color);
              iter_i!=NULL;iter_i=iter_i->Next()) {
            if(min>(v.Give_interior_simple(iter_i,v.Give_grid(),lev)))
              min = (v.Give_interior_simple(iter_i,v.Give_grid(),lev));
          }
        v.Give_grid()->Stop_calc_time_interior();
      }
      run = v.run_nearb();
      v.Give_grid()->Start_calc_time_nearbb();
      if(run) {
        // Randzellfreiheitsgrade
        for(iter_cf=v.Give_grid()->Start_P_cellpoi(lev);
            iter_cf!=NULL;iter_cf=iter_cf->Next()) {
          if(min>(v.Give_cellpoi(iter_cf,v.Give_grid())))
            min = (v.Give_cellpoi(iter_cf,v.Give_grid()));
        }
        // randnahe Punkte
        for(color=0;color<8;++color)
          for(iter_n=v.Give_grid()->Start_P_nearb(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(min>(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev)))
              min = (v.Give_nearb(iter_n,v.Give_grid(),0.0,lev));
          }
      }
      run_bo = v.run_boundary();
      if(run_bo==true) {
        // Randpunkte
        for(color=0;color<8;++color)
          for(iter_b=v.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            if(min>v.Give_Bo2p(iter_b,v.Give_grid(),lev))
              min = v.Give_Bo2p(iter_b,v.Give_grid(),lev);
          }
        // exterior grobe Randpunkte
        for(color=0;color<8;++color)
          for(iter_n=v.Give_grid()->Start_P_exteri(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(min>(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev)))
              min = (v.Give_nearb(iter_n,v.Give_grid(),0.0,lev));
          }
      }
      else if(run_bo<false) {
        // Randpunkte
        for(color=0;color<8;++color)
          for(iter_b=v.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            if(iter_b->Give_Label(-run_bo,v.Give_grid())) {
              if(min>v.Give_Bo2p(iter_b,v.Give_grid(),lev))
                min = v.Give_Bo2p(iter_b,v.Give_grid(),lev);
            }
          }
        // exterior grobe Randpunkte
        for(color=0;color<8;++color)
          for(iter_n=v.Give_grid()->Start_P_exteri(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(iter_n->Give_Label(-run_bo,lev,v.Give_grid())) {       
              if(min>(v.Give_nearb(iter_n,v.Give_grid(),0.0,lev)))
                min = (v.Give_nearb(iter_n,v.Give_grid(),0.0,lev));
            }
          }
      }
      v.Give_grid()->Stop_calc_time_nearbb();
  }
  MPI_Allreduce(&min,&total_min,1,MPI_DOUBLE,MPI_MIN,
                v.Give_grid()->Give_comm());

  return total_min;
}


// 2.: Scalarproduct


/* Binary operators */
/* product */


double product(Variable& b, Variable& a)
{
  bool run;
  int lev, run_bo, color;
  double sum;
  P_interior *iter_i;
  P_nearb    *iter_n;
  P_Bo2p     *iter_b;
  P_cellpoi  *iter_cf;
  double total_sum;
                         
  sum=0;

  a.Test_init();  // Grid storage initialization
  lev = a.Level();
  if(lev < 0)          { cout << "\n Level too coarse! "   << endl; }
  else if (lev > a.Max_Level())  { cout << "\n Level too fine! " << endl; }
  else { // Start
    run = a.run_interior();
      if(run) {
        // number of operations
        a.Give_grid()->Add_operations(2,lev);
        // innere Punkte
        a.Give_grid()->Start_calc_time_interior();
        for(color=0;color<8;++color)
          for(iter_i=a.Give_grid()->Start_P_interior(lev,color);
              iter_i!=NULL;iter_i=iter_i->Next()) {
            sum = sum + a.Give_interior_simple(iter_i,a.Give_grid(),lev) *
              b.Give_interior_simple(iter_i,a.Give_grid(),lev);
          }
        a.Give_grid()->Stop_calc_time_interior();
      }

bool ttt;
ttt =false; 
// ttt =true; 
if(ttt) {     
//if(a.Give_grid()->Give_my_rank()==0) 
cout << " sum interior: " << sum << " my_rank: "
     << a.Give_grid()->Give_my_rank() << endl;
}
      
      run = a.run_nearb();
      a.Give_grid()->Start_calc_time_nearbb();  
      if(run) {
        // Randzellfreiheitsgrade
        for(iter_cf=a.Give_grid()->Start_P_cellpoi(lev);
            iter_cf!=NULL;iter_cf=iter_cf->Next()) {
          sum = sum + a.Give_cellpoi(iter_cf,a.Give_grid()) *
            b.Give_cellpoi(iter_cf,a.Give_grid());
        } 
        // randnahe Punkte
        for(color=0;color<8;++color)
          for(iter_n=a.Give_grid()->Start_P_nearb(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            sum = sum + a.Give_nearb(iter_n,a.Give_grid(),0.0,lev) *
              b.Give_nearb(iter_n,a.Give_grid(),0.0,lev);
          }
      }
 
if(ttt) {     
//if(a.Give_grid()->Give_my_rank()==0)
cout << " sum nearb: " << sum << " my_rank: " 
     << a.Give_grid()->Give_my_rank() << endl;
}      

      run_bo = a.run_boundary();
      if(run_bo==true) {
        // Randpunkte
        for(color=0;color<8;++color) {
          for(iter_b=a.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            sum = sum + a.Give_Bo2p(iter_b,a.Give_grid(),lev) *
              b.Give_Bo2p(iter_b,a.Give_grid(),lev);
          }
        }
        // exterior grobe Randpunkte
        for(color=0;color<8;++color) {
          for(iter_n=a.Give_grid()->Start_P_exteri(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            sum = sum + a.Give_nearb(iter_n,a.Give_grid(),0.0,lev) *
              b.Give_nearb(iter_n,a.Give_grid(),0.0,lev);
          }
        }
      }
      else if(run_bo<false) {
        // Randpunkte
        for(color=0;color<8;++color)
          for(iter_b=a.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            if(iter_b->Give_Label(-run_bo,a.Give_grid())) {
              sum = sum + a.Give_Bo2p(iter_b,a.Give_grid(),lev) *
                b.Give_Bo2p(iter_b,a.Give_grid(),lev);
            }
          }
        // exterior grobe Randpunkte
        for(color=0;color<8;++color) {
          for(iter_n=a.Give_grid()->Start_P_exteri(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(iter_n->Give_Label(-run_bo,lev,a.Give_grid())) {
              sum = sum + a.Give_nearb(iter_n,a.Give_grid(),0.0,lev) *
                b.Give_nearb(iter_n,a.Give_grid(),0.0,lev);



              /*

cout << " sum exterior: " << sum << " my_rank: "
     << a.Give_grid()->Give_my_rank() 
     << " level: " << lev << endl;

      
              */









            }
          }
        }
      }
      a.Give_grid()->Stop_calc_time_nearbb();   

      
if(ttt) {
cout << " sum all: " << sum << " my_rank: " 
     << a.Give_grid()->Give_my_rank() << endl;
}
     

  }


  MPI_Allreduce(&sum,&total_sum,1,MPI_DOUBLE,MPI_SUM,
                a.Give_grid()->Give_comm());

  return total_sum;
}



double product(Variable& a1, Variable& a2, Variable& a3,
               Variable& b1, Variable& b2, Variable& b3)
{
  bool run; 
  int color;
  int run_bo;
  int lev;
  double sum, total_sum;
  P_interior *iter_i;
  P_nearb    *iter_n;
  P_Bo2p     *iter_b;
  P_cellpoi  *iter_cf;
                         
  sum=0;

  a1.Test_init();  // Grid storage initialization
  lev = a1.Level();
  if(lev < 0)          { cout << "\n Level too coarse! "   << endl; }
  else if (lev > a1.Max_Level())  { cout << "\n Level too fine! " << endl; }
  else { // Start
    run = a1.run_interior();
      if(run) {
        // number of operations
        a1.Give_grid()->Add_operations(6,lev);
        // innere Punkte
        a1.Give_grid()->Start_calc_time_interior();
        for(color=0;color<8;++color)
          for(iter_i=a1.Give_grid()->Start_P_interior(lev,color);
              iter_i!=NULL;iter_i=iter_i->Next()) {
            sum = sum + 
              a1.Give_interior_simple(iter_i,a1.Give_grid(),lev) *
              b1.Give_interior_simple(iter_i,a1.Give_grid(),lev) +
              a2.Give_interior_simple(iter_i,a1.Give_grid(),lev) *
              b2.Give_interior_simple(iter_i,a1.Give_grid(),lev) +
              a3.Give_interior_simple(iter_i,a1.Give_grid(),lev) *
              b3.Give_interior_simple(iter_i,a1.Give_grid(),lev);
          }
        a1.Give_grid()->Stop_calc_time_interior();
      }
      run = a1.run_nearb();
      a1.Give_grid()->Start_calc_time_nearbb();
      if(run) {
        // Randzellfreiheitsgrade
        for(iter_cf=a1.Give_grid()->Start_P_cellpoi(lev);
            iter_cf!=NULL;iter_cf=iter_cf->Next()) {
          sum = sum +
            a1.Give_cellpoi(iter_cf,a1.Give_grid()) *
            b1.Give_cellpoi(iter_cf,a1.Give_grid()) +
            a2.Give_cellpoi(iter_cf,a1.Give_grid()) *
            b2.Give_cellpoi(iter_cf,a1.Give_grid()) +
            a3.Give_cellpoi(iter_cf,a1.Give_grid()) *
            b3.Give_cellpoi(iter_cf,a1.Give_grid());
        }   
        // randnahe Punkte
        for(color=0;color<8;++color)
          for(iter_n=a1.Give_grid()->Start_P_nearb(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            sum = sum + 
              a1.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) *
              b1.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) +
              a2.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) *
              b2.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) +
              a3.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) *
              b3.Give_nearb(iter_n,a1.Give_grid(),0.0,lev);
          }
      }
      run_bo = a1.run_boundary();
      if(run_bo==true) {
        // Randpunkte
        a1.Give_grid()->Start_calc_time_nearbb();
        for(color=0;color<8;++color)
          for(iter_b=a1.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            sum = sum +
              a1.Give_Bo2p(iter_b,a1.Give_grid(),lev) *
              b1.Give_Bo2p(iter_b,a1.Give_grid(),lev) +
              a2.Give_Bo2p(iter_b,a1.Give_grid(),lev) *
              b2.Give_Bo2p(iter_b,a1.Give_grid(),lev) +
              a3.Give_Bo2p(iter_b,a1.Give_grid(),lev) *
              b3.Give_Bo2p(iter_b,a1.Give_grid(),lev);
          }
        // exterior grobe Randpunkte
        for(color=0;color<8;++color)
          for(iter_n=a1.Give_grid()->Start_P_exteri(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            sum = sum + 
              a1.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) *
              b1.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) +
              a2.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) *
              b2.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) +
              a3.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) *
              b3.Give_nearb(iter_n,a1.Give_grid(),0.0,lev);
          }
      }
      else if(run_bo<false) {
        // Randpunkte
        for(color=0;color<8;++color)
          for(iter_b=a1.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            if(iter_b->Give_Label(-run_bo,a1.Give_grid())) {
              sum = sum +
                a1.Give_Bo2p(iter_b,a1.Give_grid(),lev) *
                b1.Give_Bo2p(iter_b,a1.Give_grid(),lev) +
                a2.Give_Bo2p(iter_b,a1.Give_grid(),lev) *
                b2.Give_Bo2p(iter_b,a1.Give_grid(),lev) +
                a3.Give_Bo2p(iter_b,a1.Give_grid(),lev) *
                b3.Give_Bo2p(iter_b,a1.Give_grid(),lev);
            }
          }
        // exterior grobe Randpunkte
        for(color=0;color<8;++color)
          for(iter_n=a1.Give_grid()->Start_P_exteri(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            if(iter_n->Give_Label(-run_bo,lev,a1.Give_grid())) {
              sum = sum + 
                a1.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) *
                b1.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) +
                a2.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) *
                b2.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) +
                a3.Give_nearb(iter_n,a1.Give_grid(),0.0,lev) *
                b3.Give_nearb(iter_n,a1.Give_grid(),0.0,lev);
            }
          }
      }
      a1.Give_grid()->Stop_calc_time_nearbb();
  }
  MPI_Allreduce(&sum,&total_sum,1,MPI_DOUBLE,MPI_SUM,
                a1.Give_grid()->Give_comm());

  return total_sum;
}

// 3.: Print

void Print_UCD(Variable& a, ofstream *Datei) {
  if(parallel_version) {
    if(a.Level() == a.Give_grid()->Max_level())
      Print_UCD_parallel(a, Datei);
    else
      Print_UCD_parallel(a, Datei, a.Level());
  }
  else {
    if(a.Level() == a.Give_grid()->Max_level())
      a.Give_grid()->Print_Variable_AVS(Datei,a.Number_variable());
    else 
      a.Give_grid()->Print_Variable_AVS_coarse(Datei,a.Number_variable(),
                                               a.Level());
  }
};

void Print_UCD_parallel(Variable& a, ofstream *Datei) {
  Grid *grid;
  int array_num[1];
  array_num[0] = a.Number_variable();

  grid = a.Give_grid();
  
  grid->Update_Variable(Wdir,array_num,1,grid->Max_level());
  grid->Update_Variable(Edir,array_num,1,grid->Max_level());
  grid->Update_Variable(Ndir,array_num,1,grid->Max_level());
  grid->Update_Variable(Sdir,array_num,1,grid->Max_level());
  grid->Update_Variable(Tdir,array_num,1,grid->Max_level());
  grid->Update_Variable(Ddir,array_num,1,grid->Max_level());
  
  grid->Print_Variable_AVS_parallel(Datei,a.Number_variable());
};

void Print_UCD_parallel(Variable& a, ofstream *Datei, int level) {
  Grid *grid;
  int array_num[1];
  array_num[0] = a.Number_variable();

  grid = a.Give_grid();
  grid->Update_Variable(Wdir,array_num,1,level);
  grid->Update_Variable(Edir,array_num,1,level);
  grid->Update_Variable(Ndir,array_num,1,level);
  grid->Update_Variable(Sdir,array_num,1,level);
  grid->Update_Variable(Tdir,array_num,1,level);
  grid->Update_Variable(Ddir,array_num,1,level);

  a.Give_grid()->Print_Variable_AVS_parallel(Datei,a.Number_variable(),level);
};

void Print_UCD_surface_parallel(Variable& a, ofstream *Datei) {
  Grid *grid;
  int array_num[1];
  array_num[0] = a.Number_variable();

  grid = a.Give_grid();
  grid->Update_Variable(Wdir,array_num,1,grid->Max_level());
  grid->Update_Variable(Edir,array_num,1,grid->Max_level());
  grid->Update_Variable(Ndir,array_num,1,grid->Max_level());
  grid->Update_Variable(Sdir,array_num,1,grid->Max_level());
  grid->Update_Variable(Tdir,array_num,1,grid->Max_level());
  grid->Update_Variable(Ddir,array_num,1,grid->Max_level());

  a.Give_grid()->Print_surface_Variable_AVS_parallel(Datei,a.Number_variable());
};

void Print_UCD(Variable& a, Variable& b, Variable& c, ofstream *Datei) {
  if(parallel_version) {
    Print_UCD_parallel(a,b,c,Datei);
  }
  else {
    a.Give_grid()->Print_Variable_AVS(Datei,
                                      a.Number_variable(),
                                      b.Number_variable(),
                                      c.Number_variable());
  }
};

void Print_UCD_parallel(Variable& a, Variable& b, Variable& c,
                        ofstream *Datei) {
  Grid *grid;
  int array_num[3];
  array_num[0] = a.Number_variable();
  array_num[1] = b.Number_variable();
  array_num[2] = c.Number_variable();

  grid = a.Give_grid();
  grid->Update_Variable(Wdir,array_num,3,grid->Max_level());
  grid->Update_Variable(Edir,array_num,3,grid->Max_level());
  grid->Update_Variable(Ndir,array_num,3,grid->Max_level());
  grid->Update_Variable(Sdir,array_num,3,grid->Max_level());
  grid->Update_Variable(Tdir,array_num,3,grid->Max_level());
  grid->Update_Variable(Ddir,array_num,3,grid->Max_level());

  a.Give_grid()->Print_Variable_AVS_parallel(Datei,
                                             a.Number_variable(),
                                             b.Number_variable(),
                                             c.Number_variable());
};

void Print_UCD_moved(Variable& v,
                     Variable& a, Variable& b, Variable& c, ofstream *Datei) {
  if(parallel_version) {
    Print_UCD_moved_parallel(v, a, b, c, Datei);
  }
  else {
    a.Give_grid()->Print_Variable_AVS_moved(Datei,
                                            v.Number_variable(),
                                            a.Number_variable(),
                                            b.Number_variable(),
                                            c.Number_variable());
  }
}

void Print_UCD_moved_parallel(Variable& v,
                     Variable& a, Variable& b, Variable& c, ofstream *Datei) {
  Grid *grid;
  int array_num[4];
  array_num[0] = a.Number_variable();
  array_num[1] = b.Number_variable();
  array_num[2] = c.Number_variable();
  array_num[3] = v.Number_variable();

  grid = v.Give_grid();
  grid->Update_Variable(Wdir,array_num,4,grid->Max_level());
  grid->Update_Variable(Edir,array_num,4,grid->Max_level());
  grid->Update_Variable(Ndir,array_num,4,grid->Max_level());
  grid->Update_Variable(Sdir,array_num,4,grid->Max_level());
  grid->Update_Variable(Tdir,array_num,4,grid->Max_level());
  grid->Update_Variable(Ddir,array_num,4,grid->Max_level());

  a.Give_grid()->Print_Variable_AVS_moved_parallel(Datei,
                                                   v.Number_variable(),
                                                   a.Number_variable(),
                                                   b.Number_variable(),
                                                   c.Number_variable());
}

void Print_Gnuplot_moved(Variable& a, Variable& b, Variable& c,
                             ofstream *Datei) {
  a.Give_grid()->Print_Grid_Gnuplot_moved(Datei,
                                          a.Number_variable(),
                                          b.Number_variable(),
                                          c.Number_variable());
}


void Print_Dx_parallel(Variable& a, ofstream *Datei) {
  a.Give_grid()->Print_Variable_OpenDx_parallel(Datei,a.Number_variable());
};


void Print_Dx(Variable& a, ofstream *Datei) {
  if(parallel_version) {
    Print_Dx_parallel(a,Datei);
  }
  else {
    a.Give_grid()->Print_Variable_OpenDx(Datei,a.Number_variable());
  }
};


void Print_Dx_moved(Variable& v,
                     Variable& a, Variable& b, Variable& c, ofstream *Datei) {
  if(parallel_version) {
    cout << " Error: No parallel version for Print_Dx_moved !" << endl;
  }
  a.Give_grid()->Print_Variable_OpenDx_moved(Datei,
                                          v.Number_variable(),
                                          a.Number_variable(),
                                          b.Number_variable(),
                                          c.Number_variable());
}

void Print(Variable& v) {
  bool run; 
  int color;
  int run_bo;
  int lev;
  P_interior *iter_i;
  P_nearb    *iter_n;
  P_Bo2p     *iter_b;
  P_cellpoi  *iter_cf;
  P_parallel *iter_pa;

  int my_rank;
   MPI_Comm_rank(v.Give_grid()->Give_comm(), &my_rank);
                         
  v.Test_init();  // Grid storage initialization
  lev = v.Level();
  if(lev < 0)          { cout << "\n Level too coarse! "   << endl; }
  else if (lev > v.Max_Level())  { cout << "\n Level too fine! " << endl; }
  else { // Start
    run = v.run_interior();
      if(run) {
        MPI_Barrier(v.Give_grid()->Give_comm());
        if(my_rank==0) cout << " interior points: " << endl;
        MPI_Barrier(v.Give_grid()->Give_comm());
        // innere Punkte
        v.Give_grid()->Start_calc_time_interior();
        for(color=0;color<8;++color)
          for(iter_i=v.Give_grid()->Start_P_interior(lev,color);
              iter_i!=NULL;iter_i=iter_i->Next()) {
          cout << "at point: "; 
          iter_i->coordinate(v.Give_grid()).Print();
          cout << v.Give_interior_simple(iter_i,v.Give_grid(),lev) << endl;
        }
        MPI_Barrier(v.Give_grid()->Give_comm());
        if(my_rank==0) cout << " points in a boundary cell: " << endl;
        MPI_Barrier(v.Give_grid()->Give_comm());
        // Randzellfreiheitsgrade
        v.Give_grid()->Start_calc_time_nearbb();
        for(iter_cf=v.Give_grid()->Start_P_cellpoi(lev);
            iter_cf!=NULL;iter_cf=iter_cf->Next()) {
          cout << "at point: "; 
          iter_cf->coordinate(v.Give_grid()).Print();
          cout << v.Give_cellpoi(iter_cf,v.Give_grid()) << endl;
        }
      }
      run = v.run_nearb();
      if(run) {
        MPI_Barrier(v.Give_grid()->Give_comm());
        if(my_rank==0) cout << " points near the boundary: " << endl;
        MPI_Barrier(v.Give_grid()->Give_comm());
        // randnahe Punkte
        for(color=0;color<8;++color)
          for(iter_n=v.Give_grid()->Start_P_nearb(lev,color);
              iter_n!=NULL;iter_n=iter_n->Next()) {
            cout << "at point: "; 
            iter_n->coordinate(v.Give_grid()).Print();
            cout << v.Give_nearb(iter_n,v.Give_grid(),0.0,lev) << endl;
          }
      }
      run_bo = v.run_boundary();
      MPI_Barrier(v.Give_grid()->Give_comm());
      if(my_rank==0) cout << " boundary points: " << endl;
      MPI_Barrier(v.Give_grid()->Give_comm());
      if(run_bo==true) {
        // Randpunkte
        for(color=0;color<8;++color)
          for(iter_b=v.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            cout << "at point: "; 
            iter_b->coordinate(v.Give_grid()).Print();
            cout << v.Give_Bo2p(iter_b,v.Give_grid(),lev) << endl;
          }
      }
      else if(run_bo<false) {
        // Randpunkte
        for(color=0;color<8;++color)
          for(iter_b=v.Give_grid()->Start_P_Bo2p(lev,color);
              iter_b!=NULL;iter_b=iter_b->Next()) {
            if(iter_b->Give_Label(-run_bo,v.Give_grid())) {
              cout << "at point: "; 
              iter_b->coordinate(v.Give_grid()).Print();
              cout << v.Give_Bo2p(iter_b,v.Give_grid(),lev) << endl;
            }
          }
      }

      // Parallel
      MPI_Barrier(v.Give_grid()->Give_comm());
      if(my_rank==0) cout << " parallel points: " << endl;
      MPI_Barrier(v.Give_grid()->Give_comm());
      for(int r=0;r<v.Give_grid()->give_number_of_processes();++r) {
        if(r==my_rank) {
          cout << " \n rank: " << r << endl;
          for(iter_pa=v.Give_grid()->Start_P_parallel(lev);
              iter_pa!=NULL;iter_pa=iter_pa->Next()) {
            cout << "at point: "; 
            iter_pa->coordinate(v.Give_grid()).Print();
            cout << iter_pa->varM(v.Give_grid(),lev)[v.Number_variable()] 
                 << endl;
          }
        }
        MPI_Barrier(v.Give_grid()->Give_comm());
      }
  }
}


// 4.: calculate normal vector

void Normal_vector(Variable& vx, Variable& vy, Variable& vz) {
  double* h;
  Grid *grid;
  int i, num, lev, color;
  P_Bo2p *iter_b;
  BoCeData* bocedata;
  Tetraeder_storage* tet;
  D3vector normal;
  double weight;
  Index3D I;

  grid = vx.Give_grid();
  h    = new double[30];
  lev = grid->Max_level();

  // Randpunkte
  for(color=0;color<8;++color)
    for(iter_b=grid->Start_P_Bo2p(lev,color);
        iter_b!=NULL;iter_b=iter_b->Next()) {
      weight = 0.0;
      normal = D3vector(0.0,0.0,0.0);
      for(i=0;i<8;++i) {    // Durchlaufen aller Nachbarzellen
        // Info ueber i-Randzelle
        bocedata = iter_b->Give_BoData((dir_sons)i);
        I = iter_b->Ind().next((dir_sons)i,grid->Max_level()+1);
        
        if(bocedata!=NULL) {
          // Calc h
          for(num=0;num<bocedata->Give_number_points();++num) {
            if(bocedata->edge_point(num)==edge_poi_typ) {
              h[num] = grid->Give_h(I.neighbour(bocedata->corner(num)),
                                    bocedata->edge_dir(num));
            }
          }
          for(tet=bocedata->Give_boundary_tets();tet!=NULL;tet=tet->Next()) {
            normal = normal + 
              (normal_vector_of_triangle(bocedata->coord(tet->N0(),h),
                                         bocedata->coord(tet->N2(),h),
                                         bocedata->coord(tet->N1(),h))
               * tet->Give_det_surface());
          }
        }
      }
      weight = norm(normal);
      
      iter_b->varM(grid)[vx.Number_variable()] = normal.x / weight;
      iter_b->varM(grid)[vy.Number_variable()] = normal.y / weight;
      iter_b->varM(grid)[vz.Number_variable()] = normal.z / weight;
    }

  vx.Active_boundary(true);
  vx.Active_interior(false);
  vx.Active_nearb(false);
  vx.Active_Level(grid->Max_level());   

  vy.Active_boundary(true);
  vy.Active_interior(false);
  vy.Active_nearb(false);
  vy.Active_Level(grid->Max_level());   

  vz.Active_boundary(true);
  vz.Active_interior(false);
  vz.Active_nearb(false);
  vz.Active_Level(grid->Max_level());   

  delete h;
}



// 5.: member functions


Variable::Variable(Grid* gridp) {
  grid = NULL;
  Initialize(gridp);
}

Variable::Variable() {
  grid = NULL;
  number_variable = 0;
  level = 0;
  r_interior = true;
  r_boundary = true;
  r_nearb = true;
}

void Variable::Initialize(Grid* gridp) {
  if(grid!=NULL) 
    cout << " Error in Initialize. Variable is initialized! " << endl;
  grid = gridp;
  number_variable = grid->Give_number_variable();
  // all points of finest level become active
  level = gridp->Max_level();
  r_interior = true;
  r_boundary = true;
  r_nearb = true;
  if(gridp->Is_storage_initialized()==true) {
    cout << " \n Definition of variable not possible! ";
    cout << " \n Define Res_stencil_boundary before ";
    cout << " initialization of the grid! " << endl;
  }
}


void Variable::operator=(const Variable &v) {
  bool run; 
  int color;
  int  run_bo;
  int lev;
  double *vars;

  Test_init();  // Grid storage initialization
  lev = v.Level();
  if(lev < 0)          { cout << "\n Level too coarse! "   << endl; }
  else if (lev > Max_Level())  { cout << "\n Level too fine! " << endl; }

  if(grid->I_am_active()) {
    // for update type
    if(parallel_version) {
      grid->Set_type_of_update(number_variable,lev,no_update);
    }
    
    run = v.run_interior();
    if(run) {
      // innere Punkte
      for(color=0;color<8;++color)
        for(iter_i=grid->Start_P_interior(lev,color);
            iter_i!=NULL;iter_i=iter_i->Next()) {
          vars = iter_i->varM(grid,lev);
          vars[number_variable] = vars[v.Number_variable()];
        }
    }
    Active_interior(run);
    run = v.run_nearb();
    if(run) {
      // Randzellfreiheitsgrade
      for(iter_cf=grid->Start_P_cellpoi(lev);
          iter_cf!=NULL;iter_cf=iter_cf->Next()) {
        vars = iter_cf->varM(grid);
        vars[number_variable] = vars[v.Number_variable()];
      }
      // randnahe Punkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_nearb(lev,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          vars = iter_n->varM(grid,lev);
          vars[number_variable] = vars[v.Number_variable()];
        } 
    }
    Active_nearb(run);
    run_bo = v.run_boundary();
    if(run_bo==true) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=grid->Start_P_Bo2p(lev,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          vars = iter_b->varM(grid);
          vars[number_variable] = vars[v.Number_variable()];
        }
      // exterior grobe Randpunkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_exteri(lev,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          vars = iter_n->varM(grid,lev);
          vars[number_variable] = vars[v.Number_variable()];
        } 
    }
    else if(run_bo<false) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=grid->Start_P_Bo2p(lev,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          if(iter_b->Give_Label(-run_bo,grid)) {      
            vars = iter_b->varM(grid);
            vars[number_variable] = vars[v.Number_variable()];
          }
        }
      // exterior grobe Randpunkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_exteri(lev,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          if(iter_n->Give_Label(-run_bo,lev,v.Give_grid())) {
            vars = iter_n->varM(grid,lev);
            vars[number_variable] = vars[v.Number_variable()];
          }
        }
    }
    Active_boundary(run_bo);
    Active_Level(lev);   
  }
}




void Variable::operator=(double (*Formula)(double x,double y,double z)) {
  D3vector V;
  int color;
  Test_init();  // Grid storage initialization

  if(grid->I_am_active()) {
    // for update type
    if(parallel_version) {
      grid->Set_type_of_update(number_variable,level,no_update);
    }
    
    if(r_interior) {
      // innere Punkte
      for(color=0;color<8;++color)
        for(iter_i=grid->Start_P_interior(level,color);
            iter_i!=NULL;iter_i=iter_i->Next()) {
          V = iter_i->coordinate(grid);
          iter_i->varM(grid,level)[number_variable] = 
            Formula(V.x,V.y,V.z);
        }
    }
    if(r_nearb) {
      // Randzellfreiheitsgrade
      for(iter_cf=grid->Start_P_cellpoi(level);
          iter_cf!=NULL;iter_cf=iter_cf->Next()) {
        V = iter_cf->coordinate(grid);
        iter_cf->varM(grid)[number_variable] = 
          Formula(V.x,V.y,V.z);
      }
      // randnahe Punkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_nearb(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          V = iter_n->coordinate(grid);
          iter_n->varM(grid,level)[number_variable] = 
            Formula(V.x,V.y,V.z);
        }
    }
    
    if(r_boundary==true) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=grid->Start_P_Bo2p(level,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          V = iter_b->coordinate(grid);
          iter_b->varM(grid)[number_variable] = 
            Formula(V.x,V.y,V.z);
        }
      // exterior grobe Randpunkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_exteri(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          V = iter_n->coordinate(grid);
          iter_n->varM(grid,level)[number_variable] = 
            Formula(V.x,V.y,V.z);
        }
    }
    else if(r_boundary < false) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=grid->Start_P_Bo2p(level,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          if(iter_b->Give_Label(-r_boundary,grid)) {      
            V = iter_b->coordinate(grid);
            iter_b->varM(grid)[number_variable] = 
              Formula(V.x,V.y,V.z);
          }
        }
      // exterior grobe Randpunkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_exteri(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          if(iter_n->Give_Label(-r_boundary,level,grid)) {
            V = iter_n->coordinate(grid);
            iter_n->varM(grid,level)[number_variable] = 
              Formula(V.x,V.y,V.z);
          }
        }
    }
  }
}


void Variable::operator=(Input_data_object& input_object) {
  D3vector V;
  int color;
  Test_init();  // Grid storage initialization

  if(grid->I_am_active()) {
    // for update type
    if(parallel_version) {
      grid->Set_type_of_update(number_variable,level,no_update);
    }
    
    if(r_interior) {
      // innere Punkte
      for(color=0;color<8;++color)
        for(iter_i=grid->Start_P_interior(level,color);
            iter_i!=NULL;iter_i=iter_i->Next()) {
          V = iter_i->coordinate(grid);
          iter_i->varM(grid,level)[number_variable] = 
            input_object.Interpolate_value(V);
        }
    }
    if(r_nearb) {
      // Randzellfreiheitsgrade
      for(iter_cf=grid->Start_P_cellpoi(level);
          iter_cf!=NULL;iter_cf=iter_cf->Next()) {
        V = iter_cf->coordinate(grid);
        iter_cf->varM(grid)[number_variable] = 
          input_object.Interpolate_value(V);
      }
      // randnahe Punkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_nearb(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          V = iter_n->coordinate(grid);
          iter_n->varM(grid,level)[number_variable] = 
            input_object.Interpolate_value(V);
        }
    }
    
    if(r_boundary==true) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=grid->Start_P_Bo2p(level,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          V = iter_b->coordinate(grid);
          iter_b->varM(grid)[number_variable] = 
            input_object.Interpolate_value(V);
        }
      // exterior grobe Randpunkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_exteri(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          V = iter_n->coordinate(grid);
          iter_n->varM(grid,level)[number_variable] = 
            input_object.Interpolate_value(V);
        }
    }
    else if(r_boundary < false) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=grid->Start_P_Bo2p(level,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          if(iter_b->Give_Label(-r_boundary,grid)) {      
            V = iter_b->coordinate(grid);
            iter_b->varM(grid)[number_variable] = 
              input_object.Interpolate_value(V);
          }
        }
      // exterior grobe Randpunkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_exteri(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          if(iter_n->Give_Label(-r_boundary,level,grid)) {
            V = iter_n->coordinate(grid);
            iter_n->varM(grid,level)[number_variable] = 
              input_object.Interpolate_value(V);
          }
        }
    }
  }
}




void Variable::operator=(const FunctionClass& function) {
  D3vector V;
  int run_bo;
  int color;
  Test_init();  // Grid storage initialization

  if(grid->I_am_active()) {
    // for update type
    if(parallel_version) {
      grid->Set_type_of_update(number_variable,level,no_update);
    }
    
    if(function.run_interior()) {
      // innere Punkte
      for(color=0;color<8;++color)
        for(iter_i=grid->Start_P_interior(level,color);
            iter_i!=NULL;iter_i=iter_i->Next()) {
          V = iter_i->coordinate(grid);
          iter_i->varM(grid,level)[number_variable] = 
            function.Formula(V.x,V.y,V.z);
        }
    }
    Active_interior(function.run_interior());
    if(function.run_nearb()) {
      // Randzellfreiheitsgrade
      for(iter_cf=grid->Start_P_cellpoi(level);
          iter_cf!=NULL;iter_cf=iter_cf->Next()) {
        V = iter_cf->coordinate(grid);
        iter_cf->varM(grid)[number_variable] = 
          function.Formula(V.x,V.y,V.z);
      }
      // randnahe Punkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_nearb(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          V = iter_n->coordinate(grid);
          iter_n->varM(grid,level)[number_variable] = 
            function.Formula(V.x,V.y,V.z);
        }
    }
    Active_nearb(function.run_nearb());
    
    run_bo = function.run_boundary();
    if(run_bo==true) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=grid->Start_P_Bo2p(level,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          V = iter_b->coordinate(grid);
          iter_b->varM(grid)[number_variable] = 
            function.Formula(V.x,V.y,V.z);
        }
      // exterior grobe Randpunkte
       for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_exteri(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          V = iter_n->coordinate(grid);
          iter_n->varM(grid,level)[number_variable] = 
            function.Formula(V.x,V.y,V.z);
        }
    }
    else if(run_bo<false) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=grid->Start_P_Bo2p(level,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          if(iter_b->Give_Label(-run_bo,grid)) {     
            V = iter_b->coordinate(grid);
            iter_b->varM(grid)[number_variable] = 
              function.Formula(V.x,V.y,V.z);
          }
        }
      // exterior grobe Randpunkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_exteri(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          if(iter_n->Give_Label(-run_bo,level,grid)) {
            V = iter_n->coordinate(grid);
            iter_n->varM(grid,level)[number_variable] = 
              function.Formula(V.x,V.y,V.z);
          }
        }
    }
    Active_boundary(function.run_boundary());
  }
}

void Variable::operator= (double x) {
  Test_init();  // Grid storage initialization
  int color;

  if(grid->I_am_active()) {
    // for update type
    if(parallel_version) {
      grid->Set_type_of_update(number_variable,level,no_update);
      //      grid->Set_type_of_update(number_variable,level,full_update);
    }

    if(r_interior) {
      // innere Punkte
      for(color=0;color<8;++color)
        // for(color=scco;color==scco;++color)
        for(iter_i=grid->Start_P_interior(level,color);
            iter_i!=NULL;iter_i=iter_i->Next()) {
          iter_i->varM(grid,level)[number_variable] = x;
        }
    }
    if(r_nearb) {
      // Randzellfreiheitsgrade
      for(iter_cf=grid->Start_P_cellpoi(level);
          iter_cf!=NULL;iter_cf=iter_cf->Next()) {
        iter_cf->varM(grid)[number_variable] = x;
      }
      // randnahe Punkte
      for(color=0;color<8;++color)
        // for(color=scco;color==scco;++color)
        for(iter_n=grid->Start_P_nearb(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          iter_n->varM(grid,level)[number_variable] = x;
        }
    }
    if(r_boundary==true) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=grid->Start_P_Bo2p(level,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          iter_b->varM(grid)[number_variable] = x;
        }
      // exterior grobe Randpunkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_exteri(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          iter_n->varM(grid,level)[number_variable] = x;
        }
    }
    else if(r_boundary<false) {
      // Randpunkte
      for(color=0;color<8;++color)
        for(iter_b=grid->Start_P_Bo2p(level,color);
            iter_b!=NULL;iter_b=iter_b->Next()) {
          if(iter_b->Give_Label(-r_boundary,grid)) {     
            iter_b->varM(grid)[number_variable] = x;
          }
        }
      // exterior grobe Randpunkte
      for(color=0;color<8;++color)
        for(iter_n=grid->Start_P_exteri(level,color);
            iter_n!=NULL;iter_n=iter_n->Next()) {
          if(iter_n->Give_Label(-r_boundary,level,grid)) {
            iter_n->varM(grid,level)[number_variable] = x;
          }
        }
    }
  }
}
//  Parallel Ghost nodes

void Sum_ghost_nodes(Variable& v) {
  if(parallel_version) {
    Grid* grid;
    grid = v.Give_grid();
    int array_num[1];
    array_num[0] = v.Number_variable();
    grid->Sum_ghost_nodes_Variable(array_num,1,grid->Max_level());
  }
}

void Update_ghost_nodes(Variable& v) {
  if(parallel_version) {
    Grid* grid;
    grid = v.Give_grid();
    int array_num[1];
    array_num[0] = v.Number_variable();
    
    for(int d=0;d<26;++d)
      grid->Update_Variable(d,array_num,1,grid->Max_level());
  }
}

---