amroc/converters/HierarchyStorage.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding
// Brandenburgische Universitaet Cottbus
//
// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@cacr.caltech.edu

#ifndef AMROC_HIERARCHYSTORAGE_H
#define AMROC_HIERARCHYSTORAGE_H

#include <cassert>
#include <unistd.h>
#include <vector>
#include <utility>

#include "GridData3.h" 
#include "DAGHIO_hdf_ncsa.h" 
#include "BBoxList.h"

#define MAXLEV 11

template <class DataType>
class HierarchyStorage {
public:
  typedef GridData<DataType,3> data_block_type;
  typedef std::pair<data_block_type*,int> storage_type;
  typedef std::vector<storage_type> block_list_type;

public:
  HierarchyStorage(int nl) : _NLevels(nl), _RealTime(0.0) {
    _BlockList = new block_list_type*[_NLevels];
    for (int i=0; i<_NLevels; i++)
      _BlockList[i] = new block_list_type;
  }

  ~HierarchyStorage() {
    for (int i=0; i<_NLevels; i++) {
      for (typename block_list_type::iterator bit=BlockList(i).begin();
           bit!=BlockList(i).end(); bit++)
        delete (*bit).first;
      BlockList(i).clear();
      delete _BlockList[i];
    }
    delete _BlockList;
  }

