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

#ifndef ELEME_H_
#define ELEME_H_
      
// ------------------------------------------------------------
// elementary.h
//
// ------------------------------------------------------------

           // 1.:  Coordinate functions
           // 2.:  Class for Functions
           // 3.:  Restriction operator with respect to points
           // 4.:  Binary operators

// 1. Coordinate functions

// X Koordinate
class CoordinateX {
 public:
  CoordinateX() {};
  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const {
    return it_i->coordinateX(gr);
  };
  double Give_interior_coarse(const P_interior* it_i,const Grid* gr,
                              double h_mesh,
                              int l, double_stencils_in) const {
    return it_i->coordinateX(gr);
  };
  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->coordinateX(gr);
  };
  double Give_Bo2p(const P_Bo2p* it_b,const  Grid* gr, int l) const  {
    return it_b->coordinateX(gr);
  };
  double Give_cellpoi(const P_cellpoi* it_cf,const  Grid* gr,
                      const BoCeData* bocedata) const  {
    return it_cf->coordinateX(gr);
  };

  // for evaluate on cell
  double Give_interior_cell(const P_interior_cell* it_icell ,
                            const Grid* gr, int l) const {
    return it_icell->coordinateX(gr);
  };
  double Give_bo_tet(const P_boundary_tet* it_bcell, 
                      const Grid* gr) const {
    return it_bcell->coordinateX(gr);
  };

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

  // 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 true; }

  // Zeitmessung
  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 not_abstract; }
  // 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  { };
};


DExpr<CoordinateX> X();

// Y Koordinate
class CoordinateY {
 public:
  CoordinateY() {};
  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const  {
    return it_i->coordinateY(gr);
  };
  double Give_interior_coarse(const P_interior* it_i,const Grid* gr,
                              double h_mesh,
                              int l, double_stencils_in) const  {
    return it_i->coordinateY(gr);
  };
  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->coordinateY(gr);
  };
  double Give_Bo2p(const P_Bo2p* it_b,const  Grid* gr, int l)  const {
    return it_b->coordinateY(gr);
  };
  double Give_cellpoi(const P_cellpoi* it_cf,const  Grid* gr,
                      const BoCeData* bocedata)  const {
    return it_cf->coordinateY(gr);
  };

  // for evaluate on cell
  double Give_interior_cell(const P_interior_cell* it_icell ,
                            const Grid* gr, int l) const {
    return it_icell->coordinateY(gr);
  };
  double Give_bo_tet(const P_boundary_tet* it_bcell, 
                      const Grid* gr) const {
    return it_bcell->coordinateY(gr);
  };

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

  // 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 true; }

  // Zeitmessung
  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?
  bool Abstract_differential_operator()  const { return false; }
  // 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  { };
};


DExpr<CoordinateY> Y();

// Z Koordinate
class CoordinateZ {
 public:
  CoordinateZ() {};
  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in)  const {
    return it_i->coordinateZ(gr);
  };
  double Give_interior_coarse(const P_interior* it_i,const Grid* gr,
                              double h_mesh,
                              int l, double_stencils_in)  const {
    return it_i->coordinateZ(gr);
  };
  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->coordinateZ(gr);
  };
  double Give_Bo2p(const P_Bo2p* it_b,const  Grid* gr, int l) const  {
    return it_b->coordinateZ(gr);
  };
  double Give_cellpoi(const P_cellpoi* it_cf,const  Grid* gr,
                      const BoCeData* bocedata) const  {
    return it_cf->coordinateZ(gr);
  };

  // for evaluate on cell
  double Give_interior_cell(const P_interior_cell* it_icell ,
                            const Grid* gr, int l) const {
    return it_icell->coordinateZ(gr);
  };
  double Give_bo_tet(const P_boundary_tet* it_bcell, 
                      const Grid* gr) const {
    return it_bcell->coordinateZ(gr);
  };

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

  // 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 true; }

  // Zeitmessung
  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 not_abstract; }
  // 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  { };
};


DExpr<CoordinateZ> Z();



// 2.:  Class for Functions


template<class A>
class DExprFunc1 {
 private:
  double (*Formula)(double x);
  int opsf; // number of operations for Formula
  A a_;
 public:
  DExprFunc1(double (*formula)(double x), const A& a) : a_(a)
    {
      Formula = formula;
      opsf = 0;
    }

  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const { 
    return Formula(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 Formula(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 { return Formula(a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage)); }
  double Give_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    return Formula(a_.Give_Bo2p(it_b,gr,l)); }
  double Give_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                      const BoCeData* bocedata) const {
    return Formula(a_.Give_cellpoi(it_cf,gr,bocedata)); }

  // for evaluate on cell
  double Give_interior_cell(const P_interior_cell* it_icell ,
                            const Grid* gr, int l) const {
    return Formula(a_.Give_interior_cell(it_icell,gr,l));
  };
  double Give_bo_tet(const P_boundary_tet* it_bcell, 
                     const Grid* gr) const {
    return Formula(a_.Give_bo_tet(it_bcell,gr));
  };

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

  // activity of points
  int Level() const { return a_.Level(); }
  int Dominant() const { return a_.Dominant(); }
  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 opsf + a_.ops_interior(); }

  // 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; }
  // at the and of an iteration of an expression
 // 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) const {
    return a_.GS_type(var_number_left,dim);
  }; 
  void clean() const  { 
    a_.clean();  
  };
};


