Google

i/hydra

yorick banner

Home

Manual

Packages

Global Index

Keywords

Quick Reference


/*
 * hydra.i   $Id$

 * functions to access hydra-generated Silo/PDB files
 */


local hydra ;
/* DOCUMENT hydra.i
   defines several functions useful for examining and extracting

   data from hydra-generated Silo/PDB dump files:


   openh       -- use instead of openb for hydra files

   hydra_xyz   -- extracts xyz and boundary arrays

   hydra_data  -- extracts data nodal or zonal arrays

   hydra_mix   -- extracts zonal data for mixed zones

   hydra_iparm -- extracts integer parameter values

   hydra_gblk  -- extracts information relating hblks to
                  user blocks


   The function hydra_old sets the variable naming as it was for
   dump files generated with codes built before about December 1998.
   You may also set _hy_auto_enable to autodetect file version.


   SEE ALSO: openh, hydra_xyz, hydra_data, hydra_mix,

             hydra_iparm, hydra_gblk, hydra_old
 */


func hydra_old (flag)
/* DOCUMENT hydra_old, 1
            hydra_old, 0


     use old (1) or new (0) hydra dump file variable names.

   SEE ALSO: hydra
 */
{
  { extern _hydra_old; }
  _hydra_old= flag;
}

_hydra_old= 0;
_hy_auto_enable= 0;  /* set to 1 for autodetect */

func _hy_auto_version (f)
{

  _hydra_old= (_hy_auto_enable? noneof((*get_vars(f)(1))=="/Global/") :
	       _hydra_old);
}

func _hblk_var (n,suffix)
{

  if (_hydra_old) return swrite(format="/BDATA/hblk%ld/"+suffix,n);

  else            return swrite(format="/hblk%ld/bdata/"+suffix,n);
}

func _hblk_bnd (n,b)
{

  if (_hydra_old) return swrite(format="/BDATA/hblk%ld/bnd%ld/",n,b);

  else            return swrite(format="/hblk%ld/bdata/bnd%ld/",n,b);
}

func _hblk_bc (n,b)
{

  if (_hydra_old) return swrite(format="/BDATA/hblk%ld/bc%ld/",n,b);

  else            return swrite(format="/hblk%ld/bdata/bc%ld/",n,b);
}

func _hy_global (suffix)
{
  if (_hydra_old) return "/BDATA/"+suffix;
  else            return "/Global/"+suffix;
}

func _hy_gmap (n)
{

  if (_hydra_old) return swrite(format="/BDATA/gmap%ld",n);

  else            return swrite(format="/Global/Decomposition/gmap%ld/gmap",n);
}

func _hy_umap (n)
{

  if (_hydra_old) return swrite(format="/BDATA/umap%ld",n);

  else            return swrite(format="/Global/Decomposition/umap%ld/umap",n);
}

func _root_gmap (n)
{

  return swrite(format="/Decomposition/gmap%ld/gmap",n);
}

func _root_umap (n)
{

  return swrite(format="/Decomposition/umap%ld/umap",n);
}


func openh (filename)
/* DOCUMENT f= openh(filename)


     sort of properly open a hydra-generated PDB file.


     Silo (or hydra) botches the Major-Order PDB primitive; the file

     needs to be opened with open102=1.
     Unfortunately, some variables (e.g.- iparmn) are written correctly...

     Presumably meshTV would break if hydra made the mesh array dimensions
     consistent with the iparmn dimensions, although my guess is it would
     survive switching Major-Order to 102.


   SEE ALSO: hydra_xyz, hydra_data
 */
{

  return openb(filename, open102=1);
}


