src/expde/extemp/array.h

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

// ------------------------------------------------------------
//
//     array.h 
//
// ------------------------------------------------------------

#ifndef ARRAY_H_
#define ARRAY_H_


//    1. arrays of Variables
//    2. arrays of Local Variables
//    3. arrays of Restriction operators for stencils
//    4. sums and for iteration
//    5. procedure


//    1. arrays of Variables

class DExprVAR_ARR {
 private:
  // Nummer der Variablen
  int number_variable;
  int *value_;
  Grid* grid;
  Variable* start;
  int dim;
 public:
  DExprVAR_ARR(Variable* vec, int number_var, int* v, int dimension) {
    dim = dimension;
    start = vec;
    number_variable = number_var;
    value_ = v;
    grid = vec[0].Give_grid();
  }
  // special for array
  Grid* Give_grid() const {
    return grid;
  }
  int Number_variable_start() const {
    return number_variable;
  } 
  int* Give_value() const {
    return value_;
  }
  int Give_dim() const {
    return dim;
  }
  // end special
  int Number_variable() const {
    return number_variable + (*value_);
  }
  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const {
    return it_i->varM(gr,l)[number_variable+*value_];
  };
  double Give_interior_coarse(const P_interior* it_i,const Grid* gr,
                              double h_mesh, int l, double_stencils_in) const {
    return it_i->varM(gr,l)[number_variable+*value_];
  };
  double Give_nearb(const P_nearb* it_n,const Grid* gr, double h_mesh, int l,
                    const Nearb_Ablage* nearb_ablage) 
    const {
    return it_n->varM(gr,l)[number_variable+*value_];
  };
  double Give_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    return it_b->varM(gr)[number_variable+*value_];
  };  
  double Give_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                      const BoCeData* bocedata) const {
    return it_cf->varM(gr)[number_variable+*value_];
  };

  // Size of necessary stencil
  int Sice_stencil() const { return 1; }

  // activity of points
  int Level() const {
    if(*value_>=dim) *value_=0;
    return start[*value_].Level(); 
  }
  Dominace_label Dominant_lev() const { return not_dominant; }
  Dominace_label Dominant_poi() const { return not_dominant; }
  bool run_interior()  const { 
    if(*value_>=dim) *value_=0;
    return start[*value_].run_interior(); 
  }
  bool run_nearb()     const { 
    if(*value_>=dim) *value_=0;
    return start[*value_].run_nearb(); 
  }
  int run_boundary()  const {
    if(*value_>=dim) *value_=0;
    return start[*value_].run_boundary(); 
  }

  // number of operations
  int ops_interior() const { return 0; }

  // For operator ==
  void Active_Sim_Level(    int lev) const {  };
  void Active_Sim_interior(bool run) const {  };
  void Active_Sim_nearb(bool run) const {  };
  void Active_Sim_boundary( int run) const {  };
  void Active_Sim_update(Evaluation_Parallelization_object* evpar,
                         int level, type_of_update typ) const {  };

  // Put pointer to grid and activity of boundary points of test space  
  void Put_grid_rbo(Grid* gr, int r_bo)  const { };
  // For abstract differential operators:
  // --------------------------------------
  // Is the expression an abstract differential operator?
  Differential_op_typ Abstract_differential_operator()
    const { return abstract_both; }
  int Give_number_var_of_abstract_op() 
    const { return number_variable; }
  bool Give_array_variable_inserted() const {
    return true;
  }
  int Give_length_of_array_variable_inserted() const {
    return dim;
  }
  // For restriction of stencils
  // --------------------------------------
  double Give_interior_sten_element(Index3D I, double meshsize,
                                    int stelle,
                                    const Grid* grid, int l) const {
    cout << " Error in  DExprVAR_ARR! " << endl;
   return (grid->Give_variable(I.neighbour(ESTd),l)[number_variable+*value_] +
           grid->Give_variable(I.neighbour(WNDd),l)[number_variable+*value_])
      * 0.5;
  }
  double Give_boundary_sten_element(const Grid* grid, BoCeData* b_cell,
                                    double* u_ablage, int num_v) const {
    cout << " Error in  DExprVAR_ARR! " << endl;
    return b_cell->vars[num_v][number_variable+*value_];
    //    return 1.0;
  }
  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {
    for(int i=0;i<dim;++i)
      start[i].Add_variables_for_parallel(evpar);
  };
  bool GS_type(int var_number_left, int dim) const {
  return false; }; 

  // at the and of an iteration of an expression
  void clean() const { };
};


class Array_Variable {
 public:
  Array_Variable(int dimension, Grid* gridp) {
    dim  = dimension;
    start = new Variable[dim];
    for(int i=0;i<dim;++i) {
      start[i].Initialize(gridp);
    }
    number_variable = start->Number_variable();
  }
  Variable& operator[](int i) {
    return start[i];
  }
  DExpr<DExprVAR_ARR>
  operator[](Local_int& I) const {
    return DExpr<DExprVAR_ARR>(DExprVAR_ARR(start,
                                            start[0].Number_variable(),
                                            I.Give_pointer(),
                                            dim));
  }
  int Give_dim() {
    return dim;
  }

  /*
  // for procedure
  double* Give_regular_pointer(const P_interior* it_i,const Grid* gr, int l) 
    const {
    return &(it_i->varM(gr,l)[number_variable]);
  };
  double* Give_regular_pointer(const P_nearb* it_n,const Grid* gr, int l) 
    const {
    return &(it_n->varM(gr,l)[number_variable]);
  };
  double* Give_Bo2p_pointer(const P_Bo2p* it_b,const Grid* gr, int l) const {
    return &(it_b->varM(gr)[number_variable]);
  };
  double* Give_cellpoi_pointer(const P_cellpoi* it_cf,const Grid* gr) const {
    return &(it_cf->varM(gr)[number_variable]);
  };
  */