template<class A, class B, class C>
class DExprFunc3 {
 private:
  double (*Formula)(double x, double y, double z);
  int opsf; // number of operations for Formula
  A a_;
  B b_;
  C c_;
 public:
  DExprFunc3(double (*formula)(double x, double y, double z),
             const A& a, const B& b, const C& c) : a_(a), b_(b), c_(c)
    {
      Formula = formula;
      opsf = 0;
    }
  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const { 
    return Formula(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)); }
  double Give_interior_coarse(const P_interior* it_i,const Grid* gr,
                              double h_mesh,
                              int l, double_stencils_in) const { 
    return Formula(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)); }
  double Give_nearb(const P_nearb* it_n,const Grid* gr, double h_mesh, int l,
                    const Nearb_Ablage* nearb_ablage)
    const { return Formula(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)); }
  double Give_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    return Formula(a_.Give_Bo2p(it_b,gr,l),
                   b_.Give_Bo2p(it_b,gr,l),
                   c_.Give_Bo2p(it_b,gr,l)); }
  double Give_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                      const BoCeData* bocedata) const {
    return Formula(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 MAX(MAX(a_.Sice_stencil(),b_.Sice_stencil()),
               c_.Sice_stencil()); }

  // 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 a_.ops_interior()+b_.ops_interior()+c_.ops_interior()+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); 
    c_.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) const {
    if(a_.GS_type(var_number_left,dim)==true) return true;
    if(b_.GS_type(var_number_left,dim)==true) return true;
    return c_.GS_type(var_number_left,dim);
  };

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

template<class A, class B, class C, class D>
class DExprFunc4 {
 private:
  double (*Formula)(double x, double y, double z, double w);
  A a_;
  B b_;
  C c_;
  D d_;
 public:
  DExprFunc4(double (*formula)(double x, double y, double z, double w),
             const A& a, const B& b, const C& c, const D& d) 
    : a_(a), b_(b), c_(c), d_(d)
    {
      Formula = formula;
    }
  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const { 
    return Formula(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),
                   d_.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 Formula(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),
                   d_.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 { return Formula(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),
                           d_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage)); }
  double Give_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    return Formula(a_.Give_Bo2p(it_b,gr,l),
                   b_.Give_Bo2p(it_b,gr,l),
                   c_.Give_Bo2p(it_b,gr,l),
                   d_.Give_Bo2p(it_b,gr,l)); }
  double Give_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                      const BoCeData* bocedata) const {
    return Formula(a_.Give_cellpoi(it_cf,gr,bocedata),
                   b_.Give_cellpoi(it_cf,gr,bocedata),
                   c_.Give_cellpoi(it_cf,gr,bocedata),
                   d_.Give_cellpoi(it_cf,gr,bocedata)); }



  // Size of necessary stencil
  int Sice_stencil() const { 
    return MAX(MAX(a_.Sice_stencil(),b_.Sice_stencil()),
               MAX(c_.Sice_stencil(),d_.Sice_stencil())); }

  // activity of points
  int Level() const { 
    if(a_.Dominant_lev() ==  yes_dominant) return a_.Level();
    if(b_.Dominant_lev() ==  yes_dominant) return b_.Level();
    if(c_.Dominant_lev() ==  yes_dominant) return c_.Level();
    if(d_.Dominant_lev() ==  yes_dominant) return d_.Level();

    if(a_.Dominant_lev() !=  anti_dominant) return a_.Level();
    if(b_.Dominant_lev() !=  anti_dominant) return b_.Level();
    if(c_.Dominant_lev() !=  anti_dominant) return c_.Level();
    if(d_.Dominant_lev() !=  anti_dominant) return d_.Level();
    return a_.Level();
  }
  Dominace_label Dominant_lev() const {
    return (Dominace_label)(MAX(MAX(a_.Dominant_lev(),
                                    b_.Dominant_lev()),
                                MAX(c_.Dominant_lev(),
                                    d_.Dominant_lev())));
  }
  Dominace_label Dominant_poi() const {
    return (Dominace_label)(MAX(MAX(a_.Dominant_poi(),
                                    b_.Dominant_poi()),
                                MAX(c_.Dominant_poi(),
                                    d_.Dominant_poi())));
  }
  bool run_interior() const { 
    if(a_.Dominant_poi() ==  yes_dominant)  return a_.run_interior();
    if(b_.Dominant_poi() ==  yes_dominant)  return b_.run_interior();
    if(c_.Dominant_poi() ==  yes_dominant)  return c_.run_interior();
    if(d_.Dominant_poi() ==  yes_dominant)  return d_.run_interior();


    if(a_.Dominant_poi() !=  anti_dominant)  return a_.run_interior();
    if(b_.Dominant_poi() !=  anti_dominant)  return b_.run_interior();
    if(c_.Dominant_poi() !=  anti_dominant)  return c_.run_interior();
    if(d_.Dominant_poi() !=  anti_dominant)  return d_.run_interior();

    return a_.run_interior();
  }
  bool run_nearb() const { 
    if(a_.Dominant_poi() ==  yes_dominant)  return a_.run_nearb();
    if(b_.Dominant_poi() ==  yes_dominant)  return b_.run_nearb();
    if(c_.Dominant_poi() ==  yes_dominant)  return c_.run_nearb();
    if(d_.Dominant_poi() ==  yes_dominant)  return d_.run_nearb();


    if(a_.Dominant_poi() !=  anti_dominant)  return a_.run_nearb();
    if(b_.Dominant_poi() !=  anti_dominant)  return b_.run_nearb();
    if(c_.Dominant_poi() !=  anti_dominant)  return c_.run_nearb();
    if(d_.Dominant_poi() !=  anti_dominant)  return d_.run_nearb();
    return a_.run_nearb();
  }
  bool run_boundary() const { 
    if(a_.Dominant_poi() ==  yes_dominant)  return a_.run_boundary();
    if(b_.Dominant_poi() ==  yes_dominant)  return b_.run_boundary();
    if(c_.Dominant_poi() ==  yes_dominant)  return c_.run_boundary();
    if(d_.Dominant_poi() ==  yes_dominant)  return d_.run_boundary();


    if(a_.Dominant_poi() !=  anti_dominant)  return a_.run_boundary();
    if(b_.Dominant_poi() !=  anti_dominant)  return b_.run_boundary();
    if(c_.Dominant_poi() !=  anti_dominant)  return c_.run_boundary();
    if(d_.Dominant_poi() !=  anti_dominant)  return d_.run_boundary();
    return a_.run_boundary();
  }

  // number of operations
  int ops_interior() const { 
    return 0; 
  }
  // For operator ==
  void Active_Sim_Level( int lev) const { 
    a_.Active_Sim_Level(lev);
    b_.Active_Sim_Level(lev);
    c_.Active_Sim_Level(lev); 
    d_.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); 
    d_.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); 
    d_.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);     
    d_.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); 
    c_.Active_Sim_update(evpar,level,typ); 
    d_.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); 
    d_.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);
    d_.Add_variables_for_parallel(evpar);
  };
  bool GS_type(int var_number_left, int dim) const {
    if(a_.GS_type(var_number_left,dim)==true) return true;
    if(b_.GS_type(var_number_left,dim)==true) return true;
    if(c_.GS_type(var_number_left,dim)==true) return true;
    return d_.GS_type(var_number_left,dim);
  };

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