func hydra_xyz (f, ublk, i0, j0, k0, face)
/* DOCUMENT mesh= hydra_xyz(f)

         or mesh= hydra_xyz(f, ublk, i0, j0, k0, face)

         or mesh= hydra_xyz(f, ublk, i0, j0, k0)


     read a 3D mesh object from the hydra PDB/Silo file F.

     The returned mesh is _lst(xyz, bound, mbnds, blks, start).

     Note that the boundary arrays are adjusted to the hex convention

     that cells with i=1, j=1, k=1 are missing, rather than the hydra
     convention that i=imax, j=jmax, k=kmax are missing.

     In the first form, the ray entry search will start on the

     first open boundary face in the mesh.  If the actual problem
     boundary is not convex, you need to identify a surface of
     constant i, j, or k in the problem which is convex, and which

     all the rays you intend to trace intersect.
     UBLK is the user block number (starting from 0),
     I0, J0, K0 are the (1-origin) logical coordinates of a

       hydra *cell*.  Note that unlike hex cells, the hydra
       cell bounded by nodes (1,1,1) and (2,2,2) is numbered (1,1,1).
       (Hex numbers it (2,2,2).)
     FACE is the face number on cell (I0,J0,K0) which you want a
       ray to enter.  0 means the -I face, 1 the +I face, 2 the -J
       face, 3 the +J face, 4 the -K face, and 5 the +K face.
       As you step from this cell to its neighbors, then to their
       neighbors, and so on, this face must trace out a convex
       surface for the ray entry search.  Rays not intersecting
       this surface will not enter the problem; the ray trace
       will begin at this surface, not at -infinity.

     If FACE==-1 or is omitted (as in the third form), then the

     given points on the rays are assumed to lie inside the mesh,
     and a pseudo ray from the centroid of cell (I0, J0, K0) will be
     tracked to the given point on each ray; the ray will be launched
     into the cell containing that point.


     You can set a hydra_bnd_hook function before calling hydra_xyz
     if the boundary conditions for hex need to be different than
     for hydra.


   SEE ALSO: hydra_bnd_hook, hydra_data, openh
 */
{
  _hy_auto_version, f;
  { local mdims, mlens; }

  nblk= hydra_blocks(f, mdims, mlens);
  f0= f;

  /* compute strides and global offsets for all the mesh blocks */
  blo= array(0, 4, nblk);
  blo(2:4,)= mdims;

  lscratch= hydra_blks(nblk, blo);


  /* allocate the global mesh array and read it from disk

   * read bc and bnd information at same time to avoid multiple
   * loops in the case of multiple files */

  if (nblk>1) xyz= array(0., 3, sum(mlens));

  else        xyz= array(0., 3, grow([3],mdims(,1)));
  off= 0;
  binfo= array(0, 4, nblk);
  bndtmp = array(pointer, nblk);
  bctmp = array(pointer, 2, nblk);
  for (i=1 ; i<=nblk ; i++) {

    bnum= hydrap_getblk(i, f);
    off0= off+1;
    off+= mlens(i);
    xyz(1,off0:off)=

      get_member(f,swrite(format="/hblk%ld/hydro_mesh_coord0",bnum))(*);
    xyz(2,off0:off)=

      get_member(f,swrite(format="/hblk%ld/hydro_mesh_coord1",bnum))(*);
    xyz(3,off0:off)=

      get_member(f,swrite(format="/hblk%ld/hydro_mesh_coord2",bnum))(*);
    /* read [nbc, nbnd, jp, kp] */

    binfo(,i) = get_member(f,_hblk_var(bnum,"hydrodati"))([34,37,21,22]);
    n = binfo(2,i);
    strides = [1,binfo(3,i),binfo(4,i)];
    if (n > 0) {
      tmp = array(pointer, 3, n);
      bndtmp(i) = &tmp;
      for (j=1 ; j<=n ; j++) {

        prefix = _hblk_bnd(bnum,j-1);

        pn = get_member(f,prefix+"pn")(1);

        list = where(abs(pn)==strides);

        if (!numberof(list)) continue;
        pn = sign(pn)*list(1);

        tmp(1,j) = &[i, get_member(f,prefix+"len_nsend1")(1),

                     pn, get_member(f,prefix+"blk_send")(1),

                     get_member(f,prefix+"bndnbc")(1)];

        tmp(2,j) = &long(get_member(f,prefix+"ndx_send"));

        tmp(3,j) = &long(get_member(f,prefix+"ndx_recv"));
      }
    }
    n= binfo(1,i);
    if (n > 0) {
      tmp = array(pointer, n);
      tmp2 = array(0, 2, n);
      bctmp(1,i) = &tmp;
      bctmp(2,i) = &tmp2;
      for (j=1 ; j<=n ; j++) {

        prefix = _hblk_bc(bnum,j-1);

        pn = get_member(f,prefix+"pn")(1);

        list = where(abs(pn)==strides);

        if (!numberof(list)) continue;
        pn = sign(pn)*list(1);
        tmp2(1,j) = pn;

        jj = get_member(f,prefix+"len")(1);  /* unnecessary?? */

        tmp(j) = &long(get_member(f,prefix+"ndx")(1:jj));

        tmp2(2,j) = get_member(f,prefix+"rtype")(1);
      }
    }
  }
  f= f0;

  bnd_off= binfo(2,cum);
  nbnd= max(bnd_off(0),1);

  bnd_blk= bnd_len= bnd_pn= bnd_r= bnd_ri= array(0, nbnd);

  bnd_ndxs= bnd_ndxr= array(pointer, nbnd);
  jj= 0;
  for (i=1 ; i<=nblk ; i++) {
    n= binfo(2,i);
    strides= [1,binfo(3,i),binfo(4,i)];
    eq_nocopy, tmp, *bndtmp(i);
    for (j=1 ; j<=n ; j++) {
      jj++;
      eq_nocopy, tmp2, *tmp(1,j);

      if (!numberof(tmp2)) continue;
      bnd_blk(jj)= tmp2(1);
      bnd_len(jj)= tmp2(2);
      bnd_pn(jj)= tmp2(3);
      bnd_r(jj)= tmp2(4);
      bnd_ri(jj)= tmp2(5);
      bnd_ndxs(jj)= tmp(2,j);
      bnd_ndxr(jj)= tmp(3,j);
    }
  }
  bndtmp = tmp = tmp2 = [];
  bnd_r++;
  bnd_ri++;


  /* allocate the block bound array and construct it from disk data */

  bound= array(0, dimsof(xyz));
  scratch= array(0, lscratch);

  /* number of bnd nodes is sum(bnd_len), which overestimates number
   * of faces; take that as a safe length to allocate */
  n= bnd_len(sum:1:jj);

  mbnds= n? array(HX_blkbnd, n) : [];
  for (i=1,ibnd=0 ; i<=jj ; i++) {
    if (!bnd_pn(i)) continue;
    n= bnd_len(i);
    s= bnd_blk(i);
    r= bnd_r(i);
    ri= bnd_off(r) + bnd_ri(i);
    bnds= [bnd_pn(i),binfo(3,s),binfo(4,s)];
    bndr= [bnd_pn(ri),binfo(3,r),binfo(4,r)];

    jbnd= hydra_bnd(ibnd, bound, scratch, blo(,s), blo(,r),
		    bnds, bndr, n, *bnd_ndxs(i), bnd_ndxr(ri), &mbnds, r-1);
    if (jbnd<0) {

      if (hydra_bnd_check) error, "illegal bnd data";
    } else {
      ibnd = jbnd;
    }
  }
  nbnds= ibnd;
  bnd_ndxs= bnd_ndxr= [];

  /* remains only to fill in the bcs in bound */
  start= -1;
  bndr= blor= [0];  /* unused */
  for (i=1,off=0 ; i<=nblk ; off+=mlens(i++)) {
    n= binfo(1,i);
    if (n <= 0) continue;
    strides= bnds= [1,binfo(3,i),binfo(4,i)];
    blos= blo(,i);

    tbound= array(0, 3,mdims(1,i),mdims(2,i),mdims(3,i));
    tcheck= array(0, 2,n);
    eq_nocopy, tmp, *bctmp(1,i);
    eq_nocopy, tmp2, *bctmp(2,i);
    for (j=1 ; j<=n ; j++) {
      local ndxs;
      eq_nocopy, ndxs, *tmp(j);

      if (!numberof(ndxs)) continue;
      pn = tmp2(1,j);
      bnds(1) = pn;
      rtype = tmp2(2,j);

      if (rtype==0) ibnd= -1;             /* open bc */
      else if (rtype==1) ibnd= -2;  /* reflecting bc */

      else if (rtype==3) ibnd= -3;   /* zero area bc */

      else error, "unknown rtype in bc data";
      if (!is_void(hydra_bnd_hook))
	ibnd= hydra_bnd_hook(ibnd, pn, i, j, ndxs);
      jj = numberof(ndxs);

      jj= hydra_mrk(ibnd, tbound, blos, bnds, jj, ndxs);

      if (jj < 0) error, "illegal bc data";
      tcheck(1,j)= pn;
      tcheck(2,j)= jj;
    }

    jj= hydra_adj(bound, tbound, blos, n, tcheck);
    tbound= [];
    if (jj>=0 && start<0) start= jj;
  }
  bctmp = tmp = tmp2 = [];


  blks= array(HX_block, nblk);
  blks.length= blo(2:4,);
  blks.stride(1,)= 1;
  blks.stride(2:3,)= blo(2:3,);
  blks.first= blo(1,);
  blks.final= blo(1,)+blo(4,);

  if (is_void(ublk)) {

    if (start<0) error, "mesh must have at least one open bc";
  } else {
    ublks= hydra_gblk(f, nblk);
    list= where(ublks(1,)==ublk);

    if (!numberof(list)) error, "mesh has no user block #"+pr1(ublk);
    ublks= ublks(,list);

    lst= where((ublks(2,)<=i0) & (ublks(3,)>i0) &
	       (ublks(4,)<=j0) & (ublks(5,)>j0) &
	       (ublks(6,)<=k0) & (ublks(7,)>k0));
    if (!numberof(lst))

      error, "mesh ublk #"+pr1(ublk)+" has no point "+pr1([i0,j0,k0]);
    if (numberof(lst)>1)

      error, "mesh ublk #"+pr1(ublk)+" has multiple points "+pr1([i0,j0,k0]);

    /* could read ireg and check that this cell exists, but don't bother */
    lst= lst(1);
    /* get 0-origin hex cell ijk within this block, convert to 1D index */
    start= [i0,j0,k0] - ublks(2:6:2,lst) + 1;
    block= blks(list(lst));

    start= block.first + sum(start*block.stride);

    if (!is_void(face) && face!=-1) {

      if (face<0 || face>5) error, "face must be between 0 and 5 inclusive";
      start= 6*start + face;
    } else {
      start= -1 - start;
    }
  }

  if (nbnds) mbnds= mbnds(1:nbnds);

  /* form the mesh */

  return _lst(xyz, bound, mbnds, blks, start);
}


