Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

Adaptive1D.hpp

//
//: Adaptive1D.hpp
// DatenStructuren fuer 1D schnelle adaptive Verfahren
//

#ifndef ADAPTIVE1D_H
# define ADAPTIVE1D_H

#include "Wavelet.hpp"
#include <map>
#include <vector>
#include "Refine.hpp"

#define _MAX_INDEXSET_ 1025 // max number of intervals which may be stored by IndexSet

struct intlt { bool operator()(const int i,const int j) const ; } ;


struct IndexSet { 
  typedef map<int,long,intlt>  AdaptiveLE ;

  int  nr,*a,*e ;

public:
       IndexSet()          ;
       IndexSet(int Size)  ;
      ~IndexSet()          ;

  int    Init(int Size)    ;

  int    Print(char *st=NULL) ;
  void   Print(char *st,int l,int Level0) ;
  size_t Size() ;
  void   Copy(IndexSet *F)   ;
  bool   IsSame(IndexSet *F) ;

  void Union     (IndexSet *A,IndexSet *B) ;
  void Difference(IndexSet *A,IndexSet *B) ;

  bool Contains  (int i,int *r=NULL) ;

  double VecMax   (double *t) ;
  double VecMaxAbs(double *t) ;  

  void VecAdd  (double *to,double *a, double *b) ;
  void VecSub  (double *to,double *a, double *b) ;
  void VecPlus (double *to,double *a) ;
  void VecMinus(double *to,double *a) ;
  void VecMinus(double *to,double  a)  ;
  void VecMul  (double *to,const double factor)   ;
  void VecCopy (double *to,double *from)    ;
  void VecFill (double *to,const double f)  ;
  void VecClear(double *to) ;
  double VecInnerProd(double *a,double *b) ;
  void VecPrint(const char *st,double *t)  ;
  void GenerateRandom(bool boundflag , int KW , int J , bool isW , bool dorandom) ;
     
  void VecApplyWENO(double *to, double *f, double *roespeed, int *BC, int J, int k , double DX) ;

  void Load  (AdaptiveLE *al) ;
  void LoadIS(AdaptiveLE *al) ;

  void Load (AdaptiveLE *al , vector<double> *v, double *d) ;

  void Store(AdaptiveLE *al , vector<double> *v, double *d, bool add=false) ;
  //void Store(AdaptiveLE *al , vector<double> *v, double *d,int p) ;
  void LoadBoundary  (AdaptiveLE *al , vector<double> *v, double *d,int Left,int Right) ;
  void LoadBoundaryIS(AdaptiveLE *al , int Left,int Right) ;

  void StoreBoundary (AdaptiveLE *al , vector<double> *v, double *d,int Left,int Right, bool add=false) ;
  //void WriteToAdaptiveLE(AdaptiveLE *al,int *counter) ;


  void CopyWithBounds     (IndexSet *B,int l, int *BC, Wavelets *W) ;
  void CalcIndexSet       (double *d,int Size,int K0,int K1,bool Insert0,bool Insert1) ;
  void ForBoundaryWavelets(bool left, bool right, IndexSet *I, int level, Wavelets *W) ;

  // scaling functions of level l which span the given wavelets of level l
  void InducedByWavelets            (IndexSet *B,int l, int *BC, Wavelets *W) ;
  // scalings(l) influenced by the given scalings(l-1)
  void InducedByCoarseScalings      (IndexSet *S,int l, int *BC, Wavelets *W) ;  
  // scalings(l-1) required to span the given scalings(l)
  void InducedByFineScalings        (IndexSet *S,int l, int *BC, Wavelets *W) ;
  // dual scalings(l-1) required to span the given dual scalings(l)
  void InducedByFineDualScalings    (IndexSet *S,int l, int *BC, Wavelets *W) ;
  // wavelets(l) which coincide with the given scalings(l)
  void WaveletsInducedByFineScalings(IndexSet *B,int l, int *BC, Wavelets *W) ;
  // dual wavelets(l) influenced by the given dual scalings(l)
  void WaveletsInducedByFineDualScalings(IndexSet *B,int l, int *BC, Wavelets *W) ;