class Function {
 private:
  double (*Formula)(double x);
  double (*Formula3)(double x, double y,double z);
  double (*Formula4)(double x, double y,double z, double w);

  int type;  // 1 or 3 or 4
 public:
  Function(double (*formula)(double x))
    {
      Formula = formula;
      type = 1;
    }
  Function(double (*formula)(double x, double y,double z))
    {
      Formula3 = formula;
      type = 3;
    }
  Function(double (*formula)(double x, double y,double z, double w))
    {
      Formula4 = formula;
      type = 4;
    }
  template<class A>
    DExpr<DExprFunc1<A> > operator() (A a) {
    if(type != 1) cout << "\n Wrong usage of function 1! " << endl;
    typedef DExprFunc1<A>  ExprT;
    return DExpr<ExprT>(ExprT(Formula,a));
  }
  DExpr<DExprFunc1<DExprLIT> > operator() (double x) {
    if(type != 1) cout << "\n Wrong usage of function 1! " << endl;
    typedef DExprFunc1<DExprLIT>  ExprT;
    return DExpr<ExprT>(ExprT(Formula,DExprLIT(x)));
  }
  DExpr<DExprFunc1<DExprVAR> > operator() (Variable& a) {
    if(type != 1) cout << "\n Wrong usage of function 1! " << endl;
    typedef DExprFunc1<DExprVAR>  ExprT;
    return DExpr<ExprT>(ExprT(Formula,DExprVAR(a)));
  }
  // function with 3 entries
  template<class A, class B, class C>
    DExpr<DExprFunc3<A,B,C> > operator() (A a, B b, C c) {
    if(type != 3) cout << "\n Wrong usage of function 3! " << endl;
    typedef DExprFunc3<A,B,C>  ExprT;
    return DExpr<ExprT>(ExprT(Formula,a,b,c));
  }
  DExpr<DExprFunc3<DExprVAR, DExprVAR, DExprVAR> > operator() (Variable& a,
                                                               Variable& b,
                                                               Variable& c) {
    if(type != 3) cout << "\n Wrong usage of function 3! " << endl;
    typedef DExprFunc3<DExprVAR, DExprVAR, DExprVAR>  ExprT;
    return DExpr<ExprT>(ExprT(Formula3,DExprVAR(a),DExprVAR(b),DExprVAR(c)));
  }
  // function with 4 entries
  template<class A, class B, class C, class D>
    DExpr<DExprFunc4<A,B,C,D> > operator() (A a, B b, 
                                            C c, D d) {
    if(type != 4) cout << "\n Wrong usage of function 4! " << endl;
    typedef DExprFunc4<A,B,C,D>  ExprT;
    return DExpr<ExprT>(ExprT(Formula4,a,b,c,d));
  }
};



 // 3.1.:  Restriction operator with respect to points -- fast