 private:
  Variable *start;
  int dim;
  int number_variable;  // not needed
};

Macro_operator_definition(DExpr<DExprVAR_ARR>);
Macro_operator_var_definition(DExpr<DExprVAR_ARR>,DExpr<DExprVAR_ARR>);
//Macro_operator_var_definition(DExpr<DExprVAR_ARR>,Variable);


//    2. arrays of local variables

class DExprLoc_VAR_ARR {
 private:
  double* value_;
  int *Index_value_;

  int dim;
 public:
  DExprLoc_VAR_ARR(double* var_value_, int* v, int dimension) {
    value_ = var_value_;
    Index_value_ = v;
    dim = dimension;
  }
  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const {
    return value_[*Index_value_];
  }
  double Give_interior_coarse(const P_interior* it_i,const Grid* gr,
                              double h_mesh,
                              int l, double_stencils_in) const {
    return value_[*Index_value_];
  }
  double Give_nearb(const P_nearb* it_n,const Grid* gr, double h_mesh, int l,
                    const Nearb_Ablage* nearb_ablage)
    const {
    return value_[*Index_value_];
  };
  double Give_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    return value_[*Index_value_];
  };
  double Give_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                      const BoCeData* bocedata) const {
    return value_[*Index_value_];
  };
  // Size of necessary stencil
  int Sice_stencil() const { return 1; }

  // activity of points
  int Level() const { return 0; }
  Dominace_label Dominant_lev() const { return anti_dominant;  }
  Dominace_label Dominant_poi() const { return anti_dominant;  }
  bool run_interior() const { return true; }
  bool run_nearb()    const { return true; }
  int run_boundary()  const { return 1; }
  
  // number of operations
  int ops_interior() const { return 0; }
  // For operator ==
  void Active_Sim_Level( int lev) const {  };
  void Active_Sim_interior(bool run) const {  };
  void Active_Sim_nearb(bool run) const {  };
  void Active_Sim_boundary( int run) const {  };
  void Active_Sim_update(Evaluation_Parallelization_object* evpar,
                         int level, type_of_update typ) const {  };

  // Put pointer to grid and activity of boundary points of test space  
  void Put_grid_rbo(Grid* gr, int r_bo) const { };
  // Is the expression an abstract differential operator?
  Differential_op_typ Abstract_differential_operator() const
    { return abstract_both; }
  int Give_number_var_of_abstract_op() const { return -1; }
  bool Give_array_variable_inserted() const {
    return false;
  }
  int Give_length_of_array_variable_inserted() const {
    return 0;
  } 
  // For restriction of stencils
  // --------------------------------------
  double Give_interior_sten_element(Index3D I, double meshsize,
                                    int stelle,
                                    const Grid* grid, int l) const {
    if(*value_>=dim) *value_=0;
    return value_[*Index_value_];
  }
  double Give_boundary_sten_element(const Grid* grid, BoCeData* b_cell,
                                    double* u_ablage, int num_v) const {
    if(*value_>=dim) *value_=0;
    return value_[*Index_value_];
  }
  // Get and Put:
  // Set value
  void Put_interior(const P_interior* it_i, int lev, double x) { Put(x); };
  void Put_nearb(const P_nearb* it_n, int lev, double x) { Put(x); };
  void Put_Bo2p(const P_Bo2p* it_b, double x) { Put(x); };
  void Put_cellpoi(const P_cellpoi* it_cf, double x) { Put(x); };

  void Put(double x) {
    value_[*Index_value_]=x;
  };
  double Get() const {
    return value_[*Index_value_];
  };          
  double* Give_pointer() const {
    return value_;
  };
  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {};
  bool GS_type(int var_number_left, int dim) const {
  return false; }; 

  // at the and of an iteration of an expression
  void clean() const {};
};

class Array_Local_var {
 public:
  Array_Local_var(int dimension, double *startD) {
    dim = dimension;
    start_double = startD;
    start = new Local_var[dim];
    for(int i=0;i<dim;++i) {
      start[i].Initialize(&(startD[i]));
    }
  }
  Local_var& operator[](int i) const {
    return start[i];
  }
  DExpr<DExprLoc_VAR_ARR>
  operator[](Local_int& I) const {
    return DExpr<DExprLoc_VAR_ARR>(DExprLoc_VAR_ARR(start[0].Give_pointer(),
                                                    I.Give_pointer(),dim));
  }
  int Give_dim() {
    return dim;
  }

  /*
  // for procedure
  double* Give_regular_pointer(const P_interior* it_i,const Grid* gr, int l) 
    const {
    return start_double;
  };
  double* Give_regular_pointer(const P_nearb* it_n,const Grid* gr, int l) 
    const {
    return start_double;
  };
  double* Give_Bo2p_pointer(const P_Bo2p* it_b,const Grid* gr, int l) const {
    return start_double;
  };
  double* Give_cellpoi_pointer(const P_cellpoi* it_cf,const Grid* gr) const {
    return start_double;
  };
  */

 private:
   Local_var *start;
   int dim;
   double* start_double;  // not needed ?
}; 



//    3. arrays of Restriction operators for stencils


template<class A>
class DResDiff_Bo_ARR {
  A a_;

  //number of storage for boundary stencils
  int num_stencils;
  int *Index_value_;
  bool array_variable_inserted; // was an array inserted or
                                // only one variable
  int dim;                      // length of array if array was inserted

  // information about the variable put in the abstract diff. operator
  int number_variable;  // number of the variable
  bool* label_on_cell;

  // Gitter
  Grid *grid;

