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

// ------------------------------------------------------------
// extra.h
//
// ------------------------------------------------------------


#ifndef EXTRA_H_
#define EXTRA_H_
         

class DVar_Extra_Op {
 private:
  Variable v_;
  double Give_at_Index(const Index3D I, const Grid* grid) const;

 public:
  DVar_Extra_Op(const Variable& v) : v_(v) {};
  double Give_interior(const P_interior* it_i,const Grid* gr, double h_mesh,
                       int l, double_stencils_in) const {
    return v_.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 v_.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 Give_at_Index(it_n->Ind(),gr);
  };

  /*
  double Give_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    Index3D In, Inn, Innn;
    dir_3D dd;
    double un, unn, unnn, h_mesh, h_d;
    Pointtype Ptypnn, Ptypnnn;

    In = it_b->Ind();
    dd = it_b->d();
    if(dd==Wdir) { Inn = In.next_E(l); Innn = Inn.next_E(l); }
    if(dd==Edir) { Inn = In.next_W(l); Innn = Inn.next_W(l); }

    if(dd==Ndir) { Inn = In.next_S(l); Innn = Inn.next_S(l); }
    if(dd==Sdir) { Inn = In.next_N(l); Innn = Inn.next_N(l); }

    if(dd==Tdir) { Inn = In.next_D(l); Innn = Inn.next_D(l); }
    if(dd==Ddir) { Inn = In.next_T(l); Innn = Inn.next_T(l); }

    if(gr->Point_in_domain(Inn) && gr->Point_in_domain(Innn)) {
      Ptypnn  = gr->Give_type(Inn);
      Ptypnnn = gr->Give_type(Innn);
      un = Give_at_Index(In,gr);
      if(interior == Ptypnn || uniform <= Ptypnn)
        unn = gr->Give_variable(Inn,l)[v_.Number_variable()];
      else
        unn = Give_at_Index(Inn,gr);
      if(interior == Ptypnnn || uniform <= Ptypnnn)
        unnn = gr->Give_variable(Innn,l)[v_.Number_variable()];
      else
        unnn = Give_at_Index(Innn,gr);

      h_mesh = gr->H_mesh() / Zweipotenz(gr->Max_level());
      h_d    = gr->Give_h(In,dd);
      h_d = h_d / h_mesh;

      return h_d*h_d * (un-2.0*unn+unnn)/2.0 
        + h_d * (2.0*unn - un / 2.0 - 3.0/2.0 * unnn)
        + unnn;
    }

    if(gr->Point_in_domain(Inn)) {
      Ptypnn  = gr->Give_type(Inn);
      if(Ptypnn >= interior) {
        un = Give_at_Index(In,gr);
        if(interior == Ptypnn || uniform <= Ptypnn)
          unn = gr->Give_variable(Inn,l)[v_.Number_variable()];
        else
          unn = Give_at_Index(Inn,gr);
        
        h_mesh = gr->H_mesh() / Zweipotenz(gr->Max_level());
        h_d    = gr->Give_h(In,dd);
        
        return ((un - unn) / h_mesh) * h_d + un;
      }
    }
    //    return Give_at_Index(In,gr);
    return v_.Give_Bo2p(it_b,gr,l);
  };
  */
  double Give_Bo2p(const P_Bo2p* it_b,const Grid* gr, int l) const {
    Index3D In, Inn;
    dir_3D dd;
    double un, unn, h_mesh, h_d;
    Pointtype Ptypnn;

    In = it_b->Ind();
    dd = it_b->d();
    if(dd==Wdir) { Inn = In.next_E(l); }
    if(dd==Edir) { Inn = In.next_W(l); }

    if(dd==Ndir) { Inn = In.next_S(l); }
    if(dd==Sdir) { Inn = In.next_N(l); }

    if(dd==Tdir) { Inn = In.next_D(l); }
    if(dd==Ddir) { Inn = In.next_T(l); }

    if(gr->Point_in_domain(Inn)) {
      Ptypnn  = gr->Give_type(Inn);
      if(Ptypnn >= interior) {
        un = Give_at_Index(In,gr);
        if(interior == Ptypnn || uniform <= Ptypnn)
          unn = gr->Give_variable(Inn,l)[v_.Number_variable()];
        else
          unn = Give_at_Index(Inn,gr);
        
        h_mesh = gr->H_mesh() / Zweipotenz(gr->Max_level());
        h_d    = gr->Give_h(In,dd);
        
        return ((un - unn) / h_mesh) * h_d + un;
      }
    }
    //    return Give_at_Index(In,gr);
    return v_.Give_Bo2p(it_b,gr,l);
  };
  double Give_cellpoi(const P_cellpoi* it_cf,const Grid* gr,
                      const BoCeData* bocedata) const {
    return v_.Give_cellpoi(it_cf,gr,bocedata);
  };
  
  // Size of necessary stencil
  int Sice_stencil() const { return 1; }
  
  // activity of points
  int Level() const {
    return v_.Give_grid()->Max_level();
  }
  Dominace_label Dominant_lev() const { return yes_dominant; }
  Dominace_label Dominant_poi() const { return yes_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 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_abstract;
  }
  // For abstract differential operator:
  int Give_number_var_of_abstract_op() const { return 0; }

  // 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 {
    return false; 
  };

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


// Extrapolation at the boundary
DExpr<DVar_Extra_Op>
Extrapolate_near_boundary(Variable a);


#endif

        

---