Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

AdaptiveData.hpp

#ifndef _ADAPTIVE_DATA_
# define _ADAPTIVE_DATA_

#include "Refine.hpp"
#include "Adaptive1D.hpp"
#include "AdaptiveDataArray.hpp"
#include "Extensions.hpp"
#include "UniformData.hpp"
#include "Function.hpp"
#include "SparseMatrixIO.hpp"

template<int DIM>
struct AdaptiveGrid ;

template<int DIM> 
struct AdaptiveData{
 
  typedef int* BCType[DIM] ;

public:
         AdaptiveData()        ;
         AdaptiveData(AdaptiveGrid<DIM> *G) ;
        ~AdaptiveData()        ;

  void   Attach(AdaptiveGrid<DIM> *G) ;
  void   DeAttach()  ;

  size_t Size() ;

  size_t AllocatedEntries() ;

  void   SetRefine() ;

  // Refine
  void   Refine() ;
  void   Resize() ;
  // if (flag>0) (*this)-- ; else (*this)=CounterInitValue ; 
  void   UpdateRefineCounters(AdaptiveData<DIM> flags, int CounterInitValue) ;
  //
  void   Bining(double *eps, int Keps, int *count, AdaptivityCriterion *R) ; 
                                  // counts the number of waveletcoefficients c (scaled dependent on the level according to *R) in the buckets 
                                  // (c < eps[-Keps]) , (eps[-Keps]<= c < eps[-Keps+1]), .... , 
                                  //                    (eps[Keps-1]<= c < eps[Keps])  , ( eps[Keps] <= c)

  //
  //  Access to numerical values
  //

  void   FromUniform(UniformData<DIM>   *A, int *J) ;
  void   FromUniform(Matrix<double,DIM> *A, int *J) ;

  void   ToUniform  (UniformData<DIM>   *A, int *J) ;
  void   ToUniform  (Matrix<double,DIM> *A, int *J, int flag=0) ;

  // Special:
  double *GetP(int *l,int *t) ;
  int    Set(int *l,int *t,const double a) ;
  double Get(int *l,int *t) ;
  void   Set(int *l,const double a) ;
  void   Set(int dir, int k, int s, const double a) ; 
  void   Set(int  l,const double)   ;        // Set all coef of Level *k with k[0]>=l or k[1] >=l to a 
  void   Get(int *l,int *t,int dir, AdaptiveB *B, int leftrightall=0) ;
  void   Set(int *l,int *t,int dir, AdaptiveB *B, int leftrightall=0, bool add=false) ;

  void   ReadSlice (AdaptiveData<DIM+1> *X) ;
  void   WriteSlice(AdaptiveData<DIM+1> *X, bool add=false) ;

  void   SetFunction(Function *F, bool boundaryonly=false, bool notransform=false) ;
  
  void   SetBoundaryValueFunction(AdaptiveData<DIM-1> *F[DIM][2] , int BC[DIM][2] , bool ContinuousDirichletValues=false) ;

  void   SetDirichletBoundaryValues(AdaptiveData<DIM-1> *F[DIM][2] , int BC[DIM][2] , bool ContinuousDirichletValues=false) ;
  void   SetNeumannBoundaryValues  (AdaptiveData<DIM-1> *F[DIM][2] , int BC[DIM][2]) ;
  void   SetNeumannFaces           (AdaptiveData<DIM-1> *F[DIM][2] , int BC[DIM][2]) ;  // internal function, Neumann-values must have homogeneous normal derivatives
  AdaptiveData<DIM-3>* GetDi2Di1Fi0k0(AdaptiveData<DIM-1> *F[DIM][2] , int i2,int k2, int i1, int k1, int i0, int k0) ; // internal
  void   ComputeDerivatives        (AdaptiveData<DIM-1> *F[DIM][2] , int BC[DIM][2], int FD, int codim) ; // internal function
  void   InitCodimIteration        (int *t, int *n, int *f, int codim, Wavelets *W, int *tr,int *a, int *e) ;     // internal
  void   AddBiTriLinearContributionCodim2(int *l,int *t,int *n,int *f,int codim, AdaptiveData<DIM-1> *F[DIM][2],int BC[DIM][2]) ;// internal
  void   AddBiTriLinearContribution      (int *l,int *t,int *n,int *f,int codim, AdaptiveData<DIM-1> *F[DIM][2],int BC[DIM][2]) ;// internal

  void   Clear()  ;

  double Max()    ;
  double Min()    ;
  double MaxAbs() ;
  double Sum()    ;
  double SumAbs() ;
  double SumSqr() ;

