Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

UniformData.hpp

#ifndef _DATAMATRIX_
# define _DATAMATRIX_

#include "Matrix.hpp"
#include "NonAdaptive.hpp"
#include "Extensions.hpp"

#ifdef _FFTW_
  #include <rfftw.h>
#endif

struct Function ;

#define LMAXDATAARRAY 12

double _DATAMATRIX_BUFFER_s[(1<<LMAXDATAARRAY)+1],
       _DATAMATRIX_BUFFER_t[(1<<LMAXDATAARRAY)+1],
       _DATAMATRIX_BUFFER_q[(1<<LMAXDATAARRAY)+1] ;

template<int DIM>
struct UniformData {
//
public:
         UniformData()  ;
         UniformData(int *begin,int *end , Wavelets *W,double *XA=NULL, double *XE=NULL) ;
void     Init       (int *begin,int *end , Wavelets *W,double *XA=NULL, double *XE=NULL) ;

        ~UniformData() ;
  void   Clear      () ;
 
  void   SetBoundaryConditions(UniformData<DIM> *X) ;
  void   SetBoundaryConditions(int BC[DIM][2]) ;
  void   SetBoundaryConditions(int dir,int *bc) ;
  void   SetBoundaryConditions(int dir,int bc0,int bc1) ;
  bool   SameBoundaryConditions(UniformData<DIM> *X) ;
  
  // Compatibility Check for multi-array operations like A+B,A+B+C and so on
  void   CheckBC(UniformData<DIM> *X) ;
  void   CheckBC(UniformData<DIM> *X , UniformData<DIM> *Y) ;
  void   CheckBC(UniformData<DIM> *X , UniformData<DIM> *Y , UniformData<DIM> *Z) ;
  void   CheckBC(UniformData<DIM> *X , UniformData<DIM> *Y , UniformData<DIM> *Z , UniformData<DIM> *Q) ;
  void   CheckBC(UniformData<DIM> *X , UniformData<DIM> *Y , UniformData<DIM> *Z , UniformData<DIM> *Q , UniformData<DIM> *R) ;
 
  //
  void   Compatible(UniformData<DIM> *X) ;
  void   Compatible(UniformData<DIM> *X, UniformData<DIM> *Y) ;
  void   Compatible(UniformData<DIM> *X, UniformData<DIM> *Y, UniformData<DIM> *Z) ;
  void   Compatible(UniformData<DIM> *X, UniformData<DIM> *Y, UniformData<DIM> *Z, UniformData<DIM> *Q) ;
  void   Compatible(UniformData<DIM> *X, UniformData<DIM> *Y, UniformData<DIM> *Z, UniformData<DIM> *Q, UniformData<DIM> *R) ;

  void   CopyExtensions(UniformData<DIM> *X) ;

  void   SetFunction(Function *F, bool boundaryonly=false, bool notransform=false) ; 

  double Max()    ;
  double Min()    ;
  double MaxAbs() ;
  double Sum()    ;
  double SumAbs() ;
  double SumSqr() ;
 
  double LpNorm  (UniformData<DIM> *Tmp, double p) ;
  double LmaxNorm(UniformData<DIM> *Tmp) ;
  double L2Norm  (UniformData<DIM> *Tmp) ;
  double L1Norm  (UniformData<DIM> *Tmp) ;
  double Integrate() ;
  double InnerProd(UniformData<DIM> *X) ;

  void   Set(const double a) ;
  void   Copy(UniformData<DIM> *X) ;
  void   Set(UniformData<DIM> *X) ;
  void   Abs(UniformData<DIM> *X) ;
  void   Pow(UniformData<DIM> *X, double p) ;
  void   Sqr(UniformData<DIM> *X)  ;
  
  void   PPlus(UniformData<DIM> *X) ;
  void   PPlus(const double a , UniformData<DIM> *X) ;
  void   Add(UniformData<DIM> *X , UniformData<DIM> *Y) ;
  void   Add(const double a , UniformData<DIM> *X , const double b , UniformData<DIM> *Y) ;
  void   Add(const double a , UniformData<DIM> *X , const double b , UniformData<DIM> *Y , 
             const double c , UniformData<DIM> *Z) ;
  void   Add(const double a , UniformData<DIM> *X , const double b , UniformData<DIM> *Y , 
             const double c , UniformData<DIM> *Z , const double d , UniformData<DIM> *Q) ;
  void   Add(const double a , UniformData<DIM> *X , const double b , UniformData<DIM> *Y , 
             const double c , UniformData<DIM> *Z , const double d , UniformData<DIM> *Q , 
             const double f , UniformData<DIM> *R) ;

  void   MMinus(UniformData<DIM> *X) ;
  void   MMinus(const double a, UniformData<DIM> *X) ;
  void   Sub(UniformData<DIM> *X , UniformData<DIM> *Y) ;

  void   Mul(UniformData<DIM> *X , UniformData<DIM> *Y) ;
  void   Mul(UniformData<DIM> *X , UniformData<DIM> *Y , double f) ;
  void   Mul(UniformData<DIM> *X , double f) ;
  void   Mul(double f) ;
  void   XplusYmalZ (UniformData<DIM> *X , UniformData<DIM> *Y , UniformData<DIM> *Z) ;
  void   XminusYmalZ(UniformData<DIM> *X , UniformData<DIM> *Y , UniformData<DIM> *Z) ;

  // Linear operators
  void   ApplyOp(UniformData<DIM> *X,int dir,unsigned int op) ;
  void   ApplyOp(int *BCT , UniformData<DIM> *X,int dir,unsigned int op) ;
  
  void   SwitchToNoBoundaryConditions(UniformData<DIM> *X) ; // result has in any case [-1,-1] BCs
  // void    RestrictPressure() ;
  void   MRATransform       (UniformData<DIM> *X) ;
  void   MRAInverseTransform(UniformData<DIM> *X) ;
 

  void   FourierCoefficient(double *re, double *im, int *k) ;
  
