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

// ------------------------------------------------------------
// res_op.h
//
// ------------------------------------------------------------

#ifndef RESOP_H_
#define RESOP_H_
         
#include "grid/SSten.h"

//   Operators with stencils near the boundary

class Res_stencil_boundary;


// a) Abstract operator with variable
//-------------------------------------------------
template<class A>
class DResDiff_Bo {
  A a_;

  //number of storage for boundary stencils
  int num_stencils;

  // information about the variable put in the abstract diff. operator
  int number_variable;  // number of the variable
  int r_bou;            // boundary conditions of the variable 
                        // (activity of the boundary points)
  bool* label_on_cell;

  // Gitter
  Grid *grid;

  // needed for the calculation of local stiffness matrix
  double** Stencils_fine;  // eight stencils for finer cells
  double* u_ablage; // Ablage
  bool question_calc_stencil;
  bool *corner_in_domain;

 public:
  // Constructor for concrete operators
  DResDiff_Bo(const A& a, int num_s, Grid* gr, bool isRestricted);
  void Iterate_Calc_stencil(int r_bo) const;

  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 {
    static int i;
    static Celltype cetyp;
    static double sum;
    dir_sons dir_v;

    // Part 1: interior cells and boundary cells on the finest grid
    sum = a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage);

    // Part 2: boundary cells on coarser grid
    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==bo_cell) {
        sum = sum +
          Evaluate_lsm(nearb_ablage->Give_u_Recell((dir_sons)i),
                       number_variable, dir_v,
                       gr->Give_bo_stencil(it_n->Ind().next((dir_sons)i,l+1),
                                           num_stencils));
      }
    }    
    return sum;
  };
  double Give_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
                      const BoCeData* bocedata) const {
    return a_.Give_cellpoi(it_cf,gr,bocedata);
  };
  double Give_Bo2p(const P_Bo2p* it_b, const Grid* gr,int l)  const {
    return a_.Give_Bo2p(it_b,gr,l);
  };
  // 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 not_dominant; }
  Dominace_label Dominant_poi() const { return not_dominant; }
  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(); }

  // 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 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) Diagonal abstract operator
//-------------------------------------------------
template<class A>
class DResDiagDiff_Bo {
  A a_;

  //number of storage for boundary stencils
  int num_stencils;

  // Gitter
  Grid *grid;
  
  // Corresponding object for the stancils
  Res_stencil_boundary *ResS;

 public:
  // Constructor for concrete operators
  DResDiagDiff_Bo(const A& a, int num_s, Grid* gr, Res_stencil_boundary *resS) 
    : a_(a), num_stencils(num_s), grid(gr), ResS(resS) { }
  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 {
    static int i;
    static Celltype cetyp;
    static double sum;
    dir_sons dir_v;
    sum = a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage);

    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==bo_cell) {
        sum = sum + gr->Give_bo_stencil(it_n->Ind().next((dir_sons)i,l+1),
                                        num_stencils)[9*dir_v];
      }
    }

    return sum;
  };
  double Give_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
                      const BoCeData* bocedata) const {
    return a_.Give_cellpoi(it_cf,gr,bocedata);
  };
  double Give_Bo2p(const P_Bo2p* it_b, const Grid* gr,int l)  const {
    return a_.Give_Bo2p(it_b,gr,l);
  };
  // Size of necessary stencil
  int Sice_stencil() const { return a_.Sice_stencil(); }

  // 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 a_.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(Grid* gr, int r_bo) const ;    // ada: implementation further down in order to make gcc 3.4.3 happy

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


class Array_Res_stencil_boundary;
class Res_stencil_boundary {
  friend class Array_Res_stencil_boundary;
 private:
  // is stencil already used?
  // this has to be changed later!
  bool used;

  // Gitter
  Grid *grid;

  //number of storage for boundary stencils
  int num_stencils;

  Res_stencil_boundary() : grid(NULL) {};
  void Initialize(Grid* gr)  {
    grid = gr;
    used = false;
    num_stencils = gr->Give_number_bo_stencil();

    if(gr->Is_storage_initialized()==true) {
      cout << " \n Definition of Res_stencil_boundary not possible! ";
      cout << " \n Define Res_stencil_boundary before ";
      cout << " initialization of the grid! " << endl;
    }
  }
 public:
  Res_stencil_boundary(Grid* gr) {
    Initialize(gr);
  }
  
  bool Give_used() { return used;  }
  void clear()     { used = false; }

  void Calc_stencil(int r_bo);