func hydra_data (f, name)

/* DOCUMENT name_array = hydra_data(f, name)

         or pname_arrays = hydra_data(f, [name1,name2,...,nameN])

              eq_nocopy, name_array1, *pname_arrays(1)
              ...

              eq_nocopy, name_arrayN, *pname_arrays(N)


     reads variable NAME from the hydra file F.  If F is a multiblock
     file, NAME_ARRAY will be 1-D; for single block problems it will
     be 3-D.  If NAME=="matlist", you get the "Materials_matlist"
     array.  Coordinates can be obtained using the names x, y or z.


     In the second form, NAME1, ..., NAMEN are retrieved simultaneously,
     which is useful when F is a large family of files.

     Note that zone centered arrays are adjusted to the hex convention

     that cells with i=1, j=1, k=1 are missing, rather than the hydra
     convention that i=imax, j=jmax, k=kmax are missing.


   SEE ALSO: hydra_xyz, hydra_mix, hydra_array
 */
{
  { local mdims, mlens; }

  nblk = hydra_blocks(f, mdims, mlens);


  /* allocate the global data array and read it from disk */

  name = hydra_varname(name);
  dims = dimsof(name);
  nn = numberof(name);
  pdata = array(pointer, dims);
  data = [];
  off = 0;
  for (i=1 ; i<=nblk ; i++) {
    off0 = off+1;
    off += mlens(i);

    bnum = hydrap_getblk(i, f);
    for (j=1 ; j<=nn ; j++) {

      bdat = get_member(f,swrite(format="/hblk%ld/%s",bnum,name(j)));

      if (numberof(bdat)!=mlens(i)) {
        if (dimsof(bdat)(1)==1) {

          tmp = array(structof(bdat), mdims(,i)-1);
          tmp(*) = bdat;
          bdat = tmp;
        }

        tmp = array(structof(bdat), grow([3],mdims(,i)));
        tmp(2:0,2:0,2:0) = bdat;
        bdat=tmp;  tmp=[];
      }
      eq_nocopy, data, *pdata(j);
      if (is_void(data)) {

        if (nblk>1) data = array(structof(bdat), sum(mlens));

        else        data = array(structof(bdat), grow([3],mdims(,1)));
        pdata(j) = &data;
      }
      data(off0:off) = bdat(*);
    }
  }

  return dims(1)? pdata : data;
}