  double LpNorm  (AdaptiveData<DIM> *Tmp, double p) ;
  double LmaxNorm(AdaptiveData<DIM> *Tmp) ;
  double L2Norm  (AdaptiveData<DIM> *Tmp) ;
  double L1Norm  (AdaptiveData<DIM> *Tmp) ;
  double Integrate() ;  
  double InnerProd(AdaptiveData<DIM> *X) ;

  void   Set(const double a)   ;
  void   SetDiagonalMatrix(const double p, const double alpha) ; // diag(l,t)=1+alpha*sum_i 2^(l_i p)
  void   Copy(AdaptiveData<DIM> *X) ;
  void   Abs(AdaptiveData<DIM> *X)  ;
  void   Pow(AdaptiveData<DIM> *X, double p) ; // (*this)=pow(X,p) ;
  void   Sqr(AdaptiveData<DIM> *X)  ;

  void   PPlus(const double a) ;
  void   PPlus(AdaptiveData<DIM> *X  ) ;
  void   PPlus(const double a , AdaptiveData<DIM> *X) ;
  void   Add(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y) ;
  void   Add(const double a , AdaptiveData<DIM> *X , const double b , AdaptiveData<DIM> *Y) ;
  void   Add(const double a , AdaptiveData<DIM> *X , const double b , AdaptiveData<DIM> *Y , const double c , AdaptiveData<DIM> *Z) ;
  void   Add(const double a , AdaptiveData<DIM> *X , const double b , AdaptiveData<DIM> *Y , const double c , AdaptiveData<DIM> *Z , 
             const double d , AdaptiveData<DIM> *Q) ;
  void   Add(const double a , AdaptiveData<DIM> *X , const double b , AdaptiveData<DIM> *Y , const double c , AdaptiveData<DIM> *Z ,
             const double d , AdaptiveData<DIM> *Q , const double f , AdaptiveData<DIM> *R) ;

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

  void   Mul(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y) ; 
  void   Mul(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , double f) ;
  //void   Mul(AdaptiveDataArray<DIM> *X , AdaptiveDataArray<DIM> *Y)  ;
  void   Mul(AdaptiveData<DIM> *X , double f) ;
  void   Mul(double f) ;
  void   XplusYmalZ (AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z) ;
  void   XminusYmalZ(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z) ;
 
  void   Div(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y) ; // (*this):= (*X) / (*Y)  element by element
  void   Div(AdaptiveDataArray<DIM> *X , AdaptiveDataArray<DIM> *Y)  ;

  void   ApplyOp  (int *BCT , AdaptiveData<DIM> *X,int dir,unsigned int op) ;
  void   ApplyOp  (AdaptiveData<DIM> *X,int dir,unsigned int op) ;
  // multiply by diagonal matrix D_(l,t),(l,t) = \sum_i 2^(p l_i)
  void   ApplyDiagonal(AdaptiveData<DIM> *X , double p) ;
  // WENO
  void   ApplyWENO(int *BCT , AdaptiveData<DIM> *X, AdaptiveData<DIM> *RoeSpeed, int dir,unsigned int op) ;

  // Projection
  double ZeroMean() ;


  //
  //   IO
  // 

  void   WriteUDF(const char *name, UniformData<DIM> *M, AdaptiveData<DIM> *Tmp=NULL, bool CastToFloat=false) ;
  void   WriteUDF(const char *name,             int  *L, AdaptiveData<DIM> *Tmp=NULL, bool CastToFloat=false) ;
  void   WriteUDF(const char *name,             int   L, AdaptiveData<DIM> *Tmp=NULL, bool CastToFloat=false) ;

  void   WriteSparse(const char *name, AdaptiveData<DIM> **M, int N, bool CastToFloat=false, AdaptivityCriterion *R=NULL) ;
  void   WriteSparse(const char *name, AdaptivityCriterion *R=NULL) ;
  void   ReadSparse (const char *name, AdaptiveData<DIM> **M, int N, int *which) ;
  void   ReadSparse(const char *name) ;

  //
  //  Boundary Conditions
  //

  void   SetBoundaryConditions(AdaptiveData<DIM> *X) ;
  void   SetBoundaryConditions(BCType bc) ;
  void   SetBoundaryConditions(int bc[DIM][2]) ;
  // Set boundary conditions with respect to coordinate direction dir
  void   SetBoundaryConditions(int dir,int *bc) ;
  // Set boundary conditions with respect to coordinate direction dir
  void   SetBoundaryConditions(int dir,int bc0,int bc1) ;
  bool   SameBoundaryConditions(AdaptiveData<DIM> *X) ;
  