  // needed for the calculation of local stiffness matrix
  double** Stencils_fine;  // eight stencils for finer cells
  double* u_ablage; // Ablage

  bool *corner_in_domain;

 public:
  // Constructor for concrete operators
  DResDiff_Bo_ARR(const A& a, int num_s, Grid* gr, int *Ivalue_) 
    : a_(a), num_stencils(num_s), grid(gr), Index_value_(Ivalue_) {
    number_variable = a_.Give_number_var_of_abstract_op();
    array_variable_inserted = a_.Give_array_variable_inserted();
    dim = a_.Give_length_of_array_variable_inserted();
  };
  double Give_interior(const P_interior* it_i, const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const {
    return a_.Give_interior(it_i,gr,h_mesh,l, double_stencils_out);
  };
  double Give_interior_coarse(const P_interior* it_i, const Grid* gr,
                              double h_mesh,
                              int l, double_stencils_in) const {
    return a_.Give_interior_coarse(it_i,gr,h_mesh,l, double_stencils_out);
  };
  double Give_nearb(const P_nearb* it_n, const Grid* gr, double h_mesh, int l,
                    const Nearb_Ablage* nearb_ablage) const {
    static int i;
    static Celltype cetyp;
    static double sum;
    dir_sons dir_v;
    static int actual_number; 

    // Part 1: interior cells and boundary cells on the finest grid
    sum = a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage);

    if(array_variable_inserted)
      actual_number = number_variable+(*Index_value_);
    else
      actual_number = number_variable;
    // Part 2: boundary cells on coarser grid
    for(i=0;i<8;++i) {    // Durchlaufen aller Nachbarzellen
      dir_v = opposite3D((dir_sons)i);
      // Hole Info ueber i-Randzelle
      cetyp = nearb_ablage->Give_Celltyp((dir_sons)i);
//number_variable+(*Index_value_) besser machen !!!!!!!!!
      if(cetyp==bo_cell) {
        sum = sum +
          Evaluate_lsm(nearb_ablage->Give_u_Recell((dir_sons)i),
                       actual_number, dir_v,
                       gr->Give_bo_stencil(it_n->Ind().next((dir_sons)i,l+1),
                                           num_stencils+(*Index_value_)));
      }
    }    
    return sum;
  };
  double Give_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
                      const BoCeData* bocedata) const {
    return a_.Give_cellpoi(it_cf,gr,bocedata);
  };
  double Give_Bo2p(const P_Bo2p* it_b, const Grid* gr,int l)  const {
    return a_.Give_Bo2p(it_b,gr,l);
  };
  // Size of necessary stencil
  int Sice_stencil() const { return a_.Sice_stencil(); }

  // activity of points
  int Level() const { return (a_).Level(); }
  Dominace_label Dominant_lev() const { return not_dominant; }
  Dominace_label Dominant_poi() const { return not_dominant; }
  bool run_interior() const { return (a_).run_interior(); }
  bool run_nearb()    const { return (a_).run_nearb(); }
  int run_boundary() const { return (a_).run_boundary(); }

  // number of operations
  int ops_interior() const { return a_.ops_interior(); }

  // For operator ==
  void Active_Sim_Level( int lev) const {  };
  void Active_Sim_interior(bool run) const {  };
  void Active_Sim_nearb(bool run) const {  };
  void Active_Sim_boundary( int run) const {  };
  void Active_Sim_update(Evaluation_Parallelization_object* evpar,
                         int level, type_of_update typ) const {  };

  // Put pointer to grid and activity of boundary points of test space  
  void Put_grid_rbo(Grid* gr, int r_bo) const {
    a_.Put_grid_rbo(gr,r_bo);
  }; 
  // For abstract differential operators:
  // --------------------------------------
  // Is the expression an abstract differential operator?
  Differential_op_typ Abstract_differential_operator() const
    { return not_abstract; }
  //  bool Give_array_variable_inserted() const {
  //    return true;
  //  }

  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {
    if(array_variable_inserted == false)
      evpar->Variable_contained_in_expression(number_variable);
    else 
      for(int i=0;i<dim;++i)
        evpar->Variable_contained_in_expression(number_variable+i);
  };
  bool GS_type(int var_number_left, int dim_left) const {
    if(array_variable_inserted == false) {
      if(number_variable>=var_number_left &&
         number_variable<=(var_number_left+dim_left)) return true;
      else                                            return false;
    }
    else {
     if((number_variable+dim)>=var_number_left &&
        number_variable<=(var_number_left+dim_left)) return true;
     else                                            return false;
    }
  };

  // at the and of an iteration of an expression
  void clean() const { a_.clean(); };
};

class Index_Res_stencil_boundary_ARR {
 private:
  // Gitter
  Grid *grid;

  //number of storage for boundary stencils
  int num_stencils;
  int* Index_value_;
 public:
  Index_Res_stencil_boundary_ARR(int num_s, Grid* gr, int *Ivalue_) 
    : grid(gr), num_stencils(num_s),  Index_value_(Ivalue_) {
  };

  template<class A>
  DExpr<DResDiff_Bo_ARR<DExpr<A> > > 
  operator() (const DExpr<A>& a_) {
    Differential_op_typ a_abstract;
    a_abstract = a_.Abstract_differential_operator();
    if(a_abstract==not_abstract) {
      cout << " This is not an abstract differential operator! Error!" << endl;
    }
    if(a_abstract==abstract_diag)
      cout<<"\n Please use operator [] instead of operator ()! Error!" << endl;

    typedef DResDiff_Bo_ARR<DExpr<A> > ExprT;
    return DExpr<ExprT>(ExprT(a_,num_stencils,grid,Index_value_));
  }
  /*
  template<class A>
  DExpr<DResDiagDiff_Bo<DExpr<A> > > 
  operator[] (const DExpr<A>& a_) {
    Differential_op_typ a_abstract;
    a_abstract = a_.Abstract_differential_operator();
    if(a_abstract==not_abstract) {
      cout << " This is not an abstract differential operator! Error!" << endl;
    }
    if(a_abstract==abstract_with_var)
      cout<<"\n Please use operator () instead of operator []! Error!" << endl;
    typedef DResDiagDiff_Bo<DExpr<A> > ExprT;
    return DExpr<ExprT>(ExprT(a_,num_stencils,grid,this));
  }
  */
};


