Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

Operators.hpp

//
// Operators
//

#ifndef _Operators_
# define _Operators_
# include <math.h>

#include "AdaptiveData.hpp"
#include "UniformData.hpp"
#include "LevelAdaptive.hpp"

//double RealResidual(Data *) ;

extern int debugRefine ;

template<class I , int DIM>
struct Operator {
  virtual void   Apply                (I *X , I *Result)=0 ;
  virtual void   Preconditioner       (I *X , I *Result)=0 ;
  virtual void   InversePreconditioner(I *X , I *Result)=0 ;
  virtual double DiagonalEntry        (int *level , Wavelets *W)=0 ;
  virtual void   NotUsing             (I *X)=0 ;
} ;

//
// Attention: in the following classes Tmp points to to some auxiliary data used for the evaluation of the operators
// Tmp may be uninitialized. Note, this moves the memory management completely to the calling program.
// In particulary one has to assert that the calling program/subroutine does not use the Operator's (*Tmp) Data.
// It is strongly recommended to verify this condition using the NotUsing() function
// 

// scalar: par[0]*id + par[1]*dxx  + par[2]*dyy + ...
//template<int DIM>
//void Helmholtz(AdaptiveData<DIM> *X, double *par, AdaptiveData<DIM> *Result, AdaptiveData<DIM> *Tmp, AdaptiveData<DIM> *Tmp2) ;

//template<int DIM>
//void Helmholtz(UniformData<DIM>  *X, double *par, UniformData<DIM>  *Result, UniformData<DIM>  *Tmp, UniformData<DIM>  *Tmp2) ;

template<int DIM>
void HelmholtzDiagonalPreconditioner(AdaptiveData<DIM> *X , AdaptiveData<DIM> *R, Operator<AdaptiveData<DIM>,DIM> *O , bool inverse) ;

template<int DIM>
void HelmholtzDiagonalPreconditioner(UniformData<DIM> *X , UniformData<DIM> *R, Operator<UniformData<DIM>,DIM> *O , bool inverse) ;

template<class I , int DIM>
void InverseHelmholtzDiagonalPreconditioner(I *X , I *Result , Operator<AdaptiveData<DIM>,DIM> *O) ;

// To solve singular linear systems one may augment the operator
// for singular problems one may also provide a vector which spans the null space (or an approximation to)
// for the solution of Dirichlet-problems one can define faces on which Dirichlet conditions shall hold;
// to this end one has to set the attribute BC[i][k]:=0 OF THE OPERATOR ! <br>
// this forces the Apply(*X,*R) function just to copy the coefficients for boundary wavelets
// from *X to *R
template<class I , int DIM>
struct HelmholtzOperator : public Operator<I,DIM> {

                 HelmholtzOperator() ;
  virtual void   Apply         (I *X , I *Result) ;
  virtual void   Preconditioner(I *X , I *Result) ;
  virtual void   InversePreconditioner(I *X , I *Result)  ;
  virtual double DiagonalEntry (int *level , Wavelets *W) ;
  virtual void   NotUsing      (I *X) ;

  double  par[DIM+1] ;
  // (optional) vector which spanns the null space 
  I      *Null       ; 

  I     *Tmp  ;
  I     *Tmp2 ;

  bool   use_lifting_preconditioner ;

} ;

template<class I , int DIM>
struct SpaceTimeOperator : public Operator<I,DIM> {

                 SpaceTimeOperator() ;
  virtual void   Apply         (I *X , I *Result) ;
  virtual void   Preconditioner(I *X , I *Result) ;
  virtual void   InversePreconditioner(I *X , I *Result)  ;
  virtual double DiagonalEntry (int *level , Wavelets *W) ;
  virtual void   NotUsing      (I *X) ;

          void   DiagonalScale (I *X , bool inverse) ;
  I     *Tmp, *Tmp2  ;

  bool  use_lifting_preconditioner ; 
} ;



template<class I , int DIM>
struct IntegroDifferentialOperator : public Operator<I , DIM> {
  double par[3]          ;
  Operator<I , DIM> *D   ;
  LevelAdaptiveData<2*DIM>  *K   ;
  I                 *Tmp ;

                 IntegroDifferentialOperator() ;
  virtual void   Apply         (I *X , I *Result) ;
  virtual void   Preconditioner(I *X , I *Result) ;
  virtual void   InversePreconditioner(I *X , I *Result  ) ;
  virtual double DiagonalEntry (int *level  , Wavelets *W) ;
  virtual void   NotUsing      (I *X) ;
} ;


//
// IntegroDifferentialOperator 
//
// Pa[0]*u(x) + Pa[1]*(Du)(x) + Pa[2]*int_y KY(x,y)u(y)dy KX(x)
//
template<class I , int DIM>
struct TensorIntegroDifferentialOperator : public Operator<I , DIM> {
  double par[3]           ; // parameters
  Operator<I , DIM> * D   ; // Differential Operator
  LevelAdaptiveData<DIM>  *KX,*KY ; // Kernel of Integral operator
  I                  *Tmp ; // Only required if D != NULL               

                 TensorIntegroDifferentialOperator() ;
  virtual void   Apply         (I *X , I *Result) ;
  virtual void   Preconditioner(I *X , I *Result) ;
  virtual void   InversePreconditioner(I *X , I *Result  ) ;
  virtual double DiagonalEntry (int *level  , Wavelets *W) ;
  virtual void   NotUsing      (I *X) ;
} ;


template<class I , int DIM>
struct PressureOperatorI : public Operator<I,DIM> {
  I    *Tmp[2]      ;
  I    *Null        ;
  bool  do_patch    ;
  bool  use_lifting_preconditioner ;

                 PressureOperatorI() ;
  virtual void   Apply         (I *X , I *Result) ;
  virtual void   Preconditioner(I *X , I *Result) ;
  virtual void   InversePreconditioner(I *X , I *Result) ;
  virtual double DiagonalEntry (int *level , Wavelets *W) ;
  virtual void   NotUsing      (I *X) ;
} ;


//
//  Null is approximation of the (right) null space of the operator
//  In Apply: if (X->extended) then
//  R->ext = Null->InnerProd(X) ;
//  and R->PPlus(X->ext,Null)   ;
//  hence, the sattlepoint formulation of the linear operator is evaluated
//  which overcomes the difficulties in a lacking solvability constraint
//  for singular linear systems
template<class I , int DIM>
struct PressureOperatorII : public Operator<I,DIM> {
  I    *Tmp[2]      ;
  I    *Null        ;
  bool  do_restrict ;
  bool  do_patch    ;
  bool  use_lifting_preconditioner ;

                 PressureOperatorII() ;
  virtual void   Apply         (I *X , I *Result) ;
  virtual void   Preconditioner(I *X , I *Result) ;
  virtual void   InversePreconditioner(I *X , I *Result) ;
  virtual double DiagonalEntry (int *level , Wavelets *W) ;
  virtual void   NotUsing      (I *X) ;
} ;

template<class I , int DIM>
struct PressureOperatorIII : public Operator<I,DIM> {
  I    *Tmp[2]      ;
  I    *Null        ;
  bool  do_restrict ;
  bool  do_patch    ;
  bool  use_lifting_preconditioner ;

                 PressureOperatorIII() ;
  virtual void   Apply         (I *X , I *Result) ;
  virtual void   Preconditioner(I *X , I *Result) ;
  virtual void   InversePreconditioner(I *X , I *Result) ;
  virtual double DiagonalEntry (int *level , Wavelets *W) ;
  virtual void   NotUsing      (I *X) ;
} ;

