src/expde/extemp/varcon.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 VARCON_H_
#define VARCON_H_
           
// ------------------------------------------------------------
//
//     varcon.h 
//
// ------------------------------------------------------------

//
// coefficient with constant variable on cells
//


// 1. first wrapper class
/* wrapper class wich contains  ...contains the expression types 
               - ConExprVAR
               - ConExprLIT 
               - ConExprBinOp 
               - ...           */


template<class A>
class CoExpr {
 private:
  A a_;

 public:
  CoExpr(const A& x)
    : a_(x)
    {}
  // for evaluate on cell
  double Give_interior_cell(const P_interior_cell* it_icell ,
                            const Grid* gr, int l) const {
    return a_.Give_interior_cell(it_icell,gr ,l);
  };
  double Give_bo_tet(const P_boundary_tet* it_bcell, 
                      const Grid* gr) const {
    return a_.Give_bo_tet(it_bcell,gr);
  };

  // activity of points
  int Level() const { return a_.Level(); }
  Dominace_label Dominant_lev() const { return a_.Dominant_lev(); }

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

  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {
    a_.Add_variables_for_parallel(evpar);
  };
};


// 2. variable

/* class Cell_Variable */
class Cell_Variable {
 private:
  //  public:
  // Aktive Punkte:
  int  level;        // level:

  // Iterator
  P_interior_cell  *iter_ice;
  P_boundary_tet  *iter_bce;

  // Zwischenspeicher

  // Gitter
  Grid *grid;

  // Nummer der Variablen
  int number_cell_variable;

 public:
  Cell_Variable(Grid* gridp);

  int my_rank() { return grid->Give_my_rank(); };

  // Zur Speicherverwaltung:
  inline void Test_init() { grid->Test_init(); }

  int Number_cell_variable() const { return number_cell_variable; }

  // activity of level
  inline int  Level() const { return level; }
  inline void Level_up();
  inline void Level_down();
  inline void Active_Level(int Level);

  // activity of points
  inline int Max_Level() const { return grid->Max_level(); }

  inline Grid* Give_grid() const {
    if(developer_version) {
      if(grid==NULL) cout << "\n Variable is not initialized! Error!" << endl;
    }
    return grid;
  }

  // for assign
  template<class A>
    //  void operator=(const CellExpr<A>& a_);
  void operator=(const Cell_Variable &v);
  void operator=(double (*Formula)(double x,double y,double z));
  void operator=(double x);
  void operator=(Input_data_object& input_object);
  template<class A>
  void operator=(const CoExpr<A>& a_);
  template<class A>
  void operator=(const DExpr<A>& a_);

  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {  };
};

/* DExprVAR for Variable */

class CoExprVAR {
 private:
  // Aktive Punkte:
  int  level;        // level:

  // Nummer der Variablen
  int number_cell_variable;
 public:
  CoExprVAR(const Cell_Variable& vec) {
    level      = vec.Level();
    number_cell_variable =  vec.Number_cell_variable();
  }
  // for evaluate on cell
  double Give_interior_cell(const P_interior_cell* it_icell ,
                            const Grid* gr, int l) const {
    return it_icell->varM(gr)[number_cell_variable];
  };
  double Give_bo_tet(const P_boundary_tet* it_bcell, 
                      const Grid* gr) const {
    return it_bcell->varM(gr)[number_cell_variable];
  };

  // activity of points
  int Level() const { return level; }
  Dominace_label Dominant_lev() const { return not_dominant; }

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

  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {  };
};



// 3.:  double


/* DExprLIT for double */
class CoExprLIT {
 private:
  double value_;
 public:
  CoExprLIT(double value)
    { value_ = value; }
  // for evaluate on cell
  double Give_interior_cell(const P_interior_cell* it_icell ,
                            const Grid* gr, int l) const {
    return value_;
  };
  double Give_bo_tet(const P_boundary_tet* it_bcell, 
                      const Grid* gr) const {
    return value_;
  };

  // activity of points
  int Level() const { return 0; }
  Dominace_label Dominant_lev() const { return anti_dominant; }

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

  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {  };
};



// 4.:  Binary operators

/* DExprBinOp for Representing a Binary operators */
template<class A, class B, class Op>
class CoExprBinOp {
  A a_;
  B b_;
 public:
  CoExprBinOp(const A& a, const B& b)
    : a_(a), b_(b)
    {}
  // 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));
  };

  // 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()));
  }

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

  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {
    a_.Add_variables_for_parallel(evpar);
    b_.Add_variables_for_parallel(evpar);
  };
};

// 5.: transfer function

class CoExprTrans {
 private:
  // Aktive Punkte:
  int  level;        // level:

