AdaptiveGridRefine.hpp

00001 extern int debugRefine ;
00002 
00003 template<int DIM>
00004 void AdaptiveGrid<DIM>::Set(Function *F , double c) {
00005  double dl[DIM] ;
00006  int    i,anz=0 ;
00007  
00008  // save BasisPG[0] for later refinement of instances of AdaptiveData
00009  AdaptiveGrid<DIM>::ProjectionGrid  *pg=BasisPG[0] ;
00010  BasisPG[0]=OldBasisPG ;
00011  OldBasisPG=pg ;
00012  ClearPG(BasisPG[0]) ;
00013  CounterArray.Set((size_t)0) ;
00014 
00015  for (i=0; i<DIM; i++) MaxLevel[i]=LMAX ;
00016  AdaptiveDataArray<DIM>::VectorArray::iterator l(Level0 , MaxLevel) ;
00017 
00018  for (l.Set(Level0); l<=MaxLevel; ++l) {    
00019 
00020    for (i=0; i<DIM; i++) dl[i]=(double)l.i[i] ; 
00021 
00022    if (F->Eval(dl) <= c ) {
00023      int a[DIM]={0}, e[DIM] ;
00024      for (i=0; i<DIM; i++) e[i] = (l.i[i]==Level0[i]) ? (1<<l.i[i]) : (1<<(l.i[i]-1)) -1 ;
00025      Matrix<int,DIM>::iterator j(a,e) ;
00026      for (j.Set(a); j<=e; ++j) {
00027        Insert(BasisPG[0],0, l.i,j.i ) ; // insert and update CounterArray
00028        anz++ ;
00029      }
00030    }
00031  }
00032 
00033  // generate other PGs
00034  GenerateAllBasisPGs() ;
00035  was_refined=true ;
00036 
00037  TmpBasis->Resize(&CounterArray) ;
00038  for (i=0; i<ResizeGroup.count; i++) ResizeGroup.Ds[i]->ba->Resize(&CounterArray) ;
00039  for (i=0; i<RefineGroup.count; i++) RefineGroup.Ds[i]->Refine() ;
00040 
00041  SMPLoadBalance() ;
00042 
00043  // FaceGrids
00044  for (i=0;i<DIM; i++) {
00045    if (FaceGrid[i][0]) FaceGrid[i][0]->SetAsSlice(this,i,Level0[i],0) ;
00046    if (FaceGrid[i][1]) FaceGrid[i][1]->SetAsSlice(this,i,Level0[i],1<<Level0[i]) ;
00047  }
00048  
00049 }
00050 
00051 
00052 template<int DIM>
00053 void AdaptiveGrid<DIM>::SetSparse(int sl) {
00054  SparseGridFunction<DIM> F ;
00055  Set(&F,sl) ;
00056 }
00057 
00058 
00059 template<int DIM>
00060 void AdaptiveGrid<DIM>::SetSparse(int sl , int *MxLevel) {
00061  SparseGridFunction<DIM> F ;
00062  double p[DIM] ;
00063  int i ;
00064  for (i=0;i<DIM; i++) p[i]=MxLevel[i] ;
00065  F.SetParameters(p) ; 
00066  
00067  Set(&F,sl) ;
00068 }
00069 
00070 template<int DIM>
00071 void AdaptiveGrid<DIM>::SetFull(int JM) {
00072  FullGridFunction<DIM> F ;
00073  double p[DIM] ;
00074  int i ;
00075  for (i=0;i<DIM; i++) p[i]=JM ;
00076  F.SetParameters(p) ; 
00077  
00078  Set(&F,JM) ; 
00079 }
00080 
00081 template<int DIM>
00082 void AdaptiveGrid<DIM>::SetFull(int *JM) {
00083  FullGridFunction<DIM> F ;
00084  double p[DIM] ;
00085  int i ;
00086  for (i=0;i<DIM; i++) p[i]=JM[i] ;
00087  F.SetParameters(p) ; 
00088  
00089  Set(&F,MAXINT-1) ;
00090 }
00091 
00092 /***********************************************************/
00093 /*                                                         */
00094 /*  SetAsSlice of a (DIM+1)-deimensional AdaptiveGrid      */
00095 /*  The main idea is simply to run over the ProjectionGrid */
00096 /*  for the <dir>th coordiante direction and to check      */
00097 /*  whether (l,t) is contained in the AdaptiveLevels       */
00098 /*                                                         */
00099 /***********************************************************/
00100 
00101 template<int DIM>
00102 void AdaptiveGrid<DIM>::SetAsSlice(AdaptiveGrid<DIM+1> *G, int dir, int l, int t) {
00103  int i,j ;
00104  
00105  SliceParent=G   ;
00106  SliceDir   =dir ;
00107  Slicel     =l   ;
00108  Slicet     =t   ;
00109 
00110  AdaptiveGrid<DIM+1>::ProjectionGrid           *pg=G->BasisPG[dir] ;
00111  AdaptiveGrid<DIM+1>::ProjectionGrid::iterator  pgit,pgend=(*pg).end() ;
00112 
00113  // save BasisPG[0] for later refinement of instances of AdaptiveData
00114  AdaptiveGrid<DIM>::ProjectionGrid *opg=BasisPG[0] ;
00115  BasisPG[0]=OldBasisPG ;
00116  OldBasisPG=opg        ;
00117  ClearPG(BasisPG[0])   ;
00118  CounterArray.Set((size_t)0) ;
00119 
00120  for (pgit=(*pg).begin(); pgit!=pgend; ++pgit) {
00121    AdaptiveLevels::AdaptiveL *al =(*(*pgit).second)[l] ;
00122 
00123    if ((*al).find(t) != (*al).end()) { // insert
00124      AdaptiveLevelsIndex *ali=(*pgit).first  ;
00125      int ll[DIM],tt[DIM] ;
00126      for (i=0; i<DIM; i++) ll[i]=(int)ali->l[i], tt[i]=(int)ali->t[i] ;
00127      Insert(BasisPG[0],0, ll,tt ) ; // insert and update CounterArray
00128    }
00129  }
00130 
00131  GenerateAllBasisPGs() ;
00132  was_refined=true ;
00133   
00134  // set internal data
00135  WC=G->WC ;
00136  IsPeriodicValid=G->IsPeriodicValid ;
00137  j=0;
00138  for (i=0; i<DIM+1; i++)  
00139    if (i!=dir) {
00140      XA        [j]=G->XA[i] ;
00141      XE        [j]=G->XE[i] ;
00142      Level0    [j]=G->Level0[j] ;
00143      IsPeriodic[j]=G->IsPeriodic[j] ;
00144      j++ ;
00145    }
00146 
00147  // resize attached AdaptiveData
00148  TmpBasis->Resize(&CounterArray) ;
00149  for (i=0; i<ResizeGroup.count; i++) ResizeGroup.Ds[i]->ba->Resize(&CounterArray) ;
00150  for (i=0; i<RefineGroup.count; i++) RefineGroup.Ds[i]->Refine() ;
00151 
00152  SMPLoadBalance() ;
00153 
00154  // FaceGrids
00155  for (i=0;i<DIM; i++) {
00156    if (FaceGrid[i][0]) FaceGrid[i][0]->SetAsSlice(this,i,Level0[i],0) ;
00157    if (FaceGrid[i][1]) FaceGrid[i][1]->SetAsSlice(this,i,Level0[i],1<<Level0[i]) ;
00158  }
00159 
00160 }
00161 
00162 
00163 /**********************************************************/
00164 /*                                                        */
00165 /*  Grid-Refinement:                                      */
00166 /*  The main idea is simply to employ the tensor product  */
00167 /*  of the univariate refinement function                 */
00168 /*                                                        */
00169 /**********************************************************/
00170 
00171 template<int DIM>
00172 void AdaptiveGrid<DIM>::Refine(AdaptiveData<DIM> *OD , AdaptivityCriterion *R, AdaptiveData<DIM> *ActiveEntries) {
00173  int    i,ld,dir,l[DIM],MxLevel=R->MaxLevel ;
00174  size_t noe,sall ;
00175 
00176  assert(IsPeriodicValid) ;
00177 
00178  AdaptiveGrid<DIM>::ProjectionGrid           *pg,*npg    ;
00179  AdaptiveGrid<DIM>::ProjectionGrid::iterator  pgit,pgend ;
00180 
00181  AdaptivityCriterion R0(AdaptivityCriterion::ANY,-1.0), *RR ;
00182 
00183  for (dir=0; dir<DIM; dir++) {
00184 
00185    // set ANY for adaptivity in directions dir>0 
00186    // but do also refinement on higher levels 
00187    RR = (dir==0) ? R : &R0 ;
00188    RR->MaxLevel=min(MxLevel , MaxLevel[dir]) ;
00189 
00190    //if (debugRefine) printf("Refine %d %d %d %d\n",dir,Level0[dir], AB1[0]->Level0 , AB2[0]->Level0 ) ;
00191 
00192    // set BC
00193    AB2[0]->BC[0]=AB2[0]->BC[1]=AB1[0]->BC[0]=AB1[0]->BC[1]=-1 ;
00194    if (IsPeriodic[dir]) AB2[0]->BC[1]=AB1[0]->BC[1]=PERIODIC ;
00195   
00196    AB1[0]->islifting=AB2[0]->islifting=OD->Ext.islifting[dir] ;
00197 
00198    // get source/target pg    
00199     pg=BasisPG[dir] ;
00200    npg=OldBasisPG   ;
00201    assert(pg!=npg)  ;
00202 
00203    // clear target
00204    ClearPG(npg) ;
00205    CounterArray.Set((size_t)0) ;
00206 
00207    pgend=(*pg).end() ;
00208    for (pgit=(*pg).begin(); pgit!=pgend; ++pgit) {
00209      AdaptiveLevelsIndex *ali=(*pgit).first ,*nali ;
00210      AdaptiveLevels      *als=(*pgit).second,*nals ;
00211           
00212      // get level, requierd by AdptiveB::Refine
00213      for (i=0  ; i<dir  ; i++) l[i  ]=ali->l[i] ;
00214      for (i=dir; i<DIM-1; i++) l[i+1]=ali->l[i] ;
00215 
00216      // load IndexSets/Data for Basis, if dir >0 only IndexSet necessary !!
00217      for (ld=Level0[dir]; ld<=LMAX; ld++) {
00218        if (dir==0) {
00219          l[dir]=ld ;
00220          AB1[0]->L[ld].IS->Load((*als)[ld], OD->ba->d[l], AB1[0]->L[ld].d) ;
00221        } else  
00222          AB1[0]->L[ld].IS->Load((*als)[ld]) ; // no data
00223      }
00224  
00225      // Refine
00226      AB2[0]->Refine(AB1[0], RR, &l[0], dir, DIM, WC, false) ;
00227      // insert AdaptiveLevels corresponding to psgit.i in npg
00228      if (AB2[0]->Size() >0) {
00229 
00230        nals=AdaptiveLevelsAllocator     .NewItem() ;
00231        nali=AdaptiveLevelsIndexAllocator.NewItem() ;
00232 
00233        // new index (l,t)
00234        nali->dim=DIM-1 ;
00235        for (i=0;i <DIM-1; i++) { nali->l[i]=ali->l[i], nali->t[i]=ali->t[i] ;} 
00236 
00237        // write new IndexSet
00238        for (ld=Level0[dir]; ld<=LMAX; ld++) {
00239          l[dir]=ld ;
00240          IndexSetWriteToAdaptiveL((*nals)[ld], &CounterArray[l], AB2[0]->L[ld].IS) ;
00241        }
00242    
00243        // insert 
00244        pair<AdaptiveLevelsIndex* , AdaptiveLevels*> p(nali,nals) ;
00245        (*npg).insert(p) ;
00246 
00247        // set 'Active' flags in AdaptiveData associated to the OLD grid
00248        // ActiveEntries is refined to be associated to the NEW grid later !!! 
00249        if ( (ActiveEntries != NULL) && (dir==0) ) {
00250          // reload old Index Set (Flags are alreday set by 'AB2->Refine(...)') and then write to ActiveEntries
00251          for (ld=Level0[0]; ld<=LMAX; ld++) {
00252            l[0]=ld ;
00253            AB1[0]->L[ld].IS->Load ((*als)[ld]) ;
00254            AB1[0]->L[ld].IS->Store((*als)[ld], ActiveEntries->ba->d[l] , AB1[0]->L[ld].d ) ;
00255          } // next ld
00256        }
00257      }
00258  
00259    }     // next pgit 
00260 
00261    if ((R->strategy==AdaptivityCriterion::SPARSE_GRID)||(R->strategy==AdaptivityCriterion::FIXED_LEVEL)) {
00262      BasisPG[0]=npg ;
00263      OldBasisPG=pg  ;
00264      GenerateAllBasisPGs() ;
00265      goto ende ;
00266    }
00267 
00268    if (dir<DIM-1) GeneratePG(npg,dir,BasisPG[dir+1],dir+1) ;
00269      
00270  } // next dir
00271 
00272  // finish: OldBasisPG -> BasisPG[D-1] , BasisPG[0]->OldBasisPG , BasisPG[D-1] -> BasisPG[0]   
00273  pg            =BasisPG[DIM-1] ; 
00274  BasisPG[DIM-1]=OldBasisPG ;
00275  OldBasisPG    =BasisPG[0] ;
00276  BasisPG[0]    =pg ;
00277 
00278  sall=Size() ;
00279  for (dir=0; dir<DIM-1; dir++) {
00280    noe=GeneratePG(BasisPG[DIM-1], DIM-1, BasisPG[dir], dir) ;
00281    if (noe!=sall) std::cout << "Error in RefineGrid: now="<<noe<<" , Size="<<sall<<'\n' ;
00282  }
00283 
00284 ende: ;
00285 
00286  was_refined=true ;
00287  bool flag=false  ;
00288  TmpBasis->Resize(&CounterArray) ;
00289  for (i=0; i<ResizeGroup.count; i++) 
00290    if (ActiveEntries==ResizeGroup.Ds[i]) flag=true ;
00291                                    else ResizeGroup.Ds[i]->ba->Resize(&CounterArray) ;
00292 
00293  for (i=0; i<RefineGroup.count; i++) RefineGroup.Ds[i]->Refine() ;
00294 
00295  // Refine ActiveEntries even if it is in ResizeGroup only
00296  if (flag) ActiveEntries->Refine() ;
00297 
00298  R->MaxLevel=MxLevel ;
00299 
00300  SMPLoadBalance() ;
00301 
00302  // FaceGrids
00303  for (i=0;i<DIM; i++) {
00304    if (FaceGrid[i][0]) FaceGrid[i][0]->SetAsSlice(this,i,Level0[i],0) ;
00305    if (FaceGrid[i][1]) FaceGrid[i][1]->SetAsSlice(this,i,Level0[i],1<<Level0[i]) ;
00306  }
00307 
00308 }
00309 
00310 //
00311 // Refine copies the Data with respect to the old basis to another data-structure with respect to the
00312 // new basis.
00313 //
00314 template<int DIM>
00315 void AdaptiveData<DIM>::Resize() {
00316   ba->Resize(&G->CounterArray) ;
00317 }
00318 
00319 template<int DIM>
00320 void AdaptiveData<DIM>::Refine() {
00321 
00322   int    i,l0,l[DIM],ni ;
00323   size_t nj=MAXINT,j=MAXINT ;
00324 
00325   AdaptiveGrid<DIM>::ProjectionGrid             *pg  ,*npg  ;
00326   AdaptiveGrid<DIM>::ProjectionGrid::iterator    pgit,npgit ;
00327 
00328   AdaptiveLevels::AdaptiveL             *al ,*nal    ;
00329   AdaptiveLevels::AdaptiveL::iterator    alit ,nalit ;
00330   AdaptiveDataArray<DIM>::DataVector            *v,*nv ;
00331 
00332    pg=G->OldBasisPG ;
00333   npg=G->BasisPG[0] ;
00334 
00335   //check whether Grid has been refined, such that OldBasisPG is a valid structure
00336   assert(G->was_refined) ;
00337 
00338   // initialize with zero, zero fill is nesseccary
00339   G->TmpBasis->Resize(&G->CounterArray,0.0) ;
00340 
00341   //
00342   // Main loop 
00343   //
00344   for (pgit=(*pg).begin(); pgit != (*pg).end(); ++pgit) {
00345     AdaptiveLevelsIndex *ali=(*pgit).first ;
00346     AdaptiveLevels      *als=(*pgit).second,*nals ;
00347 
00348     // find the (new) 1D AdaptiveLevels structure which corresponds to the one of the old PG 
00349     npgit=(*npg).find(ali) ;
00350 
00351     // copy 1D data
00352     if (npgit!=(*npg).end()) {
00353       // prepare levels vector
00354       for (i=1; i<DIM; i++) l[i]=ali->l[i-1] ;
00355       nals=(*npgit).second ;
00356       
00357       // forall 1D Levels
00358       for (l0=G->Level0[0]; l0<= LMAX; l0++) {
00359         l[0] =l0 ;
00360         v    =ba->d[l]       ;
00361         nv   =G->TmpBasis->d[l] ;
00362         
00363         al  =(* als)[l0] ;
00364         nal =(*nals)[l0] ;
00365         
00366         if ((*al).size()>0) {
00367           alit=(*al).begin() ;
00368           i=-1 ;
00369           for (nalit=(*nal).begin(); nalit!=(*nal).end() ; nalit++) { 
00370             ni=(*nalit).first  ;
00371             nj=(*nalit).second ;
00372             while ((alit!=(*al).end()) && (i<ni)) { 
00373               i=(*alit).first  ;
00374               j=(*alit).second ;
00375               alit++ ;
00376             }
00377             if (i==ni) (*nv)[nj]=(*v)[j] ;
00378             // if (i> ni) (*nv)[nj]=0.0     ;
00379           }
00380         }
00381       }   // l0 
00382     }    // if ()     
00383   }      // pgit  
00384   
00385   //swap basis 
00386   AdaptiveDataArray<DIM> *t=ba ;
00387   ba=G->TmpBasis ;
00388   G->TmpBasis =t ;
00389  
00390   G->TmpBasis->Resize(&G->CounterArray) ;
00391 
00392 }
00393 
00394 template<int DIM>
00395 void AdaptiveGrid<DIM>::ConsistencyCheck() {
00396   int dir ;
00397 }
00398 
00399 
00400 
00401 template<int DIM>
00402 size_t AdaptiveGrid<DIM>::Size() {
00403   Matrix<size_t,DIM>::iterator  l(&CounterArray) ;
00404 
00405   size_t s=0 ;
00406   for (l=CounterArray.begin(); l<=CounterArray.end(); ++l) s += CounterArray[l] ;
00407   return s ;
00408 }
00409 
00410 
00411 template<int DIM>
00412 void AdaptiveGrid<DIM>::PrintCounterArray() {
00413   int i,n=(int)(log10(CounterArray.MaxAbs()+1)+1.5) ;
00414   Matrix<size_t,DIM>::iterator  l(&CounterArray, DIM-1 ) ;
00415   for (l=CounterArray.begin(DIM-1); l<=CounterArray.end(); ++l) {
00416     std::cout<< "CA(" ;
00417     for (i=0; i<DIM-1; i++) std::cout<<setw(2)<<l.i[i]<<',' ;
00418     std::cout<<".) " ;
00419     for (l.i[DIM-1]=CounterArray.b[DIM-1]; l.i[DIM-1]<=CounterArray.e[DIM-1]; l.i[DIM-1]++)
00420       std::cout<<setw(n)<<CounterArray[l]<<' ' ;
00421 
00422     std::cout<<'\n' ;
00423     l.i[DIM-1]=CounterArray.b[DIM-1];
00424   }
00425 }
00426 
00427 template<int DIM>
00428 void AdaptiveGrid<DIM>::PrintLargestActiveLevel() {
00429   int J[DIM],i ;
00430   LargestActiveLevel(J) ;
00431   std::cout<<"largest active levels: " ;
00432   for (i=0; i<DIM; i++) std::cout<<J[i]<<"  ";
00433   std::cout<<'\n' ;
00434 }
00435 
00436 
00437 template<int DIM>
00438 void AdaptiveGrid<DIM>::DataMemory(int *reserved , int *used) {
00439   int i,r,u ;
00440   *reserved=*used=0 ;
00441   for (i=0; i<RefineGroup.count; i++) {
00442     RefineGroup.Ds[i]->ba->Memory(&r,&u) ;
00443     (*reserved)+=r ;
00444     (*used    )+=u ;
00445   }
00446   for (i=0; i<ResizeGroup.count; i++) {
00447     ResizeGroup.Ds[i]->ba->Memory(&r,&u) ;
00448     (*reserved)+=r ;
00449     (*used    )+=u ;
00450   }
00451 }
Generated at Mon Aug 19 10:02:31 2002 for AWFD by