//
// div*grad
// if !periodic BC: Vel_BC=(-1,-1),(-1,-1) , Press_BC=(1,1),(1,1)
// Null is approximation of the (right) Null space of the operator
// In Apply: if (X->extended) then
// R->ext = Null->InnerProd(X) ;
// and R->PPlus(X->ext,Null)   ;
// hence, the sattlepoint formulation of the linear operator is evaluated
// which overcomes the difficulties in a lacking solvability constraint
// for singular linear systems
//   
template<class I , int DIM>
struct PressureOperatorIV : public Operator<I,DIM> {
  I    *Tmp[2]      ;
  I    *Null        ;
  bool  do_restrict ;
  bool  do_patch    ;
  bool  use_lifting_preconditioner ;

                 PressureOperatorIV() ;
  virtual void   Apply         (I *X , I *Result) ;
  virtual void   Preconditioner(I *X , I *Result) ;
  virtual void   InversePreconditioner(I *X , I *Result) ;
  virtual double DiagonalEntry (int *level , Wavelets *W) ;
  virtual void   NotUsing      (I *X) ;
} ;


template<class I , int DIM>
struct PressureOperatorV : public Operator<I,DIM> {
  I    *Tmp         ;
  I    *Null        ;
  bool  do_restrict ;
  bool  do_patch    ;
  bool  use_lifting_preconditioner ;

                 PressureOperatorV() ;
  virtual void   Apply         (I *X , I *Result) ;
  virtual void   Applyi        (I *X , I *Result,int i) ;
  virtual void   Preconditioner(I *X , I *Result) ;
  virtual void   InversePreconditioner(I *X , I *Result) ;
  virtual double DiagonalEntry (int *level , Wavelets *W) ;
  virtual void   NotUsing      (I *X) ;
} ;

//
// Convective term:
// Apply:                 scalar (grad X)*U
// ApplyConservative:     scalar div(UX)
// ApplyConservativeFast: scalar div(UU[i]) 
//
template<class I , int DIM>
struct ConvectionOperator : public Operator<I,DIM> {
                
                 ConvectionOperator() ;
          void   ApplyConservativeWENO(I *X  , I *Result) ;
          void   ApplyConservativeWENOUniformData(UniformData<DIM> *X  , UniformData<DIM> *Result, bool trafo=false) ;
  virtual void   Apply                (I *X  , I *Result) ;
          void   ApplyFD              (I *X  , I *Result) ;
          void   ApplyConservative    (I *X  , I *Result) ;
          void   ApplyConservativeFD  (I *X  , I *Result) ;
          void   ApplyConservativeFast(int i , I *Result) ;

  virtual void   Preconditioner       (I *X , I *Result) ;
  virtual void   InversePreconditioner(I *X , I *Result) ;
  virtual double DiagonalEntry        (int *level , Wavelets *W) ;
  virtual void   NotUsing             (I *X) ;

  I     *U[DIM]  ;
  I     *Tmp,*Tmp1 ;
} ;


//                            //
// Operators: Implementations //
//                            //

//
// Helmholtz
//
// Helmholtz evaluates the general Helmholtz-Operator
//
// Pa[0]*u + Pa[1]*uxx + Pa[2]*uyy + Pa[3]*uzz + ....
// projection on basis WITH SAME boundary conditions as X
//

template<class I,int DIM>
void Helmholtz(I *X , double *par, I *R, I *tmp, I *tmp2)
{
 assert(tmp) ;
 assert(tmp!=R) ;
 assert(tmp!=X) ;

 int  i,last,op=STIFFNESS ;
 int  BG[DIM][2]   ;
 bool change=false ;
 bool lift[DIM]    ;
 I    *tmp3 =X     ;

 last=-1 ;
 for (i=0; i<=DIM; i++) if (par[i]!=0.0) last=i ;

 if (last <0) { // all parameters vanish ...
   R->Set(0.0) ;
   R->SetBoundaryConditions(X) ;
   return ;
 }
 if (last==0) { // just the first parameter not ...
   R->Mul(X,par[0]) ;
   R->SetBoundaryConditions(X) ;
   return ;
 }

 // transform to interpolets if necessary: in any case *tmp3 points to the (modified) argument "X"
 for (i=0;i<DIM; i++) {
   if (X->Ext.islifting[i]) {
     assert(tmp2) ;
     tmp2->ApplyOp(X->Ext.BC[i], tmp3, i , LIFTING2INTERPOLET) ;
     tmp3   =tmp2 ;
     lift[i]=true ;
   } else 
     lift[i]=false ;
 }
 // here: tmp3 = X or tmp2 


 // FD Operators required, but, application of FD requires general BC (-1,-1) or (-1,10)
 // in any case tmp2 points to the modified data "X"
 if ( X->Ext.WC->InterpoletFlag && (X->Ext.WC->N<=4) ) {
   //op=FD24 ;
   op = (X->Ext.WC->N == 2) ? FD22 : FD24II ; //FD24II
   assert(tmp2!=tmp) ;
   assert(tmp2!=R) ;
   assert(tmp2!=X) ;

   for (i=0;i<DIM; i++) { 
     GetGeneralBoundaryConditions(X->Ext.BC[i],BG[i]) ;
     if ((X->Ext.BC[i][0]!=BG[i][0]) || (X->Ext.BC[i][1]!=BG[i][1])) change=true ;
   }
   if (change) {
     assert(tmp2) ;
     for (i=0;i<DIM; i++)
       tmp2->ApplyOp(BG[i], (i==0) ? tmp3 : tmp2 , i , PROJECTION) ;
   } else tmp2=tmp3 ;
 } else tmp2=tmp3 ;

 R  ->SetBoundaryConditions(tmp2->Ext.BC) ;
 tmp->SetBoundaryConditions(tmp2->Ext.BC) ;
 // all parameters are just zero:
 //if (last<=0) {
 //  R->Set(0.0) ;
 //  R->SetBoundaryConditions(X) ;
 //  return ;
 //}

 tmp->Add(0.0, tmp2 , par[0],tmp2) ;

 for (i=1;i<=DIM;i++)
   if (par[i]!=0.0) {
     R->ApplyOp(R->Ext.BC[i-1] , tmp2, i-1, op ) ;

     if (i==last)  R->Add  (1.0, tmp , par[i],R) ;
       else      tmp->PPlus(           par[i],R) ;
    }

 if (change)
   for (i=0; i<DIM; i++) R->ApplyOp(X->Ext.BC[i], R , i , PROJECTION) ;

 // transform to Lifting
 for (i=0; i<DIM; i++) 
   if (lift[i]) R->ApplyOp(R->Ext.BC[i], R , i , INTERPOLET2LIFTING) ;

}

