00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #ifndef RESOPV_H_
00023 #define RESOPV_H_
00024
00025
00027
00029
00030 class Res_stencil;
00031
00032
00033
00034 template<class A>
00035 class DResDiff {
00036 A a_;
00037
00038
00039 int num_stencils;
00040
00041
00042 int number_variable;
00043 int r_bou;
00044
00045 bool* label_on_cell;
00046
00047
00048 Grid *grid;
00049
00050
00051 double** Stencils_fine;
00052 double* u_ablage;
00053 bool question_calc_stencil;
00054 bool *corner_in_domain;
00055
00056 public:
00057
00058 DResDiff(const A& a, int num_s, Grid* gr, bool isRestricted);
00059
00060 void Iterate_Calc_stencil(int r_bov) const;
00061
00062 double Give_interior(const P_interior* it_i, const Grid* grid, double h_mesh,
00063 int lev, double_stencils_in) const {
00064 return a_.Give_interior(it_i,grid,h_mesh,lev, double_stencils_out);
00065 }
00066 double Give_interior_coarse(const P_interior* it_i, const Grid* gr,
00067 double h_mesh,
00068 int l, double_stencils_in) const {
00069 double sum;
00070 dir_sons dir_v;
00071 Nearb_Ablage* nearb_ablage;
00072 nearb_ablage = gr->Give_Nearb_Ablage();
00073 nearb_ablage->Initialize_27(gr,it_i->Ind(),l);
00074 sum = 0.0;
00075 for(int i=0;i<8;++i) {
00076 dir_v = opposite3D((dir_sons)i);
00077
00078 sum = sum +
00079 Evaluate_lsm(nearb_ablage->Give_u_Recell((dir_sons)i),
00080 number_variable, dir_v,
00081 gr->Give_stencil(it_i->Ind().next((dir_sons)i,l+1),
00082 num_stencils));
00083 }
00084 return sum;
00085 }
00086 double Give_nearb(const P_nearb* it_n, const Grid* gr, double h_mesh, int l,
00087 const Nearb_Ablage* nearb_ablage) const {
00088 static int i;
00089 static Celltype cetyp;
00090 static double sum;
00091 dir_sons dir_v;
00092
00093
00094 if(l==gr->Max_level()) {
00095 sum = a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage);
00096 }
00097 else {
00098 Nearb_Ablage* nearb_ablage2;
00099 nearb_ablage2 = gr->Give_Nearb_Ablage();
00100 nearb_ablage2->Initialize_27(gr,it_n->Ind(),l);
00101
00102
00103
00104 sum = 0.0;
00105
00106 for(i=0;i<8;++i) {
00107 dir_v = opposite3D((dir_sons)i);
00108
00109 cetyp = nearb_ablage->Give_Celltyp((dir_sons)i);
00110 if(cetyp==bo_cell || cetyp==int_cell) {
00111 sum = sum +
00112 Evaluate_lsm(nearb_ablage->Give_u_Recell((dir_sons)i),
00113 number_variable, dir_v,
00114 gr->Give_stencil(it_n->Ind().next((dir_sons)i,l+1),
00115 num_stencils));
00116 }
00117 }
00118 }
00119 return sum;
00120 };
00121 double Give_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
00122 const BoCeData* bocedata) const {
00123 return a_.Give_cellpoi(it_cf,gr,bocedata);
00124 };
00125 double Give_Bo2p(const P_Bo2p* it_b, const Grid* gr,int l) const {
00126 return a_.Give_Bo2p(it_b,gr,l);
00127 };
00128
00129 int Sice_stencil() const { return a_.Sice_stencil(); }
00130
00131
00132 int Level() const { return (a_).Level(); }
00133 Dominace_label Dominant_lev() const { return not_dominant; }
00134 Dominace_label Dominant_poi() const { return not_dominant; }
00135 bool run_interior() const { return (a_).run_interior(); }
00136 bool run_nearb() const { return (a_).run_nearb(); }
00137 int run_boundary() const { return (a_).run_boundary(); }
00138
00139
00140 int ops_interior() const { return a_.ops_interior(); }
00141
00142
00143 void Active_Sim_Level( int lev) const { };
00144 void Active_Sim_interior(bool run) const { };
00145 void Active_Sim_nearb(bool run) const { };
00146 void Active_Sim_boundary( int run) const { };
00147 void Active_Sim_update(Evaluation_Parallelization_object* evpar,
00148 int level, type_of_update typ) const { };
00149
00150
00151
00152 void Put_grid_rbo(Grid* gr, int r_bo) const;
00153
00154
00155
00156 Differential_op_typ Abstract_differential_operator() const
00157 { return not_abstract; }
00158
00159 void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar)
00160 const {
00161 a_.Add_variables_for_parallel(evpar);
00162 };
00163 bool GS_type(int var_number_left, int dim) const {
00164 return a_.GS_type(var_number_left,dim);
00165 };
00166
00167 void clean() const { a_.clean(); };
00168 };
00169
00170
00171
00172 template<class A>
00173 class DResDiagDiff {
00174 A a_;
00175
00176
00177 int num_stencils;
00178
00179
00180 Grid *grid;
00181
00182
00183 Res_stencil *ResS;
00184
00185 public:
00186
00187 DResDiagDiff(const A& a, int num_s, Grid* gr, Res_stencil *resS)
00188 : a_(a), num_stencils(num_s), grid(gr), ResS(resS) { }
00189
00190 double Give_interior(const P_interior* it_i, const Grid* gr, double h_mesh,
00191 int l, double_stencils_in) const {
00192 return a_.Give_interior(it_i,gr,h_mesh,l, double_stencils_out);
00193 }
00194 double Give_interior_coarse(const P_interior* it_i, const Grid* gr,
00195 double h_mesh,
00196 int l, double_stencils_in) const {
00197 double sum;
00198 int i;
00199 dir_sons dir_v;
00200 sum = 0.0;
00201 for(i=0;i<8;++i) {
00202 dir_v = opposite3D((dir_sons)i);
00203
00204 sum = sum + gr->Give_stencil(it_i->Ind().next((dir_sons)i,l+1),
00205 num_stencils)[9*dir_v];
00206 }
00207 return sum;
00208 };
00209 double Give_nearb(const P_nearb* it_n, const Grid* gr, double h_mesh, int l,
00210 const Nearb_Ablage* nearb_ablage) const {
00211 static int i;
00212 static Celltype cetyp;
00213 static double sum;
00214 dir_sons dir_v;
00215 if(l==gr->Max_level())
00216 sum = a_.Give_nearb(it_n,gr,h_mesh,l,nearb_ablage);
00217 else {
00218 sum = 0.0;
00219 for(i=0;i<8;++i) {
00220 dir_v = opposite3D((dir_sons)i);
00221
00222 cetyp = nearb_ablage->Give_Celltyp((dir_sons)i);
00223 if(cetyp==bo_cell || cetyp==int_cell) {
00224 sum = sum + gr->Give_stencil(it_n->Ind().next((dir_sons)i,l+1),
00225 num_stencils)[9*dir_v];
00226 }
00227 }
00228 }
00229 return sum;
00230 };
00231 double Give_cellpoi(const P_cellpoi* it_cf, const Grid* gr,
00232 const BoCeData* bocedata) const {
00233 return a_.Give_cellpoi(it_cf,gr,bocedata);
00234 };
00235 double Give_Bo2p(const P_Bo2p* it_b, const Grid* gr,int l) const {
00236 return a_.Give_Bo2p(it_b,gr,l);
00237 };
00238
00239 int Sice_stencil() const { return a_.Sice_stencil(); }
00240
00241
00242 int Level() const { return 0; }
00243 Dominace_label Dominant_lev() const { return anti_dominant; }
00244 Dominace_label Dominant_poi() const { return anti_dominant; }
00245 bool run_interior() const { return true; }
00246 bool run_nearb() const { return true; }
00247 int run_boundary() const { return true; }
00248
00249
00250 int ops_interior() const { return a_.ops_interior(); }
00251
00252
00253 void Active_Sim_Level( int lev) const { };
00254 void Active_Sim_interior(bool run) const { };
00255 void Active_Sim_nearb(bool run) const { };
00256 void Active_Sim_boundary( int run) const { };
00257 void Active_Sim_update(Evaluation_Parallelization_object* evpar,
00258 int level, type_of_update typ) const { };
00259
00260
00261 void Put_grid_rbo(Grid* gr, int r_bo) const ;
00262
00263
00264
00265
00266 Differential_op_typ Abstract_differential_operator() const
00267 { return not_abstract; }
00268
00269 void Add_variables_for_parallel(Evaluation_Parallelization_object *evpar)
00270 const {
00271 a_.Add_variables_for_parallel(evpar);
00272 };
00273 bool GS_type(int var_number_left, int dim) const {
00274 return a_.GS_type(var_number_left,dim);
00275 };
00276
00277 void clean() const { a_.clean(); };
00278 };
00279
00280
00281
00282
00283 class Res_stencil {
00284 private:
00285
00286
00287 bool used;
00288
00289
00290 Grid *grid;
00291
00292
00293 int num_stencils;
00294
00295 public:
00296 Res_stencil(Grid* gr) : grid(gr) {
00297 used = false;
00298 num_stencils = gr->Give_number_stencil();
00299
00300 if(gr->Is_storage_initialized()==true) {
00301 cout << " \n Definition of Res_stencil not possible! ";
00302 cout << " \n Define Res_stencil before ";
00303 cout << " initialization of the grid! " << endl;
00304 }
00305 }
00306
00307 bool Give_used() { return used; }
00308 void clear() { used = false; }
00309
00310 void Calc_stencil(int r_bo);
00311
00312 template<class A>
00313 DExpr<DResDiff<DExpr<A> > >
00314 operator() (const DExpr<A>& a_) {
00315 Differential_op_typ a_abstract;
00316 a_abstract = a_.Abstract_differential_operator();
00317 if(a_abstract==not_abstract) {
00318 cout << " This is not an abstract differential operator! Error!" << endl;
00319 }
00320 if(a_abstract==abstract_diag)
00321 cout<<"\n Please use operator [] instead of operator ()! Error!" << endl;
00322 typedef DResDiff<DExpr<A> > ExprT;
00323 if(used==false) {
00324 used = true;
00325 return DExpr<ExprT>(ExprT(a_,num_stencils,grid,true));
00326 }
00327 else
00328 return DExpr<ExprT>(ExprT(a_,num_stencils,grid,false));
00329 }
00330
00331
00332 template<class A>
00333 DExpr<DResDiagDiff<DExpr<A> > >
00334 operator[] (const DExpr<A>& a_) {
00335 Differential_op_typ a_abstract;
00336 a_abstract = a_.Abstract_differential_operator();
00337 if(a_abstract==not_abstract) {
00338 cout << " This is not an abstract differential operator! Error!" << endl;
00339 }
00340 if(a_abstract==abstract_with_var)
00341 cout<<"\n Please use operator () instead of operator []! Error!" << endl;
00342 typedef DResDiagDiff<DExpr<A> > ExprT;
00343 return DExpr<ExprT>(ExprT(a_,num_stencils,grid,this));
00344 }
00345 };
00346
00347
00349
00351
00352 template<class A>
00353 void DResDiagDiff<A>::Put_grid_rbo(Grid* gr, int r_bo) const {
00354 a_.Put_grid_rbo(gr,r_bo);
00355 if(ResS->Give_used()==false) {
00356 cout << "\n An abstract diagonal differential operator without";
00357 cout << "\n a variable cannot be restricted! Error! " << endl;
00358 }
00359 }
00360
00361
00362 template<class A>
00363 DResDiff<A>::DResDiff(const A& a, int num_s, Grid* gr, bool isRestricted)
00364 : a_(a), num_stencils(num_s), grid(gr) {
00365 number_variable = a_.Give_number_var_of_abstract_op();
00366 r_bou = a_.run_boundary();
00367
00368 if(isRestricted) {
00369 Stencils_fine = gr->Give_Stencils_fine();
00370 u_ablage = gr->Give_vector_ablage();
00371 label_on_cell = gr->Give_label_on_cell_ablage();
00372 corner_in_domain = gr->Give_eight_bools();
00373
00374 question_calc_stencil = true;
00375 }
00376 else {
00377 question_calc_stencil = false;
00378 }
00379 }
00380
00381 template<class A>
00382 void DResDiff<A>::Put_grid_rbo(Grid* gr, int r_bov) const {
00383 a_.Put_grid_rbo(gr,r_bov);
00384
00385 if(question_calc_stencil == true) {
00386 DResDiff<A>::Iterate_Calc_stencil(r_bov);
00387 }
00388 };
00389
00390
00391 template<class A>
00392 void DResDiff<A>::Iterate_Calc_stencil(int r_bov) const {
00393 int i, j, jj, level, lev;
00394 Celltype ce_typ;
00395 Index3D I_son, I_ecke;
00396 BoCell* b_cell;
00397 int num_v, num_u, dir_u;
00398 double* sp;
00399 double* Stencil_coarse;
00400 double mesh_size;
00401 double finest_mesh_size;
00402 double a,b, ss, hh;
00403 double weight_bocellpoint[8];
00404 dir_3D d;
00405 Index3D my_lev_index;
00406 D3vector Co;
00407 bool u_in_space, uu_in_space, v_in_space;
00408 bool cell_point_in_space;
00409 Pointtype type;
00410 int *coarse_grid_functions;
00411
00412 Index3D I;
00413
00414 coarse_grid_functions = grid->Give_coarse_grid_functions();
00415
00416
00417 int i_iter, hashtable0_leng;
00418 Point_hashtable0* point0;
00419 Point_hashtable0** hashtable0_start;
00420 hashtable0_leng = grid->Give_hashtable0_leng();
00421 hashtable0_start = grid->Give_hashtable0_start();
00422
00423
00424
00425
00426
00427 finest_mesh_size = grid->Give_finest_mesh_size();
00428 level=grid->Max_level();
00429 Iterate_hash0 {
00430
00431 I = point0->Give_Index();
00432 if(I.Cell_index()) {
00433 lev = I.Tiefe();
00434 if((point0->Give_cell_typ()==bo_cell ||
00435 point0->Give_cell_typ()==int_cell) && lev == level) {
00436
00437
00438 mesh_size = grid->H_mesh() / Zweipotenz(I.Tiefe());
00439 for(i=0;i<8;++i) {
00440
00441 I_son = I.son((dir_sons)i);
00442 ce_typ = grid->Give_cell_typ(I_son);
00443 if(ce_typ==int_cell) {
00444
00445 for(j=0;j<64;++j)
00446 Stencils_fine[i][j]
00447 = a_.Give_interior_sten_element(I_son,mesh_size,j,
00448 grid,I.Tiefe());
00449 }
00450 else if(ce_typ==fine_bo_cell) {
00451
00452
00453 for(j=0;j<64;++j) Stencils_fine[i][j] = 0.0;
00454
00455
00456 b_cell = grid->Give_Bo_cell(I_son);
00457
00458
00459 if(b_cell->Exists_bocellpoint()) {
00460 Co = b_cell->Coord_bocellpoint_normal();
00461 weight_bocellpoint[WSDd] =
00462 (1.0 - Co.x) * (1.0 - Co.y) * (1.0 - Co.z);
00463 weight_bocellpoint[ESDd] =
00464 Co.x * (1.0 - Co.y) * (1.0 - Co.z);
00465 weight_bocellpoint[WNDd] =
00466 (1.0 - Co.x) * Co.y * (1.0 - Co.z);
00467 weight_bocellpoint[ENDd] =
00468 Co.x * Co.y * (1.0 - Co.z);
00469 weight_bocellpoint[WSTd] =
00470 (1.0 - Co.x) * (1.0 - Co.y) * Co.z;
00471 weight_bocellpoint[ESTd] =
00472 Co.x * (1.0 - Co.y) * Co.z;
00473 weight_bocellpoint[WNTd] =
00474 (1.0 - Co.x) * Co.y * Co.z;
00475 weight_bocellpoint[ENTd] =
00476 Co.x * Co.y * Co.z;
00477 }
00478
00479
00480 for(dir_u=0;dir_u<8;++dir_u) {
00481 if(r_bou==1) {
00482
00483 u_in_space = true;
00484 }
00485 else if(r_bou==0) {
00486
00487 u_in_space = false;
00488 }
00489 else {
00490
00491 I_ecke = I_son.neighbour((dir_sons)dir_u);
00492 if(grid->Give_global_type(I_ecke)==multigrid) {
00493 if(grid->Give_label_bo_D(I_ecke,-r_bou,level))
00494 u_in_space = false;
00495 else
00496 u_in_space = true;
00497 }
00498 }
00499
00500 for(num_u=0;num_u<b_cell->Give_number_points();++num_u) {
00501 j = b_cell->corner(num_u);
00502 if(b_cell->edge_point(num_u) == corner_poi_typ) {
00503
00504 if(j==dir_u) u_ablage[num_u] = 1.0;
00505 else u_ablage[num_u] = 0.0;
00506 if(r_bou<=0) {
00507
00508 if(j==dir_u) u_in_space = true;
00509 }
00510 }
00511 else {
00512
00513 if(r_bou==1) {
00514
00515 uu_in_space = true;
00516 d = b_cell->edge_dir(num_u);
00517
00518 }
00519 else if(r_bou<0) {
00520
00521 d = b_cell->edge_dir(num_u);
00522 uu_in_space =
00523 grid->Give_label_bo(I_son.neighbour((dir_sons)j),
00524 d,
00525 -r_bou);
00526 }
00527 else {
00528
00529 uu_in_space = false;
00530 }
00531 if(uu_in_space) {
00532 if(j==dir_u) a = 1.0;
00533 else a = 0.0;
00534 if(Next_corner((dir_sons)j,d)==dir_u) b = 1.0;
00535 else b = 0.0;
00536 u_ablage[num_u] =
00537 (a + grid->Give_h(I_son.neighbour((dir_sons)j),d)
00538 / finest_mesh_size
00539 * (b-a));
00540 }
00541 else {
00542 u_ablage[num_u] = 0.0;
00543 }
00544 }
00545 }
00546 if(b_cell->Exists_bocellpoint()) {
00547 if(r_bou==1) {
00548
00549 u_in_space = true;
00550 }
00551 else if(r_bou<0) {
00552
00553 if(b_cell->Is_Dirichlet(-r_bou)==false)
00554 u_in_space = true;
00555 }
00556
00557 num_u = b_cell->Give_number_points();
00558 u_ablage[num_u] = weight_bocellpoint[dir_u];
00559 }
00560
00561 if(u_in_space) {
00562
00563
00564 for(num_v=0;num_v<b_cell->Give_number_points();++num_v) {
00565
00566 j = b_cell->corner(num_v);
00567 if(b_cell->edge_point(num_v) == corner_poi_typ) {
00568
00569 Stencils_fine[i][j+8*dir_u]
00570 += a_.Give_boundary_sten_element(grid,b_cell,
00571 u_ablage,num_v);
00572 }
00573 else {
00574
00575 if(r_bov==1) {
00576
00577 d = b_cell->edge_dir(num_v);
00578 v_in_space = true;
00579 }
00580 else if(r_bov<0) {
00581
00582 d = b_cell->edge_dir(num_v);
00583 v_in_space =
00584 grid->Give_label_bo(I_son.neighbour((dir_sons)j),
00585 d,
00586 -r_bov);
00587 }
00588 else {
00589
00590 v_in_space = false;
00591 }
00592 if(v_in_space) {
00593 ss = a_.Give_boundary_sten_element(grid,b_cell,
00594 u_ablage,num_v);
00595 d = b_cell->edge_dir(num_v);
00596 j = b_cell->corner(num_v);
00597
00598 hh = grid->Give_h(I_son.neighbour((dir_sons)j),d)
00599 / finest_mesh_size;
00600
00601 Stencils_fine[i][j+8*dir_u]
00602 += ss * (1.0 + hh * (0.0-1.0));
00603
00604 jj = Next_corner((dir_sons)j,d);
00605 Stencils_fine[i][jj+8*dir_u]
00606 += ss * (0.0 + hh * (1.0-0.0));
00607 }
00608 }
00609 }
00610
00611 if(b_cell->Exists_bocellpoint()) {
00612
00613 num_v = b_cell->Give_number_points();
00614 ss = a_.Give_boundary_sten_element(grid,b_cell,
00615 u_ablage,num_v);
00616 if(r_bov==1) {
00617
00618 cell_point_in_space = true;
00619 }
00620 else if(r_bov<0) {
00621
00622 if(b_cell->Is_Dirichlet(-r_bov)==false)
00623 cell_point_in_space = true;
00624 else
00625 cell_point_in_space = false;
00626 }
00627 else
00628
00629 cell_point_in_space = false;
00630 if(cell_point_in_space) {
00631
00632 for(j=0;j<8;++j)
00633 Stencils_fine[i][j+8*dir_u] += ss * weight_bocellpoint[j];
00634 }
00635 else {
00636
00637 for(num_v=0;num_v<b_cell->Give_number_points();++num_v) {
00638 if(b_cell->edge_point(num_v) == corner_poi_typ) {
00639 j = b_cell->corner(num_v);
00640 Stencils_fine[i][j+8*dir_u] += ss * weight_bocellpoint[j];
00641 }
00642 }
00643 }
00644 }
00645 }
00646 }
00647 }
00648 else {
00649 for(j=0;j<64;++j) Stencils_fine[i][j] = 0.0;
00650 }
00651 }
00652
00653
00654
00655 Stencil_coarse = grid->Give_stencil(I,num_stencils);
00656 Restrict_local_stiffness_matrix(Stencil_coarse,
00657 Stencils_fine,coarse_grid_functions);
00658
00659 if(point0->Give_cell_typ()==bo_cell) {
00660
00661 if(r_bov==0) {
00662 for(i=0;i<8;++i)
00663 if(grid->Give_global_type(I.neighbour((dir_sons)i))==multigrid) {
00664 for(j=0;j<8;++j) {
00665 Stencil_coarse[i+8*j] = 0.0;
00666 }
00667 }
00668 }
00669
00670 if(r_bou==0) {
00671 for(i=0;i<8;++i)
00672 if(grid->Give_global_type(I.neighbour((dir_sons)i))==multigrid) {
00673 for(j=0;j<8;++j) {
00674 Stencil_coarse[j+8*i] = 0.0;
00675 }
00676 }
00677 }
00678
00679
00680 if(r_bov<0) {
00681 for(i=0;i<8;++i) {
00682 type = grid->Give_global_type(I.neighbour((dir_sons)i));
00683 if(type==multigrid) {
00684 if(grid->Give_label_bo_D(I.neighbour((dir_sons)i),-r_bov,level-1)) {
00685 for(j=0;j<8;++j) {
00686 Stencil_coarse[i+8*j] = 0.0;
00687 }
00688 }
00689 }
00690 if(type==exterior) {
00691 for(j=0;j<8;++j) {
00692 Stencil_coarse[i+8*j] = 0.0;
00693 }
00694 }
00695 }
00696 }
00697
00698 if(r_bou<0) {
00699 for(i=0;i<8;++i) {
00700 type = grid->Give_global_type(I.neighbour((dir_sons)i));
00701 if(type==multigrid) {
00702 if(grid->Give_label_bo_D(I.neighbour((dir_sons)i),-r_bou,level-1)) {
00703 for(j=0;j<8;++j) {
00704 Stencil_coarse[j+8*i] = 0.0;
00705 }
00706 }
00707 }
00708 if(type==exterior) {
00709 for(j=0;j<8;++j) {
00710 Stencil_coarse[j+8*i] = 0.0;
00711 }
00712 }
00713 }
00714 }
00715 }
00716 }
00717 }
00718 }
00719
00720 if(parallel_version) {
00721 grid->Update_global_stencil(num_stencils,grid->Max_level());
00722 }
00723
00724
00725
00726 for(level=grid->Max_level()-1;level>grid->Min_level();--level) {
00727 my_lev_index = grid->Give_my_level_index(level);
00728 Iterate_hash0 {
00729
00730 I = point0->Give_Index();
00731 if(I.Cell_index()) {
00732 lev = I.Tiefe();
00733 if((point0->Give_cell_typ()==bo_cell ||
00734 point0->Give_cell_typ()==int_cell) && lev == level &&
00735 I < my_lev_index) {
00736
00737
00738
00739 mesh_size = grid->H_mesh() / Zweipotenz(I.Tiefe());
00740 for(i=0;i<8;++i) {
00741 I_son = I.son((dir_sons)i);
00742 ce_typ = grid->Give_cell_typ(I_son);
00743
00744 if(ce_typ==int_cell || ce_typ==bo_cell) {
00745 sp = grid->Give_stencil(I_son,num_stencils);
00746 for(j=0;j<64;++j) Stencils_fine[i][j] = sp[j];
00747 }
00748 else {
00749 for(j=0;j<64;++j) Stencils_fine[i][j] = 0.0;
00750 }
00751 }
00752
00753
00754 Stencil_coarse = grid->Give_stencil(I,num_stencils);
00755 Restrict_local_stiffness_matrix(Stencil_coarse,
00756 Stencils_fine,coarse_grid_functions);
00757
00758 if(point0->Give_cell_typ()==bo_cell) {
00759
00760 if(r_bov==0) {
00761 for(i=0;i<8;++i)
00762 if(grid->Give_global_type(I.neighbour((dir_sons)i))==multigrid) {
00763 for(j=0;j<8;++j) {
00764 Stencil_coarse[i+8*j] = 0.0;
00765 }
00766 }
00767 }
00768
00769 if(r_bou==0) {
00770 for(i=0;i<8;++i)
00771 if(grid->Give_global_type(I.neighbour((dir_sons)i))==multigrid) {
00772 for(j=0;j<8;++j) {
00773 Stencil_coarse[j+8*i] = 0.0;
00774 }
00775 }
00776 }
00777
00778
00779 if(r_bov<0) {
00780 for(i=0;i<8;++i) {
00781 type = grid->Give_global_type(I.neighbour((dir_sons)i));
00782 if(type==multigrid) {
00783 if(grid->Give_label_bo_D(I.neighbour((dir_sons)i),-r_bov,level-1)) {
00784 for(j=0;j<8;++j) {
00785 Stencil_coarse[i+8*j] = 0.0;
00786 }
00787 }
00788 }
00789 if(type==exterior) {
00790 for(j=0;j<8;++j) {
00791 Stencil_coarse[i+8*j] = 0.0;
00792 }
00793 }
00794 }
00795 }
00796
00797 if(r_bou<0) {
00798 for(i=0;i<8;++i) {
00799 type = grid->Give_global_type(I.neighbour((dir_sons)i));
00800 if(type==multigrid) {
00801 if(grid->Give_label_bo_D(I.neighbour((dir_sons)i),-r_bou,level-1)) {
00802 for(j=0;j<8;++j) {
00803 Stencil_coarse[j+8*i] = 0.0;
00804 }
00805 }
00806 }
00807 if(type==exterior) {
00808 for(j=0;j<8;++j) {
00809 Stencil_coarse[j+8*i] = 0.0;
00810 }
00811 }
00812 }
00813 }
00814 }
00815 }
00816 }
00817 }
00818
00819 if(parallel_version) {
00820 grid->Update_global_stencil(num_stencils,level);
00821 }
00822 }
00823 }
00824
00825 #endif
00826