func hydra_varname (name)
{

  nm = array(string, dimsof(name));

  for (i=1 ; i<=numberof(nm) ; i++) {
    nami = name(i);
    if (nami=="matlist") nm(i) = "Materials_matlist";
    else if (nami=="x")  nm(i) = "hydro_mesh_coord0";
    else if (nami=="y")  nm(i) = "hydro_mesh_coord1";
    else if (nami=="z")  nm(i) = "hydro_mesh_coord2";
    else                 nm(i) = nami+"_data";
  }
  return nm;
}


func hydra_mix (f, &matlist, name, &mixdat)

/* DOCUMENT mixdat = hydra_mix(f, matlist)

              eq_nocopy, mixn, *mixdat(1)

              eq_nocopy, mixcell, *mixdat(2)

              eq_nocopy, mixnmat, *mixdat(3)

              eq_nocopy, mixhist, *mixdat(4)

         or mix_array = hydra_mix(f, mixdat, name)

         or pmix_array = hydra_mix(f, matlist, [name1,...,nameN], mixdat)

              eq_nocopy, mix_array1, *pmix_array(1)
              ...

              eq_nocopy, mix_arrayN, *pmix_array(N)


     In first form, returns MIXDAT and MATLIST for the hydra file F.
     MIXDAT consists of two arrays: MIXN is a list of the number of

     mixed cells for each block, and MIXCELL is an index array

     into any hex global cell array (as returned by hydra_data),
     MIXNMAT is the number of mix "zones" within each cell,
     and MIXHIST is the list required in order to use the

     histogram function on a mix array.


     In the second form, reads the mix data for the variable NAME

     in the hydra file F; the MIXDAT argument must have been returned

     by a previous call to hydra_mix using the first form.

     In the third form, MATLIST and MIXDAT are both returned along
     with the set of variables NAME1, ..., NAMEN, so that a number of

     variables can be retrieved in one call (useful when F is a large
     family of files).

     For example, to compute the temperature in each cell, using
     a mass weighted average in mixed zones, you would do this:

       den = hydra_data(f,"den");

       tmat = hydra_data(f,"tmat");

       mixdat = hydra_mix(f, matlist);
       local mixcell, mixhist;

       eq_nocopy, mixcell, *mixdat(2);

       eq_nocopy, mixhist, *mixdat(4);

       denx = hydra_mix(f, mixdat, "den");

       tmatx = hydra_mix(f, mixdat, "tmat");

       vf = hydra_mix(f, mixdat, "vf");
       tavg = tmat;

       tavg(mixcell) = histogram(mixhist, tmatx*denx*vf)/den(mixcell);


   SEE ALSO: hydra_xyz, hydra_data
 */
{
  { local mdims, mlens; }

  nblk = hydra_blocks(f, mdims, mlens);
  dims = dimsof(name);

  if (is_void(name) || dims(1)) {

    /* may as well collect matlist as long as we have to read it */

    if (nblk>1) matlist = array(0., sum(mlens));

    else        matlist = array(0., grow([3],mdims(,1)));
    offs = mlens(cum);

    if (!is_void(name)) {
      nn = numberof(name);
      nm = array(string, dims);
      for (j=1 ; j<=nn ; j++) {
        nami = name(j);
        if (nami=="vf") nm(j) = "Materials_mix_vf";
        else if (nami=="mat") nm(j) = "Materials_mix_mat";
        else nm(j) = nami+"_mix";
      }
      name = nm;
      pdata = array(pointer, dims);
      qdata = array(pointer, nblk);
      data = [];
    }

    /* scan blocks for mix data */
    nmix = array(0, nblk);

    cmix = lmix = array(pointer, nblk);
    moff = 0;
    for (i=1 ; i<=nblk ; i++) {

      bnum = hydrap_getblk(i, f);

      prefix = swrite(format="/hblk%ld/",bnum);

      list = get_member(f, prefix+"Materials_matlist");
      dims = grow([3],mdims(,i));
      tmp = array(0, dims-[0,1,1,1]);
      tmp(*) = list(*);
      list = array(0, dims);
      list(2:0,2:0,2:0) = tmp;  tmp = [];
      matlist(offs(i)+1:offs(i+1)) = list(*);


      cells = where(list<0);    /* cells with mixed materials */
      if (numberof(cells)) {
	start = -list(cells);   /* initial indices into mix_next (1-origin) */

	mix_next = get_member(f, prefix+"Materials_mix_next");
	/* mix_next is a packed set of sequences of indices into vf,
	 *   with each sequence terminated by a zero

	 * - apparently, only the zero markers are used by the hydra
	 * - mix_next is the same shape as the vf or *_mix variables
	 *
	 * start is guaranteed to be a monotonically increasing

	 *   (hydra creates it by incrementing through all zones),
	 *   so each contiguous block of vf or the *_mix variables
	 *   represents the next mixed cell in zone order */
	stop = where(mix_next<=0);
	n = stop - start + 1;  /* number of materials per mixed zone */
	/* convert cells to global hex cell indices */
	cells += offs(i);
      } else {
	cells = n = [];
      }
      cmix(i) = &cells;
      lmix(i) = &n;
      nmix(i) = numberof(cells);
      if (numberof(cells)) {
        data = array(pointer, nn);
        qdata(i) = &data;
        for (j=1 ; j<=nn ; j++)

          data(j) = &get_member(f,swrite(format="/hblk%ld/%s",bnum,name(j)));
      }
    }

    n = sum(nmix);
    if (n) {
      cells = n = array(0, n);
      offs = nmix(cum);
      for (i=1 ; i<=nblk ; i++) if (nmix(i)) {
	cells(offs(i)+1:offs(i+1)) = *cmix(i);
	n(offs(i)+1:offs(i+1)) = *lmix(i);
      }
      imix = n(cum);

      marks = array(0, imix(0));
      marks(imix(1:-1)+1) = 1;

      marks = marks(psum);       /* n(1) copies of 1, n(2) copies of 2
                                  * n(3) copies of 3, etc. */
      n = imix(offs+1)(dif);
    } else {
      cells = n = marks = [];
    }
    mixdat = [&nmix, &cells, &n, &marks];

    if (is_void(qdata)) return mixdat;
    if (is_void(n)) return pdata;

  } else {
    n = *matlist(3);
  }

  if (is_void(n)) return [];
  if (is_void(qdata)) {
    if (name=="vf") name = "Materials_mix_vf";
    else if (name=="mat") name = "Materials_mix_mat";
    else name = name+"_mix";
    nn = 1;
  }
  data = pbdat = [];
  off = 0;
  for (i=1 ; i<=nblk ; i++) {
    if (!n(i)) continue;
    off0 = off+1;
    off += n(i);

    if (!is_void(qdata)) eq_nocopy, pbdat, *qdata(i);
    for (j=1 ; j<=nn ; j++) {
      if (is_void(qdata)) {

        bnum = hydrap_getblk(i, f);

        bdat = get_member(f,swrite(format="/hblk%ld/%s",bnum,name));
      } else {

        eq_nocopy, bdat, *pbdat(j);

        eq_nocopy, data, *pdata(j);
      }
      if (is_void(data)) {

        data = array(structof(bdat), sum(n));

        if (!is_void(qdata)) pdata(j) = &data;
      }
      data(off0:off) = bdat(*);
    }
  }

  return is_void(qdata)? data : pdata;
}