//
// Preconditioner::AdaptiveData<DIM>
//
template<int DIM>
void HelmholtzDiagonalPreconditioner(AdaptiveData<DIM> *X , AdaptiveData<DIM> *R, Operator<AdaptiveData<DIM>,DIM> *O , bool inverse) {
 AdaptiveDataArray<DIM>::VectorArray::iterator it(&X->ba->d),dend=X->ba->d.end() ;
 AdaptiveDataArray<DIM>::DataVector *vx,*vr ;
 size_t i,n ;
 double diag ;

 X->CheckAdaptiveGrid(R) ;
 assert(X->SameBoundaryConditions(R)) ;

 // diagonal scaling
 // iterate over all level
 for (it=X->ba->d.begin(); it<=dend; ++it) {
   vx=X->ba->d[it] ; vr=R->ba->d[it] ; n=(*vx).size() ;

   // calc diagonal entry
   diag=O->DiagonalEntry( it.i , X->G->WC) ;

   if (inverse) diag=1/diag ;

// diag=1.0 ;
   for (i=0; i<n; i++) (*vr)[i] = (*vx)[i]/diag ;
 }

 /*
 bool neumann=true ;
 for (i=0; i<DIM; i++) if ((X->Ext.BC[i][0]!=1)||(X->Ext.BC[i][1]!=1)) neumann=false ;
 if (neumann) {
   assert(0) ; 
   int    e[DIM],ld   ;
   double mx=0.0,z ;
   unsigned long j[DIM] ;

   AdaptiveGrid<DIM>::ProjectionGrid           *pg=X->G->BasisPG[0] ;
   AdaptiveGrid<DIM>::ProjectionGrid::iterator  pgit ;

   AdaptiveLevels                     *als  ;
   AdaptiveLevelsIndex                *ali  ;
   AdaptiveLevels::AdaptiveL          *al   ;
   AdaptiveLevels::AdaptiveL::iterator alit ; 
   psg=(*pg)[&X->G->Level0[1]] ;

   for (i=1; i<DIM; i++) {
     e[i]= 1<<X->G->Level0[i] ;
    }

   for (psgit=(*psg).begin(); psgit != (*psg).end(); psgit++) {

      als=(*psgit).second ;
      for (i=1; i<DIM; i++) j[i]=(*psgit).first.i[i-1] ;

      ld=X->G->Level0[0] ;
      e[0]=1<<ld ;
      al=(*als)[ld] ;

      for (alit=(*al).begin(); alit != (*al).end(); alit++) {
        j[0]=(*alit).first ;
        z=fabs( (*R->ba->d[X->G->Level0])[(*alit).second]) ;
        for (i=0;i<DIM;i++) if ((j[i]==1)||(j[i]==(unsigned long)(e[i]-1))) mx = (mx>z) ? mx : z ;
       }
     } // next psgit
  }
 */
}


//
// Preconditioner::UniformData<DIM>
// 
template<int DIM>
void HelmholtzDiagonalPreconditioner(UniformData<DIM> *X , UniformData<DIM> *R,  Operator<UniformData<DIM>,DIM> *O , bool inverse)
{
 int Level0[DIM],J[DIM],i,a[DIM],e[DIM] ;
 double diag ;

 X->ba->SameSize(R->ba) ;

 X->ba->CalcLevel(J) ;
 for (i=0;i<DIM;i++) Level0[i]=X->Ext.WC->Level0 ;

 Matrix<int,DIM>::iterator il(Level0 , J) ;

 for (il.Set(Level0); il<=J; ++il) {
   for (i=0;i<DIM;i++) {
     if (il.i[i]==Level0[i]) a[i]=0 ;
                        else a[i]=(1<<(il.i[i]-1))+1 ;
     e[i]=1<<il.i[i] ;
   }

   Matrix<double,DIM>::iterator it(a,e) ;

   diag=O->DiagonalEntry(il.i , X->Ext.WC) ;
   if (inverse) diag=1/diag ;

   for (it.Set(a); it<=e; ++it) (*R->ba)[it] = (*X->ba)[it]/diag ;
 }
}

//
// Preconditioner::LevelAdaptiveData<DIM>
// 
template<int DIM>
void HelmholtzDiagonalPreconditioner(LevelAdaptiveData<DIM> *X , LevelAdaptiveData<DIM> *R,  Operator<LevelAdaptiveData<DIM>,DIM> *O , bool inverse) { 
  Matrix< Matrix<double,DIM>* , DIM>::iterator l(X->Level0,X->J) ;
  for (l.Set(X->Level0); l<=(*X->a).end(); ++l) 
    if ((*X->A->a)[l.i]) {
      double diag=O->DiagonalEntry(l.i , X->Ext.WC) ;
      if (inverse) diag=1/diag ;

      Matrix<double,DIM>::iterator t((*X->a)[l]),e=(*(*X->a)[l]).end() ;
      Matrix<double,DIM>* al=(*X->a)[l] ;
      Matrix<double,DIM>* bl=(*R->a)[l] ;

      for (t=(*(*X->a)[l]).begin(); t<=e; ++t) (*bl)[t] = (*al)[t]/diag ; 
    }
} 

//
//
template<class I , int DIM>
void InverseHelmholtzDiagonalPreconditioner(I *X , I *Result , Operator<AdaptiveData<DIM>,DIM> *O) {
  HelmholtzDiagonalPreconditioner<DIM>(X,Result,O,true) ; 
}

template<class I , int DIM>
void InverseHelmholtzDiagonalPreconditioner(I *X , I *Result , Operator<UniformData<DIM>,DIM> *O) {
  HelmholtzDiagonalPreconditioner<DIM>(X,Result,O,true) ; 
}

template<class I , int DIM>
void InverseHelmholtzDiagonalPreconditioner(I *X , I *Result , Operator<LevelAdaptiveData<DIM>,DIM> *O) {
  HelmholtzDiagonalPreconditioner<DIM>(X,Result,O,true) ; 
}


//
//  
//
template<class I , int DIM>
HelmholtzOperator<I,DIM>::HelmholtzOperator() { 
  Tmp=Tmp2=Null=NULL ; 
  use_lifting_preconditioner=false ;
}

template<class I , int DIM>
void HelmholtzOperator<I,DIM>::NotUsing(I *X) { assert ((Tmp!=X)&&(Tmp2!=X)/*&&(Null!=X)*/) ;}


template<class I , int DIM>
void HelmholtzOperator<I,DIM>::Apply(I *X , I *R) {
  NotUsing(X) ;
  NotUsing(R) ;

  assert(Tmp) ;

  Tmp->Ext.dimext=R->Ext.dimext=X->Ext.dimext ; 
  if (Tmp2) Tmp2->Ext.dimext=X->Ext.dimext ; 

  Helmholtz<I,DIM>(X,par,R,Tmp,Tmp2 ) ;

  if (X->Ext.dimext >0 ) {
    assert(Null) ;
    assert(X->Ext.dimext==1)  ;
    R->PPlus(X->Ext.ext[0] , Null) ;
    R->Ext.ext[0]=Null->InnerProd(X) ;
    Tmp->Ext.dimext=Tmp2->Ext.dimext=0 ;
  } else {
    //  R->CopyFaces(X,BC) ;
  }

}


template<class I , int DIM>
double HelmholtzOperator<I,DIM>::DiagonalEntry(int *l , Wavelets *W) {
  int i ;
  double d=par[0] ;
  if (W->InterpoletFlag) for (i=0;i<DIM;i++) d -= par[i+1]*pow(4.0,l[i]) ; //exp(log(1<< l[i])*2.0) ;
  else                   for (i=0;i<DIM;i++) d += par[i+1]*pow(4.0,l[i]) ;
  return d;
}


template<class I , int DIM>
void HelmholtzOperator<I,DIM>::Preconditioner(I *X , I *R) {
 HelmholtzDiagonalPreconditioner<DIM>(X,R,this,false) ;
 R->Ext.CopyData(&X->Ext) ;
 
  /*
 int  i ;
 I *Y=X ;
 if (use_lifting_preconditioner) {
   for (i=0; i<DIM; i++) R->ApplyOp( (i==0) ? X : R , i , INTERPOLET2LIFTING) ;
   Y=R ;
  }

 HelmholtzDiagonalPreconditioner<DIM>(Y,R,this,false) ;

 if (use_lifting_preconditioner)
   for (i=0; i<DIM; i++) R->ApplyOp( R , i , LIFTING2INTERPOLET) ;
  */
}

template<class I , int DIM>
void HelmholtzOperator<I,DIM>::InversePreconditioner(I *X , I *R) {
 InverseHelmholtzDiagonalPreconditioner<I,DIM>(X,R,this) ;
 R->Ext.CopyData(&X->Ext) ;

  /*
 int  i ;
 I *Y=X ;
 if (use_lifting_preconditioner) {
   for (i=0; i<DIM; i++) R->ApplyOp( (i==0) ? X : R , i , INTERPOLET2LIFTING) ;
   Y=R ;
  }

 InverseHelmholtzDiagonalPreconditioner<I,DIM>(Y,R,this) ;

 if (use_lifting_preconditioner) 
   for (i=0; i<DIM; i++) R->ApplyOp( R , i , LIFTING2INTERPOLET) ;
  */
}