  // Nummer der Variablen
  int number_variable;
 public:
  CoExprTrans(const Variable& vec) {
    level      = vec.Level();
    number_variable =  vec.Number_variable();
  }
  // for evaluate on cell
  double Give_interior_cell(const P_interior_cell* it_icell ,
                            const Grid* gr, int l) const {
    return (it_icell->varENT(gr,l)[number_variable] +
            it_icell->varEND(gr,l)[number_variable] +
            it_icell->varEST(gr,l)[number_variable] +
            it_icell->varESD(gr,l)[number_variable] +
            it_icell->varWNT(gr,l)[number_variable] +
            it_icell->varWND(gr,l)[number_variable] +
            it_icell->varWST(gr,l)[number_variable] +
            it_icell->varWSD(gr,l)[number_variable])   * 0.125;
  };
  double Give_bo_tet(const P_boundary_tet* it_bcell, 
                      const Grid* gr) const {
    return (it_bcell->var0(gr)[number_variable] +
            it_bcell->var1(gr)[number_variable] +
            it_bcell->var2(gr)[number_variable] +
            it_bcell->var3(gr)[number_variable])       * 0.25;
  };

  // activity of points
  int Level() const { return level; }
  Dominace_label Dominant_lev() const { return not_dominant; }

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

  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar) 
    const {
    evpar->Variable_contained_in_expression(number_variable);
  };
};

CoExpr<CoExprTrans> Cell_Interpolation(const Variable& vec);

// 6.: Some functions

void Print_UCD(Cell_Variable& a,ofstream *Datei);
void Print_UCD_moved(Cell_Variable& v,
                     Variable& a, Variable& b, Variable& c, ofstream *Datei);
void Print_UCD_parallel(Cell_Variable& a,ofstream *Datei);

double Maximum(Cell_Variable& a);
double Minimum(Cell_Variable& a);
double L_infty(Cell_Variable& a);

// ------------------------------------------------------------
//  implementierung of the inline member functions of cell variable
// ------------------------------------------------------------

inline Cell_Variable::Cell_Variable(Grid* gridp) {
  grid = gridp;
  number_cell_variable = grid->Give_number_cell_variable();
  // all points of finest level become active
  level = gridp->Max_level();
  if(gridp->Is_storage_initialized()==true) {
    cout << " \n Definition of cell_variable not possible! ";
  }
}

inline void Cell_Variable::Level_up()   { 
  level++; 
}
inline void Cell_Variable::Level_down() { 
  level--; 
}
inline void Cell_Variable::Active_Level(int Level) { 
  level = Level; 
}

// for assign; usual
template<class A>
void Cell_Variable::operator=(const CoExpr<A>& a_) {
  int lev;
  if(grid->I_am_active()) {
    Test_init();  // Grid storage initialization
  }

  // 1.Step: level
  if(a_.Dominant_lev()==anti_dominant) lev = level;
  else                                 lev = a_.Level();
  if(lev < 0)          { cout << "\n Level too coarse! "   << endl; }
  else if (lev > Max_Level())  { cout << "\n Level too fine! " << endl; }
  
  Active_Level(lev);

  if(grid->I_am_active()) {    
    // 2.Step: Evaluation_Parallelization
    Evaluation_Parallelization_object* ev_par_obj;

    if(parallel_version) {
      ev_par_obj = grid->Give_eval_par();
      ev_par_obj->StartA_varcon_evaluation(lev);
      a_.Add_variables_for_parallel(ev_par_obj);  // this adds the variables
      // for communication
      ev_par_obj->StartB_varcon_evaluation();     // makes parallelization
    }
    // 3.Step: interior_cells
    for(iter_ice=grid->Start_P_interior_cell(level);
        iter_ice!=NULL;iter_ice=iter_ice->Next()) {
      iter_ice->varM(grid)[number_cell_variable] = 
        a_.Give_interior_cell(iter_ice,grid, level);
    }
    // 4.Step: boundary cells
    for(iter_bce=grid->Start_P_boundary_tet(level);
        iter_bce!=NULL;iter_bce=iter_bce->Next()) {
      iter_bce->varM(grid)[number_cell_variable] =
        a_.Give_bo_tet(iter_bce,grid);
    }
  }
}

// for assign; usual
template<class A>
void Cell_Variable::operator=(const DExpr<A>& a_) {
  int lev;
  Test_init();  // Grid storage initialization

  // level
  if(a_.Dominant_lev()==anti_dominant) lev = level;
  else                                 lev = a_.Level();
  if(lev < 0)          { cout << "\n Level too coarse! "   << endl; }
  else if (lev > Max_Level())  { cout << "\n Level too fine! " << endl; }

  Active_Level(lev);

  // interior_cells
  for(iter_ice=grid->Start_P_interior_cell(level);
      iter_ice!=NULL;iter_ice=iter_ice->Next()) {
    iter_ice->varM(grid)[number_cell_variable] = 
      a_.Give_interior_cell(iter_ice,grid, level);
  }
  // boundary cells
  for(iter_bce=grid->Start_P_boundary_tet(level);
      iter_bce!=NULL;iter_bce=iter_bce->Next()) {
    iter_bce->varM(grid)[number_cell_variable] =
      a_.Give_bo_tet(iter_bce,grid);
  }
}