  void InducedByOperatorMatrix      (IndexSet *S,int l, int *BC, OperatorMatrix *W) ;

  // scalings(l-1) required to span the given lifting wavelets(l)
  void ScalingsInducedByLifting     (IndexSet *S,int l, int *BC, Wavelets *W) ;


  void IndexSetForTypeIII       (IndexSet *S,int l,             int *BC, Wavelets *W) ;
  
  //
  // Let A be a XxY matrix, M subset Y, then 
  //                C(A,M):={ x in X | exists y in M: A_{xy} !=0 } 
  //
  // A is given by (*H) and (more important) by the functions First/LastNZE(P)
  //
  // CAM is used mainly to DETECT dependencies -> overestimation of intervals is valid
  //
  // It is assumed that A is a banded matrix.
  // First(j) returns the number of the smallest i such that A_{i,j}!= 0
  // Last (j) returns the number of the largest  i such that A_{i,j}!= 0
  // 
  // in case of the TransformMatrices (No Transpose) some special care in the 
  // definition of First/Last is required:
  // First: j<=2i+O  -> First=(j-O+1)/2
  // Last : 2i-U<=j  -> Last =(j+U)/2 
  // 
  void CAM(IndexSet *M, int *BC, 
           MatrixShape *H,
           int (*FirstNZEP)(int,MatrixShape *) , int (*LastNZEP )(int,MatrixShape *) ,
           int (*FirstNZE )(int,MatrixShape *) , int (*LastNZE  )(int,MatrixShape *)) ;

  //
  //                E(A,M):={x in X | forall y: A_{xy}!=0 there holds y in M}
  //
  // EAM is used mainly to find COMPLETELY dependent intervals  -> UNDERestimation of intervals 
  // is valid
  //
  // in case of the TransformMatrices (No Transpose) some special care in the 
  // definition of First/Last is required:
  // First: 2i+O<=j  -> First=(j-O  )/2
  // Last : j<=2i-U  -> Last =(j+U+1)/2
  // 
  void EAM(IndexSet *M, int *BC, 
           MatrixShape *H,
           int (*FirstNZEP)(int,MatrixShape *) , int (*LastNZEP )(int,MatrixShape *) ,
           int (*FirstNZE )(int,MatrixShape *) , int (*LastNZE  )(int,MatrixShape *)) ;

  //
  // set operations

  void SatisfyBoundCondition(int M,int N,int BL0, int BR0,bool *Insert0, bool *Insert1) ;

} ;


struct AdaptiveGS ;
struct AdaptiveBOperator ;

struct AdaptiveB {
   
   int Level0 ;
   int J ;
   int BC[2] ;
 
   double XA ;
   double XE ;
   bool ismultiscale ;
   bool islifting ;
  
   struct AdaptiveLevel {
     bool     ISowner ; 
     IndexSet *IS     ; 
     int      Size    ;
     double   *d      ;

   public:    
              AdaptiveLevel() ;
             ~AdaptiveLevel() ;

     int      Init(int size)  ;
     int      InitD(int size) ;
     int      Print()         ;
     int      Print(char *st,int level,int allflag) ;
     friend class AdaptiveGS  ;
  } L[LMAX+1] ;

public:
       AdaptiveB() ;
       AdaptiveB(AdaptiveB *B) ;
      ~AdaptiveB() ;

  int  Init(int Level0,int J,double xa=0, double xe=0) ;
  void Copy(AdaptiveB *B) ;
  void CopyIndexSets(AdaptiveB *B) ;
  bool SameIndexSets(AdaptiveB *B) ;
  void SetBoundaryConditions(int *BCs)        ;
  void SetBoundaryConditions(int bc0,int bc1) ;