//
// IntegroDifferentialOperator 
//
template<class I ,int DIM>
IntegroDifferentialOperator<I,DIM>::IntegroDifferentialOperator() {
 par[0]=par[1]=par[2]=0 ;
 D=NULL ;
 K=NULL ;
 Tmp=NULL ;
}

template<class I ,int DIM>
void IntegroDifferentialOperator<I,DIM>::Apply(I *X , I *R) {
 assert(K) ;

 if (par[2]==0) R->Set(0.0); 
          else  R->Multiply(K , X) ;

 R->Add(par[0],X , par[2],R) ;
 
 if (D) {
   NotUsing(X) ;
   NotUsing(R) ;
   D->NotUsing(X) ;
   D->NotUsing(R) ;
   D->NotUsing(Tmp) ;

   D->Apply(X,Tmp) ;
   R->Add(1.0,R , par[1],Tmp) ;
 } 
}

template<class I , int DIM>
void IntegroDifferentialOperator<I,DIM>::NotUsing(I *X) { assert (Tmp!=X) ;}

template<class I , int DIM>
double IntegroDifferentialOperator<I,DIM>::DiagonalEntry(int *l , Wavelets *W) {
  int i ;
  double d=par[0] ;
  if (D) { 
    if (W->InterpoletFlag) for (i=0;i<DIM;i++) d -= par[1]*pow(4.0,l[i]) ; //exp(log(1<< l[i])*2.0) ;
                    else   for (i=0;i<DIM;i++) d += par[1]*pow(4.0,l[i]) ;
   }

  return d;
}

template<class I , int DIM>
void IntegroDifferentialOperator<I,DIM>::Preconditioner(I *X , I *R) {
  HelmholtzDiagonalPreconditioner<DIM>(X,R,this,false) ;
}

template<class I , int DIM>
void IntegroDifferentialOperator<I,DIM>::InversePreconditioner(I *X , I *R) {
  InverseHelmholtzDiagonalPreconditioner<I,DIM>(X,R,this) ;
}


//
// TensorIntegroDifferentialOperator 
//
template<class I ,int DIM>
TensorIntegroDifferentialOperator<I,DIM>::TensorIntegroDifferentialOperator() {
 par[0]=par[1]=par[2]=0 ;
 D=NULL ;
 KX=KY=NULL ;
 Tmp=NULL ;
}

template<class I ,int DIM>
void TensorIntegroDifferentialOperator<I,DIM>::Apply(I *X , I *R) {
 assert(KX) ;
 assert(KY) ;

 double q=KY->InnerProd(X) ;
 
 R->Add(par[0],X , par[2]*q,KX) ;
 
 if (D) {
   NotUsing(X) ;
   NotUsing(R) ;
   D->NotUsing(X) ;
   D->NotUsing(R) ;
   D->NotUsing(Tmp) ;

   D->Apply(X,Tmp) ;
   R->Add(1.0,R , par[1],Tmp) ;
 }
}

template<class I , int DIM>
void TensorIntegroDifferentialOperator<I,DIM>::NotUsing(I *X) { assert (Tmp!=X) ;}

template<class I , int DIM>
double TensorIntegroDifferentialOperator<I,DIM>::DiagonalEntry(int *l , Wavelets *W) {
  int i ;
  double d=par[0] ;
  if (D) { 
    if (W->InterpoletFlag) for (i=0;i<DIM;i++) d -= par[1]*pow(4.0,l[i]) ; //exp(log(1<< l[i])*2.0) ;
                     else  for (i=0;i<DIM;i++) d += par[1]*pow(4.0,l[i]) ;
  }  
  return d;
}

template<class I , int DIM>
void TensorIntegroDifferentialOperator<I,DIM>::Preconditioner(I *X , I *R) {
  HelmholtzDiagonalPreconditioner<DIM>(X,R,this,false) ;
}

template<class I , int DIM>
void TensorIntegroDifferentialOperator<I,DIM>::InversePreconditioner(I *X , I *R) {
  InverseHelmholtzDiagonalPreconditioner<I,DIM>(X,R,this) ;
}


//
// PressureOperatorI
//
template<class I,int DIM>
PressureOperatorI<I,DIM>::PressureOperatorI() {
  Null=Tmp[0]=Tmp[1]=NULL ;
  do_patch=use_lifting_preconditioner=false ;
  // for (int i=0; i<DIM; i++) XL[i]=1.0 ;
}

template<class I,int DIM>
void PressureOperatorI<I,DIM>::NotUsing(I *X) { assert ((Tmp[0]!=X)&&(Tmp[1]!=X)&&(Null!=X)) ;}

template<class I,int DIM>
double PressureOperatorI<I,DIM>::DiagonalEntry(int *l , Wavelets *) {
 int i ;
 double d=0.0 ;
 for (i=0;i<DIM;i++) {
   double dx=(G->XE[i] - G->XA[i]) ;
   d +=1./(dx*dx)*pow(4.0 , l[i] ) ; //exp(log(1<< l[i])*2.0) ;
 } 
 return d;
}

template<class I,int DIM>
void PressureOperatorI<I,DIM>::Preconditioner(I *X , I *R) {
 int i ;
 I *Y=X ;

 if (use_lifting_preconditioner) { for (i=0; i<DIM; i++) R->ApplyOp( (i==0) ? X : R , i , INTERPOLET2LIFTING) ; Y=R ; }

 HelmholtzDiagonalPreconditioner(Y,R,this,false) ;

 if (use_lifting_preconditioner) for (i=0; i<DIM; i++) R->ApplyOp( R , i , LIFTING2INTERPOLET) ;
 
 R->Ext.CopyData(&X->Ext) ;
}

template<class I,int DIM>
void PressureOperatorI<I,DIM>::InversePreconditioner(I *X , I *R) {
 int i ;
 I *Y=X ;

 if (use_lifting_preconditioner) { for (i=0; i<DIM; i++) R->ApplyOp( (i==0) ? X : R , i , INTERPOLET2LIFTING) ; Y=R ; }

 InverseHelmholtzDiagonalPreconditioner(Y,R,this) ;

 if (use_lifting_preconditioner) for (i=0; i<DIM; i++) R->ApplyOp( R , i , LIFTING2INTERPOLET) ;

 R->Ext.CopyData(&X->Ext) ;
}

template<class I,int DIM>
void PressureOperatorI<I,DIM>::Apply(I *X , I *R) {
 NotUsing(X) ;
 NotUsing(R) ;

 int i ;
 int BG[DIM][2] ;
 int op = ( X->G->WC->InterpoletFlag && (X->G->WC->N<=4) ) ? FD24 : STIFFNESS ;

 // treat BC
 R->SetBoundaryConditions(X->Ext.BC) ;
 for (i=0; i<DIM; i++) { 
   GetGeneralBoundaryConditions(X->Ext.BC[i] , BG[i]) ;
   if ((X->Ext.BC[i][0]==0) || (X->Ext.BC[i][1]==0)) assert(!X->extended) ;
 }
 Tmp[0]->Ext.dimext=Tmp[1]->Ext.dimext=X->Ext.dimext ;

 // imbed to general BCs
 for (i=0; i<DIM; i++) Tmp[0]->ApplyOp(BG[i] , (i==0) ? X : Tmp[0] , i , PROJECTION) ;

 // Build up div*grad
 for (i=0; i<DIM; i++) {
   if (i==0) {
          R->ApplyOp(BG[i] , Tmp[0] , i , op) ;
          //R->Mul(1./(XL[0]*XL[0])) ;
    }
   else {
     Tmp[1]->ApplyOp(BG[i] , Tmp[0] , i , op) ;
     R->PPlus(1 /*./(XL[i]*XL[i])*/ , Tmp[1]) ;
    }
 } 
 // project to pressure BCs
 for (i=0; i<DIM; i++) R->ApplyOp(X->BC[i] , R , i , PROJECTION) ;

 if (X->Ext.dimext>0) {
   assert(Null) ;
   assert(X->Ext.dimext==1) ;
   R->Ext.dimext=1 ;
   R->PPlus(X->Ext.ext[0],Null) ;
   R->Ext.ext[0]=Null->InnerProd(X)    ;
   Tmp[0]->Ext.dimext=Tmp[1]->Ext.dimext=0 ;
 }
}


