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

// ------------------------------------------------------------
// mg_op.h
//
// ------------------------------------------------------------

#ifndef MGOP_H_
#define MGOP_H_
                   
// 1.:  Restriction operator for Multigrid
// 2.:  Prolongation operator for Multigrid


// 1.:  Restriction operator for Multigrid



#define near_bo_restriction_coarse(insert_I) \
var = gr->Give_variable_slow( insert_I ,lpo); \
if(var!=NULL) sum = sum + var[v_.Number_variable()];


#define near_bo_restriction_fine(insert_I) \
sum = sum + gr->Give_variable( insert_I ,lpo)[v_.Number_variable()]; 


#define near_bo_restriction_fine_bo(insert_I,d,a,b) \
var = gr->Give_variable( insert_I , d ); \
if(var!=NULL) { \
  sum = sum + ((double)a + \
               gr->Give_h( insert_I , d ) / finest_meshsize \
               * (double)(b-a)) * \
    var[v_.Number_variable()]; \
}





// preparation


class DVar_Res_Op {
  Variable v_;
  double **Stencil;
  //
  double* weight27;

  Restriction_stencil_container* restriction_stencil;

 public:
  DVar_Res_Op(const Variable& v)
    : v_(v) {
    weight27 =  v.Give_grid()->Give_weight_27();
    restriction_stencil = v.Give_grid()->Give_restriction_stencil();
  }
  double Give_interior_here(const P_interior* it_i, const Grid* grid,
                            double h_mesh,
                            int lev, double_stencils_in) const {
    return 
      it_i->var_fine_M(grid,lev)[v_.Number_variable()] +
      0.5 * (it_i->var_fine_N(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_W(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_E(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_S(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_T(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_D(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_ND(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_WN(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_WT(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_ED(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_ST(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_ES(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_EST(grid,lev)[v_.Number_variable()] +
             it_i->var_fine_WND(grid,lev)[v_.Number_variable()]);
  };
  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_nearb(const P_nearb* it_n, const Grid* gr, double h_mesh, int l,
                    const Nearb_Ablage* nearb_ablage) const {
    double sum;
    double* var;
    bool from_finest_grid;
    Index3D next, Inext;
    int lpo; // l+1
    double finest_meshsize;
    bool cell_point_in_space;

    if(l == gr->Max_level()-1)
      from_finest_grid = true;
    else 
      from_finest_grid = false;

    finest_meshsize = gr->Give_finest_mesh_size();
  
    if(mg_coefficients) {
      Point_mg_coeff_nearb* l_nearb;  // Hilfspeicher 
      Point_mg_coeff_bo*    l_bo;     // Hilfspeicher 
    
      sum = 0.0;
      for(l_nearb=it_n->Give_first_mg_coeff_nearb();
          l_nearb!=NULL;l_nearb=l_nearb->Next_res()) {
        sum = sum + l_nearb->Give_var()[v_.Number_variable()] 
          * l_nearb->Give_weight();
      }
      if(v_.run_boundary() == 1) {
        for(l_bo=it_n->Give_first_mg_coeff_bo();
            l_bo!=NULL;l_bo=l_bo->Next_res()) {
          sum = sum + l_bo->Give_var()[v_.Number_variable()] 
            * l_bo->Give_weight();
        }
      }
      else if(v_.run_boundary() < 0) {
        for(l_bo=it_n->Give_first_mg_coeff_bo();
            l_bo!=NULL;l_bo=l_bo->Next_res()) {
          if(gr->Give_label_bo(l_bo->I(),l_bo->d(),-v_.run_boundary())) {
            sum = sum + l_bo->Give_var()[v_.Number_variable()] 
              * l_bo->Give_weight();
          }
        }
      }
    }
    else {
      Index3D I, next, next_cell;
      D3vector Co;
      BoCeData* bocedata;
      int i,j,k;
      double* UU;
      
      I = it_n->Ind();

      lpo = l+1;

      // the following code was mainly generated by 
      //    void Print_bos_for_restriction(...);
      // in loc_sten.cc
      // by hand : Mittelpunkt
      if(from_finest_grid) {
        // restrict from finest grid
        sum = 0.0;

        next = I.next_W(lpo).next_W(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Edir,0,1);
          //    };

        next = I.next_W(lpo).next_W(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Edir,0,1);
          //    };

        next = I.next_W(lpo).next_W(lpo).next_N(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Edir,0,1);
          //    };
        next = I.next_W(lpo).next_W(lpo).next_N(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Edir,0,1);
          //    };
        next = I.next_W(lpo).next_S(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Edir,0,1);
          near_bo_restriction_fine_bo(next,Ndir,0,1);
          //    };      
        next = I.next_W(lpo).next_S(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Edir,0,1);
          near_bo_restriction_fine_bo(next,Ndir,0,1);
          //    };
        next = I.next_W(lpo).next_D(lpo);  
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Edir,0,1);
          near_bo_restriction_fine_bo(next,Ndir,0,1);
          near_bo_restriction_fine_bo(next,Tdir,0,1);
          //    };
        next = I.next_W(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,0);
          near_bo_restriction_fine_bo(next,Edir,1,2);
          near_bo_restriction_fine_bo(next,Ndir,1,1);
          near_bo_restriction_fine_bo(next,Sdir,1,0);
          near_bo_restriction_fine_bo(next,Tdir,1,1);
          near_bo_restriction_fine_bo(next,Ddir,1,0);
          //    };
        next = I.next_W(lpo).next_T(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,0);
          near_bo_restriction_fine_bo(next,Edir,1,1);
          near_bo_restriction_fine_bo(next,Ndir,1,0);
          near_bo_restriction_fine_bo(next,Sdir,1,0);
          near_bo_restriction_fine_bo(next,Tdir,1,0);
          near_bo_restriction_fine_bo(next,Ddir,1,1);
          //    };
        next = I.next_W(lpo).next_T(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Ddir,0,1);
          //    };
        next = I.next_W(lpo).next_N(lpo).next_D(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Tdir,0,1);
          //    };
        next = I.next_W(lpo).next_N(lpo).next_D(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,0);
          near_bo_restriction_fine_bo(next,Edir,1,1);
          near_bo_restriction_fine_bo(next,Ndir,1,0);
          near_bo_restriction_fine_bo(next,Sdir,1,0);
          near_bo_restriction_fine_bo(next,Tdir,1,1);
          near_bo_restriction_fine_bo(next,Ddir,1,0);
          //    };
        next = I.next_W(lpo).next_N(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,0);
          near_bo_restriction_fine_bo(next,Edir,1,1);
          near_bo_restriction_fine_bo(next,Ndir,1,0);
          near_bo_restriction_fine_bo(next,Sdir,1,1);
          near_bo_restriction_fine_bo(next,Tdir,1,0);
          near_bo_restriction_fine_bo(next,Ddir,1,1);
          //    };
        next = I.next_W(lpo).next_N(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Sdir,0,1);
          near_bo_restriction_fine_bo(next,Ddir,0,1);
          //    };
        next = I.next_W(lpo).next_N(lpo).next_N(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Sdir,0,1);
          //    };
        next = I.next_W(lpo).next_N(lpo).next_N(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Sdir,0,1);
          //    };
        next = I.next_S(lpo).next_S(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Ndir,0,1);
          //    };
        next = I.next_S(lpo).next_S(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Ndir,0,1);
          //    };
        next = I.next_S(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Ndir,0,1);
          near_bo_restriction_fine_bo(next,Tdir,0,1);
          //    };
        next = I.next_S(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,0);
          near_bo_restriction_fine_bo(next,Edir,1,1);
          near_bo_restriction_fine_bo(next,Ndir,1,2);
          near_bo_restriction_fine_bo(next,Sdir,1,0);
          near_bo_restriction_fine_bo(next,Tdir,1,1);
          near_bo_restriction_fine_bo(next,Ddir,1,0);
          //    };
        next = I.next_S(lpo).next_T(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,0);
          near_bo_restriction_fine_bo(next,Edir,1,1);
          near_bo_restriction_fine_bo(next,Ndir,1,1);
          near_bo_restriction_fine_bo(next,Sdir,1,0);
          near_bo_restriction_fine_bo(next,Tdir,1,0);
          near_bo_restriction_fine_bo(next,Ddir,1,1);
          //    };
        next = I.next_S(lpo).next_T(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Ddir,0,1);
          //    };
        next = I.next_D(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Tdir,0,1);
          //    };
        next = I.next_D(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,0);
          near_bo_restriction_fine_bo(next,Edir,1,1);
          near_bo_restriction_fine_bo(next,Ndir,1,1);
          near_bo_restriction_fine_bo(next,Sdir,1,0);
          near_bo_restriction_fine_bo(next,Tdir,1,2);
          near_bo_restriction_fine_bo(next,Ddir,1,0);
          //    };
        // Mittelpunkt::
        next = I;
        if(gr->Grid_point_on_finest_level(next)) {
          sum = sum + 2.0 *
            gr->Give_variable( next ,lpo)[v_.Number_variable()];
          //      near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,2,1);
          near_bo_restriction_fine_bo(next,Edir,2,1);
          near_bo_restriction_fine_bo(next,Ndir,2,1);
          near_bo_restriction_fine_bo(next,Sdir,2,1);
          near_bo_restriction_fine_bo(next,Tdir,2,1);
          near_bo_restriction_fine_bo(next,Ddir,2,1);
          //    };
        next = I.next_T(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,1);
          near_bo_restriction_fine_bo(next,Edir,1,0);
          near_bo_restriction_fine_bo(next,Ndir,1,0);
          near_bo_restriction_fine_bo(next,Sdir,1,1);
          near_bo_restriction_fine_bo(next,Tdir,1,0);
          near_bo_restriction_fine_bo(next,Ddir,1,2);
          //    };
        next = I.next_T(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Ddir,0,1);
          //    };
        next = I.next_N(lpo).next_D(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Tdir,0,1);
          //    };
        next = I.next_N(lpo).next_D(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,1);
          near_bo_restriction_fine_bo(next,Edir,1,0);
          near_bo_restriction_fine_bo(next,Ndir,1,0);
          near_bo_restriction_fine_bo(next,Sdir,1,1);
          near_bo_restriction_fine_bo(next,Tdir,1,1);
          near_bo_restriction_fine_bo(next,Ddir,1,0);
          //    };
        next = I.next_N(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,1);
          near_bo_restriction_fine_bo(next,Edir,1,0);
          near_bo_restriction_fine_bo(next,Ndir,1,0);
          near_bo_restriction_fine_bo(next,Sdir,1,2);
          near_bo_restriction_fine_bo(next,Tdir,1,0);
          near_bo_restriction_fine_bo(next,Ddir,1,1);
          //    };
        next = I.next_N(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Sdir,0,1);
          near_bo_restriction_fine_bo(next,Ddir,0,1);
          //    };
        next = I.next_N(lpo).next_N(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Sdir,0,1);
          //    };
        next = I.next_N(lpo).next_N(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Sdir,0,1);
          //    };
        next = I.next_E(lpo).next_S(lpo).next_S(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Ndir,0,1);
          //    };
        next = I.next_E(lpo).next_S(lpo).next_S(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Ndir,0,1);
          //    };
        next = I.next_E(lpo).next_S(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Ndir,0,1);
          near_bo_restriction_fine_bo(next,Tdir,0,1);
          //    };
        next = I.next_E(lpo).next_S(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,1);
          near_bo_restriction_fine_bo(next,Edir,1,0);
          near_bo_restriction_fine_bo(next,Ndir,1,1);
          near_bo_restriction_fine_bo(next,Sdir,1,0);
          near_bo_restriction_fine_bo(next,Tdir,1,1);
          near_bo_restriction_fine_bo(next,Ddir,1,0);
          //    };
        next = I.next_E(lpo).next_S(lpo).next_T(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,1);
          near_bo_restriction_fine_bo(next,Edir,1,0);
          near_bo_restriction_fine_bo(next,Ndir,1,0);
          near_bo_restriction_fine_bo(next,Sdir,1,0);
          near_bo_restriction_fine_bo(next,Tdir,1,0);
          near_bo_restriction_fine_bo(next,Ddir,1,1);
          //    };
        next = I.next_E(lpo).next_S(lpo).next_T(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Ddir,0,1);
          //    };
        next = I.next_E(lpo).next_D(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Tdir,0,1);
          //    };
        next = I.next_E(lpo).next_D(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,1);
          near_bo_restriction_fine_bo(next,Edir,1,0);
          near_bo_restriction_fine_bo(next,Ndir,1,0);
          near_bo_restriction_fine_bo(next,Sdir,1,0);
          near_bo_restriction_fine_bo(next,Tdir,1,1);
          near_bo_restriction_fine_bo(next,Ddir,1,0);
          //    };
        next = I.next_E(lpo);
        if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine(next);
        }
          near_bo_restriction_fine_bo(next,Wdir,1,2);
          near_bo_restriction_fine_bo(next,Edir,1,0);
          near_bo_restriction_fine_bo(next,Ndir,1,0);
          near_bo_restriction_fine_bo(next,Sdir,1,1);
          near_bo_restriction_fine_bo(next,Tdir,1,0);
          near_bo_restriction_fine_bo(next,Ddir,1,1);
          //    };
        next = I.next_E(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Wdir,0,1);
          near_bo_restriction_fine_bo(next,Sdir,0,1);
          near_bo_restriction_fine_bo(next,Ddir,0,1);
          //    };
        next = I.next_E(lpo).next_N(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Wdir,0,1);
          near_bo_restriction_fine_bo(next,Sdir,0,1);
          //    };
        next = I.next_E(lpo).next_N(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Wdir,0,1);
          near_bo_restriction_fine_bo(next,Sdir,0,1);
          //    };
        next = I.next_E(lpo).next_E(lpo).next_S(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Wdir,0,1);
          //    };
        next = I.next_E(lpo).next_E(lpo).next_S(lpo).next_T(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Wdir,0,1);
          //    };
        next = I.next_E(lpo).next_E(lpo).next_D(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Wdir,0,1);
          //    };
        next = I.next_E(lpo).next_E(lpo);
        //      if(gr->Grid_point_on_finest_level(next)) {
          near_bo_restriction_fine_bo(next,Wdir,0,1);
          //    };

        // bocell point (Randzellfreiheitsgrad)
        for(i=0;i<8;++i) {    // go to middle point of neighbour coarse cells
          next_cell = I.next((dir_sons)i,lpo); 
          for(j=0;j<8;++j) {    // go to middle point of neighbour fine cells
            next = next_cell.next((dir_sons)j,l+2); 
            bocedata = gr->Give_Bo_cell(next);
            if(bocedata!=NULL) {
              //if(bocedata->Exists_bocellpoint() && false) {
              if(bocedata->Exists_bocellpoint()) {

                // cout << " fine_bo_cell " << endl;
                
                UU = gr->Give_vector_ablage();

                if(v_.run_boundary()==1) {
                  // Neumann boundary condition
                  cell_point_in_space = true;
                }
                else if(v_.run_boundary()==0) {
                  // Dirichlet boundary condition
                  cell_point_in_space = false;
                }
                else {
                  // mixed boundary condition
                  if(bocedata->Is_Dirichlet(-v_.run_boundary()))
                    cell_point_in_space = false;
                  else 
                    cell_point_in_space = true;
                }

                if(cell_point_in_space) { 
                  // Neumann boundary condition
                  for(k=0;k<8;++k) {
                    UU[k] =
                      restriction_stencil->Calc_trans_weight((dir_sons)i,
                                                             (dir_sons)j,
                                                             (dir_sons)k);
                  }
                }
                else {
                  // Dirichlet boundary condition
                  for(k=0;k<8;++k) {
                    Inext = next.neighbour((dir_sons)k);
                    if(gr->Grid_point_on_finest_level(Inext))
                      UU[k] =
                        restriction_stencil->Calc_trans_weight((dir_sons)i,
                                                               (dir_sons)j,
                                                               (dir_sons)k);
                    else
                      UU[k] = 0.0;
                  }
                }

                Co = bocedata->Coord_bocellpoint_normal();

                if(developer_version) {
                  if(Co.x<0.0 || Co.x>1.0 ||
                     Co.y<0.0 || Co.y>1.0 ||
                     Co.z<0.0 || Co.z>1.0) {
                    cout << " error in DVar_Res_Op::Give_nearb! " << endl;
                  }
                }

                sum +=
                  // (UU[0] + UU[1] + UU[2] + UU[3] +
                  // UU[4] + UU[5] + UU[6] + UU[7] ) * 0.125 * 
                  (UU[WSDd] * (1.0 - Co.x) * (1.0 - Co.y) * (1.0 - Co.z) +
                   UU[ESDd] *        Co.x  * (1.0 - Co.y) * (1.0 - Co.z) +
                   UU[WNDd] * (1.0 - Co.x) *         Co.y * (1.0 - Co.z) +
                   UU[ENDd] *        Co.x  *         Co.y * (1.0 - Co.z) +
                   UU[WSTd] * (1.0 - Co.x) * (1.0 - Co.y) *        Co.z  +
                   UU[ESTd] *        Co.x  * (1.0 - Co.y) *        Co.z  +
                   UU[WNTd] * (1.0 - Co.x) *         Co.y *        Co.z  +
                   UU[ENTd] *        Co.x  *         Co.y *        Co.z)  *
                  bocedata->vars[bocedata->Give_number_points()]
                  [v_.Number_variable()];


              }
            }
          }
        }
       


        return sum * 0.5;

      }
      else {
        // restrict from coarser grids
        sum = 0.0;

        near_bo_restriction_coarse(I.next_W(lpo));
        near_bo_restriction_coarse(I.next_WT(lpo));
        near_bo_restriction_coarse(I.next_WND(lpo));
        near_bo_restriction_coarse(I.next_WN(lpo));
        near_bo_restriction_coarse(I.next_S(lpo));
        near_bo_restriction_coarse(I.next_ST(lpo));
        near_bo_restriction_coarse(I.next_D(lpo));
        near_bo_restriction_coarse(I.next_T(lpo));
        near_bo_restriction_coarse(I.next_ND(lpo));
        near_bo_restriction_coarse(I.next_N(lpo));
        near_bo_restriction_coarse(I.next_ES(lpo));
        near_bo_restriction_coarse(I.next_EST(lpo));
        near_bo_restriction_coarse(I.next_ED(lpo));
        near_bo_restriction_coarse(I.next_E(lpo));
        
        sum = sum * 0.5;

        near_bo_restriction_coarse(I);

        return sum;
      }
    }
    return sum;
  };




  double Give_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
                      const BoCeData* bocedata) const {
    cout << " Error in DVar_Res_Op::Give_cellpoi! " << endl;
    cout << " This function should never be used! " << endl;
    return 0.0;
  };
  double Give_Bo2p(const P_Bo2p* it_b, const Grid* gr,int l)  const {
    cout << " Error in DVar_Res_Op::Give_Bo2p! " << endl;
    cout << " This function should never be used! " << endl;
    return 0.0;
  };
  // Size of necessary stencil
  int Sice_stencil() const { return 15; }

  // activity of points
  int Level() const { return (v_).Level()-1; }
  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 15; }

  // 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 {
    evpar->Variable_contained_in_expression_other(v_.Number_variable());
    evpar->Other_is_fine();
  }
  bool GS_type(int var_number_left, int dim) const {
    return false; 
  }; 
  // at the and of an iteration of an expression
  void clean() const { };
};


/* Unary operator */

DExpr<DVar_Res_Op>
Restriction_FE(Variable a);


// 2.:  Prolongation operator for Multigrid




class DVar_Prol_Op {
 public:
  DVar_Prol_Op(const Variable& v)
    : v_(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* grid,
                            double h_mesh,
                            int l, double_stencils_in) const {

    Index3D I;

    I = it_i->Ind();

    if(I.I_x().Tiefe()==l) {
      if(I.I_y().Tiefe()==l) {
        if(I.I_z().Tiefe()==l) {
          // middle point
          // edge between EST-WND
          return 0.5 *
            (it_i->var_coarse_EST(grid,l)[v_.Number_variable()] +
             it_i->var_coarse_WND(grid,l)[v_.Number_variable()]);
        }
        else {
          // TD-surface
          // edge between ES-WN
          return 0.5 *
            (it_i->var_coarse_ES(grid,l)[v_.Number_variable()] +
             it_i->var_coarse_WN(grid,l)[v_.Number_variable()]);
        }
      }
      else {
        if(I.I_z().Tiefe()==l) {
          // NS-surface
          // edge between WT-ED
          return 0.5 *
            (it_i->var_coarse_WT(grid,l)[v_.Number_variable()] +
             it_i->var_coarse_ED(grid,l)[v_.Number_variable()]);
        }
        else {
          // EW-edge
          return 0.5 *
            (it_i->var_coarse_E(grid,l)[v_.Number_variable()] +
             it_i->var_coarse_W(grid,l)[v_.Number_variable()]);
        }
      }
    }
    else {
      if(I.I_y().Tiefe()==l) {
        if(I.I_z().Tiefe()==l) {
          // EW-surface
          // edge between ST-ND
          return 0.5 *
            (it_i->var_coarse_ST(grid,l)[v_.Number_variable()] +
             it_i->var_coarse_ND(grid,l)[v_.Number_variable()]);
        }
        else {
          // NS-edge
          return 0.5 *
            (it_i->var_coarse_N(grid,l)[v_.Number_variable()] +
             it_i->var_coarse_S(grid,l)[v_.Number_variable()]);
        }
      }
      else {
        if(I.I_z().Tiefe()==l) {
          // TD-edge
          return 0.5 *
            (it_i->var_coarse_T(grid,l)[v_.Number_variable()] +
             it_i->var_coarse_D(grid,l)[v_.Number_variable()]);
        }
        else {
          // idem-point
          return it_i->var_coarse_M(grid,l)[v_.Number_variable()];
        }
      }
    }

  };

  double Give_nearb(const P_nearb* it_n, const Grid* grid, double h_mesh,
                    int l, const Nearb_Ablage* nearb_ablage) const {

    // return 0.0;
    return Give_at_Index(it_n->Ind(),grid,l);
  };

  double Give_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
                      const BoCeData* bocedata) const {
    double* UU;
    Index3D Inext;
    D3vector Co;
    bool cell_point_in_space;

    // return 0.0;

    UU = gr->Give_vector_ablage();

    if(v_.run_boundary()==1) {
      // Neumann boundary condition
      cell_point_in_space = true;
    }
    else if(v_.run_boundary()==0) {
      // Dirichlet boundary condition
      cell_point_in_space = false;
    }
    else {
      // mixed boundary condition
      if(bocedata->Is_Dirichlet(-v_.run_boundary()))
        cell_point_in_space = false;
      else 
        cell_point_in_space = true;
    }
    
    if(cell_point_in_space) {
      // Neumann boundary condition
      for(int i=0;i<8;++i) {
        Inext = it_cf->Ind().neighbour((dir_sons)i);
        UU[i] = Give_at_Index(Inext,gr,it_cf->Ind().Tiefe()-1);
      }
    }
    else {
      // Dirichlet boundary condition
      for(int i=0;i<8;++i) {
        Inext = it_cf->Ind().neighbour((dir_sons)i);
        if(gr->Grid_point_on_finest_level(Inext))
          UU[i] = Give_at_Index(Inext,gr,it_cf->Ind().Tiefe()-1);
        else 
          UU[i] = 0.0;
      }
    }

    Co = bocedata->Coord_bocellpoint_normal();
        
    return
      // (UU[0] + UU[1] + UU[2] + UU[3] +
      // UU[4] + UU[5] + UU[6] + UU[7] ) * 0.125 * 
      (UU[WSDd] * (1.0 - Co.x) * (1.0 - Co.y) * (1.0 - Co.z) +
       UU[ESDd] *        Co.x  * (1.0 - Co.y) * (1.0 - Co.z) +
       UU[WNDd] * (1.0 - Co.x) *         Co.y * (1.0 - Co.z) +
       UU[ENDd] *        Co.x  *         Co.y * (1.0 - Co.z) +
       UU[WSTd] * (1.0 - Co.x) * (1.0 - Co.y) *        Co.z  +
       UU[ESTd] *        Co.x  * (1.0 - Co.y) *        Co.z  +
       UU[WNTd] * (1.0 - Co.x) *         Co.y *        Co.z  +
       UU[ENTd] *        Co.x  *         Co.y *        Co.z);
  };

  double Give_Bo2p(const P_Bo2p* it_b, const Grid* gr,int l)  const {
    double a,b;
    
    a = Give_at_Index(it_b->Ind(),gr,l);
    b = Give_at_Index(it_b->Ind().next(it_b->d(),l),gr,l);
    
    return 
      (a + 
       gr->Give_h(it_b->Ind(),it_b->d()) / gr->Give_finest_mesh_size()
       * (b-a));
  };

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

  // activity of points
  int Level() const { return (v_).Level()+1; }
  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 2; }

  // 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 {
    evpar->Variable_contained_in_expression_other(v_.Number_variable());
    evpar->Other_is_coarse();
  }
  bool GS_type(int var_number_left, int dim) const {
    return false; 
  };

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

 private:
  Variable v_;
  double **Stencil;

  double Give_at_Index(const Index3D I, const Grid* grid, int l) const;
};


/* Unary operator */

DExpr<DVar_Prol_Op>
Prolongation_FE(Variable a);


#endif

           

---