/* DExprResP for Representing a Restriction operator with respect to points */
// a) for an expression
template<class A, class R>
class DExprResP {
  A a_;
  R r_;
 public:
  DExprResP(const A& a, const R& r)
    : a_(a), r_(r)
    {}
  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 {
    return a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage);
  };
  double Give_Bo2p(const P_Bo2p* it_b, const Grid* gr, int l) const {
    return a_.Give_Bo2p(it_b,gr,l);
  };
  double Give_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
                      const BoCeData* bocedata)
    const {
    return a_.Give_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 yes_dominant; }
  bool run_interior() const { return r_.run_interior(); }
  bool run_nearb()    const { return r_.run_nearb(); }
  int run_boundary() const { return r_.run_boundary(); }

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

  // 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 a_.Abstract_differential_operator();
  }
  // 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) const {
    return a_.GS_type(var_number_left,dim); 
  };
  // at the and of an iteration of an expression
  void clean() const  { a_.clean(); };
};

// b) for a variable
template<class R>
class DExprResVarP {
 private:
  // active level
  int  level;        // level:

  // Nummer der Variablen
  int number_variable;

  // Restriction class
  R r_;
 public:
  DExprResVarP(const Variable& vec, const R& r) : r_(r) {
    level      = vec.Level();
    number_variable =  vec.Number_variable();
  }
  // number of variable
  int Number_variable() const { return number_variable; }

  // activity of points
  int Level() const { return level; }
  bool run_interior() const { return r_.run_interior(); }
  bool run_nearb()    const { return r_.run_nearb(); }
  int run_boundary()  const { return r_.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 not_both; }
  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {
    evpar->Variable_contained_in_expression(number_variable);
  };
  // at the and of an iteration of an expression
  void clean() const  { };
};


class Subgrid {
  friend  Subgrid operator&(Subgrid a, Subgrid b);
 private:
  Subgrid(bool r_i, bool r_n, int r_b, Grid* gr) : 
    static_grid(false), r_interior(r_i), 
    r_nearb(r_n), r_boundary(r_b), grid(gr) 
    { };
  Subgrid(bool static_gr, bool r_i, bool r_n, int r_b, Grid* gr) : 
    static_grid(static_gr), r_interior(r_i), 
    r_nearb(r_n), r_boundary(r_b), grid(gr) 
    { };

  void Free() { if(r_boundary<0) grid->Remove_number_label_bo(-r_boundary); }
  bool static_grid; 
    
 public:
  Subgrid(bool r_i, bool r_n, bool r_b) : 
    static_grid(true), r_interior(r_i), r_nearb(r_n), r_boundary(r_b), grid(NULL) { };
  Subgrid() : 
    static_grid(true), r_interior(false), r_nearb(false), r_boundary(false),
    grid(NULL) {};
  bool run_interior() const { return r_interior; }
  bool run_nearb()    const { return r_nearb; }
  int run_boundary() const { return r_boundary; }

  void operator=(const Subgrid &v) {
    r_interior = v.r_interior;
    r_nearb    = v.r_nearb;
    r_boundary = v.r_boundary;
    grid = v.grid;
  }
  
 protected:
  bool r_interior;
  bool r_nearb;
  int r_boundary;

  // This is used only when r_boundary is negative
  Grid* grid;
};

#define regular_points Subgrid(true,false,false)
#define interior_points Subgrid(true,true,false)
#define boundary_points Subgrid(false,false,true)
#define nearb_points Subgrid(false,true,false)
#define boundary_strip_points Subgrid(false,true,true)
#define all_points Subgrid(true,true,true)


class Sub_boundary_grid : public Subgrid {
 public:
  // Constructor for a new user defined Sub_boundary_grid
  Sub_boundary_grid(Grid* gr,
                    bool (*Formula)(double x,double y,double z)) : 
    Subgrid() {
    r_boundary = -gr->Give_number_label_bo();
    gr->Put_label_bo(Formula,-r_boundary);
    grid = gr;
  }
  Sub_boundary_grid(Variable& v) : Subgrid() {
    grid = v.Give_grid();
    r_boundary = -grid->Give_number_label_bo();
    grid->Put_label_bo(v.Number_variable(),-r_boundary);
  }
  Sub_boundary_grid operator~() {
    return Sub_boundary_grid(grid,r_boundary);
  }
 private: 
  // Constructor for complementary Sub_boundary_grid
  Sub_boundary_grid(Grid* gr, int r_boundary_org) : 
    Subgrid() {
    r_boundary = -gr->Give_number_label_bo();
    gr->Put_label_bo_complementary(-r_boundary_org,-r_boundary);
    grid = gr;
  }
};

//Sub_boundary_grid operator~(Sub_boundary_grid a) {

//}

Subgrid operator&(Subgrid a, Subgrid b);


/* Restrict on points */
template<class A>
DExpr<DExprResP<DExpr<A>, Subgrid> >
operator|(const  DExpr<A>& a, const Subgrid& r)
{
  typedef DExprResP<DExpr<A>, Subgrid> ExprT;
  return DExpr<ExprT>(ExprT(a,r));
}

DExpr<DExprResP<DExprLIT, Subgrid> >
operator|(double x, const Subgrid& r);

DExprResVarP<Subgrid>
operator|(Variable &a, const Subgrid& r);

FunctionClass
operator|(double (*formula)(double x,double y, double z), const Subgrid& r);





 // 3.2.:  Restriction operator with respect to points -- slow by a domain

// First part: DExprResDomain defined in domain.h

template<class A>
DExprResDomain<DExpr<A> >
operator|(const  DExpr<A>& a,  Domain* dom)
{
  typedef DExprResDomain<DExpr<A> > ExprT;
  return ExprT(ExprT(a,dom));
}