struct HX_block {   /* must match hex.h */
  long stride(3);
  long length(3);
  long first;
  long final;
}


struct HX_blkbnd {  /* must match hex.h */
  long block;
  long cell;
  int orient;
}


func hydra_iparm (f, name)

/* DOCUMENT value= hydra_iparm(f, name)

         or names= hydra_iparm(f)


     returns value of hydra parameter NAME from file F,
     or a list of all names in NAME is not supplied.

     If NAME is not a string, returns that parameter
     or parameters (NAME is index in the returned list of names),

     for example hydra_iparm(f,1:0) returns all parameters.


   SEE ALSO: hydra_xyz, hydra_fparm
 */
{
  _hy_auto_version, f;
  hydrap_root, f;

  if (numberof(hydrap_blocks)) hydrap_getblk, 1, f;

  iparmn0= get_member(f,_hy_global("iparmn")); /* may be dimensioned wrong :(*/

  iparmn= array(char, 8, numberof(iparmn0)/8);
  iparmn(*)= iparmn0(*);
  if (!is_void(name)) {
    if (structof(name)==string) {
      name= *pointer(name);
      n= numberof(name);
      if (n>8) name= name(1:(n=8));

      i= where(!(iparmn(1:n,)!=name)(sum,))(1);
    } else {
      i= name;
    }

    result= get_member(f,_hy_global("iparmv"))(i);
    return result;

  } else {
    n= numberof(iparmn0)/8;
    iparmn0= array(string, n);
    for (i=1 ; i<=n ; i++) iparmn0(i)= string(&iparmn(,i));
    return iparmn0;
  }
}