//
// PressureOperatorII
//
template<class I,int DIM>
PressureOperatorII<I,DIM>::PressureOperatorII() {
  Null=Tmp[0]=Tmp[1]=NULL ;
  do_patch=use_lifting_preconditioner=false ;
  //  for (int i=0; i<DIM; i++) XL[i]=1.0 ;
}

template<class I,int DIM>
void PressureOperatorII<I,DIM>::NotUsing(I *X) { assert ((Tmp[0]!=X)&&(Tmp[1]!=X)&&(Null!=X)) ;}

template<class I,int DIM>
double PressureOperatorII<I,DIM>::DiagonalEntry(int *l , Wavelets *) {
 int i ;
 double d=0.0 ;
 for (i=0;i<DIM;i++) d +=pow(4.0 , l[i]) ; //exp(log(1<< l[i])*2.0) ;
 return d;
}

template<class I,int DIM>
void PressureOperatorII<I,DIM>::Preconditioner(I *X , I *R) {
 int i ;
 I *Y=X ;

 if (use_lifting_preconditioner) { for (i=0; i<DIM; i++) R->ApplyOp( (i==0) ? X : R , i , INTERPOLET2LIFTING) ; Y=R ; }

 HelmholtzDiagonalPreconditioner(Y,R,this,false) ;

 if (use_lifting_preconditioner) for (i=0; i<DIM; i++) R->ApplyOp( R , i , LIFTING2INTERPOLET) ;
 
 if (X->extended) R->ext=X->ext ;
}

template<class I,int DIM>
void PressureOperatorII<I,DIM>::InversePreconditioner(I *X , I *R) {
 int i ;
 I *Y=X ;

 if (use_lifting_preconditioner) { for (i=0; i<DIM; i++) R->ApplyOp( (i==0) ? X : R , i , INTERPOLET2LIFTING) ; Y=R ; }

 InverseHelmholtzDiagonalPreconditioner(Y,R,this) ;

 if (use_lifting_preconditioner) for (i=0; i<DIM; i++) R->ApplyOp( R , i , LIFTING2INTERPOLET) ;

 if (X->extended) R->ext=X->ext ;
}

//
// PressureOperator for consistent Laplacian with homg Neumann/periodic boundary conditions
// NO STABILIZATION.
// !!! solvability problems for extended==false and DIVGRADFD4NOSTAB !!!
// 
template<class I,int DIM>
void PressureOperatorII<I,DIM>::Apply(I *X , I *R) {
 NotUsing(X) ;
 NotUsing(R) ;

 int i ;
 int BG[DIM][2] ;
 int op = ( X->WC->InterpoletFlag && (X->WC->N<=2) ) ? DIVGRADFD4NOSTAB : DIVGRADNOSTAB ;

 // treat BC
 R->SetBoundaryConditions(X->BC) ;
 for (i=0; i<DIM; i++) { 
   GetGeneralBoundaryConditions(X->BC[i] , BG[i]) ;
   if ((X->BC[i][0]==0) || (X->BC[i][1]==0)) assert(!X->extended) ;
 }
 Tmp[0]->extended=Tmp[1]->extended=X->extended ;

 // imbed to general BCs
 for (i=0; i<DIM; i++) Tmp[0]->ApplyOp(BG[i] , (i==0) ? X : Tmp[0] , i , PROJECTION) ;

 // Build up div*grad
 for (i=0; i<DIM; i++) {
   if (i==0) {
          R->ApplyOp(BG[i] , Tmp[0] , i , op) ;
          //R->Mul(1./(XL[0]*XL[0])) ;
    }
   else {
     Tmp[1]->ApplyOp(BG[i] , Tmp[0] , i , op) ;
     R->PPlus(1 /*./(XL[i]*XL[i])*/ , Tmp[1]) ;
    }
 } 
 // project to pressure BCs
 for (i=0; i<DIM; i++) R->ApplyOp(X->BC[i] , R , i , PROJECTION) ;

 if (X->extended) {
   assert(Null) ;
   R->extended=true ;
   R->PPlus(X->ext,Null) ;
   R->ext=Null->InnerProd(X) ;
   Tmp[0]->extended=Tmp[1]->extended=false ;
 }
}


//
// PressureOperatorIII
//
template<class I,int DIM>
PressureOperatorIII<I,DIM>::PressureOperatorIII() {
  Null=Tmp[0]=Tmp[1]=NULL ;
  do_restrict=do_patch=use_lifting_preconditioner=false ;
}

template<class I,int DIM>
void PressureOperatorIII<I,DIM>::NotUsing(I *X) { assert ((Tmp[0]!=X)&&(Tmp[1]!=X)&&(Null!=X)) ;}

template<class I,int DIM>
double PressureOperatorIII<I,DIM>::DiagonalEntry(int *l , Wavelets *) {
 int i ;
 double d=0.0 ;
 for (i=0;i<DIM;i++) d += pow(4.0 , l[i]) ;
 return d;
}

template<class I,int DIM>
void PressureOperatorIII<I,DIM>::Preconditioner(I *X , I *R) {
 int i ;
 I *Y=X ;

 if (use_lifting_preconditioner) { for (i=0; i<DIM; i++) R->ApplyOp( (i==0) ? X : R , i , INTERPOLET2LIFTING) ; Y=R ; }

 HelmholtzDiagonalPreconditioner(Y,R,this,false) ;

 if (use_lifting_preconditioner) for (i=0; i<DIM; i++) R->ApplyOp( R , i , LIFTING2INTERPOLET) ;
 
 R->Ext.CopyData(&X->Ext) ;
}

template<class I,int DIM>
void PressureOperatorIII<I,DIM>::InversePreconditioner(I *X , I *R) {
 int i ;
 I *Y=X ;

 if (use_lifting_preconditioner) { for (i=0; i<DIM; i++) R->ApplyOp( (i==0) ? X : R , i , INTERPOLET2LIFTING) ; Y=R ; }

 InverseHelmholtzDiagonalPreconditioner(Y,R,this) ;

 if (use_lifting_preconditioner) for (i=0; i<DIM; i++) R->ApplyOp( R , i , LIFTING2INTERPOLET) ;

 R->Ext.CopyData(&X->Ext) ;
}

