AMROC - Source code documentation

00001 #ifndef _included_GridFunctionBndry1_c
00002 #define _included_GridFunctionBndry1_c
00003 
00009 template <class DAGH_GFType>
00010 void GridFunction(1)<DAGH_GFType>::GF_GatherGhostBndryInfo() {
00011 
00012   const int levels = dagh.totallevels();
00013   GridBoxList** ggbl = dagh.ggbl();
00014   const int crslev = dagh.coarselevel();
00015   const int myrank = 1;
00016   int l;
00017 
00018   bndrybboxlist = new BBoxList*[levels];
00019   for (l=0; l<levels; l++) 
00020     bndrybboxlist[l] = (BBoxList*) 0;
00021   prolongbboxlist = new BBoxList*[levels];
00022   for (l=0; l<levels; l++) 
00023     prolongbboxlist[l] = (BBoxList*) 0;
00024 
00025   if (bwidth == 0) return;
00026 
00027   for (l=0; l<levels; l++) 
00028     if (l>=Max(crslev,crslev+loffset) && l<Min(levels,levels+loffset)) {
00029       if (!ggbl[l]) continue;
00030 
00031       BBox wbb(dagh.wholebbox());
00032       wbb.refine(dagh.refinedby(l));
00033       BBoxList wbbl;
00034       wbbl.add(grow(usedbbox(wbb,l,DAGH_All),bwidth));
00035 
00036       for (register GridBox* gb=ggbl[l]->first();gb;gb=ggbl[l]->next()) 
00037         wbbl -= usedbbox(*gb,l,DAGH_All);
00038 
00039       if (externalboundary()) {
00040         bndrybboxlist[l] = new BBoxList();
00041         *bndrybboxlist[l] = wbbl;
00042         if (l>0 && ggbl[l-1]) {
00043           for (register GridBox* gb=ggbl[l-1]->first();gb;gb=ggbl[l-1]->next()) {
00044             BBox bb(gb->gbBBox());
00045             bb.refine(dagh.refinefactor(l-1));
00046             *bndrybboxlist[l] -= usedbbox(bb,l,DAGH_All);
00047           }
00048           bndrybboxlist[l]->mergeboxes(overlap);
00049         }
00050       }
00051 
00052       if (adaptiveboundary()) {
00053         prolongbboxlist[l] = new BBoxList();
00054         *prolongbboxlist[l] = wbbl;
00055         if (bndrybboxlist[l]) 
00056           *prolongbboxlist[l] -= *bndrybboxlist[l];
00057         prolongbboxlist[l]->mergeboxes(overlap);
00058       }
00059 
00060       if (externalboundary()) {
00061         if (alignment != DAGH_All) 
00062           for (register BBox* bb=bndrybboxlist[l]->first();bb;
00063                bb=bndrybboxlist[l]->next()) 
00064             gdbAlignBBox(myrank, *bb, alignment);
00065 #ifdef DEBUG_PRINT_GF_COMM
00066         ( comm_service::log() << "GF_GatherGhostBndryInfo: bndrybboxlist on " << l 
00067           << *bndrybboxlist[l] << endl ).flush();
00068 #endif    
00069       }
00070       if (adaptiveboundary()) {
00071         if (alignment != DAGH_All) 
00072           for (register BBox* bb=prolongbboxlist[l]->first();bb;
00073                bb=prolongbboxlist[l]->next()) 
00074             gdbAlignBBox(myrank, *bb, alignment);
00075 #ifdef DEBUG_PRINT_GF_COMM
00076         ( comm_service::log() << "GF_GatherGhostBndryInfo: prolongbboxlist on " << l 
00077           << *prolongbboxlist[l] << endl ).flush();
00078 #endif    
00079       }
00080     }
00081 }
00082 
00083 
00084 template <class DAGH_GFType>
00085 void GridFunction(1)<DAGH_GFType>::GF_DeleteGhostBndryInfo() {
00086 
00087   const int levels = dagh.totallevels();
00088   register int l;
00089 
00090   if (bndrybboxlist) {
00091     for (l=0; l<levels; l++) 
00092       if (bndrybboxlist[l]) delete bndrybboxlist[l];
00093     delete [] bndrybboxlist;
00094   }
00095   if (prolongbboxlist) {
00096     for (l=0; l<levels; l++) 
00097       if (prolongbboxlist[l]) delete prolongbboxlist[l];
00098     delete [] prolongbboxlist;
00099   }
00100 }
00101 
00102 
00103 /*****************************************************************************/
00104 /**** Adaptive Boundary Update --- Time-Space Interpolation  ****/
00105 /*****************************************************************************/
00106 template <class DAGH_GFType>
00107 void GridFunction(1)<DAGH_GFType>::GF_AdaptiveBndryUpdate(const int time, 
00108                                                               const int level) {
00109   if (bwidth == 0) return;
00110 
00111   register int const t = dagh_timeindex(time,level);
00112   register int const l = level;
00113   
00114   if (!gdb[t] || !gdb[t][l] || l<=loffset || !adaptiveboundary()) return;
00115 
00116   GDB_Interaction* pi;
00117   int npi;
00118   int tc = GF_CurrentTime(l-1);
00119   int sc = GF_TimeStep(l-1);
00120   int tnc, tpc;
00121   double frac;
00122   
00123   if (time < tc) {
00124     frac = 1.0*(time-(tc-sc))/sc;
00125     tpc = dagh_timeindex(tc-sc,l-1);
00126     tnc = dagh_timeindex(tc,l-1);      
00127   }
00128   else if (time == tc && updatedstep == DAGHNextTime) {
00129     frac = 0.0; 
00130     tpc = dagh_timeindex(tc,l-1);
00131     tnc = dagh_timeindex(tc+sc,l-1);      
00132   }
00133   else if (time == tc && updatedstep == DAGHCurrentTime) {      
00134     frac = 1.0;
00135     tpc = dagh_timeindex(tc-sc,l-1);
00136     tnc = dagh_timeindex(tc,l-1);      
00137   }
00138   else {
00139     frac = 1.0*(time-tc)/sc;
00140     tpc = dagh_timeindex(tc,l-1);
00141     tnc = dagh_timeindex(tc+sc,l-1);
00142   }
00143   double oneminusfrac = 1.0-frac;
00144 #ifdef DEBUG_PRINT_GF_COMM
00145   ( comm_service::log() << "GF_AdaptiveBndryUpdate:" << time << "   " << tc << "   " 
00146     << frac << "   " << tpc << "   " << tnc << endl ).flush();
00147 #endif    
00148 #ifdef DEBUG_PRINT
00149   assert (frac>=0 && frac<=1.0);
00150 #endif    
00151 
00152   if (frac>0.0 && frac<1.0) {
00153 #ifdef DEBUG_PRINT
00154     assert (gdb[tnc][l-1] && gdb[tpc][l-1]);
00155 #endif    
00156     for (register int i=0; i<length[l]; i++)
00157       if (gdb[tnc][l-1] && gdb[tpc][l-1] && gdb[t][l][i]) 
00158         if ((npi=gdb[t][l][i]->gdb_parent_index)>0 && (pi=gdb[t][l][i]->gdb_parent_info) &&
00159             gdb[t][l][i]->has_adaptiveboundaries()) {
00160           GridData(1)<DAGH_GFType> &pl = gdb[t][l][i]->griddata();
00161           GridData(1)<DAGH_GFType> plpf(pl.bbox());
00162           GridData(1)<DAGH_GFType> plnf(pl.bbox());
00163           
00164           /* Interpolation on overlapping boundaries is not commutative. */
00165           /* We need a temporary work-array for pl. */
00166           GridData(1)<DAGH_GFType> work(pl.bbox()); work.copy(pl);
00167         
00168           for (register int j=0; j<npi; j++)
00169             if (gdb[tnc][l-1][pi[j].idx] && gdb[tpc][l-1][pi[j].idx]) {
00170               GridData(1)<DAGH_GFType> &palnc = gdb[tnc][l-1][pi[j].idx]->griddata();
00171               GridData(1)<DAGH_GFType> &palpc = gdb[tpc][l-1][pi[j].idx]->griddata();
00172           
00173               /* Call the prolongation function */
00174               for(register int b=0; b<2*gfrank; b++) if (gdb[t][l][i]->has_adaptiveboundary(b)) {
00175                 BBox bb = pi[j].bbox*gdb[t][l][i]->adaptivebndrybox(b);
00176                 if (!bb.empty()) {
00177 #ifdef DEBUG_PRINT_GF_COMM
00178                   ( comm_service::log() << "GF_AdaptiveBndryUpdate on :" 
00179                     << bb << endl ).flush();
00180 #endif
00181                   /* Call the prolongation function */
00182                   if (adaptiveboundary_type() == DAGHAdaptBoundaryInterp ||
00183                       adaptiveboundary_type() == DAGHAdaptBoundaryBoth) { 
00184                     if (pfunc) {
00185                       (*pfunc)(FORTRAN_ARGS(1,palpc),FORTRAN_ARGS(1,plpf),
00186                                BOUNDING_BOX(bb),myargs,&myargc);
00187                       (*pfunc)(FORTRAN_ARGS(1,palnc),FORTRAN_ARGS(1,plnf),
00188                                BOUNDING_BOX(bb),myargs,&myargc);
00189                 
00190                       /* data = (old * (1-frac)) + (new * frac) */
00191                       if (frac>0.0 && frac<1.0) 
00192                         work.lin_interp(plpf,(1.0-frac),plnf,frac,bb);
00193                       if (frac==0.0) work.copy(plpf,bb);
00194                       if (frac==1.0) work.copy(plnf,bb);                    
00195                     }
00196                   }
00197 
00198                   /* Call User-defined adaptive boundary update function */
00199                   if (adaptiveboundary_type() == DAGHAdaptBoundaryUserDef ||
00200                       adaptiveboundary_type() == DAGHAdaptBoundaryBoth) {
00201                     if (abfunc) 
00202                       (*abfunc)(FORTRAN_ARGS(1,work),
00203                                 FORTRAN_ARGS(1,palnc),&frac,
00204                                 FORTRAN_ARGS(1,palpc),&oneminusfrac,
00205                                 BOUNDING_BOX(bb),myargs,&myargc);
00206                   }
00207                   
00208                   /* Default is an error */
00209                   if (adaptiveboundary_type() != DAGHAdaptBoundaryUserDef &&
00210                       adaptiveboundary_type() != DAGHAdaptBoundaryInterp &&
00211                       adaptiveboundary_type() != DAGHAdaptBoundaryBoth) 
00212                     assert(0);
00213                 }
00214               }
00215             }
00216           pl.copy(work);
00217         }
00218   }
00219 
00220   else if (frac==0.0) {
00221 #ifdef DEBUG_PRINT
00222     assert (gdb[tpc][l-1]);
00223 #endif    
00224     for (register int i=0; i<length[l]; i++)
00225       if (gdb[tpc][l-1] && gdb[t][l][i]) 
00226         if ((npi=gdb[t][l][i]->gdb_parent_index)>0 && (pi=gdb[t][l][i]->gdb_parent_info) &&
00227             gdb[t][l][i]->has_adaptiveboundaries()) {
00228           GridData(1)<DAGH_GFType> &pl = gdb[t][l][i]->griddata();
00229           GridData(1)<DAGH_GFType> plpf(pl.bbox());
00230           
00231           /* Interpolation on overlapping boundaries is not commutative. */
00232           /* We need a temporary work-array for pl. */
00233           GridData(1)<DAGH_GFType> work(pl.bbox()); work.copy(pl);
00234         
00235           for (register int j=0; j<npi; j++)
00236             if (gdb[tpc][l-1][pi[j].idx]) {
00237               GridData(1)<DAGH_GFType> &palpc = gdb[tpc][l-1][pi[j].idx]->griddata();
00238           
00239               /* Call the prolongation function */
00240               for(register int b=0; b<2*gfrank; b++) if (gdb[t][l][i]->has_adaptiveboundary(b)) {
00241                 BBox bb = pi[j].bbox*gdb[t][l][i]->adaptivebndrybox(b);
00242                 if (!bb.empty()) {
00243 #ifdef DEBUG_PRINT_GF_COMM
00244                   ( comm_service::log() << "GF_AdaptiveBndryUpdate on :" 
00245                     << bb << endl ).flush();
00246 #endif
00247                   /* Call the prolongation function */
00248                   if (adaptiveboundary_type() == DAGHAdaptBoundaryInterp ||
00249                       adaptiveboundary_type() == DAGHAdaptBoundaryBoth) { 
00250                     if (pfunc) {
00251                       (*pfunc)(FORTRAN_ARGS(1,palpc),FORTRAN_ARGS(1,plpf),
00252                                BOUNDING_BOX(bb),myargs,&myargc);
00253                       work.copy(plpf,bb);
00254                     }
00255                   }
00256 
00257                   /* Call User-defined adaptive boundary update function */
00258                   if (adaptiveboundary_type() == DAGHAdaptBoundaryUserDef ||
00259                       adaptiveboundary_type() == DAGHAdaptBoundaryBoth) {
00260                     if (abfunc) 
00261                       (*abfunc)(FORTRAN_ARGS(1,work),
00262                                 FORTRAN_ARGS(1,palpc),&frac,
00263                                 FORTRAN_ARGS(1,palpc),&oneminusfrac,
00264                                 BOUNDING_BOX(bb),myargs,&myargc);
00265                   }
00266                   
00267                   /* Default is an error */
00268                   if (adaptiveboundary_type() != DAGHAdaptBoundaryUserDef &&
00269                       adaptiveboundary_type() != DAGHAdaptBoundaryInterp &&
00270                       adaptiveboundary_type() != DAGHAdaptBoundaryBoth) 
00271                     assert(0);
00272                 }
00273               }
00274             }
00275           pl.copy(work);
00276         }
00277   }
00278 
00279   else if (frac==1.0) {
00280 #ifdef DEBUG_PRINT
00281     assert (gdb[tnc][l-1]);
00282 #endif    
00283     for (register int i=0; i<length[l]; i++)
00284       if (gdb[tnc][l-1] && gdb[t][l][i]) 
00285         if ((npi=gdb[t][l][i]->gdb_parent_index)>0 && (pi=gdb[t][l][i]->gdb_parent_info) &&
00286             gdb[t][l][i]->has_adaptiveboundaries()) {
00287           GridData(1)<DAGH_GFType> &pl = gdb[t][l][i]->griddata();
00288           GridData(1)<DAGH_GFType> plnf(pl.bbox());
00289           
00290           /* Interpolation on overlapping boundaries is not commutative. */
00291           /* We need a temporary work-array for pl. */
00292           GridData(1)<DAGH_GFType> work(pl.bbox()); work.copy(pl);
00293         
00294           for (register int j=0; j<npi; j++)
00295             if (gdb[tnc][l-1][pi[j].idx]) {
00296               GridData(1)<DAGH_GFType> &palnc = gdb[tnc][l-1][pi[j].idx]->griddata();
00297           
00298               /* Call the prolongation function */
00299               for(register int b=0; b<2*gfrank; b++) if (gdb[t][l][i]->has_adaptiveboundary(b)) {
00300                 BBox bb = pi[j].bbox*gdb[t][l][i]->adaptivebndrybox(b);
00301                 if (!bb.empty()) {
00302 #ifdef DEBUG_PRINT_GF_COMM
00303                   ( comm_service::log() << "GF_AdaptiveBndryUpdate on :" 
00304                     << bb << endl ).flush();
00305 #endif
00306                   /* Call the prolongation function */
00307                   if (adaptiveboundary_type() == DAGHAdaptBoundaryInterp ||
00308                       adaptiveboundary_type() == DAGHAdaptBoundaryBoth) { 
00309                     if (pfunc) {
00310                       (*pfunc)(FORTRAN_ARGS(1,palnc),FORTRAN_ARGS(1,plnf),
00311                                BOUNDING_BOX(bb),myargs,&myargc);
00312                       work.copy(plnf,bb);                   
00313                     }
00314                   }
00315 
00316                   /* Call User-defined adaptive boundary update function */
00317                   if (adaptiveboundary_type() == DAGHAdaptBoundaryUserDef ||
00318                       adaptiveboundary_type() == DAGHAdaptBoundaryBoth) {
00319                     if (abfunc) 
00320                       (*abfunc)(FORTRAN_ARGS(1,work),
00321                                 FORTRAN_ARGS(1,palnc),&frac,
00322                                 FORTRAN_ARGS(1,palnc),&oneminusfrac,
00323                                 BOUNDING_BOX(bb),myargs,&myargc);
00324                   }
00325                   
00326                   /* Default is an error */
00327                   if (adaptiveboundary_type() != DAGHAdaptBoundaryUserDef &&
00328                       adaptiveboundary_type() != DAGHAdaptBoundaryInterp &&
00329                       adaptiveboundary_type() != DAGHAdaptBoundaryBoth) 
00330                     assert(0);
00331                 }
00332               }
00333             }
00334           pl.copy(work);
00335         }
00336   }
00337   /* Default is an error */
00338   else
00339     assert(0);
00340 }
00341 
00342 /*****************************************************************************/
00343 /**** External Boundary Update ****/
00344 /*****************************************************************************/
00345 template <class DAGH_GFType>
00346 void GridFunction(1)<DAGH_GFType>::GF_ExternalBndryUpdate(const int time, 
00347                                                               const int level) {
00348   if (bwidth == 0) return;
00349 
00350   register int const t = dagh_timeindex(time,level);
00351   register int const l = level;
00352   double pht = get_phystime_timevalue(time,level);
00353   
00354   if (!gdb[t] || !gdb[t][l] || !externalboundary()) return;
00355 
00356   for (register int i=0; i<length[l]; i++) if (gdb[t][l][i]) {
00357     DCoords lbc = dagh.worldCoords(gdb[t][l][i]->boundingbox().lower(), 
00358                                    gdb[t][l][i]->boundingbox().stepsize());
00359     DCoords dx(dagh.rank);
00360     for(int d=0; d<dagh.rank; d++) 
00361       dx(d) = GF_DeltaX(d,l);     
00362     
00363     GridData(1)<DAGH_GFType> &u = gdb[t][l][i]->griddata();
00364     for(register int b=0; b<2*gfrank; b++) if (gdb[t][l][i]->has_externalboundary(b)) 
00365       for (register BBox* bb=gdb[t][l][i]->externalbndrylist(b)->first();bb;
00366            bb=gdb[t][l][i]->externalbndrylist(b)->next()) {
00367 
00368 #ifdef DEBUG_PRINT_GF_COMM
00369         ( comm_service::log() << "GF_ExternalBndryUpdate on :" 
00370           << *bb << endl ).flush();
00371 #endif    
00372         /* Use User-defined External boundary function */
00373         if (boundary_type() == DAGHBoundaryUserDef && bfunc) 
00374           (*bfunc)(FORTRAN_ARGS(1,u),BOUNDING_BOX(*bb),lbc(),dx(),&b,
00375                    dagh.wholebndry(),&dagh.nbndry(),&pht,myargs,&myargc);
00376 
00377         /* H1e type shift boundaries */
00378         else if (boundary_type() == DAGHBoundaryShift) 
00379           u.equals(u,*bb,shift(*bb,b/2,(b%2 ? -1 : 1)));
00380 
00381         /* Set the boundary to user defined constant */
00382         else if (boundary_type() == DAGHBoundaryRegular) 
00383           u.equals(bvalue,*bb);
00384 
00385         /* Default is an error */
00386         else
00387           assert(0);
00388       }
00389   }
00390 }
00391 
00392 
00393 #endif
GridFunctionBndry1.h

last update: 06/01/04
