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;
}
|