  // FourierSpectrum, reliable for 0\le k < Kmin, 
  // nonzero                   for 0\le k \le Kmax
  void   FourierSpectrum(double *sp, int *Kmin, int *Kmax) ;

  
  // Calculate mean and variance of slices along the <dir> coordinate direction  
  void   MeanVariance   (double *mean, double *variance, int dir) ;


  void WriteUDF(const char *filename, UniformData<DIM> *Tmp=NULL, bool CastToFloat=false) ;
  void ReadUDF (const char *name) ;
  void ReadSparse (const char *name, int which=0) ;
  void WriteSparse(const char *filename, UniformData<DIM> **M=NULL, int N=0, bool CastToFloat=false) ;

  // special functions
  typedef double (*BoundFunction)(double *x,int d)  ; 
  void SetBoundaryFunction(BoundFunction f) ;

  inline double& operator[](const class Matrix<double,DIM>::iterator &it) { return (*(*this).ba)[it] ;} ;  
  inline double& operator[](const int *in     )                     { return (*(*this).ba)[in] ;} ;

  //internal data
  Matrix<double , DIM> *ba ;
  Extensions<DIM> Ext      ;
} ;


//                              //
// UniformData: Implementations //
//                             //

template<int DIM>
UniformData<DIM>::UniformData() {
  ba=new Matrix<double , DIM> ;
}

template<int DIM>
UniformData<DIM>::UniformData(int *A , int *E , Wavelets *W,double *XA, double *XE) {
  ba=new Matrix<double , DIM>(A,E) ;
  Ext.Init(W,XA,XE) ;
}


template<int DIM>
UniformData<DIM>::~UniformData() {
  delete ba ;
}

template<int DIM>
void UniformData<DIM>::Init(int *A , int *E , Wavelets *W,double *XA, double *XE) {
  ba->Init(A,E)     ;
  Ext.Init(W,XA,XE) ;
}

template<int DIM>
void UniformData<DIM>::Clear() {
  ba->Clear() ;
  Ext.WC=NULL ;
}

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

//
//
template<int DIM>
void UniformData<DIM>::Compatible(UniformData<DIM> *X) {
  if (!Ext.IsSame(&X->Ext)) {
    Ext.Print() ;
    X->Ext.Print() ;
    assert(0) ;
  }
}

template<int DIM>
void UniformData<DIM>::Compatible(UniformData<DIM> *X, UniformData<DIM> *Y) {
 Compatible(X) ;
 Compatible(Y) ;
}

template<int DIM>
void UniformData<DIM>::Compatible(UniformData<DIM> *X, UniformData<DIM> *Y, UniformData<DIM> *Z) {
 Compatible(X) ;
 Compatible(Y) ;
 Compatible(Z) ;
}

template<int DIM>
void UniformData<DIM>::Compatible(UniformData<DIM> *X, UniformData<DIM> *Y, UniformData<DIM> *Z, UniformData<DIM> *Q) {
 Compatible(X) ;
 Compatible(Y) ;
 Compatible(Z) ;
 Compatible(Q) ;
}

template<int DIM>
void UniformData<DIM>::Compatible(UniformData<DIM> *X, UniformData<DIM> *Y, UniformData<DIM> *Z, UniformData<DIM> *Q, UniformData<DIM> *R) {
 Compatible(X) ;
 Compatible(Y) ;
 Compatible(Z) ;
 Compatible(Q) ;
 Compatible(R) ;
}

//
// Check
template<int DIM>
void UniformData<DIM>::CheckBC(UniformData<DIM> *X) {
  assert (Ext.SameBoundaryConditions(&X->Ext)) ;
}

template<int DIM>
void UniformData<DIM>::CheckBC(UniformData<DIM> *X , UniformData<DIM> *Y) 
{ CheckBC(X) ;  CheckBC(Y) ; }

template<int DIM>
void UniformData<DIM>::CheckBC(UniformData<DIM> *X , UniformData<DIM> *Y , UniformData<DIM> *Z) 
{ CheckBC(X) ;  CheckBC(Y) ; CheckBC(Z) ; }

template<int DIM>
void UniformData<DIM>::CheckBC(UniformData<DIM> *X , UniformData<DIM> *Y , UniformData<DIM> *Z , UniformData<DIM> *Q) 
{ CheckBC(X) ;  CheckBC(Y) ; CheckBC(Z) ;  CheckBC(Q) ; }

template<int DIM>
void UniformData<DIM>::CheckBC(UniformData<DIM> *X , UniformData<DIM> *Y , UniformData<DIM> *Z , UniformData<DIM> *Q , UniformData<DIM> *R) 
{ CheckBC(X) ;  CheckBC(Y) ; CheckBC(Z) ; CheckBC(Q) ; CheckBC(R) ; }

//
// BoundaryConditions

template<int DIM>
void UniformData<DIM>::SetBoundaryConditions(UniformData<DIM> *X) {
 Ext.SetBoundaryConditions(X->Ext.BC) ;
}
template<int DIM>
void UniformData<DIM>::SetBoundaryConditions(int bc[DIM][2]) {
 Ext.SetBoundaryConditions(bc) ;
}
template<int DIM>
void UniformData<DIM>::SetBoundaryConditions(int d,int *bc) {
 Ext.SetBoundaryConditions(d,bc) ;
}
template<int DIM>
void UniformData<DIM>::SetBoundaryConditions(int d,int bc0,int bc1) {
 Ext.SetBoundaryConditions(d, bc0,bc1) ;
}
template<int DIM>
bool UniformData<DIM>::SameBoundaryConditions(UniformData<DIM> *X) {
 return Ext.SameBoundaryConditions(&X->Ext) ;
}

//
// Set

template<int DIM>
void  UniformData<DIM>::Set(const double a) { 
 ba->Set(a) ; 
 Ext.Set(0) ;   ;
}

