Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

LevelAdaptiveData.hpp

//
// Implementations for SubSpaces
//

#ifndef _SUBSPACES_
# define _SUBSPACES_

#include <assert.h>
#include "LevelAdaptive.hpp"
#include "Math.hpp"
#include "NonAdaptive.hpp"

//
// General
//
template<int DIM>
LevelAdaptiveData<DIM>::LevelAdaptiveData() {
 a=NULL ;
 G=NULL ;
 Ext.InitXAE() ;
}


template<int DIM>
LevelAdaptiveData<DIM>::LevelAdaptiveData(LevelAdaptiveGrid<DIM> *AG) {
 a=NULL  ;
 Init (AG) ;
 Ext.InitXAE() ;
}


template<int DIM>
void LevelAdaptiveData<DIM>::Init(LevelAdaptiveGrid<DIM> *AG) {
 G=AG  ;
 G->Attach(this) ;

 Ext.WC=G->W ;
 a=new Matrix< Matrix<double,DIM> *, DIM> ;
  
 int b[DIM]={0}  ;
 a->Init(b,G->L) ;
 
 Matrix< Matrix<double,DIM> *, DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) (*a)[it]=NULL ;

 Allocate() ;
}


template<int DIM>
void LevelAdaptiveData<DIM>::Clear() {
 if (a!=NULL) { 
   Matrix< Matrix<double,DIM> *, DIM>::iterator it(a) ;
   for (it=(*a).begin(); it<=(*a).end(); ++it) { delete (*a)[it] ; (*a)[it]=NULL ; }
   delete a; a=NULL  ;
 }

 G =NULL ;
}


template<int DIM>
LevelAdaptiveData<DIM>::~LevelAdaptiveData() {
 Clear() ;
}

template<int DIM>
size_t LevelAdaptiveData<DIM>::Size() {
 return G->Size() ;
}


template<int DIM>
void LevelAdaptiveData<DIM>::Allocate() {
 int i ; 
 Matrix< Matrix<double,DIM> *, DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) 
   if ( (*(G->a))[it.i] &&  ((*a)[it]==NULL) ) {
     int b[DIM],e[DIM] ;
     for (i=0; i<DIM; i++) { b[i]=0 ; e[i] = (it.i[i] == G->Level0[i]) ? 1<< G->Level0[i] : (1<<(it.i[i]-1)) -1 ; }
     (*a)[it]=new Matrix<double,DIM>(b,e) ;
    }
}


