Operations.hpp

00001 //
00002 // ApplyOp applies an univariate operator (first/second derivative) to the given Data
00003 
00004 
00005 //#include "UniformData.hpp"
00006 
00007 #ifndef _OPERATIONS_HPP_
00008 #define _OPERATIONS_HPP_
00009 # include <pthread.h>
00010 
00011 template<int DIM>
00012 void AdaptiveData<DIM>::ApplyOp(AdaptiveData<DIM> *X,int dir,unsigned int op) {
00013   ApplyOp(&Ext.BC[dir][0],X,dir,op) ;
00014 }
00015 
00016 template<int DIM>
00017 void AdaptiveData<DIM>::ApplyOp(int *BCTo, AdaptiveData<DIM> *X,int dir,unsigned int op) {
00018 
00019   int  i ;
00020   assert(G->IsPeriodicValid) ;
00021   CheckAdaptiveGrid(X) ;
00022 
00023   // valid operation
00024   if ((!G->WC->InterpoletFlag) && ((op==TRANSFORM) || (op==INVERSETRANSFORM))) {
00025     std::cout<<"AdaptiveData<DIM>::ApplyOp(.,.,.,{INVERSE}TRANSFORM) supported for Interpolets only!\n"; 
00026     std::cout<<"For other wavelets, these operations lead to valid results only for full grids\n" ; 
00027     assert(0) ;
00028   }
00029 
00030   // patch Operator
00031   if ((G->WC->InterpoletFlag) && (G->WC->N<=4) && (op==GRAD)) op=GRADFD4 ;
00032   if ((G->WC->InterpoletFlag) && (G->WC->N<=4) && (op==DIV )) op=DIVFD4  ;
00033 
00034   // simple case: nothing to do
00035   if ((op==PROJECTION) && (this==X) && (BCTo[0]==X->Ext.BC[dir][0]) && (BCTo[1]==X->Ext.BC[dir][1])) return ;
00036   if ((op==PROJECTION) && (BCTo[0]==X->Ext.BC[dir][0]) && (BCTo[1]==X->Ext.BC[dir][1]))  { Copy(X) ; return ;}
00037 
00038   // copy Flags
00039   Ext.CopyFlags(&X->Ext) ;
00040 
00041   // load boundary conditions
00042   for (i=0; i<_SMP_PROC_; i++) {
00043     G->AB1[i]->BC[0]=X->Ext.BC[dir][0] ;
00044     G->AB1[i]->BC[1]=X->Ext.BC[dir][1] ;
00045     //
00046     G->AB2[i]->BC[0]=BCTo[0] ;
00047     G->AB2[i]->BC[1]=BCTo[1] ;
00048 
00049     G->AB1[i]->islifting=Ext.islifting[dir] ;
00050     G->AB2[i]->islifting=Ext.islifting[dir] ;
00051 
00052     G->AB1[i]->XA=G->AB2[i]->XA=G->AGS[i]->XA=G->XA[dir] ;
00053     G->AB1[i]->XE=G->AB2[i]->XE=G->AGS[i]->XE=G->XE[dir] ;
00054   }
00055 
00056   Ext.BC[dir][0]=BCTo[0] ;
00057   Ext.BC[dir][1]=BCTo[1] ;
00058 
00059   if ((op==PROJECTION)&&(this!=X)) ba->Copy(X->ba) ;
00060 
00061   if (op==DIVGRAD) assert((X->Ext.BC[dir][0]==BCTo[0])&&(X->Ext.BC[dir][1]==BCTo[1])) ;
00062 
00063   // init. jobs
00064   G->smpjobcount++ ;
00065   int l[DIM]={0} ;
00066   for (i=0; i< _SMP_PROC_; i++)   
00067     G->SMPjobs[i].set(op, dir, DIM, G->smpalpstart[dir][i] , G->smpnum_alp[dir][i], &X->ba->d[l] , &ba->d[l]) ;
00068 
00069   // start threads by clearing the poll flags
00070   for (i=1; i<_SMP_PROC_; i++) G->SMPjobs[i].clrpoll() ;
00071 
00072   // start first job
00073   ProcessSMPjob((void *)&G->SMPjobs[0]) ;
00074     
00075   // set flags
00076   switch (op) {
00077     case PROJECTION: case FIRSTDERIVATIVE: case STIFFNESS: case DIV:  case GRAD: 
00078     case FD12: case FD22: case FD24: case FD14: case FD16: case GRADFD4: case DIVFD4: case RESTRICTPRESSURE: 
00079     case DIVGRAD: case DIVGRADNOSTAB: case DIVGRADFD4: case DIVGRADFD4NOSTAB: 
00080     case FD12II: case FD22II: case FD14II: case FD24II: case FD11: {
00081       assert(X->Ext.ismultiscale[dir]) ;
00082       assert(!X->Ext.islifting[dir])   ;
00083       Ext.ismultiscale[dir]=true ;
00084       Ext.islifting[dir]=false   ;
00085      }
00086     break ;
00087     case INVERSETRANSFORM: {
00088       assert(X->Ext.ismultiscale[dir]) ;
00089       assert(!X->Ext.islifting[dir])   ;
00090       Ext.ismultiscale[dir]=false ;
00091       Ext.islifting[dir]=false    ;
00092     } 
00093     break ;
00094     case TRANSFORM: {
00095       assert(!X->Ext.ismultiscale[dir]) ;
00096       assert(!X->Ext.islifting[dir])     ;
00097       Ext.ismultiscale[dir]=true ;
00098       Ext.islifting[dir]=false   ;
00099     } 
00100     break ;
00101     case INTERPOLET2LIFTING: {
00102       assert(X->Ext.ismultiscale[dir]) ;
00103       assert(!X->Ext.islifting[dir])   ;
00104       Ext.ismultiscale[dir]=true ;
00105       Ext.islifting[dir]=true    ;
00106     } 
00107     break ;
00108     case LIFTING2INTERPOLET: {
00109       assert(X->Ext.ismultiscale[dir]) ;
00110       assert(X->Ext.islifting[dir])    ;
00111       Ext.ismultiscale[dir]=true ;
00112       Ext.islifting[dir]=false   ;
00113     } 
00114     break ;
00115     default: { assert(0) ;}
00116   }
00117 
00118   // wait for end of all threads
00119   G->SMPWaitForAllThreads() ;
00120 
00121 }  
00122 
00123 
00124 template<int DIM>
00125 void AdaptiveData<DIM>::ApplyWENO(int *BCTo, AdaptiveData<DIM> *X, AdaptiveData<DIM> *RS, int dir,unsigned int op) {
00126 
00127   int  i ;
00128   assert(G->IsPeriodicValid) ;
00129   CheckAdaptiveGrid(X) ;
00130 
00131   // copy Flags
00132   Ext.CopyFlags(&X->Ext) ;
00133   for (i=0;i<DIM;i++) assert(!RS->Ext.ismultiscale[i]) ;
00134 
00135   // load boundary conditions
00136   for (i=0; i<_SMP_PROC_; i++) {
00137     G->AB1[i]->BC[0]=X->Ext.BC[dir][0] ;
00138     G->AB1[i]->BC[1]=X->Ext.BC[dir][1] ;
00139     //
00140     G->AB2[i]->BC[0]=BCTo[0] ;
00141     G->AB2[i]->BC[1]=BCTo[1] ;
00142 
00143     G->AB1[i]->islifting=Ext.islifting[dir] ;
00144     G->AB2[i]->islifting=Ext.islifting[dir] ;
00145 
00146     G->AB1[i]->XA=G->AB2[i]->XA=G->AGS[i]->XA=G->XA[dir] ;
00147     G->AB1[i]->XE=G->AB2[i]->XE=G->AGS[i]->XE=G->XE[dir] ;
00148   }
00149 
00150   Ext.BC[dir][0]=BCTo[0] ;
00151   Ext.BC[dir][1]=BCTo[1] ;
00152 
00153   // init. jobs
00154   int l[DIM]={0} ;
00155   G->smpjobcount++ ;
00156   for (i=0; i< _SMP_PROC_; i++) G->SMPjobs[i].set(op, dir, DIM, G->smpalpstart[dir][i] , G->smpnum_alp[dir][i], &X->ba->d[l] , &ba->d[l], &RS->ba->d[l]) ;
00157 
00158   // start threads
00159   for (i=1; i<_SMP_PROC_; i++) G->SMPjobs[i].clrpoll() ; 
00160 
00161   // start first thread/job
00162   ProcessSMPjobWENO((void *)&G->SMPjobs[0]) ;
00163  
00164   // wait for end of all threads
00165   G->SMPWaitForAllThreads() ;
00166 
00167 }
00168 
00169 
00170 
00171 template<int DIM>
00172 void AdaptiveData<DIM>::SetFunction(Function *F, bool boundaryonly , bool notransform ) {
00173   int    i,l[DIM],t[DIM],l0 ;
00174   double x[DIM] ;
00175   long   j ;
00176   unsigned long k ;
00177 
00178   if (!G->WC->InterpoletFlag) {
00179     std::cout<<"AdaptiveData<"<<DIM<<">::SetFunction: supported for Interpolets only!\n" ;
00180     exit(-1) ;
00181   }
00182 
00183   // set Flags
00184   if (!G->IsPeriodicValid) {
00185     std::cout<<"AdaptiveData<"<<DIM<<">::SetFunction: coordinate directions with periodic boundary conditions must be set, see AdaptiveGrid::SetPeriodicConditions\n" ;
00186     exit(-1) ;
00187   }
00188 
00189   for (i=0; i<DIM; i++) {
00190     if (G->IsPeriodic[i]) SetBoundaryConditions(i,-1,PERIODIC) ;
00191     else SetBoundaryConditions(i,-1,-1) ;
00192     
00193     Ext.islifting   [i]=false ;
00194     Ext.ismultiscale[i]=false ;
00195   }
00196   Ext.dimext=0     ;
00197 
00198   // iterate through all instances of AdaptiveLevels for the zeroth ProjectionGrid
00199   AdaptiveGrid<DIM>::ProjectionGrid          *pg=G->BasisPG[0] ;
00200   AdaptiveGrid<DIM>::ProjectionGrid::iterator pgit             ;
00201   for (pgit=(*pg).begin(); pgit!=(*pg).end(); ++pgit) { 
00202     AdaptiveLevelsIndex *ali=(*pgit).first  ;
00203     AdaptiveLevels      *als=(*pgit).second ;
00204   
00205     for (i=1; i<DIM; i++) {
00206       l[i]=ali->l[i-1] ; // collect level multiindex for access to vector array
00207       k   =ali->t[i-1] ;
00208       if (l[i] > G->Level0[i]) k=2*k+1 ;  
00209       x[i]=Ext.XA[i] + ((Ext.XE[i]-Ext.XA[i])*((double)k))/(1<<l[i]) ; // coordinates of associated gridpoint
00210     }
00211 
00212     for (l0=G->Level0[0]; l0<=LMAX; l0++) {
00213       AdaptiveLevels::AdaptiveL *al=(*als)[l0] ;
00214       AdaptiveLevels::AdaptiveL::iterator alit ;
00215       
00216       l[0]=l0 ;
00217       for (alit=(*al).begin(); alit != (*al).end(); alit++) { 
00218         k=(*alit).first ;
00219         if (l[0] > G->Level0[0]) k=2*k+1 ;
00220         x[0]=Ext.XA[0] + ((Ext.XE[0]-Ext.XA[0])*((double)k))/(1<<l0) ;
00221 
00222         j=(*alit).second ;
00223         if (j<0) j=1-j   ; 
00224         (*ba->d[l])[j]=F->Eval(x) ; // set nodal value
00225       }
00226     }
00227    }   // next pgit
00228   
00229   // transform to wavelet space
00230   if (!notransform) for (i=0; i<DIM; i++) ApplyOp(Ext.BC[i], this, i ,TRANSFORM) ;
00231   else { 
00232     if(boundaryonly) { std::cout<<"AdaptiveData<DIM>::SetFunction boundaryonly==true && notransform==true not supported yet\n"; assert(0) ;}
00233   }
00234 
00235   // boundaryonly ?
00236   if (boundaryonly)
00237     for (pgit=(*pg).begin(); pgit!=(*pg).end(); ++pgit) {
00238       AdaptiveLevelsIndex *ali=(*pgit).first  ;
00239       AdaptiveLevels      *als=(*pgit).second ;
00240      
00241       bool isboundary,isboundary0=false ; 
00242       for (i=1; i<DIM; i++) { 
00243         l[i]= ali->l[i-1] ; 
00244         t[i]= ali->t[i-1] ;
00245         if ( (l[i]==G->Level0[i]) && ((t[i]==0) || (t[i]==(1<<G->Level0[i]))) && (Ext.BC[i][1]!=PERIODIC) )  isboundary0 = true ;
00246        }
00247       
00248       for (l0=G->Level0[0]; l0<=LMAX; l0++) {
00249         AdaptiveLevels::AdaptiveL *al=(*als)[l0] ;
00250         AdaptiveLevels::AdaptiveL::iterator alit ;
00251         
00252         l[0]=l0    ;
00253         for (alit=(*al).begin(); alit != (*al).end(); alit++) { 
00254           t[0]=(*alit).first ;
00255           
00256           isboundary = isboundary0 || ( (l0==G->Level0[0]) && ((t[0]==0)||(t[0]==(1<<G->Level0[0]))) && (Ext.BC[0][1]!=PERIODIC) ) ;
00257             
00258           if (!isboundary) {
00259             j=(*alit).second ;
00260             if (j<0) j=1-j   ;
00261             (*ba->d[l])[j]=0.0 ; // clear wavelet coefficient
00262           }
00263         }
00264       } // ld
00265     }   // pgit
00266 }
00267 
00268 
00269 template<int DIM>
00270 void AdaptiveData<DIM>::SetDirichletBoundaryValues(AdaptiveData<DIM-1> *F[DIM][2] , int BC[DIM][2], bool ContinuousDirichletValues) {
00271   int d,i,k,i2,k2, nl0= 1<<G->WC->Level0 ;
00272 
00273   // initialize
00274   Set(0.0) ;
00275   for (i=0; i<DIM; i++) {
00276     if (G->IsPeriodic[i]) SetBoundaryConditions(i,-1,PERIODIC) ;
00277     else SetBoundaryConditions(i,-1,-1) ;
00278     
00279     Ext.islifting   [i]=false ;
00280     Ext.ismultiscale[i]=true  ;
00281   }
00282   Ext.dimext=0     ;
00283    
00284   // enforce continuous Dirichlet values
00285   if (!ContinuousDirichletValues) {
00286 
00287     // compute nodal values for all Dirichlet faces
00288     for (d=0; d<DIM; d++)
00289       for (k=0; k<2; k++)
00290         if (BC[d][k]==0) 
00291           for (i=0; i<DIM-1; i++) 
00292             if (F[d][k]->Ext.ismultiscale[i]) F[d][k]->ApplyOp(F[d][k] , i , INVERSETRANSFORM) ;
00293        
00294     // averaging on boundaries of the Dirichlet faces
00295     for (d=0; d<DIM; d++)   
00296       for (k=0; k<2; k++)
00297         if (BC[d][k]==0) {
00298 
00299           // start building up the global DIM-dimensional (l,t)-index for the boundary coeffcients
00300           int l[DIM],t[DIM],ll[DIM-1],kk[DIM-1],ss[DIM-1] ;
00301           l[d]=G->WC->Level0 ; 
00302           t[d]=k*nl0  ;
00303           
00304           // iterate through all indices of the face
00305           AdaptiveGrid<DIM>::ProjectionGrid          *pg=F[d][k]->G->BasisPG[0] ;
00306           AdaptiveGrid<DIM>::ProjectionGrid::iterator pgit ;
00307 
00308           for (pgit=(*pg).begin(); pgit!=(*pg).end(); ++pgit) {
00309             AdaptiveLevelsIndex *ali=(*pgit).first  ;
00310             AdaptiveLevels      *als=(*pgit).second ;
00311 
00312             // continue to build up the global (l,t)-index
00313             // d0 is the first variable coordinate of the face (w.r.t. to the DIM-dimensional cube)
00314             int j=0, d0 = (d==0) ? 1 : 0 ;
00315             for (i=0; i<DIM; i++) {
00316               if (i==d ) continue ; // l[d],t[d] already set
00317               if (i==d0) continue ; // ld0],t[d0] to be set in the loop below
00318               l[i]= ali->l[j] ;
00319               t[i]= ali->t[j] ;
00320               j++ ;
00321             }
00322 
00323             for (i=1; i<DIM-1; i++) ll[i]=ali->l[i-1] ;
00324       
00325             for (int l0=G->Level0[0]; l0<=LMAX; l0++) {
00326               AdaptiveLevels::AdaptiveL *al=(*als)[l0] ;
00327               AdaptiveLevels::AdaptiveL::iterator alit ;
00328         
00329               ll[0]=l[d0]=l0 ;
00330               for (alit=(*al).begin(); alit != (*al).end(); alit++) {
00331                 t[d0]=(*alit).first ;
00332                 // number of ~ / pointers to  values to be averaged
00333                 int     n=1 ;
00334                 double *v[DIM]={ &(*F[d][k]->ba->d[ll])[(*alit).second] } ;
00335                 
00336                 // check global index on which Dirichlet-faces it lies and do the averaging
00337                 bool doaverage=true ;
00338                 for (i=0; i<DIM; i++) {
00339                   if (i==d) continue ; // dth coordinate doesn't need to be tested
00340 
00341                   if ( (l[i]==G->WC->Level0) && ( ((t[i]==0) && (BC[i][0]==0)) || ((t[i]==nl0) && (BC[i][1]==0))) ) {
00342                     
00343                     k2 = (t[i]==0) ? 0 : 1 ; // which side 
00344                     
00345                     if (i<d) doaverage=false ; // the coefficient has been averaged already
00346                     else {
00347                       // build (DIM-1)-dimensional index for the other face[i][k2]
00348                       i2=0 ;
00349                       for (j=0; j<DIM; j++)
00350                         if (j!=i) { kk[i2]=l[j], ss[i2]=t[j], i2++; }
00351                       // find pointer to value to be averaged
00352                       v[n]=F[i][k2]->GetP(kk,ss) ;
00353                       n++ ;
00354                     }
00355                   }
00356 
00357                   if (!doaverage) break ; // no averaging needed 
00358                 }
00359              
00360                 // averaging
00361                 if (doaverage && (n>1)) {
00362                   double x=0 ;
00363                   for (i=0; i<n; i++) x += *v[i] ;
00364                   x /= n ;
00365                   for (i=0; i<n; i++) *v[i]=x ;
00366 
00367                 }              
00368               } // next alit
00369             }  // next l0             
00370           }    // next pgit
00371 
00372         }      // next (DIM-1)-dimensional face
00373         
00374   }  // ContinuousDirichletValues
00375   
00376   // compute multiscale representation for all Dirichlet faces
00377   for (d=0; d<DIM; d++) 
00378     for (k=0; k<2; k++)
00379       if (BC[d][k]==0) 
00380         for (i=0; i<DIM-1; i++)
00381           if (!F[d][k]->Ext.ismultiscale[i]) F[d][k]->ApplyOp(F[d][k] , i , TRANSFORM) ;
00382 
00383 
00384   //  set Dirichlet values  
00385   for (i=0; i<DIM; i++)
00386     for (k=0; k<2; k++)
00387       if (BC[i][k]==0)
00388         F[i][k]->WriteSlice(this) ;
00389 } 
00390 
00392 //                       //
00393 //  auxiliary functions  //
00394 //                       //
00396 template<>
00397 void AdaptiveData<2>::ComputeDerivatives(AdaptiveData<1> *F[2][2] , int BC[2][2], int FD, int ) {
00398   int i,k ;
00399   for (i=0; i<2; i++)
00400     for (k=0; k<2; k++)
00401       if (BC[i][k]==1) {
00402         // first derivatives d_j on (DIM-1)-dimensional faces
00403         assert(F[i][k]->G->tmp[1]) ;
00404         F[i][k]->G->tmp[1]->ApplyOp(F[i][k]->G->tmp[0] , 0 , FD) ;
00405       }
00406 }
00407 
00408 template<int DIM>
00409 AdaptiveData<DIM-3> *AdaptiveData<DIM>::GetDi2Di1Fi0k0(AdaptiveData<DIM-1> *F[DIM][2] , int i2,int k2, int i1, int k1, int i0, int k0) {
00410   int i1s, i2s ;
00411   i1s=(i1<i0) ? i1 : i1-1;
00412   i2s=i2 ;
00413   if ((i2>i0) && (i2<i1)) i2s=i2-1 ;
00414   if ((i2>i1) && (i2<i0)) i2s=i2-1 ;
00415   if ((i2>i0) && (i2>i1)) i2s=i2-2 ;
00416 
00417   return F[i0][k0]->G->FaceGrid[i1s][k1]->FaceGrid[i2s][k2]->tmp[0] ;
00418 }
00419 
00420 template<int DIM>
00421 void AdaptiveData<DIM>::ComputeDerivatives(AdaptiveData<DIM-1> *F[DIM][2] , int BC[DIM][2], int FD, int codim) {
00422   int i0,k0,i1,k1,i2,k2 ;
00423   for (i0=0; i0<DIM; i0++)
00424     for (k0=0; k0<2; k0++)
00425       if (BC[i0][k0]==1) {
00426         AdaptiveGrid<DIM-1> *Gi0k0=F[i0][k0]->G ;
00427         // d_i1 f_{i0,k0} on S_{i0,k0} == F[i0][k0]->G->tmp[i1^* + 1]
00428         for (i1=0; i1<DIM-1; i1++) {
00429           assert(Gi0k0->tmp[i1+1]) ;
00430           Gi0k0->tmp[i1+1]->ApplyOp(Gi0k0->tmp[0] , i1 , FD) ;
00431         }
00432         
00433         // d_i1 d_i2 f_{i0,k0} on S_{i0,k0} \bigcap S_{i1,k1} \bigcap S_{i2,k2}  ==  F[i0][k0]->G->FaceGrid[i1^*][k1]->FaceGrid[i2^**][k2]->tmp[0]
00434         if (codim==3) {
00435           for (i1=0; i1<DIM-1; i1++) {
00436             for (k1=0; k1<DIM-1; k1++) {
00437               AdaptiveGrid<DIM-2> *Gi1k1=Gi0k0->FaceGrid[i1][k1] ;
00438               Gi1k1->tmp[0]->ReadSlice(Gi0k0->tmp[i1 + 1]) ; // d_i1
00439               for (i2=0; i2<DIM-2; i2++) {
00440                 Gi1k1->tmp[i2+1]->ApplyOp( Gi1k1->tmp[0] , i2 , FD) ; // d_i2 d_i1
00441                 for (k2=0; k2<2; k2++) {
00442                   AdaptiveGrid<DIM-3> *Gi2k2=Gi1k1->FaceGrid[i2][k2] ;
00443                   Gi2k2->tmp[0]->ReadSlice(Gi1k1->tmp[i2 + 1]) ;
00444                 }
00445               }
00446             }
00447           }
00448 
00449           // averaging of mixed derivatives: 0.5*(d_i2 d_i1 + d_i1 d_i2)
00450           for (i1=0; i1<DIM; i1++)
00451             for (i2=0; i2<DIM; i2++) {
00452               if (i1==i0) continue ;
00453               if (i2==i0) continue ;
00454               if (i1==i2) continue ;
00455 
00456               for (k1=0; k1<2; k1++)
00457                 for (k2=0; k2<2; k2++) {
00458                   AdaptiveData<DIM-3> *di2i1=GetDi2Di1Fi0k0(F,i2,k2,i1,k1,i0,k0) ;
00459                   AdaptiveData<DIM-3> *di1i2=GetDi2Di1Fi0k0(F,i1,k1,i2,k2,i0,k0) ;
00460                   di2i1->Add(0.5,di2i1 , 0.5,di1i2) ;
00461                   di1i2->Copy(di2i1) ; 
00462                 }
00463             }
00464         } // codim==3
00465       } // next i0,k0
00466 }
00467 
00468 
00469 template<int DIM>
00470 void AdaptiveData<DIM>::AddBiTriLinearContributionCodim2(int *l,int *t,int *n,int *f,int codim, AdaptiveData<DIM-1> *F[DIM][2],int (*)[2] ) {
00471   int j,r,lr[DIM],tr[DIM],a[DIM],e[DIM],q,m,nl0=1<<G->WC->Level0 ;
00472  
00473   assert(codim==2) ;
00474   for (j=0; j<codim; j++)
00475     for (r=j+1; r<codim; r++) {
00476 
00477       // get derivatives
00478       q=0; for (m=0; m<DIM; m++) if (m!=n[j]) {lr[q]=l[m], tr[q]=t[m], q++;} // reduced l/t indices
00479       m= (n[r]<n[j]) ? n[r] : n[r]-1 ;                                       // derivative where to find on face ?
00480       double dnrfnjfj=F[n[j]][f[j]]->G->tmp[1+m]->Get(lr,tr) ;
00481       
00482       q=0; for (m=0; m<DIM; m++) if (m!=n[r]) {lr[q]=l[m], tr[q]=t[m], q++;}
00483       m= (n[j]<n[r]) ? n[j] : n[j]-1 ;
00484       double dnjfnrfr=F[n[r]][f[r]]->G->tmp[1+m]->Get(lr,tr) ; 
00485       
00486       // iterate over all multiindices tr
00487       InitCodimIteration(t,n,f,codim,G->WC  ,  tr,a,e) ;
00488       while (true) {
00489         
00490         double *p=GetP(l,tr) ;
00491         (*p) += 0.5*(dnrfnjfj + dnjfnrfr)*
00492           ( ((double)tr[n[j]])/nl0 - f[j] )*( ((double)tr[n[r]])/nl0 - f[r] );
00493         // advance iterator
00494         q=0 ;
00495         while (q<codim) {
00496           tr[n[q]]++ ;
00497           if (tr[n[q]]>e[q]) { tr[n[q]]=a[q], q++;}
00498           else break ;        // stop advancing the iterator
00499         }
00500         if (q==codim) break ; // stop iteration
00501       } 
00502     } // next r,j
00503 }
00504 
00505 template<>
00506 void AdaptiveData<2>::AddBiTriLinearContribution(int *l,int *t,int *n,int *f,int codim, AdaptiveData<1> *F[2][2],int BC[2][2]) {
00507   AddBiTriLinearContributionCodim2(l,t,n,f,codim,F,BC) ;
00508 }
00509 
00510 
00511 template<int DIM>
00512 void AdaptiveData<DIM>::AddBiTriLinearContribution(int *l,int *t,int *n,int *f,int codim, AdaptiveData<DIM-1> *F[DIM][2],int BC[DIM][2]) {
00513   if (codim==2) {
00514     AddBiTriLinearContributionCodim2(l,t,n,f,codim,F,BC) ;
00515     return ;
00516   }
00517   assert(codim==3) ;
00518 
00519   int j,r,lr[DIM],tr[DIM],a[DIM],e[DIM],q,m,k=0,nl0=1<<G->WC->Level0 ;
00520  
00521   // bilinear contribution
00522   for (j=0; j<codim; j++)
00523     for (r=j+1; r<codim; r++) {
00524       // codimensions n[j],n[r] vary, get the remaining codimension
00525       if ((j==0) && (r==1)) k=2 ;
00526       if ((j==0) && (r==2)) k=1 ;
00527       if  (j==1)            k=0 ;
00528       
00529       // get derivative dnrfnjfj
00530       m= (n[r]<n[j]) ? n[r] : n[r]-1 ;
00531       AdaptiveData<DIM-1> *dnrFnjfj=F[n[j]][f[j]]->G->tmp[1+m] ;
00532       //
00533       q=0; for (m=0; m<DIM; m++) if (m!=n[j]) {lr[q]=l[m], tr[q]=t[m], q++;} // reduced l/t indices
00534       m = (n[k]<n[j]) ? n[k] : n[k]-1 ;
00535       dnrFnjfj->Get(lr,tr, m , G->AB1[0] , (f[k]==0) ? 1 : 2 ) ;
00536 
00537       // get derivative dnjfnrfr
00538       m= (n[j]<n[r]) ? n[j] : n[j]-1 ;
00539       AdaptiveData<DIM-1> *dnjFnrfr=F[n[r]][f[r]]->G->tmp[1+m] ;
00540       //
00541       q=0; for (m=0; m<DIM; m++) if (m!=n[r]) {lr[q]=l[m], tr[q]=t[m], q++;} // reduced l/t indices
00542       m = (n[k]<n[r]) ? n[k] : n[k]-1 ;
00543       dnjFnrfr->Get(lr,tr, m , G->AB2[0] , (f[k]==0) ? 1 : 2 ) ;
00544 
00545       G->AB1[0]->Add(0.5,G->AB1[0] , 0.5,G->AB2[0]) ;
00546       G->AB2[0]->MMinus(G->AB1[0]) ;
00547       
00548       // iterate through tr[n[j]],tr[n[r]]
00549       InitCodimIteration(t,n,f,codim,G->WC  ,  tr,a,e) ;
00550       while (true) {
00551         double xy=( ((double)tr[n[j]])/nl0 - f[j] )*( ((double)tr[n[r]])/nl0 - f[r] ) ;
00552         G->AB2[0]->Add(0.0,G->AB2[0] , xy , G->AB1[0]) ;
00553         Set(l,tr, n[k] , G->AB2[0] , (f[k]==0) ? 1 : 2 , true) ;
00554         // advance iterator
00555         tr[n[j]]++; 
00556         if (tr[n[j]]>e[j]) {
00557           tr[n[j]]=a[j] ;
00558           tr[n[r]]++    ;
00559           if (tr[n[r]]>e[r]) break ; // stop while loop
00560         }
00561       }
00562     } // next r,j
00563   
00564   // trilinear contribution
00565   q=0; 
00566   for (m=0; m<DIM; m++) {
00567     if (m==n[0]) continue ;
00568     if (m==n[1]) continue ;
00569     if (m==n[2]) continue ;
00570     
00571     lr[q]=l[m] ;
00572     tr[q]=t[m] ;
00573     q++ ;
00574   }
00575   double c=0.0 ;
00576   for (j=0; j<codim; j++) {
00577     if (j==0) r=1, k=2 ;
00578     if (j==1) r=0, k=2 ;
00579     if (j==2) r=0, k=1 ;
00580  
00581     c += GetDi2Di1Fi0k0(F,n[k],f[k],n[r],f[r],n[j],f[j])->Get(lr,tr) ;
00582   }
00583   c=-2*c/3 ;
00584   
00585   // iterate over all multiindices tr
00586   InitCodimIteration(t,n,f,codim,G->WC  ,  tr,a,e) ;
00587   while (true) {   
00588     double *p=GetP(l,tr) ;
00589     (*p) += c*( ((double)tr[n[0]])/nl0 - f[0] )*
00590               ( ((double)tr[n[1]])/nl0 - f[1] )*
00591               ( ((double)tr[n[2]])/nl0 - f[2] );
00592     // advance iterator
00593     q=0 ;
00594     while (q<codim) {
00595       tr[n[q]]++ ;
00596       if (tr[n[q]]>e[q]) { tr[n[q]]=a[q], q++;}
00597       else break ;        // stop advancing the iterator
00598     }
00599     if (q==codim) break ; // stop iteration
00600   }
00601   
00602 
00603 }
00604 
00605 
00606 template<int DIM>
00607 void AdaptiveData<DIM>::SetNeumannBoundaryValues(AdaptiveData<DIM-1> *F[DIM][2] , int BC[DIM][2]) { 
00608   int i,j,k,m, nl0= 1<<G->WC->Level0, codim,codimmax=0, a[DIM],e[DIM] ;
00609   //int FD=FIRSTDERIVATIVE ;
00610   int FD=FD14 ;
00611     
00612   // get maximal codimension of Neumann-Neumann--faces
00613   for (i=0; i<DIM; i++)
00614     if ((BC[i][0]==1) || (BC[i][1]==1)) codimmax++ ;
00615  
00616   if (codimmax==0) return ; // no Neumann faces
00617   if (codimmax >3) {
00618     std::cout<<"AdaptiveData<"<<DIM<<">::SetBoundaryValueFunction: maximal codimension of Neumann-Neumann-...- faces must not exceed 3 for the current implementation\n" ;
00619     exit(-1) ;
00620   }
00621 
00622   // compute (copy) modified Neumann values
00623   for (i=0; i<DIM; i++) 
00624     if ((BC[i][0]==1) || (BC[i][1]==1)) {
00625       bool neig=false ; // is there a neighbouring Dirichlet face ?
00626       for (j=0; j<DIM; j++) { 
00627         if (i==j) continue ;
00628         if ((BC[j][0]==0) || (BC[j][1]==0)) {neig=true; break;}
00629       }
00630       
00631       if (!neig) { // just copy
00632         for (k=0; k<2; k++)
00633           if (BC[i][k]==1) {
00634             assert (F[i][k]->G->tmp[0]) ;
00635             F[i][k]->G->tmp[0]->Copy(F[i][k]) ;
00636           }
00637       } else {    // compute modified Neumann-values
00638         G->tmp[0]->ApplyOp(this,i,FD) ;
00639         for (k=0; k<2; k++)
00640           if (BC[i][k]==1) {
00641             F[i][k]->G->tmp[1]->ReadSlice(G->tmp[0]) ;
00642             F[i][k]->G->tmp[0]->Sub(F[i][k] , F[i][k]->G->tmp[1]) ;
00643             // Dirichlet-values have a higher priority than Neumann-values, 
00644             // but if the (tangential) derivatives of the Dirichlet-values is not consistent with
00645             // the Neumann-values, setting the Neumann values alters values on the Dirichlet faces
00646             // To circumvent this problem, the (modified) Neumann-values are projected to be zero
00647             // on Dirichlet-Neumann boundaries
00648             for (j=0; j<DIM; j++) {
00649               if (i==j) continue ;
00650 
00651               bool project=false ;
00652               int  GBC[2]={-1,-1} , DBC[2]={-1,-1} ;
00653               for (m=0; m<2; m++) 
00654                 if (BC[j][m]==0) DBC[m]=0, project=true ;
00655 
00656               if (project) {
00657                 m=(j<i) ? j : j-1 ;
00658                 F[i][k]->G->tmp[0]->ApplyOp(DBC, F[i][k]->G->tmp[0] , m , PROJECTION) ;
00659                 F[i][k]->G->tmp[0]->ApplyOp(GBC, F[i][k]->G->tmp[0] , m , PROJECTION) ;
00660               }
00661 
00662             }
00663           }
00664       }
00665     }
00666 
00667   // recursively treat dim-dimensional Neumann-Neumann-....-Neumann corners, edges, ...
00668   for (codim=codimmax; codim>=2; codim--) {
00669 
00670     G->tmp[1]->Copy(this) ;
00671     ComputeDerivatives(F,BC,FD,codim) ;
00672 
00673     // iterate through all coefficients
00674     AdaptiveGrid<DIM>::ProjectionGrid          *pg=G->BasisPG[0] ;
00675     AdaptiveGrid<DIM>::ProjectionGrid::iterator pgit ;
00676     
00677     for (pgit=(*pg).begin(); pgit!=(*pg).end(); ++pgit) {
00678       AdaptiveLevelsIndex *ali=(*pgit).first  ;
00679       AdaptiveLevels      *als=(*pgit).second ;   
00680 
00681       int l[DIM],t[DIM],l0 ;
00682       for (i=1; i<DIM; i++) {
00683         l[i]=ali->l[i-1] ;
00684         t[i]=ali->t[i-1] ;
00685       }
00686 
00687       for (l0=G->Level0[0]; l0<=LMAX; l0++) {
00688         AdaptiveLevels::AdaptiveL *al=(*als)[l0] ;
00689         AdaptiveLevels::AdaptiveL::iterator alit ;
00690         
00691         l[0]=l0 ;
00692         for (alit=(*al).begin(); alit != (*al).end(); alit++) {
00693           t[0]=(*alit).first ;
00694            
00695           // on which Neumann-faces does (l,t) live ?
00696           int n[DIM], f[DIM], codimlt=0, q, r, lr[DIM], tr[DIM] ;
00697           for (i=0; i<DIM; i++) 
00698             if (l[i]==G->Level0[i]) {
00699               if ((t[i]==  0) && (BC[i][0]==1)) {n[codimlt]=i, f[codimlt]=0, codimlt++; }
00700               if ((t[i]==nl0) && (BC[i][1]==1)) {n[codimlt]=i, f[codimlt]=1, codimlt++; }
00701             }
00702 
00703           if (codimlt==codim) {
00704   
00706             // add local linear contribution near the face //
00708             for (j=0; j<codim; j++) {
00709               // get modified boundary value of face (n[j],f[j])
00710               q=0; for (r=0; r<DIM; r++) if (r!=n[j]) {lr[q]=l[r], tr[q]=t[r], q++;}
00711               double fnk=F[n[j]][f[j]]->G->tmp[0]->Get(lr,tr) ;
00712 
00713               // iterate over all multiindices tr
00714               InitCodimIteration(t,n,f,codim,G->WC  ,  tr,a,e) ;
00715               while (true) {
00716                 double *p=GetP(l,tr) ;
00717                 (*p) += (fnk*tr[n[j]])/nl0  ;
00718                 // advance iterator
00719                 q=0 ;
00720                 while (q<codim) {
00721                   tr[n[q]]++ ;
00722                   if (tr[n[q]]>e[q]) { tr[n[q]]=a[q], q++;}
00723                   else break ;        // stop advancing the iterator
00724                 }
00725                 if (q==codim) break ; // stop iteration
00726               }
00727             }
00729             // add bi-linear contribution //
00731             AddBiTriLinearContribution(l,t,n,f,codim,F,BC) ;  
00732           }     // codim
00733 
00734         } // next alit,l0,pgit
00735       }
00736     }
00737     
00738     // compute modified Neumann-values
00739     for (i=0; i<DIM; i++)
00740       if ((BC[i][0]==1) || (BC[i][1]==1)) {
00741         G->tmp[0]->ApplyOp(this,i,FD) ;
00742         for (k=0; k<2; k++)
00743           if (BC[i][k]==1) {
00744             F[i][k]->G->tmp[1]->ReadSlice(G->tmp[0]) ;
00745             F[i][k]->G->tmp[0]->Sub(F[i][k] , F[i][k]->G->tmp[1]) ;
00746           }
00747       }
00748   } // next dim
00749 
00750   SetNeumannFaces(F,BC) ;
00751 }
00752 
00753 //
00754 template<int DIM>
00755 void AdaptiveData<DIM>::InitCodimIteration(int *t, int *n, int *f, int codim, Wavelets *W, int *tr,int *a, int *e) {
00756   int i,nl0=1<<W->Level0 ;
00757  
00758   for (i=0; i<DIM; i++) tr[i]=t[i] ; // fix components
00759   for (i=0; i<codim; i++) 
00760     if (f[i]==0) tr[n[i]]=a[i]=0          , e[i]=W->N-1;
00761     else         tr[n[i]]=a[i]=nl0-W->N+1 , e[i]=nl0   ;
00762 }
00763 
00764 //
00765 //
00766 template<int DIM>
00767 void AdaptiveData<DIM>::SetNeumannFaces(AdaptiveData<DIM-1> *F[DIM][2] , int BC[DIM][2]) { 
00768   int i,k,in[2]={0,0}, nl0= 1<<G->WC->Level0 ;
00769   double dphidx0, dphidx1 ;   
00770 
00771   // get first derivatives of left/right coarse scaling functions
00772   dphidx0=G->WC->Operators[0][0][0].AL[0][0][in] ;
00773     
00774   in[0]=G->WC->Operators[0][0][0].AR[0][0].size[0]-1 ;
00775   in[1]=G->WC->Operators[0][0][0].AR[0][0].size[1]-1 ;
00776   dphidx1=G->WC->Operators[0][0][0].AR[0][0][in] ;
00777 
00778   // set Neumann boundary conditions
00779   for (i=0; i<DIM; i++)
00780     for (k=0; k<2; k++)
00781       if (BC[i][k]==1) {
00782         F[i][k]->G->tmp[0]->Mul( ((k==1) ? 1/dphidx1 : 1/dphidx0)/nl0 ) ;
00783         F[i][k]->G->tmp[0]->WriteSlice(this,true) ;
00784       }
00785 }
00786 
00787 //
00788 //
00789 template<int DIM>
00790 void AdaptiveData<DIM>::SetBoundaryValueFunction(AdaptiveData<DIM-1> *F[DIM][2] , int BC[DIM][2], bool ContinuousDirichletValues) {
00791 
00792   if (!G->WC->InterpoletFlag) {
00793     std::cout<< "Grid::SetBoundaryValueFunction: supported for Interpolets only!\n" ;
00794     exit(-1) ;
00795   }
00796 
00797   SetDirichletBoundaryValues(F,BC,ContinuousDirichletValues) ;
00798   SetNeumannBoundaryValues  (F,BC) ;
00799 }
00800 
00801 
00802 // Integrate
00803 template<int DIM>
00804 double AdaptiveData<DIM>::Integrate() {
00805   
00806  AdaptiveGrid<DIM>::ProjectionGrid              *pg=G->BasisPG[0] ;
00807  AdaptiveGrid<DIM>::ProjectionGrid::iterator     pgit ;
00808   
00809  AdaptiveLevels::AdaptiveL          *al   ;
00810  AdaptiveLevels::AdaptiveL::iterator alit ;
00811 
00812  int    i,l[DIM] ;
00813  double integral=0.0,domain=1;
00814  
00815  // not implemented for lifting wavelets, only for interpolet WAVELETS !
00816  for (i=0; i<DIM; i++) {
00817    domain *= (Ext.XE[i] - Ext.XA[i]) ;
00818    assert(!Ext.islifting   [i]) ;
00819    assert( Ext.ismultiscale[i]) ;
00820  }
00821 
00822  //
00823  for (pgit=(*pg).begin(); pgit!=(*pg).end(); ++pgit) {
00824    AdaptiveLevelsIndex *ali=(*pgit).first  ;
00825    AdaptiveLevels      *als=(*pgit).second ;
00826     
00827    // integrals of last D-1 factors of multivariate basis functions      
00828    double c1=1.0 ;
00829    for (i=1; i<DIM; i++) {
00830      l[i]=ali->l[i-1] ;
00831      c1*=G->WC->WaveletIntegral(l[i] , ali->t[i-1] , Ext.BC[i]) ;
00832    }
00833 
00834    for (l[0]=G->Level0[0]; l[0]<=LMAX; l[0]++) {
00835      al=(*als)[l[0]] ;
00836      for (alit=(*al).begin(); alit != (*al).end(); alit++) {
00837        double  w=(*ba->d[l])[ (*alit).second ] ;
00838        integral += w*c1*G->WC->WaveletIntegral(l[0] , (int)(*alit).first , Ext.BC[0]) ;
00839      }
00840    }
00841  }   // next pgit
00842 
00843  return integral*domain ;
00844 }
00845 
00846 
00847 template<int DIM>
00848 void AdaptiveData<DIM>::ApplyDiagonal(AdaptiveData<DIM> *X, double p) {
00849 
00850  AdaptiveDataArray<DIM>::VectorArray::iterator it(&ba->d),dend=ba->d.end() ;
00851  AdaptiveDataArray<DIM>::DataVector *vx,*vr ;
00852  size_t i,n  ;
00853  double diag ;
00854  int    j    ;
00855 
00856  CheckAdaptiveGrid(X) ;
00857  CopyExtensions(X) ;
00858 
00859  // diagonal scaling
00860  // iterate over all level
00861  for (it=ba->d.begin(); it<=dend; ++it) {
00862    vx=X->ba->d[it] ; vr=ba->d[it] ; n=(*vx).size() ;
00863    assert(n==(*vr).size()) ;
00864    
00865    diag=0.0 ;
00866    for (j=0; j<DIM; j++) diag += pow(1<<it.i[j] , p) ;
00867 
00868    for (i=0; i<n; i++) (*vr)[i] = (*vx)[i] * diag ;
00869  }
00870 }
00871 
00872 template<int DIM>
00873 void   AdaptiveData<DIM>::SetDiagonalMatrix(const double p, const double alpha) {
00874  AdaptiveDataArray<DIM>::VectorArray::iterator it(&ba->d), dend=ba->d.end() ;
00875  AdaptiveDataArray<DIM>::DataVector *vr ;
00876  size_t i,n  ;
00877  double diag ;
00878  int    j    ;
00879 
00880  // diagonal scaling
00881  // iterate over all level
00882  for (it=ba->d.begin(); it<=dend; ++it) {
00883    vr=ba->d[it] ; n=(*vr).size() ;
00884    assert(n==(*vr).size()) ;
00885    
00886    diag=0.0 ;
00887    for (j=0; j<DIM; j++) diag += pow(1<<it.i[j] , p) ;
00888   
00889    for (i=0; i<n; i++) (*vr)[i] = 1 + alpha*diag ;
00890  }
00891 }
00892 
00893 
00894 
00895 
00896 
00897 
00898 //
00899 //
00900 // Adaptive->SparseMatrix Output
00901 //
00902 //
00903 
00904 template<int DIM>
00905 void AdaptiveData<DIM>::WriteSparse(const char *name, AdaptiveData<DIM> **MM, int N, bool CastToFloat, AdaptivityCriterion *R) {
00906  AdaptiveData<DIM> *M[100] ;
00907  
00908  AdaptiveGrid<DIM>::ProjectionGrid              *pg=G->BasisPG[0] ;
00909  AdaptiveGrid<DIM>::ProjectionGrid::iterator     pgit             ;
00910   
00911  AdaptiveLevels::AdaptiveL          *al   ;
00912  AdaptiveLevels::AdaptiveL::iterator alit ;
00913 
00914  int      i,l[DIM],ld ;
00915  unsigned long x[DIM],t[DIM] ;
00916  long     j ;
00917 
00918  int header[100]={0},anz=0,hs ;
00919 
00920  double wd[100] ;
00921  float  wf[100] ;
00922 
00923  bool   do_swap = UDF_CPUisBigEndian() > 0 ;
00924 
00925  FILE *fid=fopen(name,"w") ;
00926  if (fid==NULL) {
00927    std::cout<<"WriteSparse::Could not open "<<name<<'\n';
00928    exit(-1) ;
00929  }
00930 
00931  // check grids
00932  if (N>0) {
00933    assert(MM) ;
00934    for (i=0; i<N; i++) {
00935      M[i]=MM[i] ;
00936      assert(M[i])       ;
00937      assert(G==M[i]->G) ;
00938    }
00939  } else {
00940    M[0]=this; N=1 ;
00941  }
00942 
00943  // dummy write !
00944  fwrite(header,sizeof(int) , sizeof(header)/sizeof(header[0]) , fid) ;
00945   
00946  //
00947  for (pgit=(*pg).begin(); pgit!=(*pg).end(); ++pgit) {
00948    AdaptiveLevelsIndex *ali=(*pgit).first  ;
00949    AdaptiveLevels      *als=(*pgit).second ;
00950 
00951    for (i=1; i<DIM; i++) { l[i]=ali->l[i-1], t[i]=ali->t[i-1] ;}
00952          
00953    for (i=1; i<DIM; i++) code(l[i],t[i] , &x[i],LMAX+1,G->Level0[i]) ;
00954 
00955    if (do_swap) UDF_swap(&x[1],sizeof(x[0]),DIM-1) ;
00956    
00957    for (ld=G->Level0[0]; ld<=LMAX;ld++) {
00958      l[0]=ld    ;
00959      al=(*als)[ld] ;
00960      for (i=0; i<N; i++) assert (M[i]->ba->d[l]) ; // missing resize/refine ?
00961      
00962      for (alit=(*al).begin(); alit != (*al).end(); alit++) { 
00963        code(l[0],(unsigned long)(*alit).first , &x[0] , LMAX+1,G->Level0[0]) ;
00964 
00965        // collect all numerical values
00966        j=(*alit).second ;
00967        if (j<0) j=1-j   ;
00968        for (i=0; i<N; i++) wd[i]=(*M[i]->ba->d[l])[j] ;
00969 
00970        if ( (R==NULL) || (R->IsActive(wd[0], l, Ext.WC->Level0 , DIM)) ) {
00971            
00972          if (do_swap) UDF_swap(&x[0],sizeof(x[0]),1) ;
00973          
00974          fwrite(x,sizeof(x[0]),DIM,fid) ;
00975            
00976          // (cast) + write 
00977          if (CastToFloat) {
00978            for (i=0; i<N; i++) wf[i]=(float)wd[i] ;
00979            if (do_swap) UDF_swap(wf,sizeof(float),N) ;
00980            fwrite(wf , sizeof(float),N,fid) ;
00981          } else {
00982            if (do_swap) UDF_swap(wd,sizeof(double),N) ;
00983            fwrite(wd , sizeof(double),N,fid) ;
00984          }
00985          anz++ ;
00986        } // if active
00987 
00988      }
00989    }
00990  }   // next pgit
00991  fclose(fid) ;
00992  
00993  // write true header
00994  hs=header[0]=sizeof(header)/sizeof(header[0]) ;
00995  header[1]=DIM    ;
00996  header[2]=anz    ;
00997  header[3]=LMAX+1 ;
00998  
00999  float *p=(float *)(&(header[4])) ;
01000  assert(sizeof(int)==sizeof(float)) ;
01001  
01002  for (i=0; i<DIM; i++) {
01003    *p=(float)Ext.XA[i] ; p++ ;
01004    *p=(float)Ext.XE[i] ; p++ ;
01005  }
01006  
01007  header[4+2*DIM] = N ; // number of entries
01008  header[5+2*DIM] = (CastToFloat ? 1 : 0) ;
01009  
01010  // wavelets, flags, boundary conditions
01011  header[6+2*DIM]=Ext.WC->N ; // must be non zero as this is used as flag, that extended information is stored
01012  header[7+2*DIM]=Ext.WC->TypeID(Ext.WC->type) ;
01013  int st1=0,st2=0,si ;
01014  for (i=0; i<DIM; i++) {
01015    si=1<<i ;
01016    if (Ext.islifting   [i]) st1 |= si ;
01017    if (Ext.ismultiscale[i]) st2 |= si ;
01018  }
01019  header[8+2*DIM]=st1 ;
01020  header[9+2*DIM]=st2 ;
01021  for (i=0; i<DIM; i++) {
01022    header[10+2*DIM+i*2]=Ext.BC[i][0] ;
01023    header[11+2*DIM+i*2]=Ext.BC[i][1] ;
01024  }
01025 
01026 
01027  fid=fopen(name,"r+") ;
01028  if (do_swap) UDF_swap(header,sizeof(int),hs) ;
01029  fwrite(header,sizeof(int) , hs , fid) ;
01030  fclose(fid) ;
01031 }
01032 
01033 template<int DIM>
01034 void AdaptiveData<DIM>::WriteSparse(const char *name, AdaptivityCriterion *R) {
01035  AdaptiveData<DIM> *M[1]={this} ;
01036  WriteSparse(name,M,1,false,R) ;
01037 }
01038 
01039 
01040 // SparseMatrix -> Adaptive Grid
01041 template<int DIM>
01042 void AdaptiveGrid<DIM>::ReadSparse(const char *name) {
01043   
01044  int           i,l[DIM],si[DIM] ;
01045  unsigned long x[DIM],s[DIM]    ;
01046  long          j,lm[DIM]={0}    ;
01047  double        wd[100] ;
01048  float         wf[100] ;
01049 
01050  int header[200]={0},anz,jm1,N ;
01051  bool CastToFloat, do_swap ;
01052 
01053  FILE *fid=fopen(name,"r") ;
01054  if (fid==NULL) {
01055    std::cout<<"Grid::ReadSparse::Could not open "<<name<<'\n' ;
01056    exit(-1) ;
01057  }
01058 
01059  // read header
01060  ReadMATLABSparseHeader(fid,DIM, header, &anz, &jm1, &N,XA,XE,&CastToFloat,&do_swap) ;
01061  
01062  // clear Grid 
01063  Clear() ;
01064  CounterArray.Set((size_t)0) ;
01065 
01066  // Grid lesen
01067  for (i=0; i<anz; i++) {
01068    fread( x,sizeof(x[0]),DIM,fid) ;
01069    if (do_swap) UDF_swap( x,sizeof(x[0]),DIM) ;
01070 
01071    if (CastToFloat) fread(wf,sizeof(float ),N,fid) ;
01072    else             fread(wd,sizeof(double),N,fid) ;
01073    
01074    // level l ,index s decodieren
01075    for (j=0; j<DIM; j++) { 
01076      decode(&l[j],&s[j] , x[j] , jm1,Level0[j]) ;
01077      si[j]=(int)s[j] ;
01078      // max level
01079      lm[j] = (l[j]<lm[j]) ? lm[j] : l[j] ;
01080     }
01081 
01082    // insert (l,s) in Grid 
01083    long *pk=FindOrInsert(BasisPG[0] , 0 , l , si) ;
01084    // node must not be active now
01085    assert(*pk == MAXINT)     ;
01086    *pk=(long)CounterArray[l] ;
01087    CounterArray[l]++         ;
01088  }
01089  fclose(fid) ;
01090 
01091  std::cout<<"max used Level:\n" ;
01092  for (i=0;i<DIM; i++) std::cout<<lm[i]<<' ' ;
01093  std::cout<<"\n" ;
01094 
01095  // generate other PGs
01096  for (i=1; i<DIM; i++) GeneratePG(BasisPG[0], 0, BasisPG[i], i) ;
01097  
01098  // resize all Data 
01099  TmpBasis->Resize(&CounterArray) ;
01100  for (i=0; i<ResizeGroup.count; i++) ResizeGroup.Ds[i]->ba->Resize(&CounterArray) ;
01101  for (i=0; i<RefineGroup.count; i++) RefineGroup.Ds[i]->ba->Resize(&CounterArray) ;
01102 
01103  SMPLoadBalance() ;
01104 }
01105 
01106 template<int DIM>
01107 void AdaptiveData<DIM>::ReadSparse(const char *name) {
01108  AdaptiveData<DIM> *M[1]={this} ;
01109  int  which[1]={0} ;
01110  ReadSparse(name,M,1,which) ;
01111 }
01112 
01113 // SparseMatrix -> Adaptive Grid
01114 template<int DIM>
01115 void AdaptiveData<DIM>::ReadSparse(const char *name, AdaptiveData<DIM> **M, int N, int *which) {
01116   
01117  int           i,l[DIM],si[DIM],NN,waveletID,waveletN ;
01118  unsigned long x[DIM],s[DIM],st1,st2,ss ;
01119  long          j ;
01120  double        wd[100]={0},XA[DIM],XE[DIM] ;
01121  float         wf[100]={0} ;
01122 
01123  int  header[200]={0},anz,jm1 ;
01124  bool CastToFloat, do_swap ;
01125 
01126  assert(N>0) ;
01127  for (i=0; i<N; i++) {
01128    assert(M[i]) ;
01129    assert(G==M[i]->G) ;
01130  }
01131 
01132  FILE *fid=fopen(name,"r") ;
01133  if (fid==NULL) {
01134    std::cout<<"Grid::ReadSparse::Could not open "<<name<<'\n' ;
01135    exit(-1) ;
01136  }
01137 
01138  // header lesen
01139  ReadMATLABSparseHeader(fid,DIM, header, &anz, &jm1, &NN,XA,XE,&CastToFloat,&do_swap) ;
01140 
01141  // get extended information if possible
01142  waveletN=header[6+2*DIM] ;
01143  if (waveletN>0) {
01144    waveletID=header[7+2*DIM] ;
01145    if ( (waveletID!=Ext.WC->TypeID(Ext.WC->type)) || (waveletN != Ext.WC->N)) {
01146      std::cout<<"File to read contains wavelet coefficients w.r.t wavelets with ID("<<waveletN<<','<<waveletID<<") but current wavelets have ID (";
01147      std::cout<<Ext.WC->N<<','<<Ext.WC->TypeID(Ext.WC->type)<<"). This may cause wrong results\n" ;
01148      assert(0) ;
01149    }
01150    st1=header[8+2*DIM] ;
01151    st2=header[9+2*DIM] ;
01152    for (i=0; i<DIM; i++) {
01153      ss=1<<i ;
01154      Ext.islifting   [i] =  ((st1&ss) != 0);
01155      Ext.ismultiscale[i] =  ((st2&ss) != 0);
01156    }
01157 
01158    for (i=0; i<DIM; i++) {
01159      Ext.BC[i][0]=header[10+2*DIM+i*2] ;
01160      Ext.BC[i][1]=header[11+2*DIM+i*2] ;
01161    }
01162    G->SetPeriodicConditions(Ext.BC) ;
01163  }
01164 
01165  if (NN!=N) {
01166    std::cout<<"!!!!!!! WARNING  !!!!!!\n" ;
01167    std::cout<<"ReadSparse: number of components in file !=  number of Data's to be filled \n"<<NN<<','<<N<<'\n' ;
01168    if (N>NN) std::cout<<"some Data's are set to 0\n" ;
01169  }
01170 
01171  // Grid/Data lesen
01172  for (i=0; i<anz; i++) {
01173    fread( x,sizeof(x[0]),DIM,fid) ;
01174    if (do_swap) UDF_swap( x,sizeof(x[0]),DIM) ;
01175 
01176    if (CastToFloat) {
01177      fread(wf,sizeof(float),NN,fid) ;
01178      if (do_swap) UDF_swap(wf,sizeof(float),NN) ;
01179      for (j=0; j<NN; j++) wd[j]=(double)wf[j] ; 
01180    } else {
01181      fread(wd,sizeof(double),NN,fid) ;
01182      if (do_swap) UDF_swap(wd,sizeof(double),NN) ;
01183    }
01184 
01185    // level l ,index s decodieren
01186    for (j=0; j<DIM; j++) { 
01187      decode(&l[j],&s[j] , x[j] , jm1,G->Level0[j]) ;
01188      si[j]=(int)s[j] ;
01189     }
01190 
01191    // insert (l,s) in Grid 
01192    long *pk=G->FindOrInsert(G->BasisPG[0] , 0 , l , si) ;
01193    assert(*pk >=0)    ;
01194    assert(*pk < (long)G->CounterArray[l]) ;
01195    for (j=0; j<N; j++)
01196      (*M[j]->ba->d[l])[*pk]=wd[which[j]] ;
01197  }
01198 
01199  fclose(fid) ;
01200 }
01201 
01202 
01203 
01204 #endif
Generated at Mon Aug 19 10:02:32 2002 for AWFD by