  void Print(int allflag) ;
  void Print(char *st=NULL,int allflag=0) ;
  void Print(char *st) ;
  size_t Size() ;
  
  void FromFull(double *d,int *BC , Wavelets *W , double eps) ;
  void FromFull(double *d,int *BC , Wavelets *W , AdaptivityCriterion *R) ;
  void FromFull(double *d) ;
  void ToFull(double *d) ;

  void FromAdaptiveGS (AdaptiveGS *G , Wavelets *W) ;
  void FromAdaptiveGS2(AdaptiveGS *G , Wavelets *W) ;
  void FromAdaptiveGS3(AdaptiveGS *G , Wavelets *W) ;
  void FromAdaptiveGS4(AdaptiveGS *G , Wavelets *W) ;
  void AdditiveFromAdaptiveGS (AdaptiveGS *G , Wavelets *W) ;
  void HighPassFilterFromAdaptiveGS (AdaptiveGS *G , Wavelets *W) ;

  void IndexSetFromAdaptiveGS(AdaptiveGS *G , Wavelets *W) ;

  void LoadFromInterpoletGS       (AdaptiveGS *G , Wavelets *W) ;
  void LoadFromLiftingInterpoletGS(AdaptiveGS *G , Wavelets *W) ;

  //void copy(AdaptiveB *X) ;
  double Max() ;
  double MaxAbs() ;
  void Sub(AdaptiveB *X,AdaptiveB *Y) ;
  void Add(AdaptiveB *X,AdaptiveB *Y) ;
  void Add(const double a,AdaptiveB *X,const double b,AdaptiveB *Y) ;
  void Add(const double a,AdaptiveB *X,const double b,AdaptiveB *Y,const double c,AdaptiveB *Z) ;

  void PPlus(AdaptiveB *X) ;
  void PPlus(const double a,AdaptiveB *X) ;
  void PPlus(AdaptiveB *X,AdaptiveB *Y) ;
  void PPlus(const double a,AdaptiveB *X,const double b,AdaptiveB *Y) ;
  
  void MMinus(AdaptiveB *X) ;
  void MMinus(const double a,AdaptiveB *X) ;
  double InnerProd(AdaptiveB *X) ;
  void MMul(const double a) ;

  // dummy function, does nothing
  void CopyExtensions(AdaptiveB *X) ; 

  void Fill(double a) ;
  
  void Refine(AdaptiveB *B , AdaptivityCriterion *R,int *l,int dir,int Dim,Wavelets *W,bool copyflag=false) ;
  void Refine(AdaptiveB *B , AdaptivityCriterion *R,Wavelets *W) ;
  void Refine(AdaptivityCriterion *R,Wavelets *W) ;
  bool CheckBoundCondition (Wavelets *W) ;
  void InterpoletRegularity(Wavelets *W) ;

  //  
  // active operations: ApplyOp , Multiply,  ApplyFD, ApplyWENO
  // result in *this,  
  // for  ApplyOp and Multiply the  IndexSets of u are used
  // for ApplyXXX the BC of  (*this)  are the To-RB !   
  //
  // GStmp must be different from GS and is just used as temporary buffer, must be initialized, but no
  // valid values are required
  // 
  void ApplyOp  (AdaptiveB *u , int op , Wavelets *W, AdaptiveGS *GS, AdaptiveGS *GStmp) ;
  void ApplyOp  (AdaptiveB *u , int op , Wavelets *W                ) ;
  void ApplyFD  (AdaptiveB *u , int op , Wavelets *W, AdaptiveGS *GS) ;
  void ApplyWENO(AdaptiveB *u , int op , Wavelets *W, AdaptiveGS *GS, AdaptiveGS *GRoeSpeed) ;

  // GS must be initialized, but no IndexSets have to be set
  void ApplyFD2 (AdaptiveB *u , int op , Wavelets *W, AdaptiveGS *GS) ;