func hydra_fparm (f, name)

/* DOCUMENT value= hydra_fparm(f, name)

         or names= hydra_fparm(f)


     returns value of hydra parameter NAME from file F,
     or a list of all names in NAME is not supplied.

     If NAME is not a string, returns that parameter
     or parameters (NAME is index in the returned list of names),

     for example hydra_fparm(f,1:0) returns all parameters.


   SEE ALSO: hydra_xyz, hydra_iparm
 */
{
  _hy_auto_version, f;
  hydrap_root, f;

  if (numberof(hydrap_blocks)) hydrap_getblk, 1, f;

  fparmn0= get_member(f,_hy_global("fparmn")); /* may be dimensioned wrong :(*/

  fparmn= array(char, 8, numberof(fparmn0)/8);
  fparmn(*)= fparmn0(*);
  if (!is_void(name)) {
    if (structof(name)==string) {
      name= *pointer(name);
      n= numberof(name);
      if (n>8) name= name(1:(n=8));

      i= where(!(fparmn(1:n,)!=name)(sum,))(1);
    } else {
      i= name;
    }

    result= get_member(f,_hy_global("fparmv"))(i);
    return result;

  } else {
    n= numberof(fparmn0)/8;
    fparmn0= array(string, n);
    for (i=1 ; i<=n ; i++) fparmn0(i)= string(&fparmn(,i));
    return fparmn0;
  }
}


func hydra_blocks (f, &mdims, &mlens)