//
//
template<int DIM>
void LevelAdaptiveData<DIM>::Compatible(LevelAdaptiveData<DIM> *X) {
  assert(Ext.IsSame(&X->Ext)) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::Compatible(LevelAdaptiveData<DIM> *S0, LevelAdaptiveData<DIM> *S1) {
 Compatible(S0) ;
 Compatible(S1) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::Compatible(LevelAdaptiveData<DIM> *S0, LevelAdaptiveData<DIM> *S1, LevelAdaptiveData<DIM> *S2) {
 Compatible(S0) ;
 Compatible(S1) ;
 Compatible(S2) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::Compatible(LevelAdaptiveData<DIM> *S0, LevelAdaptiveData<DIM> *S1, LevelAdaptiveData<DIM> *S2, LevelAdaptiveData<DIM> *S3) {
 Compatible(S0) ;
 Compatible(S1) ;
 Compatible(S2) ;
 Compatible(S3) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::CopyExtensions(LevelAdaptiveData<DIM> *X) {
  Ext.CopyFlags(&X->Ext) ;
}

//
//
template<int DIM>
void LevelAdaptiveData<DIM>::SetBoundaryConditions(int B[DIM][2]) {
 Ext.SetBoundaryConditions(B) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::SetBoundaryConditions(LevelAdaptiveData<DIM> *S) {
 Ext.SetBoundaryConditions(&S->Ext) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::CheckBC(int BC0[DIM][2] , int BC1[DIM][2]) {
 int i ;
 for (i=0; i<DIM; i++) { assert(BC0[i][0]==BC1[i][0]) ; assert(BC0[i][1]==BC1[i][1]) ;}
}

//
// Copy
//
template<int DIM>
void LevelAdaptiveData<DIM>::Copy(LevelAdaptiveData<DIM> *S) {
 assert (A==S->A) ;
 
 int i ;
 for (i=0; i<DIM; i++) { BC[i][0]=S->BC[i][0] ; BC[i][1]=S->BC[i][1] ;}  
 
 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) 
    if ((*(G->a))[it.i]) (*a)[it]->Copy( (*S->a)[it] ) ;
 
 Ext.CopyFlags(&S->Ext) ;
 Ext.CopyData (&S->Ext) ;
}

//
// linear algebra
//
template<int DIM>
void LevelAdaptiveData<DIM>::PPlus(const double c0,LevelAdaptiveData<DIM> *S0) {
 Compatible(S0) ;

 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) 
    if ((*(G->a))[it.i]) (*a)[it]->PPlus(c0,(*S0->a)[it]) ;

 Ext.PPlus(c0,&S0->Ext) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::Add(const double c0,LevelAdaptiveData<DIM> *S0 , const double c1,LevelAdaptiveData<DIM> *S1) {
 S0->Compatible(S1) ;
 CopyExtensions(S0) ;

 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) if ((*(G->a))[it.i]) (*a)[it]->Add(c0,(*(S0->a))[it] , c1,(*(S1->a))[it]) ;

 Ext.Add(c0,&S0->Ext , c1,&S1->Ext) ;
}


template<int DIM>
void LevelAdaptiveData<DIM>::Add(const double c0,LevelAdaptiveData<DIM> *S0 , const double c1,LevelAdaptiveData<DIM> *S1 , const double c2,LevelAdaptiveData<DIM> *S2) {
 S0->Compatible(S1,S2) ;
 CopyExtensions(S0) ;

 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it)
   if ((*(G->a))[it.i]) (*a)[it]->Add(c0,(*(S0->a))[it] , c1,(*(S1->a))[it]  , c2,(*(S2->a))[it]) ;
 
 Ext.Add(c0,&S0->Ext , c1,&S1->Ext , c2,&S2->Ext) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::Add(const double c0,LevelAdaptiveData<DIM> *S0 , const double c1,LevelAdaptiveData<DIM> *S1 , const double c2,LevelAdaptiveData<DIM> *S2 ,
                         const double c3,LevelAdaptiveData<DIM> *S3) {
 S0->Compatible(S1,S2,S3) ;

 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it)
   if ((*(G->a))[it.i]) (*a)[it]->Add(c0,(*(S0->a))[it] , c1,(*(S1->a))[it]  , c2,(*(S2->a))[it] , c3,(*(S3->a))[it]) ;

 Ext.Add(c0,&S0->Ext , c1,&S1->Ext , c2,&S2->Ext , c3,&S3->Ext) ;
}


template<int DIM>
void LevelAdaptiveData<DIM>::Sub(LevelAdaptiveData<DIM> *S0 , LevelAdaptiveData<DIM> *S1) {
 S0->Compatible(S1) ;
 CopyExtensions(S0) ;

 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) if ((*(G->a))[it.i]) (*a)[it]->Sub((*(S0->a))[it] , (*(S1->a))[it]) ;

 Ext.Sub(c0,&S0->Ext , c1,&S1->Ext) ;
}


template<int DIM>
void LevelAdaptiveData<DIM>::Mul(double f) {

 Matrix<Matrix<double,DIM> * , DIM>::iterator l(a) ;
 for (l=(*a).begin(); l<=(*a).end(); ++l) if ((*(G->a))[l.i]) (*a)[l]->Mul(f) ;

 Ext.Mul(f);
}


template<int DIM>
double LevelAdaptiveData<DIM>::InnerProd(LevelAdaptiveData<DIM> *S) {
 double s=0 ;
 Compatible(S) ;

 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) if ((*(G->a))[it.i])  s+= (*a)[it]->InnerProd((*(S->a))[it]) ;

 return s + Ext.InnerProd(&S->Ext) ;
}

template<int DIM>
double LevelAdaptiveData<DIM>::Max() {
 double s=-1e+100,r ;
 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) 
   if ((*(G->a))[it.i]) { r=(*a)[it]->Max() ; s = max(s,r) ; }
 return max(s,Ext.Max()) ;
}

template<int DIM>
double LevelAdaptiveData<DIM>::Min() {
 double s=+1e+100,r ;
 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) 
   if ((*(G->a))[it.i]) { r=(*a)[it]->Max() ; s = min(s,r) ; }
 return  min(s,Ext.Min()) ;
}