class Array_Res_stencil_boundary {
 public:
  Array_Res_stencil_boundary(int dimension, Grid* gr) {
    dim  = dimension;
    start = new Res_stencil_boundary[dim];
    for(int i=0;i<dim;++i) {
      start[i].Initialize(gr);
    }
  }
  Res_stencil_boundary& operator[](int i) const {
    return start[i];
  }
  Index_Res_stencil_boundary_ARR
  operator[](Local_int& I) const {
    for(int i=0;i<dim;++i) {
      if(start[i].Give_used()==false) {
        cout << " Error: Use Array_Res_stencil_boundary first time " 
             << " with operator operator[](int i) !! " << endl;
      }
    }
    return Index_Res_stencil_boundary_ARR(start[0].num_stencils,
                                          start[0].grid,
                                          I.Give_pointer());
  }
  int Give_dim() {
    return dim;
  }
 private:
  Res_stencil_boundary *start;
  int dim;
};






//    4. sums and for iteration



template<class A>
class Expr_for_iterator {
 private:
  A a_;
  int* value_;
  int begin;
  int end;
  
  double dummy;
 public:
  Expr_for_iterator(const A& x, int b, int e, int* v)
    : a_(x), value_(v), begin(b), end(e)
    {}
  // for evaluate
  void Run_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const {
    for((*value_)=begin;(*value_)<=end;(*value_)++)
      a_.Run_interior(it_i,gr,h_mesh,l, double_stencils_out);
  };
  void Run_interior_coarse(const P_interior* it_i,const Grid* gr,
                           double h_mesh,
                           int l, double_stencils_in) const {
    for((*value_)=begin;(*value_)<=end;(*value_)++)
      a_.Run_interior_coarse(it_i,gr,h_mesh,l, double_stencils_out);
  };
  void Run_nearb(const P_nearb* it_n,const Grid* gr, double h_mesh, int l,
                    const Nearb_Ablage* nearb_ablage)
    const {
    for((*value_)=begin;(*value_)<=end;(*value_)++)
      a_.Run_nearb(it_n,gr,h_mesh,l,nearb_ablage);
  };
  void Run_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    for((*value_)=begin;(*value_)<=end;(*value_)++)
      a_.Run_Bo2p(it_b,gr,l);
  };
  void Run_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                      const BoCeData* bocedata) const {
    for((*value_)=begin;(*value_)<=end;(*value_)++)
      a_.Run_cellpoi(it_cf,gr,bocedata);
  };

  // Size of necessary stencil
  int Sice_stencil() const { return a_.Sice_stencil(); }

  // activity of points
  int Level() const { return a_.Level(); }
  Dominace_label Dominant_lev() const { return a_.Dominant_lev(); }
  Dominace_label Dominant_poi() const { return a_.Dominant_poi(); }
  bool run_interior() const { return a_.run_interior(); }
  bool run_nearb()    const { return a_.run_nearb(); }
  int run_boundary() const { return a_.run_boundary(); }

  // number of operations
  int ops_interior() const { return a_.ops_interior()*(end-begin+1); }   

  // For operator ==
  void Active_Sim_Level( int lev) const { a_.Active_Sim_Level(lev); }
  void Active_Sim_interior( bool run) const { a_.Active_Sim_interior(run); }
  void Active_Sim_nearb( bool run) const { a_.Active_Sim_nearb(run); }
  void Active_Sim_boundary(int run) const { a_.Active_Sim_boundary(run); }
  void Active_Sim_update(Evaluation_Parallelization_object* evpar,
                         int level, type_of_update typ) const 
    {  a_.Active_Sim_update(evpar,level,typ); };

  // Put pointer to grid and activity of boundary points of test space  
  void Put_grid_rbo(Grid* gr, int r_bo) const {
    a_.Put_grid_rbo(gr,r_bo); }
  // Is the expression an abstract differential operator?
   Differential_op_typ Abstract_differential_operator() const
    { return not_abstract; }
  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {
    a_.Add_variables_for_parallel(evpar);
  };
  bool GS_type(int var_number_left, int dim_left) const {
    return a_.GS_type(var_number_left, dim_left);
  };

  // at the and of an iteration of an expression
  void clean() const { a_.clean(); };
};

template<class A>
class Expr_for_sum {
 private:
  A a_;
  int* value_;
  int begin;
  int end;
  
  double* sum;
 public:
  Expr_for_sum(const A& x, int b, int e, int* v) 
    : a_(x), value_(v), begin(b), end(e) {
    if(e<b) cout << "\n Empty Sum! Error! " << endl;
    sum = new double;
  };

