src/expde/extemp/opera.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 OPERA_H_
#define OPERA_H_
  
#include "grid/SSten.h"
        
// ------------------------------------------------------------
// operator.h
//
// ------------------------------------------------------------

// Remark: this file uses the formulas in source/formulas/diffopc.h

//   Constant coefficient Differential operators

 
// a.: Operator




/* DExprDiff_15S_Op for Differential operators */
template<class V, class DiffOp>
class DVarDiff_15S_Op {
  V v_;
  double *Stencil;      // stencil independent of meshsize

 public:
  // Constructor for concrete operators
  DVarDiff_15S_Op(const V& v)
    : v_(v) {
    DiffOp opp;
    if(opp.typ_u()==Poisson_type ||
       opp.typ_v()==Poisson_type) 
      Stencil = v.Give_grid()->Give_Loc_stencil_Poisson();
    else {
      Stencil = v.Give_grid()->Give_Loc_stencil(opp.typ_u(),opp.typ_v());
    }
  }
  DVarDiff_15S_Op(const DExpr<V>& v)
    : v_(v) {
    DiffOp opp;
    if(opp.typ_u()==Poisson_type ||
       opp.typ_v()==Poisson_type) 
      Stencil = v.Give_grid()->Give_Loc_stencil_Poisson();
    else {
      Stencil = v.Give_grid()->Give_Loc_stencil(opp.typ_u(),opp.typ_v());
    }
  }
  double Give_interior(const P_interior* it_i, const Grid* grid, double h_mesh,
                       int lev, double_stencils_in) const {
    return Give_interior_here(it_i,grid,h_mesh,lev, double_stencils_out);
  }
  double Give_interior_coarse(const P_interior* it_i, const Grid* grid,
                              double h_mesh,
                              int lev, double_stencils_in) const {
    return Give_interior_here(it_i,grid,h_mesh,lev, double_stencils_out);
  }
  double Give_interior_here(const P_interior* it_i, const Grid* gr,
                            double h_mesh,
                            int l, double_stencils_in) const {
   return DiffOp::stencil(LuN[v_.Number_variable()],
                          LuW[v_.Number_variable()],
                          LuM[v_.Number_variable()],
                          LuE[v_.Number_variable()],
                          LuS[v_.Number_variable()],
                          LuT[v_.Number_variable()],
                          LuD[v_.Number_variable()],
                          LuND[v_.Number_variable()],
                          LuWN[v_.Number_variable()],
                          LuWT[v_.Number_variable()],
                          LuED[v_.Number_variable()],
                          LuST[v_.Number_variable()],
                          LuES[v_.Number_variable()],
                          LuEST[v_.Number_variable()],
                          LuWND[v_.Number_variable()],
                          h_mesh);
  };
  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 BoCeData* bocedata;
    static Celltype cetyp;
    static double sum;
    static int num_v;
    dir_sons dir_v;