  // Compatibility Check for multi-array operations like A+B,A+B+C and so on 
  void   CheckBC(AdaptiveData<DIM> *X) ;
  void   CheckBC(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y) ;
  void   CheckBC(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z) ;
  void   CheckBC(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z , AdaptiveData<DIM> *Q) ;
  void   CheckBC(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z , AdaptiveData<DIM> *Q , AdaptiveData<DIM> *R) ;

  //
  void   Compatible(AdaptiveData<DIM> *X) ;
  void   Compatible(AdaptiveData<DIM> *X, AdaptiveData<DIM> *Y) ;
  void   Compatible(AdaptiveData<DIM> *X, AdaptiveData<DIM> *Y, AdaptiveData<DIM> *Z) ;
  void   Compatible(AdaptiveData<DIM> *X, AdaptiveData<DIM> *Y, AdaptiveData<DIM> *Z, AdaptiveData<DIM> *Q) ;
  void   Compatible(AdaptiveData<DIM> *X, AdaptiveData<DIM> *Y, AdaptiveData<DIM> *Z, AdaptiveData<DIM> *Q, AdaptiveData<DIM> *R) ;

  //
  void   CheckAdaptiveGrid(AdaptiveData<DIM> *X) ;
  void   CheckAdaptiveGrid(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y) ;
  void   CheckAdaptiveGrid(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z) ;
  void   CheckAdaptiveGrid(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z , AdaptiveData<DIM> *Q) ;
  void   CheckAdaptiveGrid(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z , AdaptiveData<DIM> *Q , AdaptiveData<DIM> *R) ;

  void   CopyExtensions(AdaptiveData<DIM> *X) ;

  // numerical data
  AdaptiveDataArray<DIM>  *ba ;
  // adaptive grid
  AdaptiveGrid<DIM>       *G  ;
  // boundary conditions, filter coefficients and augmentions for singular linear systems
  Extensions<DIM>          Ext ;
} ;


//
//
//
template<>
struct AdaptiveData<0> {

   AdaptiveData()                    {G=NULL;}
   AdaptiveData(AdaptiveGrid<0> *Gr) {G=Gr;}
  ~AdaptiveData() {;}
 
  void   ReadSlice (AdaptiveData<1> *X)                 ;
  void   WriteSlice(AdaptiveData<1> *X, bool add=false) ;

  void   Add (const double a , AdaptiveData<0> *X , const double b , AdaptiveData<0> *Y) {data=a*X->data+b*Y->data;}
  void   Copy(AdaptiveData<0> *X) {data=X->data;}

  double*GetP(int * ,int * )                {return &data;}
  int    Set (int * ,int * ,const double a) {data=a; return 0;}
  double Get (int * ,int * )                {return data;}

  AdaptiveGrid<0> *G ;
  // numerical data
  double           data ;
} ;


#include "AdaptiveGrid.hpp"

//                       //
// Data::Implementations //
//                       //

// for specialization
void AdaptiveData<0>::ReadSlice (AdaptiveData<1> *X) { 
  assert(G->SliceParent==X->G); 
  data=X->Get(&G->Slicel,&G->Slicet); 
}

void AdaptiveData<0>::WriteSlice(AdaptiveData<1> *X, bool add) { 
  assert(G->SliceParent==X->G); 
  double *p=X->GetP(&G->Slicel,&G->Slicet); 
  if (add) (*p)+=data; 
  else     (*p) =data;
}

//
// general
//
template<int DIM>
AdaptiveData<DIM>::AdaptiveData() { 
  G=NULL ;
  ba=new AdaptiveDataArray<DIM> ; 
}


template<int DIM>
AdaptiveData<DIM>::AdaptiveData(AdaptiveGrid<DIM> *Gr) {
  G=NULL ;
  ba=new AdaptiveDataArray<DIM> ; 
  Attach(Gr) ;
}


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


template<int DIM>
void AdaptiveData<DIM>::Clear(){ ba->Clear() ;}


template<int DIM>
void AdaptiveData<DIM>::Attach(AdaptiveGrid<DIM> *Gr) {
  assert((G==NULL)&&(Gr!=NULL)) ;
  G=Gr ;
  G->ResizeGroup.Insert(this)  ;
  G->RefineGroup.Remove(this)  ;
  ba->Resize(&G->CounterArray) ;
  Ext.WC=G->WC ;
  Ext.XA=G->XA ;
  Ext.XE=G->XE ;
}


template<int DIM>
void AdaptiveData<DIM>::DeAttach() {
  assert(G!=NULL) ;
  G->ResizeGroup.Remove(this) ;
  G->RefineGroup.Remove(this) ;
  G=NULL  ;
  Ext.WC=NULL ;
  Ext.XA=Ext.XE=NULL ;
}

