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@amroc.net

#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"

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 Symmetry[], int Periodic[], 
             int Eliminate, DataType EliminateValue, const int& MaxProcs) {

    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 {
            BlockList(level).push_back(storage_type(new data_block_type(i,j,k,ii,jj,kk,s,ss,sss),proc));  
            AMRreadData(filename,rank,dims,BlockList(level).back().first->data()); 
          }
          NoData = false;
        } 
   
      } while(result>=0);
      SDSclose(filename);
      if (!parData) break;
    }
    if (NoData) {
      std::cerr << "Nothing read for " << filename << " !" << std::endl;
      return -1; 
    }
    
    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

---