DExprResDomain<DExprLIT>
operator|(double x, const Domain& dom);

     /*
DExprResVarP<Subgrid>
operator|(Variable &a, const Subgrid& r);
     */

 // 4.:  Binary operators

/* DExprBinOp for Representing a Binary operators */
template<class A, class B, class Op>
class DExprBinOp {
  A a_;
  B b_;
 public:
  DExprBinOp(const A& a, const B& b)
    : a_(a), b_(b)
    {}
  // for evaluation of expression
  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const {
    return Op::apply(a_.Give_interior(it_i,gr,h_mesh,l, double_stencils_out),
                     b_.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 Op::apply(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));
  }
  double Give_nearb(const P_nearb* it_n,const Grid* gr, double h_mesh, int l,
                    const Nearb_Ablage* nearb_ablage)
    const {
    return Op::apply(a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage),
                     b_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage));
  };
  double Give_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    return Op::apply(a_.Give_Bo2p(it_b,gr,l),
                     b_.Give_Bo2p(it_b,gr,l));
  };
  double Give_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                      const BoCeData* bocedata) const {
    return Op::apply(a_.Give_cellpoi(it_cf,gr,bocedata),
                     b_.Give_cellpoi(it_cf,gr,bocedata));
  };

  // for evaluate on cell
  double Give_interior_cell(const P_interior_cell* it_icell ,
                            const Grid* gr, int l) const {
    return Op::apply(a_.Give_interior_cell(it_icell,gr,l),
                     b_.Give_interior_cell(it_icell,gr,l));
  };
  double Give_bo_tet(const P_boundary_tet* it_bcell, 
                      const Grid* gr) const {
    return Op::apply(a_.Give_bo_tet(it_bcell,gr),
                     b_.Give_bo_tet(it_bcell,gr));
  };

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

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

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

  // For operator ==
  void Active_Sim_Level( int lev) const { 
    a_.Active_Sim_Level(lev); b_.Active_Sim_Level(lev); };
  void Active_Sim_interior(bool run) const { 
    a_.Active_Sim_interior(run); b_.Active_Sim_interior(run); };
  void Active_Sim_nearb(bool run) const { 
    a_.Active_Sim_nearb(run); b_.Active_Sim_nearb(run); };
  void Active_Sim_boundary( int run) const {
    a_.Active_Sim_boundary(run); b_.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); };

 // 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); 
  }
 // For abstract differential operators:
  // --------------------------------------
  // Is the expression an abstract differential operator?
  Differential_op_typ Abstract_differential_operator()  const {
    Differential_op_typ a_abstract, b_abstract;
    a_abstract = a_.Abstract_differential_operator();
    b_abstract = b_.Abstract_differential_operator();
    if(a_abstract==not_abstract ||
       b_abstract==not_abstract) 
      return not_abstract;
    if(a_abstract==abstract_both)
      return b_abstract;
    if(b_abstract==abstract_both)
      return a_abstract;
    if(a_.Give_number_var_of_abstract_op() !=
       b_.Give_number_var_of_abstract_op()) return not_abstract;
    return a_abstract;

    /* old
    if(a_.Abstract_differential_operator()==not_abstract ||
       b_.Abstract_differential_operator()==not_abstract) 
      return not_abstract;
    if(a_.Abstract_differential_operator()==abstract_both)
      return b_.Abstract_differential_operator();
    if(b_.Abstract_differential_operator()==abstract_both)
      return a_.Abstract_differential_operator();
    if(a_.Give_number_var_of_abstract_op() !=
       b_.Give_number_var_of_abstract_op()) return not_abstract;
    return a_.Abstract_differential_operator();
    */
    // cout << " It is not allowed to construct such an expression! " << endl;
  }
  // For abstract differential operator:
  int Give_number_var_of_abstract_op() const { 
    return MAX(a_.Give_number_var_of_abstract_op(),
               b_.Give_number_var_of_abstract_op()); 
  };
  bool Give_array_variable_inserted() const {
    return  (a_.Give_array_variable_inserted() ||
             b_.Give_array_variable_inserted());
  }
  int Give_length_of_array_variable_inserted() const {
    return MAX(a_.Give_length_of_array_variable_inserted(),
               b_.Give_length_of_array_variable_inserted());
  }
  // For restriction of stencils
  // --------------------------------------
  double Give_interior_sten_element(Index3D I, double meshsize,
                                    int stelle,
                                    const Grid* grid, int l) const {
    return Op::apply(a_.Give_interior_sten_element(I,meshsize,stelle,grid,l),
                     b_.Give_interior_sten_element(I,meshsize,stelle,grid,l));
  }
  double Give_boundary_sten_element(const Grid* grid, BoCeData* b_cell,
                                    double* u_ablage, int num_v) const {
    return Op::apply(a_.Give_boundary_sten_element(grid,b_cell,u_ablage,num_v),
                     b_.Give_boundary_sten_element(grid,b_cell,u_ablage,num_v));
  }
  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {
    a_.Add_variables_for_parallel(evpar);
    b_.Add_variables_for_parallel(evpar);
  };
  bool GS_type(int var_number_left, int dim) const {
    if(a_.GS_type(var_number_left,dim)==true) return true;
    return b_.GS_type(var_number_left,dim); 
  };
  // at the and of an iteration of an expression
  void clean() const  {
    a_.clean(); 
    b_.clean(); 
  };
};