template<int DIM>
size_t AdaptiveData<DIM>::Size() {
  if (G==NULL) std::cout <<"AdaptiveData not attached to a grid\n";
  else return G->Size() ;
}

template<int DIM>
size_t AdaptiveData<DIM>::AllocatedEntries() {
  int reserved,used ;
  ba->Memory(&reserved,&used) ;
  return (size_t) reserved    ;
}

template<int DIM>
void AdaptiveData<DIM>::SetRefine() {
  assert(G!=NULL) ;
  G->ResizeGroup.Remove(this) ;
  G->RefineGroup.Insert(this) ;
}

template<int DIM>
void AdaptiveData<DIM>::SetBoundaryConditions(AdaptiveData<DIM> *X) {
  Ext.SetBoundaryConditions(X->Ext.BC) ;
  G->SetPeriodicConditions(Ext.BC) ;
}

template<int DIM>
void AdaptiveData<DIM>::SetBoundaryConditions(int bc[DIM][2]) {
  Ext.SetBoundaryConditions(bc) ;
  G->SetPeriodicConditions(Ext.BC) ;
}

template<int DIM>
void AdaptiveData<DIM>::SetBoundaryConditions(int d,int *bc) {
  Ext.SetBoundaryConditions(d,bc) ;
}
template<int DIM>
void AdaptiveData<DIM>::SetBoundaryConditions(int d,int bc0,int bc1) {
  Ext.SetBoundaryConditions(d,bc0,bc1) ;
}


template<int DIM>
bool AdaptiveData<DIM>::SameBoundaryConditions(AdaptiveData<DIM> *X) {
  return Ext.SameBoundaryConditions(&X->Ext) ;
}

template<int DIM>
void AdaptiveData<DIM>::CheckBC(AdaptiveData<DIM> *X) {
  assert(SameBoundaryConditions(X)) ;
}

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

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

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

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

template<int DIM>
void AdaptiveData<DIM>::CheckAdaptiveGrid(AdaptiveData<DIM> *X) 
{ assert (G==X->G) ;}

template<int DIM>
void AdaptiveData<DIM>::CheckAdaptiveGrid(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y) 
{ assert ((G==X->G)&&(G==Y->G)) ;}

template<int DIM>
void AdaptiveData<DIM>::CheckAdaptiveGrid(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z) 
{ assert ((G==X->G)&&(G==Y->G)&&(G==Z->G)) ;}

template<int DIM>
void AdaptiveData<DIM>::CheckAdaptiveGrid(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z , AdaptiveData<DIM> *Q) 
{ assert ((G==X->G)&&(G==Y->G)&&(G==Z->G)&&(G==Q->G)) ;}

template<int DIM>
void AdaptiveData<DIM>::CheckAdaptiveGrid(AdaptiveData<DIM> *X , AdaptiveData<DIM> *Y , AdaptiveData<DIM> *Z , AdaptiveData<DIM> *Q , AdaptiveData<DIM> *R) 
{ assert ((G==X->G)&&(G==Y->G)&&(G==Z->G)&&(G==Q->G)&&(G==R->G)) ;}


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

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

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

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

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

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


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

template<int DIM>
void  AdaptiveData<DIM>::Set(int *l , const double a) { 
  ba->Set(l,a) ; 
}

template<int DIM>
void AdaptiveData<DIM>::Set(int l , const double a) {
  AdaptiveDataArray<DIM>::VectorArray::iterator it(&ba->d) ;
  for (it=ba->d.begin(); it<= ba->d.end(); ++it) {
    bool set=false ;
    for (int i=0; i<DIM; i++) if (it.i[i]>=l) set=true ;
    if (set) ba->Set(it.i,a) ; 
  }
}

// Get pointer to a particular value 
template<int DIM>
double *AdaptiveData<DIM>::GetP(int *l,int *t) {
  AdaptiveDataArray<DIM>::DataVector *v=ba->d[l] ;
  long second = *(G->FindOrInsert(G->BasisPG[0],0, l , t )) ;
  if ((second)==MAXINT) return NULL ;
  return &(*v)[second] ;
}

// Set particular value 
template<int DIM>
int AdaptiveData<DIM>::Set(int *l,int *t,const double a) {
  double *p=GetP(l,t) ;
  if (p==NULL) return -1 ;
  (*p)=a ;
  return 0 ;
}

// Get particular value 
template<int DIM>
double AdaptiveData<DIM>::Get(int *l,int *t) {
  double *p=GetP(l,t) ;
  if (p==NULL) exit(-1) ;
  return *p ;
}