//
// PressureOperator for consistent Laplacian with homg Neumann/periodic boundary conditions
template<class I,int DIM>
void PressureOperatorIII<I,DIM>::Apply(I *X , I *R) {
 NotUsing(X) ;
 NotUsing(R) ;

 int i ;
 int BG[DIM][2] ;
 int op = ( X->Ext.WC->InterpoletFlag && (X->Ext.WC->N<=2) ) ? DIVGRADFD4 : DIVGRAD ;

 // treat BC
 R->SetBoundaryConditions(X) ;
 for (i=0; i<DIM; i++) { 
   GetGeneralBoundaryConditions(X->Ext.BC[i] , BG[i]) ;
   if ((X->Ext.BC[i][0]==0) || (X->Ext.BC[i][1]==0)) assert(X->Ext.dimext==0) ;
 }
 Tmp[0]->Ext.dimext=Tmp[1]->Ext.dimext=X->Ext.dimext ;

 // imbed to general BCs
 for (i=0; i<DIM; i++) Tmp[0]->ApplyOp(BG[i] , (i==0) ? X : Tmp[0] , i , PROJECTION) ;

 // build up div*grad
 for (i=0; i<DIM; i++) {
   if (i==0) R->ApplyOp(BG[i] , Tmp[0] , i , op) ;
   else {
     Tmp[1]->ApplyOp(BG[i] , Tmp[0] , i , op) ;
     R->PPlus( Tmp[1] ) ;
    }
 }

 // project to pressure BCs
 for (i=0; i<DIM; i++) R->ApplyOp(X->Ext.BC[i] , R , i , PROJECTION) ;

 if (X->Ext.dimext) {
   assert(Null) ;
   R->Ext.dimext=1 ;
   R->PPlus(X->Ext.ext[0],Null) ;
   R->Ext.ext[0]=Null->InnerProd(X) ;
   Tmp[0]->Ext.dimext=Tmp[1]->Ext.dimext=0 ;
 }


}


//
// PressureOperatorIV
//
template<class I,int DIM>
PressureOperatorIV<I,DIM>::PressureOperatorIV() {
  Null=Tmp[0]=Tmp[1]=NULL ;
  do_restrict=do_patch=use_lifting_preconditioner=false ;
  //  for (int i=0; i<DIM; i++) XL[i]=1.0 ;
}

template<class I,int DIM>
void PressureOperatorIV<I,DIM>::NotUsing(I *X) { assert ((Tmp[0]!=X)&&(Tmp[1]!=X)&&(Null!=X)) ;}

template<class I,int DIM>
double PressureOperatorIV<I,DIM>::DiagonalEntry(int *l , Wavelets *) {
 int i ;
 double d=0.0 ;
 for (i=0;i<DIM;i++) d +=/*1./(XL[i]*XL[i])**/pow(4.0 , l[i]) ; //exp(log(1<< l[i])*2.0) ;
 return d;
}

template<class I,int DIM>
void PressureOperatorIV<I,DIM>::Preconditioner(I *X , I *R) {
 int i ;
 I *Y=X ;

 HelmholtzDiagonalPreconditioner(Y,R,this,false) ;
 
 if (X->extended) R->ext=X->ext ;
}

template<class I,int DIM>
void PressureOperatorIV<I,DIM>::InversePreconditioner(I *X , I *R) {
 int i ;
 I *Y=X ;

 InverseHelmholtzDiagonalPreconditioner(Y,R,this) ;

 if (X->extended) R->ext=X->ext ;
}

//
// PressureOperator for consistent Laplacian with homg Neumann/periodic boundary conditions
template<class I,int DIM>
void PressureOperatorIV<I,DIM>::Apply(I *X , I *R) {
 NotUsing(X) ;
 NotUsing(R) ;

 int i,j ;
 int BC[DIM][2] ;

 // treat BC
 R->SetBoundaryConditions(X->BC) ;
 for (i=0; i<DIM; i++) {
   assert(X->BC[i][0]==-1) ;
   assert((X->BC[i][1]==-1) || (X->BC[i][1]==PERIODIC)) ;

   if (X->BC[i][1]==PERIODIC) { BC[i][0]=-1, BC[i][1]=PERIODIC ;}
                              else  { BC[i][0]=BC[i][1]=0 ;}
 }

 Tmp[0]->extended=Tmp[1]->extended=X->extended ;

 // Build up div*grad
 R->Set(0.0) ; 
 R->ext=0    ;
 for (i=0; i<DIM; i++) {
   Tmp[0]->ApplyOp(X->BC[i] , X , i , FD14) ;
   for (j=0; j<DIM; j++) {
     Tmp[0]->ApplyOp(   BC[j] , Tmp[0] , j , PROJECTION) ; 
     Tmp[0]->ApplyOp(X->BC[j] , Tmp[0] , j , PROJECTION) ;
    }

   /*
   int A[2]={0,0},J[2]={6,6},E[2]={64,64} ;
   Matrix<double,2> MM(A,E) ;
   Tmp[0]->ToUniform(&MM,J) ;
   if ((i==1)&& (debugRefine!=0)) { MM.WriteMATLAB("R1") ;exit(-1) ;}
   */

   Tmp[0]->ApplyOp(X->BC[i] , Tmp[0] , i , FD14) ;

   R->PPlus(1/*./(XL[i]*XL[i])*/ , Tmp[0]) ;
 } 

 if (X->extended) {
   assert(Null) ;
   R->extended=true ;
   R->PPlus(X->ext,Null) ;
   R->ext=Null->InnerProd(X) ;
   Tmp[0]->extended=Tmp[1]->extended=false ;
 }

}

//
// PressureOperatorV
//
template<class I,int DIM>
PressureOperatorV<I,DIM>::PressureOperatorV() {
  Null=Tmp=NULL ;
  do_restrict=do_patch=use_lifting_preconditioner=false ;
  // for (int i=0; i<DIM; i++) XL[i]=1.0 ;
}

template<class I,int DIM>
void PressureOperatorV<I,DIM>::NotUsing(I *X) { assert ((Tmp!=X)&&(Null!=X)) ;}

template<class I,int DIM>
double PressureOperatorV<I,DIM>::DiagonalEntry(int *l , Wavelets *) {
 return 1.0;
}

template<class I,int DIM>
void PressureOperatorV<I,DIM>::Preconditioner(I *X , I *R) {
 R->Copy(X) ;
 R->ext=X->ext ;
}

template<class I,int DIM>
void PressureOperatorV<I,DIM>::InversePreconditioner(I *X , I *R) {
 R->Copy(X)    ;
 R->ext=X->ext ;
}

/*
// PressureOperator for consistent Laplacian with homg Neumann/periodic boundary conditions
template<class I,int DIM>
void PressureOperatorV<I,DIM>::Apply(I *X , I *R) {
 NotUsing(X) ;
 NotUsing(R) ;

 int i,j ;
 int BC[DIM][2] ;

 for (i=0; i<DIM; i++) BC[i][0]=BC[i][1]=-1 ;

 //Tmp[0]->extended=Tmp[1]->extended=X->extended ;

 // Build up div*grad
 R->Set(0.0) ;
 for (i=0; i<DIM; i++) {
   int op[DIM]={NOOPERATION , NOOPERATION},re[DIM] ;
   op[i]=FD14 ;
   Tmp->ba->TensorMatrixOperation(BC,X->ba,BC,X->WC,applyOp,op,re) ;
   Tmp->extended=X->extended ;

   Matrix<double,DIM>::iterator s(Tmp->ba) ;
   for (s=(*Tmp->ba).begin(); s<=(*Tmp->ba).end(); ++s) 
     for (j=0; j<DIM; j++) if ((s.i[j]==0) || (s.i[j]==Tmp->ba->e[j])) (*Tmp->ba)[s]=0.0 ;

   op[i]=FD14 ;
   Tmp->ba->TensorMatrixOperation(BC,Tmp->ba,BC,X->WC,applyOp,op,re) ;

   R->PPlus(1 , Tmp) ;
 } 

 if (X->extended) {
   assert(Null) ;
   R->extended=true ;
   R->PPlus(X->ext,Null) ;
   R->ext=Null->InnerProd(X) ;
   Tmp->extended=false ;
 }

}

template<class I,int DIM>
void PressureOperatorV<I,DIM>::Applyi(I *X , I *R,int i) {
 NotUsing(X) ;
 NotUsing(R) ;

 int j ;
 int BC[DIM][2] ;

 for (j=0; j<DIM; j++) BC[j][0]=BC[j][1]=-1 ;

 //Tmp[0]->extended=Tmp[1]->extended=X->extended ;

 // Build up div*grad
 R->Set(0.0) ;

   int op[DIM]={NOOPERATION , NOOPERATION},re[DIM] ;
   op[i]=FD12 ;
   Tmp->ba->TensorMatrixOperation(BC,X->ba,BC,X->WC,applyOp,op,re) ;
   Tmp->extended=X->extended ;

   Matrix<double,DIM>::iterator s(Tmp->ba) ;
   for (s=(*Tmp->ba).begin(); s<=(*Tmp->ba).end(); ++s) 
     for (j=0; j<DIM; j++) if ((s.i[j]==0) || (s.i[j]==Tmp->ba->e[j])) (*Tmp->ba)[s]=0.0 ;

   Tmp->ba->TensorMatrixOperation(BC,Tmp->ba,BC,X->WC,applyOp,op,re) ;

   R->PPlus(1 , Tmp) ;


 if (X->extended) {
   assert(Null) ;
   R->extended=true ;
   R->PPlus(X->ext,Null) ;
   R->ext=Null->InnerProd(X) ;
   Tmp->extended=false ;
 }
}
*/

