FileVisualGridToSilo.h

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Santiago Lombeyda

# ifndef strequals
# define strequals(src,dst) (!strcmp(src,dst))
# endif

  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  //                                                   WRITE SILO FILE
  //
  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  // min/max bounds mark lower front left corner of cells (not vertex!)
  // nvertex*=(dims[i]+1)
  // stepsize marks the number of finest mesh that could fit in 
  //          current cell size
  // dx marks the conversion from integer/double cell numbering to real
  //    world units (without these the coords make a cube)
  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  void WriteSiloFile(char *silofilename)  {
#   if DEBUG
#    ifdef OBNOXIOUS_DEBUG
      bool verbose=true;
#    else
      bool verbose=false;
#    endif
#   else
    bool verbose=false;
#   endif

    using std::malloc;
    using std::toupper;

    bool matchCoordAndVarDims=true;
    int ndims=3;
    int nvardims;
    int nodespercell;

    DimShape(nvardims,nodespercell);

    if (matchCoordAndVarDims) ndims=nvardims;

    if (silofilename==NULL || silofilename[0]==0) {
      cerr << "Abort: currentl cannot write SILO to stdout. "
                     "Please provide FileName in your display_file.in" << endl;
      return;
    }

    if (verbose) {
      cerr << endl << LEN_TKEYS << endl;
      cerr << "Starting convesion to Silo." << endl;
    }


    SiloCreator * silofptr= new SiloCreator(silofilename);

    if (silofptr==NULL) {
      cerr << "Abort: could not create SILO file." << endl;
      return;
    }

    if (verbose) {
      cerr << "Volume bounds " <<  geom[0] << ","
                               <<  geom[2] << ","
                               <<  geom[4] << " - "
                               <<  geom[1] << ","
                               <<  geom[3] << ","
                               <<  geom[5] << "  ::  "
                               <<  dx[0] << ","
                               <<  dx[1] << ","
                               <<  dx[2] << endl;
    }

    char   * meshblockname;
    char   * meshblockdirname;
    char   * title;
    double   min[3],max[3];
    int      dims[3]={0,0,0};
    float    fmin[3];
    float    fmax[3];
    double   spacing[3];
    int      idims[3];
    int      vardims[3];
    int      nchunk=0;
    int      nblock=0;
    int      maxlevel=-1;

    float     **blockptrs= (float **)malloc(nkeys*BlockList().size()*sizeof(float *));
    int        *blocklevel=(int    *)malloc(BlockList().size()*sizeof(int));
    char      **vartitles= (char  **)malloc(sizeof(char *)*nkeys);

    SiloMesh     **meshptrs= (SiloMesh **)malloc(BlockList().size()*sizeof(SiloMesh *));
    SiloVariable **varptrs=  (SiloVariable  **)malloc(BlockList().size()*sizeof(SiloVariable  *)*nkeys);
    int            vardatacentering=DB_NODECENT;

    for (typename block_list_type::iterator bit=BlockList().begin();
                  bit!=BlockList().end();
                  bit++, nblock++) {
      int length;

      if (FKeys() == 1) {
        vardatacentering=DB_NODECENT;
        length = (**bit).NNodes();
      } else {
        vardatacentering=DB_ZONECENT;
        length = (**bit).NCells();
      }

      float* dat = new float[length*nkeys];

      if (dat==NULL) {
        cerr<<"Abort: could not allocate memory for block no." << nblock << endl; 
        return;
      }

      if (verbose) {
        cerr << "Encoding block #" << nblock << " of size " << length << endl;
      }

      for (int i=0; i<ndims; i++) {

        min[i]=(**bit).DB().bbox().lower(i);
        max[i]=(**bit).DB().bbox().upper(i);

        // if dumping as cells
        dims[i]=(**bit).DB().bbox().extents(i);
        dims[i]++;

        // spacing in terms of min refinement
        spacing[i]=(**bit).DB().bbox().stepsize(i);
        max[i]+=spacing[i];

        // debug check by ralf (?)
        if (verbose) {
          if (max[i]>48) {
            cerr << "for axis " << i << " got max " << max[i] << endl;
          }
        }

        fmin[i]=(float)(min[i]*(double)dx[i]+(double)geom[i*2]);
        fmax[i]=(float)(max[i]*(double)dx[i]+(double)geom[i*2]);
        idims[i]=(int)dims[i];
        vardims[i]=idims[i];
        if (vardatacentering==DB_ZONECENT) vardims[i]--;
        if (vardims[i]<=0) vardims[i]=1;

        // debug check by ralf (?)
        if (verbose) {
          if (fmax[i]>geom[i*2+1]) {
            cerr << "for axis " << i << " got fmax " << fmax[i] << endl;
          }
        }
      }

      if (verbose) {
        cerr << "   [" << min[0] << "," <<
                          min[1] << "," <<
                          min[2] << "] - [" <<
                          max[0] << "," <<
                          max[1] << "," <<
                          max[2] << "]    steps of (" <<
                          spacing[0] << "," <<
                          spacing[1] << "," <<
                          spacing[2] << ")" << endl;
      }

      // creating a directory per block
      meshblockdirname= new char[LEN_TKEYS];

      sprintf(meshblockdirname,"BLOCK%04i",nblock);
      silofptr->SetToRoot();
      silofptr->MakeAndSetDirectory(meshblockdirname);

      meshblockname= new char[LEN_TKEYS];
      sprintf(meshblockname,"block%04i",nblock);

      if (verbose) {
        cerr << "* ";
        cerr << meshblockname << ":" << idims[0] << "x" << idims[1] << "x" << idims[2] << endl;
        cerr << "    ";
        cerr << fmin[0] << "," << fmin[1] << "," << fmin[2] << " - ";
        cerr << fmax[0] << "," << fmax[1] << "," << fmax[2] << endl;
      }

#     ifdef USE_PRESET_AMROC_BLOCK_COORDS
      // this code is for flattened unstructured grids
      // : not current valid
      float* xyz = new float[3*(**bit).DB().bbox().size()];
      SetGridNodes((float (*)[3]) xyz);
      if (ndims<3)
        for (int idx=0,cidx=0; cidx<NNodes();cidx++) 
          for (int icidx=0; icidx<ndims; icidx++)
            xyz[idx++]=xyz[cidx*3+icidx];
      meshptrs[nblock]= new SiloMesh((char *)meshblockname,ndims,(int *)idims,(float *)xyz);
#     else
      meshptrs[nblock]= new SiloMesh((char *)meshblockname,ndims,(int *)idims,(float *)fmin,(float *)fmax);
#     endif
      silofptr->SetMesh(meshptrs[nblock]);
      silofptr->DrawMesh();

      int offset=0;
      for (register int k=0; k<nkeys; k++) {

        blockptrs[nchunk]=dat+offset;

        int poffset=offset;
        if (FKeys() == 1) {
          // Node-centered
          (**bit).SetScalNodes(&(keys[k]),dat,offset,base::CompRead);
        } else  {
          // Cell-centered
          (**bit).SetScalCells(&(keys[k]),dat,offset,base::CompRead); 
        }
        if (offset==poffset) {
          std::cerr << "Warning: SetScal::CompRead of key #" << k << " created no data offset. Variable may not be defined." << std::endl;
        }


        // get variable name
        title = new char[LEN_TKEYS];
          strcpy(title,TheEvaluator().TitleKeys(keys[k]-1));
          int i=strlen(title)-1;
          while (title[i]==' ' && i>1) { title[i] = '\0'; i--; }
          int n=strlen(title);
          for (i=0; i<n; i++) {
            if (title[i]==' ') title[i]='_';
            if (title[i]=='-') title[i]='_';
            if (title[i]=='+') title[i]='_';
            if (title[i]=='/') title[i]='_';
          }

          // Replace GNU extension strdup() with standard function strcpy()
          // char * varTITLE=strdup(title);
          char* varTITLE = (char*) malloc(strlen(title)+1);
          strcpy(varTITLE,title);
        for (unsigned int i=0; i<strlen(title); i++) varTITLE[i]=toupper(title[i]);

        if (nblock==0) {
          // Replace GNU extension strdup() with standard function strcpy()
          // vartitles[nchunk]=strdup(varTITLE);
          vartitles[nchunk] = (char*) malloc(strlen(varTITLE)+1);
          strcpy(vartitles[nchunk],varTITLE);
        }

        if (strequals(varTITLE,"LEVEL") || strequals(varTITLE,"LEVELS")) {
          blocklevel[nblock]=(**bit).Level();
          if (blocklevel[nblock]>maxlevel) maxlevel=blocklevel[nblock];
        } 
#       ifdef DEBUG_REPLACE_VALUES
        if (nvardims==3) {
          int idx=0;
          for (int idx_k=0; idx_k<idims[2]; idx_k++)
            for (int idx_j=0; idx_j<idims[1]; idx_j++)
              for (int idx_i=0; idx_i<idims[0]; idx_i++,idx++)
                blockptrs[nchunk][idx]=(((k%2)==1)?k:-k)*(min[0]+min[1]+min[2]+idx_i+idx_j+idx_k);
                //(blockptrs[nchunk])[idx]=k;
        }
#       endif
        varptrs[nchunk]=new SiloVariable(varTITLE,nvardims,vardims,meshptrs[nblock],1,blockptrs+nchunk,vardatacentering);
        silofptr->SetVariable(varptrs[nchunk]);
        silofptr->DrawVariable();
        if (verbose) cerr << "    +" << varTITLE << endl;
        nchunk++;
      }
      //for (int idx=0; idx<nkeys; idx++) cerr << "    +" << dat[idx*length+1] << endl;
    }
    
#   ifndef SILO_NO_GROUPING
    silofptr->SetToRoot();
    silofptr->CreateMultiMesh("/ALL_LEVELS");
      SiloMesh **imeshptrs= (SiloMesh **)malloc(BlockList().size()*sizeof(SiloMesh *));
      for (int i=0; i<nblock; i++) imeshptrs[i]=meshptrs[nblock-i-1];
      silofptr->SetMultiMesh(nblock,imeshptrs);
    silofptr->WriteMultiMesh();

    silofptr->SetToRoot();
    silofptr->MakeAndSetDirectory("LEVELS");
 
    if (verbose) cerr << "Creating per level collections." << endl;
    for (int i=maxlevel; i>=0; i--) {
      int ninlevel=0;
      for (int idx=0; idx<nblock; idx++) if (blocklevel[idx]==i) ninlevel++;
      if (ninlevel>0) {
        if (verbose) cerr << "Creating mesh collection for level " << i << " [" << ninlevel << "] ";
        SiloMesh ** blocksinlevel= (SiloMesh **)malloc(sizeof(SiloMesh *)*ninlevel);
        for (int idx=0, ni=0; idx<nblock; idx++)
          if (blocklevel[idx]==i)
            blocksinlevel[ni++]=meshptrs[idx];
        char * levelmesh = new char[LEN_TKEYS];
        sprintf(levelmesh,"/LEVELS/LEVEL_%03i",i);
        if (verbose) cerr << levelmesh << endl;
        silofptr->CreateMultiMesh(levelmesh);
        silofptr->SetMultiMesh(ninlevel,blocksinlevel);
        silofptr->WriteMultiMesh();
      }
    }

    int nvars= nkeys;
    for (int i=0; i<nvars; i++) {
      SiloVariable ** blocksinvar= (SiloVariable **)malloc(sizeof(SiloVariable *)*nblock);
      //for (int ni=0, idx=i; idx<nchunk; idx+=nvars) blocksinvar[ni++]=varptrs[idx];
      for (int ni=nblock, idx=i; idx<nchunk; idx+=nvars) blocksinvar[--ni]=varptrs[idx];

      char * varmeshname = new char[LEN_TKEYS];
      silofptr->SetToRoot();
      sprintf(varmeshname,"VAR_%s",vartitles[i]);
      silofptr->MakeAndSetDirectory(varmeshname);

      sprintf(varmeshname,"/VAR_%s/ALL_LEVELS",vartitles[i]);

      silofptr->CreateMultiVar(varmeshname);
      silofptr->SetMultiVar(nblock,blocksinvar);
      silofptr->WriteMultiVar();

      for (int l=maxlevel; l>=0; l--) {
        int ninlevel=0;
        for (int idx=0; idx<nblock; idx++)
          if (blocklevel[idx]==l) ninlevel++;
        if (ninlevel>0) {
          SiloVariable ** blocksinlevel= (SiloVariable **)malloc(sizeof(SiloVariable *)*ninlevel);
        for (int idx=0, ni=0; idx<nblock; idx++)
          if (blocklevel[idx]==l)
            blocksinlevel[ni++]=varptrs[idx*nvars+i];
          sprintf(varmeshname,"/VAR_%s/LEVEL_%03i",vartitles[i],l);
          if (verbose)
            cerr << "Writing collection " << varmeshname << " " << ninlevel << "/" << nblock << endl;
          silofptr->CreateMultiVar(varmeshname);
          silofptr->SetMultiVar(ninlevel,blocksinlevel);
          silofptr->WriteMultiVar();
        }
      }
    }
#   endif

    //for (int idx=0; idx<nchunk; idx++) fprintf(stderr,"%8i %f\n",idx,blockptrs[idx][1]);

    if (verbose) cerr << "Closing silo file " << silofilename << endl ;
    silofptr->Close();

    if (verbose) cerr << "DONE" << endl;

    //delete silofptr;

    cerr << endl;

  }

---