template<int DIM>
void AdaptiveData<DIM>::Set(int dir, int k, int s, const double a) {
  AdaptiveGrid<DIM>::ProjectionGrid           *pg=G->BasisPG[dir] ;
  AdaptiveGrid<DIM>::ProjectionGrid::iterator  pgit               ;
  
  int    i,j,l[DIM] ;

  for (pgit=(*pg).begin(); pgit!=(*pg).end(); ++pgit) {
    AdaptiveLevelsIndex *ali=(*pgit).first  ;
    AdaptiveLevels      *als=(*pgit).second ;
   
    j=0 ;
    for (i=0; i<DIM; i++)
      if (i!=dir) {
        l[i]=ali->l[j] ;
        j++ ;
      }
    l[dir]=k ;

    AdaptiveLevels::AdaptiveL *al=(*als)[k] ;
    AdaptiveLevels::AdaptiveL::iterator alit=(*al).find(s) ;
    assert(alit!=(*al).end()) ;

    j=(*alit).second ;
    int rl[DIM],rt[DIM] ;
    for (i=0; i<DIM; i++) rl[i]=ali->l[i], rt[i]=ali->t[i] ;
    (*ba->d[l])[j] = a ;

   }   // next pgit
}

template<int DIM>
void AdaptiveData<DIM>::Get(int *l,int *t,int dir, AdaptiveB *B, int leftrightall) {
  int K0,K1,ld,i ;
  
  AdaptiveLevelsIndex *lt=AdaptiveLevelsIndexAllocator.NewItem() ;
  lt->dim=DIM-1 ;
  for (i=0    ; i<dir; i++) lt->l[i  ]=l[i], lt->t[i  ]=t[i] ;
  for (i=dir+1; i<DIM; i++) lt->l[i-1]=l[i], lt->t[i-1]=t[i] ;

  AdaptiveGrid<DIM>::ProjectionGrid          *pg=G->BasisPG[dir]  ;
  AdaptiveGrid<DIM>::ProjectionGrid::iterator plt=(*pg).find(lt) ;
 
  // if no such 1-dimensional line, return empty AdaptiveB
  if (plt==(*pg).end()) {
    for (ld=G->Level0[dir]; ld<=LMAX; ld++) B->L[ld].IS->nr=0 ;
    std::cout<< "AdaptiveData<"<<DIM<<">::Get:  attempt to load AdaptiveB* from non-existent AdaptiveLevels\n" ;
    lt->Print() ;
    assert(0) ;
    return ;
  }
  AdaptiveLevelsIndexAllocator.DeleteItem(lt) ;

  // otherwise load values
  AdaptiveLevels *als=(*plt).second ;
  for (ld=G->Level0[dir]; ld<=LMAX; ld++) {
    l[dir]=ld ;

    K0=-1, K1=(1<<ld)+1 ;
    if (ld > B->Level0) {
      if (leftrightall & 1) K0=                 G->WC->KW0 ; 
      if (leftrightall & 2) K1=(1<<(ld-1)) -1 - G->WC->KW1 ; 
    } else { 
      if (leftrightall & 1) K0=                 G->WC->K0 ;
      if (leftrightall & 2) K1=(1<<ld)        - G->WC->K1 ; 
    }
    
    if (leftrightall==0) B->L[ld].IS->Load        ((*als)[ld] , ba->d[l] , B->L[ld].d ) ;
    else                 B->L[ld].IS->LoadBoundary((*als)[ld] , ba->d[l] , B->L[ld].d , K0 , K1 ) ;
  }
}