  // for evaluate
  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const {
    (*value_)=begin;
    (*sum)=a_.Give_interior(it_i,gr,h_mesh,l, double_stencils_out);
    for((*value_)=begin+1;(*value_)<=end;(*value_)++)
      (*sum) = (*sum)+a_.Give_interior(it_i,gr,h_mesh,l, double_stencils_out);
    return (*sum);
  };
  double Give_interior_coarse(const P_interior* it_i,const Grid* gr, 
                              double h_mesh, int l, double_stencils_in) const {
    (*value_)=begin;
    (*sum) = a_.Give_interior_coarse(it_i,gr,h_mesh,l, double_stencils_out);
    for((*value_)=begin+1;(*value_)<=end;(*value_)++)
      (*sum)=(*sum)+
        a_.Give_interior_coarse(it_i,gr,h_mesh,l, double_stencils_out);
    return (*sum);
  };
  double Give_nearb(const P_nearb* it_n,const Grid* gr, double h_mesh, int l,
                    const Nearb_Ablage* nearb_ablage)
    const {
    (*value_)=begin;
    (*sum) = a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage);
    for((*value_)=begin+1;(*value_)<=end;(*value_)++)
      (*sum) = (*sum) + a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage);
    return (*sum);
  };
  double Give_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    (*value_)=begin;
    (*sum) = a_.Give_Bo2p(it_b,gr,l);
    for((*value_)=begin+1;(*value_)<=end;(*value_)++)
      (*sum) = (*sum) + a_.Give_Bo2p(it_b,gr,l);
    return (*sum);
  };
  double Give_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                      const BoCeData* bocedata) const {
    (*value_)=begin;
    (*sum) = a_.Give_cellpoi(it_cf,gr,bocedata);
    for((*value_)=begin+1;(*value_)<=end;(*value_)++)
      (*sum) = (*sum) + a_.Give_cellpoi(it_cf,gr,bocedata);
    return (*sum);
  };

  // Size of necessary stencil
  int Sice_stencil() const { return a_.Sice_stencil(); }

  // activity of points
  int Level() const { return a_.Level(); }
  Dominace_label Dominant_lev() const { return a_.Dominant_lev(); }
  Dominace_label Dominant_poi() const { return a_.Dominant_poi(); }
  bool run_interior() const { return a_.run_interior(); }
  bool run_nearb()    const { return a_.run_nearb(); }
  int run_boundary() const { return a_.run_boundary(); }

  // number of operations
  int ops_interior() const { return a_.ops_interior()*(end-begin+1)
                               +(end-begin); }   

  // For operator ==
  void Active_Sim_Level( int lev) const { a_.Active_Sim_Level(lev); }
  void Active_Sim_interior( bool run) const { a_.Active_Sim_interior(run); }
  void Active_Sim_nearb( bool run)    const { a_.Active_Sim_nearb(run); }
  void Active_Sim_boundary(int run) const { a_.Active_Sim_boundary(run); }
  void Active_Sim_update(Evaluation_Parallelization_object* evpar,
                         int level, type_of_update typ) const 
    {  a_.Active_Sim_update(evpar,level,typ); };

  // Put pointer to grid and activity of boundary points of test space
  void Put_grid_rbo(Grid* gr, int r_bo)  const { a_.Put_grid_rbo(gr,r_bo); }
  // For abstract differential operators:
  // --------------------------------------
  // Is the expression an abstract differential operator?
  Differential_op_typ Abstract_differential_operator() const  {
    return not_abstract; }
  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {
    a_.Add_variables_for_parallel(evpar);
  };
  bool GS_type(int var_number_left, int dim_left) const {
    return a_.GS_type(var_number_left, dim_left);
  };

  // at the and of an iteration of an expression
  void clean() const { delete(sum); a_.clean(); };
};



class Schleife {
 public:
  int* value_;
  int begin;
  int end;
  Schleife(int* v, int b, int e) : value_(v), begin(b), end(e)
    {};
  template<class A>
    DWrapSim<Expr_for_iterator<DWrapSim<A> > >
    operator()(const DWrapSim<A> & a)
    {
      typedef Expr_for_iterator<DWrapSim<A> >  ExprT;
      return DWrapSim<ExprT>(ExprT(a,begin,end,value_));
    }
};




class SchleifeSum {
 public:
  int* value_;
  int begin;
  int end;
  SchleifeSum(int* v, int b, int e) : value_(v), begin(b), end(e)
    {};
  template<class A>
    DExpr<Expr_for_sum<DExpr<A> > >
    operator()(const DExpr<A> & a)
    {
      typedef Expr_for_sum<DExpr<A> >  ExprT;
      return DExpr<ExprT>(ExprT(a,begin,end,value_));
    }
};





//  for iterator function


Schleife For(Local_int& I, int begin, int end);


SchleifeSum Sum(Local_int& I, int begin, int end);

// 5.:  Class for procedure

