Adaptive1D.cc

00001 //: Adaptive1D.cc
00002 
00003 #include<stdlib.h>
00004 #include<math.h>
00005 #include<assert.h>
00006 
00007 #include "Adaptive1D.hpp"
00008 #include "NonAdaptive.hpp"
00009 #include "Buffer.hpp"
00010 extern int debugRefine ;
00011  
00012 bool intlt::operator()(const int i,const int j) const { return i<j ;}
00013 
00014 AdaptiveB::AdaptiveLevel::AdaptiveLevel()  { d=NULL  ; IS=NULL ; Size=0 ; ISowner=false ; } 
00015 AdaptiveB::AdaptiveLevel::~AdaptiveLevel() { delete []d; d=NULL; if (ISowner) delete IS ; IS=NULL ; }
00016 
00017 int AdaptiveB::AdaptiveLevel::InitD(int size) {
00018   Size=size     ;
00019   ISowner=false ;
00020   if(d==NULL) {
00021     d = new double[size+1] ;
00022   } else {;}
00023   for (int i=0; i<=size; i++) d[i]=MNAN ;
00024   return 1  ;
00025 }
00026 
00027 int AdaptiveB::AdaptiveLevel::Init(int size) {
00028   if(!InitD(size)) return 0 ;
00029   IS=new IndexSet(size) ;
00030   ISowner=true ;
00031   return 1 ;
00032 }
00033 
00034 AdaptiveGS::AdaptiveLevel::AdaptiveLevel()  { d=NULL   ; Size=0 ;}
00035 AdaptiveGS::AdaptiveLevel::~AdaptiveLevel() { delete []d; }
00036 
00037 int AdaptiveGS::AdaptiveLevel::Init(int size) {
00038   Size=size ;
00039   
00040   if(!I.Init(size))   return 0  ;
00041   if(!II.Init(size))  return 0  ;
00042   if(!III.Init(size)) return 0  ;
00043   if(!IV.Init(size))  return 0  ;
00044 
00045   if(d==NULL) { 
00046     d = new double[size+1] ;
00047     for (int i=0; i<size; i++) d[i]=MNAN ;
00048    }
00049   return 1  ;
00050 }
00051 
00052 int AdaptiveB::AdaptiveLevel::Print(char *st,int level,int allflag) {
00053   int c,i,r ;
00054   c=IS->Print() ; 
00055   if (allflag &&(IS->nr>0)) for (i=0;i<min(Size,1065);i++) std::cout << st << '(' << level<<','<<i << "): "<< d[i] <<'\n';
00056   else {
00057     for (r=0; r<IS->nr; r++)
00058       for (i=IS->a[r]; i<=IS->e[r]; i++) std::cout << st << '(' << level<<','<<i << "): "<< d[i] <<'\n';
00059   }
00060   return c ;
00061 }
00062 
00063 int AdaptiveGS::AdaptiveLevel::Print(char *st,int level,int allflag) {
00064 int c ;
00065 
00066   c=I.Print() ; 
00067   II.Print() ; 
00068   III.Print() ; 
00069   IV.Print() ; 
00070   if (allflag &&(I.nr>0)) for (int i=0;i<min(Size,1065);i++) std::cout << st << '(' << level<<','<<i << "): "<< d[i] <<'\n';
00071   return c ;
00072 }
00073 
00074 /***************************************/
00075 /*                                     */
00076 /*   Funktionen fuer AdaptiveB::       */
00077 /*                                     */
00078 /***************************************/
00079 
00080 AdaptiveB::AdaptiveB () { 
00081   Level0=J=0 ;
00082   BC[0]=BC[1]=-MAXINT ; 
00083   XA=0.0 ; XE=1.0 ; 
00084   islifting   =false ; 
00085   ismultiscale=true ;
00086   for (unsigned int i=0; i<sizeof(Buffers)/sizeof(Buffers[0]); i++) Buffers[i]=NULL ;
00087 }
00088 
00089 AdaptiveB::AdaptiveB (AdaptiveB *B) {
00090   Level0=B->Level0 ; J=B->J ;
00091   BC[0]=B->BC[0]  ; BC[1]=B->BC[1] ;
00092 
00093   XA=0.0 ; XE=1.0 ;
00094   islifting   =false ;
00095   ismultiscale=true  ;
00096 
00097   int c,i;
00098   c=L[Level0].InitD((1<<Level0)+1) ; assert(c) ;
00099 
00100   for (i=Level0+1;i<=J;i++) { 
00101     c=L[i].InitD(1<<(i-1)) ;
00102     if (c==0) {std::cout << "AdaptiveB::AdaptiveB(AdaptiveB *B)  Could not allocate memory\n" ; exit(-1) ;}
00103   }
00104 
00105   for (i=Level0; i<=J; i++) { L[i].IS=B->L[i].IS ; assert(L[i].IS) ;} 
00106 
00107   for (i=0; i<(int)(sizeof(Buffers)/sizeof(Buffers[0])); i++) Buffers[i]=B->Buffers[i] ;
00108 }
00109 
00110 AdaptiveB::~AdaptiveB() { 
00111   Level0=J=0 ;
00112 }
00113 
00114 int AdaptiveB::Init(int Lev0,int j, double xa, double xe) {
00115   XA=xa, XE=xe ;
00116   if (Level0) { assert((Lev0==Level0)&&(j==J)) ; return 1 ;} // schon irgendwannmal initialisiert worden
00117 
00118   Level0=Lev0 ; J=j ; 
00119   if(L[Level0].Init((1<<Level0)+1)==0) return 0 ;
00120 
00121   for (int i=Level0+1;i<=J;i++)
00122     if(L[i].Init(1<<(i-1))==0) return 0 ;
00123 
00124   return 1 ;
00125 }
00126 
00127 void AdaptiveB::Copy(AdaptiveB *B) {
00128  int l ;
00129  BC[0] =B->BC[0] ;
00130  BC[1] =B->BC[1] ;
00131  Level0=B->Level0 ;
00132  J     =B->J ; 
00133  XA    =B->XA ;
00134  XE    =B->XE ;
00135    
00136  CopyIndexSets(B) ;
00137  for (l=Level0; l<=J; l++) L[l].IS->VecCopy(L[l].d , B->L[l].d) ;
00138 }
00139 
00140 
00141 void AdaptiveB::SetBoundaryConditions(int *BCs)        { BC[0]=BCs[0] ; BC[1]=BCs[1] ;}
00142 void AdaptiveB::SetBoundaryConditions(int bc0,int bc1) { BC[0]=bc0    ; BC[1]=bc1    ;} 
00143 
00144 //void AdaptiveB::Print() { Print(NULL,0) ;}
00145 void AdaptiveB::Print(int allflag) { Print("Bd",allflag) ;}
00146 void AdaptiveB::Print(char *st,int allflag) {
00147   int c=0 ;
00148   std::cout << "Level0: "<<Level0<<" J: "<<J<<" BC "<<BC[0]<<','<<BC[1]<<'\n' ;
00149   for (int l=Level0;l<=J;l++) {
00150     std::cout << "Level "<<l<<'\n'; 
00151     c+=L[l].Print(st,l,allflag) ;
00152   }
00153   std::cout << "NZE: "<<c<<'\n';
00154 }
00155 
00156 //
00157 // print only the active IndexSets in the following fashion:
00158 // an output with three columns is generated:   *st , x(l,i) , l
00159 // where x(l,i)=i*2^(-l) 
00160 void AdaptiveB::Print(char *st) {
00161   for (int l=Level0;l<=J;l++) L[l].IS->Print(st,l,Level0) ;
00162 }
00163 
00164 
00165 size_t AdaptiveB::Size() {
00166   size_t s=0 ;
00167   for (int l=Level0; l<=J; l++) s += L[l].IS->Size() ; 
00168   return s ;
00169 }
00170 
00171 // Calculates IndexSet of nonvanishing coefficients 
00172 // the coefficients are stored in d[0...Size-1] 
00173 // if any coefficient with index < K0        is active, the complete interval [0..K-1] is activated
00174 // if any coefficient with index > Size-1-K1 is active, the complete interval [size-K1..size-1] is act.
00175 //
00176 // if Insert0/1 ist true then the interval [0..K-1] / [size-K1..size-1] is always activated
00177 // 
00178 void IndexSet::CalcIndexSet(double *d,int Size,int K0,int K1,bool Insert0,bool Insert1) {
00179  int i,aa=0,ee=-1 ;
00180 
00181  nr=0 ;
00182  // condition for K0
00183  i=0 ;
00184  while (i<K0) {
00185    if (d[i]!=0) Insert0=true ;
00186    i++ ;
00187  }
00188 
00189  i=max(K0,0) ;
00190  if (Insert0) {
00191      a[0]=0 ;
00192      while ((i<Size) && (d[i]!=0)) i++ ;
00193      e[0]=i-1 ;
00194      nr=1 ;
00195  }
00196    
00197  while (i<Size) {
00198    while ((i<Size) && (d[i]==0)) i++; 
00199    aa=i   ;
00200    while ((i<Size) && (d[i]!=0)) i++; 
00201    ee=i-1 ;
00202    if (aa<=ee) a[nr]=aa , e[nr]=ee , nr++ ;
00203  }
00204 
00205  // condition for K1
00206  i=nr-1 ;
00207  while ((i>=0) && (e[i]>Size-1-K1)) i-- ;
00208 
00209  // the intervals [i+1..nr-1] match the condition: cleap them up 
00210  if (i+1<=nr-1) {
00211    if (a[i+1]>Size-1-K1) a[i+1]=Size-K1 ;
00212    e[i+1]=Size-1 ;
00213    nr=i+2 ;
00214  } else {
00215    if (Insert1) { a[nr]=Size-K1 ; e[nr++]=Size-1 ;} // insert a new interval
00216  }
00217 
00218  // Patch: if the last interval and it's predecessor touch, merge them
00219  if (nr>1)
00220    if (a[nr-1]==e[nr-2]+1) { e[nr-2]=e[nr-1] ; nr--;}
00221 }
00222 
00223 
00224 void AdaptiveB::CopyIndexSets(AdaptiveB *S)
00225 {
00226  int l ;
00227  for (l=Level0;l<=J;l++) L[l].IS->Copy(S->L[l].IS) ;
00228 }
00229 
00230 bool AdaptiveB::SameIndexSets(AdaptiveB *S)
00231 {
00232  int l ;
00233  for (l=Level0;l<=J;l++) if (!L[l].IS->IsSame(S->L[l].IS)) return false ;
00234  return true ;
00235 }
00236 
00237 void AdaptiveB::FromFull(double *d, int *BCs, Wavelets *W, double eps) {
00238   class AdaptivityCriterion R(AdaptivityCriterion::L2_THRESHOLD,eps) ;
00239   FromFull(d,BCs,W,&R) ;
00240 }
00241 
00242 void AdaptiveB::FromFull(double *d, int *BCs, Wavelets *W , AdaptivityCriterion *R) {
00243   int    i,l,size ;
00244   double x     ;
00245   int    K0=-1,K1=-1,KW0=-1,KW1=-1 ;
00246   bool   Insert0,Insert1,NotPeriodic=(BCs[1]!=PERIODIC) ;
00247   
00248   if (NotPeriodic) {
00249     K0=W->K0 ;
00250     K1=W->K1 ;
00251     KW0=W->KW0 ;
00252     KW1=W->KW1 ;
00253   }
00254 
00255   BC[0]=BCs[0] , BC[1]=BCs[1] ;
00256 
00257   bool BoundC=false ;   
00258   // this condition is neccessary if also in the adaptive case, the 
00259   // set of adaptive wavelets shall be large enough to represent P_(V_0)(u)
00260   // and also the imbedding P_V(u) 
00261   // this not trivial !!!   
00262 
00263   Insert0=Insert1=false ;
00264   for (l=J; l>Level0; l--) {
00265     size=1<<(l-1) ;
00266     for (i=0;i<size;i++) x=d[i+size+1] , L[l].d[i]= (R->IsActive(x,&l,W->Level0,1)) ? x : 0 ;
00267     L[l].IS->CalcIndexSet(L[l].d,size,KW0,KW1,Insert0&BoundC,Insert1&BoundC) ;
00268     
00269     // Check whether boundary wavelets on lower levels must be activated
00270     if (L[l].IS->nr>0) {
00271       Insert0 =Insert0 || ((L[l].IS->a[0            ]==0     ) && NotPeriodic) ;
00272       Insert1 =Insert1 || ((L[l].IS->e[L[l].IS->nr-1]==size-1) && NotPeriodic) ;
00273     }
00274   }
00275   
00276   // copy Coefficients for Wavelet Basis 
00277   size=(1<<Level0) ;
00278   if (NotPeriodic) size++ ;
00279   for (i=0;i<size;i++)   x=d[i] , L[Level0].d[i]= (R->IsActive(x,&Level0,W->Level0,1))  ? x : 0 ;
00280   L[Level0].IS->CalcIndexSet(L[Level0].d,size,K0,K1,Insert0,Insert1) ;
00281   /*
00282   L[Level0].IS->nr=1      ;
00283   L[Level0].IS->a[0]=0    ;
00284   L[Level0].IS->e[0]=size ;
00285   */
00286 }
00287 
00288 
00289 void AdaptiveB::FromFull(double *d) {
00290   int l,*a,*e,r,i,be ;
00291   for (l=Level0; l<=J; l++) {
00292     a=L[l].IS->a;
00293     e=L[l].IS->e;
00294     be = (l==Level0) ? 0 : (1<<(l-1))+1 ;
00295     for (r=0; r<L[l].IS->nr; r++)
00296       for (i=a[r]; i<=e[r]; i++) 
00297         L[l].d[i]=d[be+i] ;
00298   }
00299 }
00300 
00301 
00302 void AdaptiveB::ToFull(double *d) {
00303   int i,l,o  ; 
00304   double *d0 ;
00305 // copy Coefficients for Wavelet Basis 
00306   
00307   for (l=Level0; l<=J; l++) {
00308     d0= L[l].d ;
00309     o= (l==Level0) ? 0 : L[l].Size+1 ; 
00310     for (i=0;i<L[l].Size;i++) d[i+o]=0.0 ;
00311     for (int r=0; r<L[l].IS->nr; r++)
00312       for (int j=L[l].IS->a[r];j<=L[l].IS->e[r];j++) d[j+o]=d0[j] ;
00313    }  
00314 }
00315 
00316 void AdaptiveB::FromAdaptiveGS(AdaptiveGS *GS , Wavelets *W) // Erzeugt Basis aus Adaptivem Erzeugenden System
00317 {                                                                       // Erzeugenden Sytem wird veraendert !!!
00318  int l ;
00319 
00320  BC[0]=GS->BC[0] ,  BC[1]=GS->BC[1] ;
00321  XA=GS->XA ; XE=GS->XE ;
00322 
00323  // calculate all levels from up to down 
00324  for (l=J; l>Level0 ;l--) {
00325    W->GtJ.MatVec(    L[l  ].d,     L[l  ].IS   ,  GS->L[l].d,&GS->L[l].IV,BC ,l,1, Buffers[0]) ;
00326    W->HtJ.MatVec(GS->L[l-1].d,&GS->L[l-1].III  ,  GS->L[l].d,&GS->L[l].I ,BC ,l,1, Buffers[0]) ; // fuer Level0 eigentlich nur fuer III geschnitten IV
00327  }
00328  L[Level0].IS->VecCopy(L[Level0].d , GS->L[Level0].d) ;
00329  
00330  BoundaryCorrection(BC,L[Level0].d,Level0) ;
00331 }
00332 
00333 //
00334 // Test-Funktion
00335 // 
00336 void AdaptiveB::FromAdaptiveGS2(AdaptiveGS *GS , Wavelets *W) // Erzeugt Basis aus Adaptivem Erzeugenden System
00337 {                                                                       // Erzeugenden Sytem wird veraendert !!!
00338  int l ;
00339  
00340  BC[0]=GS->BC[0] ,  BC[1]=GS->BC[1] ;
00341  XA=GS->XA ; XE=GS->XE ;
00342 
00343  // additional values of the generating system are calculated by the  
00344  // canonic interpolation (one step of inverse wavelet transform) 
00345 
00346  for (l=Level0;l<J;l++) {
00347    GS->L[l+1].II.InducedByCoarseScalings(&GS->L[l].I,l+1,BC,W) ;
00348    GS->L[l+1].III.Difference(&GS->L[l+1].II , &GS->L[l+1].I)   ;
00349    W->HtJ.MatTVec(GS->L[l+1].d,&GS->L[l+1].III  ,  GS->L[l].d,&GS->L[l].I,BC, l ,1, Buffers[0]) ;
00350  } 
00351 
00352  // Berechne alle Level von oben nach unten 
00353  for (l=J; l>Level0 ;l--) {
00354    W->GtJ.MatVec(    L[l  ].d,     L[l  ].IS   ,  GS->L[l].d,&GS->L[l].I ,BC ,l,1, Buffers[0]) ;
00355    W->HtJ.MatVec(GS->L[l-1].d,&GS->L[l-1].I    ,  GS->L[l].d,&GS->L[l].II,BC ,l,1, Buffers[0]) ; // fuer Level0 eigentlich nur fuer III geschnitten IV
00356  }
00357  L[Level0].IS->VecCopy(L[Level0].d , GS->L[Level0].d) ;
00358 
00359  BoundaryCorrection(BC,L[Level0].d,Level0) ;
00360 }
00361 
00362 // AdaptiveGS contains both the index sets S(l) for the adaptive generating systems as well as the nodal values u^S(l)
00363 // starting from these values we compute the wavelet coefficients by successive low and high pass filter operations.
00364 // If (islifting==true) lifting steps are included
00365 // note: the values in *GS are modified
00366 void AdaptiveB::FromAdaptiveGS3(AdaptiveGS *GS, Wavelets *W)
00367 {
00368  int l ;
00369  
00370  BC[0]=GS->BC[0] ,  BC[1]=GS->BC[1] ; 
00371  XA=GS->XA ; XE=GS->XE ;
00372 
00373  // start from high levels to low levels 
00374  for (l=J; l>Level0 ;l--) {
00375 
00376    W->GtJ.MatVec(    L[l  ].d,     L[l  ].IS   ,  GS->L[l].d,&GS->L[l].I ,BC ,l,1, Buffers[0]) ; // high pass filter
00377 
00378    GS->L[l-1].III.IndexSetForTypeIII(&GS->L[l].I,l, BC, W) ;                                     // low pass filter
00379    W->HtJ.MatVec(GS->L[l-1].d,&GS->L[l-1].III  ,  GS->L[l].d,&GS->L[l].I ,BC ,l,1, Buffers[0]) ;
00380 
00381    if (islifting) W->Q.MatTVec(GS->L[l-1].d,&GS->L[l-1].I  ,  L[l].d,L[l].IS , BC , l-1 , -1, Buffers[0]) ;
00382  }
00383 
00384  L[Level0].IS->VecCopy(L[Level0].d , GS->L[Level0].d) ;
00385  BoundaryCorrection(BC,L[Level0].d , Level0) ;
00386 
00387  if (islifting) for (l=Level0; l<=J; l++) L[l].IS->VecMul(L[l].d , 1./sqrt((double)(1<<l))) ;
00388 }
00389 
00390 
00391 void AdaptiveB::AdditiveFromAdaptiveGS(AdaptiveGS *GS, Wavelets *W) // Erzeugt Basis aus Adaptivem Erzeugenden System
00392 {                                                                       // Erzeugenden Sytem wird veraendert !!!
00393  int l ;
00394 
00395  BC[0]=GS->BC[0] ,  BC[1]=GS->BC[1] ;
00396  
00397  // Berechne alle Level von oben nach unten 
00398  for (l=J; l>Level0 ;l--) {
00399    W->GJ.MatVec(    L[l  ].d,     L[l  ].IS ,  GS->L[l].d,&GS->L[l].I ,BC ,l,1, Buffers[0]) ;
00400    W->HJ.MatVec(GS->L[l-1].d,&GS->L[l-1].I  ,  GS->L[l].d,&GS->L[l].I ,BC ,l,0, Buffers[0]) ; // fuer Level0 eigentlich nur fuer III geschnitten IV
00401  }
00402  L[Level0].IS->VecCopy(L[Level0].d , GS->L[Level0].d) ;
00403 }
00404 
00405 // HighPassFilterFromAdaptiveGS generates a wavelet basis representation by simply  
00406 // performing the high pass filtering
00407 void AdaptiveB::HighPassFilterFromAdaptiveGS(AdaptiveGS *GS, Wavelets *W) {
00408  int l ;
00409 
00410  BC[0]=GS->BC[0] ,  BC[1]=GS->BC[1] ;
00411  
00412  for (l=J; l>Level0 ;l--) {
00413    W->GJ.MatVec( L[l].d,L[l].IS ,  GS->L[l].d,&GS->L[l].I ,BC ,l,1, Buffers[0]) ;
00414  }
00415  L[Level0].IS->VecCopy(L[Level0].d , GS->L[Level0].d) ;
00416 }
00417 
00418 
00419 void AdaptiveB::IndexSetFromAdaptiveGS(AdaptiveGS *G, Wavelets *WC)
00420 {
00421  int l ;
00422  for (l=J;l>=Level0;l--) 
00423    L[l].IS->WaveletsInducedByFineScalings(&G->L[l].I,l,G->BC,WC) ;
00424 }
00425 
00426 
00427 //
00428 // generate an adaptive basis, which fulfills the regularity condition
00429 // for interpolet generating systems
00430 // This amounts fills up from fine to coarse levels
00431 // 
00432 void AdaptiveB::InterpoletRegularity(Wavelets *W) 
00433 {
00434  int  l,r,su,so ;
00435  IndexSet I1(4000),I2(4000) , IE3(4000) ;
00436  bool NotPeriodic,Insert0,Insert1 ;
00437  
00438  if (debugRefine) { std::cout << "IP "<<Level0<<'\n'; L[Level0].IS->Print("IP") ; }
00439 
00440  NotPeriodic=(BC[1]!=PERIODIC) ;
00441  Insert0=Insert1=false ;
00442  
00443  IE3.nr=0 ; //empty
00444  for (l=J-1; l>=Level0; l--) {
00445    // compute scaling functions on level l induced by Lifting-Wavelets on level l+1  and by  scaling functions on level l+1 
00446    if (islifting) {
00447      I1.ScalingsInducedByLifting(L[l+1].IS , l+1 , BC , W) ;
00448      I2.InducedByFineScalings   (&IE3      , l+1 , BC , W) ;
00449      IE3.Union(&I1 , &I2) ;
00450    }
00451 
00452    // compute scaling functions on level l+1 which span the wavelets on level l+1;
00453    // these quite immediately lead to the fathers of the level l+1 wavelets on level l
00454    I2.InducedByWavelets(L[l+1].IS,l+1,BC,W) ;
00455 
00456    // fathers on level l ...
00457    for (r=0; r<I2.nr; r++) {
00458      su=I2.a[r]/2 ; // index set of associated gridpoints on next coarser level
00459      so=I2.e[r]/2 ;
00460      
00461      if (l>Level0) {
00462        su=(su&1) ? su : su+1 ;
00463        so=(so&1) ? so : so-1 ;
00464 
00465        I1.a[r]=su/2 ;
00466        I1.e[r]=so/2 ;
00467      } else {
00468        I1.a[r]=su ;
00469        I1.e[r]=so ;
00470      }
00471    }
00472    I1.nr=I2.nr ;
00473 
00474    // I1 may contain overlapping intervals, merge these
00475    if (I1.nr==0) I2.nr=0 ;
00476    else {
00477      I2.a[0]=I1.a[0]; 
00478      I2.e[0]=I1.e[0];  
00479      I2.nr=1 ;
00480      for (r=1; r<I1.nr; r++) {
00481        if (I2.e[I2.nr-1]>=I1.a[r]-1) { I2.e[I2.nr-1]=I1.e[r]; } // touch or overlap
00482        else {
00483          I2.a[I2.nr]=I1.a[r], // copy interval
00484          I2.e[I2.nr]=I1.e[r], 
00485          I2.nr++ ;
00486         }
00487      }
00488    }
00489 
00490    // union of: old basis , fathers , RequiredForScalings on Level l
00491    if (IE3.nr>0) {
00492      I1.Union(L[l].IS , &I2) ;
00493      I2.WaveletsInducedByFineDualScalings(&IE3 , l , BC , W) ;
00494    } else 
00495      I1.Copy(L[l].IS) ;
00496 
00497    L[l].IS->Union(&I1 , &I2) ;
00498 
00499    // BoundCondition
00500    if (NotPeriodic) {
00501      int sl,K0,K1 ;
00502      if (l>Level0) { K0=W->KW0-1; K1=(1<<(l-1))-W->KW1 ; sl=(1<<(l-1))-1; }
00503      else          { K0=W->K0-1 ; K1=(1<<l) +1 -W->K1  ; sl= 1<<l       ; }
00504      L[l].IS->SatisfyBoundCondition( sl , -MAXINT , K0 , K1 ,&Insert0,&Insert1) ; 
00505    }
00506  }
00507 
00508 }
00509 
00510 
00511 //
00512 // Load only the important coefficients of the GS in a basis data structure
00513 // this requires the regularity condition from above
00514 //
00515 void AdaptiveB::LoadFromInterpoletGS(AdaptiveGS *G, Wavelets *)
00516 {
00517  int l,r,i,*a,*e,n ;
00518  double *d,*d2 ;
00519  
00520  BC[0]=G->BC[0] ;
00521  BC[1]=G->BC[1] ;
00522 
00523  L[Level0].IS->VecCopy(L[Level0].d , G->L[Level0].d) ;
00524  for (l=Level0+1; l<=J; l++) {
00525    a=L[l].IS->a ;
00526    e=L[l].IS->e ;
00527    n=L[l].IS->nr ;
00528    d=L[l].d      ;
00529    d2=G->L[l].d ;
00530    for (r=0; r<n; r++)
00531      for (i=a[r];i<=e[r];i++) d[i]=d2[2*i+1] ;
00532 
00533    if (BC[0]==1) d[0]=d2[2] ;
00534    if (BC[1]==1) d[(1<<(l-1))-1]=d2[(1<<l)-2] ;
00535  }
00536 }
00537 
00538 //
00539 // Load only the important coefficients of the GS in a basis data structure
00540 // this requires the regularity condition from above
00541 // We can not assume that the GS entries of intermediate Levels are still 
00542 // interpolating, therefore we have to go really from the finest to the coarsest
00543 // Level
00544 //
00545 
00546 void AdaptiveB::LoadFromLiftingInterpoletGS(AdaptiveGS *G, Wavelets *)
00547 {
00548  int l,r,i,*a,*e,n  ;
00549  double *d,*ds,*dt ;
00550  
00551  BC[0]=G->BC[0] ;
00552  BC[1]=G->BC[1] ;
00553 
00554  for (l=J; l>Level0; l--) {
00555    a =G->L[l].I.a  ;
00556    e =G->L[l].I.e  ;
00557    n =G->L[l].I.nr ;
00558    d =L[l].d       ;
00559    ds=G->L[l  ].d  ;
00560    dt=G->L[l-1].d  ;
00561 
00562    for (r=0; r<n; r++)
00563    for (i=a[r]; i<=e[r]; i++)
00564      if (i&1) d [i>>1]=ds[i] ;
00565         else  dt[i>>1]=ds[i] ; 
00566      
00567    if (BC[0]==1) d[0]=ds[2] ;
00568    if (BC[1]==1) d[(1<<(l-1))-1]=ds[(1<<l)-2] ;
00569  }
00570 
00571  L[Level0].IS->VecCopy(L[Level0].d , G->L[Level0].d) ;
00572 
00573 }
00574 
00575 
00576 //
00577 // RestoreFromInterpoletBasis restores the complete generating system from the compact
00578 // representation in the basis data structure
00579 //
00580 void AdaptiveGS::RestoreFromInterpoletBasis(AdaptiveB *B, Wavelets *) {
00581   RestoreFromInterpoletBasis(B, this, NULL) ; 
00582 }
00583 
00584 /*
00585 void AdaptiveGS::RestoreFromInterpoletBasis(class AdaptiveB *B,struct Wavelets *)
00586 {
00587  int l,r,i,*a,*e,n,er ;
00588  double *d,*d2 ; 
00589 
00590  BC[0]=B->BC[0] ;
00591  BC[1]=B->BC[1] ;
00592 
00593  // 
00594  // restore the representatives
00595  L[Level0].I.VecFill(L[Level0].d , 0.0) ;
00596  B->L[Level0].IS->VecCopy(L[Level0].d , B->L[Level0].d) ;
00597 
00598  for (l=Level0+1; l<=J; l++) {
00599    a=B->L[l].IS->a ;
00600    e=B->L[l].IS->e ;
00601    n=B->L[l].IS->nr ;
00602    d=B->L[l].d ;
00603    d2=L[l].d   ;
00604    for (r=0; r<n; r++) {
00605      er=e[r];
00606      for (i=a[r];i<=er;i++) d2[2*i+1]=d[i] ;
00607    }
00608 
00609    if (BC[0]==1) {d2[2]=d[0]                   ;d2[1]=0.0 ;}
00610    if (BC[1]==1) {d2[(1<<l)-2]=d[(1<<(l-1))-1] ;d2[(1<<l)-1]=0 ;}
00611  } 
00612 
00613  //
00614  // restores the complete rest
00615  for (l=Level0+1; l<=J; l++) {
00616    a=L[l].I.a  ;
00617    e=L[l].I.e  ;
00618    n=L[l].I.nr ;
00619    d=L[l].d    ;
00620    d2=L[l-1].d ;
00621 
00622    double w0=d[2],w1=d[(1<<l)-2] ;
00623 
00624    for (r=0; r<n; r++) {
00625      er=e[r]&(0xfffffe) ;
00626      for (i=(a[r]+1)&(0xfffffe) ; i<=er ; i+=2) {
00627        d[i]=d2[i/2] ;
00628      }
00629    }
00630 
00631    if (BC[0]==1) d[2]=w0 ;
00632    if (BC[1]==1) d[(1<<l)-2]=w1 ;
00633 
00634  }
00635 }
00636 */
00637 
00638 //
00639 //
00640 void AdaptiveGS::RestoreFromInterpoletBasis(AdaptiveB *B , AdaptiveGS * G, Wavelets *W) {
00641  int l,r,i,*a,*e,n,er ;
00642  double *d,*d2 ; 
00643 
00644  BC[0]=B->BC[0] ;
00645  BC[1]=B->BC[1] ;
00646 
00647  if (this != G) for (l=Level0+1; l<=J; l++) L[l].I.VecFill(L[l].d , 1e+133) ;
00648 
00649  // 
00650  // restore the representatives
00651  L[Level0].I.VecFill(L[Level0].d , 0.0) ;
00652  B->L[Level0].IS->VecCopy(L[Level0].d , B->L[Level0].d) ;
00653 
00654  for (l=Level0+1; l<=J; l++) {
00655    a=B->L[l].IS->a ;
00656    e=B->L[l].IS->e ;
00657    n=B->L[l].IS->nr ;
00658    d=B->L[l].d ;
00659    d2=L[l].d   ;
00660    for (r=0; r<n; r++) {
00661      er=e[r];
00662      for (i=a[r];i<=er;i++) d2[2*i+1]=d[i] ;
00663    }
00664 
00665    if (BC[0]==1) {d2[2]=d[0]                   ;d2[1]=0.0 ;}
00666    if (BC[1]==1) {d2[(1<<l)-2]=d[(1<<(l-1))-1] ;d2[(1<<l)-1]=0 ;}
00667  } 
00668 
00669  //
00670  // restores the rest in the Basis-generating system 
00671  for (l=Level0+1; l<=J; l++) {
00672    a=G->L[l].I.a  ;
00673    e=G->L[l].I.e  ;
00674    n=G->L[l].I.nr ;
00675    d =L[l  ].d    ;
00676    d2=L[l-1].d    ;
00677 
00678    double w0=d[2],w1=d[(1<<l)-2] ;
00679 
00680    for (r=0; r<n; r++) {
00681      er=e[r]&(0xfffffe) ;
00682      for (i=(a[r]+1)&(0xfffffe) ; i<=er ; i+=2) {
00683        d[i]=d2[i/2] ;
00684      }
00685    }
00686 
00687    if (BC[0]==1) d[2]=w0 ;
00688    if (BC[1]==1) d[(1<<l)-2]=w1 ;
00689 
00690  }
00691 
00692  // 
00693  // now the additional entries are computed by interpolet interpolation
00694  if (this == G) return ; // nothing to do
00695  
00696  double *t=Buffers[0]->lock(); assert(t) ;
00697  IndexSet I1(4000) ;
00698  for (l=Level0+1; l<=J ; l++) {
00699    W->HJ.MatTVec(t,&L[l].I,  L[l-1].d,   &L[l-1].I ,BC , l-1 ,1, Buffers[1]) ;
00700    
00701    a=L[l].I.a ;
00702    e=L[l].I.e ;
00703    n=L[l].I.nr ;
00704    d=L[l].d    ;
00705   
00706    for (r=0; r<n ; r++) 
00707      for (i=a[r]; i<=e[r]; i++) if (d[i]==1e133) d[i]=t[i] ;
00708 
00709   }
00710 
00711  Buffers[0]->unlock() ;
00712 }
00713 
00714 
00715 /****************************************/
00716 /*                                      */
00717 /*   Funktionen fuer AdaptiveGS::       */
00718 /*                                      */
00719 /****************************************/
00720 
00721 AdaptiveGS::AdaptiveGS()  { 
00722   Level0=J=0 ; 
00723   BC[0]=BC[1]=-MAXINT ;
00724   XA=0.0 ; XE=1.0 ;
00725   for (unsigned int i=0; i<sizeof(Buffers)/sizeof(Buffers[0]); i++) Buffers[i]=NULL;
00726 }
00727 
00728 AdaptiveGS::~AdaptiveGS() { Level0=J=0 ;}
00729 
00730 int AdaptiveGS::Init(int Lev0,int j) {
00731   if (Level0) {assert ((Level0==Lev0)&&(J==j)) ; return 1 ;}
00732   Level0=Lev0 ; J=j ;
00733 
00734   for (int i=Level0;i<=J;i++) if(L[i].Init((1<<i)+1)==0) return 0  ;
00735   return 1  ;
00736 }
00737 
00738 
00739 void AdaptiveGS::CopyIndexSets(AdaptiveGS *S)
00740 {
00741  int l ;
00742  for (l=Level0;l<=J;l++) {
00743    L[l].I  .Copy(&S->L[l].I  ) ;
00744    L[l].II .Copy(&S->L[l].II ) ; 
00745    L[l].III.Copy(&S->L[l].III) ;
00746    L[l].IV .Copy(&S->L[l].IV ) ;
00747  }
00748 }
00749 
00750 bool AdaptiveGS::SameIndexSetsI(AdaptiveGS *S)
00751 {
00752  int l ;
00753  for (l=Level0; l<=J; l++) if (!L[l].I.IsSame(&S->L[l].I)) return false ;
00754  return true ;
00755 }
00756 bool AdaptiveGS::SameIndexSetsII(AdaptiveGS *S)
00757 {
00758  int l ;
00759  for (l=Level0; l<=J; l++) if (!L[l].II.IsSame(&S->L[l].II)) return false ;
00760  return true ;
00761 }
00762 bool AdaptiveGS::SameIndexSetsIII(AdaptiveGS *S)
00763 {
00764  int l ;
00765  for (l=Level0; l<=J; l++) if (!L[l].III.IsSame(&S->L[l].III)) return false ;
00766  return true ;
00767 }
00768 bool AdaptiveGS::SameIndexSetsIV(AdaptiveGS *S)
00769 {
00770  int l ;
00771  for (l=Level0; l<=J; l++) if (!L[l].IV.IsSame(&S->L[l].IV)) return false ;
00772  return true ;
00773 }
00774 
00775 bool AdaptiveGS::SameIndexSets(AdaptiveGS *S) {
00776  return SameIndexSetsI  (S) && SameIndexSetsII(S) && 
00777         SameIndexSetsIII(S) && SameIndexSetsIV(S) ;
00778 }
00779 
00780 
00781 
00782 
00783 
00784 void AdaptiveGS::Copy(AdaptiveGS *S)
00785 {
00786  int l ;
00787  for (l=Level0;l<=J;l++) {
00788    L[l].I  .Copy(&S->L[l].I  ) ;
00789    L[l].II .Copy(&S->L[l].II ) ; 
00790    L[l].III.Copy(&S->L[l].III) ;
00791    L[l].IV .Copy(&S->L[l].IV ) ;
00792 
00793    L[l].I.VecCopy(L[l].d , S->L[l].d) ;
00794  }
00795 }
00796 
00797 
00798 void AdaptiveGS::Print() { Print(NULL,0) ;}
00799 void AdaptiveGS::Print(int allflag) { Print("GS",allflag) ;}
00800 void AdaptiveGS::Print(char *st,int allflag) {
00801   int c=0 ;
00802   std::cout << "Level0: "<<Level0<<" J: "<<J<<" BC "<<BC[0]<<','<<BC[1]<<'\n';
00803   for (int l=Level0;l<=J;l++) {
00804     std::cout <<"Level "<<l<<'\n'; 
00805     c+=L[l].Print(st,l,allflag) ;
00806    }
00807   std::cout <<"NZE "<<c<<'\n';
00808 }
00809 
00810 void AdaptiveGS::IndexSetForAdaptiveGS(AdaptiveB *B, Wavelets *WC)
00811 {
00812  L[J].I.InducedByWavelets(B->L[J].IS,J,B->BC,WC) ;
00813  L[J].IV.Copy(&L[J].I) ;
00814 
00815  for (int l=J-1;l>=Level0;l--)  {
00816    L[l].IV.InducedByWavelets    (B->L[l].IS,l  ,B->BC,WC) ; // R(GtlT,Tau(l))
00817    L[l].II.InducedByFineScalings(&L[l+1].I ,l+1,B->BC,WC) ; // R(Hl+1,S(l+1))
00818    L[l].I.Union(&L[l].II,&L[l].IV) ;                        // S(l)=..
00819 
00820    L[l].III.IndexSetForTypeIII  (&L[l+1].I,l+1,B->BC,WC) ;
00821  }
00822 }
00823 
00824 void AdaptiveGS::IndexSetForAdaptiveGS2(AdaptiveB *B , Wavelets *WC)
00825 {
00826  int l ;
00827  for (l=Level0;l<=J;l++) L[l].IV.InducedByWavelets(B->L[l].IS,l,B->BC,WC) ;
00828 
00829  L[J].II.nr=0 ;
00830  L[J].I.Union(&L[J].II,&L[J].IV) ;
00831  for (l=J-1;l>=Level0;l--)  {
00832    L[l].II.InducedByFineScalings(&L[l+1].I,l+1,B->BC,WC) ;
00833    L[l].I.Union(&L[l].II,&L[l].IV) ;
00834    L[l].III.nr=1 ;          // jeden Zugriff auf III ordentlich versauern
00835    L[l].III.a[0]=-MAXINT ;
00836    L[l].III.e[0]=-MAXINT ;
00837  }
00838 }
00839 
00840 
00841 void AdaptiveGS::FromAdaptiveB(AdaptiveB *B , Wavelets *W)
00842 {
00843  int l ;
00844 
00845  BC[0]=B->BC[0] , BC[1]=B->BC[1] ; 
00846  XA=B->XA ; XE=B->XE ;
00847 
00848  if (B->islifting) for (l=Level0; l<=J; l++) B->L[l].IS->VecMul(B->L[l].d , sqrt((double)(1<<l))) ;
00849 
00850  // copy Level0 , it is not assumed that L[Level0].I == B->L[Level0].IS
00851  L[Level0].I.VecFill(L[Level0].d,0.0) ;
00852  B->L[Level0].IS->VecCopy(L[Level0].d , B->L[Level0].d) ;
00853  
00854  // compute all other levels from coarse to fine
00855  for (l=Level0+1; l<=J ;l++) {
00856    if (B->islifting) W->Q.MatTVec(L[l-1].d,&L[l-1].I  ,  B->L[l].d,B->L[l].IS , BC , l-1 , +1, Buffers[0]) ;
00857 
00858    W->HJ.MatTVec(L[l].d,&L[l].I  ,    L[l-1].d,   &L[l-1].I ,BC , l-1 ,1, Buffers[0]) ;
00859    W->GJ.MatTVec(L[l].d,&L[l].I  , B->L[l  ].d, B->L[l  ].IS,BC , l-1 ,0, Buffers[0]) ;
00860 
00861    BoundaryCorrection(BC,L[l].d,l) ;
00862  }
00863 
00864 }
00865 
00866 
00867 void AdaptiveGS::FromAdaptiveB1(AdaptiveB *B, Wavelets *W) {
00868 
00869  int l,i, *BCF=B->BC , BCT[2] ;
00870  static IndexSet   I[30],*IP[30] ;
00871  static AdaptiveGS GS    ; // auxiliary GS 
00872 
00873  bool up0,up1,down0,down1 ;
00874 
00875  BCT[0]=BC[0] , BCT[1]=BC[1] ; 
00876  XA=B->XA ; XE=B->XE ;
00877 
00878  //
00879  // check, what to do
00880  down0 = (BCF[0]==-1) && (BCT[0]==0 ) ;
00881  down1 = (BCF[1]==-1) && (BCT[1]==0 ) ;
00882  up0   = (BCF[0]== 0) && (BCT[0]==-1) ;
00883  up1   = (BCF[1]== 0) && (BCT[1]==-1) ;
00884  assert (!((down0&&up1)||(down1&&up0))) ; // on one side up and on the other down is not allowed
00885  assert (!(down0||down1)) ;               // down not allowed
00886  assert (!((BCF[1]==PERIODIC)&&(BCT[1]!=PERIODIC))) ;
00887  assert (!((BCF[1]!=PERIODIC)&&(BCT[1]==PERIODIC))) ;
00888 
00889  //
00890  // first calculation of ordinary generating system with the same boundary
00891  // conditions than the basis
00892  //
00893   
00894  FromAdaptiveB(B,W) ;
00895  if (!(up0||up1)) return  ;
00896 
00897  // 
00898  // modification of boundary scaling functions
00899  //
00900  
00901  // imbedding of gs
00902  Projection(BCT,W) ;
00903 
00904  //   
00905  // first calculation of indexset for boundary wavelets and corresponding adaptive gs
00906  for (i=0;i<30;i++) I[i].Init(10) ;
00907  GS.Init(Level0,J) ;
00908 
00909  for (l=Level0+1;l<=J;l++) {
00910    I[l].ForBoundaryWavelets(down0||up0,down1||up1,B->L[l].IS,l,W) ;
00911    // tricky: replace the original IndexSet by the boundary indexset
00912    IP[l]=B->L[l].IS ; B->L[l].IS=&I[l] ;
00913   }
00914  // calculate boundary adapted gs
00915  GS.IndexSetForAdaptiveGS2(B,W) ;
00916 
00917  // represent the wavelets of all levels by the appropriate scaling functions
00918  for (l=Level0+1; l<=J ;l++)
00919    W->GJ.MatTVec(GS.L[l].d,&GS.L[l].I  ,  B->L[l].d,B->L[l].IS,BCF, l-1 ,1 , Buffers[0]) ;
00920 
00921  // imbedding of only the boundary scaling functions
00922  GS.BC[0]=BCF[0] , GS.BC[1]=BCF[1] ;
00923  GS.Projection(BCT,W) ;
00924 
00925  //
00926  // now calc. contributions of W^i_0 (i>j) to V^j :overwrite GS 
00927  double *t=Buffers[1]->lock() ;assert(t) ;
00928  for (l=J-1;l>=Level0;l--) {
00929    W->HJ.MatVec(&t[0],&GS.L[l].I ,  GS.L[l+1].d,&GS.L[l+1].I ,BCT,l+1,1 , Buffers[0]) ;
00930    GS.L[l].I.VecPlus(L[l].d,&t[0]    ) ;
00931    GS.L[l].I.VecPlus(GS.L[l].d, &t[0]) ;
00932   } 
00933  Buffers[1]->unlock() ;
00934 
00935  //
00936  // undo replacement
00937  for (l=Level0+1;l<=J;l++) B->L[l].IS=IP[l] ;
00938  
00939  BC[0]=BCT[0] , BC[1]=BCT[1] ;
00940 
00941 }
00942 
00943 
00944 
00945 //
00946 // Project performs the projection V->V(_0) 
00947 // if the BC of BCT and this->BC agree nothing is done
00948 // Project fails, if e.g. BCT[0]==-1 but this->BC[0]==0 :: this operation is not supported
00949 void AdaptiveGS::Projection(int *BCT,Wavelets *W) {
00950  int  l,i,n ;
00951  bool down0,down1,up0,up1 ;
00952  bool do0,do1 ; 
00953 
00954  //
00955  // check, what to do
00956  down0 = (BC[0]==-1) && (BCT[0]>=0 ) ;
00957  down1 = (BC[1]==-1) && (BCT[1]>=0 ) ;
00958  up0   = (BC[0]>= 0) && (BCT[0]==-1) ;
00959  up1   = (BC[1]>= 0) && (BCT[1]==-1) ;
00960 
00961  if (!(down0||up0||down1||up1)) return ;
00962 
00963  // projection of gs
00964  for (l=Level0; l<=J; l++) {
00965    i=1<<l ;   
00966 
00967    do0=do1=false ;
00968    if ((n=L[l].I.nr-1)>=0) {
00969      do0 = (L[l].I.a[0]==0) ;
00970      do1 = (L[l].I.e[n]==i) ;
00971     }
00972  
00973    if (!W->InterpoletFlag) {
00974      if (down0&do0) W->OP0[BCT[0]+1].MatVec (L[l].d,L[l].d) ;
00975      if (down1&do1) W->OP1[BCT[1]+1].MatVec (&L[l].d[i+1-W->K1],&L[l].d[i+1-W->K1]) ;
00976 
00977      if (up0&do0)   W->OP0[BC[0]+1].MatTVec (L[l].d,L[l].d) ;
00978      if (up1&do1)   W->OP1[BC[1]+1].MatTVec (&L[l].d[i+1-W->K1],&L[l].d[i+1-W->K1]) ;
00979    } else {
00980    
00981      if (down0&do0) L[l].d[  BCT[0]]=0 ;
00982      if (down1&do1) L[l].d[i-BCT[1]]=0 ;
00983 
00984      if (up0&do0)   W->OP0[BC[0]].MatVec (L[l].d,L[l].d) ;
00985      if (up1&do1)   W->OP1[BC[1]].MatVec (&L[l].d[i+1-W->K1],&L[l].d[i+1-W->K1]) ;
00986    }
00987 
00988  } // next l
00989 
00990  BC[0]=BCT[0] ;
00991  BC[1]=BCT[1] ;
00992 }
00993 
00994 
00995 //
00996 // BoundaryQuadrature(which) performs the polynomial exact (inverse) qudrature
00997 // for only the boundary scaling functions
00998 //  which >0 ==> quadrature
00999 //  which <0 ==> inverse quadrature
01000 //
01001 void AdaptiveGS::BoundaryQuadrature(Wavelets *W, int which) {
01002  int  l,i,n ;
01003  bool do0,do1 ;
01004  Matrix2D *W0,*W1 ;
01005 
01006  if ((BC[1]==PERIODIC) || (which==0)) return ;
01007  
01008  // 
01009  // select quadrature matrices
01010   
01011  if (which>0) {
01012    W0=&W->W.AL[0] ;
01013    W1=&W->W.AR[0] ;
01014  } else {
01015    W0=&W->IW.AL[0] ;
01016    W1=&W->IW.AR[0] ;
01017  }
01018 
01019  //
01020  // quadrature of gs
01021  for (l=Level0; l<=J; l++) {
01022    i=1<<l ;   
01023 
01024    do0=do1=false ;
01025    if ((n=L[l].I.nr-1)>=0) {
01026      do0 = (L[l].I.a[0]==0) ;
01027      do1 = (L[l].I.e[n]==i) ;
01028     }
01029  
01030    if (do0) W0->MatVec (L[l].d,L[l].d) ;
01031    if (do1) W1->MatVec (&L[l].d[i+1-W->K1],&L[l].d[i+1-W->K1]) ;
01032   }
01033 }
01034
Generated at Mon Aug 19 10:02:31 2002 for AWFD by