template<int DIM>
void AdaptiveData<DIM>::Set(int *l,int *t,int dir, AdaptiveB *B, int leftrightall, bool add) {
  int K0,K1,ld,i ;
  
  AdaptiveLevelsIndex *lt=AdaptiveLevelsIndexAllocator.NewItem() ;
  lt->dim=DIM-1 ;
  for (i=0    ; i<dir; i++) lt->l[i  ]=l[i], lt->t[i  ]=t[i] ;
  for (i=dir+1; i<DIM; i++) lt->l[i-1]=l[i], lt->t[i-1]=t[i] ;

  AdaptiveGrid<DIM>::ProjectionGrid          *pg=G->BasisPG[dir]  ;
  AdaptiveGrid<DIM>::ProjectionGrid::iterator plt=(*pg).find(lt) ;
 
  // if no such 1-dimensional line, return empty AdaptiveB
  if (plt==(*pg).end()) {
    std::cout<< "AdaptiveData<"<<DIM<<">::Get:  attempt to load AdaptiveB* from non-existent AdaptiveLevels\n" ;
    lt->Print() ;
    assert(0) ;
  }
  AdaptiveLevelsIndexAllocator.DeleteItem(lt) ;

  // otherwise load values
  AdaptiveLevels *als=(*plt).second ;
  int ldold=l[dir] ;
  for (ld=G->Level0[dir]; ld<=LMAX; ld++) {
    l[dir]=ld ;

    K0=-1, K1=(1<<ld)+1 ;
    if (ld > G->Level0[dir]) {
      if (leftrightall & 1) K0=                 G->WC->KW0 ; 
      if (leftrightall & 2) K1=(1<<(ld-1)) -1 - G->WC->KW1 ; 
    } else { 
      if (leftrightall & 1) K0=                 G->WC->K0 ;
      if (leftrightall & 2) K1=(1<<ld)        - G->WC->K1 ; 
    }

    if (leftrightall==0) B->L[ld].IS->Store        ((*als)[ld] , ba->d[l] , B->L[ld].d , add) ;
    else                 B->L[ld].IS->StoreBoundary((*als)[ld] , ba->d[l] , B->L[ld].d , K0 , K1 ,add) ;
  }
  l[dir]=ldold ;
}





template<int DIM>
void AdaptiveData<DIM>::ReadSlice(AdaptiveData<DIM+1> *X) {
  // checks
  assert (G->SliceParent=X->G) ;

  AdaptiveGrid<DIM>::ProjectionGrid              *pg=X->G->BasisPG[G->SliceDir] ;
  AdaptiveGrid<DIM>::ProjectionGrid::iterator     pgit             ;
  
  int    i,l[DIM+1] ;
  long   j ;

  // set Flags
  j=0 ;
  for (i=0; i<DIM+1; i++)
    if (i!=G->SliceDir) {
      Ext.islifting   [j]=X->Ext.islifting   [i] ;
      Ext.ismultiscale[j]=X->Ext.ismultiscale[i] ;
      Ext.BC[j][0]       =X->Ext.BC[i][0] ;
      Ext.BC[j][1]       =X->Ext.BC[i][1] ;
      j++ ;
    }

  for (pgit=(*pg).begin(); pgit!=(*pg).end(); ++pgit) {
    AdaptiveLevelsIndex *ali=(*pgit).first  ;
    AdaptiveLevels      *als=(*pgit).second ;
   
    j=0 ;
    for (i=0; i<DIM+1; i++)
      if (i!=G->SliceDir) {
        l[i]=ali->l[j] ;
        j++ ;
      }
    l[G->SliceDir]=G->Slicel ;

    AdaptiveLevels::AdaptiveL *al=(*als)[G->Slicel] ;
    AdaptiveLevels::AdaptiveL::iterator alit=(*al).find(G->Slicet) ;
    if (alit ==(*al).end()) continue;

    j=(*alit).second ;
    int rl[DIM],rt[DIM] ;
    for (i=0; i<DIM; i++) rl[i]=ali->l[i], rt[i]=ali->t[i] ;
    Set(rl , rt, (*X->ba->d[l])[j]) ;

   }   // next pgit
}  


template<int DIM>
void AdaptiveData<DIM>::WriteSlice(AdaptiveData<DIM+1> *X, bool add) {
  // checks
  assert(G->SliceParent=X->G) ;

  int  i ;
  long j=0 ;
  for (i=0; i<DIM+1; i++)
    if (i!=G->SliceDir) {
      assert(Ext.islifting   [j]==X->Ext.islifting   [i]) ;
      assert(Ext.ismultiscale[j]==X->Ext.ismultiscale[i]) ;
      assert(Ext.BC[j][0]       ==X->Ext.BC[i][0]) ;
      assert(Ext.BC[j][1]       ==X->Ext.BC[i][1]) ;
      j++ ;
    }

  AdaptiveGrid<DIM+1>::ProjectionGrid          *pg=X->G->BasisPG[G->SliceDir] ;
  AdaptiveGrid<DIM+1>::ProjectionGrid::iterator pgit ;
  
  int l[DIM+1] ;
  for (pgit=(*pg).begin(); pgit!=(*pg).end(); ++pgit) {
    AdaptiveLevelsIndex *ali=(*pgit).first  ;
    AdaptiveLevels      *als=(*pgit).second ;
   
    j=0 ;
    for (i=0; i<DIM+1; i++)
      if (i!=G->SliceDir) {
        l[i]=ali->l[j] ;
        j++ ;
      }
    l[G->SliceDir]=G->Slicel ;
   
    AdaptiveLevels::AdaptiveL *al=(*als)[G->Slicel] ;
    AdaptiveLevels::AdaptiveL::iterator alit=(*al).find(G->Slicet) ;
    if (alit==(*al).end()) continue;

    j=(*alit).second ;
    int rl[DIM],rt[DIM] ;
    for (i=0; i<DIM; i++) rl[i]=ali->l[i], rt[i]=ali->t[i] ;
    if (add) (*X->ba->d[l])[j] = (*X->ba->d[l])[j] + Get(rl , rt) ;
    else     (*X->ba->d[l])[j] =                     Get(rl , rt) ;
   }   // next pgit
}  