/* Applicative Template classes */
class DApAdd {
 public:
  DApAdd() { }
  static inline double apply(double a, double b)
    { return a+b; }
};


class DApDivide {
 public:
  DApDivide() { }
  static inline double apply(double a, double b)
    { return a/b; }
};


class DApMul {
 public:
  DApMul() { }
  static inline double apply(double a, double b)
    { return a*b; }
};

class DApSub {
 public:
  DApSub() { }
  static inline double apply(double a, double b)
    { return a-b; }
};


 // 4.:  Binary operators



/* Mul */
template<class A, class B>
DExpr<DExprBinOp<DExpr<A>, DExpr<B>, DApMul> >
operator*(const DExpr<A>& a,const  DExpr<B>& b)
{
  typedef DExprBinOp<DExpr<A>, DExpr<B>, DApMul> ExprT;
  return DExpr<ExprT>(ExprT(a,b));
}


template<class A>
DExpr<DExprBinOp<DExprLIT, DExpr<A>, DApMul> >
operator*(double x,const  DExpr<A>& a)
{
  typedef DExprBinOp<DExprLIT, DExpr<A>, DApMul> ExprT;
  return DExpr<ExprT>(ExprT(DExprLIT(x),a));
}

template<class A>
DExpr<DExprBinOp<Local_var, DExpr<A>, DApMul> >
operator*(Local_var& x,const  DExpr<A>& a)
{
  typedef DExprBinOp<Local_var, DExpr<A>, DApMul> ExprT;
  return DExpr<ExprT>(ExprT(x,a));
}

template<class A>
DExpr<DExprBinOp<DExpr<A>, DExprLIT, DApMul> >
operator*(const   DExpr<A>& a, double x)
{
  typedef DExprBinOp<DExpr<A>, DExprLIT, DApMul> ExprT;
  return DExpr<ExprT>(ExprT(a,DExprLIT(x)));
}

template<class A>
DExpr<DExprBinOp<DExpr<A>, Local_var, DApMul> >
operator*(const   DExpr<A>& a, Local_var& x)
{
  typedef DExprBinOp<DExpr<A>, Local_var, DApMul> ExprT;
  return DExpr<ExprT>(ExprT(a,x));
}

template<class A>
DExpr<DExprBinOp<DExprVAR, DExpr<A>, DApMul> >
operator*(  Variable& v,const  DExpr<A>& a)
{
  typedef DExprBinOp<DExprVAR, DExpr<A>, DApMul> ExprT;
  return DExpr<ExprT>(ExprT(DExprVAR(v),a));
}



template<class A>
DExpr<DExprBinOp<DExpr<A>, DExprVAR, DApMul> >
operator*(const  DExpr<A>& a, Variable& v)
{
  typedef DExprBinOp<DExpr<A>, DExprVAR, DApMul> ExprT;
  return DExpr<ExprT>(ExprT(a,DExprVAR(v)));
}


DExpr<DExprBinOp<DExprLIT, DExprVAR, DApMul> >
operator*(double x, Variable& a);

DExpr<DExprBinOp<Local_var, DExprVAR, DApMul> >
operator*( Local_var& x,Variable& a);


DExpr<DExprBinOp<DExprVAR, DExprLIT, DApMul> >
operator*(Variable& a,double x);

DExpr<DExprBinOp<DExprVAR, Local_var, DApMul> >
operator*(Variable& a,Local_var& x);

DExpr<DExprBinOp<DExprVAR, DExprVAR, DApMul> >
operator*(Variable& b, Variable& a);


DExpr<DExprBinOp<Local_var, DExprLIT, DApMul> >
operator*(Local_var&  a,double x);

DExpr<DExprBinOp<DExprLIT, Local_var, DApMul> >
operator*(double x, Local_var&  a);


DExpr<DExprBinOp<Local_var, Local_var, DApMul> >
operator*(Local_var& b, Local_var& a);



/* Add */
template<class A, class B>
DExpr<DExprBinOp<DExpr<A>, DExpr<B>, DApAdd> >
operator+(const DExpr<A>& a,const  DExpr<B>& b)
{
  typedef DExprBinOp<DExpr<A>, DExpr<B>, DApAdd> ExprT;
  return DExpr<ExprT>(ExprT(a,b));
}


template<class A>
DExpr<DExprBinOp<DExprLIT, DExpr<A>, DApAdd> >
operator+(double x,const  DExpr<A>& a)
{
  typedef DExprBinOp<DExprLIT, DExpr<A>, DApAdd> ExprT;
  return DExpr<ExprT>(ExprT(DExprLIT(x),a));
}

template<class A>
DExpr<DExprBinOp<Local_var, DExpr<A>, DApAdd> >
operator+(Local_var& x,const  DExpr<A>& a)
{
  typedef DExprBinOp<Local_var, DExpr<A>, DApAdd> ExprT;
  return DExpr<ExprT>(ExprT(x,a));
}

template<class A>
DExpr<DExprBinOp<DExpr<A>, DExprLIT, DApAdd> >
operator+(const   DExpr<A>& a, double x)
{
  typedef DExprBinOp<DExpr<A>, DExprLIT, DApAdd> ExprT;
  return DExpr<ExprT>(ExprT(a,DExprLIT(x)));
}