/* DOCUMENT gnblk= hydra_blocks(f, mdims, mlens)

     returns number of blocks GNBLK, block dimensions MDIMS, and

     block lengths MLENS for the hydra mesh in file F.
     MDIMS is 3-by-NBLK, MLENS is GNBLK elements.


   SEE ALSO: hydra_xyz, hydra_iparm
 */
{
  hydrap_root, f;

  /* get number of blocks */
  if (!numberof(hydrap_blocks))

    gnblk= hydra_iparm(f, "gnblk");
  else
    gnblk= numberof(hydrap_blocks);

  /* get shape of each block */
  mlens= array(0, gnblk);
  mdims= array(0, 3, gnblk);
  if (hydra_rootfile(f)) {

    ndims = hydra_iparm(f, "ndims");
    for (i=1 ; i<=gnblk ; i++) {

      name = swrite(format="/Decomposition/gmap%ld/gmap", i-1);
      gl = get_member(f, name);
      mdims(1,i) = len = gl(7) - gl(4) + 1;
      if (ndims >= 2) {
        mdims(2,i) = d = gl(8) - gl(5) + 1;
        len *= d;
        if (ndims == 3) {
          mdims(3,i) = d = gl(9) - gl(6) + 1;
          len *= d;
        }
      }
      mlens(i) = len;
    }
  } else {
    for (i=1 ; i<=gnblk ; i++) {

      bnum= hydrap_getblk(i, f);

      name= swrite(format="/hblk%ld/hydro_mesh_coord0", bnum);

      mlens(i)= numberof(get_member(f,name));

      dims= dimsof(get_member(f,name));
      d= dims(1);
      mdims(1:d,i)= dims(2:d+1);
    }
  }

  return gnblk;
}


func hydra_ublk (f, nblk)
/* DOCUMENT ublk= hydra_ublk(f)


     return user block information from the hydra PDB/Silo file F.

     Each hblk in the mesh corresponds to a particular imin:imax,
     jmin:jmax, kmin:kmax in a particular ublk.  The return value is

     a 2D long array 7-by-numberof(h blocks):

     ublk(1,)=   user block number for this hblk
     ublk(2:3,)= ublk [imin,imax] of this hblk
     ublk(4:5,)= ublk [jmin,jmax] of this hblk
     ublk(6:7,)= ublk [kmin,kmax] of this hblk


   SEE ALSO: hydra_xyz, hydra_data, openh
 */
{
  _hy_auto_version, f;


  unblk = hydra_iparm(f,"unblk")

  ublk= array(0, 7, unblk);
  list= [1,4,7,5,8,6,9];
  if (hydra_rootfile(f)) {
    for (i=0 ; i<unblk ; i++) {

      umap = get_member(f,_root_umap(i))(1:9);
      ublk(,i+1) = umap(list);
    }
  } else {
    for (i=0 ; i<unblk ; i++) {

      bnum= hydrap_getblk(1, f);

      umap= get_member(f,_hy_umap(i))(1:9);
      ublk(,i+1)= umap(list);
    }
  }
  ublk(2:7,)-= 1;  /* adjust for ghost nodes and 0-origin vs 1-origin */

  return ublk;
}


func hydra_gblk (f, nblk)
/* DOCUMENT gblk= hydra_gblk(f)


     return global block information from the hydra PDB/Silo file F.

     Each hblk in the mesh corresponds to a particular imin:imax,
     jmin:jmax, kmin:kmax in a particular gblk.  The return value is

     a 2D long array 7-by-numberof(h blocks):

     gblk(1,)=   user block number for this hblk
     gblk(2:3,)= gblk [imin,imax] of this hblk
     gblk(4:5,)= gblk [jmin,jmax] of this hblk
     gblk(6:7,)= gblk [kmin,kmax] of this hblk


   SEE ALSO: hydra_xyz, hydra_data, openh
 */
{
  _hy_auto_version, f;
  if (is_void(nblk)) {
    { local mdims, mlens; }

    nblk= hydra_blocks(f, mdims, mlens);
  }

  gblk= array(0, 7, nblk);
  list= [1,4,7,5,8,6,9];
  if (hydra_rootfile(f)) {
    for (i=1 ; i<=nblk ; i++) {

      gmap = get_member(f,_root_gmap(i-1))(1:9);
      gnblk = gmap(2);
      gblk(,i) = gmap(list);
    }
  } else {
    for (i=1 ; i<=nblk ; i++) {

      bnum= hydrap_getblk(i, f);

      gnblk= get_member(f,_hblk_var(bnum,"hydrodati"))(36);

      gmap= get_member(f,_hy_gmap(gnblk))(1:9);
      gblk(,i)= gmap(list);
    }
  }
  gblk(2:7,)-= 1;  /* adjust for ghost nodes and 0-origin vs 1-origin */

  return gblk;
}