    sum = 0.0;

    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);

      if(cetyp==fine_bo_cell) {
        bocedata = nearb_ablage->Give_BoData((dir_sons)i);
        // Berechnung der Nummer der Testfunktion auf i-Randzelle
        for(num_v=0;(!(bocedata->corner(num_v)==dir_v)) ||
              bocedata->edge_point(num_v);) { ++num_v; }
        // Auswertung
        sum = sum + DiffOp::F_bo_cell(bocedata,num_v,v_.Number_variable());
      }
      else if(cetyp==int_cell) {
        sum = sum +
          DiffOp::F_reg_cell(nearb_ablage->Give_u_Recell((dir_sons)i),
                             v_.Number_variable(),
                             dir_v,h_mesh,Stencil);
        //      if(now) cout << " sum r " << sum << endl;
      }
    }
    return sum;
  };
  double Give_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
                      const BoCeData* bocedata) const {
    // Auswertung 
    return DiffOp::F_bo_cell(bocedata,bocedata->Give_number_points(),
                             v_.Number_variable());
  };
  double Give_Bo2p(const P_Bo2p* it_b, const Grid* gr,int l)  const {
    static int i;
    static BoCeData* bocedata;
    static double sum;
    static int num_v;
    dir_sons dir_v;
    bool stop;


    // return 0.0;

   sum = 0.0;
    for(i=0;i<8;++i) {    // Durchlaufen aller Nachbarzellen
      dir_v = opposite3D((dir_sons)i);
      // Info ueber i-Randzelle
      bocedata = it_b->Give_BoData((dir_sons)i); 
      if(bocedata!=NULL) {

        // Berechnung der Nummer der Testfunktion auf i-Randzelle
        stop = false;
        for(num_v=0;!stop && bocedata->Give_number_points() > num_v;) {
          stop = bocedata->corner(num_v)   == dir_v      &&
                 bocedata->edge_dir(num_v) == it_b->d()  &&
                 bocedata->edge_point(num_v);
          if(!stop) ++num_v;
        }

        if(bocedata->Give_number_points() > num_v) {
          // Auswertung
            sum = sum + 
              DiffOp::F_bo_cell(bocedata,num_v,v_.Number_variable());
        }
      }
    }

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

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

  // number of operations
  int ops_interior() const { return DiffOp::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(const 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_with_var; }
  int Give_number_var_of_abstract_op() const { 
    return v_.Number_variable_start();
  }
  bool Give_array_variable_inserted() const {
    return v_.Give_array_variable_inserted();
  }
  int Give_length_of_array_variable_inserted() const {
    return v_.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 level) const {
    return Stencil[stelle] * DiffOp::Factor_reg_cell(meshsize);
  }
  double Give_boundary_sten_element(const Grid* grid, BoCeData* b_cell,
                                    double* u_ablage, int num_v) const {
    return DiffOp::Local_matrix_bo_cell(b_cell,u_ablage,num_v);
  }

  // for Parallelization
  void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar)
    const {
    evpar->Variable_contained_in_expression(v_.Number_variable());
  }
  bool GS_type(int var_number_left, int dim) const {
    if(var_number_left<=v_.Number_variable() &&
       var_number_left+dim>=v_.Number_variable())
      return true;
    return false; 
  };

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


/* Unary operators */
// Laplace

class DExprVAR_ARR;

#define Macro_operator_definition(insertA) \
DExpr<DVarDiff_15S_Op<insertA, laplace_FE_const> > \
Laplace_FE(const insertA& a); \
/* DxDx */ \
DExpr<DVarDiff_15S_Op<insertA, dxdx_FE_const> > \
DxDx_FE(insertA a); \
/* Dydy */ \
DExpr<DVarDiff_15S_Op<insertA, dydy_FE_const> > \
DyDy_FE(insertA a); \
/* Dzdz */ \
DExpr<DVarDiff_15S_Op<insertA, dzdz_FE_const> > \
DzDz_FE(insertA a); \
/* Helm */ \
DExpr<DVarDiff_15S_Op<insertA, helm_FE_const> > \
Helm_FE(insertA a); \
/* ///////////////////////// \
// Nichsymmetrische Terme \
/* Dxdy */ \
DExpr<DVarDiff_15S_Op<insertA, dxdy_FE_const> > \
DxDy_FE(insertA a); \
/* Dydx */ \
DExpr<DVarDiff_15S_Op<insertA, dydx_FE_const> > \
DyDx_FE(insertA a); \
/* Dxdz */ \
DExpr<DVarDiff_15S_Op<insertA, dxdz_FE_const> > \
DxDz_FE(insertA a); \
/* Dzdx */ \
DExpr<DVarDiff_15S_Op<insertA, dzdx_FE_const> > \
DzDx_FE(insertA a); \
/* Dydz */ \
DExpr<DVarDiff_15S_Op<insertA, dydz_FE_const> > \
DyDz_FE(insertA a); \
/* Dzdy */ \
DExpr<DVarDiff_15S_Op<insertA, dzdy_FE_const> > \
DzDy_FE(insertA a); \
DExpr<DVarDiff_15S_Op<insertA, dxhelm_FE_const> > \
DX_FE(insertA a); \
DExpr<DVarDiff_15S_Op<insertA, helmdx_FE_const> > \
DX_FE_t(insertA a); \
DExpr<DVarDiff_15S_Op<insertA, dyhelm_FE_const> > \
DY_FE(insertA a); \
DExpr<DVarDiff_15S_Op<insertA, helmdy_FE_const> > \
DY_FE_t(insertA a); \
DExpr<DVarDiff_15S_Op<insertA, dzhelm_FE_const> > \
DZ_FE(insertA a); \
DExpr<DVarDiff_15S_Op<insertA, helmdz_FE_const> > \
DZ_FE_t(insertA a); \
/* ///////////////////////// \
// Rand-operatoren \
DExpr<DVarDiff_15S_Op<insertA, L2boundary_const> > \
Int_boundary(insertA a);

Macro_operator_definition(Variable)
// other usage of Macro_operator_definition in array.h
 

// b.: Diagonal operator


// We need this, since expression is const!
class Storage_for_one_Stencil {
  public :
    double *Stencil;
};


/* DExprDiff_15S_Op for Differential operators */
template<class DiffOp>
class DiagonalDiff_Op {
 private:
  Grid* grid;
  Storage_for_one_Stencil* OS; 
 public:
  DiagonalDiff_Op()  {
    OS = new Storage_for_one_Stencil;
  }
  double Give_interior(const P_interior* it_i, const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const {
    return DiffOp::diag_stencil(h_mesh);
  };
  double Give_interior_coarse(const P_interior* it_i, const Grid* gr,
                              double h_mesh,
                              int l, double_stencils_in) const {
    return DiffOp::diag_stencil(h_mesh);
  };
  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 BoCeData* bocedata;
    static double sum;
    static Celltype cetyp;
    static int num_v;
    dir_sons dir_v;

    sum = 0.0;
    for(i=0;i<8;++i) {    // Durchlaufen aller Nachbarzellen
      dir_v = opposite3D((dir_sons)i);
      cetyp = nearb_ablage->Give_Celltyp((dir_sons)i);
      if(cetyp==fine_bo_cell) {
        // Hole Info ueber i-Randzelle
        bocedata = nearb_ablage->Give_BoData((dir_sons)i); 
        // Berechnung der Nummer der Testfunktion auf i-Randzelle
        for(num_v=0;(!(bocedata->corner(num_v)==dir_v)) ||
              bocedata->edge_point(num_v);) { ++num_v; }
        // Auswertung
        sum = sum + DiffOp::diag_F_bo_cell(bocedata,num_v);     
      }
      else if(cetyp==int_cell) {
        sum = sum +
          DiffOp::diag_F_reg_cell(dir_v,h_mesh,OS->Stencil);
      }
    }
    return sum;
  };
  double Give_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
                      const BoCeData* bocedata) const {
    // Auswertung 
    return DiffOp::diag_F_bo_cell(bocedata,bocedata->Give_number_points());
  };
  double Give_Bo2p(const P_Bo2p* it_b, const Grid* gr,int l)  const {
    static int i;
    static BoCeData* bocedata;
    static double sum;
    static int num_v;
    dir_sons dir_v;
    bool stop;

    sum = 0.0;
    for(i=0;i<8;++i) {    // Durchlaufen aller Nachbarzellen
      dir_v = opposite3D((dir_sons)i);
      // Info ueber i-Randzelle
      //  wegg bo:    bocedata = bo2p_ablage->Give_BoData((dir_sons)i); 
      bocedata = it_b->Give_BoData((dir_sons)i); 
      if(bocedata!=NULL) {
        // Berechnung der Nummer der Testfunktion auf i-Randzelle
        stop = false;
        for(num_v=0;!stop && bocedata->Give_number_points() > num_v;) {
          stop = bocedata->corner(num_v)   == dir_v      &&
                 bocedata->edge_dir(num_v) == it_b->d()  &&
                 bocedata->edge_point(num_v);
          if(!stop) ++num_v;
        }

        if(bocedata->Give_number_points() > num_v) {
          // Auswertung
          sum = sum + 
            DiffOp::diag_F_bo_cell(bocedata,num_v);
        }
      }
    }
    return sum;
  };
  // Size of necessary stencil
  int Sice_stencil() const { return 15; }

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

  // number of operations
  int ops_interior() const { return DiffOp::ops_diag_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 { 
    DiffOp opp;
    if(opp.typ_u()==Poisson_type ||
       opp.typ_v()==Poisson_type) 
      OS->Stencil = gr->Give_Loc_stencil_Poisson();
    else {
      OS->Stencil = gr->Give_Loc_stencil(opp.typ_u(),opp.typ_v());
    }
  };
  // Is the expression an abstract differential operator?
  Differential_op_typ Abstract_differential_operator() const
    { return abstract_diag; }
  int Give_number_var_of_abstract_op() const { return -1; }
  // 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 { delete(OS); };
};


/* Unary operators */

// Laplace
DExpr<DiagonalDiff_Op<laplace_FE_const> >
Diag_Laplace_FE();

// DxDx
DExpr<DiagonalDiff_Op<dxdx_FE_const> >
Diag_DxDx_FE();

// DyDy
DExpr<DiagonalDiff_Op<dydy_FE_const> >
Diag_DyDy_FE();

// DzDz
DExpr<DiagonalDiff_Op<dzdz_FE_const> >
Diag_DzDz_FE();

// Helm
DExpr<DiagonalDiff_Op<helm_FE_const> >
Diag_Helm_FE();

// Nichsymmetrische Terme

// DxDy
DExpr<DiagonalDiff_Op<dxdy_FE_const> >
Diag_DxDy_FE();

// DyDx
DExpr<DiagonalDiff_Op<dydx_FE_const> >
Diag_DyDx_FE();

// DxDz
DExpr<DiagonalDiff_Op<dxdz_FE_const> >
Diag_DxDz_FE();

// DzDx
DExpr<DiagonalDiff_Op<dzdx_FE_const> >
Diag_DzDx_FE();

// DyDz
DExpr<DiagonalDiff_Op<dydz_FE_const> >
Diag_DyDz_FE();

// DzDy
DExpr<DiagonalDiff_Op<dzdy_FE_const> >
Diag_DzDy_FE();

// DX
DExpr<DiagonalDiff_Op<dxhelm_FE_const> >
Diag_DX_FE();

// DX transpose
DExpr<DiagonalDiff_Op<helmdx_FE_const> >
Diag_DX_FE_t();

// DY
DExpr<DiagonalDiff_Op<dyhelm_FE_const> >
Diag_DY_FE();

// DY transpose
DExpr<DiagonalDiff_Op<helmdy_FE_const> >
Diag_DY_FE_t();

// DZ
DExpr<DiagonalDiff_Op<dzhelm_FE_const> >
Diag_DZ_FE();

// DZ transpose
DExpr<DiagonalDiff_Op<helmdz_FE_const> >
Diag_DZ_FE_t();


// Rand-operatoren

DExpr<DiagonalDiff_Op<L2boundary_const> >
Diag_Int_boundary();

#endif
        

---