//
// ConvectiveOperator
//

template<class I , int DIM>
ConvectionOperator<I,DIM>::ConvectionOperator() {
  Tmp=NULL ;
  for (int i=0;i<DIM;i++) { U[i]=NULL /*, XL[i]=1.0*/; }
}

template<class I , int DIM>
void ConvectionOperator<I,DIM>::NotUsing(I *X) { assert (Tmp!=X); assert (Tmp1!=X); }

template<class I , int DIM>
double ConvectionOperator<I,DIM>::DiagonalEntry(int * , Wavelets *) { assert(0) ; return 0 ;}

template<class I , int DIM>
void ConvectionOperator<I,DIM>::Preconditioner(I * , I * ) { assert(0) ; }
template<class I , int DIM>
void ConvectionOperator<I,DIM>::InversePreconditioner(I * , I * ) { assert(0) ; }

template<class I , int DIM>
void ConvectionOperator<I,DIM>::Apply(I *X , I *R) {
  int i,j ;

  NotUsing(X) ;
  NotUsing(R) ;

  for (i=0; i<DIM; i++) {
    for (j=0; j<DIM; j++)
      Tmp->ApplyOp( U[0]->Ext.BC[j] , (j==0) ? X : Tmp , j , (j==i) ? FIRSTDERIVATIVE : PROJECTION) ;

    for (j=0; j<DIM; j++)
      Tmp->ApplyOp( Tmp , j , INVERSETRANSFORM) ;

    if (i==0)
      R->Mul(Tmp , U[i]) ;
    else {
      Tmp->Mul  (Tmp , U[i]) ;
        R->PPlus(Tmp) ;
      }
  }

  for (j=0; j<DIM; j++) R->ApplyOp( R ,j , TRANSFORM) ;
  for (j=0; j<DIM; j++) R->ApplyOp( X->Ext.BC[j] , R , j , PROJECTION) ;

} 

//
template<class I , int DIM>
void ConvectionOperator<I,DIM>::ApplyFD(I *X , I *R) {
  int i,j,BG[DIM][2] ;

  NotUsing(X) ;
  NotUsing(R) ;

  for (j=0; j<DIM; j++) GetGeneralBoundaryConditions(X->Ext.BC[j],BG[j]) ;

  for (j=0; j<DIM; j++) Tmp1->ApplyOp( BG[j] , (j==0) ? X : Tmp1 , j , PROJECTION) ;

  for (i=0; i<DIM; i++) {
    
    Tmp->ApplyOp( BG[i] , Tmp1 , i , FD14) ;

    for (j=0; j<DIM; j++) Tmp->ApplyOp( Tmp , j , INVERSETRANSFORM) ;

    if (i==0)
      R->Mul(Tmp , U[i]) ;
    else {
      Tmp->Mul  (Tmp , U[i]) ;
        R->PPlus(Tmp) ;
      }
  }

  for (j=0; j<DIM; j++) R->ApplyOp( R ,j , TRANSFORM) ;
  for (j=0; j<DIM; j++) R->ApplyOp( X->Ext.BC[j] , R , j , PROJECTION) ;

}


template<class I , int DIM>
void ConvectionOperator<I,DIM>::ApplyConservativeFast(int dir , I *R) {
  int i,j ;

  for (i=0; i<DIM; i++) 
    for (j=0; j<DIM; j++) 
      assert(!U[j]->Ext.islifting[i]) ;

  NotUsing(X) ;
  NotUsing(R) ;

  for (i=0; i<DIM; i++) {
    Tmp->Mul(U[dir] , U[i]) ;
    for (j=0; j<DIM; j++) Tmp->ApplyOp( Tmp , j , TRANSFORM) ;

    if (i==0)    R->ApplyOp( U[0]->Ext.BC[i] , Tmp , i , FIRSTDERIVATIVE) ;
       else {  Tmp->ApplyOp( U[0]->Ext.BC[i] , Tmp , i , FIRSTDERIVATIVE) ;
                 R->Add(Tmp) ; }
  }

}


template<class I , int DIM>
void ConvectionOperator<I,DIM>::ApplyConservative(I *X , I *R) {
  int i,j,BG[DIM][2] ;

  for (i=0; i<DIM; i++) assert(!X->Ext.islifting[i]) ;

  NotUsing(X) ;
  NotUsing(R) ;

  assert(Tmp)  ;
  assert(Tmp1) ;

  for (j=0; j<DIM; j++) GetGeneralBoundaryConditions(X->Ext.BC[j],BG[j]) ;
  for (j=0; j<DIM; j++) Tmp1->ApplyOp(BG[j] , (j==0) ? X : Tmp1 , j , PROJECTION) ;

  for (j=0; j<DIM; j++) Tmp1->ApplyOp(BG[j], Tmp1, j, INVERSETRANSFORM) ;

  for (i=0; i<DIM; i++) {
    Tmp->Mul(U[i] , Tmp1 ) ;
    for (j=0; j<DIM; j++) Tmp->ApplyOp( Tmp , j , TRANSFORM) ;

    if (i==0)    R->ApplyOp(BG[i] , Tmp , i , FIRSTDERIVATIVE) ;
       else {  Tmp->ApplyOp(BG[i] , Tmp , i , FIRSTDERIVATIVE) ;
                 R->PPlus(Tmp) ; }
  }

  for (j=0; j<DIM; j++) R->ApplyOp(X->Ext.BC[j] , R , j , PROJECTION) ;
}


//
// ApplyConservativeFD was instable in ModonsAdap/Merging problem
// on the other hand, it worked well for the simple transport problem ....
//
template<class I , int DIM>
void ConvectionOperator<I,DIM>::ApplyConservativeFD(I *X , I *R) {
  int i,j,BG[DIM][2] ;
  bool change=false  ;

  for (i=0; i<DIM; i++) assert(!X->Ext.islifting[i]) ;

  NotUsing(X) ;
  NotUsing(R) ;

  assert(Tmp) ;
  assert(Tmp1) ;

  for (j=0; j<DIM; j++) {
    GetGeneralBoundaryConditions(X->Ext.BC[j],BG[j]) ;
    Tmp1->ApplyOp(BG[j] , (j==0) ? X : Tmp1 , j , PROJECTION) ;
  }

  for (j=0; j<DIM; j++) Tmp1->ApplyOp(BG[j] , Tmp1 , j , INVERSETRANSFORM) ;

  for (i=0; i<DIM; i++) {
    Tmp->Mul(U[i] , Tmp1 ) ;
    for (j=0; j<DIM; j++) Tmp->ApplyOp(BG[j] ,  Tmp , j , TRANSFORM) ;

    if (i==0)    R->ApplyOp( BG[i] , Tmp , i , FD14) ;
       else {  Tmp->ApplyOp( BG[i] , Tmp , i , FD14) ;
                 R->PPlus(Tmp) ; }
  }

  for (j=0; j<DIM; j++) R->ApplyOp(X->Ext.BC[j] , R , j , PROJECTION) ;
}



