Extensions.hpp

00001 //
00002 // data structure for extension information 
00003 // used by AdaptiveData<DIM> as well as UniformData<DIM> and LevelAdaptiveGrid.hpp
00004 //
00005 
00006 #ifndef _EXTENSIONS_
00007 # define _EXTENSIONS_
00008 
00009 # define _EXTENSIONS_MAXDIM_ 100
00010 
00011 # include "Math.hpp"
00012 
00019 template<int DIM>
00020 struct Extensions {
00021 public: 
00023          Extensions() ;
00024         ~Extensions() ;
00025  
00027   void   SetBoundaryConditions(int bc[DIM][2])  ;
00029   void   SetBoundaryConditions(int dir,int *bc) ;
00031   void   SetBoundaryConditions(int dir,int bc0,int bc1) ;
00033   bool   SameBoundaryConditions(Extensions<DIM> *E) ;
00035   bool   IsSame(Extensions<DIM> *E) ;
00037   void   CopyFlags(Extensions<DIM> *E) ;
00039   void   CopyData (Extensions<DIM> *E) ;
00044   void   Init(Wavelets *W, double *xa=NULL, double *xe=NULL) ;
00046   void   InitXAE() ; 
00047 
00049   void   Print() ;
00050 
00051   void   WriteUDF(FILE *f) ;
00052 
00053 
00054 // data
00056   bool   ismultiscale[DIM] ;
00058   bool   islifting[DIM]    ;
00059 
00060   Wavelets *WC ;
00061  
00062   //double XA[DIM],XE[DIM] ; // grid corresponds to box 
00063                            // default: XA_i=0, XE_i=1
00065   double *XA ;
00067   double *XE ;
00068   bool   ismyXAE ;         // true if XA/XE were allocated using InitXAE
00070   int    BC[DIM][2] ;      // boundary conditions
00071 
00075   int    dimext   ;
00077   double ext[_EXTENSIONS_MAXDIM_] ;
00079   double Max()    ;
00081   double Min()    ;
00083   double MaxAbs() ;
00085   double Sum()    ;
00087   double SumAbs()    ;
00089   double SumSqr() ;
00091   double InnerProd(Extensions<DIM> *X) ;
00092 
00097   void   PPlus(const double a)       ; 
00099   void   PPlus(Extensions<DIM> *X  ) ;
00101   void   PPlus(const double a , Extensions<DIM> *X)   ;
00103   void   Add(Extensions<DIM> *X , Extensions<DIM> *Y) ;
00105   void   Add(const double a , Extensions<DIM> *X , const double b , Extensions<DIM> *Y) ;
00107   void   Add(const double a , Extensions<DIM> *X , const double b , Extensions<DIM> *Y , const double c , Extensions<DIM> *Z) ;
00109   void   Add(const double a , Extensions<DIM> *X , const double b , Extensions<DIM> *Y , const double c , Extensions<DIM> *Z , 
00110              const double d , Extensions<DIM> *Q) ;
00112   void   Add(const double a , Extensions<DIM> *X , const double b , Extensions<DIM> *Y , const double c , Extensions<DIM> *Z ,
00113              const double d , Extensions<DIM> *Q , const double f , Extensions<DIM> *R) ;
00114 
00116   void   MMinus(Extensions<DIM> *X) ;
00118   void   MMinus(const double a, Extensions<DIM> *X)   ;
00120   void   Sub(Extensions<DIM> *X , Extensions<DIM> *Y) ;
00122   void   Mul(Extensions<DIM> *X , Extensions<DIM> *Y) ;
00124   void   Mul(Extensions<DIM> *X , Extensions<DIM> *Y , double f) ;
00126   void   Mul(Extensions<DIM> *X , double f) ;
00128   void   Mul(double f) ;
00129 
00131   void   XplusYmalZ (Extensions<DIM> *X , Extensions<DIM> *Y , Extensions<DIM> *Z) ;
00133   void   XminusYmalZ(Extensions<DIM> *X , Extensions<DIM> *Y , Extensions<DIM> *Z) ;
00135   void   Div(Extensions<DIM> *X , Extensions<DIM> *Y) ;
00137   void   Abs(Extensions<DIM> *X) ;
00139   void   Sqr(Extensions<DIM> *X) ;
00141   void   Pow(Extensions<DIM> *X, double p) ;
00143   void   Set(const double a) ;
00144 } ;
00145 
00146 
00147 template<int DIM>
00148 Extensions<DIM>::Extensions() {
00149   int i ;
00150   for (i=0; i<DIM; i++) {
00151     ismultiscale[i]=true  ;
00152     islifting   [i]=false ;
00153     BC[i][0]=BC[i][1]=-MAXINT ; // no correct initialization causes lots of sigsegv's  
00154   }
00155   
00156   WC=NULL ;
00157   ismyXAE=false ;
00158   XA=XE=NULL    ;
00159   
00160   dimext=0 ;
00161   for (i=0; i<_EXTENSIONS_MAXDIM_; i++) ext[i]=0 ; 
00162 }
00163 
00164 template<int DIM>
00165 Extensions<DIM>::~Extensions() {
00166   if (ismyXAE) {
00167     delete []XA; 
00168     delete []XE ;
00169   }
00170 }
00171 
00172 template<int DIM>
00173 void Extensions<DIM>::Init(Wavelets *W, double *A, double *E) {
00174   WC=W ;
00175   if (A!=NULL) {
00176     assert(E!=NULL) ;
00177     XA=A, XE=E ; 
00178   } else InitXAE() ;
00179 }
00180 
00181 template<int DIM>
00182 void Extensions<DIM>::InitXAE() {
00183   if (XA!=NULL) {
00184     assert(XE) ;
00185   } else {
00186     ismyXAE=true ;
00187     XA=new double[DIM];
00188     XE=new double[DIM];
00189     for (int i=0; i<DIM; i++) { XA[i]=0.0 ; XE[i]=1.0 ; }
00190   }
00191 }
00192 
00193 template<int DIM>
00194 void Extensions<DIM>::SetBoundaryConditions(int bc[DIM][2]) {
00195   int i ;
00196   for (i=0; i<DIM; i++) {
00197     BC[i][0]=bc[i][0] ;
00198     BC[i][1]=bc[i][1] ;
00199   }
00200 }
00201 
00202 template<int DIM>
00203 void Extensions<DIM>::SetBoundaryConditions(int dir,int *bc) {
00204   BC[dir][0]=bc[0] ;
00205   BC[dir][1]=bc[1] ;
00206 }
00207 
00208 template<int DIM>
00209 void Extensions<DIM>::SetBoundaryConditions(int dir,int bc0,int bc1) {
00210   BC[dir][0]=bc0 ;
00211   BC[dir][1]=bc1 ;
00212 }
00213 
00214 template<int DIM>
00215 bool Extensions<DIM>::SameBoundaryConditions(Extensions<DIM> *E) {
00216   int i ;
00217   for (i=0; i<DIM; i++) 
00218     if ( (BC[i][0] != E->BC[i][0]) || (BC[i][1] != E->BC[i][1])) return false ;
00219  
00220   return true ;
00221 }
00222 
00223 template<int DIM>
00224 void Extensions<DIM>::CopyFlags(Extensions<DIM> *E) {
00225   int i ;
00226   dimext=E->dimext ;
00227   WC    =E->WC     ;
00228   if (!ismyXAE) {
00229     XA    =E->XA ;
00230     XE    =E->XE ;
00231   }
00232  
00233   for (i=0; i<DIM; i++) {
00234     islifting   [i]=E->islifting[i]    ;
00235     ismultiscale[i]=E->ismultiscale[i] ;
00236     
00237     BC[i][0]=E->BC[i][0] ;
00238     BC[i][1]=E->BC[i][1] ;
00239    
00240     if (ismyXAE) {
00241       XA[i]=E->XA[i] ;
00242       XE[i]=E->XE[i] ;
00243     }
00244   }
00245 
00246   //
00247   // ext =E->ext ; 
00248   // this line of code would be WRONG --> dont copy ANY numerical data like ext !!!
00249   //
00250 }
00251 
00252 template<int DIM>
00253 void Extensions<DIM>::CopyData(Extensions<DIM> *E) {
00254  for (int i=0; i<dimext; i++) ext[i]=E->ext[i] ;
00255 }
00256 
00257 
00258 template<int DIM>
00259 bool Extensions<DIM>::IsSame(Extensions<DIM> *E) {
00260   int i ;
00261   if (dimext!=E->dimext) goto print ;
00262   if (WC    !=E->WC)     goto print ;
00263   
00264   for (i=0; i<DIM; i++) {
00265     if (islifting   [i]!=E->islifting[i])    goto print ;
00266     if (ismultiscale[i]!=E->ismultiscale[i]) goto print ;
00267     
00268     if (BC[i][0]!=E->BC[i][0]) goto print ;
00269     if (BC[i][1]!=E->BC[i][1]) goto print ;
00270 
00271     if (XA[i]!=E->XA[i]) goto print ;
00272     if (XE[i]!=E->XE[i]) goto print ;
00273   }
00274 
00275   return true ;
00276 
00277  print:
00278   E->Print() ;
00279   Print() ;
00280   return false ;
00281 }
00282 
00283 template<int DIM>
00284 void Extensions<DIM>::Print() {
00285   std::cout<<"WC="<<WC<<" dimext="<<dimext<<'\n' ;
00286   std::cout<<"i: lift ismult |  BC[0] BC[1] |  XA XE\n" ;
00287   for (int i=0; i<DIM; i++)
00288     std::cout<<i<<":  "<<islifting[i]<<"   "<<ismultiscale[i]<<"  |   "<<BC[i][0]<<' '<<BC[i][1]<<"  |  "<<XA[i]<<' '<<XE[i]<<'\n' ;
00289 }
00290 
00291 
00292 template<int DIM>
00293 void Extensions<DIM>::WriteUDF(FILE *f) {
00294   int  i ;
00295   unsigned int si,st1=0, st2=0;
00296   
00297   for (i=0; i<DIM; i++) {
00298     si=1<<i ;
00299     if (islifting   [i]) st1 |= si ;
00300     if (ismultiscale[i]) st2 |= si ;
00301   }
00302   
00303   fprintf(f," %d %d  %d %d ",WC->N,WC->TypeID(WC->type),st1,st2) ;
00304 
00305   for (i=0; i<DIM; i++) fprintf(f," %d %d ",BC[i][0],BC[i][1]) ;
00306 
00307   fprintf(f,"\n") ;
00308 }
00309 
00310 
00311 // algebraic functions on ext[]
00312 template<int DIM>
00313 double Extensions<DIM>::Max() {
00314   double r=-HUGE_VAL ;
00315   for (int i=0; i<dimext; i++)  r=max(r,ext[i]) ;
00316   return r ;
00317 }
00318 
00319 template<int DIM>
00320 double Extensions<DIM>::Min() {
00321   double r=+HUGE_VAL ;
00322   for (int i=0; i<dimext; i++)  r=min(r,ext[i]) ;
00323   return r ;
00324 }
00325 
00326 template<int DIM>
00327 double Extensions<DIM>::MaxAbs() {
00328   double r=-HUGE_VAL ;
00329   for (int i=0; i<dimext; i++)  r=max(r,fabs(ext[i])) ;
00330   return r ;
00331 }
00332 
00333 template<int DIM>
00334 double Extensions<DIM>::Sum() {
00335   double r=0 ;
00336   for (int i=0; i<dimext; i++)  r += ext[i] ;
00337   return r ;
00338 }
00339 
00340 template<int DIM>
00341 double Extensions<DIM>::SumAbs() {
00342   double r=0 ;
00343   for (int i=0; i<dimext; i++)  r += fabs(ext[i]) ;
00344   return r ;
00345 }
00346 
00347 template<int DIM>
00348 double Extensions<DIM>::SumSqr() {
00349   double r=0 ;
00350   for (int i=0; i<dimext; i++)  r += ext[i]*ext[i] ;
00351   return r ;
00352 }
00353 
00354 
00355 template<int DIM>
00356 double Extensions<DIM>::InnerProd(Extensions<DIM> *E) {
00357   assert(dimext==E->dimext) ;
00358   double r=0 ;
00359   for (int i=0; i<dimext; i++)  r += ext[i] * E->ext[i] ;
00360   return r ;
00361 }
00362 
00363 // PPlus / Add / MMinus / Sub /Mul /Div / Abs /Pow
00364 template<int DIM>
00365 void Extensions<DIM>::PPlus(const double a) {
00366   for (int i=0; i<dimext; i++) ext[i] += a ;
00367 }
00368 template<int DIM>
00369 void Extensions<DIM>::PPlus(Extensions<DIM> *X) {
00370   assert(dimext==X->dimext) ;
00371   for (int i=0; i<dimext; i++) ext[i] += X->ext[i] ;
00372 }
00373 template<int DIM>
00374 void Extensions<DIM>::PPlus(const double a , Extensions<DIM> *X) {
00375   assert(dimext==X->dimext) ;
00376   for (int i=0; i<dimext; i++) ext[i] += a*X->ext[i] ; 
00377 }
00378 
00379 // ADD
00380 template<int DIM>
00381 void Extensions<DIM>::Add(Extensions<DIM> *X , Extensions<DIM> *Y) {
00382   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00383   for (int i=0; i<dimext; i++) ext[i] = X->ext[i] + Y->ext[i] ; 
00384 }
00385 
00386 template<int DIM>
00387 void Extensions<DIM>::Add(const double a , Extensions<DIM> *X , 
00388                           const double b , Extensions<DIM> *Y) {
00389   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00390   for (int i=0; i<dimext; i++) ext[i] = a*X->ext[i] + b*Y->ext[i] ; 
00391 }
00392 
00393 template<int DIM>
00394 void Extensions<DIM>::Add(const double a , Extensions<DIM> *X , 
00395                           const double b , Extensions<DIM> *Y , 
00396                           const double c , Extensions<DIM> *Z ) {
00397   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00398   assert(dimext==Z->dimext) ;
00399   for (int i=0; i<dimext; i++) ext[i] = a*X->ext[i] + b*Y->ext[i] + c*Z->ext[i] ; 
00400 }
00401 
00402 template<int DIM>
00403 void Extensions<DIM>::Add(const double a , Extensions<DIM> *X , 
00404                           const double b , Extensions<DIM> *Y , 
00405                           const double c , Extensions<DIM> *Z ,
00406                           const double d , Extensions<DIM> *Q ) {
00407   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00408   assert(dimext==Z->dimext) ; assert(dimext==Q->dimext) ;
00409   for (int i=0; i<dimext; i++) 
00410     ext[i] = a*X->ext[i] + b*Y->ext[i] + c*Z->ext[i] + d*Q->ext[i] ;
00411 }
00412 
00413 template<int DIM>
00414 void Extensions<DIM>::Add(const double a , Extensions<DIM> *X , 
00415                           const double b , Extensions<DIM> *Y , 
00416                           const double c , Extensions<DIM> *Z ,
00417                           const double d , Extensions<DIM> *Q ,
00418                           const double e , Extensions<DIM> *R ) {
00419   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00420   assert(dimext==Z->dimext) ; assert(dimext==Q->dimext) ;
00421   assert(dimext==R->dimext) ;
00422   for (int i=0; i<dimext; i++) 
00423     ext[i] = a*X->ext[i] + b*Y->ext[i] + c*Z->ext[i] + d*Q->ext[i] + e*R->ext[i] ;
00424 }
00425 
00426 // MMINUS
00427 template<int DIM>
00428 void Extensions<DIM>::MMinus(Extensions<DIM> *X) {
00429   assert(dimext==X->dimext) ;
00430   for (int i=0; i<dimext; i++) ext[i] -= X->ext[i] ;
00431 }
00432 template<int DIM>
00433 void Extensions<DIM>::MMinus(const double a , Extensions<DIM> *X) {
00434   assert(dimext==X->dimext) ;
00435   for (int i=0; i<dimext; i++) ext[i] -= a*X->ext[i] ; 
00436 }
00437 // SUB
00438 template<int DIM>
00439 void Extensions<DIM>::Sub(Extensions<DIM> *X , Extensions<DIM> *Y) {
00440   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00441   for (int i=0; i<dimext; i++) ext[i] = X->ext[i] - Y->ext[i] ; 
00442 }
00443 // MUL
00444 template<int DIM>
00445 void Extensions<DIM>::Mul(Extensions<DIM> *X , Extensions<DIM> *Y) {
00446   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00447   for (int i=0; i<dimext; i++) ext[i] = X->ext[i] * Y->ext[i] ; 
00448 }
00449 template<int DIM>
00450 void Extensions<DIM>::Mul(Extensions<DIM> *X , Extensions<DIM> *Y, double f) {
00451   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00452   for (int i=0; i<dimext; i++) ext[i] = X->ext[i] * Y->ext[i] * f ; 
00453 }
00454 template<int DIM>
00455 void Extensions<DIM>::Mul(Extensions<DIM> *X , double f) {
00456   assert(dimext==X->dimext) ;
00457   for (int i=0; i<dimext; i++) ext[i] *= X->ext[i] * f ; 
00458 }
00459 template<int DIM>
00460 void Extensions<DIM>::Mul(double f) {
00461   for (int i=0; i<dimext; i++) ext[i] *= f ; 
00462 }
00463 // DIV
00464 template<int DIM>
00465 void Extensions<DIM>::Div(Extensions<DIM> *X , Extensions<DIM> *Y) {
00466   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00467   for (int i=0; i<dimext; i++) ext[i] = X->ext[i] / Y->ext[i] ; 
00468 }
00469 // XplusYmalZ
00470 template<int DIM>
00471 void Extensions<DIM>::XplusYmalZ(Extensions<DIM> *X , Extensions<DIM> *Y , 
00472                                  Extensions<DIM> *Z) {
00473   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00474   assert(dimext==Z->dimext) ;
00475   for (int i=0; i<dimext; i++) ext[i] = X->ext[i] + (Y->ext[i]*Z->ext[i]) ;
00476 }
00477 template<int DIM>
00478 void Extensions<DIM>::XminusYmalZ(Extensions<DIM> *X , Extensions<DIM> *Y , 
00479                                   Extensions<DIM> *Z) {
00480   assert(dimext==X->dimext) ; assert(dimext==Y->dimext) ;
00481   assert(dimext==Z->dimext) ;
00482   for (int i=0; i<dimext; i++) ext[i] = X->ext[i] - (Y->ext[i]*Z->ext[i]) ;
00483 }
00484 // ABS/POW/Set
00485 template<int DIM>
00486 void Extensions<DIM>::Abs(Extensions<DIM> *X) {
00487   assert(dimext==X->dimext) ;
00488   for (int i=0; i<dimext; i++) ext[i] = fabs(X->ext[i]) ;
00489 }
00490 
00491 template<int DIM>
00492 void Extensions<DIM>::Sqr(Extensions<DIM> *X) {
00493   assert(dimext==X->dimext) ;
00494   for (int i=0; i<dimext; i++) ext[i] = sqr(X->ext[i]) ;
00495 }
00496 
00497 template<int DIM>
00498 void Extensions<DIM>::Pow(Extensions<DIM> *X , double p) {
00499   assert(dimext==X->dimext) ;
00500   for (int i=0; i<dimext; i++) ext[i] = pow(X->ext[i], p) ;
00501 }
00502 template<int DIM>
00503 void Extensions<DIM>::Set(const double p) {
00504   for (int i=0; i<dimext; i++) ext[i] = p ;
00505 }
00506 
00507 #endif
Generated at Mon Aug 19 10:02:31 2002 for AWFD by