template<int DIM>
void  UniformData<DIM>::SetFunction(Function *F, bool boundaryonly, bool notransform) {
 int i,ii ;

 // test/set valid boundary conditions
 for (i=0; i<DIM; i++) {
   for (int k=0; k<2; k++)
     if ( (Ext.BC[i][k]!=-1) && (Ext.BC[i][k]!=0) && (Ext.BC[i][k]!=1) && (Ext.BC[i][k]!=PERIODIC) ) {
       std::cout<< "UniformData<"<<DIM<<">::SetFunction():: you must set valid boundary conditions first, to define the periodic directions\n" ;  
       assert(0) ;
     }
   Ext.BC[i][0]=-1 ;
   if (Ext.BC[i][1]!=PERIODIC) Ext.BC[i][1]=-1 ;
 }

 Matrix<double,DIM>::iterator s(ba) ;

 for (s.Set(ba->b); s<=ba->e; ++s) {
   double x[DIM] ;
   for (i=0; i<DIM ;i++) {
     ii=s.i[i] - ba->b[i] ;
     x[i]=Ext.XA[i]+((Ext.XE[i]-Ext.XA[i])*ii) / (ba->e[i]-ba->b[i]) ;
    }
   (*ba)[s]=F->Eval(x) ;  
 }
 
 Ext.Set(0)     ;
 for (i=0; i<DIM; i++) Ext.ismultiscale[i]=Ext.islifting[i]=false ; 
 if (!notransform) for (i=0; i<DIM; i++) ApplyOp(Ext.BC[i], this, i, TRANSFORM) ;
 else { 
   if(boundaryonly) { std::cout<< "UniformData<"<<DIM<<">::SetFunction boundaryonly==true && notransform==true not supported yet\n"; assert(0) ;}
  }

 if (boundaryonly) {
   int L[DIM],l ;
   unsigned long t ;
   ba->CalcLevel(L) ;
   
   for (s.Set(ba->b); s<=ba->e; ++s) {
     bool boundary=false ;
     for (i=0; i<DIM; i++) {
       decode(&l,&t,s.i[i],L[i],Ext.WC->Level0) ;
       if ( (l==Ext.WC->Level0) && ((t==0)||(t==(unsigned long)(1<<Ext.WC->Level0))) && (Ext.BC[i][1]!=PERIODIC)) boundary=true ;
     }

     if (!boundary) (*ba)[s]=0 ;
   }  
 }
}


//
// Linear Algebra