func hydra_rootfile (f)
{

  return strpart(print(f)(1),-4:0) == ".root";
}


func hydrap_root (f)
{
  { extern hydrap_pf, hydrap_unblk, hydrap_names, hydrap_blocks; }
  /* for now, use name -- other possibility is /_fileinfo == "dump file" */
  pf = print(f);
  if (pf(1) == hydrap_pf) return;
  hydrap_pf = pf(1);
  if (hydra_rootfile(f)) {

    hydrap_unblk= get_member(f,"/Decomposition/NumBlocks")(1);

    gnblk= get_member(f,"/Decomposition/NumDomains")(1);

    hydrap_names= array(string, gnblk);

    toks= [string(0), string(&get_member(f,"/hydro_mesh_meshnames"))];
    for (i=1 ; i<=gnblk ; i++) {
      toks= strtok(toks(2),";");

      hydrap_names(i)= strtok(toks(1),":")(1);
    }

    hydrap_names= strpart(pf(2),17:0) + hydrap_names; /* add directory part */

    hydrap_blocks= get_member(f,"/DomainFiles")(1:gnblk);
  } else {
    hydrap_unblk= hydrap_names= hydrap_blocks= hydrap_curname= [];
  }
}


func hydrap_getblk (i, &f)
{
  { extern hydrap_curname; }
  if (numberof(hydrap_blocks)) {
    if (i==1) hydrap_curname= string(0);
    if (hydrap_curname != hydrap_names(i)) {
      hydrap_curname= hydrap_names(i);
      f= openh(hydrap_curname);
    }
    return hydrap_blocks(i);
  } else {
    return i-1;
  }
}


func hydra_array (f, ublk, name)

/* DOCUMENT name_array = hydra_array(f, ublk, name)

         or pname_arrays = hydra_array(f, ublk, [name1,name2,...,nameN])

              eq_nocopy, name_array1, *pname_arrays(1)
              ...

              eq_nocopy, name_arrayN, *pname_arrays(N)


     reads variable array NAME for user block UBLK from the hydra file F.  
     If NAME=="matlist", you get the "Materials_matlist" array.
     Coordinates can be obtained using the names x, y or z.

     Ublk numbering starts at 0.


     Note that here zone centered arrays are given using the hydra convention
     so that i=imax, j=jmax, k=kmax are missing.  Thus in order to use the 

     Yorick plc and plf functions correctly you should index the plotted
     variable i.e. for a 2D array.
     plf, den(1:-1,1:-1), y, x


   SEE ALSO: hydra_xyz, hydra_mix
 */
{
  { local mdims, mlens; }

  gnblk= hydra_blocks(f, mdims, mlens);

  unblk= hydra_iparm(f,"unblk");

  ublk++;

  if (ublk > unblk) {

    print,"Error bad ublk number specified - only %d user blocks,"+
      " numbering starts from 0",unblk;
    return ;
  }


  /* allocate the global data array and read it from disk */

  name = hydra_varname(name);
  dims = dimsof(name);
  nn = numberof(name);
  pdata = array(pointer, dims);
  data = [];
  off= 0;
  umap = hydra_ublk(f,gnblk);
  gmap = hydra_gblk(f,gnblk);

  ndims = hydra_iparm(f,"ndims");
  if (ndims == 2)

    udims= grow([2],umap(3,ublk),umap(5,ublk)) ;
  else

    udims= grow([3],umap(3,ublk),umap(5,ublk),umap(7,ublk)) ;
  for (i=1 ; i<=gnblk ; i++) {
    if ((gmap(1,i)+1) != ublk) continue;

    bnum= hydrap_getblk(i, f);
    for (j=1 ; j<=nn ; j++) {

      bdat= get_member(f,swrite(format="/hblk%ld/%s",bnum,name(j)));

      eq_nocopy, data, *pdata(j);
      if (is_void(data)) {

        data = array(structof(bdat), udims);
        pdata(j) = &data;
      }

      dm = dimsof(bdat);
      ix = dm(2);
      jx = dm(3);
      ioff = gmap(2,i)-1;
      joff = gmap(4,i)-1;
      if (ndims == 2) {
        data(ioff+1:ioff+ix,joff+1:joff+jx) = bdat;
      } else {
        kx = dm(4);
        koff = gmap(6,i)-1;
        data(ioff+1:ioff+ix,joff+1:joff+jx,koff+1:koff+kx) = bdat;
      }
    }
  }

  return dims(1)? pdata : data;
}