template<class A>
DExpr<DExprBinOp<DExpr<A>, Local_var, DApAdd> >
operator+(const   DExpr<A>& a, Local_var& x)
{
  typedef DExprBinOp<DExpr<A>, Local_var, DApAdd> ExprT;
  return DExpr<ExprT>(ExprT(a,x));
}

template<class A>
DExpr<DExprBinOp<DExprVAR, DExpr<A>, DApAdd> >
operator+(  Variable& v,const  DExpr<A>& a)
{
  typedef DExprBinOp<DExprVAR, DExpr<A>, DApAdd> ExprT;
  return DExpr<ExprT>(ExprT(DExprVAR(v),a));
}



template<class A>
DExpr<DExprBinOp<DExpr<A>, DExprVAR, DApAdd> >
operator+(const  DExpr<A>& a, Variable& v)
{
  typedef DExprBinOp<DExpr<A>, DExprVAR, DApAdd> ExprT;
  return DExpr<ExprT>(ExprT(a,DExprVAR(v)));
}



DExpr<DExprBinOp<DExprLIT, DExprVAR, DApAdd> >
operator+(double x,Variable& a);

DExpr<DExprBinOp<Local_var, DExprVAR, DApAdd> >
operator+( Local_var& x,Variable& a);

DExpr<DExprBinOp<DExprVAR, DExprLIT, DApAdd> >
operator+(Variable& a,double x);

DExpr<DExprBinOp<DExprVAR, Local_var, DApAdd> >
operator+(Variable& a,Local_var& x);

DExpr<DExprBinOp<DExprVAR, DExprVAR, DApAdd> >
operator+(Variable& b, Variable& a);

DExpr<DExprBinOp<Local_var, DExprLIT, DApAdd> >
operator+(Local_var&  a,double x);

DExpr<DExprBinOp<DExprLIT, Local_var, DApAdd> >
operator+(double x, Local_var&  a);

DExpr<DExprBinOp<Local_var, Local_var, DApAdd> >
operator+(Local_var& b, Local_var& a);

/* Sub */
template<class A, class B>
DExpr<DExprBinOp<DExpr<A>, DExpr<B>, DApSub> >
operator-(const DExpr<A>& a,const  DExpr<B>& b)
{
  typedef DExprBinOp<DExpr<A>, DExpr<B>, DApSub> ExprT;
  return DExpr<ExprT>(ExprT(a,b));
}


template<class A>
DExpr<DExprBinOp<DExprLIT, DExpr<A>, DApSub> >
operator-(double x,const  DExpr<A>& a)
{
  typedef DExprBinOp<DExprLIT, DExpr<A>, DApSub> ExprT;
  return DExpr<ExprT>(ExprT(DExprLIT(x),a));
}

template<class A>
DExpr<DExprBinOp<Local_var, DExpr<A>, DApSub> >
operator-(Local_var& x,const  DExpr<A>& a)
{
  typedef DExprBinOp<Local_var, DExpr<A>, DApSub> ExprT;
  return DExpr<ExprT>(ExprT(x,a));
}

template<class A>
DExpr<DExprBinOp<DExpr<A>, DExprLIT, DApSub> >
operator-(const   DExpr<A>& a, double x)
{
  typedef DExprBinOp<DExpr<A>, DExprLIT, DApSub> ExprT;
  return DExpr<ExprT>(ExprT(a,DExprLIT(x)));
}

template<class A>
DExpr<DExprBinOp<DExpr<A>, Local_var, DApSub> >
operator-(const   DExpr<A>& a, Local_var& x)
{
  typedef DExprBinOp<DExpr<A>, Local_var, DApSub> ExprT;
  return DExpr<ExprT>(ExprT(a,x));
}

template<class A>
DExpr<DExprBinOp<DExprVAR, DExpr<A>, DApSub> >
operator-(  Variable& v,const  DExpr<A>& a)
{
  typedef DExprBinOp<DExprVAR, DExpr<A>, DApSub> ExprT;
  return DExpr<ExprT>(ExprT(DExprVAR(v),a));
}




template<class A>
DExpr<DExprBinOp<DExpr<A>, DExprVAR, DApSub> >
operator-(const  DExpr<A>& a, Variable& v)
{
  typedef DExprBinOp<DExpr<A>, DExprVAR, DApSub> ExprT;
  return DExpr<ExprT>(ExprT(a,DExprVAR(v)));
}



DExpr<DExprBinOp<DExprLIT, DExprVAR, DApSub> >
operator-(double x,Variable& a);

DExpr<DExprBinOp<Local_var, DExprVAR, DApSub> >
operator-( Local_var& x,Variable& a);

DExpr<DExprBinOp<DExprVAR, DExprLIT, DApSub> >
operator-(Variable& a,double x);


DExpr<DExprBinOp<DExprVAR, Local_var, DApSub> >
operator-(Variable& a,Local_var& x);

DExpr<DExprBinOp<DExprVAR, DExprVAR, DApSub> >
operator-(Variable& b, Variable& a);

DExpr<DExprBinOp<Local_var, DExprLIT, DApSub> >
operator-(Local_var&  a,double x);


DExpr<DExprBinOp<DExprLIT, Local_var, DApSub> >
operator-(double x, Local_var&  a);

DExpr<DExprBinOp<Local_var, Local_var, DApSub> >
operator-(Local_var& b, Local_var& a);