  template<class A>
  DExpr<DResDiff_Bo<DExpr<A> > > 
  operator() (const DExpr<A>& a_) {
    Differential_op_typ a_abstract;
    a_abstract = a_.Abstract_differential_operator();
    if(a_abstract==not_abstract) {
      cout << " This is not an abstract differential operator! Error!" << endl;
    }
    if(a_abstract==abstract_diag)
      cout<<"\n Please use operator [] instead of operator ()! Error!" << endl;

    typedef DResDiff_Bo<DExpr<A> > ExprT;
    if(used==false) { // restrict stencil
      used = true;
      return DExpr<ExprT>(ExprT(a_,num_stencils,grid,true));
    }
    else              // restrict stencil not necessary
      return DExpr<ExprT>(ExprT(a_,num_stencils,grid,false));
  }


  template<class A>
  DExpr<DResDiagDiff_Bo<DExpr<A> > > 
  operator[] (const DExpr<A>& a_) {
    Differential_op_typ a_abstract;
    a_abstract = a_.Abstract_differential_operator();
    if(a_abstract==not_abstract) {
      cout << " This is not an abstract differential operator! Error!" << endl;
    }
    if(a_abstract==abstract_with_var)
      cout<<"\n Please use operator () instead of operator []! Error!" << endl;
    typedef DResDiagDiff_Bo<DExpr<A> > ExprT;
    return DExpr<ExprT>(ExprT(a_,num_stencils,grid,this));
  }
};


// Implementation of some member functions

template<class A>
void DResDiagDiff_Bo<A>::Put_grid_rbo(Grid* gr, int r_bo) const {
  a_.Put_grid_rbo(gr,r_bo);
  if(ResS->Give_used()==false) {
    cout << "\n An abstract diagonal differential operator without";
    cout << "\n a variable cannot be restricted! Error! " << endl;
  }
}



template<class A>
DResDiff_Bo<A>::DResDiff_Bo(const A& a, int num_s, Grid* gr, bool isRestricted)
  : a_(a), num_stencils(num_s), grid(gr) {

  number_variable = a_.Give_number_var_of_abstract_op();
  r_bou           = a_.run_boundary();

  if(isRestricted) {
    Stencils_fine = gr->Give_Stencils_fine();
    u_ablage = gr->Give_vector_ablage();
    label_on_cell = gr->Give_label_on_cell_ablage();
    corner_in_domain = gr->Give_eight_bools();

    question_calc_stencil = true;
  }
  else {
    question_calc_stencil = false;
  }
}

template<class A>
void DResDiff_Bo<A>::Put_grid_rbo(Grid* gr, int r_bov) const {
  a_.Put_grid_rbo(gr,r_bov);

  if(question_calc_stencil == true) {
    DResDiff_Bo<A>::Iterate_Calc_stencil(r_bov);
  }
};