//
//
//
template<class I , int DIM>
void ConvectionOperator<I,DIM>::ApplyConservativeWENO(I *X , I *R) {
  int i,j,BG[DIM][2] ;

  for (i=0; i<DIM; i++) assert(!X->Ext.islifting[i]) ;

  NotUsing(X) ;
  NotUsing(R) ;

  assert(Tmp) ;
  assert(Tmp1) ;

  for (j=0; j<DIM; j++) {
    GetGeneralBoundaryConditions(X->Ext.BC[j],BG[j]) ;
    Tmp1->ApplyOp(BG[j] , (j==0) ? X : Tmp1 , j , PROJECTION) ;
  }

  for (j=0; j<DIM; j++) Tmp1->ApplyOp(BG[j] , Tmp1 , j , INVERSETRANSFORM) ;

  for (i=0; i<DIM; i++) {
    Tmp->Mul(U[i] , Tmp1 ) ;
    for (j=0; j<DIM; j++) Tmp->ApplyOp(BG[j] ,  Tmp , j , TRANSFORM) ;

    if (i==0)    R->ApplyWENO(BG[i], Tmp, U[i], i , WENO5) ;
       else {  Tmp->ApplyWENO(BG[i], Tmp, U[i], i , WENO5) ;
                 R->PPlus(Tmp) ; }
  }

  for (j=0; j<DIM; j++) R->ApplyOp(X->Ext.BC[j] , R , j , PROJECTION) ;
 
}



//
//
//
template<class I , int DIM>
void ConvectionOperator<I,DIM>::ApplyConservativeWENOUniformData(UniformData<DIM> *X , UniformData<DIM> *R, bool trafo) {
  int i,j,dir,BG[DIM][2],J[DIM] ;

  NotUsing(X) ;
  NotUsing(R) ;

  assert(Tmp)  ;
  assert(Tmp1) ;

  X->ba->CalcLevel(J) ;

  for (j=0; j<DIM; j++) {
    assert(X->Ext.ismultiscale[j]==  false) ;
    assert(X->Ext.BC[j][1]==PERIODIC) ;
    GetGeneralBoundaryConditions(X->Ext.BC[j],BG[j]) ;
  }
  Tmp1->Copy(X) ;

  // *this=0 ;
  R->CopyExtensions(Tmp1) ;
  R->Set(0.0) ;

  for (dir=0; dir<DIM; dir++) {
    Tmp->Mul(U[dir] , Tmp1 ) ;

    if (trafo) {
      for (i=0; i<DIM; i++) 
        if (i!=dir) Tmp->ApplyOp(BG[i], Tmp, i, TRANSFORM) ;
    }

    Matrix<double,DIM>::iterator s, Xend=(*X->ba).end() ;
    int b=X->ba->b[dir] , e=X->ba->e[dir] ;

    IndexSet IS(2) ; IS.a[0]=b; IS.e[0]=e; IS.nr=1 ;

    for (s=(*X->ba).begin(dir) ; s<=Xend ; ++s) {
      double *Fs=_DATAMATRIX_BUFFER_s, *Ft=_DATAMATRIX_BUFFER_t, *Fr=_DATAMATRIX_BUFFER_q ;

      int *i=&s.i[dir] ;
      for ((*i)=b; (*i)<=e; (*i)++) { 
        Fs[(*i)-b]=(*Tmp->ba)   [s] ; // F
        Fr[(*i)-b]=(*U[dir]->ba)[s] ; // roespeed
      }

      IS.VecApplyWENO(Ft,Fs, Fr, BG[dir], J[dir], 3 , X->Ext.XE[dir]-X->Ext.XA[dir]) ;

      for ((*i)=b; (*i)<=e; (*i)++) (*Tmp->ba)[s] = Ft[(*i)-b] ;
    }
    
    if (trafo) {
      for (i=0; i<DIM; i++) 
        if (i!=dir) Tmp->ApplyOp(BG[i], Tmp, i, INVERSETRANSFORM) ;
    }
   
    R->PPlus(Tmp) ;

  }

  //for (j=0; j<DIM; j++) R->ApplyOp(BG[j] ,  R , j , TRANSFORM) ;
}


template<class I,int DIM>
SpaceTimeOperator<I,DIM>::SpaceTimeOperator() {
 Tmp=Tmp2=NULL ;
 use_lifting_preconditioner=false ;
}

template<class I,int DIM>
void SpaceTimeOperator<I,DIM>::Apply(I *X , I *Result) {

  int i,BG[DIM][2] ;
  for (i=0; i<DIM; i++) GetGeneralBoundaryConditions(X->Ext.BC[i],BG[i]) ;
  for (i=0; i<DIM; i++) Tmp2->ApplyOp(BG[i], (i==0) ? X : Tmp2, i , PROJECTION) ; 
  Result->ApplyOp(BG[0], Tmp2, 0, FD11) ;
  for (i=1; i<DIM; i++) {
    Tmp   ->ApplyOp(BG[i], Tmp2, i, FD24) ;
    Result->MMinus(1.0 , Tmp) ;
  }
  for (i=0; i<DIM; i++) Result->ApplyOp(X->Ext.BC[i], Result, i , PROJECTION) ;

}

template<class I,int DIM>
void SpaceTimeOperator<I,DIM>::DiagonalScale(I *X, bool inverse) {

 // iterate over all level
 Matrix< vector<double> * , DIM>::iterator it(&X->ba->d),dend=X->ba->d.end() ;
 AdaptiveDataArray<DIM>::DataVector *vx ;
 size_t i,n ;
 for (it=X->ba->d.begin(); it<=dend; ++it) {
   vx=X->ba->d[it] ; n=(*vx).size() ;

   // calc diagonal entry
   double diag=pow(2.0 , it.i[0]) ;
   for (i=1; i<DIM; i++) diag += pow(4.0 , it.i[i]) ;
   if  (!inverse) diag=1.0/diag ;
   for (i=0; i<n; i++) (*vx)[i] *= diag ;
 } 
 
}

template<class I,int DIM>
void SpaceTimeOperator<I,DIM>::Preconditioner(I *X , I *R) {
  int i;
  R->Copy(X) ;
  if (use_lifting_preconditioner) for (i=1; i<DIM; i++) R->ApplyOp(R->Ext.BC[i], R, i, INTERPOLET2LIFTING) ;
  DiagonalScale(R,false) ;
  if (use_lifting_preconditioner) for (i=1; i<DIM; i++) R->ApplyOp(R->Ext.BC[i], R, i, LIFTING2INTERPOLET) ;
}

template<class I,int DIM>
void SpaceTimeOperator<I,DIM>::InversePreconditioner(I *X , I *R) {
  int i;
  R->Copy(X) ;
  if (use_lifting_preconditioner) for (i=1; i<DIM; i++) R->ApplyOp(R->Ext.BC[i], R, i, INTERPOLET2LIFTING) ;
  DiagonalScale(R,true) ;
  if (use_lifting_preconditioner) for (i=1; i<DIM; i++) R->ApplyOp(R->Ext.BC[i], R, i, LIFTING2INTERPOLET) ;
}

template<class I,int DIM>
void SpaceTimeOperator<I,DIM>::NotUsing(I *X) { assert ( (Tmp!=X) && (Tmp2!=X)) ;}

template<class I,int DIM>
double SpaceTimeOperator<I,DIM>::DiagonalEntry(int * , Wavelets *) { return 1; }
#endif

---