Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

BlendingData.hpp

//
// implementations for BlendingSpaces
//


#ifndef _BLENDINGSPACES_
# define _BLENDINGSPACES_

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

//
// General
//

template<int DIM>
BlendingData<DIM>::BlendingData() {
 a=NULL  ;
 A=NULL  ;
 Ext.InitXAE() ;
}


template<int DIM>
BlendingData<DIM>::BlendingData(BlendingWeights<DIM> *Ac) {
 a=NULL  ;
 Init(Ac) ;
 Ext.InitXAE() ;
}


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


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


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


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

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

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

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

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

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

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

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


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

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

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


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

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

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

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

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


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

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

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

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

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

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


//
// Set
//
template<int DIM>
void BlendingData<DIM>::SetFunction(Function *F) {

 Matrix<Matrix<double,DIM> * , DIM>::iterator it(a) ;
 for (it=(*a).begin(); it<=(*a).end(); ++it) if ((*(A->a))[it.i]  !=0 ) (*a)[it]->Set(F) ;
 
 for (int i=0; i<DIM; i++) {
   Ext.BC[i][0]=Ext.BC[i][1]=-1 ;
   Ext.ismultiscale[i]=false ;
   Ext.islifting   [i]=false ;
 }
 Ext.dimext=0 ;
}


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

 // nasty implementation: we allocate an auxiliary, dummy UniformData
 // and replace the actual array of numerical values by the arrays of *this
 // This way we may use UniformData::ApplyOp()
 
 int bg[DIM]={0},e[DIM]={0} ;
 UniformData<DIM> X(bg,e,Ext.WC,Ext.XA,Ext.XE),
                  Y(bg,e,Ext.WC,Ext.XA,Ext.XE) ;
 Matrix<double,DIM> *Xba=X.ba, *Yba=Y.ba ; // keep old arrays

 Y.Ext.CopyFlags(&S->Ext) ;

 // forall ProjectionSubSpaces
 Matrix<Matrix<double,DIM> * , DIM>::iterator il(A->G->Level0,A->G->L) ;
 for (il.Set(A->G->Level0); il<=(*a).end(); ++il) {
   if ((*A->a)[il.i] != 0) {
     // fake the arrays
     X.ba=(   *a)[il] ;
     Y.ba=(*S->a)[il] ;
     X.ApplyOp(BCT, &Y, dir, op) ;
   } 
 }

 Ext.CopyFlags(&X.Ext) ;
 Ext.CopyData (&X.Ext) ;
 // set old arrays to be deleted by ~X and ~Y
 X.ba=Xba, Y.ba=Yba ;


}


#endif

---