// ------------------------------------------------------------
//  Addition, Multi, ...
// ------------------------------------------------------------


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


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


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

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


#define Macro_expressions_for_cells(insertA,insertB) \
template<class A, class B> \
CoExpr<CoExprBinOp<CoExpr<A>, CoExpr<B>, insertA> > \
insertB(const CoExpr<A>& a,const  CoExpr<B>& b) \
{ \
  typedef CoExprBinOp<CoExpr<A>, CoExpr<B>, insertA> ExprT; \
  return CoExpr<ExprT>(ExprT(a,b)); \
} \
CoExpr<CoExprBinOp<CoExprLIT, CoExprVAR, insertA> > \
insertB(double x, Cell_Variable& a); \
template<class A> \
CoExpr<CoExprBinOp<CoExprLIT, CoExpr<A>, insertA> > \
insertB(double x,const  CoExpr<A>& a) \
{ \
  typedef CoExprBinOp<CoExprLIT, CoExpr<A>, insertA> ExprT; \
  return CoExpr<ExprT>(ExprT(CoExprLIT(x),a)); \
} \
template<class A> \
CoExpr<CoExprBinOp<CoExpr<A>, CoExprLIT, insertA> > \
insertB(const   CoExpr<A>& a, double x) \
{ \
  typedef CoExprBinOp<CoExpr<A>, CoExprLIT, insertA> ExprT; \
  return CoExpr<ExprT>(ExprT(a,CoExprLIT(x))); \
} \
template<class A> \
CoExpr<CoExprBinOp<CoExprVAR, CoExpr<A>, insertA> > \
insertB(  Cell_Variable& v,const  CoExpr<A>& a) \
{ \
  typedef CoExprBinOp<CoExprVAR, CoExpr<A>, insertA> ExprT; \
  return CoExpr<ExprT>(ExprT(CoExprVAR(v),a)); \
} \
template<class A> \
CoExpr<CoExprBinOp<CoExpr<A>, CoExprVAR, insertA> > \
insertB(const  CoExpr<A>& a, Cell_Variable& v) \
{ \
  typedef CoExprBinOp<CoExpr<A>, CoExprVAR, insertA> ExprT; \
  return CoExpr<ExprT>(ExprT(a,CoExprVAR(v))); \
} \
CoExpr<CoExprBinOp<CoExprVAR, CoExprLIT, insertA> > \
insertB(Cell_Variable& a,double x); \
CoExpr<CoExprBinOp<CoExprVAR, CoExprVAR, insertA> > \
insertB(Cell_Variable& b, Cell_Variable& a); \
template<class A, class B> \
CoExpr<CoExprBinOp<DExpr<A>, CoExpr<B>, insertA> > \
insertB(const DExpr<A>& a,const  CoExpr<B>& b) \
{ \
  typedef CoExprBinOp<DExpr<A>, CoExpr<B>, insertA> ExprT; \
  return CoExpr<ExprT>(ExprT(a,b)); \
} \
template<class A, class B> \
CoExpr<CoExprBinOp<CoExpr<A>, DExpr<B>, insertA> > \
insertB(const CoExpr<A>& a,const  DExpr<B>& b) \
{ \
  typedef CoExprBinOp<CoExpr<A>, DExpr<B>, insertA> ExprT; \
  return CoExpr<ExprT>(ExprT(a,b)); \
} \
template<class A> \
CoExpr<CoExprBinOp<CoExprVAR, DExpr<A>, insertA> > \
insertB(  Cell_Variable& v,const  DExpr<A>& a) \
{ \
  typedef CoExprBinOp<CoExprVAR, DExpr<A>, insertA> ExprT; \
  return CoExpr<ExprT>(ExprT(CoExprVAR(v),a)); \
} \
template<class A> \
CoExpr<CoExprBinOp<DExpr<A>, CoExprVAR, insertA> > \
insertB(const  DExpr<A>& a, Cell_Variable& v) \
{ \
  typedef CoExprBinOp<DExpr<A>, CoExprVAR, insertA> ExprT; \
  return CoExpr<ExprT>(ExprT(a,CoExprVAR(v))); \
}



Macro_expressions_for_cells(CoApMul,operator*)
Macro_expressions_for_cells(CoApAdd,operator+)
Macro_expressions_for_cells(CoApDivide,operator/)
Macro_expressions_for_cells(CoApSub,operator-)

#endif
              

---