  void Multiply  (AdaptiveB *u , AdaptiveB *v, Wavelets *W,int which) ;
  void Projection(AdaptiveB *u , Wavelets *W , AdaptiveGS *GS, AdaptiveGS *GStmp) ;

  int  Solve(AdaptiveB *RHS, AdaptiveBOperator  *Op) ;

  // special for div*grad operators
  void RestrictPressure() ;

  friend class AdaptiveGS ;

  doubleBuffer *Buffers[2] ;
} ;
#ifdef _TYPEDEF_STRUCTS_
typedef struct AdaptiveB AdaptiveB;
#endif

struct AdaptiveGS {

  int    Level0,J,BC[2] ;
  double XA,XE ;

  struct AdaptiveLevel {
    IndexSet I,II,III,IV ;
    int      Size        ;
    double   *d          ;

public:
          AdaptiveLevel() ;
         ~AdaptiveLevel() ;

     int  Init(int size) ;
     int  Print() ;
     int  Print(char *st,int level,int allflag) ;
  } L[LMAX+1] ;
  
public:
       AdaptiveGS() ;
      ~AdaptiveGS() ;
  int  Init(int Lev0,int J) ;
  void Print()              ;
  void Print(int allflag)   ;
  void Print(char *st,int allflag) ;

  void IndexSetForAdaptiveGS  (AdaptiveB *B , Wavelets *W) ; //with    III, using IV
  void IndexSetForAdaptiveGS2 (AdaptiveB *B , Wavelets *W) ; //without III
  void IndexSetForAdaptiveGSOp(AdaptiveB *X , Wavelets *W) ;
  void IndexSetForAdaptiveGSFD(AdaptiveB *X , OperatorMatrix *FD , Wavelets *W) ;

  //
  // Get adaptive generating system.
  //

  void FromAdaptiveB         (AdaptiveB *B , Wavelets *W) ; 
  void FromAdaptiveB1        (AdaptiveB *B , Wavelets *W) ;
  void Projection            (int *BCTo , Wavelets *W)    ; // project V->(V_0) 

  void BoundaryQuadrature    (Wavelets  *W , int which)   ;
  void ApplyFD               (int *BCTo , AdaptiveGS *G , Wavelets *W , int op) ;
  void ApplyWENO             (int *BCTo , AdaptiveGS *G , AdaptiveGS *GRoeSpeed, int op) ;

  void RestoreFromInterpoletBasis(AdaptiveB *B , Wavelets *W) ;

  void RestoreFromInterpoletBasis(AdaptiveB *B , AdaptiveGS *G , Wavelets *W) ;

  void CopyIndexSets(AdaptiveGS *G) ;
  void Copy         (AdaptiveGS *G) ;

  bool SameIndexSetsI  (AdaptiveGS *G) ;
  bool SameIndexSetsII (AdaptiveGS *G) ;
  bool SameIndexSetsIII(AdaptiveGS *G) ;
  bool SameIndexSetsIV (AdaptiveGS *G) ;
  bool SameIndexSets   (AdaptiveGS *G) ; // I && II && III && IV

  // some linear Algebra Stuff
  void   MMinus   (AdaptiveGS *G) ;
  double InnerProd(AdaptiveGS *G) ;

  doubleBuffer *Buffers[2] ;
} ;
#ifdef _TYPEDEF_STRUCTS_
typedef struct AdaptiveGS AdaptiveGS;
#endif

// generalized Helmholtz-operator
struct AdaptiveBOperator {

       AdaptiveBOperator(double aa0,double aa1,Wavelets *W) ;

  void apply         (AdaptiveB *X, AdaptiveB *R) ;
  void Preconditioner(AdaptiveB *X, AdaptiveB *R) ;
  // data
  double a0,a1 ; // (op u,v) == a0*(u,v) + a1*(grad u,grad v)
  Wavelets *WC ;
} ;

#endif

---