AdaptiveGrid.hpp

00001 #ifndef _ADAPTIVE_GRID_
00002 # define _ADAPTIVE_GRID_
00003 
00004 # include<assert.h>
00005 # include<vector>
00006 # include<map>
00007 # include<pthread.h>
00008 
00009 # include "Wavelet.hpp"
00010 # include "Adaptive1D.hpp"
00011 # include "NonAdaptive.hpp"
00012 # include "UDF.h"
00013 # include "Refine.hpp"
00014 # include "Buffer.hpp"
00015 
00016 # include "Matrix.hpp"
00017 # include "AdaptiveLevels.hpp"
00018 # include "BlockAllocator.hpp"
00019 # include "DataGroup.hpp"
00020 
00021 BlockAllocator<AdaptiveLevels     > AdaptiveLevelsAllocator      ;
00022 BlockAllocator<AdaptiveLevelsIndex> AdaptiveLevelsIndexAllocator ;
00023 
00024 //
00025 //  THE GRID 
00026 //
00027 # define _MAX_ADAPTIVELEVELPOINTERS_ 100000
00028 
00043 template<int DIM>
00044 struct AdaptiveGrid {
00045  
00047   typedef map<AdaptiveLevelsIndex* , AdaptiveLevels* , AdaptiveLevelsIndexCmp > ProjectionGrid ;
00048 
00049   //public:
00051          AdaptiveGrid() ;
00053          AdaptiveGrid(Wavelets *W,bool Init_tmp = true) ;
00054   //  see Init(.,.,...)
00055          AdaptiveGrid(Wavelets *W,bool Init_tmp, 
00056                       AdaptiveB **B1, AdaptiveB **B2, AdaptiveGS **GS1, AdaptiveGS **GS2, SMPjob *SMP) ;
00058         ~AdaptiveGrid() ;
00061   void   Init(Wavelets *W,bool Init_tmp=true) ;
00062   //  Init-function with extended functionality, if B1,...,SMP are not NULL (all of them must be !=NULL then)
00063   //  the resources AB1,AB2 and in particular the threads are not allocated/invoked but shared
00064   //  This is necessary for initializing lower-dimensional subgrids ...  
00065   void   Init(Wavelets *W,bool Init_tmp, 
00066               AdaptiveB **B1, AdaptiveB **B2, AdaptiveGS **GS1, AdaptiveGS **GS2, SMPjob *SMP) ;
00067 
00070   void   SetPeriodicConditions(int  BC[DIM][2]) ;
00071   void   Clear() ;
00072 
00075   void   Set(Function *F, double c) ;
00077   void   SetSparse(int sl) ;
00078   void   SetSparse(int sl,int *MaxLevel) ;
00080   void   SetFull  (int *L) ;
00082   void   SetFull  (int  L) ;
00083   
00087   void   SetAsSlice(AdaptiveGrid<DIM+1> *G, int dir, int l, int t) ;
00088 
00102   void   Refine(AdaptiveData<DIM> *TestData , AdaptivityCriterion *R , AdaptiveData<DIM> *ActiveEntries = NULL) ;
00103   
00104   // check wether RefineGrid failed to produce consistent pointers
00105   void   ConsistencyCheck() ;
00106 
00108   size_t Size() ;
00110   void   PrintCounterArray()    ;
00112   void   DataMemory(int *reserved , int *used) ; // return the size of memory (in bytes) allocated for all attached Data<DIM>s.
00114   void   LargestActiveLevel(int *L) ;
00116   void   PrintLargestActiveLevel() ;
00117 
00118   //
00119   AdaptiveLevels::AdaptiveL* FindOrInsertAdaptiveL(ProjectionGrid *PG,int dirTo, int *l,int *t) ;
00122   long  *FindOrInsert(ProjectionGrid *PG,int dirTo, int *l,int *t) ;
00125   long  Insert       (ProjectionGrid *PG,int dirTo, int *l,int *t) ;
00128   void   Insert      (ProjectionGrid *PG,int dirTo, int *l,int *t , long  vectorindex) ;
00131   long   Find        (ProjectionGrid *PG,int dirTo, int *l,int *t) ;
00132 
00133   // Get the numerical value for index (l,t) and the given AdaptiveDataArray. <br>
00134   //double Read        (ProjectionGrid *PG,int dirTo, int *l,int *t , AdaptiveDataArray<DIM> *Dat) ;
00135   // If (l,t) is in the adaptive index set associated to *PG, then the corresponding component of vector (*Dat)[l] := x <br>
00136   // Otherwise insert (l,t) in the adaptive index set associated to *PG, make sure that the vector (*Dat)[l] is large enough
00137   // to store the new entry and set this entry to x, adjust the sizes (CounterArray[l])  
00138   //void   Write       (ProjectionGrid *PG,int dTo, int *l,int *t , AdaptiveDataArray<DIM> *Dat,double x) ;
00139 
00141   size_t GeneratePG         (ProjectionGrid *PGFrom,int dirFrom , ProjectionGrid *PGTo  ,int dirTo) ;
00142   void   GenerateAllBasisPGs() ;
00143   void   DumpPG             (ProjectionGrid *PG, int dir) ;
00144   void   ClearPG            (ProjectionGrid *PG, int dir=-1) ; 
00145 
00147   void   SMPLoadBalance() ;
00148 
00149   //void   WriteMatLab      (int *J,const char *name, Wavelets *W) ;
00150   void   WriteSparse(const char *name) ;
00153   void   ReadSparse (const char *name) ;
00154  
00155   //internal data:
00156   int    Level0  [DIM]   ; //default: 0
00157   int    MaxLevel[DIM]   ; //default: {LMAX-1,...}
00158   Wavelets        *WC    ;
00159   bool   IsPeriodic[DIM] ;
00160   bool   IsPeriodicValid ;
00161   double XA[DIM],XE[DIM] ; // box associated to the grid: [XA_0,XE_0] x ... x [XA_.,XE_.], default XA=0.0 , XE=1.0 ;
00162  
00163   // if *this is just a slice of a (DIM+1)-dimensional AdaptiveGrid we store some information about the
00164   // relationship
00165   AdaptiveGrid<DIM+1> *SliceParent ;
00166   int                  SliceDir    ;
00167   int                  Slicel      ;
00168   int                  Slicet      ;
00169 
00170   // ProjectionGrids: each is a  DIM-1 dimensional Array of HashTables of AdaptiveLevels
00171   ProjectionGrid PGS[DIM+1] ;
00172   ProjectionGrid *BasisPG[DIM] ; // (*BasisPG[i]) allows for the access to lines of the Basis along the i-th coordinate
00173   ProjectionGrid *OldBasisPG ;
00174   bool was_refined ; 
00175 
00177   AdaptiveGrid<DIM-1> *FaceGrid[DIM][2] ;
00178 
00179   enum {NUM_TMP=13} ;
00181   AdaptiveData<DIM>   *tmp[NUM_TMP] ;
00182 
00183   AdaptiveDataArray<DIM> *TmpBasis  ; // memory for new data in AdaptiveData::Refine() ;
00184   Matrix<size_t,DIM>  CounterArray ;  // numbers of active coefficients for each level 
00185   DataGroup<DIM>      RefineGroup , ResizeGroup ;
00186 
00187   //
00188   // SMP stuff
00189   //
00190   // Alp[DIM] will contain pointers to all the AdaptiveLevels which have to  
00191   // be processed in ApplyOp(..,dir,...)
00192   // smpJobs[dir][_SMP_PROC_] is a set of increasing numbers with smpJobs[dir][0]=0 
00193   // such that thread i will process all AdaptiveLevels numbered in Alp[dir]
00194   //
00195   AdaptiveLevelsPointer *smpalp[DIM], *smpalptmp,*smpalpstart[DIM][_SMP_PROC_]  ;
00196   size_t                 smpnum_alp[DIM][_SMP_PROC_] ;
00197 
00198   SMPjob                 *SMPjobs ; // data structure with all the information
00199 
00200   pthread_t              smpthreads [_SMP_PROC_] ;
00201   size_t                 smpjobcount ; // 
00202 
00203   // 'accumulators'  for the univariate operator evaluation
00204   AdaptiveB             **AB1, **AB2; 
00205   AdaptiveGS            **AGS, **AG2;
00206   bool AB1owner ;
00207 
00208   void                   SMPWaitForAllThreads() ; // waits for all threads to be finished 
00209 
00210   /*
00216   struct iterator {
00218     iterator() {;}
00220     void Init(AdaptiveGrid<DIM> *Gr) ;
00222     bool isend() ;
00224     iterator& operator++() ;
00226     int l[DIM] ;
00228     int t[DIM] ;
00230     size_t j ;
00232     AdaptiveLevels *als ;
00234     AdaptiveL      *al  ;
00235   private:
00236     ProjectionGrid::iterator pgit   ;
00237     AdaptiveL::iterator      alit   ;
00238     AdaptiveGrid<DIM>*       G      ;
00239   } ;
00240 
00241 
00244   struct boundaryiterator {
00246     boundaryiterator() ;
00248     void Init(AdaptiveGrid<DIM> *G) ;
00250     bool isend() ;
00252     boundaryiterator& operator++() ;
00254     int l[DIM] ;
00256     int t[DIM] ;
00258     size_t j ;
00260     AdaptiveLevels *als ;
00262     AdaptiveLevelsIndex *ali ;
00264     AdaptiveL      *al  ;
00265   private:
00266     ProjectionGrid::iterator pgit ;
00267     AdaptiveL::iterator      alit ;
00268     bool  ende ;
00269   } ;
00270   */
00271 
00272 } ;
00273 
00274 
00275 
00276 
00277 // Specialization for 0-dimensional AdaptiveGrids which contain exactly one point
00278 // Here, just dummy functions with the right interfaces are needed
00279 template<>
00280 struct AdaptiveGrid<0> {
00281          AdaptiveGrid() {Init(NULL,false);}
00282          AdaptiveGrid(Wavelets *  ,bool ) {Init(NULL,false);}
00283          AdaptiveGrid(Wavelets *  ,bool , 
00284                       AdaptiveB **, AdaptiveB **, AdaptiveGS **, AdaptiveGS **, SMPjob *) {Init(NULL,false);}
00285         ~AdaptiveGrid() { for (int i=0; i<100; i++)  delete tmp[i]; }
00286   void   Init(Wavelets *,bool Init_tmp=true) {Init(NULL,Init_tmp,NULL,NULL,NULL,NULL,NULL);}
00287   void   Init(Wavelets *,bool Init_tmp, AdaptiveB **, AdaptiveB **, AdaptiveGS **, AdaptiveGS **, SMPjob *) ;
00288   void   SetAsSlice(AdaptiveGrid<1> *G, int dir, int l, int t) {SliceParent=G; SliceDir=dir; Slicel=l; Slicet=t; }
00289 
00290   AdaptiveGrid<1> *SliceParent ;
00291   int              SliceDir    ;
00292   int              Slicel      ;
00293   int              Slicet      ;
00294 
00295   AdaptiveData<0> *tmp[100] ;
00296 } ;
00297 
00298 
00300 //                   //
00301 //  Implementations  //
00302 //                   //
00304 
00305 /*
00307 // iterators //
00309 template<int DIM>
00310 void AdaptiveGrid<DIM>::Init(AdaptiveGrid<DIM> *Gr) {
00311   G=Gr ;
00312   pgit  =(*G->BasisPG[0]).begin() ;
00313   l[0]  =  G->Level0[0] ;
00314   if (pgit != (*G->BasisPG[0]).end()) {
00315     ali=(*pgit).first  ;
00316     als=(*pgit).second ;
00317     for (int i=1; i<DIM; i++) { l[i]=ali->l[i-1], t[i]=ali->t[i-1]; }
00318     
00319     al  =(*als)[l[0]] ;
00320     alit=(*al).begin();
00321     if (alit != (*al).end()) {
00322       t[0]=(*alit).first  ;
00323       j   =(*alit).second ;
00324     } else ende=true ;
00325 
00326   } else ende=true ;
00327 }
00328 
00329 template<int DIM>
00330 void AdaptiveGrid<DIM>::isend() { return ende ; }
00331 
00332 template<class I,int DIM>
00333 AdaptiveGrid<DIM>::iterator& AdaptiveGrid<DIM>::iterator::operator++() {
00334   alit++ ;
00335  setalit:; // new spatial index
00336   if (alit != (*al).end()) {
00337     t[0]=(*alit).first  ;
00338     j   =(*alit).second ;
00339   } else { // new level
00340 incl0: ;
00341     l[0]++ ;
00342     if (l[0]<=LMAX) {
00343       al  =(*als)[l[0]] ;
00344       alit=(*al).begin();
00345       goto setalit ;
00346     } else { // new AdaptiveLevels
00347       pgit++ ;
00348       if (pgit!=(*G->BasisPG[0]).end()) {
00349         ali=(*pgit).first  ;
00350         als=(*pgit).second ;   
00351         for (int i=1; i<DIM; i++) { l[i]=ali->l[i-1], t[i]=ali->t[i-1]; }
00352         
00353         l[0]=G->Level0[0] ;
00354         al  =(*als)[l[0]] ;
00355         alit=(*al).begin();
00356         goto setalit ;
00357       } else { ende=true ; // stop }
00358     }
00359   }
00360   return *this ;
00361 } 
00362 */
00363 
00364 void AdaptiveGrid<0>::Init(Wavelets *,bool , AdaptiveB **, AdaptiveB **, AdaptiveGS **, AdaptiveGS **, SMPjob *) {
00365   SliceParent=NULL ;
00366   for (int i=0; i<100; i++) tmp[i]=new AdaptiveData<0>(this) ;
00367 }
00368 
00369 template<int DIM>
00370 AdaptiveGrid<DIM>::AdaptiveGrid() {
00371   size_t i ;
00372 
00373   IsPeriodicValid=false ;
00374   for (i=0; i<DIM; i++) Level0[i]=0 ;
00375   Init((Wavelets*)NULL) ;
00376 }
00377 
00378 template<int DIM>
00379 AdaptiveGrid<DIM>::AdaptiveGrid(Wavelets *WC,bool InitData) { 
00380   size_t i ;
00381   IsPeriodicValid=false ;
00382   for (i=0; i<DIM; i++) Level0[i]=0 ;
00383   Init(WC,InitData) ;
00384 }
00385 
00386 template<int DIM>
00387 AdaptiveGrid<DIM>::AdaptiveGrid(Wavelets *WC,bool InitData,
00388                                 AdaptiveB **B1, AdaptiveB **B2, AdaptiveGS **GS1, AdaptiveGS **GS2, SMPjob *SMP) { 
00389   size_t i ;
00390   IsPeriodicValid=false ;
00391   for (i=0; i<DIM; i++) Level0[i]=0 ;
00392   Init(WC,InitData,B1,B2,GS1,GS2,SMP) ;
00393 }
00394 
00395 template<int DIM>
00396 void AdaptiveGrid<DIM>::Init(Wavelets *W,bool InitData) {
00397   Init(W,InitData,NULL,NULL,NULL,NULL,NULL) ;
00398 }
00399 
00400 template<int DIM>
00401 void AdaptiveGrid<DIM>::Init(Wavelets *W,bool InitData,
00402                              AdaptiveB **B1, AdaptiveB **B2, AdaptiveGS **GS1, AdaptiveGS **GS2, SMPjob *SMP) 
00403 { 
00404   int k,i,j,err ;
00405 
00406   if (debugRefine) std::cout<<"Grid::Init _SMP_PROC_="<<_SMP_PROC_<<" , LMAX="<<LMAX<<'\n' ;
00407 
00408   if (Level0[0]!=0) {
00409     std::cout<< "Grid initialized before\n Whats going on here ?\n" ;
00410     exit(-1) ;
00411   }
00412 
00413   //security check for AdaptiveLevelsPointer
00414   if (DMAX < DIM) {
00415     std::cout<< "Instantiation of AdaptiveGrid with dimension "<<DIM<<" larger than DMAX="<<DMAX<<" not possible!\n See/edit $(HOME)/Makefile\n" ;
00416     exit(-1) ;
00417   }
00418 
00419   // set simple values
00420   for (k=0; k<DIM  ; k++) { 
00421     Level0    [k]= (W!=NULL) ? W->Level0 : 0 ;
00422     MaxLevel  [k]= LMAX-1 ;
00423     IsPeriodic[k]= false  ;
00424     XA[k]=0.0 ;
00425     XE[k]=1.0 ;
00426    }
00427 
00428   WC=W ;
00429   was_refined=false ;
00430 
00431   // Distribution of PGS
00432   for (i=0;i<DIM;i++) BasisPG[i]=&PGS[i] ;
00433   OldBasisPG=&PGS[DIM] ;
00434 
00435   // set Grid for tmp-Data
00436   for (i=0; i<NUM_TMP; i++) {
00437     if (InitData) tmp[i]=new AdaptiveData<DIM>(this) ;
00438              else tmp[i]=NULL ;
00439    }
00440 
00441   // allocate temp. memory for refine
00442   TmpBasis=new AdaptiveDataArray<DIM> ;
00443 
00444   // init. CounterArray
00445   int b[DIM],e[DIM] ;
00446   for (i=0; i<DIM; i++) b[i]=0 , e[i]=LMAX ;
00447   CounterArray.Init(&b[0],&e[0]) ;
00448   CounterArray.Set((size_t)0) ;
00449 
00450   //
00451   // allocation of the 1D AdaptiveBasis/GS Datastructures
00452   //
00453   if (B1!=NULL) {
00454 
00455     //check for consistency
00456     assert(B2) ;
00457     assert(GS1) ;
00458     assert(GS2) ;
00459     assert(SMP) ;
00460  
00461     // set resources
00462     AB1owner=false ;
00463     AB1=B1 ;
00464     AB2=B2 ;
00465     AGS=GS1 ;
00466     AG2=GS2 ;
00467     SMPjobs=SMP ;
00468   } else {
00469     AB1owner=true ;
00470     AB1=new AdaptiveB* [_SMP_PROC_] ;
00471     AB2=new AdaptiveB* [_SMP_PROC_] ;
00472     AGS=new AdaptiveGS*[_SMP_PROC_] ;
00473     AG2=new AdaptiveGS*[_SMP_PROC_] ;
00474     SMPjobs=new SMPjob [_SMP_PROC_] ;
00475 
00476     for (i=0; i<_SMP_PROC_; i++) {
00477       err=0 ;
00478       AB1[i]=new AdaptiveB() ;
00479       if (AB1[i]->Init(Level0[0],LMAX)==0) err |= 1; 
00480       for (j=0; j<2; j++) AB1[i]->Buffers[j]=new doubleBuffer(LMAX) ;
00481 
00482       AB2[i]=new AdaptiveB(AB1[i]) ; // AB2 and AB1 share the same IndexSets
00483       AGS[i]=new AdaptiveGS() ;
00484       AG2[i]=new AdaptiveGS() ;
00485     
00486       if (AGS[i]->Init(Level0[0],LMAX)==0) err |= 2 ;
00487       if (AG2[i]->Init(Level0[0],LMAX)==0) err |= 4 ;
00488 
00489       for (j=0; j<2; j++) { 
00490         AGS[i]->Buffers[j]=AB1[i]->Buffers[j] ;
00491         AG2[i]->Buffers[j]=AB1[i]->Buffers[j] ;
00492       }
00493 
00494       if (err) { 
00495         std::cout << "AdaptiveGrid<"<<DIM<<">::Init(..) Could allocate AdaptiveB/GS err="<<err<<'\n' ;
00496         exit(-1) ;
00497       }
00498  
00499       SMPjobs[i].AB1 = AB1[i] ;
00500       SMPjobs[i].AB2 = AB2[i] ;
00501       SMPjobs[i].AGS = AGS[i] ;
00502       SMPjobs[i].AG2 = AG2[i] ;
00503       SMPjobs[i].WC  = WC     ;
00504       SMPjobs[i].DIM = MAXINT ;
00505       SMPjobs[i].dir = MAXINT ;      
00506  
00507 #if _SMP_PROC_ > 1
00508       if (i>0) {
00509         err=pthread_create(&smpthreads[i], NULL , PollProcessSMPjob , (void*) &SMPjobs[i] );
00510         if (err) { std::cout << "failed start thread "<<i<<" error "<< err<<'\n' ; exit(-1) ; }
00511       }
00512 #endif
00513     }
00514 
00515   }
00516  
00517   // init the linear storage for SMP parallelization
00518   err=0 ;
00519   for (i=0; i<DIM; i++)
00520     if ((smpalp[i]=(AdaptiveLevelsPointer *)calloc(_MAX_ADAPTIVELEVELPOINTERS_,
00521                                                    sizeof(AdaptiveLevelsPointer))) == NULL) err |= 1 ;
00522   
00523   if ((smpalptmp=(AdaptiveLevelsPointer *)calloc(_MAX_ADAPTIVELEVELPOINTERS_,
00524                                                  sizeof(AdaptiveLevelsPointer))) == NULL) err |= 2 ;
00525   
00526   if (err) {
00527     std::cout << "AdaptiveGrid<" <<DIM<< ">::Init  Could not allocate smpalp/tmp "<<err<<'\n' ;
00528     exit(-1) ;
00529   }
00530 
00531   smpjobcount=0 ;
00532 
00533   //Slice stuff  
00534   SliceParent=NULL ;
00535   SliceDir   =-1   ;
00536   Slicel     =0    ;
00537   Slicet     =-1   ;
00538 
00539  // grids for faces
00540   for (i=0; i<DIM; i++) {
00541     FaceGrid[i][0]=new AdaptiveGrid<DIM-1>(W,InitData,AB1,AB2,AGS,AG2,SMPjobs) ;
00542     FaceGrid[i][1]=new AdaptiveGrid<DIM-1>(W,InitData,AB1,AB2,AGS,AG2,SMPjobs) ;
00543   }
00544 
00545 }
00546 
00547 
00548 template<int DIM>
00549 AdaptiveGrid<DIM>::~AdaptiveGrid() {
00550   int  d,i ;
00551 
00552   // remove FaceGrids 
00553   for (i=0; i<DIM; i++) { 
00554     delete FaceGrid[i][0] ;
00555     delete FaceGrid[i][1] ;
00556   }
00557 
00558   // clear some of the resources
00559   for (i=0; i<DIM; i++) free(smpalp[i]) ;
00560                         free(smpalptmp) ;
00561 
00562   // removal of all AdaptiveLevels
00563   for (d=0; d<DIM+1; d++) ClearPG(&PGS[d]) ;
00564 
00565   for (d=0;d<NUM_TMP;d++) delete tmp[d] ;
00566   delete TmpBasis ;
00567   
00568   // clear all Data associated to (*this)
00569   for (i=0; i<ResizeGroup.count+RefineGroup.count; i++) {
00570     AdaptiveData<DIM> *X ;
00571     if (i < ResizeGroup.count) X=ResizeGroup.Ds[i] ;
00572                           else X=RefineGroup.Ds[i-ResizeGroup.count] ;
00573     delete X->ba ;
00574     X->ba=NULL   ;
00575     X->G =NULL   ;
00576   }
00577 
00578   // clear shared resources
00579   if (AB1owner) {
00580     // send final signal to threads
00581     for (i=1; i<_SMP_PROC_; i++) { 
00582       SMPjobs[i].terminate=true ; 
00583       SMPjobs[i].clrpoll() ; 
00584     }
00585   
00586     for (size_t i=0; i<_SMP_PROC_; i++) {
00587       if (AB1[i]) { 
00588         for (int j=0; j<2; j++) delete AB1[i]->Buffers[j]; 
00589         delete AB1[i] ; 
00590       }
00591       if (AB2[i]) delete AB2[i] ;
00592       if (AGS[i]) delete AGS[i] ;
00593       if (AG2[i]) delete AG2[i] ;
00594     }
00595     delete []AB1 ;
00596     delete []AB2 ;
00597     delete []AGS ;
00598     delete []AG2 ; 
00599     delete []SMPjobs ;
00600   }
00601 
00602 }
00603 
00604 
00605 template<int DIM>
00606 void AdaptiveGrid<DIM>::Clear() {
00607   for (int d=0; d<DIM+1; d++) ClearPG(&PGS[d]) ;
00608   CounterArray.Set((size_t)0) ;
00609 }
00610 
00611 
00612 template<int DIM>
00613 void AdaptiveGrid<DIM>::SetPeriodicConditions(int BC[DIM][2]) {
00614   int i ;
00615   if (IsPeriodicValid) {
00616     for (i=0; i<DIM; i++) 
00617       if (IsPeriodic[i] != (BC[i][1]==PERIODIC)) { 
00618         std::cout << "AdaptiveGrid::SetPeriodicConditions: you attach AdaptiveDatas with different periodic directions to the same AdaptiveGrid\n" ; 
00619         assert(0) ;
00620        }
00621   } else {
00622     for (i=0;i<DIM;i++) IsPeriodic[i]=(BC[i][1]==PERIODIC) ;
00623     IsPeriodicValid=true ;
00624   }
00625 }
00626 
00627 
00628 template<int DIM>
00629 long AdaptiveGrid<DIM>::Find(ProjectionGrid *PG,int dir, int *l,int *t) {
00630   int  i ;
00631   long r = -MAXINT ;
00632 
00633   // generate key
00634   AdaptiveLevelsIndex *lt=AdaptiveLevelsIndexAllocator.NewItem() ;
00635   lt->dim=DIM-1 ;
00636   for (i=0    ; i<dir; i++) lt->l[i  ]=l[i], lt->t[i  ]=t[i];
00637   for (i=dir+1; i<DIM; i++) lt->l[i-1]=l[i], lt->t[i-1]=t[i];
00638  
00639   // find key
00640   ProjectionGrid::iterator  plt=(*PG).find(lt) ;
00641   
00642   if (plt!=(*PG).end()) {
00643     AdaptiveLevels::AdaptiveL           *al =(*(*plt).second)[l[dir]] ; 
00644     AdaptiveLevels::AdaptiveL::iterator  ait=(*al).find(t[dir]) ;
00645     if (ait!=(*al).end()) r=(*ait).second ;
00646    }
00647 
00648   // remove auxiliary key
00649   AdaptiveLevelsIndexAllocator.DeleteItem(lt) ;
00650   return r ;
00651 }
00652 
00653 
00654 template<int DIM>
00655 AdaptiveLevels::AdaptiveL* AdaptiveGrid<DIM>::FindOrInsertAdaptiveL(ProjectionGrid *PG,int dir, int *l,int *t) {
00656   int i ;
00657   bool ins=false ;
00658 
00659   // generate key
00660   AdaptiveLevelsIndex *lt=AdaptiveLevelsIndexAllocator.NewItem() ;
00661   lt->dim=DIM-1 ;
00662   for (i=0    ; i<dir; i++) lt->l[i  ]=l[i], lt->t[i  ]=t[i] ;
00663   for (i=dir+1; i<DIM; i++) lt->l[i-1]=l[i], lt->t[i-1]=t[i] ;
00664 
00665   // find key
00666   ProjectionGrid::iterator  plt=(*PG).find(lt) ;
00667   AdaptiveLevels           *als ;
00668 
00669   // insert if necessary
00670   if (plt==(*PG).end()) {
00671     ins=true ;
00672     als=AdaptiveLevelsAllocator.NewItem() ;
00673     pair<AdaptiveLevelsIndex* , AdaptiveLevels*> p(lt , als) ;
00674     (*PG).insert(p) ;
00675   } else
00676     als=(*plt).second ; 
00677 
00678   // remove auxiliary key
00679   if (!ins) AdaptiveLevelsIndexAllocator.DeleteItem(lt) ;
00680 
00681   return (*als)[l[dir]] ;
00682 }
00683 
00684 template<int DIM>
00685 long AdaptiveGrid<DIM>::Insert(ProjectionGrid *PG,int dir, int *l,int *t) {
00686   long vi;
00687   AdaptiveLevels::AdaptiveL           *al =FindOrInsertAdaptiveL(PG,dir,l,t) ;
00688   AdaptiveLevels::AdaptiveL::iterator alit=(*al).find(t[dir]) ;
00689 
00690   if (alit==(*al).end()) vi=(*al)[t[dir]]=(long)(CounterArray[l]++) ;
00691   else                   vi=(*alit).second ;
00692   return vi;
00693 }
00694 
00695 template<int DIM>
00696 void AdaptiveGrid<DIM>::Insert(ProjectionGrid *PG,int dir, int *l,int *t , long second) {
00697   AdaptiveLevels::AdaptiveL *al = FindOrInsertAdaptiveL(PG,dir,l,t) ;
00698   (*al)[t[dir]]=second ;
00699 }
00700 
00701 
00702 template<int DIM>
00703 long *AdaptiveGrid<DIM>::FindOrInsert(ProjectionGrid *PG,int dir, int *l,int *t) {
00704   AdaptiveLevels::AdaptiveL           *al = FindOrInsertAdaptiveL(PG,dir,l,t) ;
00705   AdaptiveLevels::AdaptiveL::iterator alit= (*al).find(t[dir]) ;
00706 
00707   if (alit==(*al).end()) {
00708     (*al)[t[dir]]=MAXINT ;
00709     alit=(*al).find(t[dir]) ; 
00710   }
00711   return &((*alit).second) ;
00712 }
00713 
00714 
00715 /*
00716 template<int DIM>
00717 double AdaptiveGrid<DIM>::Read(AdaptiveGrid<DIM>::ProjectionGrid *PG,int dTo, int *l,int *t , AdaptiveDataArray<DIM> *Dat) {
00718   long *i=FindOrInsert(PG,dTo,l,t) ;
00719   assert(*i!=MAXINT) ;
00720   return (*(*Dat).d[l])[*i] ;
00721 }
00722 
00723 template<int DIM>
00724 void AdaptiveGrid<DIM>::Write(AdaptiveGrid<DIM>::ProjectionGrid *PG,int dTo, int *l,int *t , AdaptiveDataArray<DIM> *Dat,double x) {
00725   long *i=FindOrInsert(PG,dTo,l,t) ;
00726 
00727   AdaptiveDataArray<DIM>::DataVector *v=(*Dat).d[l] ;
00728   if (*i==MAXINT) {
00729     *i=(long)(*v).size() ;
00730     (*v).push_back(x) ;
00731     CounterArray[l]=(size_t)((*i)+1) ;
00732   } else {
00733     (*v)[*i]=x ;
00734   } 
00735 }
00736 */
00737 
00738 
00739 template<int DIM>
00740 void AdaptiveGrid<DIM>::LargestActiveLevel(int *J) {
00741   int i ;
00742   for (i=0; i<DIM; i++) J[i]=0 ;
00743 
00744   Matrix<size_t,DIM>::iterator it(&CounterArray) ;
00745   for (it=CounterArray.begin(); it<=CounterArray.end(); ++it)
00746     if (CounterArray[it]>0) 
00747       for (i=0; i<DIM; i++) if (J[i]<it.i[i]) J[i]=it.i[i] ;
00748 }
00749 
00750 //
00751 // ProjectionGrid-Stuff
00752 //
00753 
00754 template<int DIM>
00755 void AdaptiveGrid<DIM>::ClearPG(AdaptiveGrid<DIM>::ProjectionGrid *PG, int dir) {
00756   ProjectionGrid::iterator pgit ;
00757 
00758   for (pgit=(*PG).begin(); pgit!=(*PG).end(); ++pgit) {
00759     AdaptiveLevelsIndex *ali=(*pgit).first  ;
00760     AdaptiveLevels      *als=(*pgit).second ;
00761 
00762     if (dir>=0) {
00763       DumpPG(PG,dir) ;
00764       std::cout<< "ClearPG "<<this<<' '<<DIM<<' '<<(*PG).size()<<' '<<als<<'\n' ;
00765     }
00766 
00767     AdaptiveLevelsIndexAllocator.DeleteItem(ali) ;
00768     AdaptiveLevelsAllocator.     DeleteItem(als) ;
00769   }
00770   (*PG).clear() ;
00771 }
00772 
00773 //
00774 // GeneratePG generates the a PG Datastructure for accessing the Data
00775 // with respect to the 'dTo' direction. 
00776 // The Base is the PG Datastructure for accessing with respect to the 'dFrom' 
00777 // direction
00778 //  
00779 
00780 template<int DIM>
00781 size_t AdaptiveGrid<DIM>::GeneratePG(ProjectionGrid *PGF,int dFrom , ProjectionGrid *PGT,int dTo) {
00782   size_t  numberofentries=0 ;
00783   int     i, l[DIM], t[DIM], ld ; 
00784 
00785   ClearPG(PGT) ;
00786   ProjectionGrid::iterator pgit , pgend=(*PGF).end() ; 
00787 
00788   for (pgit=(*PGF).begin(); pgit!=pgend; ++pgit) {
00789     AdaptiveLevelsIndex       *ali=(*pgit).first  ;
00790     AdaptiveLevels            *als=(*pgit).second ;
00791     AdaptiveLevels::AdaptiveL *al ;
00792     AdaptiveLevels::AdaptiveL::iterator alit,alend ;
00793 
00794     // calculate level and index
00795     for (i=0    ; i<dFrom; i++) l[i  ]=ali->l[i], t[i  ]=ali->t[i] ;
00796     for (i=dFrom; i<DIM-1; i++) l[i+1]=ali->l[i], t[i+1]=ali->t[i] ;
00797 
00798     for (ld=Level0[dFrom]; ld<=LMAX; ld++) {
00799       l[dFrom]=ld ;
00800       al      =(*als)[ld]  ;
00801       alend   =(*al).end() ;
00802       size_t alsize=0 ;
00803       for (alit=(*al).begin(); alit!=alend; alit++) {
00804         // complete index
00805         t[dFrom]=(*alit).first ;
00806         // insert   (l,t) with known index to numerical data (*alit).second
00807         Insert(PGT,dTo,l,t,(*alit).second) ;
00808         alsize++ ;
00809       }
00810       // some consistency checks
00811       if (alsize!=(*al).size()) std::cout<< "Error in GeneratePG: alsize="<<alsize<<" , al.size="<<(*al).size()<<'\n' ;
00812       numberofentries += alsize ;
00813     }
00814   }
00815 
00816   return numberofentries ;
00817 }
00818 
00819 
00820 template<int DIM>
00821 void AdaptiveGrid<DIM>::GenerateAllBasisPGs() {
00822   size_t noe,sall=Size() ;
00823 
00824   for(int d=1;d<DIM;d++) {
00825     noe=GeneratePG(BasisPG[0],0,BasisPG[d],d) ;
00826     if (noe != sall) {
00827       std::cout<< "Error in GenerateAllBasisPGs: noe="<<noe<<" , Size=",sall<<'\n' ;
00828       PrintCounterArray() ;
00829     }
00830   }
00831 
00832 }
00833 
00834 
00835 template<int DIM>
00836 void AdaptiveGrid<DIM>::DumpPG(ProjectionGrid *pg,int dir) {
00837   int     i, l[DIM], t[DIM], ld, n=0 ;
00838 
00839   ProjectionGrid::iterator pgit , pgend=(*pg).end() ; 
00840 
00841   for (pgit=(*pg).begin(); pgit!=pgend; ++pgit) {
00842     AdaptiveLevelsIndex       *ali=(*pgit).first  ;
00843     AdaptiveLevels            *als=(*pgit).second ;
00844     AdaptiveLevels::AdaptiveL *al ;
00845     AdaptiveLevels::AdaptiveL::iterator alit,alend ;
00846 
00847     std::cout<< "DumpPG("<<dir<<") AL(S/I) "<<als<<' '<<ali<<'\n' ;
00848 
00849     // calculate level and index
00850     for (i=0  ; i<dir  ; i++) l[i  ]=ali->l[i], t[i  ]=ali->t[i] ;
00851     for (i=dir; i<DIM-1; i++) l[i+1]=ali->l[i], t[i+1]=ali->t[i] ;
00852 
00853     for (ld=Level0[dir]; ld<=LMAX; ld++) {
00854       l[dir]=ld ;
00855       al    =(*als)[ld]  ;
00856       alend =(*al).end() ;
00857       //size_t alsize=0 ;
00858       for (alit=(*al).begin(); alit!=alend; alit++) {
00859         n++ ;
00860         // complete index
00861         t[dir]   =(*alit).first ;
00862         size_t si=(*alit).second ;  
00863         std::cout<< "DumpPG("<<dir<<"): (" ;
00864         if (DIM>0) { 
00865           std::cout<<l[0] ;
00866           for (i=1; i<DIM; i++) std::cout<<','<<l[i] ;
00867         }
00868         std::cout<<"),(" ;
00869         if (DIM>0) { 
00870           std::cout<<t[0] ;
00871           for (i=1; i<DIM; i++) std::cout<<','<<t[i] ;
00872         }
00873         std::cout<<") "<<si<<'\n' ;
00874       }
00875     }
00876   }
00877 
00878   std::cout<<"DumpPG("<<dir<<") Size "<<n<<'\n' ;
00879 }
00880 
00881 
00882 template<int DIM>
00883 void AdaptiveGrid<DIM>::WriteSparse(const char *name) {
00884   tmp[0]->Set(1.0) ;
00885   tmp[0]->WriteSparse(name) ;
00886 }
00887 
00888 #include"SMP.hpp"
00889 
00890 #endif
Generated at Mon Aug 19 10:02:31 2002 for AWFD by