template<int DIM>
double LevelAdaptiveData<DIM>::MaxAbs() {
 double s=-1e+100,r ;
 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) 
   if ((*(G->a))[it.i]) { r=(*a)[it]->MaxAbs() ; s = max(s,r) ; }
 return max(s,Ext.MaxAbs()) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::Set(const double c) {
 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) 
   if ((*(G->a))[it.i]) (*a)[it]->Set(c) ;
 Ext.Set(0) ;
}


//
// Matrix and IO
//

template<int DIM>
void LevelAdaptiveData<DIM>::FromUniform(UniformData<DIM> *M) {
 FromUniform(M->ba) ;
 Ext.CopyFlags(&M->Ext) ;
 Ext.CopyData (&M->Ext) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::FromUniform(Matrix<double,DIM> *M) { 
 int i,JM[DIM] ;

 M->CalcLevel(JM) ;

 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) {
   bool inner=true ;
   for (i=0; i<DIM; i++) if (it.i[i]>JM[i]) inner=false ;
   if ( (*G->a)[it.i] && inner ) {
     // copy submatrices
     int im[DIM] ;
     Matrix<double,DIM> *ait=(*a)[it] ;
     Matrix<double,DIM>::iterator in( ait ) ;
     for (in=(*ait).begin(); in <= (*ait).end(); ++in) { 
       for (i=0; i<DIM; i++) im[i]=in.i[i] + ( (it.i[i]==Level0[i]) ?  0 : (1<<(it.i[i]-1))+1 ) ;
       (*ait)[in] = (*M)[im] ;
      }
    }
  }
}