template<class A>
void DResDiff_Bo<A>::Iterate_Calc_stencil(int r_bov) const {
  int i, j, jj, level, lev;
  Celltype ce_typ;
  Index3D I_son, I_ecke;
  BoCell* b_cell;
  int num_v, num_u, dir_u;
  double* sp;
  double* Stencil_coarse;
  double mesh_size;
  double finest_mesh_size;
  double a,b, ss, hh;
  double weight_bocellpoint[8];
  dir_3D d;
  Index3D my_lev_index;
  D3vector Co;
  bool u_in_space, uu_in_space, v_in_space;
  bool cell_point_in_space;
  Pointtype type;
  int *coarse_grid_functions;

  Index3D I;

  coarse_grid_functions = grid->Give_coarse_grid_functions();

  // for iteration
  int i_iter, hashtable0_leng;
  Point_hashtable0* point0;
  Point_hashtable0** hashtable0_start;  //Zeiger auf Tabelle
  hashtable0_leng  = grid->Give_hashtable0_leng();
  hashtable0_start = grid->Give_hashtable0_start();


  // Remark: By recursion cell type of I is: bo_cell

  // I. First step: cells on one level under finest level
  finest_mesh_size = grid->Give_finest_mesh_size();
  level=grid->Max_level();
  Iterate_hash0 {
   // a) search for boundary cells    
   I = point0->Give_Index();
   if(I.Cell_index()) {
    lev = I.Tiefe();
    if(point0->Give_cell_typ()==bo_cell && lev == level) {
      // b) Get the finer stencils
      // mesh_size = grid->H_mesh() / Zweipotenz(I.Tiefe()+1);
      mesh_size = grid->H_mesh() / Zweipotenz(I.Tiefe());
      for(i=0;i<8;++i) {
        // iterate 8 finer cells
        I_son = I.son((dir_sons)i);
        ce_typ = grid->Give_cell_typ(I_son);
        if(ce_typ==int_cell) {
          // if int_cell
          for(j=0;j<64;++j)
            Stencils_fine[i][j]
              = a_.Give_interior_sten_element(I_son,mesh_size,j,
                                              grid,I.Tiefe());
        }
        else if(ce_typ==fine_bo_cell) {
          // if fine_bo_cell
          // first: 1): set zero
          for(j=0;j<64;++j) Stencils_fine[i][j] = 0.0;

          // second: 2): set bocell
          b_cell = grid->Give_Bo_cell(I_son);
          
          // third: 3): calc weights for bocellpoint
          if(b_cell->Exists_bocellpoint()) {
            Co = b_cell->Coord_bocellpoint_normal();
            weight_bocellpoint[WSDd] = 
              (1.0 - Co.x) * (1.0 - Co.y) * (1.0 - Co.z);
            weight_bocellpoint[ESDd] =
              Co.x  * (1.0 - Co.y) * (1.0 - Co.z);
            weight_bocellpoint[WNDd] = 
              (1.0 - Co.x) *         Co.y * (1.0 - Co.z);
            weight_bocellpoint[ENDd] =
              Co.x  *         Co.y * (1.0 - Co.z);
            weight_bocellpoint[WSTd] = 
              (1.0 - Co.x) * (1.0 - Co.y) *        Co.z; 
            weight_bocellpoint[ESTd] =
              Co.x  * (1.0 - Co.y) *        Co.z; 
            weight_bocellpoint[WNTd] =
              (1.0 - Co.x) *         Co.y *        Co.z; 
            weight_bocellpoint[ENTd] =
              Co.x  *         Co.y *        Co.z;
          }

          // forth: 4): calculate local stiffness matrix
          for(dir_u=0;dir_u<8;++dir_u) { // iter. on u
            if(r_bou==1) {
              // case pure Neumann for u
              u_in_space = true;
            }
            else if(r_bou==0) {
              // case pure Dirichlet for u
              u_in_space = false;
            }
            else {
              // case mixed for u
              I_ecke = I_son.neighbour((dir_sons)dir_u);
              if(grid->Give_global_type(I_ecke)==multigrid) {
                if(grid->Give_label_bo_D(I_ecke,-r_bou,level))
                  u_in_space = false;
                else
                  u_in_space = true;
              }
            }

            for(num_u=0;num_u<b_cell->Give_number_points();++num_u) {
              j = b_cell->corner(num_u);
              if(b_cell->edge_point(num_u) == corner_poi_typ) {
                // num_u is corner point
                if(j==dir_u) u_ablage[num_u] = 1.0;
                else         u_ablage[num_u] = 0.0;
                if(r_bou<=0) {
                  // case pure Dirichlet and mixed for u
                  if(j==dir_u) u_in_space = true;
                }
              }
              else {
                // num_u is edge point          
                if(r_bou==1) {
                  // case pure Neumann for u
                  uu_in_space = true;
                  d = b_cell->edge_dir(num_u);
                  
                }
                else if(r_bou<0) {
                  // case mixed for u
                  d = b_cell->edge_dir(num_u);
                  uu_in_space = 
                    grid->Give_label_bo(I_son.neighbour((dir_sons)j), // "I"
                                        d,                            // "d"
                                        -r_bou);
                }
                else {
                  // case pure Dirichlet for u
                  uu_in_space = false;
                }
                if(uu_in_space) {
                  if(j==dir_u) a = 1.0;
                  else         a = 0.0;
                  if(Next_corner((dir_sons)j,d)==dir_u) b = 1.0;
                  else                                  b = 0.0;
                  u_ablage[num_u] =
                    (a + grid->Give_h(I_son.neighbour((dir_sons)j),d)
                     / finest_mesh_size
                     * (b-a));
                }
                else {
                  u_ablage[num_u] = 0.0;
                }
              }
            }
            if(b_cell->Exists_bocellpoint()) {
              if(r_bou==1) {
                // case pure Neumann
                u_in_space = true;
              }
              else if(r_bou<0) {
                // case mixed for u
                if(b_cell->Is_Dirichlet(-r_bou)==false)
                  u_in_space = true;
              }
              // num_u is bocellpoint
              num_u = b_cell->Give_number_points();
              u_ablage[num_u] = weight_bocellpoint[dir_u];
            }

            if(u_in_space) { // go on only when u is in space
              // v: cases of boundary conditions for v!
              // corners of boundary cell
              for(num_v=0;num_v<b_cell->Give_number_points();++num_v) {
                // iter. on v
                j = b_cell->corner(num_v);
                if(b_cell->edge_point(num_v) == corner_poi_typ) {
                  // num_v is corner point
                  Stencils_fine[i][j+8*dir_u] 
                    += a_.Give_boundary_sten_element(grid,b_cell,
                                                    u_ablage,num_v);
                }
                else {
                  // num_v is edge point
                  if(r_bov==1) {
                    // case pure Neumann for v
                    d = b_cell->edge_dir(num_v);
                    v_in_space = true;
                  }
                  else if(r_bov<0) {
                    // case mixed for v
                    d = b_cell->edge_dir(num_v);
                    v_in_space = 
                      grid->Give_label_bo(I_son.neighbour((dir_sons)j), // "I"
                                          d,                            // "d"
                                          -r_bov);
                  }
                  else {
                    // case Dirichlet for v
                    v_in_space = false;
                  }
                  if(v_in_space) {
                    ss = a_.Give_boundary_sten_element(grid,b_cell,
                                                       u_ablage,num_v);
                    d = b_cell->edge_dir(num_v);
                    j = b_cell->corner(num_v);
                    
                    hh = grid->Give_h(I_son.neighbour((dir_sons)j),d)
                      / finest_mesh_size;
                    // corner at num_v 
                    Stencils_fine[i][j+8*dir_u] 
                      += ss * (1.0 + hh * (0.0-1.0));
                    // opposite corner of num_v 
                    jj = Next_corner((dir_sons)j,d);
                    Stencils_fine[i][jj+8*dir_u]
                      += ss * (0.0 + hh * (1.0-0.0));
                  }
                }
              }

              if(b_cell->Exists_bocellpoint()) {
                // num_v is bocellpoint
                num_v = b_cell->Give_number_points();
                ss = a_.Give_boundary_sten_element(grid,b_cell,
                                                   u_ablage,num_v);
                if(r_bov==1) {
                  // case pure Neumann
                  cell_point_in_space = true;
                }
                else if(r_bov<0) {
                  // case mixed for u
                  if(b_cell->Is_Dirichlet(-r_bov)==false)
                    cell_point_in_space = true;
                  else 
                    cell_point_in_space = false;
                }
                else
                  // case pure Dirichlet
                  cell_point_in_space = false;
                if(cell_point_in_space) {
                  // case pure Neumann  or mixed for v
                  for(j=0;j<8;++j)
                    Stencils_fine[i][j+8*dir_u] += ss * weight_bocellpoint[j];
                }
                else {
                  // case pure Dirichlet for v
                  for(num_v=0;num_v<b_cell->Give_number_points();++num_v) {
                    if(b_cell->edge_point(num_v) == corner_poi_typ) {
                      j = b_cell->corner(num_v);
                      Stencils_fine[i][j+8*dir_u] += ss * weight_bocellpoint[j];
                    }
                  }
                }
              }
            }
          }
        }
        else {
          for(j=0;j<64;++j) Stencils_fine[i][j] = 0.0;
        }
      }


      // c) Coarse the finer stencils
      Stencil_coarse = grid->Give_bo_stencil(I,num_stencils);
      Restrict_local_stiffness_matrix(Stencil_coarse,
                                      Stencils_fine,coarse_grid_functions);

      // Dirichlet for v
      if(r_bov==0) {
        for(i=0;i<8;++i) 
          if(grid->Give_global_type(I.neighbour((dir_sons)i))==multigrid) {
            for(j=0;j<8;++j) {
              Stencil_coarse[i+8*j] = 0.0;
            }
          }
      }
      // Dirichlet for u
      if(r_bou==0) {
        for(i=0;i<8;++i) 
          if(grid->Give_global_type(I.neighbour((dir_sons)i))==multigrid) {
            for(j=0;j<8;++j) {
              Stencil_coarse[j+8*i] = 0.0;
            }
          }
      }
      
      // mixed for v
      if(r_bov<0) {
        for(i=0;i<8;++i) {
          type = grid->Give_global_type(I.neighbour((dir_sons)i));
          if(type==multigrid) {
            if(grid->Give_label_bo_D(I.neighbour((dir_sons)i),-r_bov,level-1)) {
              for(j=0;j<8;++j) {
                Stencil_coarse[i+8*j] = 0.0;
              }
            }
          }
          if(type==exterior) {
            for(j=0;j<8;++j) {
              Stencil_coarse[i+8*j] = 0.0;
            }
          }
        }
      }
      // mixed for u
      if(r_bou<0) {
        for(i=0;i<8;++i) {
          type = grid->Give_global_type(I.neighbour((dir_sons)i));
          if(type==multigrid) {
            if(grid->Give_label_bo_D(I.neighbour((dir_sons)i),-r_bou,level-1)) {
              for(j=0;j<8;++j) {
                Stencil_coarse[j+8*i] = 0.0;
              }
            }
          }
          if(type==exterior) {
            for(j=0;j<8;++j) {
              Stencil_coarse[j+8*i] = 0.0;
            }
          }
        }
      }
    }
   }
  }
  // d) Parallel:
  if(parallel_version)
    grid->Update_boundary_stencil(num_stencils,grid->Max_level());

  // II. Second step: all coarser level
  for(level=grid->Max_level()-1;level>grid->Min_level();--level) {
    my_lev_index = grid->Give_my_level_index(level);
    Iterate_hash0 {
     // a) search for boundary cells
     I = point0->Give_Index();
     if(I.Cell_index()) {
      lev = I.Tiefe();
      if(point0->Give_cell_typ()==bo_cell && lev == level &&
         I < my_lev_index) {

        // b) Get the finer stencils
        // mesh_size = grid->H_mesh() / Zweipotenz(I.Tiefe()+1);
        mesh_size = grid->H_mesh() / Zweipotenz(I.Tiefe());
        for(i=0;i<8;++i) {
          I_son = I.son((dir_sons)i);
          ce_typ = grid->Give_cell_typ(I_son);

          if(ce_typ==int_cell) {
            for(j=0;j<64;++j)
              Stencils_fine[i][j]
                = a_.Give_interior_sten_element(I_son,mesh_size,j,
                                                grid,I.Tiefe());
          }
          else if(ce_typ==bo_cell) {
            sp = grid->Give_bo_stencil(I_son,num_stencils);
            for(j=0;j<64;++j) Stencils_fine[i][j] = sp[j];
          }
          else {
            for(j=0;j<64;++j) Stencils_fine[i][j] = 0.0;
          }
        }
  
        // c) Coarse the finer stencils
        Stencil_coarse = grid->Give_bo_stencil(I,num_stencils);
        Restrict_local_stiffness_matrix(Stencil_coarse,
                                        Stencils_fine,coarse_grid_functions);

        // Dirichlet for v
        if(r_bov==0) {
          for(i=0;i<8;++i) 
            if(grid->Give_global_type(I.neighbour((dir_sons)i))==multigrid) {
              for(j=0;j<8;++j) {
                Stencil_coarse[i+8*j] = 0.0;
              }
            }
        }
        // Dirichlet for u
        if(r_bou==0) {
          for(i=0;i<8;++i) 
            if(grid->Give_global_type(I.neighbour((dir_sons)i))==multigrid) {
              for(j=0;j<8;++j) {
                Stencil_coarse[j+8*i] = 0.0;
              }
            }
        }

        // mixed for v
        if(r_bov<0) {
          for(i=0;i<8;++i) {
            type = grid->Give_global_type(I.neighbour((dir_sons)i));
            if(type==multigrid) {
              if(grid->Give_label_bo_D(I.neighbour((dir_sons)i),-r_bov,level-1)) {
                for(j=0;j<8;++j) {
                  Stencil_coarse[i+8*j] = 0.0;
                }
              }
            }
            if(type==exterior) {
              for(j=0;j<8;++j) {
                Stencil_coarse[i+8*j] = 0.0;
              }
            }
          }
        }
        // mixed for u
        if(r_bou<0) {
          for(i=0;i<8;++i) {
            type = grid->Give_global_type(I.neighbour((dir_sons)i));
            if(type==multigrid) {
              if(grid->Give_label_bo_D(I.neighbour((dir_sons)i),-r_bou,level-1)) {
                for(j=0;j<8;++j) {
                  Stencil_coarse[j+8*i] = 0.0;
                }
              }
            }
            if(type==exterior) {
              for(j=0;j<8;++j) {
                Stencil_coarse[j+8*i] = 0.0;
              }
            }
          }
        }
        
      }
     }
    }
    // d) Parallel:
    if(parallel_version) {
      grid->Update_boundary_stencil(num_stencils,level);
    }
  }
}




#endif

           

---