  int ReadIn(char* fname, int RefineFactor[], Coords DimUsed, int Symmetry[], int Periodic[], 
             int Eliminate, DataType EliminateValue, const int& MaxProcs, 
             BBox *_where, BBox *_to, BBox *_from, const int& RestrictData) {

    char dataname[256];
    char filename[256];
    int32 level,gridID,step,rank,proc;
    int result;
    int32 ori[3],dims[3],res[4];
    int i,j,k,ii,jj,kk,s,ss,sss,maxlevel=0;
    int maxsize[3] = { 0, 0, 0};
    short olap[3] = { 0, 0, 0};
    bool parData=false, NoData=true;

    std::sprintf(filename,"%s.hdf",fname);
    if (access(filename, F_OK)) parData = true;

    for (int me=0; me<MaxProcs; me++) {
      if (parData) {
        std::sprintf(filename,"%s.%d.hdf",fname,me);
        if (access(filename, F_OK)) continue;
      }
      if (access(filename, R_OK)) {
        std::cerr << "No read permission on " << filename << " !" << std::endl;
        return -1; 
      }
      if (access(filename, W_OK)) {
        std::cerr << "No write permission on " << filename << " !"
                  << std::endl;
        return -1; 
      }
      
      do{
      
        result=AMRreadAttribs(filename,dataname,&proc,&level,&gridID,&step,
                              &_RealTime,&rank,dims,ori,res);

        if(result>=0){
          if (SDSgetNT(filename)!=DFNT_FLOAT32 && sizeof(DataType)==4) {
            std::cerr << "Numbertype must be float!" << std::endl;
            return -1; 
          }
          if (SDSgetNT(filename)!=DFNT_FLOAT64 && sizeof(DataType)==8) {
            std::cerr << "Numbertype must be double!" << std::endl;
            return -1; 
          }
          
          i=ori[0];
          j=ori[1];
          k=ori[2];
          s=res[0];
          ss=res[1]; 
          sss=res[2];
          ii=i+s*(dims[0]-1);
          jj=j+ss*(dims[1]-1);
          kk=k+sss*(dims[2]-1);
          
          if(level>maxlevel) maxlevel=level;
          if(ii+s>maxsize[0]) maxsize[0]=ii+s;
          if(jj+s>maxsize[1]) maxsize[1]=jj+s;
          if(kk+s>maxsize[2]) maxsize[2]=kk+s;
          
          if (Eliminate) {
            BBox tmpbb(3,i,j,k,ii,jj,kk,s,ss,sss);
            Coords step(3,s,ss,sss);
            data_block_type tmpblock(tmpbb);
            AMRreadData(filename,rank,dims,tmpblock.data()); 
            
            BBoxList cbbl, obbl; 
            BeginFastIndex3(tmp, tmpblock.bbox(), tmpblock.data(), DataType);
            for_3 (n,m,l,tmpblock.bbox(),step)
              if (FastIndex3(tmp,n,m,l)==EliminateValue)
                cbbl.add(BBox(3,n,m,l,n,m,l,s,ss,sss));
            end_for
            EndFastIndex3(tmp); 

            obbl.add(tmpblock.bbox());
            obbl -= cbbl;
            obbl.mergeboxes(olap);
            for (BBox *bb=obbl.first();bb;bb=obbl.next()) {
              BlockList(level).push_back(storage_type(new data_block_type(*bb),proc));
              BlockList(level).back().first->equals(tmpblock,*bb);
            }
          }
          else {
            BBox tmpbb(3,i,j,k,ii,jj,kk,s,ss,sss);
            Coords step(3,s,ss,sss);
            BlockList(level).push_back(storage_type(new data_block_type(tmpbb),proc));  
            AMRreadData(filename,rank,dims,BlockList(level).back().first->data()); 
            BeginFastIndex3(tmp, tmpbb, BlockList(level).back().first->data(), DataType);
            for_3 (n,m,l,tmpbb,step)
              if (FastIndex3(tmp,n,m,l)==EliminateValue)
                FastIndex3(tmp,n,m,l) = FLT_MAX;
            end_for
            EndFastIndex3(tmp);             
          }
          NoData = false;

          if (BlockList(level).size()>0) {
            if (!_where[level].empty()) {
              BBox where(_where[level]);
              Coords ss(BlockList(level).back().first->stepsize());
              for (int i=0;i<3;i++) {
                if (where.stepsize(i)>ss(i))
                  where.refine(i,where.stepsize(i)/ss(i),0);
                else where.coarsen(i,ss(i)/where.stepsize(i));
              }
              BBox tmpbb = BlockList(level).back().first->bbox()*where;
              std::cerr << ".";
              data_block_type* tmpblock = new data_block_type(tmpbb);
              tmpblock->equals(*(BlockList(level).back().first));
              delete BlockList(level).back().first;
              BlockList(level).back().first = tmpblock;
            }
            else {
              unsigned int cp=BlockList(level).size()-1;
              Coords ss(BlockList(level).at(cp).first->stepsize());
              for (register int lev=0; lev<MAXLEV; lev++) {         
                if (!_to[lev].empty() && !_from[lev].empty() && ss==_from[lev].stepsize()) {
                  BBox to(_to[lev]), from(_from[lev]);
                  BBox bb(BlockList(level).at(cp).first->bbox());
                  Coords cs;
                  for (int i=0;i<3;i++)
                    if (from.extents(i)/to.extents(i)>1)
                      cs(i) = from.extents(i)/to.extents(i);
                    else cs(i) = 1;
                  
                  data_block_type* tmpblock = new data_block_type(coarsen(bb,cs));
                  BeginFastIndex3(src, bb, BlockList(level).at(cp).first->data(), DataType);
                  BeginFastIndex3(tmp, tmpblock->bbox(), tmpblock->data(), DataType);
                  for_3 (n,m,l,tmpblock->bbox(),tmpblock->stepsize())
                    FastIndex3(tmp,n,m,l)=0.;
                    BBox bbadd(3,n,m,l,n,m,l,tmpblock->stepsize(0),tmpblock->stepsize(1),tmpblock->stepsize(2));
                    bbadd.refine(cs); bbadd*=bb;
                    for_3 (nn,mm,ll,bbadd,bbadd.stepsize())
                      FastIndex3(tmp,n,m,l)+=FastIndex3(src,nn,mm,ll);
                    end_for
                    FastIndex3(tmp,n,m,l)/=bbadd.size();
                  end_for
                  EndFastIndex3(tmp); 
                  EndFastIndex3(src); 
                
                  BBox newbb(tmpblock->bbox());
                  Coords ex = newbb.extents();
                  newbb.stepsize() = to.stepsize();
                  newbb.lower() += to.lower()-from.lower();
                  newbb.lower() /= (from.extents()*from.stepsize())/(to.extents()*to.stepsize());
                  newbb.upper() = newbb.lower()-newbb.stepsize()+ex*newbb.stepsize();
                  std::cerr << ".";
                  DataType *data;
                  tmpblock->deallocate(data);
                  delete tmpblock;
                  delete BlockList(level).at(cp).first;
                  BlockList(level).erase(BlockList(level).begin()+cp);
                  data_block_type* newblock = new data_block_type(newbb,data);
                  BlockList(lev).push_back(storage_type(newblock,proc));
                  if (lev>maxlevel) maxlevel=lev;
                }
              }
            }
          }
        } 
   
      } while(result>=0);
      SDSclose(filename);
      if (!parData) break;
    }
    if (NoData) {
      std::cerr << "Nothing read for " << filename << " !" << std::endl;
      return -1; 
    }
    
    if (RestrictData != 0) {
      for (level=maxlevel;level>0; level--){    
        std::cerr << "Restrict level " << level;
        typename block_list_type::iterator bit;
        if (BlockList(level-1).empty()) {
          BBoxList bblc;
          for (bit=BlockList(level).begin();bit!=BlockList(level).end();bit++)
            bblc.add((*bit).first->bbox());
          short olap[DAGHMaxRank];
          for (int i=0;i<DAGHMaxRank;i++) olap[i]=0;
          bblc.mergeboxes(olap);
          for (BBox *_b = bblc.first();_b;_b=bblc.next()) {
            for (int i=0;i<3;i++)
              if (DimUsed(i)) _b->coarsen(i,RefineFactor[level-1]);
            BlockList(level-1).push_back(storage_type(new data_block_type(*_b),0));  
          }
        }  
        for (bit=BlockList(level).begin();bit!=BlockList(level).end();bit++) {
          BBox bb = (*bit).first->bbox();
          typename block_list_type::iterator bitc;
          for (bitc=BlockList(level-1).begin();bitc!=BlockList(level-1).end();bitc++) {
            BBox bbc = (*bitc).first->bbox();
            Coords cs = bbc.stepsize()/bb.stepsize();
            BBox bbwork(bbc * coarsen(bb,cs));
            if (!bbwork.empty()) {
              std::cerr << ".";
              BeginFastIndex3(src, bb, (*bit).first->data(), DataType);
              BeginFastIndex3(dst, bbc, (*bitc).first->data(), DataType);
              for_3 (n,m,l,bbwork,bbwork.stepsize())
                FastIndex3(dst,n,m,l)=0.;
                BBox bbadd(3,n,m,l,n,m,l,bbwork.stepsize(0),bbwork.stepsize(1),bbwork.stepsize(2));
                bbadd.refine(cs); bbadd*=bb;
                for_3 (nn,mm,ll,bbadd,bbadd.stepsize())
                  FastIndex3(dst,n,m,l)+=FastIndex3(src,nn,mm,ll);
                end_for
                FastIndex3(dst,n,m,l)/=bbadd.size();
              end_for
              EndFastIndex3(dst); 
              EndFastIndex3(src); 
            }
          }
        }
        std::cerr << std::endl;
      }
    }

    for (int d=0; d<3; d++) 
      if (Symmetry[d] != 0) 

        for (level=0;level<=maxlevel; level++){      
          block_list_type TempStorage;          
          typename block_list_type::iterator bit;
          for (bit=BlockList(level).begin();bit!=BlockList(level).end();bit++)
          {
            if (Symmetry[d] > 0) {
              BBox bbwork = (*bit).first->bbox();
              bbwork.lower(d)=2*maxsize[d]-(*bit).first->bbox().upper(d)-
                (*bit).first->bbox().stepsize(d);
              bbwork.upper(d)=2*maxsize[d]-(*bit).first->bbox().lower(d)-
                (*bit).first->bbox().stepsize(d); 
              TempStorage.push_back(storage_type(new data_block_type(bbwork),(*bit).second));
         
              BeginFastIndex3(source,(*bit).first->bbox(),
                              (*bit).first->data(), DataType); 
              BeginFastIndex3(target, TempStorage.back().first->bbox(),
                              TempStorage.back().first->data(), DataType);
              switch(d) {
              case 0:
                for_3 (n, m, l,(*bit).first->bbox(), (*bit).first->bbox().stepsize())
                  FastIndex3(target, 2*maxsize[0]-n-(*bit).first->bbox().stepsize(0),
                             m,l) = FastIndex3(source,n,m,l);  
                end_for
                break;
              case 1:
                for_3 (n, m, l,(*bit).first->bbox(), (*bit).first->bbox().stepsize())
                  FastIndex3(target,n, 2*maxsize[1]-m-(*bit).first->bbox().stepsize(1),
                             l) = FastIndex3(source,n,m,l);  
                end_for
              break;
              default:
                for_3 (n, m, l,(*bit).first->bbox(), (*bit).first->bbox().stepsize())
                  FastIndex3(target,n,m,2*maxsize[2]-l-(*bit).first->bbox().stepsize(2)
                             ) = FastIndex3(source,n,m,l);  
                end_for
              break;
              }
              EndFastIndex3(source); 
              EndFastIndex3(target);  
            }

            else {
              BBox bbwork = (*bit).first->bbox();
              bbwork.lower(d)=maxsize[d]+(*bit).first->bbox().lower(d);
              bbwork.upper(d)=maxsize[d]+(*bit).first->bbox().upper(d); 
              TempStorage.push_back(storage_type(new data_block_type(bbwork),(*bit).second));
              TempStorage.back().first->copy(
                *((*bit).first),TempStorage.back().first->bbox(),
                (*bit).first->bbox());

              bbwork = (*bit).first->bbox();
              bbwork.lower(d)=maxsize[d]-(*bit).first->bbox().upper(d)-
                (*bit).first->bbox().stepsize(d);
              bbwork.upper(d)=maxsize[d]-(*bit).first->bbox().lower(d)-
                (*bit).first->bbox().stepsize(d); 
              TempStorage.push_back(storage_type(new data_block_type(bbwork),(*bit).second));

              BeginFastIndex3(source,(*bit).first->bbox(),
                              (*bit).first->data(), DataType); 
              BeginFastIndex3(target, TempStorage.back().first->bbox(),
                              TempStorage.back().first->data(), DataType);
              switch(d) {
              case 0:
                for_3 (n, m, l, (*bit).first->bbox(),(*bit).first->bbox().stepsize())
                  FastIndex3(target, maxsize[0]-n-(*bit).first->bbox().stepsize(0), 
                             m,l) = FastIndex3(source,n,m,l);  
                end_for
              break;
              case 1:
                for_3 (n, m, l, (*bit).first->bbox(),(*bit).first->bbox().stepsize())
                  FastIndex3(target,n, maxsize[1]-m-(*bit).first->bbox().stepsize(1),
                             l) = FastIndex3(source,n,m,l);  
                end_for
              break;
              default:
                for_3 (n, m, l, (*bit).first->bbox(),(*bit).first->bbox().stepsize())
                  FastIndex3(target,n,m,maxsize[2]-l-(*bit).first->bbox().stepsize(2)
                             ) = FastIndex3(source,n,m,l);  
              end_for
              break;
              }
              EndFastIndex3(source); 
              EndFastIndex3(target);  
            }
          }
          
          if (Symmetry[d] < 0) {
            for (bit=BlockList(level).begin();bit!=BlockList(level).end(); bit++) 
              delete (*bit).first;
            BlockList(level).clear();
          }

          for(bit=TempStorage.begin();bit!=TempStorage.end();bit++) 
            BlockList(level).push_back(storage_type((*bit).first,(*bit).second));
          TempStorage.clear();           
        }

      else if (Periodic[d] > 0)

        for (level=0;level<=maxlevel; level++){      
          block_list_type TempStorage;          
          typename block_list_type::iterator bit;
          for (bit=BlockList(level).begin();bit!=BlockList(level).end();bit++)
            for (int pc=1; pc<=Periodic[d]; pc++) {
              BBox bbwork = (*bit).first->bbox();
              bbwork.lower(d)=pc*maxsize[d]+(*bit).first->bbox().lower(d);
              bbwork.upper(d)=pc*maxsize[d]+(*bit).first->bbox().upper(d); 
              TempStorage.push_back(
                storage_type(new data_block_type(bbwork),(*bit).second));
              TempStorage.back().first->copy(
                *(*bit).first, bbwork, (*bit).first->bbox());
            }
                             
          for (bit=TempStorage.begin();bit!=TempStorage.end();bit++) 
            BlockList(level).push_back(
              storage_type((*bit).first,(*bit).second));
          TempStorage.clear();           
        }
    return (parData ? 1 : 0);
  }


  block_list_type& BlockList(int Level) {
    assert(Level >=0 && Level < _NLevels);
    return (*_BlockList[Level]);
  }
  
  inline const int& NLevels() const { return _NLevels; }
  inline int NLevels() { return _NLevels; }
  inline float64 RealTime() const { return _RealTime; }
  
private:
  int _NLevels;
  float64 _RealTime;
  block_list_type** _BlockList;
};


#endif

---