// 5.1.a. procedure with 3  parameters, and per reference
template<class A, class B, class C>
class DExprProcedure3 {
 private:
  void (*Procedure)(double* x, double* y, double* z);
  int opsf; // number of operations for Formula
  A a_;
  B b_;
  C c_;
 public:
  DExprProcedure3(void (*procedure)(double* x, double* y, double* z),
             const A& a, const B& b, const C& c) : a_(a), b_(b), c_(c)
    {
      Procedure = procedure;
      opsf = 0;
    }
  // for evaluate
  void Run_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const {
    Procedure(a_.Give_regular_pointer(it_i,gr,l),
              b_.Give_regular_pointer(it_i,gr,l),
              c_.Give_regular_pointer(it_i,gr,l));
  };
  void Run_interior_coarse(const P_interior* it_i,const Grid* gr,
                           double h_mesh,
                           int l, double_stencils_in) const {
    Procedure(a_.Give_regular_pointer(it_i,gr,l),
              b_.Give_regular_pointer(it_i,gr,l),
              c_.Give_regular_pointer(it_i,gr,l));
  };
  void Run_nearb(const P_nearb* it_n,const Grid* gr, double h_mesh, int l,
                    const Nearb_Ablage* nearb_ablage)
    const {
    Procedure(a_.Give_regular_pointer(it_n,gr,l),
              b_.Give_regular_pointer(it_n,gr,l),
              c_.Give_regular_pointer(it_n,gr,l));
  };
  void Run_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    Procedure(a_.Give_Bo2p_pointer(it_b,gr,l),
              b_.Give_Bo2p_pointer(it_b,gr,l),
              c_.Give_Bo2p_pointer(it_b,gr,l));
  };
  void Run_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
                   const BoCeData* bocedata) const {
    Procedure(a_.Give_cellpoi_pointer(it_cf,gr),
              b_.Give_cellpoi_pointer(it_cf,gr),
              c_.Give_cellpoi_pointer(it_cf,gr));
  };

  // Size of necessary stencil
  int Sice_stencil() const { return 1; }

  // activity of points
  int Level() const { 
    if(c_.Dominant_lev() ==  yes_dominant || 
       (b_.Dominant_lev() == anti_dominant &&
        a_.Dominant_lev() == anti_dominant))
      return c_.Level();
    if(b_.Dominant_lev() ==  yes_dominant || 
       (c_.Dominant_lev() == anti_dominant &&
        a_.Dominant_lev() == anti_dominant))
      return b_.Level();
    if(a_.Dominant_lev() ==  yes_dominant || 
       (b_.Dominant_lev() == anti_dominant &&
        c_.Dominant_lev() == anti_dominant))
      return a_.Level();
   if(c_.Dominant_lev() == anti_dominant)
     return MAX(a_.Level(),b_.Level()); 
   if(b_.Dominant_lev() == anti_dominant)
     return MAX(a_.Level(),c_.Level()); 
   if(a_.Dominant_lev() == anti_dominant)
     return MAX(c_.Level(),b_.Level());
   return MAX(a_.Level(),b_.Level(),c_.Level()); 
  }
  Dominace_label Dominant_lev() const {
    return (Dominace_label)(MAX(a_.Dominant_lev(),
                                b_.Dominant_lev(),
                                c_.Dominant_lev()));
  }
  Dominace_label Dominant_poi() const {
    return (Dominace_label)(MAX(a_.Dominant_poi(),
                                b_.Dominant_poi(),
                                c_.Dominant_poi()));
  }
  bool run_interior() const { 
    if(c_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return c_.run_interior();
    if(b_.Dominant_poi() ==  yes_dominant || 
       (c_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return b_.run_interior();
    if(a_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        c_.Dominant_poi() == anti_dominant))
      return a_.run_interior();
    if(c_.Dominant_poi() == anti_dominant)
      return b_.run_interior();
    return c_.run_interior();
  }
  bool run_nearb() const { 
    if(c_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return c_.run_nearb();
    if(b_.Dominant_poi() ==  yes_dominant || 
       (c_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return b_.run_nearb();
    if(a_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        c_.Dominant_poi() == anti_dominant))
      return a_.run_nearb();
    if(c_.Dominant_poi() == anti_dominant)
      return b_.run_nearb();
    return c_.run_nearb();
  }
  bool run_boundary() const { 
    if(c_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return c_.run_boundary();
    if(b_.Dominant_poi() ==  yes_dominant || 
       (c_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return b_.run_boundary();
    if(a_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        c_.Dominant_poi() == anti_dominant))
      return a_.run_boundary();
    if(c_.Dominant_poi() == anti_dominant)
      return b_.run_boundary();
    return c_.run_boundary();
  }

  // number of operations
  int ops_interior() const { 
    return opsf;
  }
  // For operator ==
  void Active_Sim_Level( int lev) const { 
    a_.Active_Sim_Level(lev);
    b_.Active_Sim_Level(lev);
    c_.Active_Sim_Level(lev); 
  };
  void Active_Sim_interior(bool run) const { 
    a_.Active_Sim_interior(run);
    b_.Active_Sim_interior(run);
    c_.Active_Sim_interior(run); 
  };
  void Active_Sim_nearb(bool run) const { 
    a_.Active_Sim_nearb(run); 
    b_.Active_Sim_nearb(run); 
    c_.Active_Sim_nearb(run); 
  };
  void Active_Sim_boundary( int run) const {
    a_.Active_Sim_boundary(run); 
    b_.Active_Sim_boundary(run); 
    c_.Active_Sim_boundary(run); 
  };
  void Active_Sim_update(Evaluation_Parallelization_object* evpar,
                         int level, type_of_update typ) const {
    a_.Active_Sim_update(evpar,level,typ);
    b_.Active_Sim_update(evpar,level,typ); 
    b_.Active_Sim_update(evpar,level,typ); 
  };

  // Put pointer to grid and activity of boundary points of test space  
  void Put_grid_rbo(Grid* gr, int r_bo) const {
    a_.Put_grid_rbo(gr,r_bo);
    b_.Put_grid_rbo(gr,r_bo); 
    c_.Put_grid_rbo(gr,r_bo); 
  }
 // For abstract differential operators:
  // --------------------------------------
  // Is the expression an abstract differential operator?
  Differential_op_typ Abstract_differential_operator()  const {
    return not_abstract;
  };
 // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar)
    const {
    a_.Add_variables_for_parallel(evpar);
    b_.Add_variables_for_parallel(evpar);
    c_.Add_variables_for_parallel(evpar);
  };
  bool GS_type(int var_number_left, int dim_left) const {
   if(c_.GS_type(var_number_left, dim_left)==true) return true;
   if(b_.GS_type(var_number_left, dim_left)==true) return true;
    return a_.GS_type(var_number_left, dim_left);
  };

  // at the and of an iteration of an expression
  void clean() const  {
    a_.clean(); 
    b_.clean(); 
    c_.clean(); 
  };
};

/*
// 5.1.b. procedure with 3  parameters, and not per reference
template<class A, class B, class C>
class DExprProcedure3_nr {
 private:
  void (*Procedure)(double x, double y, double z);
  int opsf; // number of operations for Formula
  A a_;
  B b_;
  C c_;
 public:
  DExprProcedure3(void (*procedure)(double x, double y, double z),
             const A& a, const B& b, const C& c) : a_(a), b_(b), c_(c)
    {
      Procedure = procedure;
      opsf = 0;
    }
  // for evaluate
  void Run_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                      int l, double_stencils_in) const { 
    Procedure(a_.Give_interior(it_i,gr,h_mesh,l, double_stencils_out),
                   b_.Give_interior(it_i,gr,h_mesh,l, double_stencils_out),
                   c_.Give_interior(it_i,gr,h_mesh,l, double_stencils_out)); 
  }
  void Run_interior_coarse(const P_interior* it_i,const Grid* gr,
                           double h_mesh,
                           int l, double_stencils_in) const { 
    Procedure(a_.Give_interior_coarse(it_i,gr,h_mesh,l,
                                      double_stencils_out),
              b_.Give_interior_coarse(it_i,gr,h_mesh,l,
                                      double_stencils_out),
              c_.Give_interior_coarse(it_i,gr,h_mesh,l,
                                      double_stencils_out)); 
  }
  void Run_nearb(const P_nearb* it_n,const Grid* gr, double h_mesh, int l,
                 const Nearb_Ablage* nearb_ablage) const { 
    Procedure(a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage),
              b_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage),
              c_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage)); 
  }
  void Run_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    Procedure(a_.Give_Bo2p(it_b,gr,l),
              b_.Give_Bo2p(it_b,gr,l),
              c_.Give_Bo2p(it_b,gr,l)); 
  }
  void Run_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                   const BoCeData* bocedata) const {
    Procedure(a_.Give_cellpoi(it_cf,gr,bocedata),
              b_.Give_cellpoi(it_cf,gr,bocedata),
              c_.Give_cellpoi(it_cf,gr,bocedata)); 
  }
  
  // Size of necessary stencil
  int Sice_stencil() const { return 1; }

  // activity of points
  int Level() const { 
    if(c_.Dominant_lev() ==  yes_dominant || 
       (b_.Dominant_lev() == anti_dominant &&
        a_.Dominant_lev() == anti_dominant))
      return c_.Level();
    if(b_.Dominant_lev() ==  yes_dominant || 
       (c_.Dominant_lev() == anti_dominant &&
        a_.Dominant_lev() == anti_dominant))
      return b_.Level();
    if(a_.Dominant_lev() ==  yes_dominant || 
       (b_.Dominant_lev() == anti_dominant &&
        c_.Dominant_lev() == anti_dominant))
      return a_.Level();
   if(c_.Dominant_lev() == anti_dominant)
     return MAX(a_.Level(),b_.Level()); 
   if(b_.Dominant_lev() == anti_dominant)
     return MAX(a_.Level(),c_.Level()); 
   if(a_.Dominant_lev() == anti_dominant)
     return MAX(c_.Level(),b_.Level());
   return MAX(a_.Level(),b_.Level(),c_.Level()); 
  }
  Dominace_label Dominant_lev() const {
    return (Dominace_label)(MAX(a_.Dominant_lev(),
                                b_.Dominant_lev(),
                                c_.Dominant_lev()));
  }
  Dominace_label Dominant_poi() const {
    return (Dominace_label)(MAX(a_.Dominant_poi(),
                                b_.Dominant_poi(),
                                c_.Dominant_poi()));
  }
  bool run_interior() const { 
    if(c_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return c_.run_interior();
    if(b_.Dominant_poi() ==  yes_dominant || 
       (c_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return b_.run_interior();
    if(a_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        c_.Dominant_poi() == anti_dominant))
      return a_.run_interior();
    if(c_.Dominant_poi() == anti_dominant)
      return b_.run_interior();
    return c_.run_interior();
  }
  bool run_nearb() const { 
    if(c_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return c_.run_nearb();
    if(b_.Dominant_poi() ==  yes_dominant || 
       (c_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return b_.run_nearb();
    if(a_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        c_.Dominant_poi() == anti_dominant))
      return a_.run_nearb();
    if(c_.Dominant_poi() == anti_dominant)
      return b_.run_nearb();
    return c_.run_nearb();
  }
  bool run_boundary() const { 
    if(c_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return c_.run_boundary();
    if(b_.Dominant_poi() ==  yes_dominant || 
       (c_.Dominant_poi() == anti_dominant &&
        a_.Dominant_poi() == anti_dominant))
      return b_.run_boundary();
    if(a_.Dominant_poi() ==  yes_dominant || 
       (b_.Dominant_poi() == anti_dominant &&
        c_.Dominant_poi() == anti_dominant))
      return a_.run_boundary();
    if(c_.Dominant_poi() == anti_dominant)
      return b_.run_boundary();
    return c_.run_boundary();
  }

  // number of operations
  int ops_interior() const { 
    return opsf;
  }
  // For operator ==
  void Active_Sim_Level( int lev) const { 
    a_.Active_Sim_Level(lev);
    b_.Active_Sim_Level(lev);
    c_.Active_Sim_Level(lev); 
  };
  void Active_Sim_interior(bool run) const { 
    a_.Active_Sim_interior(run);
    b_.Active_Sim_interior(run);
    c_.Active_Sim_interior(run); 
  };
  void Active_Sim_nearb(bool run) const { 
    a_.Active_Sim_nearb(run); 
    b_.Active_Sim_nearb(run); 
    c_.Active_Sim_nearb(run); 
  };
  void Active_Sim_boundary( int run) const {
    a_.Active_Sim_boundary(run); 
    b_.Active_Sim_boundary(run); 
    c_.Active_Sim_boundary(run); 
  };
  void Active_Sim_update(Evaluation_Parallelization_object* evpar,
                         int level, type_of_update typ) const {
    a_.Active_Sim_update(evpar,level,typ);
    b_.Active_Sim_update(evpar,level,typ); 
    b_.Active_Sim_update(evpar,level,typ); 
  };

  // Put pointer to grid and activity of boundary points of test space  
  void Put_grid_rbo(Grid* gr, int r_bo) const {
    a_.Put_grid_rbo(gr,r_bo);
    b_.Put_grid_rbo(gr,r_bo); 
    c_.Put_grid_rbo(gr,r_bo); 
  }
 // For abstract differential operators:
  // --------------------------------------
  // Is the expression an abstract differential operator?
  Differential_op_typ Abstract_differential_operator()  const {
    return not_abstract;
  };
 // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar)
    const {
    a_.Add_variables_for_parallel(evpar);
    b_.Add_variables_for_parallel(evpar);
    c_.Add_variables_for_parallel(evpar);
  };
  bool GS_type(int var_number_left, int dim_left) const {
   if(c_.GS_type(var_number_left, dim_left)==true) return true;
   if(b_.GS_type(var_number_left, dim_left)==true) return true;
    return a_.GS_type(var_number_left, dim_left);
  };

  // at the and of an iteration of an expression
  void clean() const  {
    a_.clean(); 
    b_.clean(); 
    c_.clean(); 
  };
};
*/


// 5.3. procedure definition
class Procedure {
  private:
  void (*Procedure3)(double* x, double* y,double* z);
  void (*Procedure4)(double* x, double* y,double* z, double* u);

  int type;  // 3 or 4 

  int length_x, length_y, length_z; // only for array
 public:

  Procedure(void (*procedure)(double* x, double* y,double* z),
            int len_x, int len_y, int len_z)
    {
      Procedure3 = procedure;
      type = 3;
      length_x = len_x;
      length_y = len_y;
      length_z = len_z;
    }
  Procedure(void (*procedure)(double* x, double* y,double* z))
    {
      Procedure3 = procedure;
      type = 3;
      length_x = 0;
      length_y = 0;
      length_z = 0;
    }
  // procedure with 3 entries
  // a) for array
  DWrapSim<DExprProcedure3<DExprVAR, DExprVAR, DExprVAR> > operator() 
     (Array_Variable& a, Array_Variable& b, Array_Variable& c) {
    if(a.Give_dim()<length_x &&
       b.Give_dim()<length_y &&
       c.Give_dim()<length_z) {
      cout << "\n Wrong length of vectors! " << endl;
    }
    if(type != 3) cout << "\n Wrong usage of procedure 3! " << endl;
    typedef DExprProcedure3<DExprVAR, DExprVAR, DExprVAR>  ExprT;
    return DWrapSim<ExprT>(ExprT(Procedure3,
      DExprVAR(a[0]),DExprVAR(b[0]),DExprVAR(c[0])));
  }

  DWrapSim<DExprProcedure3<Local_var, DExprVAR, DExprVAR> > operator() 
     (Array_Local_var& a, Array_Variable& b, Array_Variable& c) {
    if(a.Give_dim()<length_x &&
       b.Give_dim()<length_y &&
       c.Give_dim()<length_z) {
      cout << "\n Wrong length of vectors! " << endl;
    }
    if(type != 3) cout << "\n Wrong usage of procedure 3! " << endl;
    typedef DExprProcedure3<Local_var, DExprVAR, DExprVAR>  ExprT;
    return DWrapSim<ExprT>(ExprT(Procedure3,
      a[0],DExprVAR(b[0]),DExprVAR(c[0])));
  }

  DWrapSim<DExprProcedure3<DExprVAR, Local_var, DExprVAR> > operator() 
     (Array_Variable& a, Array_Local_var& b, Array_Variable& c) {
    if(a.Give_dim()<length_x &&
       b.Give_dim()<length_y &&
       c.Give_dim()<length_z) {
      cout << "\n Wrong length of vectors! " << endl;
    }
    if(type != 3) cout << "\n Wrong usage of procedure 3! " << endl;
    typedef DExprProcedure3<DExprVAR, Local_var, DExprVAR>  ExprT;
    return DWrapSim<ExprT>(ExprT(Procedure3,
      DExprVAR(a[0]),b[0],DExprVAR(c[0])));
  }

  DWrapSim<DExprProcedure3<DExprVAR, DExprVAR, Local_var> > operator() 
     (Array_Variable& a, Array_Variable& b, Array_Local_var& c) {
    if(a.Give_dim()<length_x &&
       b.Give_dim()<length_y &&
       c.Give_dim()<length_z) {
      cout << "\n Wrong length of vectors! " << endl;
    }
    if(type != 3) cout << "\n Wrong usage of procedure 3! " << endl;
    typedef DExprProcedure3<DExprVAR, DExprVAR, Local_var>  ExprT;
    return DWrapSim<ExprT>(ExprT(Procedure3,
      DExprVAR(a[0]),DExprVAR(b[0]),c[0]));
  }
  /*
  // b) for scalars
  template<class A, class B, class C>
    DWrapSim<DExprProcedure3_nr<A,B,C> > operator() (A a, B b, C c) {
    if(type != 3) cout << "\n Wrong usage of procedure 3! " << endl;
    typedef DExprProcedure3_nr<A, B, C>  ExprT;
    return DWrapSim<ExprT>(ExprT(Procedure3,a,b,c));
  }
  */
};


#endif
        

---