template<int DIM>
void  UniformData<DIM>::Copy(UniformData<DIM> *X) {
 SetBoundaryConditions(X) ;
 ba->Copy(X->ba)       ;
 Ext.CopyFlags(&X->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::Set(UniformData<DIM> *X) {
 SetBoundaryConditions(X) ;
 ba->Set(X->ba)        ;
 Ext.CopyFlags(&X->Ext) ;

 // patch grid extensions
 double *XA=X->Ext.XA, *Xa=Ext.XA ,
        *XE=X->Ext.XE, *Xe=Ext.XE ;
 int    *A=X->ba->b, *a=ba->b ,
        *E=X->ba->e, *e=ba->e ;
 assert(XA!=Xa) ;
 assert(XE!=Xe) ; 

 for (int i=0; i<DIM; i++) {
   Xa[i]=XA[i]+((a[i] - A[i])*(XE[i]-XA[i]))/(E[i]-A[i]) ;
   Xe[i]=XA[i]+((e[i] - A[i])*(XE[i]-XA[i]))/(E[i]-A[i]) ;
 }
}

template<int DIM>
void  UniformData<DIM>::Abs(UniformData<DIM> *X) {
 SetBoundaryConditions(X) ;
 ba->Abs(X->ba) ;
 Ext.CopyFlags(&X->Ext) ;
 Ext.Abs(&X->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::Pow(UniformData<DIM> *X, double p) {
 SetBoundaryConditions(X) ;
 ba->Pow(X->ba,p) ;
 Ext.CopyFlags(&X->Ext) ;
 Ext.Pow(&X->Ext,p)    ;
}

template<int DIM>
void  UniformData<DIM>::Sqr(UniformData<DIM> *X) {
 SetBoundaryConditions(X) ;
 ba->Sqr(X->ba) ;
 Ext.CopyFlags(&X->Ext) ;
 Ext.Sqr      (&X->Ext) ;
}


template<int DIM>
double UniformData<DIM>::Max() { 
 return max(ba->Max() , Ext.Max()) ; 
}

template<int DIM>
double UniformData<DIM>::Min() { 
 return min(ba->Min() , Ext.Min()) ;
}

template<int DIM>
double UniformData<DIM>::MaxAbs() {
 return max(ba->MaxAbs() , Ext.MaxAbs()) ; 
}

template<int DIM>
double UniformData<DIM>::Sum() {
 return ba->Sum() + Ext.Sum() ;
}

template<int DIM>
double UniformData<DIM>::SumAbs() {
 return ba->SumAbs() + Ext.SumAbs() ;
}

template<int DIM>
double UniformData<DIM>::SumSqr() {
 return ba->SumSqr() + Ext.SumSqr() ;
}

template<int DIM>
double UniformData<DIM>::InnerProd(UniformData<DIM> *X) { 
 Compatible(X) ;
 return ba->InnerProd(X->ba) + Ext.InnerProd(&X->Ext) ;
}

/*
//
// Norms
template<int DIM>
double UniformData<DIM>::LpNorm(UniformData<DIM> *Tmp, double p) {
 int i,GBC[2] ;
 for (i=0; i<DIM; i++) {
   GetGeneralBoundaryConditions(BC[i],GBC) ;
   Tmp->ApplyOp(GBC, (i==0) ? this : Tmp , i , PROJECTION) ;
 }

 for (i=0; i<DIM; i++) Tmp->ApplyOp(GBC, Tmp , i , INVERSETRANSFORM) ;
 Tmp->Abs(Tmp) ;

 // Maximum Norm
 if (p==HUGE_VAL) return Tmp->Max() ;

 Tmp->Pow(Tmp,p) ;
 for (i=0; i<DIM; i++) Tmp->ApplyOp(GBC, Tmp , i , TRANSFORM) ;

 return pow(Tmp->Integrate(),1/p) ;
}

template<int DIM>
double UniformData<DIM>::LmaxNorm(UniformData<DIM> *Tmp) { return LpNorm(Tmp,HUGE_VAL) ; }

template<int DIM>
double UniformData<DIM>::L2Norm(UniformData<DIM> *Tmp)   { return LpNorm(Tmp,2.0) ; }

template<int DIM>
double UniformData<DIM>::L1Norm(UniformData<DIM> *Tmp)   { return LpNorm(Tmp,1.0) ; }
*/


template<int DIM>
void  UniformData<DIM>::PPlus(UniformData<DIM> *X) { 
 Compatible(X) ;
 ba->PPlus( X->ba)  ;
 Ext.PPlus(&X->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::PPlus(const double a,UniformData<DIM> *X) { 
 Compatible(X) ;
 ba->PPlus(a, X->ba ) ;
 Ext.PPlus(a,&X->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::Add(UniformData<DIM> *X,UniformData<DIM> *Y) { 
 X->Compatible(Y) ;
 CopyExtensions(X) ;
 ba->Add( X->ba  ,  Y->ba) ; 
 Ext.Add(&X->Ext , &Y->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::Add(const double a,UniformData<DIM> *X,const double b,UniformData<DIM> *Y) { 
 X->Compatible(Y) ;
 CopyExtensions(X) ;
 ba->Add(a, X->ba  , b, Y->ba) ;
 Ext.Add(a,&X->Ext , b,&Y->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::Add(const double a,UniformData<DIM> *X,const double b,UniformData<DIM> *Y,
                const double c,UniformData<DIM> *Z) 
{ 
 X->Compatible(Y,Z) ;
 CopyExtensions(X) ;
 ba->Add(a, X->ba  , b, Y->ba  , c, Z->ba) ;
 Ext.Add(a,&X->Ext , b,&Y->Ext , c,&Z->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::Add(const double a,UniformData<DIM> *X,const double b,UniformData<DIM> *Y,
                const double c,UniformData<DIM> *Z,const double d,UniformData<DIM> *Q) 
{ 
 X->Compatible(Y,Z,Q) ;
 CopyExtensions(X) ;
 ba->Add(a, X->ba  , b, Y->ba  , c, Z->ba , d, Q->ba) ;
 Ext.Add(a,&X->Ext , b,&Y->Ext , c,&Z->Ext, d,&Q->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::Add(const double a,UniformData<DIM> *X,const double b,UniformData<DIM> *Y,
                const double c,UniformData<DIM> *Z,const double d,UniformData<DIM> *Q,
                const double e,UniformData<DIM> *R) 
{ 
 X->Compatible(Y,Z,Q,R) ;
 CopyExtensions(X) ;
 ba->Add(a,X->ba,b,Y->ba,c,Z->ba,d,Q->ba,e,R->ba) ;
 Ext.Add(a,&X->Ext , b,&Y->Ext , c,&Z->Ext, d,&Q->Ext , e,&R->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::MMinus(UniformData<DIM> *X) { 
 Compatible(X) ;
 ba->MMinus( X->ba)  ;
 Ext.MMinus(&X->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::MMinus(const double a,UniformData<DIM> *X) { 
 Compatible(X) ;
 ba->MMinus(a,X->ba) ;
 Ext.MMinus(a,&X->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::Sub(UniformData<DIM> *X,UniformData<DIM> *Y) { 
 X->Compatible(Y) ;
 CopyExtensions(X) ;
 ba->Sub( X->ba  ,  Y->ba)  ;
 Ext.Sub(&X->Ext , &Y->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::Mul(UniformData<DIM>     *X ,UniformData<DIM>      *Y) { 
 X->Compatible(Y) ;
 CopyExtensions(X) ;
 ba->Mul( X->ba  ,  Y->ba)  ;
 Ext.Mul(&X->Ext , &Y->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::Mul(UniformData<DIM>      *X ,UniformData<DIM>      *Y , double f) { 
 X->Compatible(Y) ;
 CopyExtensions(X) ;
 ba->Mul( X->ba  ,  Y->ba , f) ;
 Ext.Mul(&X->Ext , &Y->Ext, f) ;
}

template<int DIM>
void  UniformData<DIM>::Mul(UniformData<DIM>      *X , double f) { 
 CopyExtensions(X) ;
 ba->Mul( X->ba , f) ;
 Ext.Mul(&X->Ext, f) ;
}

template<int DIM>
void  UniformData<DIM>::Mul(double f) { 
 ba->Mul(f) ;
 Ext.Mul(f) ;
}

template<int DIM>
void  UniformData<DIM>::XplusYmalZ (UniformData<DIM> *X , UniformData<DIM> *Y , UniformData<DIM> *Z) { 
 X->Compatible(Y,Z) ;
 CopyExtensions(X) ;
 ba->XplusYmalZ(X->ba,Y->ba,Z->ba) ;
 Ext.XplusYmalZ(&X->Ext,&Y->Ext,&Z->Ext) ;
}

template<int DIM>
void  UniformData<DIM>::XminusYmalZ(UniformData<DIM> *X , UniformData<DIM> *Y , UniformData<DIM> *Z) { 
 X->Compatible(Y,Z) ;
 CopyExtensions(X) ;
 ba->XminusYmalZ(X->ba,Y->ba,Z->ba) ;
 Ext.XminusYmalZ(&X->Ext,&Y->Ext,&Z->Ext) ;
}


//
// linear Operators

template<int DIM>
void UniformData<DIM>::ApplyOp(UniformData<DIM> *X,int dir,unsigned int op) {
  ApplyOp(Ext.BC[dir] , X , dir , op) ; 
}

template<int DIM>
void UniformData<DIM>::ApplyOp(int *BCTo , UniformData<DIM> *X,int dir,unsigned int operation) {
 int op,res,J[DIM],dmy ;
 Matrix<double,DIM>::UniDirOperation uniop=NULL ;
 bool Trafo=false ;

 op=NOOPERATION ;
 ba->CalcLevel(J) ;
 ba->SameSize(X->ba) ;

#ifdef _FFTW_
 int  fftN=1<<J[dir]  ;
 fftw_real *fftIn = (fftw_real *)malloc(sizeof(fftw_real)*fftN);
 fftw_real *fftOut= (fftw_real *)malloc(sizeof(fftw_real)*fftN);
 rfftw_plan fftp=NULL ;
#endif

 CopyExtensions(X) ;
 double DX=Ext.XE[dir]-Ext.XA[dir] ;

 // the fully consistent AFD2 algorithms coincide with standard algorithms in the non-adaptive case.
 if (operation==FD14II) operation=FD14 ;
 if (operation==FD24II) operation=FD24 ;

 switch (operation) {
 case FIRSTDERIVATIVE: case STIFFNESS: case DIV: case GRAD: case FD12: case FD14: case FD22: case FD24: 
 case DIVGRADFD4NOSTAB: case DIVGRAD: {
   assert(!Ext.islifting[dir]) ;
   uniop=applyOp               ;
   op   =operation             ;
   Trafo=Ext.ismultiscale[dir] ;
 }
 break ;
 case PROJECTION: {
   uniop=Projection ;
   op   =PROJECTION ;
   Trafo=Ext.ismultiscale[dir] ;
 }
 break ;
 case TRANSFORM: {
   assert(!Ext.ismultiscale[dir]) ;
   Ext.ismultiscale[dir]=true ;
   uniop=bwt ;
   op   =J[dir]-Ext.WC->Level0 ;
 }
 break ;
 case INVERSETRANSFORM: {
   assert(Ext.ismultiscale[dir]) ;
   Ext.ismultiscale[dir]=false   ;
   uniop=ibwt ;
   op   =J[dir]-Ext.WC->Level0 ;
 }
 break ;
 case INTERPOLET2LIFTING: {
   assert(!Ext.islifting[dir]) ;
   assert(Ext.WC->InterpoletFlag)  ;
   if (!Ext.ismultiscale[dir]) { Copy(X); Ext.islifting[dir]=true; return; }
 }
 break ;
 case LIFTING2INTERPOLET: {
   assert(Ext.islifting[dir])  ;
   assert(Ext.WC->InterpoletFlag)  ;
   if (!Ext.ismultiscale[dir]) { Copy(X); Ext.islifting[dir]=false; return; }
 }
 break ;
 case POLYQUAD: case INVERSEPOLYQUAD: case NODALQUAD: case INVERSENODALQUAD: case ONEPOINTQUAD: case INVERSEONEPOINTQUAD: case SHORTQUAD: case INVERSESHORTQUAD: {
   assert(!Ext.ismultiscale[dir]) ;
   uniop=Quadrature ;
   op   =operation  ;
 }
 break ;
 case SMOOTH0: case SMOOTH1: case SMOOTH2: {
   assert(!Ext.ismultiscale[dir]) ;
   uniop=applyOp   ;
   op   =operation ;
 }
 break ;
 case SORTHB: {
   uniop=SortHB     ;
   op   =Ext.WC->Level0 ;
 }
 break ;
 case FFT:  {
   assert(!Ext.ismultiscale[dir]) ;
   assert(Ext.BC[dir][1]==PERIODIC) ;
   Ext.ismultiscale[dir]=true ;
#ifdef _FFTW_
   fftp = rfftw_create_plan(fftN, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
#endif
 }
 break ;
 case IFFT: {
   assert(Ext.ismultiscale[dir]) ;
   assert(Ext.BC[dir][1]==PERIODIC) ;
   Ext.ismultiscale[dir]=false   ;
#ifdef _FFTW_
   fftp = rfftw_create_plan(fftN, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
#endif
 }
 break ;
 default: { std::cout<< "UniformData<"<<DIM<<">::ApplyOp(): unsupported operation "<<operation<<"\n"; assert (0) ;}
 }
 
 //
 // main loop along coordinate <dir>
 //
 Matrix<double,DIM>::iterator it,Xend=(*ba).end() ;
 int b=ba->b[dir] , e=ba->e[dir] ;
 double *s=_DATAMATRIX_BUFFER_s, *t=_DATAMATRIX_BUFFER_t, *q=_DATAMATRIX_BUFFER_q ;
 
 for (it=(*ba).begin(dir) ; it<=Xend ; ++it) {

   int *i=&it.i[dir] ;
   for ((*i)=b; (*i)<=e; (*i)++) s[(*i)-b]=(*X->ba)[it] ;

   switch (operation) {
   case FFT: case IFFT: {
#ifdef _FFTW_
     for (int j=0; j<fftN; j++) fftIn[j]=s[j] ;
     rfftw_one(fftp, fftIn, fftOut);
     for (int j=0; j<fftN; j++) t[j]=fftOut[j] ;
#endif
   }
   break ;
   case INTERPOLET2LIFTING: case LIFTING2INTERPOLET: {
     Ext.WC->LiftingFlag=Ext.islifting[dir] ;
     ibwt (q,X->Ext.BC[dir], NULL , s,X->Ext.BC[dir] , Ext.WC, J[dir] , J[dir]-Ext.WC->Level0, DX) ;
     Ext.WC->LiftingFlag = !Ext.WC->LiftingFlag ;
     bwt  (t,BCTo , &dmy , q,BCTo , Ext.WC, J[dir] , J[dir]-Ext.WC->Level0, DX) ;
   }
   break ;
   default: { // all other operations 
     Ext.WC->LiftingFlag=Ext.islifting[dir] ;
     if (!Trafo) {
       uniop(t,BCTo , &res , s,X->Ext.BC[dir] , Ext.WC, J[dir] , op , DX) ;
     } else {
       ibwt (t,X->Ext.BC[dir], NULL , s,X->Ext.BC[dir] , Ext.WC, J[dir] , J[dir]-Ext.WC->Level0, DX) ;
       if (op==DIVGRADFD4NOSTAB) {
         applyOp(q,BCTo , &res , t,X->Ext.BC[dir] , Ext.WC, J[dir] , FD14, DX) ;
         applyOp(s,BCTo , &res , q,X->Ext.BC[dir] , Ext.WC, J[dir] , FD14, DX) ;
       }  
       else {
         uniop(s,BCTo , &res , t,X->Ext.BC[dir] , Ext.WC, J[dir] , op, DX) ;
       }
       bwt  (t,BCTo , &dmy , s,BCTo , Ext.WC, J[dir] , J[dir]-Ext.WC->Level0, DX) ;
     }
   } 
   break ;
   }
   
   for ((*i)=b; (*i)<=e; (*i)++) (*ba)[it]=t[(*i)-b] ;
 } 
 
 Ext.SetBoundaryConditions(dir, BCTo) ;

 //
 // clean up
 //
 switch (operation) {
 case FFT: case IFFT: {
   Ext.ismultiscale[dir] = (operation==FFT);
#ifdef _FFTW_
   free(fftIn) ;
   free(fftOut) ;
   rfftw_destroy_plan(fftp);
#endif
 }
 break ;
 case INTERPOLET2LIFTING: case LIFTING2INTERPOLET: {
   Ext.islifting[dir] = !Ext.islifting[dir] ;
 }
 break ;
 default: ;
 }
 
}

template<int DIM>
void UniformData<DIM>::SwitchToNoBoundaryConditions(UniformData<DIM> *X) {
  int dir,i,s1[DIM],s2[DIM] ;
  if (this != X) Copy(X) ;
  
  for (dir=0; dir<DIM; dir++) {
    Matrix<double,DIM>::iterator s(ba , dir) ;
    if (Ext.BC[dir][0]==0) {
      for (s=(*ba).begin(); s<=(*ba).end(); ++s) {
        for (i=0; i<DIM; i++) s1[i]=s.i[i] ;
        s1[dir]=0 ;
        (*ba)[s1]=0 ;
      }
    }
    
    if (Ext.BC[dir][0]==1) {
      for (s=(*ba).begin(); s<=(*ba).end(); ++s) {
        for(i=0; i<DIM; i++) s1[i]=s2[i]=s.i[i] ;
        s1[dir]=1 ; 
        s2[dir]=0 ;
        (*ba)[s1]=(*ba)[s2] ;
      }
    }

    if (Ext.BC[dir][1]==0) {
      for (s=(*ba).begin(); s<=(*ba).end(); ++s) {
        for (i=0; i<DIM; i++) s1[i]=s.i[i] ;
        s1[dir]=ba->e[dir] ;
        (*ba)[s1]=0 ;
       }
    }

    if (Ext.BC[dir][1]==1) {
      for (s=(*ba).begin(); s<=(*ba).end(); ++s) {
        for(i=0; i<DIM; i++) s1[i]=s2[i]=s.i[i] ;
        s1[dir]=ba->e[dir]-1 ; 
        s2[dir]=ba->e[dir] ;
        (*ba)[s1]=(*ba)[s2] ;
      }
    }

    if (Ext.BC[dir][1]==PERIODIC) {
      for (s=(*ba).begin(); s<=(*ba).end(); ++s) {
        for(i=0; i<DIM; i++) s1[i]=s2[i]=s.i[i] ;
        s1[dir]=ba->e[dir]; s2[dir]=0 ;
        (*ba)[s1]=(*ba)[s2] ;
      }
    }

   Ext.BC[dir][0]=Ext.BC[dir][1]=-1 ;
  }
}



template<int DIM>
double UniformData<DIM>::Integrate() {
 // Integrate of multi / single scale function ?
 int i,J[DIM] ;
 bool ms=Ext.ismultiscale[0] ;
  
 // only complete multi- or single scale allowed
 for (i=1; i<DIM; i++) assert (ms==Ext.ismultiscale[i]) ;

 ba->CalcLevel(J) ;
 double r=0.0 ;
 if (ms) { // multiscale

   int A[DIM],a[DIM],e[DIM] ;
   for (i=0; i<DIM; i++) A[i]=Ext.WC->Level0 ;
  
   Matrix<int,DIM>::iterator l(A,J) ;
   for (l.Set(A); l<=J; ++l) {
     for (i=0; i<DIM; i++) {
       a[i] = (l.i[i]==Ext.WC->Level0) ? 0 : (1<<(l.i[i]-1))+1 ;
       e[i] = 1<<l.i[i] ;
     }
     Matrix<double,DIM>::iterator s(a,e) ;
     for (s.Set(a); s<=e; ++s) {
       double t=1.0 ;
       for (i=0; i<DIM; i++) t *=Ext.WC->WaveletIntegral(l.i[i],s.i[i]-a[i],Ext.BC[i]) ;
       r +=(*ba)[s]*t ;
     } // s
   } // l

 } else {  // single scale
 
   Matrix<double,DIM>::iterator s(ba) ;
   for (s=(*ba).begin(); s<= (*ba).end(); ++s) {
     double t=1.0 ;
     for (i=0; i<DIM; i++) t *=Ext.WC->ScalingFunctionIntegral(J[i],s.i[i],Ext.BC[i]) ;

     r +=(*ba)[s]*t ;
   }
 }

return r ;
}


//
// Norms
template<int DIM>
double UniformData<DIM>::LpNorm(UniformData<DIM> *Tmp, double p) {
 int i,GBC[2] ;
 for (i=0; i<DIM; i++) {
   GetGeneralBoundaryConditions(Ext.BC[i],GBC) ;
   Tmp->ApplyOp(GBC, (i==0) ? this : Tmp , i , PROJECTION) ;
 }

 for (i=0; i<DIM; i++)
   if (Tmp->Ext.islifting[i]) Tmp->ApplyOp(GBC, Tmp , i , LIFTING2INTERPOLET) ;

 for (i=0; i<DIM; i++) Tmp->ApplyOp(GBC, Tmp , i , INVERSETRANSFORM) ;
 Tmp->Abs(Tmp) ;

 // Maximum Norm
 if (p==HUGE_VAL) return Tmp->Max() ;

 Tmp->Pow(Tmp,p) ;
 return pow(Tmp->Integrate(),1/p) ;
}

template<int DIM>
double UniformData<DIM>::LmaxNorm(UniformData<DIM> *Tmp) { return LpNorm(Tmp,HUGE_VAL) ; }

template<int DIM>
double UniformData<DIM>::L2Norm(UniformData<DIM> *Tmp)   { return LpNorm(Tmp,2.0) ; }

template<int DIM>
double UniformData<DIM>::L1Norm(UniformData<DIM> *Tmp)   { return LpNorm(Tmp,1.0) ; }


//
//  isotropic transforms
//
template<int DIM>
void  UniformData<DIM>::MRATransform(UniformData<DIM> *X) {
  int L[DIM],l,dir,dmy ;

  ba->SameSize (X->ba) ;
  for (dir=0; dir<DIM; dir++) assert(!X->Ext.ismultiscale[dir]) ;

  ba->CalcLevel(L) ;
  CopyExtensions(X) ;
  for (dir=0; dir<DIM; dir++)  Ext.ismultiscale[dir] = true ;
  Ext.SetBoundaryConditions(X->Ext.BC) ;

  for (l=L[0]; l>=Ext.WC->Level0+1; l--) {
    int e[DIM] ;
    for (dir=0; dir<DIM; dir++) e[dir]=1<<l ;

    for (dir=0; dir<DIM; dir++) {
      double DX=Ext.XE[dir]-Ext.XA[dir] ;
      Ext.WC->LiftingFlag=Ext.islifting[dir] ;
      
      Matrix<double,DIM>::iterator it ;
      int b=ba->b[dir] ;

      UniformData<DIM> *Y = ((l==L[0])&&(dir==0)) ? X : this ;

      for (it=(*ba).begin(dir); it<=e; ++it) {
        double *s=_DATAMATRIX_BUFFER_s, *t=_DATAMATRIX_BUFFER_t ;

        int *i=&it.i[dir] ;
        for ((*i)=b; (*i)<=e[dir]; (*i)++) s[(*i)-b]=(*Y->ba)[it] ;
        bwt  (t,Ext.BC[dir], &dmy , s,Ext.BC[dir], Ext.WC, l , 1 , DX) ;
        for ((*i)=b; (*i)<=e[dir]; (*i)++) (*ba)[it]=t[(*i)-b] ;
      } // it
    } // dir
  } // l 
}


template<int DIM>
void  UniformData<DIM>::MRAInverseTransform(UniformData<DIM> *X) {
  int L[DIM],l,dir,dmy ;

  ba->SameSize (X->ba) ;
  for (dir=0; dir<DIM; dir++) assert(X->Ext.ismultiscale[dir]) ;

  ba->CalcLevel(L) ;
  CopyExtensions(X) ;
  for (dir=0; dir<DIM; dir++)  Ext.ismultiscale[dir] = false ;
  Ext.SetBoundaryConditions(X->Ext.BC) ;

  for (l=Ext.WC->Level0+1; l<=L[0]; l++) {
    int e[DIM] ;
    for (dir=0; dir<DIM; dir++) e[dir]=1<<l ;

    for (dir=0; dir<DIM; dir++) {
      double DX=Ext.XE[dir]-Ext.XA[dir] ;
      Ext.WC->LiftingFlag=Ext.islifting[dir] ;
      
      Matrix<double,DIM>::iterator it ;
      int b=ba->b[dir] ;

      UniformData<DIM> *Y = ((l==L[0])&&(dir==0)) ? X : this ;

      for (it=(*ba).begin(dir); it<=e; ++it) {
        double *s=_DATAMATRIX_BUFFER_s, *t=_DATAMATRIX_BUFFER_t ;

        int *i=&it.i[dir] ;
        for ((*i)=b; (*i)<=e[dir]; (*i)++) s[(*i)-b]=(*Y->ba)[it] ;
        ibwt (t,Ext.BC[dir], &dmy , s,Ext.BC[dir], Ext.WC, l , 1 , DX) ;
        for ((*i)=b; (*i)<=e[dir]; (*i)++) (*ba)[it]=t[(*i)-b] ;
      } // it
    } // dir
  } // l 
}



// 
//  IO & UDF
//

template<int DIM>
void UniformData<DIM>::ReadSparse(const char *filename, int which) {
  ba->ReadSparse(filename,which,Ext.XA,Ext.XE,Ext.WC->Level0) ;
}

template<int DIM>
void UniformData<DIM>::WriteSparse(const char *filename, UniformData<DIM> **M, int N, bool CastToFloat) {
  int i ;
  Matrix<double,DIM> *MM[100] ;
  for (i=0; i<N ; i++) MM[i]=M[i]->ba ;
  
  ba->WriteSparse(filename,Ext.XA, Ext.XE , Ext.WC->Level0, MM, N, CastToFloat) ;
}


template<int DIM>
void UniformData<DIM>::WriteUDF(const char *name, UniformData<DIM> *Tmp, bool Cast2Float) {
  UniformData<DIM> *A=this ;
  if (Tmp!=NULL) { // transform to nodal values whenever necessary
    for (int i=0; i<DIM; i++)
      if (Ext.ismultiscale[i]) {
        Tmp->ApplyOp(Ext.BC[i] , A , i , INVERSETRANSFORM) ;
        A=Tmp ;
      }
  }
  A->ba->WriteUDF(name,&A->Ext, "scalar" , NULL , 1, Cast2Float, false) ;
}

template<int DIM>
void UniformData<DIM>::ReadUDF(const char *name) {
  ba->ReadUDF(name,&Ext, "scalar") ;
}


//
//  FOURIER STUFF
//

template<int DIM>
void UniformData<DIM>::FourierSpectrum(double *sp, int *Kmin,int *Kmax) {
  int i,j,N[DIM],anz=1 ;
  
  // dimensions, number of elements, size of domain, aspect ratio 
  double ar[DIM],mdx=1e+100,DX[DIM],vol=1.0 ;
  for (i=0; i<DIM; i++) {
    N[i] =(ba->e[i] - ba->b[i]) ;
    anz *=N[i] ;
    
    DX[i]=Ext.XE[i]-Ext.XA[i] ;
    vol *=DX[i] ;
    mdx  =min(mdx,DX[i]) ;
  }

  for (i=0; i<DIM; i++) ar[i]=DX[i]/mdx ;
  
  // find smallest reliable wave number
  j=10000 ;
  for (i=0; i<DIM; i++) j=min(j,(int)(N[i]/ar[i])) ;
  *Kmin=j ;

  // now the FourierSpectrum 
  int a[DIM],e[DIM] ;
  double mxf=0 ;
  for (i=0; i<DIM; i++) {
    a[i] =-(N[i]/2 -1) ; // lowest  frequency
    e[i] =  N[i]/2     ; // highesy frq
    mxf += sqr(e[i]/ar[i]) ;
  }
  mxf=sqrt(mxf)  ;
  *Kmax=(int)mxf ;
  for (i=0; i<=*Kmax; i++) sp[i]=0 ;
  
  Matrix<double,DIM>::iterator s(a,e) ;
  for (s.Set(a); s<=e; ++s) {
    // wavenumber
    double k=0,re,im ;
    for (i=0; i<DIM; i++) k += sqr(s.i[i]/ar[i]) ;

    FourierCoefficient(&re,&im,s.i)     ;
    sp[(int)sqrt(k)] += (re*re + im*im) ;
  }
 
  // remove wrong scaling
  for (i=0; i<=*Kmax; i++) {
   sp[i] *= vol/((double)anz*anz) ;
  }
}


//
// access to Fouriercoefficients is a little bit complicated, due to
// the special storage scheme of the real->complex FFT of the fftw
//
// since real an imaginary parts are stored separately from each other and
// then are seperately transformed with respect to the other directions
// one has to collect the real an imaginary parts from all 2^D qudrants
// plus: one has to care about different signs with wich these 
// terms contribute to the final value
//

template<int DIM>
void UniformData<DIM>::FourierCoefficient(double *re, double *im, int *k) {
  int kk[DIM],N[DIM],q=0,i,d2=1<<DIM ;
  
  // determine quadrant and index for quadrant 0 
  for (i=0; i<DIM; i++) {
    N[i]=(ba->e[i] - ba->b[i]) ;
    if (k[i] < 0) { kk[i]= -k[i], q |= (1<<i) ;}
    else          { kk[i]=  k[i] ; }
  }
   
  // collect all entries from the various subquadrants
  double rj[1<<DIM] ;
  int rjq ;
  for (rjq=0; rjq<d2; rjq++) {
    int kq[DIM] ;
    for (i=0; i<DIM; i++) {
      if (rjq & (1<<i) ) kq[i]=N[i]-kk[i] ;
      else               kq[i]=     kk[i] ;
    }
    rj[rjq]=(*ba)[kq] ;
  } 
 
  // delete spurious entries
  for (i=0 ; i< DIM; i++) {
    if ((kk[i]==0)||(kk[i]==N[i]/2)) 
      for (rjq=0; rjq<d2; rjq++)
        if (rjq & (1<<i)) rj[rjq]=0 ;
  }

  // sum the different entries 
  *re=*im=0 ;
  for (rjq=0; rjq<d2; rjq++) {

    // nj = number of active bits in rjq == power of sqrt(-1) -> if nj is even: rj[rjq] is real
    //                                                           if nj is odd : rj[rjq] is imaginary
    // nj/2 gives the power of (-1) with wich rj[rjq] contributes to the real/imaginary part
    //   
    // number of coinciding bits in rjq and q == number of multiplys by -1 

    int nj=0,nq=0 ;
    for (i=0; i<DIM; i++) {
      if ( rjq &      (1<<i)) nj++ ;
      if ((rjq & q) & (1<<i)) nq++ ;
    }

    // get signs    
    double sg,sq ;
    i=nj/2 ;
    sg = (i  & 1) ? -1 : 1 ;
    sq = (nq & 1) ? -1 : 1 ;

    if (nj&1) { // nj odd -> rj[rjq] contributes to imaginary part
       *im -= rj[rjq]*sg*sq ;
    } else {  // nj even -> rj[rjq] contributes to real part 
       *re += rj[rjq]*sg*sq ;
    }
  }
}



//
//  MEAN and VARIANCE
//
template<int DIM>
void UniformData<DIM>::MeanVariance(double *mean, double *variance, int dir) {
  int i,j,e=(*ba).size(dir),size ;
    // calculate mean values along a x-z plane
  for (j=0; j<e; j++) {
    double u=0,u2=0,x;

    Matrix<double,DIM>::iterator s(ba , dir) ;
    size=0 ;
    for (s=(*ba).begin(dir); s<=(*ba).end(); ++s) {
      int t[DIM] ;
      for (i=0; i<DIM; i++) t[i]=s.i[i] ;
      t[dir]=j  ;
      x   = (*this)[t] ;
      u  += x   ;
      u2 += x*x ;
      size++ ;
    }
    u  /= size ;
    u2 /= size ;
    mean    [j]  = u ;
    variance[j]  = u2 - u*u ;
  }
}

#endif

---