VisualGrid.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_VISUALGRID_H
#define AMROC_VISUALGRID_H

#include "VisualGridBase.h"
#include "HierarchyStorage.h" 
#include "VisualGridData.h"
#include "Evaluator.h"

#include "BBoxList.h"
#include "DAGHDefaults.h"

#include <vector>
#include <utility>
#include <map>
#include <functional>
#include <algorithm>
#include <iostream>
#include <string>
#include <cstdio>
#include <cfloat>

template <class VectorType>
class VisualGrid : public VisualGridBase {
  typedef typename VectorType::InternalDataType DataType;

public:
  typedef VisualGridData<VectorType> data_block_type;
  typedef std::vector<data_block_type*> block_list_type;
  typedef HierarchyStorage<DataType> hierarchy_storage_type;
  typedef typename hierarchy_storage_type::data_block_type hdata_block_type;
  typedef typename hierarchy_storage_type::block_list_type hblock_list_type;
  typedef typename hblock_list_type::iterator hblock_list_iterator;
  typedef typename data_block_type::direction direction;
  typedef std::pair<int, int> node_pair_type;
  typedef std::multimap<int, int> node_pair_list_type;
  typedef std::vector<float> point_list_type;
  typedef Evaluator<VectorType> evaluator_type;
  typedef VisualGridBase vgridbase_type;

  VisualGrid(evaluator_type* ev) : 
    vgridbase_type(ev), _FKeys(1), _MinLevel(0), _MaxLevel(0), _NLevels(1), _MaxProcs(1024),
    StepFormat(0), StepDirectories(0), _CutOut(1), _Eliminate(1), _EliminateValue(-1.e37), _parData(0) {
    PlotLevel = new int[NLevels()];
    int i;
    for (i=0; i<NLevels(); i++)
      PlotLevel[i] = 0;

    VectorType dummy;
    VectorLength = dummy.length();
    CompName = new std::string[VectorLength];
    CompDir = "-";
    CompRead = new bool[VectorLength];
    for (i=0; i<VectorLength; i++) 
      CompRead[i] = false;
    for (i=0; i<3; i++) {
      Symmetry[i] = 0; Periodic[i] = 0; 
      show[2*i] = 0; show[2*i+1] = -1.0; cut[i] = FLT_MAX;
      show_exact[i] = 0;
    }
    DimUsed = Coords(3,1);
    Make3d = false;
    for (i=0; i<DAGHMaxLevels; i++)
      _RefineFactor[i] = DAGHDefaultRefineFactor;
  }

  ~VisualGrid() {
    delete PlotLevel;
    delete [] CompName;
    delete [] CompRead;
  }
  
  //******************************************************************************
  // Abstract class interface
  //******************************************************************************
  virtual void BlockListAppend(hdata_block_type& workdata, BBox& bb, int& comp, int& Level,
                               int& proc, float* dx, float* geom) = 0;
  //******************************************************************************

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    ControlDevice GHCtrl = Ctrl.getSubDevice("GridHierarchy");
    char VariableName[32];
    for (int d=0; d<3; d++) {
      std::sprintf(VariableName,"Cells(%d)",d+1);
      RegisterAt(GHCtrl,VariableName,shape[d]);
      std::sprintf(VariableName,"GeomMin(%d)",d+1);
      RegisterAt(GHCtrl,VariableName,geom[2*d]);
      std::sprintf(VariableName,"GeomMax(%d)",d+1);
      RegisterAt(GHCtrl,VariableName,geom[2*d+1]);     
    }   
    RegisterAt(GHCtrl,"MaxLevels",_NLevels);
    for (register int i=0; i<DAGHMaxLevels; i++) {
      std::sprintf(VariableName,"RefineFactor(%d)",i+1);
      RegisterAt(GHCtrl,VariableName,_RefineFactor[i]);
    }