template<int DIM>
void LevelAdaptiveData<DIM>::ToUniform(UniformData<DIM> *M) {
 ToUniform(M->ba) ;
 M->Ext.CopyFlags(&Ext) ;
 M->Ext.CopyData (&Ext) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::ToUniform(Matrix<double,DIM> *M) {
 int i,JM[DIM] ;

 M->CalcLevel(JM) ;
 M->Set(0.0) ;

 Matrix<Matrix<double,DIM> * , DIM>::iterator il(a) ;
 for (il=(*a).begin(); il<=(*a).end(); ++il) {
   bool inner=true ;
   for (i=0; i<DIM; i++) if (il.i[i]>JM[i]) inner=false ;
   if ( (*G->a)[il.i] && inner ) {
     // copy submatrices
     int is[DIM] ;
     Matrix<double,DIM> *ail=(*a)[il] ;
     Matrix<double,DIM>::iterator it( ail ) ;
     for (it=(*ail).begin(); it <= (*ail).end(); ++it) { 
       for (i=0; i<DIM; i++) is[i]=it.i[i] + ( (il.i[i]==G->Level0[i]) ?  0 : (1<<(il.i[i]-1)) + 1 ) ;
       (*M)[is] = (*ail)[it] ;
      }
    }
  }
}


// WriteUDF
template<int DIM>
void LevelAdaptiveData<DIM>::WriteUDF(const char *name, int  *L, LevelAdaptiveData<DIM> *Tmp, bool CastToFloat) {
  int a[DIM]={0},e[DIM],i ;
  for (i=0; i<DIM; i++) e[i]=1<<L[i] ;
  UniformData <DIM> M(a,e,Ext.WC) ;
  WriteUDF(name,&M,Tmp,CastToFloat) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::WriteUDF(const char *name, int  L, LevelAdaptiveData<DIM> *Tmp, bool CastToFloat) {
  int LL[DIM],i ;
  for (i=0; i<DIM; i++) LL[i]=L ;
  WriteUDF(name,LL,Tmp,CastToFloat) ;
}


template<int DIM>
void LevelAdaptiveData<DIM>::WriteUDF(const char *name, UniformData<DIM> *M, LevelAdaptiveData<DIM> *Tmp, bool CastToFloat) {
  int  i,BG[DIM][2],L[DIM] ;
  LevelAdaptiveData<DIM> *X=this ;

  // change of basis to Interpolets without boundary conditions
  if (Tmp) {
    for (i=0; i<DIM; i++) {
      assert(Ext.ismultiscale[i]) ;
      if (Ext.islifting[i]) {
        Tmp->ApplyOp( Ext.BC[i] , X , i, LIFTING2INTERPOLET) ;
        X=Tmp ;
      }
    }
    for (i=0; i<DIM; i++) {
      GetGeneralBoundaryConditions(Ext.BC[i],BG[i]) ;    
      Tmp->ApplyOp( BG[i] , X , i, PROJECTION) ;
      X=Tmp ;
    }
  }

  ToUniform(M) ;
  if (Tmp) for (i=0; i<DIM; i++) M->ApplyOp(BG[i], M, i, INVERSETRANSFORM) ;

  M.WriteUDF(name,NULL,CastToFloat) ;
}



//
// Operators
//
template<int DIM>
void LevelAdaptiveData<DIM>::ApplyOp(  LevelAdaptiveData<DIM> *S , int dir , unsigned int op) {
 ApplyOp( BC[dir] , S, dir , op) ;
}

template<int DIM>
void LevelAdaptiveData<DIM>::ApplyOp(int *BCto , LevelAdaptiveData<DIM> *S , int dir , unsigned int op) {
 // check
 assert(A==S->A) ;
 assert(dir<DIM) ;

 int i,m=0 ;
 // allocate aux. vectors ;
 for (i=0; i<DIM; i++) m=max(m , G->J[i]) ;
 double s[(1<<m)+1] , t[(1<<m)+1] ;

 // forall ProjectionSubSpaces
 Matrix<Matrix<double,DIM> * , DIM>::iterator il(G->Level0,G->J,dir) ;
 for (il.Set(Level0); il<=(*a).end(); ++il) {

   // forall entries in each ProjectionSubSpaces 
   Matrix<double,DIM> *ail=(*a)[il] ;
   if ((*G->a)[il.i]) {
     assert (ail) ;
     Matrix<double,DIM>::iterator it( ail , dir ) ;

     for(it=(*ail).begin(dir); it <= (*ail).end(); ++it) {

       // copy one line
       int is[DIM],ik[DIM],JJ=0 ;
       int n=0     ;
       for (i=0; i<DIM; i++) { is[i]=it.i[i] ; ik[i]=il.i[i] ; }
       for (ik[dir]=0; ik[dir]<=G->J[dir]; ik[dir]++)
         if ((*G->a)[ik]) {
           Matrix<double,DIM> *aik=(*S->a)[ik] ;
           assert(aik) ;
           int e = (ik[dir]==Level0[dir]) ? (1<<Level0[dir]) : (1<<(ik[dir]-1)) -1 ;
           for (is[dir]=0; is[dir]<=e; is[dir]++) s[n++] = (*aik)[is] ;
           JJ=ik[dir] ;
         }

       // apply operation
       switch (op) {
         case STIFFNESS: case FIRSTDERIVATIVE: case FD12: case FD14: case FD22: case FD24: case FD16: 
         case DIV: case GRAD: case PROJECTION: {
           assert(!WC->InterpoletFlag) ;
           ibwt   (t,S->BC[dir] , &m  , s, S->BC[dir] , WC , JJ , JJ-G->Level0[dir] , 0) ;
           applyOp(s,BCto       , NULL, t, S->BC[dir] , WC , JJ , op                , Ext.XE[dir]-Ext.XA[dir]) ;
           bwt    (t,BCto       , &m  , s, BCto       , WC , JJ , JJ-G->Level0[dir] , 0) ;
          } break ;
         
         default: assert(0) ;
       }

       // write back result
       n=0 ;
       for (ik[dir]=0; ik[dir]<=G->J[dir]; ik[dir]++)
         if ((*G->a)[ik]) {
           Matrix<double,DIM> *aik=(*a)[ik] ;
           int e = (ik[dir]==G->Level0[dir]) ? (1<<G->Level0[dir]) : (1<<(ik[dir]-1)) -1 ;
           for (is[dir]=0; is[dir]<=e; is[dir]++) (*aik)[is] = t[n++] ;
         }

       } // next line 
     }   // if ()
   }     // next ProjectionSubSpace

 // Boundary Conditions 
 Ext.CopyFlags(&S->Ext) ;
 Ext.SetBoundaryConditions(dir,BCto[0] , BCto[1]) ;
}


template<int DIM>
void LevelAdaptiveData<DIM>::Blend(BlendingData<DIM> *B) {
 int i,la[DIM],sa[DIM],se[DIM],ta[DIM],te[DIM] ;

 // check
 for (i=0; i<DIM; i++) {
   la[i]=Ext.WC->Level0 ;
   assert(B->Ext.ismultiscale[i]) ;
 }

 // init. 
 Set(0.0) ;
 double w ;

 // forall active BlendingSpaces
 Matrix<Matrix<double,DIM> * , DIM>::iterator k(B->a) ;
 for (k=(*B->a).begin(); k<=(*B->a).end(); ++k) {
   if ( (w=(*B->A->a)[k.i]) != 0) {
     Matrix<double,DIM> *ak=(*B->a)[k] ;
     assert (ak) ;
     
     Matrix<Matrix<double,DIM> * , DIM>::iterator l(la, k.i) ;
     for (l.Set(la);  l<=k.i; ++l) {
       for (i=0; i<DIM; i++)
         if (l.i[i]==la[i]) {
           sa[i]=ta[i]=0        ;
           se[i]=te[i]=1<<la[i] ;
         } else {
           sa[i]=(1<<(l.i[i]-1))+1 ;
           se[i]= 1<<l.i[i]        ;
           ta[i]=0                 ;
           te[i]=se[i]-sa[i]       ;
         }
       
       Matrix<double , DIM>::iterator s(sa, se) ;
       Matrix<double , DIM>::iterator t(ta, te) ;
       Matrix<double , DIM> *al= (*a)[l] ;
       assert(al) ;
       for (s.Set(sa), t.Set(ta); s<=se; ++s,++t) (*al)[t] += w * (*ak)[s] ;

     } // l
   }   // if ()
 }     // next BlendingSpace

 Ext.CopyFlags(&B->Ext) ;
 Ext.CopyData (&B->Ext) ;
}


template<int DIM>
void LevelAdaptiveData<DIM>::Multiply(LevelAdaptiveData<2*DIM> *K , LevelAdaptiveData<DIM> *S) {

  if (Ext.WC->InterpoletFlag) {
    std::cout <<"LevelAdaptiveData<DIM>::Multiply supported for orthogonal Daubechies wavelets only!\n" ;
    assert(0) ;
  }

 int i  ;
 double w ;
 for (i=0; i<DIM; i++) {
   assert(S->Ext.BC[i][0] == K->Ext.BC[DIM+i][0]) ;
   assert(S->Ext.BC[i][1] == K->Ext.BC[DIM+i][1]) ;

   assert(S->Ext.ismultiscale[i]) ;
   assert(K->Ext.ismultiscale[i]) ;
   assert(K->Ext.ismultiscale[DIM+i]) ;

   Ext.ismultiscale[i]=true ;
   Ext.BC[i][0]=K->Ext.BC[i][0] ;
   Ext.BC[i][1]=K->Ext.BC[i][1] ;
 }

 // forall SubSpaces of (*this)
 Matrix<Matrix<double,DIM> * , DIM>::iterator l(a) ;
 for (l=(*a).begin(); l<=(*a).end(); ++l) {
   // forall entries in each SubSpace 
   if ( (*G->a)[l.i] ) {
     Matrix<double,DIM> *al=(*a)[l] ;
     Matrix<double,DIM>::iterator s(al) ;
     for (s=(*al).begin(); s<= (*al).end(); ++s) {
       w=0 ;
       // forall SubSpaces in S with active K(.,S)
       Matrix<Matrix<double,DIM> * , DIM>::iterator Sl(S->a) ;
       for (Sl=(*S->a).begin(); Sl<=(*S->a).end(); ++Sl) { 
         Matrix<double,DIM> *Sal=(*S->a)[Sl] ;
         // joint level for K
         int                 Kl[2*DIM] , Ks[2*DIM] ;
         for (i=0; i<DIM; i++) {
           Kl[i    ]= l.i[i] ;
           Kl[DIM+i]=Sl.i[i] ;
           Ks[i    ]= s.i[i] ;
         }

         if ( (*S->G->a)[Sl.i] && (*K->G->a)[Kl] ) {

           Matrix<double,2*DIM> *Kal=(*K->a)[Kl] ;
           Matrix<double,DIM>::iterator Ss(Sal) ;
           for (Ss=(*Sal).begin(); Ss<= (*Sal).end(); ++Ss) {
              // joint spatial index  
              for (i=0; i<DIM; i++) Ks[DIM+i]=Ss.i[i] ;
              w +=(*Kal)[Ks]*(*Sal)[Ss] ;
           }  // next Ss

         }    // if 
       }      // next Sl

       // Der Lohn der Muehen:
       (*al)[s]=w ;
     }        // next s
   }          // if
 }            // next l

}

#endif

---