/* Divide */
template<class A, class B>
DExpr<DExprBinOp<DExpr<A>, DExpr<B>, DApDivide> >
operator/(const DExpr<A>& a,const  DExpr<B>& b)
{
  typedef DExprBinOp<DExpr<A>, DExpr<B>, DApDivide> ExprT;
  return DExpr<ExprT>(ExprT(a,b));
}


template<class A>
DExpr<DExprBinOp<DExprLIT, DExpr<A>, DApDivide> >
operator/(double x,const  DExpr<A>& a)
{
  typedef DExprBinOp<DExprLIT, DExpr<A>, DApDivide> ExprT;
  return DExpr<ExprT>(ExprT(DExprLIT(x),a));
}

template<class A>
DExpr<DExprBinOp<Local_var, DExpr<A>, DApDivide> >
operator/(Local_var& x,const  DExpr<A>& a)
{
  typedef DExprBinOp<Local_var, DExpr<A>, DApDivide> ExprT;
  return DExpr<ExprT>(ExprT(x,a));
}

template<class A>
DExpr<DExprBinOp<DExpr<A>, DExprLIT, DApDivide> >
operator/(const   DExpr<A>& a, double x)
{
  typedef DExprBinOp<DExpr<A>, DExprLIT, DApDivide> ExprT;
  return DExpr<ExprT>(ExprT(a,DExprLIT(x)));
}

template<class A>
DExpr<DExprBinOp<DExpr<A>, Local_var, DApDivide> >
operator/(const   DExpr<A>& a, Local_var& x)
{
  typedef DExprBinOp<DExpr<A>, Local_var, DApDivide> ExprT;
  return DExpr<ExprT>(ExprT(a,x));
}

template<class A>
DExpr<DExprBinOp<DExprVAR, DExpr<A>, DApDivide> >
operator/(  Variable& v,const  DExpr<A>& a)
{
  typedef DExprBinOp<DExprVAR, DExpr<A>, DApDivide> ExprT;
  return DExpr<ExprT>(ExprT(DExprVAR(v),a));
}


template<class A>
DExpr<DExprBinOp<DExpr<A>, DExprVAR, DApDivide> >
operator/(const  DExpr<A>& a, Variable& v)
{
  typedef DExprBinOp<DExpr<A>, DExprVAR, DApDivide> ExprT;
  return DExpr<ExprT>(ExprT(a,DExprVAR(v)));
}


DExpr<DExprBinOp<DExprLIT, DExprVAR, DApDivide> >
operator/(double x,Variable& a);

DExpr<DExprBinOp<Local_var, DExprVAR, DApDivide> >
operator/( Local_var& x,Variable& a);


DExpr<DExprBinOp<DExprVAR, DExprLIT, DApDivide> >
operator/(Variable& a,double x);


DExpr<DExprBinOp<DExprVAR, Local_var, DApDivide> >
operator/(Variable& a,Local_var& x);

DExpr<DExprBinOp<DExprVAR, DExprVAR, DApDivide> >
operator/(Variable& b, Variable& a);

DExpr<DExprBinOp<Local_var, DExprLIT, DApDivide> >
operator/(Local_var&  a,double x);


DExpr<DExprBinOp<DExprLIT, Local_var, DApDivide> >
operator/(double x, Local_var&  a);


DExpr<DExprBinOp<Local_var, Local_var, DApDivide> >
operator/(Local_var& b, Local_var& a);


 // 5.:  Unary operator

/* DExprBinOp for Representing a Binary operators */
template<class A>
class DExprMinusUnaOp {
  A a_;
 public:
  DExprMinusUnaOp(const A& a)
    : a_(a)
    {}
  // for evaluation of expression
  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 {
    return - a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage);
  };
  double Give_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    return - a_.Give_Bo2p(it_b,gr,l);
  };
  double Give_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                      const BoCeData* bocedata) const {
    return - a_.Give_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() + 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);
  }
 // For abstract differential operators:
  // --------------------------------------
  // Is the expression an abstract differential operator?
  Differential_op_typ Abstract_differential_operator()  const {
    return a_.Abstract_differential_operator();
  }
  // For abstract differential operator:
  int Give_number_var_of_abstract_op() const { 
    return a_.Give_number_var_of_abstract_op();
  };
  // For restriction of stencils
  // --------------------------------------
  double Give_interior_sten_element(Index3D I, double meshsize,
                                    int stelle,
                                    const Grid* grid, int level) const {
    return - a_.Give_interior_sten_element(I,meshsize,stelle,grid,level);
  }
  double Give_boundary_sten_element(const Grid* grid, BoCeData* b_cell,
                                    double* u_ablage, int num_v) const {
    return - a_.Give_boundary_sten_element(grid,b_cell,u_ablage,num_v);
  }
  // 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) const {
    return a_.GS_type(var_number_left,dim); 
  };
  void clean() const  { a_.clean(); }
};



/* Sub */
template<class A>
DExpr<DExprMinusUnaOp<DExpr<A> > >
operator-(const DExpr<A>& a)
{
  typedef DExprMinusUnaOp<DExpr<A> > ExprT;
  return DExpr<ExprT>(ExprT(a));
}


DExpr<DExprMinusUnaOp<Local_var> >
operator-(Local_var& x);

DExpr<DExprMinusUnaOp<DExprVAR> >
operator-(Variable& v);



#endif
      

---