    ControlDevice FilesCtrl = Ctrl.getSubDevice("Files");
    for (int comp=0; comp<VectorLength; comp++) {
      char ComponentName[16];
      std::sprintf(ComponentName,"-");
      CompName[comp] = ComponentName; 
      std::sprintf(ComponentName,"OutputName(%d)",comp+1);
      RegisterAt(FilesCtrl,ComponentName,CompName[comp]);
    }    
    RegisterAt(FilesCtrl,"OutputDirectory",CompDir);
    RegisterAt(FilesCtrl,"StepFormat",StepFormat);
    RegisterAt(FilesCtrl,"StepDirectories",StepDirectories);    
  }
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, ""); }
  
  virtual void update() {    
    std::string iname;
    if (UpdateName()) iname = *UpdateName();
    else iname = "display.in";
    ControlDevice DisplayCtrl(GetFileControlDevice(iname.c_str(),""));
    RegisterAt(DisplayCtrl,"Type",_FKeys); 
    RegisterAt(DisplayCtrl,"DisplayMinLevel",_MinLevel); 
    RegisterAt(DisplayCtrl,"DisplayMaxLevel",_MaxLevel);
    RegisterAt(DisplayCtrl,"CutOut",_Eliminate);
    RegisterAt(DisplayCtrl,"CutOutValue",_EliminateValue);
    RegisterAt(DisplayCtrl,"MaxProcs",_MaxProcs);
    delete PlotLevel;
    PlotLevel = new int[NLevels()];
    for (int lev=0; lev<NLevels(); lev++) {
      char LevelName[16];      
      std::sprintf(LevelName,"PlotGrid(%d)",lev);
      RegisterAt(DisplayCtrl,LevelName,PlotLevel[lev]);
    }    
    for (int d=0; d<3; d++) {
      char VariableName[16];
      std::sprintf(VariableName,"ShowMin(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,show[2*d]);
      std::sprintf(VariableName,"ShowMax(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,show[2*d+1]);     
      std::sprintf(VariableName,"ShowCut(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,cut[d]);     
      std::sprintf(VariableName,"Symmetry(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,Symmetry[d]);     
      std::sprintf(VariableName,"Periodic(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,Periodic[d]);     
    }   
    TheEvaluator().register_at(DisplayCtrl,"");
    DisplayCtrl.update();  

    TheEvaluator().update();
  }

  virtual int SetUp(int step, std::string** filenames, int Nfilenames) {
    vgridbase_type::SetUp(step,filenames,Nfilenames);
    bool initialize = true;

    if (initialize) {
      for (register int i=0; i<NLevels()-1; i++)
        if (_RefineFactor[i]<DAGHDefaultRefineFactor)
          _RefineFactor[i] = DAGHDefaultRefineFactor;
      _RefineFactor[NLevels()-1] = 1;
      int d, ii;
      for (ii=1; ii<NLevels(); ii++)
        for (d=0; d<3; d++) 
          shape[d] *= _RefineFactor[ii-1];
      
      float dm = 1.e37;
      for (d=0; d<3; d++) {
        dx[d] = (geom[2*d+1]-geom[2*d]) / shape[d]; 
        if (dx[d]>0. && dx[d]<dm) dm = dx[d];
      } 
      for (d=0; d<3; d++) {
        if (dx[d] <= 0.0) {
          DimUsed(d) = 0;
          dx[d] = dm;
          shape[d] = 1;
          for (int comp=VectorLength-1; comp>d; comp--) 
            CompName[comp] = CompName[comp-1];
          if (VectorLength >= 4)
            CompName[d+1] = "-";
        }
      }

      for (d=0; d<3; d++) { 
        if (Symmetry[d] != 0) {
          shape[d] += shape[d];
          if (Symmetry[d] > 0)
            geom[2*d+1] += geom[2*d+1]-geom[2*d];
          else
            geom[2*d] -= geom[2*d+1]-geom[2*d];
        }
        else if (Periodic[d] > 0) {
          shape[d] += Periodic[d]*shape[d];
          geom[2*d+1] += Periodic[d]*(geom[2*d+1]-geom[2*d]);
        }
      }      

      for (d=0; d<3; d++) {
        if (DimUsed(d) && cut[d]>=geom[2*d] && cut[d]<=geom[2*d+1]) {
          show[2*d]   = cut[d];
          show[2*d+1] = cut[d];
          show_exact[d] = 1;
        }
        if (DimUsed(d)==0 || show[2*d]>show[2*d+1]) {
          show_shape[2*d]   = 0;
          show_shape[2*d+1] = shape[d];
        }
        else {
          show_shape[2*d]   = int((Max(show[2*d],geom[2*d])-
                                   geom[2*d]) / dx[d]);
          show_shape[2*d+1] = int((Min(show[2*d+1],geom[2*d+1])-
                                   geom[2*d]+dx[d]/2.0) / dx[d]);
        }
      }
    }

    int i;
    for (i=0; i<Nfilenames; i++) 
      CompName[i] = *filenames[i];
    for (i=0; i<VectorLength; i++) {
      if (CompName[i].c_str()[0] == '-')
        continue;

      char ioname[DAGHBktGFNameWidth+16];
      OutputFileName(ioname,CompName[i].c_str(),step);
      hierarchy_storage_type* HS = new hierarchy_storage_type(NLevels());
      _parData=HS->ReadIn(ioname, Symmetry, Periodic, _Eliminate, 
                          _EliminateValue, MaxProcs());
      if (_parData >= 0) {
        CompRead[i] = true;
        if (HS->RealTime() >= _RealTime) _RealTime = HS->RealTime();
        std::cerr << "Reading Component " << i << " : " << CompName[i]
                  << " at " << RealTime() << std::endl;
      }

      int* NCells = new int[MaxProcs()*NLevels()];
      int pc, Level;
      for (pc=0; pc<MaxProcs()*NLevels(); pc++) NCells[pc] = 0;
      for (Level=0; Level<NLevels(); Level++) 
        for (hblock_list_iterator hsbit=HS->BlockList(Level).begin();
             hsbit!=HS->BlockList(Level).end(); hsbit++) 
          NCells[MaxProcs()*Level+(*hsbit).second] +=
            (*hsbit).first->bbox().size();
      int* WCells = new int[MaxProcs()]; 
      int* Work = new int[MaxProcs()]; 
      int WholeCells = 0, WholeWork = 0, NProcs = 0;
      for (pc=0; pc<MaxProcs(); pc++) {
        WCells[pc] = NCells[pc]; 
        Work[pc] = NCells[pc];
        int fac = 1;
        for (Level=1; Level<NLevels(); Level++) {
          fac *= RefineFactor(Level-1);
          WCells[pc] += NCells[MaxProcs()*Level+pc];
          Work[pc] += fac*NCells[MaxProcs()*Level+pc];
        }
        WholeCells += WCells[pc];       
        WholeWork += Work[pc];
        if (Work[pc] != 0) 
          if (pc > NProcs) NProcs = pc;
      }
      for (pc=0; pc<MaxProcs(); pc++) 
        if (WCells[pc] != 0) 
          std::cerr << "Processor " << pc << ": " << WCells[pc]
                    << " Cells, Load: "
                    << float(Work[pc]) / (float(WholeWork) / float(NProcs+1))
                    << std::endl;         
      std::cerr << "On all Procs: " << WholeCells  << " Cells" << std::endl;
      delete [] NCells;
      delete [] WCells;
      delete [] Work;
 
      if (initialize) {
        if (_MinLevel < 0) 
          _MinLevel = 0;
        if (_MinLevel >=  HS->NLevels()) 
          _MinLevel = HS->NLevels()-1;
        if (_MaxLevel < 0) 
          _MaxLevel = 0;
        if (_MaxLevel >=  HS->NLevels()) 
          _MaxLevel = HS->NLevels()-1;
      }

      typename block_list_type::iterator bit=BlockList().begin();
      for (Level=MaxLevel(); Level>=MinLevel(); Level--) {
        int sh[3];
        int d;
        for (d=0; d<3; d++) {
          sh[d] = show_shape[2*d+1];
          if (sh[d]>show_shape[2*d]) sh[d] -= 1;
        }
        globalbb = BBox(3,show_shape[0],show_shape[2],show_shape[4],
                        sh[0],sh[1],sh[2],1);   

        for (d=0; d<3; d++) {
          int maxrf=1;
          for (int ii=(show_exact[d] ? Level : 0); ii<NLevels()-1; ii++) 
            maxrf *= _RefineFactor[ii];   
          globalbb.coarsen(d,maxrf);
          globalbb.refine(d,maxrf,0);
        }

        int cnt = 0;
        int cells = 0;
        for (hblock_list_iterator hsbit=HS->BlockList(Level).begin();
             hsbit!=HS->BlockList(Level).end(); hsbit++, cnt++) {
          if (initialize) {
            cells += (*hsbit).first->bbox().size();
            BBox bb((*hsbit).first->bbox());
            bb *= globalbb;
            bb.lower() = (bb.lower()/bb.stepsize())*bb.stepsize();
            bb.upper() = (bb.upper()/bb.stepsize())*bb.stepsize();
            BBoxList bbl;
            bbl.add(bb);
            if (_CutOut) {
              BBoxList bblall;
              for (typename block_list_type::iterator bitall=BlockList().begin();
                   bitall!=BlockList().end(); bitall++)
                bblall.add(coarsen((**bitall).DB().bbox(),
                                   (*hsbit).first->bbox().stepsize()/
                                   (**bitall).DB().bbox().stepsize()));
              bbl -= bblall;
            }
            if (!bbl.isempty()) 
              for (BBox* bbp = bbl.first(); bbp; bbp=bbl.next()) 
                if (!bbp->empty()) 
                  BlockListAppend(*((*hsbit).first), *bbp, i, Level, (*hsbit).second, dx, geom);
          }
          else {
            if (bit!=BlockList().end()) {
              bool inside;
              do {
                BBox bb = (*hsbit).first->bbox()*(**bit).DB().bbox();
                if (!bb.empty() &&
                    bb.stepsize()==(**bit).DB().bbox().stepsize()) {
                  BeginFastIndex3(source, (*hsbit).first->bbox(), (*hsbit).first->data(),
                                  DataType);
                  BeginFastIndex3(target, (**bit).DB().bbox(),
                                (**bit).DB().data(), VectorType);
                  for_3 (n, m, l, bb, bb.stepsize())
                    FastIndex3(target,n,m,l)(i) = FastIndex3(source,n,m,l); 
                  end_for
                  EndFastIndex3(source);
                  EndFastIndex3(target);
                  inside=true;
                  bit++;
                }
                else inside=false;
              } while (inside && bit!=BlockList().end());
            }
          }
        }
        if (initialize) 
          std::cerr << "On Level " << Level << ": " << cnt << " Grids with " 
                    << cells << " Cells read in." << std::endl;
      }
      
      initialize = false;
      delete HS;
    }
    
    if (initialize)
      return -1;

    for (i=0; i<VectorLength; i++) {
      if (CompName[i].c_str()[0] == '-') {
        std::cerr << "Setting Component " << i << " to zero." << std::endl;
        for (typename block_list_type::iterator bit=BlockList().begin(); 
             bit!=BlockList().end(); bit++) {
          BeginFastIndex3(target, (**bit).DB().bbox(), 
                          (**bit).DB().data(), VectorType);
          for_3 (n, m, l, (**bit).DB().bbox(), (**bit).DB().bbox().stepsize())
            FastIndex3(target,n,m,l)(i) = DataType(0.0); 
          end_for
          EndFastIndex3(target);
        }
      } 
    }
    return 0;
  }
  
  void OutputFileTime(char *time, const int& Time) {
    if (StepFormat<=0)
      std::sprintf(time,"%d",Time); 
    else if (StepFormat==1) {
      if (Time<=999999) 
        std::sprintf(time,"%6.6d",Time);
      else if (Time<=99999999) 
        std::sprintf(time,"%8.8d",Time);
      else
        std::sprintf(time,"%d",Time);
    }
  }    

  void OutputFileName(char *ioname, const char* name, const int& Time) {
    char time[20], help[DAGHBktGFNameWidth+256];
    OutputFileTime(time,Time);
    if (CompDir.c_str()[0] != '-' && CompDir.length()>0) 
      std::sprintf(ioname,"%s/",CompDir.c_str());
    else
      ioname[0]='\0';
    if (StepDirectories>0) {
      std::sprintf(help,"%s/",time); 
      strcat(ioname,help);
    }
    std::sprintf(help,"%s_%s",name,time); 
    strcat(ioname,help);
  }

  virtual int NCells() {
    int cells=0; 
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++)
      cells += (*bit)->NCells();
    return cells;
  }  
  virtual int NNodes() {
    int nodes=0; 
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++)
      nodes += (*bit)->NNodes();
    return nodes;
  }
  
  virtual void SetBlocks(int blocks[]) {
    int offset=0; int noffset=0;
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++)
      (*bit)->SetBlocks(blocks,offset,noffset);
  }

  virtual void SetGridCells(float xyz[][3]) {
    int offset=0;
    std::cerr << "SetGridCells()";
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++) 
      (*bit)->SetGridCells(xyz,offset);
    std::cerr << std::endl;
  }
  virtual void SetGridNodes(float xyz[][3]) {
    int offset=0;
    std::cerr << "SetGridNodes()";
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++) 
      (*bit)->SetGridNodes(xyz,offset);
    std::cerr << std::endl;
  }

  virtual void SetGridConns(int con[][8]) {
    int offset=0;
    std::cerr << "SetGridConns()";
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++) 
      (*bit)->SetGridConns(con,offset);
    std::cerr << std::endl;
  }

  virtual void SetScalCells(int *key,float v[]) {
    int offset=0;
    std::cerr << "SetScalCells()";
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++) 
      (*bit)->SetScalCells(key,v,offset,CompRead);
    std::cerr << std::endl;
  }
  virtual void SetScalNodes(int *key,float v[]) {
    int offset=0;
    std::cerr << "SetScalNodes()";
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++) 
      (*bit)->SetScalNodes(key,v,offset,CompRead);
    std::cerr << std::endl;
  }
  
  virtual void SetPointsCells(int *key,point_list_type &p,float v[]) {
    std::cerr << "SetPointsCells()";
    point_list_type p_tmp(p);
    for (register unsigned int i=0;i<p.size()/3;i++) v[i]=-1.e37;
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++) 
      (*bit)->SetPointsCells(key,p_tmp,v,CompRead);
    std::cerr << std::endl;
  }
  virtual void SetPointsNodes(int *key,point_list_type &p,float v[]) {
    std::cerr << "SetPointsNodes()";
    point_list_type p_tmp(p);
    for (register unsigned int i=0;i<p.size()/3;i++) v[i]=-1.e37;
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++) 
      (*bit)->SetPointsNodes(key,p_tmp,v,CompRead);
    std::cerr << std::endl;
  }
  
  virtual void VectNodes(int *key,float v[][3]) {
    int offset=0;
    std::cerr << "VectNodes()";
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++) 
      (*bit)->VectNodes(key,v,offset,CompRead);
    std::cerr << std::endl;
  } 
  virtual void VectCells(int *key,float v[][3]) {
    int offset=0;
    std::cerr << "VectCells()";
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++) 
      (*bit)->VectCells(key,v,offset,CompRead);
    std::cerr << std::endl;
  } 
      
  virtual int Size() const { return BlockList().size(); }
  virtual const int& MinLevel() const { return _MinLevel; }
  virtual const int& MaxLevel() const { return _MaxLevel; }
  virtual const int& NLevels() const { return _NLevels; }
  virtual const int& MaxProcs() const { return _MaxProcs; }
  virtual const int& RefineFactor(const int l) const
    { return _RefineFactor[l]; }
  virtual const BBox& GlobalBBox() const { return globalbb; }
  virtual const int& FKeys() const { return _FKeys; }

  virtual evaluator_type& TheEvaluator() const 
    { return (evaluator_type &) vgridbase_type::TheEvaluator(); }
  virtual evaluator_type& TheEvaluator() 
    { return (evaluator_type &) vgridbase_type::TheEvaluator(); }

  block_list_type& BlockList() { return *_BlockList; }
  block_list_type& BlockList() const { return *_BlockList; }

  const Coords Dimensions() const { return ( Make3d ? Coords(3,1): DimUsed) ; }

protected:
  block_list_type* _BlockList;
  int _FKeys;
  int _MinLevel;
  int _MaxLevel;
  int _NLevels;
  int _MaxProcs;
  int _RefineFactor[DAGHMaxLevels];
  int* PlotLevel;
  int shape[3];
  BBox globalbb;
  float geom[2*3];
  float show[2*3];  
  float cut[3];
  int show_shape[2*3];
  int show_exact[3];
  float dx[3];
public:
  Coords DimUsed;
  bool Make3d;
protected:
  int VectorLength;
  std::string* CompName;
  std::string CompDir;
  int StepFormat, StepDirectories;
  bool* CompRead;
  int Symmetry[3];
  int Periodic[3];
  int _CutOut;
  int _Eliminate;
  DataType _EliminateValue;
  int _parData;
};

#endif

---