Refine.cc

00001 #include "Refine.hpp"
00002 #include "Adaptive1D.hpp"
00003 #include "Buffer.hpp"
00004 
00005 AdaptivityCriterion::AdaptivityCriterion()                             { SetStrategy(0,0)   ;}
00006 AdaptivityCriterion::AdaptivityCriterion(const int s)                  { SetStrategy(s,0)   ;}
00007 AdaptivityCriterion::AdaptivityCriterion(const int s,const double eps) { SetStrategy(s,eps) ;}
00008 
00009 void AdaptivityCriterion::SetStrategy(const int s,const double e) {
00010  strategy=s    ;
00011  eps     =e    ;
00012  MaxLevel=LMAX ;
00013  q1  =0 ;
00014  qinf=0 ;
00015 }
00016 
00017 void AdaptivityCriterion::SetStrategy(const int    s) {SetStrategy(s,0) ;} 
00018 void AdaptivityCriterion::SetEps     (const double e) {eps=e ;}
00019 void AdaptivityCriterion::Print() {
00020   if (strategy==L2_THRESHOLD) std::cout<<"strategy=L2_THRESHOLD\n" ; 
00021   if (strategy==H1_THRESHOLD) std::cout<<"strategy=H1_THRESHOLD\n" ; 
00022   if (strategy==Hq_THRESHOLD) std::cout<<"strategy=Hq_THRESHOLD q1="<<q1<<" qinf="<<qinf<<'\n' ;
00023   if (strategy==SPARSE_GRID ) std::cout<<"strategy=SPARSE_GRID\n" ; 
00024   if (strategy==ANY         ) std::cout<<"strategy=ANY\n" ;
00025   std::cout<<"MaxLevel="<<MaxLevel<<" eps="<<eps<<'\n' ;  
00026 }
00027 
00028 
00029 bool AdaptivityCriterion::IsActive(const double x,const int *l,const int Level0,const int Dim) {
00030   double y=fabs(x),w=0 ;
00031   bool  l0=true    ;
00032   
00033   int i,l1=0,linf=0 ;
00034   for (i=0; i<Dim; i++) {
00035     l1   += l[i] ;
00036     linf  =max(linf,l[i]) ;
00037     if (l[i]>Level0) l0=false ;
00038   }
00039  
00040   if (linf>MaxLevel) return false ;
00041   if (l0           ) return true  ;
00042   
00043   if (strategy==L2_THRESHOLD) return (y>=eps) ;
00044   if (strategy==H1_THRESHOLD) { for (i=0;i<Dim;i++) w += 1 << l[i] ;
00045                                 return y*w>=eps ;}
00046 
00047   if (strategy==Hq_THRESHOLD) { for (i=0;i<Dim;i++) w += pow(2.0,l[i]*qinf) ;
00048                                 return y*pow(2.0,l1*q1)*w>=eps ;}
00049      
00050   if (strategy==SPARSE_GRID) return (l1   <= eps) ;
00051   if (strategy==FIXED_LEVEL) return (linf <= eps) ;
00052   if (strategy==ANY        ) return true ;
00053   return true ;
00054 }
00055 
00056 extern int debugRefine ;
00057 void IndexSet::SatisfyBoundCondition(int M,int  ,int BL0, int BR0,bool *Insert0, bool *Insert1) {
00058 int r ;
00059 //
00060 // left boundary
00061 
00062   if ((nr>0)&&(a[0]<=BL0)) *Insert0=true ;
00063 
00064   if (*Insert0) {
00065     if (nr==0) { a[0]=0; e[0]=BL0; nr=1 ; }                 // insert one interval
00066     else { 
00067       if (a[0]>0) {
00068         if (a[0]<=BL0+1) { a[0]=0 ;e[0]=max(e[0],BL0) ; }   // expand first interval
00069         else {                                              //  insert new interval
00070           for (r=nr; r>0; r--) { a[r]=a[r-1]; e[r]=e[r-1]; }
00071           a[0]=0; e[0]=BL0;
00072           nr++ ;
00073         }
00074       } else { // a[0]==0 ;
00075         e[0]=max(e[0],BL0) ;
00076         if ((nr>1)&&(e[0]+1>=a[1])) {
00077           e[1]=max(e[0],e[1]) ;
00078           for (r=1;r<nr;r++) { a[r-1]=a[r]; e[r-1]=e[r] ;}
00079           nr-- ;
00080           a[0]=0 ;
00081         }     
00082 
00083       }
00084     }   // if (n>0)
00085   } 
00086 
00087   //
00088   // right boundary
00089 
00090   if ((nr>0)&&(e[nr-1]>=BR0)) *Insert1=true ;
00091 
00092   if (*Insert1) {
00093     if (nr==0) { a[0]=BR0; e[0]=M; nr=1 ; }                            // insert one interval
00094     else { 
00095       if (e[nr-1]<M) {
00096         if (e[nr-1]>=BR0-1) { a[nr-1]=min(a[nr-1],BR0); e[nr-1]=M ; }  // expand last interval
00097         else {                                                         // insert new interval
00098           a[nr]=BR0 ; e[nr]=M ;
00099           nr++ ;
00100          }
00101       } else {
00102         a[nr-1]=min(a[nr-1],BR0) ;
00103         if ((nr>1)&&(e[nr-2]+1>=a[nr-1])) { e[nr-2]=M ;nr-- ;}
00104       }
00105      }   // if (n>0)
00106    } 
00107 }
00108 
00109 
00110 //
00111 // checks whether the given adaptive basis satisfies the bound-cone-condition
00112 //
00113 
00114 bool AdaptiveB::CheckBoundCondition(Wavelets *W) {
00115  if (BC[1]==PERIODIC) return true ;
00116   
00117  bool Insert0=false, Insert1=false ;
00118  int  l,K0,K1,si=-MAXINT,*a,*e,n ;
00119 
00120  K0=W->KW0 ; K1=W->KW1 ;
00121  for (l=J;l>Level0;l--) {
00122    if (l >Level0) { K0=W->KW0 ; K1=W->KW1 ; si=1<<(l-1) ; }
00123    if (l==Level0) { K0=W->K0  ; K1=W->K1  ; si=(1<<l)+1 ; }
00124    a=L[l].IS->a ; e=L[l].IS->e ; n=L[l].IS->nr ;
00125 
00126    if (n>0) {
00127      if (Insert0||(a[0  ]<K0    )) assert((a[0  ]==0   )&&(e[0  ]+1>=K0)) ;
00128      if (Insert1||(e[n-1]>=si-K1)) assert((e[n-1]==si-1)&&(a[n-1]<=si-K1)) ;
00129     
00130      Insert0 = Insert0 || (a[0]<K0) ;
00131      Insert1 = Insert1 || (e[n-1]>=si-K1) ;
00132    }
00133   }
00134 
00135  return true ;
00136 }
00137 
00138 
00139 int LeftP   (int j,MatrixShape *S) { return j-S->U ; }
00140 int RightP  (int j,MatrixShape *S) { return j+S->O ; }
00141 int Left    (int j,MatrixShape *S) { int i=j-S->U ; 
00142                                      if(i<=S->BL0) i=0 ;
00143                                      return i ; }
00144 int Right   (int j,MatrixShape *S) { int i=j+S->O ;
00145                                      if (i>=S->BR0) i=S->M ;
00146                                      return i ; }
00147  
00148 int UpLeftP (int j,MatrixShape *S) { return 2*j-S->U ; }
00149 int UpRightP(int j,MatrixShape *S) { return 2*j+S->O ; }
00150 int UpLeft  (int j,MatrixShape *S) { int i=2*j-S->U ; 
00151                                      if(i<=S->BL0) i=0 ;
00152                                      return i ;  }
00153 int UpRight (int j,MatrixShape *S) { int i=2*j+S->O ;
00154                                      if (i>=S->BR0) i=S->M ;
00155                                      return i ; }
00156 
00157 // Refine 1 dim. case with data copy
00158 void AdaptiveB::Refine(AdaptiveB *B , AdaptivityCriterion *R,Wavelets *W) {
00159   int l[1] ;
00160   Refine(B,R,&l[0],0,1,W,true) ;
00161 }
00162 // Refine 1 dim. case with data copy
00163 void AdaptiveB::Refine(AdaptivityCriterion *R,Wavelets *W) {
00164   int l[1] ;
00165   Refine(this,R,&l[0],0,1,W,true) ; 
00166 }
00167 
00168 //
00169 // Refine refines an instance of AdaptiveB
00170 // since this function is called also in multivariat codes, 
00171 // the dimension and also the multivariate level must be be known 
00172 // 
00173 // general function
00174 
00175 void AdaptiveB::Refine(AdaptiveB *B , AdaptivityCriterion *R,int *l,int dir,int Dim,Wavelets *W,bool copyflag)
00176 {
00177  static IndexSet I[LMAX+1],II[LMAX+1],III ;
00178  int ld,r,i,*a,*e,n,*al,*el,nl,N,M ;
00179  double *d ;
00180  bool Insert0,Insert1,NotPeriodic = (B->BC[1]!=PERIODIC) ;
00181 
00182  // reserve auxiliary indexsets
00183  for (ld=W->Level0; ld<=J; ld++) { 
00184     I[ld].Init(1<<(ld-1)) ;
00185    II[ld].Init(1<<(ld-1)) ;
00186   }
00187  
00188  i=III.Init(1<<(LMAX-1)) ;
00189  assert(i) ;
00190  // 
00191  // coarse on each level:   IS -> I
00192  //
00193  for (ld=J;ld>=Level0;ld--) {
00194 
00195    // get indexset pointers
00196    a=B->L[ld].IS->a ;
00197    e=B->L[ld].IS->e ;
00198    n=B->L[ld].IS->nr ;
00199    d=B->L[ld].d ;
00200 
00201    al=I[ld].a ;
00202    el=I[ld].e ;
00203    l[dir]=ld ;
00204 
00205    nl=-1 ;
00206    for (r=0;r<n;r++) 
00207    for (i=a[r]; i<=e[r]; i++) {
00208      // insert ?
00209      if(R->IsActive(d[i],l,Level0,Dim)) {
00210        //yes. 
00211        //still in an active interval-> set new possible end
00212        int ee=(nl>=0) ? el[nl]+1 :  -10000 ;
00213        if (ee==i) el[nl]=i ;
00214        else {nl++ ; al[nl]=el[nl]=i ; } //new interval
00215      }
00216     }
00217    // finished main loop with active interval-> finish interval now
00218    // set number of intervals
00219    I[ld].nr=nl+1 ;
00220 
00221    // Set 'Active' Flags in old Basis 
00222    if (!copyflag) {
00223      B->L[ld].IS->VecFill(d, 0.0) ;
00224      I[ld].       VecFill(d, 1.0) ;
00225     }
00226   }
00227 
00228  //
00229  // refine on each level:  I -> II
00230  //  
00231  struct MatrixShape S ;
00232  S.U=S.O=AdaptivityCriterion::REFINEBOUNDBOX ;
00233  for (ld=Level0;ld<=J;ld++) {
00234    // get dimensions of current level   
00235    // and set shape of index-inducing matrix
00236    N=(ld==Level0) ? 1<<ld : (1<<(ld-1))-1 ; // max interval [0..N]
00237    S.M=S.N=N ;
00238    if (ld==Level0) {
00239      S.BL1=W->K0-1   ;
00240      S.BR1=N+1-W->K1 ;
00241    } else {
00242      S.BL1=W->KW0-1   ;
00243      S.BR1=N+1-W->KW1 ;
00244    }
00245    S.BL0=S.BL1 ;
00246    S.BR0=S.BR1 ;
00247    
00248    // get refined index set  
00249    II[ld].CAM(&I[ld],BC,&S, LeftP,RightP, Left,Right) ;
00250   }
00251 
00252  //
00253  // refine on higher level: II -> IS
00254  //
00255  Insert0=Insert1=false ;
00256  S.U=S.O=AdaptivityCriterion::REFINEBOUNDBOX ;
00257  
00258  for (ld=min(J-1 , R->MaxLevel-1)+1; ld<J; ld++) L[ld+1].IS->nr=0 ;
00259 
00260  for (ld=min(J-1 , R->MaxLevel-1); ld>Level0; ld--) {
00261    // get dimensions of current level   
00262    // and set shape of index-inducing matrix
00263    M=(1<< ld   )-1 ;
00264    N=(1<<(ld-1))-1 ; // max interval [0..N]
00265    S.M=M ; S.N=N   ; // ACHTUNG: Fuer Randwavelets werden auschliesslich die 
00266    S.BL1=W->KW0-1  ; // Randwavelets auf hoherem Level aktiviert - keine anderen
00267    S.BL0=S.BL1     ;
00268    S.BR1=N+1-W->KW1 ;
00269    S.BR0=M+1-W->KW1 ;
00270    
00271    // get refined index set  
00272    if ((R->strategy!=AdaptivityCriterion::FIXED_LEVEL)&&
00273        (R->strategy!=AdaptivityCriterion::SPARSE_GRID)
00274     /* && (R->strategy!=AdaptivityCriterion::ANY)*/   ) {
00275 
00276      III.CAM(&I[ld],BC,&S,
00277              UpLeftP, //TransformMatrixTransposeFirstNZEP,
00278              UpRightP,//TransformMatrixTransposeLastNZEP,
00279              UpLeft,  //TransformMatrixTransposeFirstNZE,
00280              UpRight);//TransformMatrixTransposeLastNZE) ;
00281    } else III.nr=0 ;
00282   
00283    // complete new level
00284    L[ld+1].IS->Union(&II[ld+1],&III) ;
00285 
00286    // Check boundary condition     
00287    if (NotPeriodic)  L[ld+1].IS->SatisfyBoundCondition(M,N,S.BL0,S.BR0,&Insert0,&Insert1) ;
00288 
00289  }   // next ld
00290 
00291 
00292  //
00293  // refinement from first to second level requires special treatment
00294  //
00295  ld=Level0 ;
00296  M=(1<<ld)-1   ; // max interval [0..M]
00297  N=(1<<ld)     ;
00298  S.M=M ; S.N=N ; 
00299  S.BL1=W->K0-1 ;
00300  S.BL0=W->KW0-1   ;
00301  S.BR1=N+1-W->K1  ;
00302  S.BR0=M+1-W->KW1 ;
00303 
00304  // get refined index set  
00305  if ((R->strategy!=AdaptivityCriterion::FIXED_LEVEL)&&
00306      (R->strategy!=AdaptivityCriterion::SPARSE_GRID)) {
00307    III.CAM(&I[ld],BC,&S,
00308            LeftP,  //OperatorMatrixFirstNZEP,
00309            RightP, //OperatorMatrixLastNZEP,
00310            Left,   //OperatorMatrixFirstNZE,
00311            Right); //OperatorMatrixLastNZE) ;
00312    } else III.nr=0 ;
00313 
00314  // complete new level
00315  L[ld+1].IS->Union(&II[ld+1],&III) ;
00316 
00317  // Check boundary condition     
00318  if (NotPeriodic)  L[ld+1].IS->SatisfyBoundCondition(M,N,S.BL0,S.BR0,&Insert0,&Insert1) ; 
00319  
00320  // treat Level0
00321  L[ld].IS->Copy(&II[ld]) ;
00322  M=(1<<ld)     ; // max interval [0..M]
00323  N=(1<<ld)     ;
00324 
00325  S.BL1=W->K0-1 ; // Randwavelets auf hoherem Level aktiviert - keine anderen
00326  S.BL0=S.BL1   ;
00327  S.BR1=N+1-W->K1  ;
00328  S.BR0=S.BR1      ;  
00329 
00330  if (NotPeriodic)  L[ld].IS->SatisfyBoundCondition(M,N,S.BL0,S.BR0,&Insert0,&Insert1) ;
00331 
00332  // cone + boundcondition
00333  if (W->InterpoletFlag) InterpoletRegularity(W) ;
00334 
00335  // Copy data  OR   set flags for active entries
00336  if (copyflag) { 
00337    for (ld=Level0; ld<=J; ld++) {
00338      assert(L[ld].IS != B->L[ld].IS) ;
00339      III.Difference(L[ld].IS,B->L[ld].IS) ;
00340      III.VecClear  (L[ld].d) ;
00341      if (this!=B) B->L[ld].IS->VecCopy(L[ld].d , B->L[ld].d) ;
00342     }
00343  }
00344 }
00345 
00346 
00347 
00348 
00349
Generated at Mon Aug 19 10:02:32 2002 for AWFD by