FileVisualGrid.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
// Santiago Lombeyda (VTU, SILO)

#ifndef AMROC_FILEVISUALGRID_H
#define AMROC_FILEVISUALGRID_H

#include "VisualGrid.h" 
#include "HDFToFile/FileVisualGridBase.h"
#include "EncodeToBase64.h"
#include <strings.h>
#include <cctype>
#include <hlimits.h>

#if defined(HAVE_LIBSILO) || defined(HAVE_LIBSILOH5)
#include <silo.h>
#endif

#define PROBE_TOLERANCE 1.e-6

typedef float points_list_type[][3];
typedef int   cons_list_type[][8];

template <class VectorType>
class FileVisualGrid : public FileVisualGridBase, 
                       public VisualGrid<VectorType> {
  typedef typename VectorType::InternalDataType DataType;
  typedef VisualGrid<VectorType> base;
public:
  typedef typename base::data_block_type data_block_type;
  typedef typename base::block_list_type block_list_type;
  typedef typename base::hierarchy_storage_type hierarchy_storage_type;
  typedef typename base::hdata_block_type hdata_block_type;
  typedef typename base::evaluator_type evaluator_type;
  typedef typename base::point_list_type point_list_type;

  typedef std::multimap<float, int> index_ordered_map_type;
  typedef std::pair<float, int> ordered_index_type;

public:
  FileVisualGrid(evaluator_type* ev) : FileVisualGridBase(ev), 
                                       base(ev), OutputType(0) {
    base::_BlockList = (block_list_type *) new block_list_type;
    nkeys = 0;
    keys = (int *) 0;
    OutputName = "-";
    for (int d=0; d<3; d++) {
      LinePoint[d]  = 0.0;   
      LineVector[d] = 0.0;   
      RegData[d] = 10;
    } 
    RegLim[0] = FLT_MIN;
    RegLim[1] = FLT_MAX;
    RegMax = 1;
    LevelCutOut = 0;
    PointOutput = 0;
    _Step = -1;
    FileMerge = 1;
    PointFile = "";
  }

  ~FileVisualGrid() {
    if (keys) delete [] keys;
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++)
      delete (*bit);
    BlockList().clear();      
    delete base::_BlockList;
  }

  virtual void update() {    
    std::string iname;
    if (UpdateName()) iname = *UpdateName();
    else iname = "display_file.in";

    // Generating display_file.in if file is missing; 
    std::ifstream infile;
    std::ofstream outfile;
    infile.open(iname.c_str(), std::ios::in);
    if (!infile) {
      if (!UpdateName()) {
        outfile.open(iname.c_str(), std::ios::out);
        outfile << "Type  6\n\nDisplayMinLevel    0\nDisplayMaxLevel    10\n\n";
        outfile << "FileName convert\nFileType 6\nKeys p,t,s,i\n"; 
        outfile.close();
      }
    }
    else{
      infile.close();
    }

    ControlDevice DisplayCtrl(GetFileControlDevice(iname.c_str(),""));
    RegisterAt(DisplayCtrl,"Type",base::_FKeys); 
    RegisterAt(DisplayCtrl,"FileName",OutputName);     
    RegisterAt(DisplayCtrl,"FileType",OutputType);   
    RegisterAt(DisplayCtrl,"Keys",StrKeys);     
    RegisterAt(DisplayCtrl,"DisplayMinLevel",base::_MinLevel); 
    RegisterAt(DisplayCtrl,"DisplayMaxLevel",base::_MaxLevel);
    RegisterAt(DisplayCtrl,"CutOut",base::_Eliminate);
    RegisterAt(DisplayCtrl,"CutOutValue",base::_EliminateValue);
    RegisterAt(DisplayCtrl,"RegularMax",RegMax);
    RegisterAt(DisplayCtrl,"RegularLimit(1)",RegLim[0]);
    RegisterAt(DisplayCtrl,"RegularLimit(2)",RegLim[1]);
    RegisterAt(DisplayCtrl,"PointOutput",PointOutput);
    RegisterAt(DisplayCtrl,"LevelCutOut",LevelCutOut);
    RegisterAt(DisplayCtrl,"FileMerge",FileMerge);
    RegisterAt(DisplayCtrl,"PointFile",PointFile);
    for (int d=0; d<3; d++) {
      char VariableName[16];
      std::sprintf(VariableName,"ShowMin(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,base::show[2*d]);
      std::sprintf(VariableName,"ShowMax(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,base::show[2*d+1]);     
      std::sprintf(VariableName,"ShowCut(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,base::cut[d]);     
      std::sprintf(VariableName,"Symmetry(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,base::Symmetry[d]);     
      std::sprintf(VariableName,"Periodic(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,base::Periodic[d]);     
      std::sprintf(VariableName,"LinePoint(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,LinePoint[d]);     
      std::sprintf(VariableName,"LineVector(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,LineVector[d]);          
      std::sprintf(VariableName,"RegularData(%d)",d+1);
      RegisterAt(DisplayCtrl,VariableName,RegData[d]);
    }   
    TheEvaluator().register_at(DisplayCtrl,"");
    DisplayCtrl.update();  

    TheEvaluator().update();
  }

  virtual int SetUp(int step, std::string** filenames, int Nfilenames) { 
    int res;
    if ((OutputType==0 || OutputType==1) && !PointFile.empty())
      base::_CutOut = 0;
    if (OutputType==6 || OutputType==7) { 
      base::_CutOut = 0;
      _Step = step;
    }
    if (OutputType==15)
      if (LevelCutOut<=0 && FKeys()!=1)
        base::_CutOut = 0;
      
    if ((res=base::SetUp(step,filenames,Nfilenames)) >= 0) {
      char append_str[32], time[20];
      base::OutputFileTime(time, step);
      if (OutputName != "-") {
        if (OutputType>=0 && OutputType<=2) 
          std::sprintf(append_str,"_%s.txt",time);
        else if (OutputType>=3 && OutputType<=5) 
          std::sprintf(append_str,"_%s.dx",time);
        else if (OutputType==8) 
          std::sprintf(append_str,"_%s.dat",time);
        else if (OutputType==10 || OutputType==11) 
          std::sprintf(append_str,"_%s.vtk",time);
        else if (OutputType==12 || OutputType==13) 
          std::sprintf(append_str,"_%s.vtu",time);
        else if (OutputType==15) 
          std::sprintf(append_str,"_%s.silo",time);
        if (OutputType<=5 || OutputType>=8) OutputName += append_str;
      }

      int* work = new int[StrKeys.length()];
      for (unsigned int sk=0; sk<StrKeys.length(); sk++) 
        for (int i=0; i<TheEvaluator().NKeys(); i++) 
          if (TheEvaluator().KeyboardKeys()[i] == StrKeys.c_str()[sk]) {
            work[nkeys] = i+1; nkeys++; 
          }
      
      if (keys) delete [] keys;
      keys = new int[nkeys];
      for (int k=0; k<nkeys; k++)
        keys[k] = work[k];

      delete work;
    }
    return res;
  }


  virtual void Write() {
    if (OutputType>=6 && OutputType<=7 && OutputName=="-")
      OutputName = "out";
    if (OutputName != "-") {
      if (OutputType <= 5 || ( OutputType >= 8 && OutputType < 15 )) {
        std::ofstream outfile;
        outfile.open(OutputName.c_str(), std::ios::out);
        std::ostream* out = new std::ostream(outfile.rdbuf());
        if (OutputType == 0) 
          WriteASCIITabular(0,0,*out);
        else if (OutputType == 1) 
          WriteASCIITabular(1,0,*out); 
        else if (OutputType == 2) 
          WriteRegularTabular(*out);
        else if (OutputType == 3) 
          WriteDXFile(*out);
        else if (OutputType == 4) 
          WriteDXASCIIFile(*out);
        else if (OutputType == 5) 
          WriteDXExternalFilter(*out);
        else if (OutputType == 8) 
          WriteTecplotASCII(*out);
        else if (OutputType == 10) 
          WriteVTKASCIIFile(*out);
        else if (OutputType == 11) 
          WriteVTKFile(*out);
        else if (OutputType == 12)
          WriteVTUASCIIFile(*out);
        else if (OutputType == 13)
          WriteVTUFile(*out);
        outfile.close();
        delete out;
      }
      else if (OutputType == 6) 
        WriteHDFConvert(0);
      else if (OutputType == 7) 
        WriteHDFConvert(1);
      else if (OutputType == 15)
        WriteSiloFile((char *)OutputName.c_str());
    }
    else {
      if (OutputType == 0) 
        WriteASCIITabular(0,0,std::cout);
      else if (OutputType == 1) 
        WriteASCIITabular(1,0,std::cout); 
      else if (OutputType == 2) 
        WriteRegularTabular(std::cout);
      else if (OutputType == 3) 
        WriteDXFile(std::cout);
      else if (OutputType == 4) 
        WriteDXASCIIFile(std::cout);
      else if (OutputType == 5) 
        WriteDXExternalFilter(std::cout);
      else if (OutputType == 8) 
        WriteTecplotASCII(std::cout);
      else if (OutputType == 10)
        WriteVTKASCIIFile(std::cout);
      else if (OutputType == 11)
        WriteVTKFile(std::cout);
      else if (OutputType == 12) 
        WriteVTUASCIIFile(std::cout);
      else if (OutputType == 13) 
        WriteVTUFile(std::cout);
      else if (OutputType == 15)
        WriteSiloFile(0);
    }
  }

  void WriteRegularTabular(std::ostream& out) {
    if (!keys) return;
    int length;
    if (FKeys() == 1) 
      length = NNodes();
    else
      length = NCells();
    
    float* xyz = new float[3*length];
    SetGrid((float (*)[3]) xyz); 
    float* dat = new float[nkeys*length];
    float* val = new float[nkeys];
    for (int i=0; i<nkeys*length; i++) *(dat+i) = 0.;
    for (int k=0; k<nkeys; k++) 
      SetScal(&(keys[k]),&(dat[k*length]));

    float regdx[3], reggeom[6];
    int xd = -1, yd = -1, fd = -1;
    for (int d=0; d<3; d++) {
      if (base::DimUsed(d)==0 || base::show[2*d]>=base::show[2*d+1]) {
        reggeom[2*d] = base::geom[2*d]; reggeom[2*d+1] = base::geom[2*d+1]; 
      }
      else {
        reggeom[2*d] = Max(base::show[2*d],base::geom[2*d]); 
        reggeom[2*d+1] = Min(base::show[2*d+1],base::geom[2*d+1]); 
      }
      if (RegData[d]>0) regdx[d] = (reggeom[2*d+1] - reggeom[2*d])/RegData[d];
      else regdx[d] = 0.0;
      if (yd<0 && xd>=0 && regdx[d]>0.0) yd = d;
      if (xd<0 && regdx[d]>0.0) xd = d;
      if (fd<0 && regdx[d]==0.0) fd = d;
    }
    float xs = reggeom[2*xd] + (FKeys()==1 ? 0.0 : regdx[xd]/2.0);
    while (xs<=reggeom[2*xd+1] && regdx[xd]>0.0) {
      float ys = reggeom[2*yd] + (FKeys()==1 ? 0.0 : regdx[yd]/2.0);
      while (ys<=reggeom[2*yd+1] && regdx[yd]>0.0) {
        float mindist=FLT_MAX;
        float fddist=(RegMax ? FLT_MIN : FLT_MAX);
        int inearest = -1;
        for (register int i=0; i<length; i++) {
          float sum=std::sqrt((xs-xyz[3*i+xd])*(xs-xyz[3*i+xd])+
                              (ys-xyz[3*i+yd])*(ys-xyz[3*i+yd]));
          if (sum<mindist-PROBE_TOLERANCE*(mindist>0?mindist:1.0) ||
              inearest<0) {
            mindist = sum;
            inearest = i;
          }
          if (fd>=0 && sum<mindist+PROBE_TOLERANCE*(mindist>0?mindist:1.0)) {
            if (RegMax && xyz[3*i+fd]>fddist &&
                dat[i]>=RegLim[0] && dat[i]<=RegLim[1]) {
              fddist = xyz[3*i+fd]; inearest = i;
            }
            if (!RegMax && xyz[3*i+fd]<fddist &&
                dat[i]>=RegLim[0] && dat[i]<=RegLim[1]) {
              fddist = xyz[3*i+fd]; inearest = i;
            }
          }
        }
        out << xs << " " << ys;
        for (register int k=0; k<nkeys; k++) 
          out << " " << dat[k*length+inearest];
        out << std::endl;
        std::cerr << ".";
        ys += regdx[yd];
      }
      out << std::endl;
      std::cerr << std::endl;
      xs += regdx[xd];
    }

    delete [] val;
    delete [] dat;
    delete [] xyz;
  }

  void WriteASCIITabular(int ordered, int tecplot, std::ostream& out) {
    if (!keys) return;

    if (PointFile.empty()) {
      int length;
      if (FKeys() == 1)
        length = NNodes();
      else
        length = NCells();
      
      float* xyz = new float[3*length];
      SetGrid((float (*)[3]) xyz);  
      int nkl = nkeys*length;
      float* dat = new float[nkl];
      for (register int i=0; i<nkl; i++) *(dat+i) = 0.;
      register int k;
      for (k=0; k<nkeys; k++) 
        SetScal(&(keys[k]),&(dat[k*length]));
      
      if (!tecplot) {
        out << RealTime() << ";";
        if (!LineProbe()) {
          if (!ordered) {
            if (base::DimUsed(0)) out << "x ;";
            if (base::DimUsed(1)) out << "y ;";
            if (base::DimUsed(2)) out << "z ;";
          }
          else
            if (base::DimUsed(0)) out << "x ;";
            else if (base::DimUsed(1)) out << "y ;";
            else if (base::DimUsed(2)) out << "z ;";
        }
        else
          out << "r ;";
        for (k=0; k<nkeys; k++) 
          out << TheEvaluator().TitleKeys(keys[k]-1) << ";";
      }
      else {
        out << "title = \"t = " << RealTime() << "\"" << std::endl;
        out << "variables = ";
        if (!LineProbe()) out << "\"x\"";
        else out << "\"r\"";
        for (k=0; k<nkeys; k++) 
          out << ", \"" << TheEvaluator().TitleKeys(keys[k]-1) << "\"";
      }
      out << std::endl;
   
      if (!ordered) {
        for (register int i=0; i<length; i++) {
          if (!LineProbe()) {
            for (register int d=0; d<3; d++) 
              if (base::DimUsed(d))
                out << xyz[3*i+d] << " ";        
            for (k=0; k<nkeys; k++) 
                out << dat[k*length+i] << " ";
            out << std::endl;
          }
          else {
            int lp = 1, d;
            float di = 0.0, xs = 0.0, c = 0.0;
            for (d=0; d<3; d++) 
              if (base::DimUsed(d)) {
                di += (xyz[3*i+d]-LinePoint[d])*(xyz[3*i+d]-LinePoint[d]);
                if (LineVector[d] != 0.0) {
                  xs += (xyz[3*i+d]-LinePoint[d])/LineVector[d]; c += 1.0;
                }
              }   
            c = xs/c;
            for (d=0; d<3; d++) 
              if (base::DimUsed(d) && 
                  std::fabs(c*LineVector[d]-xyz[3*i+d]+LinePoint[d])>
                  PROBE_TOLERANCE) {
                lp = 0; break;
              }
            if (lp) {
              out << std::sqrt(di) << " ";        
              for (k=0; k<nkeys; k++) 
                out << dat[k*length+i] << " ";
              out << std::endl;
            }
          }
        }
      }
      else {
        index_ordered_map_type imap;
        for (register int i=0; i<length; i++) 
          if (!LineProbe()) {
            for (register int d=0; d<3; d++)
                if (base::DimUsed(d)) {
                  imap.insert(ordered_index_type(xyz[3*i+d],i));
                  break;
                }
          }
          else {
            int lp = 1, d;
            float di = 0.0, xs = 0.0, c = 0.0;
            for (d=0; d<3; d++) 
              if (base::DimUsed(d)) {
                di += (xyz[3*i+d]-LinePoint[d])*(xyz[3*i+d]-LinePoint[d]);
                if (LineVector[d] != 0.0) {
                  xs += (xyz[3*i+d]-LinePoint[d])/LineVector[d]; c += 1.0;
                }
              }
            c = xs/c;
            for (d=0; d<3; d++) 
              if (base::DimUsed(d) && 
                  std::fabs(c*LineVector[d]-xyz[3*i+d]+LinePoint[d])>
                  PROBE_TOLERANCE) {
                lp = 0; break;
              }
            if (lp) 
              imap.insert(ordered_index_type(std::sqrt(di),i));
          }
        
        for (index_ordered_map_type::iterator it = imap.begin();
             it != imap.end(); ++it) {
          out << (*it).first << " ";
          for (k=0; k<nkeys; k++) 
            out << dat[k*length+(*it).second] << " ";
          out << std::endl;
        }
      }
    
      delete [] dat;
      delete [] xyz;
    }    
    else {
      std::ifstream in(PointFile.c_str(), std::ios::in);
      if (!in) {
        std::cerr << "Cannot open " << PointFile << " for reading!" << std::endl;
        return;
      }
      float SP[3];
      int d, k;
      point_list_type p;
      while (!in.eof()) {
        for (d=0; d<3; d++) 
          if (base::DimUsed(d))
            in >> SP[d];
          else SP[d]=0.0;
        if (in.eof()) continue;
        p.insert(p.end(), SP, SP+3);
      }
      in.close();
      if (!p.empty()) {
        unsigned int Np=p.size()/3;
        float *v=new float[Np*nkeys];
        for (k=0; k<nkeys; k++) 
          if (FKeys() == 1)
            base::SetPointsNodes(&(keys[k]),p,&(v[k*Np]));
          else
            base::SetPointsCells(&(keys[k]),p,&(v[k*Np]));
        
        for (register unsigned int i=0;i<Np;i++) {
          bool print=true;
          for (k=0; k<nkeys; k++) 
            if (v[k*Np+i]<=float(-1.e37)) { print=false; break; }
          if (print) {
            for (d=0; d<3; d++) 
              if (base::DimUsed(d))
                out << p[3*i+d] << " ";
            for (k=0; k<nkeys; k++) 
              out << v[k*Np+i] << " ";
            out << std::endl;
          }
          else {
            std::cerr << "Point ( ";
            for (d=0; d<3; d++) 
              if (base::DimUsed(d))
                std::cerr << p[3*i+d] << " ";
            std::cerr << ") outside of domain. Skipping..." << std::endl; 
          }
        }
        delete [] v;
      }
    }
  }
  
  void DimShape(int& dimshape, int& conshape) {
    if (base::Make3d) {
      dimshape = 3;
      conshape = 8;
    }
    else {
      dimshape = 0;
      conshape = 1;
      for (int d=0; d<3; d++) {
        dimshape += base::DimUsed(d);
        conshape *= (base::DimUsed(d) ? 2 : 1);
      }
    }
  }

  //  
  // All nodes are dumped out. Equal nodes in different subgrids are NOT deleted, but
  // DX does not seem to be disturbed.
  //
  void WriteDXFile(std::ostream& out) {
    int dxfile_dataoffset=0;
    int dimshape, conshape;
    DimShape(dimshape,conshape);

    int nnodes = NNodes(), ncells = NCells();    
    float* xyz = new float[3*nnodes];
    out << "object 1 class array type float rank 1 shape " << dimshape;
    out << " items " << nnodes;
#if defined DX_LSB
    out << " lsb";
#else
    out << " msb";
#endif   
    out << " binary data " << dxfile_dataoffset << std::endl;
    SetGridNodes((float (*)[3]) xyz);
    dxfile_dataoffset+=dimshape*nnodes*sizeof(float);

    int* con = new int[8*ncells];    
    out << "object 2 class array type int rank 1 shape " << conshape;
    out << " items " << ncells;
#if defined DX_LSB
    out << " lsb";
#else
    out << " msb";
#endif   
    out << " binary data " << dxfile_dataoffset << std::endl;
    SetGridConns((int (*)[8]) con);  
    out << "attribute \"element type\" string ";
    if (dimshape==3) out << "\"cubes\"";
    else if (dimshape==2) out << "\"quads\"";
    else if (dimshape==1) out << "\"lines\"";

    out << std::endl << "attribute \"ref\" string \"positions\"" << std::endl;
    dxfile_dataoffset+=conshape*ncells*sizeof(int);

    int length;
    if (FKeys() == 1)
      length = nnodes;
    else
      length = ncells;
    int nkl = nkeys*length;
    float* dat = new float[nkl];
    for (register int i=0; i<nkl; i++) *(dat+i) = 0.;
    out << "object 3 class array type float";
    if (nkeys>1) out << " rank 1 shape " << nkeys;
    out << " items " << length;
#if defined DX_LSB
    out << " lsb";
#else
    out << " msb";
#endif   
    out << " binary data " << dxfile_dataoffset << std::endl;
    register int k;
    for (k=0; k<nkeys; k++) 
      SetScal(&(keys[k]),&(dat[k*length]));
    if (FKeys() == 6)
      out << "attribute \"dep\" string \"connections\"" << std::endl;
    dxfile_dataoffset+= nkeys*length*sizeof(float);

    out << "object \"default\" class field" << std::endl;
    out << "component \"positions\" value 1" << std::endl;
    out << "component \"connections\" value 2" << std::endl;
    out << "component \"data\" value 3" << std::endl;
    out << "end" << std::endl;
    
    register int n,t;
    if (dimshape==1) {
      for (n=0,t=0; n<nnodes; n++)
        xyz[t++] = xyz[3*n];
    }
    else if (dimshape==2) {
      for (n=0,t=0; n<nnodes; n++) {
        xyz[t++] = xyz[3*n]; xyz[t++] = xyz[3*n+1];
      }
    }
    out.write((char*) xyz, dimshape*nnodes*sizeof(float));

    if (conshape==2) {
      for (n=0,t=0; n<ncells; n++) {
        con[t++] = con[8*n]; con[t++] = con[8*n+1];
      }
    }
    if (conshape==4) {
      for (n=0,t=0; n<ncells; n++) {
        con[t++] = con[8*n  ]; con[t++] = con[8*n+1]; 
        con[t++] = con[8*n+2]; con[t++] = con[8*n+3];
      }
    }
    out.write((char*) con, conshape*ncells*sizeof(int));

    if (nkeys>1) {
      float* datout = new float[nkeys*length];
      int nn = 0;
      for (n=0; n<length; n++)
        for (k=0; k<nkeys; k++) {
          datout[nn] = dat[k*length+n]; nn++;
        }
      out.write((char*) datout, sizeof(float)*nkeys*length);
      delete [] datout;
    }
    else 
      out.write((char*) dat, sizeof(float)*length);
      
    delete [] xyz;
    delete [] con;
    delete [] dat;
  }    

  void WriteDXASCIIFile(std::ostream& out) {
    int dimshape, conshape;
    DimShape(dimshape,conshape);
    
    int nnodes = NNodes(), ncells = NCells();
    float* xyz = new float[3*nnodes];
    out << "object 1 class array type float rank 1 shape " << dimshape;
    out << " items " << nnodes << " data follows" << std::endl;
    SetGridNodes((float (*)[3]) xyz);
    register int i;
    for (i=0; i<nnodes; i++) {
      for (int d=0; d<dimshape; d++) 
        out << xyz[3*i+d] << " ";
      out << std::endl;
    }

    int* con = new int[8*ncells];    
    out << "object 2 class array type int rank 1 shape " << conshape;
    out << " items " << ncells << " data follows" << std::endl;
    SetGridConns((int (*)[8]) con);  
    for (i=0; i<ncells; i++) {
      for (int d=0; d<conshape; d++) 
        out << con[8*i+d] << " ";
      out << std::endl;
    }

    out << "attribute \"element type\" string ";
    if (dimshape==3) out << "\"cubes\"";
    else if (dimshape==2) out << "\"quads\"";
    else if (dimshape==1) out << "\"lines\"";

    out << std::endl << "attribute \"ref\" string \"positions\"" << std::endl;
    int length;
    if (FKeys() == 1)
      length = nnodes;
    else
      length = ncells;
    float* dat = new float[nkeys*length];
    for (i=0; i<nkeys*length; i++) *(dat+i) = 0.;
    out << "object 3 class array type float";
    if (nkeys>1) out << " rank 1 shape " << nkeys;
    out << " items " << length << " data follows" << std::endl;
    register int k;
    for (k=0; k<nkeys; k++) 
      SetScal(&(keys[k]),&(dat[k*length]));
    for (i=0; i<length; i++) {
      for (k=0; k<nkeys; k++) 
        out << " " << dat[k*length+i];
      out << std::endl;
    }
    if (FKeys() == 6)
      out << "attribute \"dep\" string \"connections\"" << std::endl;

    out << "object \"default\" class field" << std::endl;
    out << "component \"positions\" value 1" << std::endl;
    out << "component \"data\" value 3" << std::endl;
    out << "component \"connections\" value 2" << std::endl;
    out << "end" << std::endl;
        
    delete [] xyz;
    delete [] con;
    delete [] dat;
  }    

  void WriteDXExternalFilter(std::ostream& out) {
    int dimshape, conshape;
    DimShape(dimshape,conshape);
    
    int nnodes = NNodes(), ncells = NCells();
    float* xyz = new float[3*nnodes];
    out << "object 1 class array type float rank 1 shape " << dimshape;
    out << " items " << nnodes;
#if defined DX_LSB
    out << " lsb";
#else
    out << " msb";
#endif   
    out << " binary data follows" << std::endl;
    SetGridNodes((float (*)[3]) xyz);
    register int n,t;
    if (dimshape==1) {
      for (n=0,t=0; n<nnodes; n++)
        xyz[t++] = xyz[3*n];
    }
    else if (dimshape==2) {
      for (n=0,t=0; n<nnodes; n++) {
        xyz[t++] = xyz[3*n]; xyz[t++] = xyz[3*n+1];
      }
    }
    out.write((char*) xyz, dimshape*nnodes*sizeof(float));
    delete [] xyz;

    int* con = new int[8*ncells];    
    out << "object 2 class array type int rank 1 shape " << conshape;
    out << " items " << ncells;
#if defined DX_LSB
    out << " lsb";
#else
    out << " msb";
#endif   
    out << " binary data follows" << std::endl;
    SetGridConns((int (*)[8]) con);  
    if (conshape==2) {
      for (n=0,t=0; n<ncells; n++) {
        con[t++] = con[8*n]; con[t++] = con[8*n+1];
      }
    }
    if (conshape==4) {
      for (n=0,t=0; n<ncells; n++) {
        con[t++] = con[8*n  ]; con[t++] = con[8*n+1]; 
        con[t++] = con[8*n+2]; con[t++] = con[8*n+3];
      }
    }
    out.write((char*) con, conshape*ncells*sizeof(int));
    delete [] con;

    out << "attribute \"element type\" string ";
    if (dimshape==3) out << "\"cubes\"";
    else if (dimshape==2) out << "\"quads\"";
    else if (dimshape==1) out << "\"lines\"";

    out << "attribute \"ref\" string \"positions\"" << std::endl;

    int length;
    if (FKeys() == 1)
      length = nnodes;
    else
      length = ncells;
    float* dat = new float[nkeys*length];
    for (int i=0; i<nkeys*length; i++) *(dat+i) = 0.;
    out << "object 3 class array type float ";
    if (nkeys>1) out << "rank 1 shape " << nkeys << " ";
    if (FKeys() == 6)
      out << "attribute \"dep\" string \"connections\"" << std::endl;
    out << "items " << length;
#if defined DX_LSB
    out << " lsb";
#else
    out << " msb";
#endif   
    out << " binary data follows" << std::endl;
    register int k;
    for (k=0; k<nkeys; k++) 
      SetScal(&(keys[k]),&(dat[k*length]));
    if (nkeys>1) {
      float* datout = new float[nkeys*length];
      int nn = 0;
      for (n=0; n<length; n++)
        for (k=0; k<nkeys; k++) {
          datout[nn] = dat[k*length+n]; nn++;
        }
      out.write((char*) datout, sizeof(float)*nkeys*length);
      delete [] datout;
    }
    else 
      out.write((char*) dat, sizeof(float)*length);      
    delete [] dat;

    out << "object \"default\" class field" << std::endl;
    out << "component \"positions\" value 1" << std::endl;
    out << "component \"connections\" value 2" << std::endl;
    out << "component \"data\" value 3" << std::endl;
    out << "end" << std::endl;  
  }    

  void WriteHDFConvert(int DimReduce) {
    if (!keys || _Step<0 || OutputName == "-") return;
    float64 realtime = RealTime();
    int32 ori[3],dims[3],res[4];
    int32 proc,level,gridID;
    int32 step = _Step;
    int32 rank = 3;
    char target[256];
    char dname[16];
    
    if (BlockList().size()>MAX_NC_VARS && FileMerge>=0) {
      std::cerr << "Maximal block count in HDF4 file exceeded. Setting FileMerge to 0." << std::endl;
      FileMerge = 0;
    }

    for (register int k=0; k<nkeys; k++) {
      std::string TargetName = OutputName;
      dname[0] = TheEvaluator().KeyboardKeys()[keys[k]-1]; dname[1] = 0;
      char time[20];
      base::OutputFileTime(time, _Step);
      if (nkeys>1)
        std::sprintf(target,"%s_%s_%s.hdf",OutputName.c_str(),dname,time);   
      else
        std::sprintf(target,"%s_%s.hdf",OutputName.c_str(),time);   
      
      std::cerr << "SetScalCells()";
      gridID=0;
      int last_proc = -1;
      for (typename block_list_type::reverse_iterator bit=BlockList().rbegin();
           bit!=BlockList().rend(); bit++, gridID++) {

        proc = (**bit).Processor();
        level = (**bit).Level();
        
        if (base::_parData==1 && FileMerge<=0) 
          if (proc!=last_proc) {
            if (last_proc>=0) {
              SDSflush(target);
              SDSclose(target);
            }
            if (nkeys>1)
              std::sprintf(target,"%s_%s_%s.%d.hdf",OutputName.c_str(),dname,time,(int)proc);   
            else
              std::sprintf(target,"%s_%s.%d.hdf",OutputName.c_str(),time,(int)proc); 
            last_proc = proc;
          }

        int offset=0;
        hdata_block_type tmpblock((**bit).DB().bbox());
        float* work = new DataType[(**bit).DB().bbox().size()];

        (**bit).SetScalCells(&(keys[k]),work,offset,base::CompRead);    
        BeginFastIndex3(tmp, tmpblock.bbox(), tmpblock.data(), DataType);
        int idx = 0;
        for_3 (n,m,l,tmpblock.bbox(),tmpblock.bbox().stepsize())
          FastIndex3(tmp,n,m,l) = work[idx];
          idx++;              
        end_for
        EndFastIndex3(tmp); 
        delete [] work;

        register int d;
        for (d=0; d<3; d++) {
          ori[d] = (int32) tmpblock.lower(d);
          if (!base::show_exact[d]) res[d] = (int32) tmpblock.stepsize(d);
          else res[d] = 1;
          dims[d] = (int32) tmpblock.extents(d);
        }
        res[3] = 1;
        if (!base::DimUsed(2)) {
          dims[2] = 1;
          ori[2] = 0;
          res[2] = 1;
          res[3] = 1;
          if (!base::DimUsed(1)) {
            dims[1] = 1;
            ori[1] = 0;
            res[1] = 1;
          }
        }
        if (DimReduce) {
          int dc = 0;
          for (d=0; d<3; d++) 
            if (base::show_exact[d]) {
              for (register int dd=d+1; dd<3; dd++) {
                ori[dd-1] = ori[dd];
                res[dd-1] = res[dd];
                dims[dd-1] = dims[dd]; 
              }
              dc++;
              ori[3-dc] = 0; res[3-dc] = 1; dims[3-dc] = 1; 
            }
        }       

        if (sizeof(DataType) == SIZE_FLOAT32)
          SDSsetNT(target,DFNT_FLOAT32);    
        else if (sizeof(DataType) == SIZE_FLOAT64)
          SDSsetNT(target,DFNT_FLOAT64);    
        else if (sizeof(DataType) == SIZE_FLOAT128)
          SDSsetNT(target,DFNT_FLOAT128); 
        else
          std::cerr << "SDSsetNT() not called!" << std::endl;
        AMRwriteData(target,dname,proc,level,gridID,
                     step,realtime,rank,dims,ori,res,
                     tmpblock.data());
      }
      SDSflush(target);
      SDSclose(target);
      std::cerr << std::endl;
    }
  }    
  
  void WriteTecplotASCII(std::ostream& out) {
    int dimshape, conshape;
    DimShape(dimshape,conshape);

    if (dimshape==1) { 
      WriteASCIITabular(1,1,out); 
      return;
    }
    
    int length;
    if (FKeys() == 1)
      length = NNodes();
    else
      length = NCells();
    out << "title = \"t = " << RealTime() << "\"" << std::endl;
    float* xyz = new float[3*length];
    SetGrid((float (*)[3]) xyz);

    out << "variables = \"x\", \"y\"";
    if (dimshape>2) out << ", \"z\"";
    register int i,k;
    for (k=0; k<nkeys; k++) {
      char* title = new char[LEN_TKEYS];
      std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
      i=std::strlen(title)-1;
      while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }      
      out << ", \"" << title << "\"";
    }
    out << std::endl;
    float* dat = new float[nkeys*length];
    for (i=0; i<nkeys*length; i++) *(dat+i) = 0.;
    for (k=0; k<nkeys; k++) 
      SetScal(&(keys[k]),&(dat[k*length]));

    int cnt=0,cbit=0,blength;
    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++, cbit++) {

      Coords ext = (*bit)->DB().extents();
      if (FKeys() == 1) {
        ext += base::DimUsed; 
        blength = (*bit)->NNodes();
      }
      else
        blength = (*bit)->NCells();

      out << "zone t=\"B-" << cbit << " L-" << (*bit)->Level() << " P-" 
          << (*bit)->Processor() << "\", i=" << ext(0) << ", j=" << ext(1);
      if (dimshape>2) out << ", k=" << ext(2);

      if (PointOutput<=0) {
        out << ", f=block" << std::endl;
        for (register int d=0; d<dimshape; d++)
          for (i=0; i<blength; i++) {
            out << xyz[3*(cnt+i)+d] << " ";
            if ((i+1)%16==0) out << std::endl; 
          }
        out << std::endl << std::endl;
        for (k=0; k<nkeys; k++) {
          for (i=0; i<blength; i++) {
            out << dat[k*length+cnt+i] << " ";
            if ((i+1)%16==0) out << std::endl; 
          }
          out << std::endl << std::endl;
        }
      }
      else {
        out << ", f=point" << std::endl;
        for (i=0; i<blength; i++) {
          for (register int d=0; d<dimshape; d++) 
            out << xyz[3*(cnt+i)+d] << " ";
          for (k=0; k<nkeys; k++) 
            out << " " << dat[k*length+cnt+i];
          out << std::endl;
        }
        out << std::endl;
      }
      cnt += blength;
    }
        
    delete [] xyz;
    delete [] dat;
  }


  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  //                                            WRITE VTU ASCII FILE
  //
  // this function converts the amroc dumped hdf files (in memory)
  // into an xml based VTK unstructured grid format. notice that the
  // levels have been flattened. it is quite inefficient.
  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# define VTK_USED_DIMS 3
  void WriteVTUASCIIFile(std::ostream& out)  {
    // VTI ASCI File writer. based of VTKASCII writer.
    // http://www.vtf.website/~slombey/asci/vtk/
    int dimshape, conshape;
    DimShape(dimshape,conshape);

    out << "<?xml version=\"1.0\"?>" << std::endl;
    out << "<?meta name=\"source\" content=\"HDF data created with AMROC's hdf2file\"?>" 
        << std::endl;
    out << "<?meta name=\"timestamp\" content=\"" << RealTime() << "\"?>" << std::endl;

    // VTU XML HEADER
    int nnodes = NNodes(), ncells = NCells();
    float* xyz = new float[3*nnodes];
    out << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" 
        << std::endl;
    out << "<UnstructuredGrid>" << std::endl;
    out << "<Piece NumberOfPoints=\"" << nnodes << "\" NumberOfCells=\"" << ncells << "\">" 
        << std::endl;

    // flatten points (i.e. remove duplicates) and print
    out << "<Points>" << std::endl;
    out << "<DataArray type=\"Float32\" NumberOfComponents=\""<< VTK_USED_DIMS 
        << "\" format=\"ascii\" Name=\"vertex\">" << std::endl;
    SetGridNodes((float (*)[3]) xyz);
    register int i;
    for (i=0; i<nnodes; i++) {
      for (register int d=0; d<VTK_USED_DIMS; d++) 
        out << xyz[3*i+d] << " ";
      out << std::endl;
    }
    out << "</DataArray>" << std::endl;
    out << "</Points>" << std::endl;

    // cells
    int* con = new int[8*ncells];    
    out << "<Cells>" << std::endl;
    out << "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl;
    // print cells->connectivities
    SetGridConns((int (*)[8]) con);  
    for (i=0; i<ncells; i++) {
      for (register int d=0; d<conshape; d++) 
        out << con[8*i+d] << " ";
      out << std::endl;
    }
    out << "</DataArray>" << std::endl;

    // print cells->element offsets
    out << "<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl;
    for (i=1; i<=ncells; i++)
      out << i*conshape << std::endl;
    out << "</DataArray>" << std::endl;
    
    // print cells->types
    int cell_type = 1;
    if (dimshape==3) cell_type = 11;
    else if (dimshape==2) cell_type = 8;
    else if (dimshape==1) cell_type = 3;

    out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">" << std::endl;
    for (i=0; i<ncells; i++) 
      out << cell_type << std::endl;
    out << "</DataArray>" << std::endl;

    out << "</Cells>" << std::endl;
    
    int length;
    if (FKeys() == 1) {
      length = nnodes;
      out << "<PointData>" << std::endl;
    }
    else {
      length = ncells;
      out << "<CellData>" << std::endl;
    }
    float* dat = new float[nkeys*length];
    for (i=0; i<nkeys*length; i++) *(dat+i) = 0.;
    register int k;
    for (k=0; k<nkeys; k++) 
      SetScal(&(keys[k]),&(dat[k*length]));
    for (k=0; k<nkeys; k++) {
      char* title = new char[LEN_TKEYS];
      // clean up title string
        std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
        i=std::strlen(title)-1;
        while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
        int n=std::strlen(title);
        for (i=0; i<n; i++) if (title[i]==' ') title[i]='_';
      out << "<DataArray type=\"Float32\" Name=\"" << title 
          << "\" format=\"ascii\" NumberOfComponents=\"1\">" << std::endl;
      for (i=0; i<length; i++) 
        out << dat[k*length+i] << std::endl; 
      out << "</DataArray>" << std::endl;
    }

    if (FKeys() == 1)
      out << "</PointData>" << std::endl;
    else 
      out << "</CellData>" << std::endl;
    out << "</Piece>" << std::endl;
    out << "</UnstructuredGrid>" << std::endl;
    out << "</VTKFile>" << std::endl;
 
    delete [] xyz;
    delete [] con;
    delete [] dat;
  }

  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  //                                                     WRITE VTU FILE
  //
  // this function converts the amroc dumped hdf files (in memory)
  // into an xml based VTK unstructured grid format. the data is dumped
  // in binary after being 64base converted. (see mime). an extra integer
  // is added at the beginning of each binary array with the length.
  // + some loops have been flattened for efficiency.
  // + the data is dumped to a buffer (2megs) to make more efficient
  //   writes. the size of this buffer can be altered (hardwired). see
  //   encodeToBase64.h
  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  void WriteVTUFile(std::ostream& out)  {
    // VTI ASCI File writer. based of VTKASCII writer.
    // http://www.vtf.website/~slombey/asci/vtk/
    int dimshape, conshape;
    DimShape(dimshape,conshape);
    unsigned long appendedDataOffset=0;

     // hack or needed? not in vtk specs, but seems necessary
    bool addHeaderSizeAsInt32=true;

    out << "<?xml version=\"1.0\"?>" << std::endl;
    out << "<?meta name=\"source\" content=\"HDF data created with AMROC's hdf2file\"?>" 
        << std::endl;
    out << "<?meta name=\"timestamp\" content=\"" << RealTime() << "\"?>" << std::endl;

    // VTU XML HEADER
    int nnodes = NNodes(), ncells = NCells();
    float* xyz = new float[3*nnodes];
    out << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=";
    if (isBigEndian())
        out << "\"BigEndian\">" << std::endl;
    else
        out << "\"LittleEndian\">" << std::endl;
    out << "<UnstructuredGrid>" << std::endl;
    out << "<Piece NumberOfPoints=\"" << nnodes << "\" NumberOfCells=\"" << ncells << "\">" 
        << std::endl;

    // flatten points (i.e. remove duplicates) and print
    out << "<Points>" << std::endl;
    out << "<DataArray type=";
    if (sizeof(float)==4)
      out << "\"Float32\"";
    else if (sizeof(float)==8)
      out << "\"Float64\"";
    else out << "\"Float" << sizeof(float)*8 <<"\""; // usupported!
    out << " NumberOfComponents=\""<< VTK_USED_DIMS 
        << "\" Name=\"vertex\" format=\"appended\" offset=\"" << appendedDataOffset 
        << "\" />" << std::endl;

    unsigned long deltaDataOffset= VTK_USED_DIMS*nnodes;
    deltaDataOffset*=sizeof(float);
    if (addHeaderSizeAsInt32) deltaDataOffset+=4;
    if (deltaDataOffset%3>0) deltaDataOffset+=(3-(deltaDataOffset%3));
    deltaDataOffset=deltaDataOffset/3*4;
    appendedDataOffset+=deltaDataOffset;

    out << "</Points>" << std::endl;
    out << "<Cells>" << std::endl;
    out << "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" 
        << appendedDataOffset << "\" />" << std::endl;

    deltaDataOffset= conshape*ncells;
    deltaDataOffset*=sizeof(int);
    if (addHeaderSizeAsInt32) deltaDataOffset+=4;
    if (deltaDataOffset%3>0) deltaDataOffset+=(3-(deltaDataOffset%3));
    deltaDataOffset=deltaDataOffset/3*4;
    appendedDataOffset+=deltaDataOffset;

    out << "<DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" 
        << appendedDataOffset << "\" />" << std::endl;
    deltaDataOffset= ncells;
    deltaDataOffset*=sizeof(int);
    if (addHeaderSizeAsInt32) deltaDataOffset+=4;
    if (deltaDataOffset%3>0) deltaDataOffset+=(3-(deltaDataOffset%3));
    deltaDataOffset=deltaDataOffset/3*4;
    appendedDataOffset+=deltaDataOffset;

    out << "<DataArray type=\"UInt8\" Name=\"types\" format=\"appended\" offset=\"" 
        << appendedDataOffset << "\" />" << std::endl;
    deltaDataOffset= ncells;
    deltaDataOffset*=sizeof(char);
    if (addHeaderSizeAsInt32) deltaDataOffset+=4;
    if (deltaDataOffset%3>0) deltaDataOffset+=(3-(deltaDataOffset%3));
    deltaDataOffset=deltaDataOffset/3*4;
    appendedDataOffset+=deltaDataOffset;

    out << "</Cells>" << std::endl;
    
    int length;
    if (FKeys() == 1) {
      length = nnodes;
      out << "<PointData>" << std::endl;
    } else {
      length = ncells;
      out << "<CellData>" << std::endl;
    }

    for (int k=0; k<nkeys; k++) {
      char* title = new char[LEN_TKEYS];
      // clean up title string
        std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
        int i=std::strlen(title)-1;
        while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
        int n=std::strlen(title);
        for (i=0; i<n; i++) if (title[i]==' ') title[i]='_';

      out << "<DataArray type=";
      if (sizeof(float)==4) out << "\"Float32\"";
      else if (sizeof(float)==8) out << "\"Float64\"";
      out << " Name=\"" << title << "\" NumberOfComponents=\"1\" format=\"appended\" offset=\"" 
          << appendedDataOffset << "\" />" << std::endl;

      deltaDataOffset=length;
      deltaDataOffset*=sizeof(float);
      if (addHeaderSizeAsInt32) deltaDataOffset+=4;
      if (deltaDataOffset%3>0) deltaDataOffset+=(3-(deltaDataOffset%3));
      deltaDataOffset=deltaDataOffset/3*4;
      appendedDataOffset+=deltaDataOffset;
    }

    if (FKeys() == 1)
      out << "</PointData>" << std::endl;
    else 
      out << "</CellData>" << std::endl;
    out << "</Piece>" << std::endl;
    out << "</UnstructuredGrid>" << std::endl;

    // APPENDED DATA 

    out << "<AppendedData encoding=\"base64\">" << std::endl;
    out << "_";

    SetGridNodes((float (*)[3]) xyz);

    // encode x,y,z (go in sets of 3, and encode leftovers at the end)
    register int i;
    int dataBlockSize;
    dataBlockSize=nnodes*VTK_USED_DIMS*sizeof(float);
    if (addHeaderSizeAsInt32)
      encodeToBase64AndWrite((unsigned char *)&dataBlockSize,sizeof(int),out);
    encodeToBase64AndWrite((unsigned char *)xyz,dataBlockSize,out,true);

    // cells
    int* con = new int[8*ncells];    

    // print cells->connectivities
    SetGridConns((int (*)[8]) con);  

    // encode conns(go in sets of 3, and encode leftovers at the end)
    dataBlockSize=ncells*conshape*sizeof(int);
    if (addHeaderSizeAsInt32)
      encodeToBase64AndWrite((unsigned char *)&dataBlockSize,sizeof(int),out);
    int deltamaxcon=8*sizeof(int);
    int deltacon   =conshape*sizeof(int);
    unsigned char *srccon=(unsigned char *)con;
    unsigned char *dstcon=(unsigned char *)con;
    if (deltamaxcon!=deltacon && conshape<8) {
      for (i=1; i<ncells; i++) {
         srccon+=deltamaxcon;
         dstcon+=deltacon;
         memcpy(dstcon,srccon,deltacon);
      }
    }
    encodeToBase64AndWrite((unsigned char *)con,dataBlockSize,out,true);

    // print cells->element offsets
    // encode element offsets (go in sets of 3, and encode leftovers at the end)
    dataBlockSize=ncells*sizeof(int);
    if (addHeaderSizeAsInt32)
      encodeToBase64AndWrite((unsigned char *)&dataBlockSize,sizeof(int),out);
    for (i=0; i<ncells; i++) con[i]= (i+1)*conshape;
    encodeToBase64AndWrite((unsigned char *)con,dataBlockSize,out,true);


    // print cells->types
    unsigned char cell_type = 1;
    if (dimshape==3) cell_type = 11;
    else if (dimshape==2) cell_type = 8;
    else if (dimshape==1) cell_type = 3;


    // encode element celltypes (go in sets of 3, and encode leftovers at the end)
    dataBlockSize=ncells;
    if (addHeaderSizeAsInt32) {
      encodeToBase64AndWrite((unsigned char *)&dataBlockSize,sizeof(int),out);
    }
    int four_cell_types;
    ((char *)(&four_cell_types))[0]=cell_type;
    ((char *)(&four_cell_types))[1]=cell_type;
    ((char *)(&four_cell_types))[2]=cell_type;
    ((char *)(&four_cell_types))[3]=cell_type;
    for (i=0; i<(ncells+3)/4; i++) con[i]= four_cell_types;
    encodeToBase64AndWrite((unsigned char *)con,dataBlockSize,out,true);

    float* dat = new float[nkeys*length];
    bzero(dat,nkeys*length*sizeof(float));
    register int k;
    for (k=0; k<nkeys; k++) 
      SetScal(&(keys[k]),&(dat[k*length]));
    for (k=0; k<nkeys; k++) {
      // encode each of the data keys [vars] (go in sets of 3, and encode leftovers at the end)
      dataBlockSize=length*sizeof(float);
      if (addHeaderSizeAsInt32) {
        encodeToBase64AndWrite((unsigned char *)&dataBlockSize,sizeof(int),out);
      }
      encodeToBase64AndWrite((unsigned char *)&(dat[k*length]),dataBlockSize,out,true);
    }
    out <<  std::endl;

    out << "<? meta name=\"datasize\" size=\"" << appendedDataOffset << "\" ?>" << std::endl;

    out << std::endl;
    out << "</AppendedData>" << std::endl;
    out << "</VTKFile>" << std::endl;
 
    delete [] xyz;
    delete [] con;
    delete [] dat;
  }

  void WriteVTKASCIIFile(std::ostream& out)  {
    int dimshape, conshape;
    DimShape(dimshape,conshape);

    int nnodes = NNodes(), ncells = NCells();
    float* xyz = new float[3*nnodes];
    out << "# vtk DataFile Version 3.0" << std::endl;
    out << "HDF data created with AMROC's hdf2file" << std::endl;
    out << "ASCII" << std::endl;
    out << "DATASET UNSTRUCTURED_GRID" << std::endl;
    out << "POINTS " << nnodes << " float" << std::endl;
    SetGridNodes((float (*)[3]) xyz);
    register int i;
    for (i=0; i<nnodes; i++) {
      for (int d=0; d<VTK_USED_DIMS; d++) 
        out << xyz[3*i+d] << " ";
      out << std::endl;
    }

    int* con = new int[8*ncells];    
    out << "CELLS " << ncells << " " << ncells*(conshape+1) << std::endl;
    SetGridConns((int (*)[8]) con);  
    for (i=0; i<ncells; i++) {
      out << conshape << " ";
      for (int d=0; d<conshape; d++) 
        out << con[8*i+d] << " ";
      out << std::endl;
    }

    int cell_type = 1;
    if (dimshape==3) cell_type = 11;
    else if (dimshape==2) cell_type = 8;
    else if (dimshape==1) cell_type = 3;

    out << "CELL_TYPES " << ncells << std::endl;
    for (i=0; i<ncells; i++) 
      out << cell_type << std::endl;
    
    int length;
    if (FKeys() == 1) {
      length = nnodes;
      out << "POINT_DATA ";
    }
    else {
      length = ncells;
      out << "CELL_DATA ";
    }
    out << length << std::endl;    
    float* dat = new float[nkeys*length];
    bzero(dat,nkeys*length*sizeof(float));

    register int k;
    for (k=0; k<nkeys; k++) 
      SetScal(&(keys[k]),&(dat[k*length]));
    out << "FIELD Data " << nkeys << std::endl;     
    for (k=0; k<nkeys; k++) {
      char* title = new char[LEN_TKEYS];
      std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
      i=std::strlen(title)-1;
      while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
      int n=std::strlen(title);
      for (i=0; i<n; i++) 
        if (title[i]==' ') title[i]='_';
      out << title << " 1 " << length << " float" << std::endl;
      for (i=0; i<length; i++) 
        out << dat[k*length+i] << std::endl; 
    }
    
    delete [] xyz;
    delete [] con;
    delete [] dat;
  }   


  void flipBytes(char *c, int n) {
     char *ptr=c;
     char x;
     n/=4;
     for (register int i=0; i<n; i++) {
       x=ptr[3]; ptr[3]=ptr[0]; ptr[0]=x;
       x=ptr[2]; ptr[2]=ptr[1]; ptr[1]=x;
       ptr+=4;
     }
  }
  void WriteVTKFile(std::ostream& out) {
    int dimshape, conshape;
    DimShape(dimshape,conshape);
    bool writeSizeAsHeader=false;       // equivalent to addHeaderSizeAsInt32 on binary VTU
    bool flipdatabytes=((char *)(&conshape))[3]==0;

    int nnodes = NNodes(), ncells = NCells();
    float* xyz = new float[3*nnodes];
    out << "# vtk DataFile Version 3.0" << std::endl;
    out << "HDF data created with AMROC's hdf2file" << std::endl;
    out << "BINARY" << std::endl;
    out << "DATASET UNSTRUCTURED_GRID" << std::endl;
    out << "POINTS " << nnodes << " float" << std::endl;
    SetGridNodes((float (*)[3]) xyz);
    register int n;
    if (writeSizeAsHeader) {
      int binaryBlockSize= VTK_USED_DIMS*nnodes;
      out.write((char*)&binaryBlockSize,sizeof(int));
    }
    if (flipdatabytes) flipBytes((char*) xyz, VTK_USED_DIMS*nnodes*sizeof(float));
    out.write((char*) xyz, VTK_USED_DIMS*nnodes*sizeof(float));
    out << std::endl;
    delete [] xyz;

    int* con = new int[8*ncells];    
    int* con2 = new int[9*ncells];  
    int c2shape = conshape+1;
    out << "CELLS " << ncells << " " 
        << ncells*(conshape+1) << std::endl;
    SetGridConns((int (*)[8]) con);  
    for (n=0; n<ncells; n++) 
      con2[c2shape*n] = conshape;
    for (register int i=0; i<conshape; i++)
      for (n=0; n<ncells; n++) 
        con2[c2shape*n+i+1] = con[8*n+i];
    if (writeSizeAsHeader) {
      int binaryBlockSize= (conshape+1)*ncells;
      out.write((char*)&binaryBlockSize,sizeof(int));
    }
    if (flipdatabytes) flipBytes((char*) con2, (conshape+1)*ncells*sizeof(int));
    out.write((char*) con2, (conshape+1)*ncells*sizeof(int));
    out << std::endl;
    delete [] con; delete [] con2;
    
    int cell_type = 1;
    if (dimshape==3) cell_type = 11;
    else if (dimshape==2) cell_type = 8;
    else if (dimshape==1) cell_type = 3;
    int* type = new int[ncells];    
    for (n=0; n<ncells; n++)
      type[n] = cell_type;
    out << "CELL_TYPES " << ncells << std::endl;
    if (writeSizeAsHeader) {
      int binaryBlockSize=  ncells;
      out.write((char*)&binaryBlockSize,sizeof(int));
    }
    if (flipdatabytes) flipBytes((char*) type, ncells*sizeof(int));
    out.write((char*) type, ncells*sizeof(int));
    out << std::endl;
    delete [] type;

    int length;
    if (FKeys() == 1) {
      length = nnodes;
      out << "POINT_DATA ";
    }
    else {
      length = ncells;
      out << "CELL_DATA ";
    }
    out << length << std::endl;    
    float* dat = new float[nkeys*length];
    for (n=0; n<nkeys*length; n++) *(dat+n) = 0.;
    register int k;
    for (k=0; k<nkeys; k++) 
      SetScal(&(keys[k]),&(dat[k*length]));
    out << "FIELD Data " << nkeys << std::endl;  
    for (k=0; k<nkeys; k++) {
      char* title = new char[LEN_TKEYS];
      std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
      int i=std::strlen(title)-1;
      while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
      int n=std::strlen(title);
      for (i=0; i<n; i++) 
        if (title[i]==' ') title[i]='_';
      out << title << " 1 " << length << " float" << std::endl;
      if (writeSizeAsHeader) {
        int binaryBlockSize=  length;
        out.write((char*)&binaryBlockSize,sizeof(int));
      }

      if (flipdatabytes) flipBytes((char*) &(dat[k*length]), length*sizeof(float));
      out.write((char*) &(dat[k*length]), length*sizeof(float));
      out << std::endl;
    }
    delete [] dat;
  }    

  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  //                                                   WRITE SILO FILE
  //
  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#if defined(HAVE_LIBSILO) || defined(HAVE_LIBSILOH5)
  void WriteSiloFile(char *silofilename)  {
    if (silofilename==0 || silofilename[0]==0) {
      std::cerr << "Abort: Cannot write SILO to stdout. "
                << "Please provide FileName in your display_file.in" << std::endl;
      return;
    }
#if defined(HAVE_LIBSILOH5)
    DBfile * silofptr= DBCreate(silofilename, DB_CLOBBER, DB_LOCAL, NULL, DB_HDF5);
#else
    DBfile * silofptr= DBCreate(silofilename, DB_CLOBBER, DB_LOCAL, NULL, DB_PDB);
#endif
    if (!silofptr) {
      std::cerr << "Abort: could not create SILO file." << std::endl;
      return;
    }

    int dimshape, conshape;
    DimShape(dimshape,conshape);

    int nnodes = NNodes(), ncells = NCells();
    int vardatacentering, length;
    if (FKeys() == 1) {
      vardatacentering=DB_NODECENT;
      length = nnodes;
    } else {
      vardatacentering=DB_ZONECENT;
      length = ncells;
    }

    std::cerr << std::endl;

    float * xyz = new float[3*nnodes];
    char  * coordnames [3]= {"X","Y","Z"};
    float * coordinates[3];

    coordinates[0]=xyz;
    coordinates[1]=coordinates[0]+nnodes;
    coordinates[2]=coordinates[1]+nnodes;

#   ifdef USE_PRESET_AMROC_BLOCK_COORDS
    SetGridNodes((float (*)[3]) xyz);
    register int n,t;
    if (dimshape==1) {
      for (n=0,t=0; n<nnodes; n++)
        xyz[t++] = xyz[3*n];
    }
    else if (dimshape==2) {
      for (n=0,t=0; n<nnodes; n++) {
        xyz[t++] = xyz[3*n]; xyz[t++] = xyz[3*n+1];
      }
    }
#   endif

    float* dat = new float[nkeys*length];
    register int i,k;
    for (i=0; i<nkeys*length; i++) *(dat+i) = 0.;
    for (k=0; k<nkeys; k++) 
      SetScal(&(keys[k]),&(dat[k*length]));

    int    *blocklevel = new int[BlockList().size()];
    char  **vartitles  = new char* [nkeys];
    int     leveltitle = -1;

    // get variable name
    for (k=0; k<nkeys; k++) {
      char *title = new char[LEN_TKEYS];
      std::strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
      i=std::strlen(title)-1;
      while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
      int n=std::strlen(title);
      for (i=0; i<n; i++) {
        if (title[i]==' ' ||
            title[i]=='.' ||
            title[i]=='-' ||
            title[i]=='+' ||
            title[i]=='/'    ) title[i]='_';
      }
      vartitles[k] = new char[std::strlen(title)+1];
      for (i=0; i<n; i++) 
        vartitles[k][i]=std::toupper(title[i]);
      vartitles[k][i]='\0';
      delete [] title;
      if (!std::strcmp(vartitles[k],"LEVEL") || !std::strcmp(vartitles[k],"LEVELS"))
      leveltitle=k;
    }

    int offset=0, maxlevel=-1, nblock=0;
    char **meshptrs= new char *[BlockList().size()];
    char **varptrs=  new char *[BlockList().size()];

    for (typename block_list_type::iterator bit=BlockList().begin();
         bit!=BlockList().end(); bit++, nblock++) {
      float fmin[3], fmax[3];
      int idims[3], vardims[3];
      for (i=0; i<dimshape; i++) {
        fmin[i] = (**bit).DB().bbox().lower(i)*base::dx[i]+base::geom[i*2];
        fmax[i] = ((**bit).DB().bbox().upper(i)+
                   (**bit).DB().bbox().stepsize(i))*base::dx[i]+base::geom[i*2];
        idims[i]= (**bit).DB().bbox().extents(i)+1;
        vardims[i]=idims[i];
        if (vardatacentering==DB_ZONECENT) vardims[i]--;
        if (vardims[i]<=0) vardims[i]=1;
      }

      // creating a directory per block
      char *meshblockdirname= new char[LEN_TKEYS];
      meshblockdirname[0]='/';
      meshblockdirname++;
      sprintf(meshblockdirname,"BLOCK%04i",nblock);
      DBSetDir(silofptr,"/");
      DBMkDir(silofptr,meshblockdirname);
      DBSetDir(silofptr,meshblockdirname);

      char *meshblockname= meshblockdirname + std::strlen(meshblockdirname)+1;
      meshblockname[-1]='/';
      sprintf(meshblockname,"block%04i",nblock);

      meshptrs[nblock]=meshblockdirname;
      varptrs[nblock]= meshblockname;

#     ifdef USE_PRESET_AMROC_BLOCK_COORDS
      std::cerr<<"using USE_PRESET_AMROC_BLOCK_COORDS" << std::endl;
      for (i=0; i<dimshape; i++)
        coordinates[i]=&(xyz[dimshape*offset+i]);
      DBPutQuadmesh(silofptr, meshblockname ,coordnames, coordinates, idims, dimshape, DB_FLOAT, DB_COLLINEAR, NULL);
#     else
      for (i=0; i<dimshape; i++) {
        for (register int idx=0; idx<idims[i]; idx++)
          coordinates[i][idx]=fmin[i]+(float)idx*(fmax[i]-fmin[i])/(float)(idims[i]-1);
      }
      DBPutQuadmesh(silofptr, meshblockname ,coordnames, coordinates, idims, dimshape, DB_FLOAT, DB_COLLINEAR, NULL);
#     endif

      for (k=0; k<nkeys; k++) {
        if (k==leveltitle) {
          blocklevel[nblock]=(**bit).Level();
          if (blocklevel[nblock]>maxlevel) maxlevel=blocklevel[nblock];
        }
        DBPutQuadvar1(silofptr, vartitles[k], meshblockname, &(dat[offset+k*length]),
                      vardims, dimshape, NULL, 0,DB_FLOAT,  vardatacentering, NULL);
      }
      if (FKeys()==1) offset += (**bit).NNodes();
      else offset += (**bit).NCells();
    }

    bool   dump_sets=true;
    if (dump_sets) {
    DBSetDir(silofptr,"/");
    char * levelmesh = new char[LEN_TKEYS];
    char * varmeshname = new char[LEN_TKEYS];
 
    int   * meshtypes=new int[BlockList().size()];
    char ** imeshptrs=new char *[nblock];
    for (int i=0; i<nblock; i++) {
      meshtypes[i]= DB_QUAD_RECT;
      imeshptrs[i]=meshptrs[nblock-i-1];
    }
    DBPutMultimesh(silofptr, "/ALL_LEVELS", nblock, meshptrs, meshtypes, NULL);

    DBSetDir(silofptr,"/");
    DBMkDir(silofptr,"LEVELS");
    for (int l=maxlevel; l>=0; l--) {
      int ninlevel=0;
      for (register int idx=0; idx<nblock; idx++) {
        if (blocklevel[idx]==l) {
          imeshptrs[ninlevel++]=meshptrs[idx]-1;
        }
      }  
      if (ninlevel>0) {
        sprintf(levelmesh,"/LEVELS/LEVEL_%03i",l);
        DBPutMultimesh(silofptr, levelmesh, ninlevel, imeshptrs, meshtypes, NULL);
      }
    }

    for (register int idx=0; idx<nblock; idx++) meshtypes[idx]= DB_QUADVAR;
    for (k=0; k<nkeys; k++) {
      sprintf(varmeshname,"VAR_%s",vartitles[k]);
      for (register int idx=0; idx<nblock; idx++) std::strcpy(varptrs[idx],vartitles[k]);

      DBSetDir(silofptr,"/");
      DBMkDir(silofptr,varmeshname);
      DBSetDir(silofptr,varmeshname);
      sprintf(varmeshname,"/VAR_%s/ALL_LEVELS",vartitles[k]);
      for (register int idx=0; idx<nblock; idx++) imeshptrs[idx]=meshptrs[idx]-1;
      DBPutMultivar(silofptr, varmeshname, nblock, imeshptrs, meshtypes, NULL);
      DBSetDir(silofptr,"/");

      for (int l=maxlevel; l>=0; l--) {
        int ninlevel=0;
        for (register int idx=0; idx<nblock; idx++)
          if (blocklevel[idx]==l)
            imeshptrs[ninlevel++]=meshptrs[idx]-1;
        if (ninlevel>0) {
          sprintf(varmeshname,"/VAR_%s/LEVEL_%03i",vartitles[k],l);
          DBPutMultivar(silofptr, varmeshname, ninlevel, imeshptrs, meshtypes, NULL);
        }
      }
    }
    delete [] levelmesh;
    delete [] varmeshname;
    }
    DBClose(silofptr);


    //for (i=0; i<nblock; i++) if (meshptrs[i]) delete meshptrs[i];
    delete [] meshptrs;
    delete [] varptrs;
    delete [] xyz;
    delete [] dat;
    delete [] blocklevel;
    delete [] vartitles;

  }
#else
  void WriteSiloFile(char *silofilename)  {
     std::cerr << "hdf2file SILO package not compiled. "
               << "Please reconfigure with SILO_DIR defined and recompile." << std::endl;
  }  
#endif

  virtual void SetGrid(float xyz[][3]) {
    if (FKeys() == 1)
      SetGridNodes(xyz);
    else
      SetGridCells(xyz);
  }

  virtual void SetScal(int *key,float v[]) {
    if (FKeys() == 1)
      SetScalNodes(key, v);
    else
      SetScalCells(key, v);
  }

  virtual void Vect(int *key,float v[][3]) {
    if (FKeys() == 1)
      VectNodes(key, v);
    else
      VectCells(key, v);
  }

  virtual void BlockListAppend(hdata_block_type& workdata, BBox& bb, int& comp, int& Level,
                               int& proc, float* dx, float* geom) {
    BlockList().push_back(new data_block_type(*this, workdata, bb, comp,
                                              TheEvaluator(), Level, proc, dx, geom));
  }

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

  int LineProbe() const 
    { return (LinePoint[0]!=LineVector[0] || LinePoint[1]!=LineVector[1] || 
              LinePoint[2]!=LineVector[2]); }
 
  virtual void register_at(ControlDevice& Ctrl,const std::string& prefix) {
    base::register_at(Ctrl,prefix); }
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, ""); }
  virtual int NCells() { return base::NCells(); }
  virtual int NNodes() { return base::NNodes(); }
  virtual void SetBlocks(int blocks[]) { base::SetBlocks(blocks); }
  virtual void SetGridCells(float xyz[][3]) {
    base::SetGridCells(xyz); }
  virtual void SetGridNodes(float xyz[][3]) {
    base::SetGridNodes(xyz); }
  virtual void SetGridConns(int con[][8]) {
    base::SetGridConns(con); }
  virtual void SetScalCells(int *key,float v[]) {
    base::SetScalCells(key,v); }
  virtual void SetScalNodes(int *key,float v[]) {
    base::SetScalNodes(key,v); }
  virtual void VectNodes(int *key,float v[][3]) {
    base::VectNodes(key,v); }
  virtual void VectCells(int *key,float v[][3]) {
    base::VectCells(key,v); }
  virtual int Size() const { return base::Size(); }
  virtual const int& MinLevel() const { return base::MinLevel(); }
  virtual const int& MaxLevel() const { return base::MaxLevel(); }
  virtual const int& NLevels() const { return base::NLevels(); }
  virtual const int& RefineFactor(const int l) const {
    return base::RefineFactor(l); }
  virtual const BBox& GlobalBBox() const {
    return base::GlobalBBox(); }
  virtual const int& FKeys() const { return base::FKeys(); }
  virtual evaluator_type& TheEvaluator() const {
    return base::TheEvaluator(); }
  virtual evaluator_type& TheEvaluator() {
    return base::TheEvaluator(); }
  virtual float64 RealTime() const { return base::RealTime(); }
  virtual void SetUpdateName(std::string* iname) {
    base::SetUpdateName(iname); }
  virtual std::string* UpdateName() { return base::UpdateName(); }
  virtual std::string* UpdateName() const {
    return base::UpdateName(); }

private:
  std::string OutputName;
  int OutputType;   
  std::string StrKeys;   
  int* keys;
  int nkeys;
  float LinePoint[3], LineVector[3];
  int RegData[3];
  float RegLim[2];
  int RegMax;
  int LevelCutOut;
  int PointOutput;
  int _Step;
  int FileMerge;
  std::string PointFile;
};

#endif

---