template<int DIM>
void  AdaptiveData<DIM>::Copy(AdaptiveData<DIM> *X) { 
  CopyExtensions(X) ;
  Ext.CopyData(&X->Ext) ;
  ba->Copy(X->ba) ;
}

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

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

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

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

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

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

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

template<int DIM>
void  AdaptiveData<DIM>::PPlus(const double a) {
  ba->PPlus(a) ;
  Ext.PPlus(a) ;
}


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

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

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

template<int DIM>
void  AdaptiveData<DIM>::Add(const double a,AdaptiveData<DIM> *X,const double b,AdaptiveData<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  AdaptiveData<DIM>::Add(const double a,AdaptiveData<DIM> *X,const double b,AdaptiveData<DIM> *Y,
                     const double c,AdaptiveData<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  AdaptiveData<DIM>::Add(const double a,AdaptiveData<DIM> *X,const double b,AdaptiveData<DIM> *Y,
                             const double c,AdaptiveData<DIM> *Z,const double d,AdaptiveData<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  AdaptiveData<DIM>::Add(const double a,AdaptiveData<DIM> *X,const double b,AdaptiveData<DIM> *Y,
                             const double c,AdaptiveData<DIM> *Z,const double d,AdaptiveData<DIM> *Q,
                             const double e,AdaptiveData<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  AdaptiveData<DIM>::MMinus(AdaptiveData<DIM> *X) { 
  Compatible(X) ;
  ba->MMinus( X->ba ) ;
  Ext.MMinus(&X->Ext) ;
}

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

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

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

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


template<int DIM>
void  AdaptiveData<DIM>::Mul(AdaptiveData<DIM>      *X ,AdaptiveData<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  AdaptiveData<DIM>::Mul(AdaptiveDataArray<DIM> *X ,AdaptiveDataArray<DIM> *Y) {
  ba->Mul(X , Y) ;
}
*/

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

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

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

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


//
// Norms
template<int DIM>
double AdaptiveData<DIM>::LpNorm(AdaptiveData<DIM> *Tmp, double p) {
  int i,GBC[2] ;
  AdaptiveData<DIM> *X=this ;
  
  for (i=0; i<DIM; 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],GBC) ;
    Tmp->ApplyOp(GBC, (i==0) ? X : 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 AdaptiveData<DIM>::LmaxNorm(AdaptiveData<DIM> *Tmp) { return LpNorm(Tmp,HUGE_VAL) ; }

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

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

// WriteUDF
template<int DIM>
void AdaptiveData<DIM>::WriteUDF(const char *name, int  *L, AdaptiveData<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 AdaptiveData<DIM>::WriteUDF(const char *name, int  L, AdaptiveData<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 AdaptiveData<DIM>::WriteUDF(const char *name, UniformData<DIM> *M, AdaptiveData<DIM> *Tmp, bool CastToFloat) {
  int  i,BG[DIM][2],L[DIM] ;
  AdaptiveData<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 ; 
    }
  }
   
  M->ba->CalcLevel(L) ;
  X->ToUniform(M,L) ;
  // inverse transform ?
  if (Tmp) 
    for (i=0; i<DIM; i++) M->ApplyOp(BG[i], M, i, INVERSETRANSFORM) ;

  M->WriteUDF(name,NULL,CastToFloat) ;
}


template<int DIM>
void  AdaptiveData<DIM>::Bining(double *eps,int Keps, int *count, AdaptivityCriterion *R) {

  AdaptiveDataArray<DIM>::VectorArray::iterator it(&ba->d),dend=ba->d.end() ;
  AdaptiveDataArray<DIM>::DataVector *v ;
  
  int i,n,j,a,e ;
  double c,eps0=R->eps ;

  for (i=-(Keps+1); i<=Keps; i++) count[i]=0 ;

  // loop over all levels
  for (it=ba->d.begin(); it<=dend; ++it) {
    v=ba->d[it]    ;
    n=(*v).size() ;

    // loop over coefficients of one level
    for (i=0; i<n; i++) {
      c = (*v)[i] ;

      // binary search for c in eps[]
      R->eps=eps[-Keps] ;
      if (!R->IsActive(c, it.i, G->Level0, DIM)) count[-(Keps+1)]++ ; // c<eps[-Keps]
      else {
        R->eps=eps[Keps] ;
        if (R->IsActive(c, it.i, G->Level0, DIM)) count[Keps]++ ;     // c>=eps[Keps]
        else { // eps[-Keps]<= c < eps[Keps]
          a=-Keps, e=Keps ;
          while (a<e-1) {
            j=(a+e)/2 ;
            R->eps=eps[j] ;
            if (!R->IsActive(c, it.i, G->Level0, DIM)) e=j ; // c<eps[j]
            else                                       a=j ;
          }
          // if e>a: eps[a] <= c < eps[e];  otherwise c=eps[a]=eps[e]
          count[a]++ ;
        } // else
      } // else
    } // next coefficient
  }   // next subspace

  R->eps=eps0 ;
}


//                       // 
//  from/to UniformData  //
//                       //
template<int DIM>
void AdaptiveData<DIM>::FromUniform(UniformData<DIM> *A,int *J) {
  Ext.CopyFlags(&A->Ext) ;
  Ext.CopyData (&A->Ext) ;
  G->SetPeriodicConditions(Ext.BC) ;
  for (int i=0; i<DIM; i++)
    if (!Ext.ismultiscale[i]) {
      std::cout<< "AdaptiveData<"<<DIM<<">::FromUniform: UniformData must be in multiscale representation\n" ;
      assert(0);
    }
  FromUniform(A->ba , J) ;
}


template<int DIM>
void AdaptiveData<DIM>::FromUniform(Matrix<double,DIM>* A , int *J) {
  int  i,a[DIM]={0},e[DIM],t1[DIM] ;
 
  G->SetFull(J) ;

  Matrix<double,DIM>::iterator l(G->Level0,J) ;
  for (l.Set(G->Level0); l<=J; ++l) {
    for (i=0; i<DIM; i++)
      e[i]=(l.i[i]==G->Level0[i]) ? 1<<l.i[i] : (1<<(l.i[i]-1))-1 ;
    
    Matrix<double,DIM>::iterator t(a,e) ;
    for (t.Set(a) ; t<=e ;++t) {
      for (i=0; i<DIM; i++) t1[i] = ( (l.i[i]==G->Level0[i]) ? 0 : (1<<(l.i[i]-1))+1 )  +  t.i[i] ;
      Set(l.i, t.i, (*A)[t1]) ;
    }
  }
}


template<int DIM>
void AdaptiveData<DIM>::ToUniform(UniformData<DIM> *A , int *J)  { 
 ToUniform(A->ba,J,0) ;
 A->Ext.CopyFlags(&Ext) ;
 A->Ext.CopyData (&Ext) ;
}

template<int DIM>
void AdaptiveData<DIM>::ToUniform(Matrix<double,DIM> *A , int *J, int flag)
{
 int i,l[DIM],dj[DIM] ; /*,sel=1*/ ; // sel / flag only for debug purposes
 size_t k,t[DIM],a[DIM] ; 

 A->Set(0.0) ;

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

 for (pgit=(*pg).begin(); pgit!=(*pg).end(); ++pgit) {
   AdaptiveLevelsIndex *ali=(*pgit).first  ;
   AdaptiveLevels      *als=(*pgit).second ;
   bool ins=true ;
   for (i=1; i<DIM; i++) { 
     l[i]=ali->l[i-1] ;
     t[i]=ali->t[i-1] ;
     if (l[i]>J[i]) ins=false ;
    }

   if (ins) 
     for (l[0]=G->Level0[0]; l[0]<=J[0]; l[0]++) {
       AdaptiveLevels::AdaptiveL           *al =(*als)[l[0]] ;
       AdaptiveLevels::AdaptiveL::iterator alit ;
       
       AdaptiveDataArray<DIM>::DataVector *v=ba->d[l] ;
       
       for (i=0;i<DIM;i++) a[i]=(l[i]==G->Level0[i]) ? 0 : (1<<(l[i]-1))+1 ;
       for (alit=(*al).begin(); alit!=(*al).end(); alit++) {
         t[0]=(*alit).first  ;
         k   =(*alit).second ;
         for (i=0;i<DIM;i++) dj[i]=(int)(t[i]+a[i]) ;
         (*A)[dj]= (flag) ? (double)(k+1) : (*v)[k] ;
       }
     }
 }
}




#include "AdaptiveGridRefine.hpp"
#include "Operations.hpp"

#endif

---