Adaptive1DOperations.cc

00001 //: Adaptive1DOperations.cc
00002 
00003 #include "Adaptive1D.hpp"
00004 #include "NonAdaptive.hpp"
00005 #include "Buffer.hpp"
00006 
00007 extern int debugRefine ;
00008 AdaptiveBOperator::AdaptiveBOperator(double aa0,double aa1,Wavelets *W) {
00009  a0=aa0 ;
00010  a1=aa1 ;
00011  WC=W   ;
00012 }
00013 
00014 void AdaptiveBOperator::apply(AdaptiveB *X,AdaptiveB *R) {
00015  assert(X->SameIndexSets(R)) ;
00016  R->ApplyOp(X,2,WC) ;
00017  R->Add(a0,X,a1,R)  ;
00018 }
00019 
00020 void AdaptiveBOperator::Preconditioner(AdaptiveB *X,AdaptiveB *R) {
00021  int l,i,r,*a,*e,n ;
00022  double *d,*dr,q   ;
00023  for (l=X->Level0; l<=X->J; l++) {
00024    q=a1*(1<<(2*l))+a0 ;
00025    a=X->L[l].IS->a  ;
00026    e=X->L[l].IS->e  ;
00027    n=X->L[l].IS->nr ; 
00028    d=X->L[l].d ;
00029    dr=R->L[l].d ;
00030    for (r=0; r<n; r++) for (i=a[r]; i<=e[r]; i++) dr[i]=d[i]/q ;
00031  }
00032 }
00033 
00034 // 
00035 // ApplyOp(int op,int *BCT,int *BCF, Wavelets *W) berechnet Galerkin-Projektion 
00036 // eines Diff-Operators bz.w. Matrix-Vektor Product mit SteifigkeitsMatrix
00037 // op==1 1. Ableitung
00038 // op==2 SteifigkeitsMatrix
00039 //
00040 
00041 //extern int debugRefine ;
00042 
00043 void AdaptiveB::ApplyOp(AdaptiveB *X, int op, Wavelets *W)
00044 { static AdaptiveGS G1, G2 ;
00045   G1.Init(Level0,J)    ;
00046   G1.IndexSetForAdaptiveGSOp(X,W) ;
00047   G2.Init(Level0,J)    ;
00048   ApplyOp(X,op,W,&G1,&G2) ;
00049 }
00050 
00051 //
00052 // Calculate IndexSet for adaptive GS required for evaluation of operators
00053 // 
00054 void AdaptiveGS::IndexSetForAdaptiveGSOp(AdaptiveB *X, Wavelets *W) 
00055 {
00056  int l , *BCF=X->BC ;
00057 
00058  for (l=Level0;l<=J;l++)
00059    L[l].III.InducedByWavelets( X->L[l].IS , l, BCF , W) ;  
00060    
00061  for (l=J; l>=Level0; l--) {
00062    if (l==J) {
00063      L[l].IV.InducedByOperatorMatrix(&L[l].III , l , BCF , &W->Operators[0][0][0]) ;
00064    } else {
00065      L[l].II.InducedByOperatorMatrix(&L[l].III , l , BCF , &W->Operators[0][0][0]) ;
00066      L[l].I.InducedByFineScalings  (&L[l+1].IV ,l+1, BCF , W) ;
00067      L[l].IV.Union(&L[l].II , &L[l].I) ;
00068    }
00069  }
00070 
00071  //
00072 
00073  L[J-1].II.InducedByOperatorMatrix( X->L[J].IS, J-1, BCF, &W->Operators[0][0][1]) ;
00074  for (l=J-1; l>Level0; l--) {
00075    L[l-1].III.InducedByOperatorMatrix    ( X->L[l].IS , l-1, BCF, &W->Operators[0][0][1]) ;
00076    L[l-1].I  .InducedByFineDualScalings  (&L[l].II , l  , BCF, W) ;
00077 
00078    if (l>Level0+1) L[l-1].II.Union(&L[l-1].III , &L[l-1].I ) ;
00079               else L[l-1].I .Union(&L[l-1].III , &L[l-1].II) ;
00080  }
00081  l=Level0 ;
00082  L[l].II.Union(&L[l].I , X->L[l].IS) ;
00083 
00084  for (l=Level0; l<=J; l++) L[l].I.Union(&L[l].II , &L[l].IV) ;
00085 }
00086 
00087 
00088 void AdaptiveB::ApplyOp(AdaptiveB *X, int op, Wavelets *W, AdaptiveGS *G1, AdaptiveGS *G2)
00089 {
00090  int     l,op0=NOOPERATION,opcode=-MAXINT ;
00091  int    *BCT=BC , *BCF=X->BC ;
00092  double  DX=XE-XA ; 
00093 
00094  BoundaryCorrection(BCF,X->L[Level0].d,Level0) ;
00095 
00096  // finite difference Operators
00097  if (W->N==2) {
00098    if (op==FIRSTDERIVATIVE) ApplyFD(X, FD12 , W, G1) ;
00099    if (op==STIFFNESS      ) ApplyFD(X, FD22 , W, G1) ;
00100    return ;
00101   }
00102  
00103  if ((op==FD12)||(op==FD14)||(op==FD16)||(op==FD22)||(op==FD24)||
00104      (op==GRADFD4)||(op==DIVFD4)||(op==GRADFD6)||(op==DIVFD6)) { 
00105    ApplyFD(X, op , W, G1) ; 
00106    return ;
00107  }
00108 
00109  // 
00110  if (op==DIV ) {op=FIRSTDERIVATIVE ; op0=DIV  ;}
00111  if (op==GRAD) {op=FIRSTDERIVATIVE ; op0=GRAD ;}
00112 
00113  if (op==FIRSTDERIVATIVE) opcode=0 ;
00114  if (op==STIFFNESS      ) opcode=1 ;
00115 
00116  assert(X!=this) ;
00117  CopyIndexSets(X) ;
00118  
00119  // compute b(.,.) via generating system
00120  l=Level0 ;
00121  W->Operators[opcode][0][0].MatVec(L[l].d, L[l].IS,BCT , X->L[l].d, X->L[l].IS,BCF ,l,1 ,DX, Buffers[0]) ;
00122  
00123  G1->L[Level0].I.VecFill  (G1->L[l].d, 0.0) ;
00124   X->L[Level0].IS->VecCopy(G1->L[l].d, X->L[l].d) ;
00125  
00126  for (l=Level0+1; l<=J; l++) {
00127    W->Operators[opcode][1][1].MatVec(L[l].d, L[l].IS,BCT,  X->L[l  ].d,   X->L[l  ].IS, BCF, l-1, 1 , DX , Buffers[0]) ;
00128    W->Operators[opcode][1][0].MatVec(L[l].d, L[l].IS,BCT, G1->L[l-1].d, &G1->L[l-1].I , BCF, l-1, 0 , DX , Buffers[0]) ;
00129 
00130    if (l<J) {
00131      W->HJ.MatTVec(G1->L[l].d, &G1->L[l].I, G1->L[l-1].d, &G1->L[l-1].I , BCF, l-1, 1 , Buffers[0]) ;
00132      W->GJ.MatTVec(G1->L[l].d, &G1->L[l].I,  X->L[l  ].d,   X->L[l  ].IS, BCF, l-1, 0 , Buffers[0]) ;
00133    }
00134  }
00135 
00136 
00137  // compute a(j,.) 
00138  assert(J>Level0) ;
00139 
00140  double *t=Buffers[0]->lock() ;assert(t) ;
00141 
00142  W->Operators[opcode][0][1].MatVec(t,&G1->L[J-1].I, BCT  ,  X->L[J].d, X->L[J].IS,BCF, J-1,1 , DX , Buffers[1]) ;
00143  for (l=J-1 ; l>Level0; l--) {
00144    W->GtJ.MatVec(    L[l  ].d,      L[l  ].IS, t, &G1->L[l].I, BCT, l, 0 , Buffers[1]) ;
00145    W->HtJ.MatVec(G2->L[l-1].d, &G1->L[l-1].I , t, &G1->L[l].I, BCT, l, 1 , Buffers[1]) ;
00146    W->Operators[opcode][0][1].MatVec(t, &G1->L[l-1].I, BCT, X->L[l].d, X->L[l].IS, BCF, l-1, 1, DX , Buffers[1]);
00147    
00148    G1->L[l-1].I.VecPlus(t, G2->L[l-1].d) ;
00149  }
00150 
00151  G1->L[Level0].I.VecPlus(L[Level0].d,t) ;
00152  Buffers[0]->unlock() ;
00153 
00154  //
00155  // DivGrad
00156  // 
00157  if (op0!=NOOPERATION) {
00158    assert(X->BC[0]==-1) ;
00159    assert((X->BC[1]==-1)||(X->BC[1]==PERIODIC)) ;
00160    assert((op0==GRAD)||(op0==DIV)) ;
00161 
00162    int FD=-1 ;
00163    double alpha=1e+100 ;
00164    if (W->N==6) { alpha=300.0; FD = (op0==GRAD) ? FD801 : FD80 ;}
00165    if (W->N==4) { alpha= 15.0; FD = (op0==GRAD) ? FD601 : FD60 ;}
00166 
00167    G1->FromAdaptiveB(X,W) ;
00168    G1->ApplyFD(BCT , G1 , W , FD) ;
00169 
00170    G1->IndexSetForAdaptiveGS(X,W) ;
00171    X->FromAdaptiveGS3(G1 , W) ;
00172    if (op0==GRAD) PPlus(+alpha, X) ;
00173           else    PPlus(-alpha, X) ;
00174  }
00175 
00176  BoundaryCorrection(BCT,L[Level0].d,Level0) ;
00177 }
00178 
00179 //
00180 //  Projection via the adaptive GS
00181 //
00182 void AdaptiveB::Projection(AdaptiveB *X, Wavelets *W, AdaptiveGS *G1, AdaptiveGS *G2) {  
00183  IndexSet *IP[30] ; 
00184  int    l   ;
00185  int    *BCT=BC , *BCF=X->BC ;
00186  bool   down0,down1,up0,up1 ;
00187 
00188  if (X==this) return ;
00189 
00190  // copy basis
00191  CopyIndexSets(X) ;
00192  for (l=Level0; l<=J; l++) L[l].IS->VecCopy(L[l].d , X->L[l].d) ; 
00193 
00194  //
00195  // check, what to do
00196  down0 = (BCF[0]==-1) && (BCT[0]>=0 ) ;
00197  down1 = (BCF[1]==-1) && (BCT[1]>=0 ) ;
00198  up0   = (BCF[0]>= 0) && (BCT[0]==-1) ;
00199  up1   = (BCF[1]>= 0) && (BCT[1]==-1) ;
00200 
00201  assert (!((down0&&up1)||(down1&&up0))) ; // on one side up and on the other down is not allowed
00202 
00203  if (!(down0||up0||down1||up1)) return ;  // copy was just enough
00204 
00205  //
00206  // first calculation of IndexSet for boundary wavelets
00207 
00208  for (l=Level0+1; l<=J; l++) {
00209    // dirty hack: for piecewise linear interpolets we build the complete genertaing system, 
00210    // since KW0=KW1=0 leads to problems in the optimized code
00211 
00212    if (W->InterpoletFlag && (W->N==2)) {
00213      IP[l]=X->L[l].IS ;
00214    } else {
00215      // 
00216      // general case: mark just the boundary wavelets und use these only to build the generating system
00217      // the required index set is stored in GS2.L[l].I which is not used for other purposes
00218      //
00219      G2->L[l].I.ForBoundaryWavelets(down0||up0,down1||up1,X->L[l].IS,l,W) ;
00220      // tricky: replace the original IndexSet by the boundary indexset
00221      IP[l]=X->L[l].IS       ; // will be used to restore the original Indexsets 
00222      X->L[l].IS=&G2->L[l].I ;
00223    }
00224   }
00225 
00226  // calculate boundary adapted gs
00227  G1->IndexSetForAdaptiveGS2(X,W) ;
00228  
00229  //
00230  // projection of V->V_0
00231  if (down0|down1) {
00232 
00233    // calc. generating system 
00234    G1->BC[0]=BCF[0] ;
00235    G1->BC[1]=BCF[1] ;
00236    G1->FromAdaptiveB(X,W) ;
00237 
00238    // projection of gs  
00239    G1->Projection(BCT,W) ;
00240  } 
00241   
00242  //
00243  // change of basis
00244  if (up0|up1) {
00245 
00246    //
00247    // calculate imbedding W^j_0 in V^j :stored in GS
00248 
00249    // first copy Level0
00250    G1->L[Level0].I.VecClear (G1->L[Level0].d) ;
00251     X->L[Level0].IS->VecCopy(G1->L[Level0].d, X->L[Level0].d) ;
00252 
00253    // represent the wavelets of all levels by the appropriate scaling functions
00254    for (l=Level0+1; l<=J ;l++) 
00255       W->GJ.MatTVec(G1->L[l].d, &G1->L[l].I, X->L[l].d,X->L[l].IS,BCF, l-1 ,1, Buffers[0]) ;
00256 
00257    // change basis for gs
00258    G1->BC[0]=BCF[0] ;
00259    G1->BC[1]=BCF[1] ;   
00260    G1->Projection(BCT,W) ;
00261 
00262    //
00263    // now calc. contributions of W^i_0 (i>j) to V^j :stored in GS2 
00264    //
00265 
00266    double *t=Buffers[0]->lock(); assert(t) ;
00267 
00268    G1->L[J].I.VecClear(G2->L[J].d) ;
00269    for (l=J-1;l>=Level0;l--) {
00270      G1->L[l+1].I.VecAdd( t , G2->L[l+1].d , G1->L[l+1].d) ;
00271      W->HtJ.MatVec(G2->L[l].d,&G1->L[l].I ,  t, &G1->L[l+1].I ,BCT,l+1,1, Buffers[1]) ;
00272     }
00273 
00274    Buffers[0]->unlock(); 
00275 
00276    //
00277    // now calc. all contributions of W^i_0 (i<=j) to V_j  :overwrite GS
00278    for (l=Level0+1; l<=J; l++)
00279      W->HJ.MatTVec (G1->L[l].d,&G1->L[l].I  , G1->L[l-1].d, &G1->L[l-1].I, BCT , l-1 ,0, Buffers[0]) ;
00280 
00281    // finally sum up both contributions :overwrite GS
00282    for (l=Level0; l<=J; l++)
00283      G1->L[l].I.VecPlus(G1->L[l].d , G2->L[l].d) ;
00284  }
00285 
00286  // 
00287  // Projection of the generating system to wavelet spaces
00288  for (l=J; l>Level0; l--) 
00289    W->GtJ.MatVec(L[l].d, X->L[l].IS , G1->L[l].d ,&G1->L[l].I,BCT , l,1, Buffers[0]) ;
00290    // W->GtJ.MatVec(L[l].d, &I[l] , GS.L[l].d ,&GS.L[l].I,BCT , l,1) ; // alte Version 7.2. 2000
00291 
00292  l=Level0 ;
00293  L[Level0].IS->VecCopy(L[Level0].d,G1->L[l].d) ;
00294  //BoundaryCorrection(BCT,L[l].d,l) ;
00295 
00296  //
00297  // restore complete old basis
00298  for (l=Level0+1; l<=J; l++) X->L[l].IS=IP[l] ;
00299 
00300 } 
00301 
00302 
00303 //
00304 //
00305 void AdaptiveB::RestrictPressure() {
00306  int    l,*a,*e,*a1,*e1,n,n1,i,r,r1,cnt ;
00307  double *d ;
00308  
00309  // clear each 20th coefficient c[l,i] , which satisfies 
00310  //  (*) c[l+1,2*i] and c[l+1,2*i+1] are NOT active
00311  // 
00312  for (l=Level0+1; l<J; l++) {
00313    n =L[l  ].IS->nr ;
00314    a =L[l  ].IS->a ;
00315    e =L[l  ].IS->e ;
00316    d =L[l  ].d ;
00317 
00318    n1=L[l+1].IS->nr ;
00319    a1=L[l+1].IS->a;
00320    e1=L[l+1].IS->e;
00321 
00322    cnt=r1=0 ;
00323    for (r=0; r<n; r++) 
00324      for (i=a[r]; i<=e[r]; i++) {
00325        int i21=2*i+1 ;
00326        // find r1 such that (r1>=n1)   or  (2*i+1 <= L[l+1].e[r1])
00327        while ((r1<n1) && (i21)>e1[r1]) r1++ ;
00328        if ((r1>=n1) || (2*i<a1[r1])) { 
00329          if ( (cnt%20)== 0) d[i]=0 ;
00330          cnt++ ;
00331         }
00332       }
00333  }
00334 
00335  // clear finest coefficients in any case
00336  n=L[J].IS->nr ;
00337  a=L[J].IS->a ;
00338  e=L[J].IS->e ;
00339  d=L[J].d  ;
00340  for (r=0; r<n; r++) 
00341    for (i=a[r]; i<=e[r]; i++) d[i]=0 ;
00342  
00343 }
00344 
00345 // 
00346 // FD Operators
00347 
00348 //
00349 //
00350 void AdaptiveGS::ApplyFD(int *BCT , AdaptiveGS *G , Wavelets *W , int op) {
00351   int l ;
00352 
00353   double DX=XE-XA, *f, *t=Buffers[0]->lock(); assert(t) ;
00354   
00355   for (l=Level0; l<=J; l++) {
00356     if (this!=G) f=G->L[l].d ;
00357            else { 
00358              G->L[l].I.VecCopy(t,G->L[l].d) ; 
00359              f=t ;
00360            }
00361     switch (op) {
00362       case FD12: case FD14: case FD16: case FD22: case FD24:  case FD40: case FD60: case FD80: 
00363     case GRADFD4: case DIVFD4: case FD11://case GRADFD6: case DIVFD6:
00364         W->FDOperators[W->FDOp(op)].MatVec(L[l].d,&L[l].I,BCT  ,  f,&G->L[l].I,G->BC, l , 1 , DX , Buffers[1]) ;
00365        break ;
00366       case FD401: {
00367          W->FDOperators[W->FDOp(FD40)].MatVec(L[l].d,&L[l].I,BCT  ,  f,&G->L[l].I,G->BC, l , 1 , DX, Buffers[1] ) ;
00368          if (BCT[1]!=PERIODIC) L[l].d[0]=L[l].d[1<<l]=0 ;
00369         }
00370        break ;
00371       case FD601: {
00372          W->FDOperators[W->FDOp(FD60)].MatVec(L[l].d,&L[l].I,BCT  ,  f,&G->L[l].I,G->BC, l , 1 , DX, Buffers[1]) ;
00373          if (BCT[1]!=PERIODIC) L[l].d[0]=L[l].d[1<<l]=0 ;
00374         }
00375        break ;
00376       case FD801: {
00377          W->FDOperators[W->FDOp(FD80)].MatVec(L[l].d,&L[l].I,BCT  ,  f,&G->L[l].I,G->BC, l , 1 , DX, Buffers[1]) ;
00378          if (BCT[1]!=PERIODIC) L[l].d[0]=L[l].d[1<<l]=0 ;
00379         }
00380        break ;
00381      default: assert(0) ;
00382     }
00383   }
00384 
00385   // BC[0]=BCT[0]; BC[1]=BCT[1] ; ??
00386 
00387  Buffers[0]->unlock() ;
00388 }
00389 
00390 
00391 void AdaptiveGS::ApplyWENO(int *  , AdaptiveGS *G, AdaptiveGS *GRoeSpeed, int op) {
00392   int l ;
00393   
00394   double DX=XE-XA, *f, *t=Buffers[0]->lock(); assert(t) ;
00395 
00396   assert(GRoeSpeed != this) ;
00397   assert(GRoeSpeed != G   ) ;
00398 
00399   for (l=Level0; l<=J; l++) {
00400     if (this!=G) f=G->L[l].d ;
00401            else {
00402              G->L[l].I.VecCopy(t,G->L[l].d) ; 
00403              f=t ;
00404            }
00405     switch (op) {
00406       case WENO3: G->L[l].I.VecApplyWENO(L[l].d, f, GRoeSpeed->L[l].d, BC, l, 2 , DX) ;
00407        break ;
00408       case WENO5: G->L[l].I.VecApplyWENO(L[l].d, f, GRoeSpeed->L[l].d, BC, l, 3 , DX) ;
00409        break ;
00410       default:    assert(0) ;
00411        break ;
00412     }
00413   }
00414   Buffers[0]->unlock() ;
00415 }
00416 
00417 //
00418 //
00419 void AdaptiveGS::IndexSetForAdaptiveGSFD(AdaptiveB *X , OperatorMatrix *FD , Wavelets *W) {
00420  int l , *BCF=X->BC   ;
00421  struct MatrixShape H ;
00422  
00423  assert ((BCF[0]==-1)&& ((BCF[1]==-1)||(BCF[1]==PERIODIC))) ;
00424 
00425  //
00426  // calc index sets: IV is the index set required for the wavelet transform in the final step of this alg.
00427 
00428  for (l=J; l>=Level0; l--) {
00429    L[l].IV .InducedByWavelets(X->L[l].IS, l , BCF , W) ;
00430 
00431    FD->GetShape(l,&H) ;
00432    H.Transpose()      ;// Attention U,O are not interchanged, thus works only if U==O
00433    if (l==J)
00434      L[l].I.CAM(&L[l].IV, BCF , &H ,
00435                 OperatorMatrixFirstNZEP , OperatorMatrixLastNZEP ,
00436                 OperatorMatrixFirstNZE  , OperatorMatrixLastNZE ) ;
00437    else {
00438      L[l].III.CAM(&L[l].IV, BCF , &H ,  
00439                   OperatorMatrixFirstNZEP , OperatorMatrixLastNZEP ,  
00440                   OperatorMatrixFirstNZE  , OperatorMatrixLastNZE ) ;
00441 
00442      L[l].II.InducedByFineScalings   (&L[l+1].I,l+1, BCF , W)   ;
00443      L[l].I .Union                   (&L[l].II , &L[l].III) ;
00444    }
00445  }
00446 
00447  // now, subindex sets where result of FD matrix-vector product is exact
00448  // ???
00449  for (l=Level0; l<=J; l++) {
00450    FD->GetShape(l,&H) ;
00451    L[l].II.EAM(&L[l].I , BCF , &H , 
00452                OperatorMatrixFirst2NZEP , OperatorMatrixLast2NZEP ,  
00453                OperatorMatrixFirst2NZE  , OperatorMatrixLast2NZE ) ;
00454   }
00455 
00456  // now, type III scaling functions for Delta-correction, 
00457  for (l=Level0; l<J; l++) {
00458    W->HtJ.GetShape(l+1,&H) ;
00459    L[l].III.EAM(&L[l+1].II , BCF , &H ,
00460                 TransformMatrixFirst2NZEP , TransformMatrixLast2NZEP ,
00461                 TransformMatrixFirst2NZE  , TransformMatrixLast2NZE ) ;
00462   }
00463 }
00464 
00465 //
00466 //
00467 void AdaptiveB::ApplyFD(AdaptiveB *X, int op, Wavelets *W, AdaptiveGS *GS)
00468 {
00469  BoundaryCorrection(X->BC,X->L[Level0].d,Level0) ;
00470 
00471  GS->FromAdaptiveB(X , W) ;
00472  GS->ApplyFD(X->BC , GS , W , op) ;
00473  FromAdaptiveGS(GS , W) ;
00474 
00475  BoundaryCorrection(BC,L[Level0].d,Level0) ;
00476 }
00477 
00478 
00479 
00480 double AFD2Multiply(int i,double *d,double *fd,IndexSet *I) {
00481 
00482   bool ok = I->Contains(i-3) && I->Contains(i-1) && I->Contains(i-0) && I->Contains(i+1) && I->Contains(i+3) ;
00483   if (!ok) { std::cout<< "AFD2 "<<i<<'\n'; I->Print(); printf("%d\n",1/(i-1)); assert(0); }
00484 
00485   return fd[0]*d[i-3] + fd[2]*d[i-1] +fd[3]*d[i] + fd[4]*d[i+1] + fd[6]*d[i+3] ;
00486 }
00487 
00488 
00489 double AFD2MultiplyP(int i,double *d,double *fd,IndexSet *I, int l2) {
00490   int l1=l2-1, im3=(i-3+l2)&l1, im1=(i-1+l2)&l1, ip1=(i+1+l2)&l1, ip3=(i+3+l2)&l1 ; 
00491   
00492   bool ok = I->Contains(im3) && I->Contains(im1) && I->Contains(i) && I->Contains(ip1) && I->Contains(ip3) ;
00493   if (!ok) { std::cout<< "AFD2 "<<i<<'\n'; I->Print(); printf("%d\n",1/(i-1)); assert(0); }
00494   
00495   return fd[0]*d[im3] + fd[2]*d[im1] +fd[3]*d[i] + fd[4]*d[ip1] + fd[6]*d[ip3] ;
00496 }
00497 
00498 
00499 void AdaptiveB::ApplyFD2(AdaptiveB *X, int op, Wavelets *W, AdaptiveGS *GS)
00500 {
00501  int  l,*a,*e,r,i,nr,*a2,*e2,K,l2,l1,op2,ra,re,nr2,NK2,N2,j ;
00502  double fd24[7]={-1./72, 0, 81./72, -160./72, 81./72,0 , -1./72},
00503         fd14[7]={ 1./48, 0,-27./48,        0, 27./48,0 , -1./48},*fd,q,p, *d,*d3, DX=XE-XA ;
00504 
00505  double *f=Buffers[1]->lock(); 
00506  assert(f) ;
00507 
00508  //bool dbg=false ;// (X->L[3].IS->nr==1) && (X->L[4].IS->nr==1) && (X->L[5].IS->nr==0) && (X->L[3].IS->a[0]==0) && (X->L[3].IS->e[0]==4) && (X->L[4].IS->a[0]==0) && (X->L[4].IS->e[0]==2) ; 
00509 
00510  switch (op) {
00511  case FD24II: op2=FD24, K=4 ; p=2; fd=fd24; break ; 
00512  case FD14II: op2=FD14, K=4 ; p=1; fd=fd14; break ; 
00513  default: {
00514    std::cout <<"AdaptiveB::ApplyFD2 : operation "<<op<<"not yet supported\n" ;
00515    exit(-1) ;
00516   }
00517  }
00518 
00519  NK2=W->N-2+K/2 ;
00520  N2 =W->N-2     ;
00521 
00522  BoundaryCorrection(X->BC,X->L[Level0].d,Level0) ;
00523 
00524  // compute the generating system
00525  //
00526  // patch for Level0: since not all points of the coarsest level are true gridpoints it may happen
00527  // that some of our default stencils can not be evaluated. (The rationale is rather technical and involved.) 
00528  // Therefore, we virtually fill up the coarsest level. To this end, the values in the 
00529  // non-active gridpoints are set zero.
00530  //  
00531  for (i=0; i<=(1<<Level0); i++) GS->L[Level0].d[i]=0 ;
00532  GS->IndexSetForAdaptiveGS(X , W) ;
00533  GS->FromAdaptiveB(X , W) ;
00534 
00535  // virtual index set for Level0
00536  IndexSet I0(2) ;
00537  I0.a[0]=0 ;
00538  I0.e[0]=(BC[1]==PERIODIC) ? (1<<Level0)-1 : (1<<Level0) ;
00539  I0.nr  =1 ;
00540  
00541  for (l=Level0; l<=J; l++) {
00542    l2=1<<l ;
00543    l1=l2-1 ;
00544    q =pow(l2,p) ;
00545 
00546    // non-periodic BCs: left/right bound of first/last interval are 0 or l2 ?
00547    bool closea0=false, closee0=false, closean=false, closeen=false ; 
00548    // periodic BCs: 
00549    bool complete=false, periodicinterval=false, onedisappears=false ;
00550  
00551    // calculate index set \subset S for which standard FD stencils may be used
00552    // Except the case that the interval boundaries collapse with the boundaries of [0,...,2^l]
00553    // some special FD stencils must be used at the boundaries. Therefore we seperate the points where
00554    // a special treatment is necessary. 
00555    a =GS->L[l].IV.a; e =GS->L[l].IV.e; nr=GS->L[l].IV.nr; // S
00556    a2=GS->L[l].II.a; e2=GS->L[l].II.e;                    // subset of S 
00557    if (l==Level0) { a =I0.a; e =I0.e; nr=I0.nr; } 
00558 
00559    // nothing to do in this level means nothing to do also for all higher levels
00560    if (nr==0) break ; 
00561 
00562    ra=1; re=nr-2;
00563    if (BC[1]!=PERIODIC) {
00564      if (a[0]==0) a2[0]=0, closea0=true ;
00565      else         a2[0]=a[0]+NK2 ;
00566 
00567      if (e[0]==l2) e2[0]=l2, closee0=true ;
00568      else          e2[0]=e[0]-NK2 ;
00569 
00570      if (a[nr-1]==0) a2[nr-1]=0, closean=true ;
00571      else            a2[nr-1]=a[nr-1]+NK2 ;
00572 
00573      if (e[nr-1]==l2) e2[nr-1]=l2, closeen=true ;
00574      else             e2[nr-1]=e[nr-1]-NK2 ;
00575  
00576      nr2=nr ;
00577 
00578    } else { // periodic case
00579 
00580      if ((a[0]==0) && (e[nr-1]>=l2-1)) { // a real periodic interval! this is the complicated case
00581 
00582        if (nr==1) {  // just one big interval
00583          a2[0]=0; e2[0]=l2-1; nr2=1 ;
00584          complete=true ;
00585          re=0 ;      // do no other interval calculations later on
00586          goto later;
00587        }
00588      
00589        if (e[0]-NK2<0) { // the first interval will disappear
00590          e2[nr-2]=e[0   ]-NK2+l2 ;
00591          a2[nr-2]=a[nr-1]+NK2    ;
00592          for (r=1; r<nr-1; r++) { a2[r-1]=a[r]+NK2; e2[r-1]=e[r]-NK2; }
00593          nr2=nr-1 ;
00594          periodicinterval=onedisappears=true ;
00595          re=0 ;
00596          goto later ;
00597        }
00598 
00599        if (a[nr-1]+NK2>=l2) { // last interval will disapear
00600          a2[0]=a[nr-1]+NK2-l2 ;
00601          e2[0]=e[0]-NK2       ;
00602          nr2=nr-1  ;
00603          periodicinterval=onedisappears=true ;
00604          goto later ; 
00605        }
00606      
00607        // last and first interval survive the restriction operation
00608        a2[0   ]=0           ;e2[0   ]=e[0]-NK2 ;
00609        a2[nr-1]=a[nr-1]+NK2 ;e2[nr-1]=l2-1     ;
00610        nr2=nr    ;
00611        periodicinterval=true ;
00612        goto later ;
00613 
00614      } else { // no periodic interval !
00615        ra=0; re=nr-1     ;
00616        nr2=nr ;
00617      }
00618    } //
00619 
00620  later:
00621    GS->L[l].II.nr=nr2 ;
00622    for (r=ra; r<=re; r++) {a2[r]=a[r]+NK2; e2[r]=e[r]-NK2; }   
00623 
00624    // test 
00625    for (r=0; r<nr2; r++) {
00626     assert(a2[r]<=e2[r]); 
00627     assert(a2[r]>=0) ;
00628     assert(e2[r]>=0) ;
00629    }
00630 
00631    // apply standard FD stencil
00632    d =GS->L[l  ].d ;
00633    d3=GS->L[l-1].d ;
00634 
00635    //if (debugRefine) 
00636    //GS->L[l].IV.VecFill(f,1e+100) ; // debug
00637    W->FDOperators[W->FDOp(op2)].MatVec(f,&GS->L[l].II,BC , d,&GS->L[l].IV,BC  , l , 1 , DX , Buffers[0]) ;
00638 
00639    // apply special f.d. stencil at the boundaries of inner intervals
00640    if (BC[1]!=PERIODIC) {
00641      if (!closea0) { 
00642        i=a2[0]-1; f[i]=AFD2Multiply(i,d,fd,&GS->L[l].IV)*q ;
00643        for (i=a[0]; i<=a[0]+N2; i +=2) f[i]=d3[i/2] ;
00644      }
00645 
00646      if (!closee0) { 
00647        i=e2[0]+1; f[i]=AFD2Multiply(i,d,fd,&GS->L[l].IV)*q ;
00648        for (i=e[0]-N2; i<=e[0]; i +=2) f[i]=d3[i/2] ;
00649      }
00650 
00651      if (!closean) { 
00652        i=a2[nr-1]-1; f[i]=AFD2Multiply(i,d,fd,&GS->L[l].IV)*q ;
00653        for (i=a[nr-1]; i<=a[nr-1]+N2; i +=2) f[i]=d3[i/2] ;
00654      }
00655      if (!closeen) { 
00656        i=e2[nr-1]+1; f[i]=AFD2Multiply(i,d,fd,&GS->L[l].IV)*q ;
00657        for (i=e[nr-1]-N2; i<=e[nr-1]; i +=2) f[i]=d3[i/2] ;
00658      }
00659      for (r=1; r<nr-1; r++) { 
00660        i=a2[r]-1; f[i]=AFD2Multiply(i,d,fd,&GS->L[l].IV)*q ;
00661        i=e2[r]+1; f[i]=AFD2Multiply(i,d,fd,&GS->L[l].IV)*q ;
00662        for (i=a[r]   ; i<=a[r]+N2; i +=2) f[i]=d3[i/2] ;
00663        for (i=e[r]-N2; i<=e[r]   ; i +=2) f[i]=d3[i/2] ;
00664      }
00665   
00666    } else { // periodic 
00667 
00668      if (!complete) {
00669 
00670        //
00671        // apply boundary stencils
00672        //
00673        ra=0; re=nr2-1 ; // by default all intervals are inner intervals
00674        if ((periodicinterval==true) && (onedisappears==false)) { // last and first need special treatment
00675          i=(e2[0   ]+1   )&l1; f[i]=AFD2MultiplyP(i,d,fd,&GS->L[l].IV,l2)*q ;
00676          i=(a2[nr-1]-1+l2)&l1; f[i]=AFD2MultiplyP(i,d,fd,&GS->L[l].IV,l2)*q ;
00677          ra=1; re=nr-2 ; // standard treatment for inner intervals only
00678        }
00679        // inner intervals 
00680        for (r=ra; r<=re; r++) {
00681          i=(a2[r]-1+l2)&l1; f[i]=AFD2MultiplyP(i,d,fd,&GS->L[l].IV,l2)*q ;
00682          i=(e2[r]+1+l2)&l1; f[i]=AFD2MultiplyP(i,d,fd,&GS->L[l].IV,l2)*q ;
00683        }
00684        
00685        //
00686        // copy from lower levels
00687        //
00688        ra=0; re=nr-1 ;
00689        if (periodicinterval==true) { // last and first intervall need special treatment
00690          for (i=e[0   ]-N2; i<=e[0   ]   ; i +=2) { j=(i+l2)&l1; f[j]=d3[j/2] ; }
00691          for (i=a[nr-1]   ; i<=a[nr-1]+N2; i +=2) { j=(i   )&l1; f[j]=d3[j/2] ; }
00692          ra=1; re=nr-2 ;
00693        }
00694 
00695        for (r=ra; r<=re; r++) {
00696          for (i=e[r]-N2; i<=e[r]   ; i +=2) { j=(i+l2)&l1; f[j]=d3[j/2] ; }
00697          for (i=a[r]   ; i<=a[r]+N2; i +=2) { j=(i   )&l1; f[j]=d3[j/2] ; }
00698        }
00699 
00700      } // !complete
00701 
00702    }
00703 
00704    GS->L[l].IV.VecCopy(d,f) ;
00705 
00706  } // l
00707  Buffers[1]->unlock() ;
00708 
00709  // get Basis representation of result using I, III ,IV
00710  FromAdaptiveGS(GS , W) ;
00711  BoundaryCorrection(BC,L[Level0].d,Level0) ;
00712 }
00713 
00714 
00715 
00716 //
00717 //
00718 void AdaptiveB::ApplyWENO(AdaptiveB *X, int op, Wavelets *W, AdaptiveGS *GS, AdaptiveGS *GRoeSpeed)
00719 {
00720  BoundaryCorrection(X->BC,X->L[Level0].d,Level0) ;
00721 
00722  // Now, get generating system
00723  GS->FromAdaptiveB(X , W) ;
00724  
00725  // Apply Operator
00726  GS->ApplyWENO(X->BC, GS, GRoeSpeed, op) ;
00727  
00728  // get Basis representation of result using I, III ,IV
00729  FromAdaptiveGS(GS , W) ;
00730 
00731  BoundaryCorrection(BC,L[Level0].d,Level0) ;
00732 }
00733 
00734 //
00735 //  Multiplikations-Routine ueber Erzeugenden-System
00736 //
00737 
00738 void AdaptiveB::Multiply(AdaptiveB *u , AdaptiveB *v, Wavelets *W,int which)
00739 {
00740   /* static */  AdaptiveGS U,V ;
00741  int l , *uBC=u->BC ;  
00742  
00743  assert( (which==ONE_POINT_MULTIPLY) || (which==BOUND_MULTIPLY) ) ;
00744  assert(u->SameIndexSets(v)) ;
00745  assert(BC[0]==v->BC[0])  ;
00746  assert(BC[1]==v->BC[1])  ;
00747   
00748  CopyIndexSets(u) ;
00749  BC[0]=uBC[0] ;
00750  BC[1]=uBC[1] ;
00751 
00752  U.Init(Level0,J) ; // wird tatsachlich nur das aller erstemal aufgerufen
00753  V.Init(Level0,J) ;
00754 
00755  // Erzeugenden-System aufstellen
00756 
00757  U.IndexSetForAdaptiveGS(u,W) ;
00758  V.CopyIndexSets(&U) ;
00759 
00760  U.FromAdaptiveB(u ,W) ;
00761  V.FromAdaptiveB(v ,W) ;
00762 
00763  //
00764  // switch to nodal values
00765  if (which==BOUND_MULTIPLY) {
00766    U.BoundaryQuadrature(W,-1) ; 
00767    V.BoundaryQuadrature(W,-1) ;
00768  }
00769 
00770  //
00771  // Multipliy of generating systems
00772  for (l=Level0;l<=J;l++) {
00773    double q=1<<l ;
00774    double *du =U.L[l].d    ;
00775    double *dv =V.L[l].d    ;
00776    int    *a  =U.L[l].I.a  ;
00777    int    *e  =U.L[l].I.e  ;
00778    int     n  =U.L[l].I.nr ;
00779    int     r,i ;
00780 
00781    q=sqrt(q) ;
00782 
00783    for (r=0;r<n;r++) for (i=a[r];i<=e[r];i++) du[i]*=(dv[i]*q) ;
00784     
00785    if ((BC[1]!=PERIODIC)&&(which==ONE_POINT_MULTIPLY)) {
00786      int s,nJ=1<<l ;
00787      Matrix2D *A   ;   
00788   
00789      A=&W->ATau0[BC[0]+1] ;
00790      s=A->size[1] ;
00791      for (i=0;i<s;i++) du[i]*= A->a[0][i] ;
00792 
00793      A=&W->ATau1[BC[1]+1] ; 
00794      s=A->size[1] ;
00795      for (i=0;i<s;i++)  du[nJ-s+1+i]*= A->a[0][i] ;
00796    } 
00797  }
00798 
00799  // back to scaling function coefficients
00800  if ((BC[1]!=PERIODIC)&&(which==BOUND_MULTIPLY)) U.BoundaryQuadrature(W,1) ;
00801 
00802  // back to wavelet coefficients
00803  FromAdaptiveGS2(&U ,W) ;
00804 }
Generated at Mon Aug 19 10:02:31 2002 for AWFD by