[Matrix-SIG] Errors in Gist files in LLNLPython2
Paul F. Dubois
Paul F. Dubois" <dubois1@llnl.gov
Tue, 9 Jun 1998 14:19:49 -0700
This is a multi-part message in MIME format.
------=_NextPart_000_0011_01BD93B1.A9ACA3C0
Content-Type: text/plain;
charset="iso-8859-1"
Content-Transfer-Encoding: 7bit
I accidentally missed the one Python file in the LLNL distribution that used
the routines whose names were changed to bitwise_and, bitwise_or, and
bitwise_xor. File Graphics/OOG/Lib/slice3.py needs all "boolean_and" changed
to "bitwise_and".
The corrected version is attached for you convenience. Since this is
something that one can fix in ten seconds with a text editor I won't make a
new release just to fix it (as I expect to make a new release relatively
soon anyway). I apologize for any inconvenience.
Apparently some previous version of NumPy actually had these functions
working. I think we probably lost them on the "Hinsen change" where the
repairs were made to enable NumPy to be a shared library on more platforms.
If your software like Gist was using them, I just broke your software with
this name change. Since nobody who used these functions evidently upgraded
to the previous NumPy, there can't be too many of you...but my apologies
anyway. I wanted this change in before we wrote a lot of documentation.
The second file is one that Zane had announced a bug fix for, and I thought
he had checked in, but he hadn't. It is Graphics/Gist/Lib/gist.py.
------=_NextPart_000_0011_01BD93B1.A9ACA3C0
Content-Type: application/octet-stream;
name="slice3.py"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="slice3.py"
# Copyright (c) 1996, 1997, The Regents of the University of California.
# All rights reserved. See LEGAL.LLNL for full text and disclaimer.
# SLICE3.PY
# find 2D slices of a 3D hexahedral mesh
# $Id: slice3.py,v 1.9 1997/11/12 21:51:52 motteler Exp $
#
# Change (ZCM 12/4/96) Apparently _draw3_list, which is global in
# pl3d.py, can be fetched from there, but once this has been done,
# assignments to it over there are not reflected in the copy here.
# This has been fixed by creating an access function.
from Numeric import *
from shapetest import *
from types import *
from pl3d import *
from arrayfns import *
from gistC import *
#
# Caveats:
# (A) Performance is reasonably good, but may still be a factor of
# several slower than what could be achieved in compiled code.
# (B) Only a simple in-memory mesh model is implemented here.
# However, hooks are supplied for more interesting possibilities,
# such as a large binary file resident mesh data base.
# (C) There is a conceptual difficulty with _walk3 for the case
# of a quad face all four of whose edges are cut by the slicing
# plane. This can only happen when two opposite corners are
# above and the other two below the slicing plane. There are
# three possible ways to connect the four intersection points in
# two pairs: (1) // (2) \\ and (3) X There is a severe problem
# with (1) and (2) in that a consistent decision must be made
# when connecting the points on the two cells which share the
# face - that is, each face must carry information on which way
# it is triangulated. For a regular 3D mesh, it is relatively
# easy to come up with a consistent scheme for triangulating faces,
# but for a general unstructured mesh, each face itself must carry
# this information. This presents a huge challenge for data flow,
# which I don't believe is worthwhile. Because the X choice is
# unique, and I don't see why we shouldn't use it here.
# For contouring routines, we reject the X choice on esthetic
# grounds, and perhaps that will prove to be the case here as
# well - but I believe we should try the simple way out first.
# In this case, we are going to be filling these polygons with
# a color representing a function value in the cell. Since the
# adjacent cells should have nearly the same values, the X-traced
# polygons will have nearly the same color, and I doubt there will
# be an esthetic problem. Anyway, the slice3 implemented
# below produces the unique X (bowtied) polygons, rather than
# attempting to choose between // or \\ (non-bowtied) alternatives.
# Besides, in the case of contours, the trivial alternating
# triangulation scheme is just as bad esthetically as every
# zone triangulated the same way!
def plane3 (normal, point) :
# plane3(normal, point)
# or plane3([nx,ny,nz], [px,py,pz])
# returns [nx,ny,nz,pp] for the specified plane.
# the normal doesn't really need to be normalized, but this
# has the desirable side effect of blowing up if normal=3D=3D0
newnorm =3D zeros (4, Float)
newnorm [0:3] =3D normal / sqrt (sum (normal*normal))
newnorm [3] =3D sum (multiply (normal, point))
return newnorm
_Mesh3Error =3D "Mesh3Error"
def mesh3 (x, y =3D None, z =3D None, ** kw) :
# mesh3(x,y,z)
# or mesh3(x,y,z, funcs =3D [f1,f2,...])
# or mesh3(xyz, funcs =3D [f1,f2,...])
# or mesh3(nxnynz, dxdydz, x0y0z0, funcs =3D [f1,f2,...])
# make mesh3 argument for slice3, xyz3, getv3, etc., functions.
# X, Y, and Z are each 3D coordinate arrays. The optional F1, F2,
# etc. are 3D arrays of function values (e.g. density, temperature)
# which have one less value along each dimension than the coordinate
# arrays. The "index" of each zone in the returned mesh3 is
# the index in these cell-centered Fi arrays, so every index from
# one through the total number of cells indicates one real cell.
# The Fi arrays can also have the same dimensions as X, Y, or Z
# in order to represent point-centered quantities.
# If X has four dimensions and the length of the first is 3, then
# it is interpreted as XYZ (which is the quantity actually stored
# in the returned cell list).
# If X is a vector of 3 integers, it is interpreted as [nx,ny,nz]
# of a uniform 3D mesh, and the second and third arguments are
# [dx,dy,dz] and [x0,y0,z0] respectively. (DXDYDZ represent the
# size of the entire mesh, not the size of one cell, and NXNYNZ are
# the number of cells, not the number of points.)
# Added by ZCM 1/13/97: if x, y, and z are one-dimensional of
# the same length and if the keyword verts exists and yields
# an NCELLS by 8 integer array, then we have an unstructured
# rectangular mesh, and the subscripts of cell i's vertices
# are verts[i, 0:8].
# Added by ZCM 10/10/97: if x, y, and z are one-dimensional
# of the same length or not, and verts does not exist, then
# we have a structured reectangular mesh with unequally spaced
# nodes.
# other sorts of meshes are possible -- a mesh which lives
# in a binary file is an obvious example -- which would need
# different workers for xyz3, getv3, getc3, and iterator3
# iterator3_rect may be more general than the other three;
# as long as the cell dimensions are the car of the list
# which is the 2nd car of m3, it will work=20
dims =3D shape (x)
if len (dims) =3D=3D 1 and y !=3D None and len (x) =3D=3D len (y) \
and z !=3D None and len(x) =3D=3D len (z) and kw.has_key ("verts") =
:
virtuals =3D [xyz3_irreg, getv3_irreg,
getc3_irreg, iterator3_irreg]
dims =3D kw ["verts"]
if type (dims) !=3D ListType :
m3 =3D [virtuals, [dims, array ( [x, y, z])], []]
else : # Irregular mesh with more than one cell type
sizes =3D ()
for nv in dims :
sizes =3D sizes + (shape (nv) [0],) # no. cells of this type
totals =3D [sizes [0]]
for i in range (1, len (sizes)) :
totals.append (totals [i - 1] + sizes [i]) #total cells so =
far
m3 =3D [virtuals, [dims, array ( [x, y, z]), sizes, totals], =
[]]
if kw.has_key ("funcs") :
funcs =3D kw ["funcs"]
else :
funcs =3D []
i =3D 0
for f in funcs:
if len (f) !=3D len (x) and len (f) !=3D shape (dims) [0] :
# if vertex-centered, f must be same size as x.
# if zone centered, its length must match number of cells.
raise _Mesh3Error, "F" + `i` + " is not a viable 3D cell =
value"
m3 [2] =3D m3 [2] + [f]
i =3D i + 1
return m3
virtuals =3D [xyz3_rect, getv3_rect, getc3_rect, iterator3_rect]
if len (dims) =3D=3D 4 and dims [0] =3D=3D 3 and min (dims) >=3D 2 :
xyz =3D x
dims =3D dims [1:4]
elif len (dims) =3D=3D 1 and len (x) =3D=3D 3 and type (x [0]) =3D=3D =
IntType \
and y !=3D None and z !=3D None and len (y) =3D=3D len (z) =3D=3D =
3 :
xyz =3D array ([y, z])
dims =3D (1 + x [0], 1 + x [1], 1 + x [2])
virtuals [0] =3D xyz3_unif
elif len (dims) =3D=3D 1 and y !=3D None and z !=3D None and len =
(y.shape) =3D=3D 1 \
and len (z.shape) =3D=3D 1 and x.typecode () =3D=3D y.typecode () =
=3D=3D \
z.typecode () =3D=3D Float :=20
# regular mesh with unequally spaced points
dims =3D array ( [len (x), len (y), len (z)], Int)
xyz =3D [x, y, z] # has to be a list since could be different =
lengths
virtuals [0] =3D xyz3_unif
else :
if len (dims) !=3D 3 or min (dims) < 2 or \
y =3D=3D None or len (shape (y)) !=3D 3 or shape (y) !=3D dims =
or \
z =3D=3D None or len (shape (z)) !=3D 3 or shape (z) !=3D dims:
raise _Mesh3Error, "X,Y,Z are not viable 3D coordinate mesh =
arrays"
xyz =3D array ( [x, y, z])
dim_cell =3D (dims [0] - 1, dims [1] - 1, dims [2] - 1)
m3 =3D [virtuals, [dim_cell, xyz], []]
if kw.has_key ("funcs") :
funcs =3D kw ["funcs"]
else :
funcs =3D []
i =3D 0
for f in funcs:
if len (f.shape) =3D=3D 3 and \
( (f.shape [0] =3D=3D dims [0] and f.shape [1] =3D=3D dims [1] =
and
f.shape [2] =3D=3D dims [2]) or (f.shape [0] =3D=3D dim_cell =
[0] and
f.shape [1] =3D=3D dim_cell [1] and f.shape [2] =3D=3D =
dim_cell [2])) :
m3 [2] =3D m3 [2] + [f]
i =3D i + 1
else :
raise _Mesh3Error, "F" + `i` + " is not a viable 3D cell value"
return m3
# Ways that a list of polygons can be extracted:
# Basic idea:
# (1) At each *vertex* of the cell list, a function value is defined.
# This is the "slicing function", perhaps the equation of a =
plane,
# perhaps some other vertex-centered function.
# (2) The slice3 routine returns a list of cells for which the
# function value changes sign -- that is, cells for which some
# vertices have positive function values, and some negative.
# The function values and vertex coordinates are also returned.
# (3) The slice3 routine computes the points along the *edges*
# of each cell where the function value is zero (assuming linear
# variation along each edge). These points will be vertices of
# the polygons. The routine also sorts the vertices into cyclic
# order.
# (4) A "color function" can be used to assign a "color" or other
# value to each polygon. If this function depends only on the
# coordinates of the polygon vertices (e.g.- 3D lighting), then
# the calculation can be done elsewhere. There are two other
# possibilities: The color function might be a cell-centered
# quantity, or a vertex-centered quantity (like the slicing
# function) on the mesh. In these cases, slice3 already
# has done much of the work, since it "knows" cell indices,
# edge interpolation coefficients, and the like.
#
# There are two particularly important cases:
# (1) Slicing function is a plane, coloring function is either a
# vertex or cell centered mesh function. Coloring function
# might also be a *function* of one or more of the predefined
# mesh functions. If you're eventually going to sweep the whole
# mesh, you want to precalculate it, otherwise on-the-fly might
# be better.
# (2) Slicing function is a vertex centered mesh function,
# coloring function is 3D shading (deferred).
#
# fslice(m3, vertex_list)
# vertex_list_iterator(m3, vertex_list, mesh3)
# fcolor(m3, vertex_list, fslice_1, fslice_2)
# the coloring function may need the value of fslice at the vertices
# in order to compute the color values by interpolation
# two "edge functions": one to detect edges where sign of fslice =
changes,
# second to interpolate for fcolor
# second to interpolate for fcolor
#
# slice3(m3, fslice, &nverts, &xyzverts, <fcolor>)
_Slice3Error =3D "Slice3Error"
def slice3 (m3, fslice, nverts, xyzverts, * args, ** kw) :
# slice3 (m3, fslice, nverts, xyzverts)
# or color_values=3D slice3(m3, fslice, nverts, xyzverts, fcolor)
# or color_values=3D slice3(m3, fslice, nverts, xyzverts, fcolor, =
1)
# slice the 3D mesh M3 using the slicing function FSLICE, returning
# the list [NVERTS, XYZVERTS, color]. Note that it is impossible to
# pass arguments as addresses, as yorick does in this routine.
# NVERTS is the number of vertices in each polygon of the slice, and
# XYZVERTS is the 3-by-sum(NVERTS) list of polygon vertices. If the
# FCOLOR argument is present, the values of that coloring function on
# the polygons are returned as the value of the slice3 function
# (numberof(color_values) =3D=3D numberof(NVERTS) =3D=3D number of =
polygons).
# If the slice function FSLICE is a function, it should be of the
# form:
# func fslice(m3, chunk)
# returning a list of function values on the specified chunk of the
# mesh m3. The format of chunk depends on the type of m3 mesh, so
# you should use only the other mesh functions xyz3 and getv3 which
# take m3 and chunk as arguments. The return value of fslice should
# have the same dimensions as the return value of getv3; the return
# value of xyz3 has an additional first dimension of length 3.
# N. B. (ZCM 2/24/97) I have eliminated the globals iso_index
# and _value, so for isosurface_slicer only, the call must be
# of the form fslice (m3, chunk, iso_index, _value).
# Likewise, I have eliminated normal and projection, so
# for plane slicer only, we do fslice (m3, chunk, normal, =
projection).
# If FSLICE is a list of 4 numbers, it is taken as a slicing plane
# with the equation FSLICE(+:1:3)*xyz(+)-FSLICE(4), as returned by
# plane3.
# If FSLICE is a single integer, the slice will be an isosurface for
# the FSLICEth variable associated with the mesh M3. In this case,
# the keyword value=3D must also be present, representing the value
# of that variable on the isosurface.
# If FCOLOR is nil, slice3 returns nil. If you want to color the
# polygons in a manner that depends only on their vertex coordinates
# (e.g.- by a 3D shading calculation), use this mode.
# If FCOLOR is a function, it should be of the form:
# func fcolor(m3, cells, l, u, fsl, fsu, ihist)
# returning a list of function values on the specified cells of the
# mesh m3. The cells argument will be the list of cell indices in
# m3 at which values are to be returned. l, u, fsl, fsu, and ihist
# are interpolation coefficients which can be used to interpolate
# from vertex centered values to the required cell centered values,
# ignoring the cells argument. See getc3 source code.
# The return values should always have dimsof(cells).
# If FCOLOR is a single integer, the slice will be an isosurface for
# the FCOLORth variable associated with the mesh M3.
# If the optional argument after FCOLOR is non-nil and non-zero,
# then the FCOLOR function is called with only two arguments:
# func fcolor(m3, cells)
# The keyword argument NODE, if present and nonzero, is a signal
# to return node-centered values rather than cell-centered
# values. (ZCM 4/16/97)
global _poly_permutations
iso_index =3D None
if type (fslice) !=3D FunctionType :
if not kw.has_key ("value") and not is_scalar (fslice) and \
len (shape (fslice)) =3D=3D 1 and len (fslice) =3D=3D 4 :
normal =3D fslice [0:3]
projection =3D fslice [3]
fslice =3D _plane_slicer
elif is_scalar (fslice) and type (fslice) =3D=3D IntType :
if not kw.has_key ("value") :
raise _Slice3Error, \
"value=3D keyword required when FSLICE is mesh variable"
_value =3D kw ["value"]
iso_index =3D fslice
fslice =3D _isosurface_slicer
else :
raise _Slice3Error, \
"illegal form of FSLICE argument, try help,slice3"
if kw.has_key ("node") :
node =3D kw ["node"]
else :
node =3D 0
# will need cell list if fcolor function to be computed
need_clist =3D len (args) > 0
if len (args) > 1 :
nointerp =3D args [1]
else :
nointerp =3D None
if need_clist :
fcolor =3D args [0]
if fcolor =3D=3D None :
need_clist =3D 0
else :
fcolor =3D None
=20
# test the different possibilities for fcolor
if need_clist and type (fcolor) !=3D FunctionType :
if not is_scalar (fcolor) or type (fcolor) !=3D IntType :
raise _Slice3Error, \
"illegal form of FCOLOR argument, try help,slice3"
# chunk up the m3 mesh and evaluate the slicing function to
# find those cells cut by fslice=3D=3D0
# chunking avoids potentially disastrously large temporaries
got_xyz =3D 0
ntotal =3D 0
# The following are used only for an irregular mesh, to
# give the sizes of each portion of the mesh.
ntotal8 =3D 0
ntotal6 =3D 0
ntotal5 =3D 0
ntotal4 =3D 0
# The following are used only for an irregular mesh, to
# give the indices of the different types of chunk in the
# results list.
i8 =3D []
i6 =3D []
i5 =3D []
i4 =3D []
itot =3D [i4, i5, i6, i8]
nchunk =3D 0
results =3D []
chunk =3D iterator3 (m3)
cell_offsets =3D [0, 0, 0, 0]
while chunk !=3D None :
# get the values of the slicing function at the vertices of
# this chunk
if fslice =3D=3D _isosurface_slicer :
fs =3D fslice (m3, chunk, iso_index, _value)
# an isosurface slicer brings back a list [vals, None]
# where vals is simply an array of the values of the
# iso_index'th function on the vertices of the specified
# chunk, or a triple, consisting of the array of
# values, an array of relative cell numbers in the
# chunk, and an offset to add to the preceding to
# get absolute cell numbers.
elif fslice =3D=3D _plane_slicer :
fs =3D fslice (m3, chunk, normal, projection)
# In the case of a plane slice, fs is a list [vals, _xyz3]
# (or [ [vals, clist, cell_offset], _xyz3] in the irregular =
case)
# where _xyz3 is the array of vertices of the chunk. _xyz3
# is ncells by 3 by something (in the irregular case),
# ncells by 3 by 2 by 2 by 2 in the regular case,
# and 3 by ni by nj by nk otherwise. vals will be
# the values of the projections of the corresponding
# vertex on the normal to the plane, positive if in
# front, and negative if in back.
else :
fs =3D fslice (m3, chunk)
if node =3D=3D 1 and fcolor !=3D None and fcolor !=3D FunctionType =
:
# need vertex-centered data
col =3D getv3 (fcolor, m3, chunk)
if type (col) =3D=3D ListType :
col =3D col [0]
else :
col =3D None
# ZCM 2/24/97 Elimination of _xyz3 as a global necessitates the =
following:
# (_xyz3 comes back as the last element of the list fs)
_xyz3 =3D fs [1]
fs =3D fs [0]
irregular =3D type (fs) =3D=3D ListType
if irregular :
cell_offset =3D fs [2]
# will need cell list if fslice did not compute xyz
got_xyz =3D _xyz3 !=3D None
need_clist =3D need_clist or not got_xyz
# If the m3 mesh is totally unstructured, the chunk should be
# arranged so that fslice returns an ncells-by-2-by-2-by-2
# (or ncells-by-3-by-2 or ncells-by-5 or ncells-by-4) array
# of vertex values of the slicing function. Note that a
# chunk of an irregular mesh always consists of just one
# kind of cell.
# On the other hand, if the mesh vertices are arranged in a
# rectangular grid (or a few patches of rectangular grids), the
# chunk should be the far less redundant rectangular patch.
if (irregular) :
# fs is a 2-sequence, of which the first element is an =
ncells-by-
# 2-by-2-by-2 (by-3-by-2, by-5, or by-4) array, and the second
# is the array of corresponding cell numbers.
# here is the fastest way to generate the required cell list
dims =3D shape (fs [0])
dim1 =3D dims [0]
slice3_precision =3D 0.0
if len (dims) =3D=3D 4 : # hex case
# Note that the sum below will be between 1 and 7
# precisely if f changes sign in the cell.
critical_cells =3D bitwise_and (add.reduce \
(reshape (ravel (transpose (less (fs [0], =
slice3_precision))), \
(8, dim1))), 7)
if (sum (critical_cells) !=3D 0) :
clist =3D take (fs [1], nonzero (critical_cells))
ntotal8 =3D ntotal8 + len (clist)
else :
clist =3D None
i8.append (len (results))
cell_offsets [3] =3D cell_offset
elif len (dims) =3D=3D 3 : # prism case
# Note that the sum below will be between 1 and 5
# precisely if f changes sign in the cell.
critical_cells =3D add.reduce \
(reshape (ravel (transpose (less (fs [0], =
slice3_precision))), \
(6, dim1)))
critical_cells =3D logical_and (greater (critical_cells, 0),
less (critical_cells, 6))
if (sum (critical_cells) !=3D 0) :
clist =3D take (fs [1], nonzero (critical_cells))
ntotal6 =3D ntotal6 + len (clist)
else :
clist =3D None
i6.append (len (results))
cell_offsets [2] =3D cell_offset
elif dims [1] =3D=3D 5 : # pyramid case
# Note that the sum below will be between 1 and 4
# precisely if f changes sign in the cell.
critical_cells =3D add.reduce \
(reshape (ravel (transpose (less (fs [0], =
slice3_precision))), \
(5, dim1)))
critical_cells =3D logical_and (greater (critical_cells, 0),
less (critical_cells, 5))
if (sum (critical_cells) !=3D 0) :
clist =3D take (fs [1], nonzero (critical_cells))
ntotal5 =3D ntotal5 + len (clist)
else :
clist =3D None
i5.append (len (results))
cell_offsets [1] =3D cell_offset
else : # tet case
# Note that the sum below will be between 1 and 3
# precisely if f changes sign in the cell.
critical_cells =3D bitwise_and (add.reduce \
(reshape (ravel (transpose (less (fs [0], =
slice3_precision))), \
(4, dim1))), 3)
if (sum (critical_cells) !=3D 0) :
clist =3D take (fs [1], nonzero (critical_cells))
ntotal4 =3D ntotal4 + len (clist)
else :
clist =3D None
i4.append (len (results))
cell_offsets [0] =3D cell_offset
else :
dims =3D shape (fs)
# fs is an ni-by-nj-by-nk array
# result of the zcen is 0, 1/8, 2/8, ..., 7/8, or 1
# slice3_precision =3D max (ravel (abs (fs))) * (-1.e-12)
slice3_precision =3D 0
clist1 =3D ravel (zcen_ (zcen_ (zcen_
(array (less (fs, slice3_precision), Float), 0), 1), 2))
clist1 =3D logical_and (less (clist1, .9), greater (clist1, =
.1))
if sum (clist1) > 0 :
clist =3D nonzero (clist1)
ntotal =3D ntotal + len (clist)
else :
clist =3D None
i8.append (len (results)) # Treat regular case as hex
if clist !=3D None :
# we need to save:
# (1) the absolute cell indices of the cells in clist
# (2) the corresponding ncells-by-2-by-2-by-2 (by-3-by-2,
# by-5, or by-4) list of fslice
# values at the vertices of these cells
if (irregular) :
# extract the portions of the data indexed by clist
fs =3D take (fs [0], clist)
if got_xyz :
_xyz3 =3D take (_xyz3, clist)
if col :
col =3D take (col, clist)
else :
# extract the to_corners portions of the data indexed by =
clist
indices =3D to_corners3 (clist, dims [1], dims [2])
no_cells =3D shape (indices) [0]
indices =3D ravel (indices)
fs =3D reshape (take (ravel (fs), indices),\
(no_cells, 2, 2, 2))
if got_xyz :
new_xyz3 =3D zeros ( (no_cells, 3, 2, 2, 2), Float )
new_xyz3 [:, 0, ...] =3D reshape (take (ravel (_xyz3 [0, =
...]),\
indices), (no_cells, 2, 2, 2))
new_xyz3 [:, 1, ...] =3D reshape (take (ravel (_xyz3 [1, =
...]),\
indices), (no_cells, 2, 2, 2))
new_xyz3 [:, 2, ...] =3D reshape (take (ravel (_xyz3 [2, =
...]),\
indices), (no_cells, 2, 2, 2))
_xyz3 =3D new_xyz3
del new_xyz3
if col !=3D None :
col =3D reshape (take (ravel (col), indices), (no_cells, =
2, 2, 2))
# NB: col represents node colors, and is only used
# if those are requested.
# here, the iterator converts to absolute cell indices without
# incrementing the chunk
if (need_clist) :
clist =3D iterator3 (m3, chunk, clist)
else :
clist =3D None
nchunk =3D nchunk + 1
need_vert_col =3D col !=3D None
results.append ( [clist, fs, _xyz3, col])
else :
results.append ( [None, None, None, None])
chunk =3D iterator3 (m3, chunk)
# endwhile chunk !=3D None
# collect the results of the chunking loop
if not ntotal and not (ntotal8 + ntotal6 + ntotal5 + ntotal4) :
return None
if ntotal : # (regular mesh, but can be handled same as hex)
ntotal8 =3D ntotal
i8 =3D range (len (results))
itot [3] =3D i8
ntot =3D [ntotal4, ntotal5, ntotal6, ntotal8]
new_results =3D []
for i in range (len (ntot)) :
# This loop processes each kind of cell independently,
# the results to be combined at the end.
if ntot [i] =3D=3D 0 : # No cells of type i
continue
if need_clist :
clist =3D zeros (ntot [i], Int)
fs =3D zeros ( (ntot [i], _no_verts [i]), Float )
if got_xyz :
xyz =3D zeros ( (ntot [i], 3, _no_verts [i]), Float )
else :
xyz =3D None
if need_vert_col :
col =3D zeros ( (ntot [i], _no_verts [i]), Float )
else :
col =3D None
k =3D 0
# collect the results of the chunking loop
for j in range (len (itot [i])) :
l =3D k
k =3D k + len (results [itot [i] [j]] [0])
if need_clist :
clist [l:k] =3D results [itot [i] [j]] [0]
fs [l:k] =3D reshape (results [itot [i] [j]] [1], (k - l, =
_no_verts [i]))
if xyz !=3D None :
xyz [l:k] =3D reshape (results [itot [i] [j]] [2],
(k - l, 3, _no_verts [i]))
if col !=3D None :
col [l:k] =3D reshape (results [itot [i] [j]] [3],
(k - l, _no_verts [i]))
if not got_xyz :
# zcm 2/4/97 go to absolute cell list again
if i > 0 and len (m3 [1]) > 2 :
adder =3D m3 [1] [3] [i - 1]
else :
adder =3D 0
xyz =3D reshape (xyz3 (m3, clist + adder), (ntot [i], 3, =
_no_verts [i]))
# produce the lists of edge intersection points
# -- generate (nsliced)x12 (9, 8, 6) array of edge mask values
# (mask non-zero if edge is cut by plane)
below =3D less (fs, 0.0)
# I put the following into C for speed
mask =3D find_mask (below, _node_edges [i])
list =3D nonzero (mask)
edges =3D array (list, copy =3D 1)
cells =3D edges / _no_edges [i]
edges =3D edges % _no_edges [i]
# construct edge endpoint indices in fs, xyz arrays
# the numbers are the endpoint indices corresponding to
# the order of the _no_edges [i] edges in the mask array
lower =3D take (_lower_vert [i], edges) + _no_verts [i] * cells
upper =3D take (_upper_vert [i], edges) + _no_verts [i] * cells
fsl =3D take (ravel (fs), lower)
fsu =3D take (ravel (fs), upper)
# following denominator guaranteed non-zero
denom =3D fsu - fsl
fsu =3D fsu / denom
fsl =3D fsl / denom
new_xyz =3D zeros ( (len (lower), 3), Float )
new_xyz [:, 0] =3D reshape ( (take (ravel (xyz [:, 0]), lower) * =
fsu - \
take (ravel (xyz [:, 0]), upper) * fsl), (len (lower),))
new_xyz [:, 1] =3D reshape ( (take (ravel (xyz [:, 1]), lower) * =
fsu - \
take (ravel (xyz [:, 1]), upper) * fsl), (len (lower),))
new_xyz [:, 2] =3D reshape ( (take (ravel (xyz [:, 2]), lower) * =
fsu - \
take (ravel (xyz [:, 2]), upper) * fsl), (len (lower),))
xyz =3D new_xyz
del new_xyz
if col !=3D None :
# Extract subset of the data the same way
col =3D take (ravel (col), lower) * fsu - \
take (ravel (col), upper) * fsl
# The xyz array is now the output xyzverts array,
# but for the order of the points within each cell.
# give each sliced cell a "pattern index" between 0 and 255
# (non-inclusive) representing the pattern of its 8 corners
# above and below the slicing plane
p2 =3D left_shift (ones (_no_verts [i], Int) , array (
[0, 1, 2, 3, 4, 5, 6, 7], Int) [0: _no_verts [i]])
pattern =3D transpose (sum (transpose (multiply (below, p2))))
# broadcast the cell's pattern onto each of its sliced edges
pattern =3D take (pattern, list / _no_edges [i])
# Let ne represent the number of edges of this type of cell,
# and nv the number of vertices.
# To each pattern, there corresponds a permutation of the
# twelve edges so that they occur in the order in which the
# edges are to be connected. Let each such permuation be
# stored as a list of integers from 0 to ne - 1 such that
# sorting the integers into increasing order rearranges the edges =
at
# the corresponding indices into the correct order. (The position
# of unsliced edges in the list is arbitrary as long as the sliced
# edges are in the proper order relative to each other.)
# Let these permutations be stored in a ne-by-2**nv - 2 array
# _poly_permutations (see next comment for explanation of 4 * ne):
pattern =3D take (ravel (transpose (_poly_permutations [i])),=20
_no_edges [i] * (pattern - 1) + edges) + 4 * _no_edges [i] * =
cells
order =3D argsort (pattern)
xyz1 =3D zeros ( (len (order), 3), Float )
xyz1 [:,0] =3D take (ravel (xyz [:,0]), order)
xyz1 [:,1] =3D take (ravel (xyz [:,1]), order)
xyz1 [:,2] =3D take (ravel (xyz [:,2]), order)
xyz =3D xyz1
if col !=3D None :
col =3D take (col, order)
edges =3D take (edges, order)
pattern =3D take (pattern, order)
# cells(order) is same as cells by construction */
# There remains only the question of splitting the points in
# a single cell into multiple disjoint polygons.
# To do this, we need one more precomputed array: poly_splits
# should be another ne-by-2**nv - 2 array with values between 0 =
and 3
# 0 for each edge on the first part, 1 for each edge on the
# second part, and so on up to 3 for each edge on the fourth
# part. The value on unsliced edges can be anything, say 0.
# With a little cleverness poly_splits can be combined with
# _poly_permutations, by putting _poly_permutations =3D
# _poly_permutations(as described above) + _no_edges =
[i]*poly_splits
# (this doesn't change the ordering of _poly_permutations).
# I assume this has been done here:
pattern =3D pattern / _no_edges [i]
# now pattern jumps by 4 between cells, smaller jumps within cells
# get the list of places where a new value begins, and form a
# new pattern with values that increment by 1 between each plateau
pattern =3D dif_ (pattern, 0)
nz =3D nonzero (pattern)
list =3D zeros (len (nz) + 1, Int)
list [1:] =3D nz + 1
newpat =3D zeros (len (pattern) + 1, Int)
newpat [0] =3D 1
newpat [1:] =3D cumsum (not_equal (pattern, 0)) + 1
pattern =3D newpat
nverts =3D histogram (pattern) [1:]
xyzverts =3D xyz
# finally, deal with any fcolor function
if fcolor =3D=3D None :
new_results.append ( [nverts, xyzverts, None])
continue
# if some polys have been split, need to split clist as well
if len (list) > len (clist) :
clist =3D take (clist, take (cells, list))
if col =3D=3D None :
if nointerp =3D=3D None :
if type (fcolor) =3D=3D FunctionType :
col =3D fcolor (m3, clist + cell_offsets [i], lower, =
upper, fsl,
fsu, pattern - 1)
else :
col =3D getc3 (fcolor, m3, clist + cell_offsets [i], =
lower, upper,
fsl, fsu, pattern - 1)
else :
if type (fcolor) =3D=3D FunctionType :
col =3D fcolor (m3, clist + cell_offsets [i])
else :
col =3D getc3 (fcolor, m3, clist + cell_offsets [i])
new_results.append ( [nverts, xyzverts, col])
# New loop to consolidate the return values
nv_n =3D 0
xyzv_n =3D 0
col_n =3D 0
for i in range (len (new_results)) :
nv_n =3D nv_n + len (new_results [i] [0])
xyzv_n =3D xyzv_n + shape (new_results [i] [1]) [0]
if new_results [i] [2] !=3D None :
col_n =3D col_n + len (new_results [i] [2])
nverts =3D zeros (nv_n, Int)
xyzverts =3D zeros ( (xyzv_n, 3), Float )
if col_n !=3D 0 :
col =3D zeros (col_n, Float )
else :
col =3D None
nv_n1 =3D 0
xyzv_n1 =3D 0
col_n1 =3D 0
for i in range (len (new_results)) :
nv_n2 =3D len (new_results [i] [0])
xyzv_n2 =3D shape (new_results [i] [1]) [0]
nverts [nv_n1:nv_n1 + nv_n2] =3D new_results [i] [0]
xyzverts [xyzv_n1:xyzv_n1 + xyzv_n2] =3D new_results [i] [1]
if new_results [i] [2] !=3D None :
col_n2 =3D len (new_results [i] [2])
col [col_n1:col_n1 + col_n2] =3D new_results [i] [2]
col_n1 =3D col_n1 + col_n2
nv_n1 =3D nv_n1 + nv_n2
xyzv_n1 =3D xyzv_n1 + xyzv_n2
return [nverts, xyzverts, col]
# The iterator3 function combines three distinct operations:
# (1) If only the M3 argument is given, return the initial
# chunk of the mesh. The chunk will be no more than
# _chunk3_limit cells of the mesh.
# (2) If only M3 and CHUNK are given, return the next CHUNK,
# or [] if there are no more chunks.
# (3) If M3, CHUNK, and CLIST are all specified, return the
# absolute cell index list corresponding to the index list
# CLIST of the cells in the CHUNK.
# Do not increment the chunk in this case.
#
# The form of the CHUNK argument and return value for cases (1)
# and (2) is not specified, but it must be recognized by the
# xyz3 and getv3 functions which go along with this iterator3.
# (For case (3), CLIST and the return value are both ordinary
# index lists.)
_Slice3MeshError =3D "Slice3MeshError"
def slice3mesh (xyz, * args, ** kw) :
# slice3mesh (z, color =3D None, smooth =3D 0)
# slice3mesh (nxny, dxdy, x0y0, z, color =3D None, smooth =3D 0)
# slice3mesh (x, y, z, color =3D None, smooth =3D 0)
#
# slice3mesh returns a triple [nverts, xyzverts, color]
# nverts is no_cells long and the ith entry tells how many
# vertices the ith cell has.
# xyzverts is sum (nverts) by 3 and gives the vertex
# coordinates of the cells in order.
# color, if present, is len (nverts) long and contains
# a color value for each cell in the mesh if smooth =3D=3D 0;
# sum (nverts) long and contains a color value for each
# node in the mesh if smooth =3D=3D 1.
# There are a number of ways to call slice3mesh:
# slice3mesh (z, color =3D None, smooth =3D 0)
# z is a two dimensional array of function values, assumed
# to be on a uniform mesh nx by ny nodes (assuming z is nx by ny)
# nx being the number of nodes in the x direction, ny the number
# in the y direction.
# color, if specified, is either an nx - 1 by ny - 1 array
# of cell-centered values by which the surface is to
# be colored, or an nx by ny array of vertex-
# centered values, which will be averaged over each
# cell to give cell-centered values if smooth =3D=3D 0, or
# returned as a node-centered array sum (nverts) long if
# smooth =3D=3D 1.
# slice3mesh (nxny, dxdy, x0y0, z, color =3D None, smooth =3D 0)
# In this case, slice3mesh accepts the specification for
# a regular 2d mesh: nxny is the number of cells in the
# x direction and the y direction (i. e., its two components
# are nx - 1 and ny - 1, nx by ny being the node size;
# x0y0 are the initial
# values of x and y; and dxdy are the increments in the
# two directions. z is the height of a surface above
# the xy plane and must be dimensioned nx by ny.=20
# color, if specified, is as above.
# slice3mesh (x, y, z, color =3D None, smooth =3D 0)
# z is as above, an nx by ny array of function values
# on a mesh of the same dimensions. There are two choices
# for x and y: they can both be one-dimensional, dimensioned
# nx and ny respectively, in which case they represent a
# mesh whose edges are parallel to the axes; or else they
# can both be nx by ny, in which case they represent a
# general quadrilateral mesh.
# color, if specified, is as above.
two_d =3D 0
if kw.has_key ("smooth") :
smooth =3D kw ["smooth"]
else :
smooth =3D 0
if len (args) =3D=3D 0 :
# Only the z argument is present
if len (shape (xyz)) !=3D 2 :
raise _Slice3MeshError, \
"z must be two dimensional."
else :
z =3D xyz
ncx =3D shape (xyz) [0]
ncy =3D shape (xyz) [1]
x =3D arange (ncx, typecode =3D Float )
y =3D arange (ncy, typecode =3D Float )
elif len (args) =3D=3D 3 :
# must be the (nxny, dxdy, x0y0, z...) form
ncx =3D xyz [0] + 1
ncy =3D xyz [1] + 1
x =3D arange (ncx, typecode =3D Float ) * args [0] [0] + args [1] =
[0]
y =3D arange (ncy, typecode =3D Float ) * args [0] [1] + args [1] =
[1]
z =3D args [2]
if (ncx, ncy) !=3D shape (z) :
raise _Slice3MeshError, \
"The shape of z must match the shape of x and y."
elif len (args) =3D=3D 2 :
# must be the x, y, z format
x =3D xyz
y =3D args [0]
z =3D args [1]
dims =3D shape (x)
if len (dims) =3D=3D 2 :
two_d =3D 1
if dims !=3D shape (y) or dims !=3D shape (z) :
raise _Slice3MeshError, \
"The shapes of x, y, and z must match."
ncx =3D dims [0]
ncy =3D dims [1]
elif len (dims) =3D=3D 1 :
ncx =3D dims [0]
ncy =3D len (y)
if (ncx, ncy) !=3D shape (z) :
raise _Slice3MeshError, \
"The shape of z must match the shape of x and y."
else :
raise _Slice3MeshError, \
"Unable to decipher arguments to slice3mesh."
else :
raise _Slice3MeshError, \
"Unable to decipher arguments to slice3mesh."
nverts =3D ones ( (ncx - 1) * (ncy - 1), Int) * 4
ncxx =3D arange (ncx - 1, typecode =3D Int) * (ncy)
ncyy =3D arange (ncy - 1, typecode =3D Int)
if kw.has_key ("color") :
color =3D kw ["color"]
else :
color =3D None
if color !=3D None :
# col =3D array (len (nverts), Float )
if shape (color) =3D=3D (ncx - 1, ncy - 1) :
col =3D color
elif shape (color) =3D=3D (ncx, ncy) and smooth =3D=3D 0 :
col =3D ravel (color)
# Lower left, upper left, upper right, lower right
col =3D 0.25 * (take (col, ravel (add.outer ( ncxx, ncyy))) +
take (col, ravel (add.outer ( ncxx, ncyy + 1))) +
take (col, ravel (add.outer ( ncxx + ncy, ncyy + 1))) +
take (col, ravel (add.outer ( ncxx + ncy, ncyy))))
elif shape (color) =3D=3D (ncx, ncy) and smooth !=3D 0 :
# Node-centered colors are wanted (smooth plots)
col =3D ravel (color)
col =3D ravel (transpose (array ( [
take (col, ravel (add.outer ( ncxx, ncyy))),
take (col, ravel (add.outer ( ncxx, ncyy + 1))),
take (col, ravel (add.outer ( ncxx + ncy, ncyy + 1))),
take (col, ravel (add.outer ( ncxx + ncy, ncyy)))])))
else :
raise _Slice3MeshError, \
"color must be cell-centered or vertex centered."
else :
col =3D None
xyzverts =3D zeros ( (4 * (ncx -1) * (ncy -1), 3), Float )
if not two_d :
x1 =3D multiply.outer (ones (ncy - 1, Float), x [0:ncx - 1])
x2 =3D multiply.outer (ones (ncy - 1, Float), x [1:ncx])
xyzverts [:, 0] =3D ravel (transpose (array ([x1, x1, x2, x2])))
del x1, x2
y1 =3D multiply.outer (y [0:ncy - 1], ones (ncx - 1))
y2 =3D multiply.outer (y [1:ncy], ones (ncx - 1))
xyzverts [:, 1] =3D ravel (transpose (array ([y1, y2, y2, y1])))
del y1, y2
else :
newx =3D ravel (x)
xyzverts [:, 0] =3D ravel (transpose (array ( [
take (newx, ravel (add.outer ( ncxx, ncyy))),
take (newx, ravel (add.outer ( ncxx, ncyy + 1))),
take (newx, ravel (add.outer ( ncxx + ncy, ncyy + 1))),
take (newx, ravel (add.outer ( ncxx + ncy, ncyy)))])))
newy =3D ravel (y)
xyzverts [:, 1] =3D ravel (transpose (array ( [
take (newy, ravel (add.outer ( ncxx, ncyy))),
take (newy, ravel (add.outer ( ncxx, ncyy + 1))),
take (newy, ravel (add.outer ( ncxx + ncy, ncyy + 1))),
take (newy, ravel (add.outer ( ncxx + ncy, ncyy)))])))
newz =3D ravel (z)
xyzverts [:, 2] =3D ravel (transpose (array ( [
take (newz, ravel (add.outer ( ncxx, ncyy))),
take (newz, ravel (add.outer ( ncxx, ncyy + 1))),
take (newz, ravel (add.outer ( ncxx + ncy, ncyy + 1))),
take (newz, ravel (add.outer ( ncxx + ncy, ncyy)))])))
=20
return [nverts, xyzverts, col]
=20
def iterator3 (m3 , chunk =3D None, clist =3D None) :
return m3 [0] [3] (m3, chunk, clist)
# biggest temporary is 3 doubles times this,
# perhaps 4 or 5 doubles times this is most at one time
_chunk3_limit =3D 10000
def iterator3_rect (m3, chunk, clist) :
# Note: if you look at the yorick version of this routine, you
# will see that the significance of the subscripts is reversed.
# This is because we do things in row-major order.
global _chunk3_limit
=20
if chunk =3D=3D None :
dims =3D m3 [1] [0] # [ni,nj,nk] cell dimensions
[ni, nj, nk] =3D [dims [0], dims [1], dims [2]]
njnk =3D nj * nk
if _chunk3_limit <=3D nk :
# stuck with 1D chunks
ck =3D (nk - 1) / _chunk3_limit + 1
cj =3D ci =3D 0
elif _chunk3_limit <=3D njnk :
# 2D chunks
ci =3D ck =3D 0
cj =3D (njnk - 1) / _chunk3_limit + 1
else :
# 3D chunks
cj =3D ck =3D 0
ci =3D (njnk * ni - 1) / _chunk3_limit + 1
chunk =3D array ( [[ci =3D=3D 0, cj =3D=3D 0, ck =3D=3D 0],
[not ci, nj * (ci !=3D 0) + (ck !=3D 0),
nk * ( (cj + ci) !=3D 0)],
[ci, cj, ck], [ni, nj, nk]])
else :
ni =3D chunk [3,0]
nj =3D chunk [3,1]
nk =3D chunk [3,2]
njnk =3D nj * nk
offsets =3D array ( [njnk, nj, 1], Int)
if clist !=3D None :
# add offset for this chunk to clist and return
return sum (offsets * ( chunk [0] - 1)) + clist
# increment to next chunk
xi =3D chunk [1, 0]
xj =3D chunk [1, 1]
xk =3D chunk [1, 2]
np =3D chunk [2, 2]
if (np) :
# 1D chunks
if xk =3D=3D nk :
if xj =3D=3D nj :
if xi =3D=3D ni : return None
xi =3D xi + 1
xj =3D 1;
else :
xj =3D xj + 1
xk =3D 0
ck =3D xk + 1
step =3D ck / np=20
frst =3D ck % np # first frst steps are step+1
if (xk < (step + 1) * frst) : step =3D step + 1
xk =3D xk + step
chunk [0] =3D array ( [xi, xj, ck])
chunk [1] =3D array ( [xi, xj, xk])
else :
np =3D chunk [2, 1]
if (np) :
if (xj =3D=3D nj) :
if (xi =3D=3D ni) : return None
xi =3D xi + 1
xj =3D 0
cj =3D xj + 1
step =3D nj / np
frst =3D nj % np # first frst steps are step+1
if (xj < (step + 1) * frst) : step =3D step + 1
xj =3D xj + step
chunk [0, 0:2] =3D array ( [xi, cj])
chunk [1, 0:2] =3D array ( [xi, xj])
else :
if xi =3D=3D ni : return None
ci =3D xi + 1
np =3D chunk [2, 0]
step =3D ni / np
frst =3D ni % np # first frst steps are step+1
if (xi < (step + 1) * frst) : step =3D step + 1
xi =3D xi + step
chunk [0, 0] =3D ci
chunk [1, 0] =3D xi
return chunk
def iterator3_irreg (m3, chunk, clist) :
# Does the same thing as iterator3_rect only for an irregular
# rectangular mesh. It simply splits a large mesh into smaller
# parts. Whether this is necessary I am not sure.
# Certainly it makes it easier in the irregular case to handle
# the four different types of cells separately.
# if clist is present, in the irregular case it is already
# the list of absolute cell indices, so it is simply returned.
# This and other routines to do with irregular meshes return a
# chunk which is a 2-list. The first item delimits the chunk;
# the second gives a list of corresponding cell numbers.
global _chunk3_limit
if clist !=3D None:
return clist
dims =3D m3 [1] [0] # ncells by _no_verts array of subscripts
# (or a list of from one to four of same)
if type (dims) !=3D ListType :
if chunk =3D=3D None: # get the first chunk
return [ [0, min (shape (dims) [0], _chunk3_limit)],
arange (0, min (shape (dims) [0], _chunk3_limit),
typecode =3D Int)]
else : # iterate to next chunk
start =3D chunk [0] [1]
if start >=3D shape(dims) [0] :
return None
else :
return [ [start, min (shape (dims) [0], start + =
_chunk3_limit)],
arange (start, min (shape (dims) [0],
start + _chunk3_limit),
typecode =3D Int)]
else :
totals =3D m3 [1] [3] # cumulative totals of numbers of cells
if chunk =3D=3D None :
return [ [0, min (totals [0], _chunk3_limit)],
arange (0, min (totals [0], _chunk3_limit),
typecode =3D Int)]
else : # iterate to next chunk
start =3D chunk [0] [1]
if start >=3D totals [-1] :
return None
else :
for i in range (len (totals)) :
if start < totals [i] :
break
return [ [start, min (totals [i], start + _chunk3_limit)],
arange (start,
min (totals [i], start + _chunk3_limit),=20
typecode =3D Int)]
def getv3 (i, m3, chunk) :
# getv3(i, m3, chunk)
# return vertex values of the Ith function attached to 3D mesh M3
# for cells in the specified CHUNK. The CHUNK may be a list of
# cell indices, in which case getv3 returns a 2x2x2x(dimsof(CHUNK))
# list of vertex coordinates. CHUNK may also be a mesh-specific data
# structure used in the slice3 routine, in which case getv3 may
# return a (ni)x(nj)x(nk) array of vertex values. For meshes which
# are logically rectangular or consist of several rectangular
# patches, this is up to 8 times less data, with a concomitant
# performance advantage. Use getv3 when writing slicing functions
# for slice3.
return m3 [0] [1] (i, m3, chunk)
_Getv3Error =3D "Getv3Error"
def getv3_rect (i, m3, chunk) :
fi =3D m3 [2]
i =3D i - 1
if i < 0 or is_scalar (fi) or i >=3D len (fi) :
raise _Getv3Error, "no such mesh function as F" + `i`
dims =3D m3 [1] [0]
if dims =3D=3D shape (fi [i]) :
raise _Getv3Error, "mesh function F" + `i` + " is not =
vertex-centered"
if len (shape (chunk)) !=3D 1 :
c =3D chunk
# The difference here is that our arrays are 0-based, while
# yorick's are 1-based; and the last element in a range is not
# included in the result array.
return fi [i] [c [0, 0] - 1:1 + c [1, 0], c [0, 1] - 1:1 + c [1, =
1] ,
c [0, 2] - 1:1 + c [1, 2]]
else :
# Need to create an array of fi values the same size and shape
# as what to_corners3 returns.
# To avoid exceedingly arcane calculations attempting to
# go backwards to a cell list, this branch returns the list
# [<function values>, chunk]
# Then it is trivial for slice3 to find a list of cell
# numbers in which fi changes sign.
indices =3D to_corners3 (chunk, dims [0] + 1, dims [1] + 1)
no_cells =3D shape (indices) [0]
indices =3D ravel (indices)
retval =3D reshape (take (ravel (fi [i]), indices), (no_cells, 2, =
2, 2))
return [retval, chunk]
def getv3_irreg (i, m3, chunk) :
# for an irregular mesh, returns a 3-list whose elements are:
# (1) the function values for the ith function on the vertices of the
# given chunk. (The function values must have the same dimension
# as the coordinates; there is no attempt to convert zone-centered
# values to vertex-centered values.)
# (2) an array of relative cell numbers within the list of cells
# of this type.
# (3) a number that can be added to these relative numbers to gives
# the absolute cell numbers for correct access to their coordinates
# and function values.
=20
fi =3D m3 [2]
i =3D i - 1
if i < 0 or is_scalar (fi) or i >=3D len (fi) :
raise _Getv3Error, "no such function as F" + `i`
# len (fi [i]) and the second dimension of m3 [1] [1] (xyz) should
# be the same, i. e., there is a value associated with each =
coordinate.
if len (fi [i]) !=3D len (m3 [1] [1] [0]) :
raise _Getv3Error, "mesh function F" + `i` + " is not =
vertex-centered."
verts =3D m3 [1] [0]
oldstart =3D chunk [0] [0]
oldfin =3D chunk [0] [1]
no_cells =3D oldfin - oldstart
if type (verts) !=3D ListType : # Only one kind of cell in mesh
indices =3D ravel (verts [oldstart:oldfin])
else : # A list of possibly more than one kind
sizes =3D m3 [1] [2]
totals =3D m3 [1] [3]
for j in range (len (totals)) :
if oldfin <=3D totals [j] :
break
verts =3D verts [j]
if j > 0 :
start =3D oldstart - totals [j - 1]
fin =3D oldfin - totals [j - 1]
else :
start =3D oldstart=20
fin =3D oldfin
indices =3D ravel (verts [start:fin])
tc =3D shape (verts) [1]
# ZCM 2/4/97 the array of cell numbers must be relative
if tc =3D=3D 8 : # hex cells
return [ reshape (take (fi [i], indices), (no_cells, 2, 2, 2)),
arange (0, no_cells, typecode =3D Int), oldstart]
elif tc =3D=3D 6 : # pyramids
return [ reshape (take (fi [i], indices), (no_cells, 3, 2)),
arange (0, no_cells, typecode =3D Int), oldstart]
else : # tetrahedron or pyramid
return [ reshape (take (fi [i], indices), (no_cells, tc)),
arange (0, no_cells, typecode =3D Int), oldstart]
_Getc3Error =3D "Getc3Error"
def getc3 (i, m3, chunk, *args) :
# getc3(i, m3, chunk)
# or getc3(i, m3, clist, l, u, fsl, fsu, cells)
# return cell values of the Ith function attached to 3D mesh M3
# for cells in the specified CHUNK. The CHUNK may be a list of
# cell indices, in which case getc3 returns a (dimsof(CHUNK))
# list of vertex coordinates. CHUNK may also be a mesh-specific data
# structure used in the slice3 routine, in which case getc3 may
# return a (ni)x(nj)x(nk) array of vertex values. There is no
# savings in the amount of data for such a CHUNK, but the gather
# operation is cheaper than a general list of cell indices.
# Use getc3 when writing colorng functions for slice3.
# If CHUNK is a CLIST, the additional arguments L, U, FSL, and FSU
# are vertex index lists which override the CLIST if the Ith attached
# function is defined on mesh vertices. L and U are index lists into
# the (dimsof(CLIST))x2x2x2 vertex value array, say vva, and FSL
# and FSU are corresponding interpolation coefficients; the zone
# centered value is computed as a weighted average of involving these
# coefficients. The CELLS argument is required by histogram to do
# the averaging. See the source code for details.
# By default, this conversion (if necessary) is done by averaging
# the eight vertex-centered values.
if len (args) =3D=3D 0 :
l =3D None
u =3D None
fsl =3D None
fsu =3D None
cells =3D None
elif len (args) =3D=3D 5 :
l =3D args [0]
u =3D args [1]
fsl =3D args [2]
fsu =3D args [3]
cells =3D args [4]
else :
raise _Getc3Error, "getc3 requires either three or eight =
parameters."
return m3 [0] [2] (i, m3, chunk, l, u, fsl, fsu, cells)
def getc3_rect (i, m3, chunk, l, u, fsl, fsu, cells) :
fi =3D m3 [2]
m3 =3D m3 [1]
if ( i < 1 or i > len (fi)) :
raise _Getc3Error, "no such mesh function as F" + `i - 1`
dims =3D m3 [0]
if shape (fi [i - 1]) =3D=3D dims :
# it is a cell-centered quantity
if len (shape (chunk)) !=3D 1 :
c =3D chunk
# The difference here is that our arrays are 0-based, while
# yorick's are 1-based; and the last element in a range is not
# included in the result array.
return fi [i - 1] [c [0, 0] - 1:1 + c [1, 0],
c [0, 1] - 1:1 + c [1, 1] ,
c [0, 2] - 1:1 + c [1, 2]]
else :
[k, l. m] =3D dims
return reshape (take (ravel (fi [i - 1]), chunk),
(len (chunk), k, l, m))
else :
# it is vertex-centered, so we take averages to get cell quantity
if len (shape (chunk)) !=3D 1 :
c =3D chunk
# The difference here is that our arrays are 0-based, while
# yorick's are 1-based; and the last element in a range is not
# included in the result array.
return zcen_ (zcen_( zcen_ (
(fi [i - 1] [c [0, 0] - 1:1 + c [1, 0],
c [0, 1] - 1:1 + c [1, 1] ,
c [0, 2] - 1:1 + c [1, 2]]), 0), 1), 2)
else :
indices =3D to_corners3 (chunk, dims [1] + 1, dims [2] + 1)
no_cells =3D shape (indices) [0]
indices =3D ravel (indices)
corners =3D take (ravel (fi [i - 1]), indices)
if l =3D=3D None :
return 0.125 * sum (transpose (reshape (corners, (no_cells, =
8))))
else :
# interpolate corner values to get edge values
corners =3D (take (corners, l) * fsu -
take (corners, u) * fsl) / (fsu -fsl)
# average edge values (vertex values of polys) on each poly
return histogram (cells, corners) / histogram (cells)
def getc3_irreg (i, m3, chunk, l, u, fsl, fsu, cells) :
# Same thing as getc3_rect, i. e., returns the same type of
# data structure, but from an irregular rectangular mesh.
# m3 [1] is a 2-list; m3[1] [0] is an array whose ith element
# is an array of coordinate indices for the ith cell,
# or a list of up to four such arrays.
# m3 [1] [1] is the 3 by nverts array of coordinates.
# m3 [2] is a list of arrays of vertex-centered or cell-centered
# data.
# chunk may be a list, in which case chunk [0] is a 2-sequence
# representing a range of cell indices; or it may be a =
one-dimensional
# array, in which case it is a nonconsecutive set of cell indices.
# It is guaranteed that all cells indexed by the chunk are the
# same type.
fi =3D m3 [2]
if i < 1 or i > len (fi) :
raise _Getc3Error, "no such mesh function as F" + `i - 1`
verts =3D m3 [1] [0]
if type (verts) =3D=3D ListType :
sizes =3D m3 [1] [2]
totals =3D m3 [1] [3]
if type (verts) =3D=3D ListType and totals [-1] =3D=3D len (fi [i - =
1]) or \
type (verts) !=3D ListType and shape (verts) [0] =3D=3D len (fi [i =
- 1]) :
# cell-centered case
if type (chunk) =3D=3D ListType :
return fi [i - 1] [chunk [0] [0]:chunk [0] [1]]
elif type (chunk) =3D=3D ArrayType and len (shape (chunk)) =3D=3D =
1 :
return take (fi [i - 1], chunk)
else :
raise _Getc3Error, "chunk argument is incomprehensible."
if len (fi [i - 1]) !=3D shape (m3 [1] [1]) [1] :
raise _Getc3Error, "F" + `i - 1` + " has the wrong size to be " \
"either zone-centered or node-centered."
# vertex-centered case
# First we need to pick up the vertex subscripts, which are
# also the fi [i - 1] subscripts.
if type (verts) !=3D ListType :
if type (chunk) =3D=3D ListType :
indices =3D verts [chunk [0] [0]:chunk [0] [1]]
elif type (chunk) =3D=3D ArrayType and len (shape (chunk)) =3D=3D =
1 :
indices =3D take (verts, chunk)
else :
raise _Getc3Error, "chunk argument is incomprehensible."
else :
# We have a list of vertex subscripts, each for a different
# type of cell; need to extract the correct list:
if type (chunk) =3D=3D ListType :
start =3D chunk [0] [0]
fin =3D chunk [0] [1]
for j in range (len (totals)) :
if fin <=3D totals [j] :
break
verts =3D verts [j]
if j > 0 :
start =3D start - totals [j - 1]
fin =3D fin - totals [j - 1]
indices =3D verts [start:fin]
elif type (chunk) =3D=3D ArrayType and len (shape (chunk)) =3D=3D =
1 :
for j in range (len (totals)) :
if chunk [-1] <=3D totals [j] :
break
verts =3D verts [j]
ch =3D chunk
if j > 0 :
ch =3D chunk - totals [j - 1]
indices =3D take (verts, ch)
else :
raise _Getc3Error, "chunk argument is incomprehensible."
shp =3D shape (indices)
no_cells =3D shp [0]
indices =3D ravel (indices)
corners =3D take (fi [i - 1], indices)
if l =3D=3D None :
return (1. / shp [1]) * transpose ((sum (transpose (reshape =
(corners,
(no_cells, shp [1]))) [0:shp [1]])))
else :
# interpolate corner values to get edge values
corners =3D (take (corners, l) * fsu -
take (corners, u) * fsl) / (fsu -fsl)
# average edge values (vertex values of polys) on each poly
return histogram (cells, corners) / histogram (cells)
_no_verts =3D array ( [4, 5, 6, 8])
_no_edges =3D array ( [6, 8, 9, 12])
# Lower and upper vertex subscripts for each edge
_lower_vert4 =3D array ( [0, 0, 0, 1, 2, 3], Int)
_lower_vert5 =3D array ( [0, 0, 0, 0, 1, 2, 3, 4], Int)
_lower_vert6 =3D array ( [0, 1, 0, 1, 2, 3, 0, 2, 4], Int)
_lower_vert8 =3D array ( [0, 1, 2, 3, 0, 1, 4, 5, 0, 2, 4, 6], Int)
_lower_vert =3D [_lower_vert4, _lower_vert5, _lower_vert6, _lower_vert8]
_upper_vert4 =3D array ( [1, 2, 3, 2, 3, 1], Int)
_upper_vert5 =3D array ( [1, 2, 3, 4, 2, 3, 4, 1], Int)
_upper_vert6 =3D array ( [4, 5, 2, 3, 4, 5, 1, 3, 5], Int)
_upper_vert8 =3D array ( [4, 5, 6, 7, 2, 3, 6, 7, 1, 3, 5, 7], Int)
_upper_vert =3D [_upper_vert4, _upper_vert5, _upper_vert6, _upper_vert8]
_node_edges8_s =3D array ( [ [1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0],
[0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0],
[0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0],
[1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
[0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1],
[0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1]], Int)
_node_edges8 =3D array ( [ [0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1],
[0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1],
[0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
[1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0],
[0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0],
[0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0],
[0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0],
[1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0]], Int)
_node_edges6_s =3D array ( [ [1, 0, 1, 0, 0, 0, 1, 0, 0],
[0, 1, 0, 1, 0, 0, 1, 0, 0],
[0, 0, 1, 0, 1, 0, 0, 1, 0],
[0, 0, 0, 1, 0, 1, 0, 1, 0],
[1, 0, 0, 0, 1, 0, 0, 0, 1],
[0, 1, 0, 0, 0, 1, 0, 0, 1]], Int)
_node_edges6 =3D array ( [ [0, 1, 0, 0, 0, 1, 0, 0, 1],
[1, 0, 0, 0, 1, 0, 0, 0, 1],
[0, 0, 0, 1, 0, 1, 0, 1, 0],
[0, 0, 1, 0, 1, 0, 0, 1, 0],
[0, 1, 0, 1, 0, 0, 1, 0, 0],
[1, 0, 1, 0, 0, 0, 1, 0, 0]], Int)
_node_edges4_s =3D array ( [ [1, 1, 1, 0, 0, 0],
[1, 0, 0, 1, 0, 1],
[0, 1, 0, 1, 1, 0],
[0, 0, 1, 0, 1, 1]], Int)
_node_edges4 =3D array ( [ [0, 0, 1, 0, 1, 1],
[0, 1, 0, 1, 1, 0],
[1, 0, 0, 1, 0, 1],
[1, 1, 1, 0, 0, 0]], Int)
_node_edges5_s =3D array ( [ [1, 1, 1, 1, 0, 0, 0, 0],
[1, 0, 0, 0, 1, 0, 0, 1],
[0, 1, 0, 0, 1, 1, 0, 0],
[0, 0, 1, 0, 0, 1, 1, 0],
[0, 0, 0, 1, 0, 0, 1, 1]], Int)
_node_edges5 =3D array ( [ [0, 0, 0, 1, 0, 0, 1, 1],
[0, 0, 1, 0, 0, 1, 1, 0],
[0, 1, 0, 0, 1, 1, 0, 0],
[1, 0, 0, 0, 1, 0, 0, 1],
[1, 1, 1, 1, 0, 0, 0, 0]], Int)
_node_edges =3D [_node_edges4_s, _node_edges5_s, _node_edges6_s, =
_node_edges8_s]
_node_edges3 =3D [_node_edges4, _node_edges5, _node_edges6, =
_node_edges8]
def _construct3 (itype) :
global _node_edges
global _no_verts
global _no_edges
i =3D arange (1, 2**_no_verts [itype] - 1, typecode =3D Int)
if itype =3D=3D 0 :
below =3D transpose (not_equal (array ( [bitwise_and (i, 8),
bitwise_and (i, 4),
bitwise_and (i, 2),
bitwise_and (i, 1)]), 0))
elif itype =3D=3D 1 :
below =3D transpose (not_equal (array ( [bitwise_and (i, 16),
bitwise_and (i, 8),
bitwise_and (i, 4),
bitwise_and (i, 2),
bitwise_and (i, 1)]), 0))
elif itype =3D=3D 2 :
below =3D transpose (not_equal (array ( [bitwise_and (i, 32),
bitwise_and (i, 16),
bitwise_and (i, 8),
bitwise_and (i, 4),
bitwise_and (i, 2),
bitwise_and (i, 1)]), 0))
elif itype =3D=3D 3 :
below =3D transpose (not_equal (array ( [bitwise_and (i, 128),
bitwise_and (i, 64),
bitwise_and (i, 32),
bitwise_and (i, 16),
bitwise_and (i, 8),
bitwise_and (i, 4),
bitwise_and (i, 2),
bitwise_and (i, 1)]), 0))
# For some reason the node edges for a cell need to be in different =
order
# here than in slice3 to get the correct results. Hence _node_edges3.
mask =3D find_mask (below, _node_edges3 [itype])
return construct3 (mask, itype)
# =
------------------------------------------------------------------------
_poly_permutations4 =3D _construct3 (0)
_poly_permutations5 =3D _construct3 (1)
_poly_permutations6 =3D _construct3 (2)
_poly_permutations8 =3D _construct3 (3)
_poly_permutations =3D [_poly_permutations4, _poly_permutations5,=20
_poly_permutations6, _poly_permutations8]
_ContourError =3D "ContourError"
# =
------------------------------------------------------------------------
def plzcont (nverts, xyzverts, contours =3D 8, scale =3D "lin", clear =
=3D 1,
edges =3D 0, color =3D None, cmin =3D None, cmax =3D None,=20
zaxis_min =3D None, zaxis_max =3D None, split =3D 0) :
# plzcont (nverts, xyzverts, contours =3D 8, scale =3D "lin", clear =3D =
1,
# edges =3D 0, color =3D None, cmin =3D None, cmax =3D None, split =3D =
0
# zaxis_min =3D None, zaxis_max =3D None, )
# Plot filled z contours on the specified surface. NVERTS and
# XYZVERTS arrays specify the polygons for the surface being
# drawn. CONTOURS can be one of the following:
# N, an integer: Plot N contours (therefore, N+1 colored
# components of the surface)
# CVALS, a vector of floats: draw the contours at the
# specified levels.
# SCALE can be "lin", "log", or "normal" specifying the
# contour scale. (Only applicable if contours =3D N, of course).
# If CLEAR =3D 1, clear the display list first.
# If EDGES =3D 1, plot the edges.
# The algorithm is to apply slice2x repeatedly to the surface.
# If color =3D=3D None, then bytscl the palette into N + 1 colors
# and send each of the slices to pl3tree with the appropriate color.
# If color =3D=3D "bg", will plot only the edges.
# If CMIN is given, use it instead of the minimum z actually
# being plotted in the computation of contour levels. If CMAX is =
given,
# use it instead of the maximum z actually being plotted in the
# computation of contour levels. This is done so that a component
# of a larger graph will have the same colors at the same levels
# as every other component, rather than its levels being based
# on its own max and min, which may lie inside those of the
# rest of the graph.
# ZAXIS_MIN and ZAXIS_MAX represent axis limits on z as expressed
# by the user. If present, ZAXIS_MIN will inhibit plotting of all
# lesser z values, and ZAXIS_MAX will inhibit the plotting of all
# greater z values.
# =
------------------------------------------------------------------------
# 1. Get contour colors
if type (contours) =3D=3D IntType :
n =3D contours
if cmin !=3D None :
vcmin =3D cmin
minz =3D min (xyzverts [:, 2])
else :
vcmin =3D min (xyzverts [:, 2])
minz =3D vcmin
if cmax !=3D None :
vcmax =3D cmax
maxz =3D max (xyzverts [:, 2])
else :
vcmax =3D max (xyzverts [:, 2])
maxz =3D vcmax
if scale =3D=3D "lin" :
vc =3D vcmin + arange (1, n + 1, typecode =3D Float) * \
(vcmax - vcmin) / (n + 1)
elif scale =3D=3D "log" :
vc =3D vcmin + exp (arange (1, n + 1, typecode =3D Float) * =
\
log (vcmax - vcmin) / (n + 1))
elif scale =3D=3D "normal" :
zlin =3D xyzverts [:, 2]
lzlin =3D len (zlin)
zbar =3D add.reduce (zlin) / lzlin
zs =3D sqrt ( (add.reduce (zlin ** 2) - lzlin * zbar ** 2) /
(lzlin - 1))
z1 =3D zbar - 2. * zs
z2 =3D zbar + 2. * zs
diff =3D (z2 - z1) / (n - 1)
vc =3D z1 + arange (n) * diff
else :
raise _ContourError, "Incomprehensible scale parameter."
elif type (contours) =3D=3D ArrayType and contours.typecode () =
=3D=3D Float :
n =3D len (contours)
vc =3D sort (contours)
else :
raise _ContourError, "Incorrect contour specification."
if split =3D=3D 0 :
colors =3D (arange (n + 1, typecode =3D Float) * (199. / =
n)).astype ('b')
else :
colors =3D (arange (n + 1, typecode =3D Float) * (99. / =
n)).astype ('b')
# 2. Loop through slice2x calls
nv =3D array (nverts, copy =3D 1)
xyzv =3D array (xyzverts, copy =3D 1)
if clear =3D=3D 1 :
clear3 ( ) # Clear any previous plot or we're in trouble
# find imin--contours below this number need not be computed,
# and imax--contours at this level and above need not be computed.
imin =3D imax =3D 0
for i in range (n) :
if vc [i] <=3D minz :
imin =3D i + 1
if vc [i] >=3D maxz :
imax =3D i
break
if i =3D=3D n - 1 :
imax =3D n
# now make sure that the minimum and maximum contour levels =
computed
# are not outside the axis limits.
if zaxis_min !=3D None and zaxis_min > vc [imin] :
for i in range (imin, imax) :
if i + 1 < imax and zaxis_min > vc [i + 1] :
imin =3D i + 1
else :
break
vc [imin] =3D zaxis_min
if zaxis_max !=3D None and zaxis_max < vc [imax - 1] :
for i in range (imax - imin) :
if imax - 2 >=3D imin and zaxis_max < vc [imax - 2] :
imax =3D imax - 1
else :
break
vc [imax - 1] =3D zaxis_max
for i in range (imin, imax) :
[nv, xyzv, d1, nvb, xyzvb, d2] =3D \
slice2x (array ( [0., 0., 1., vc [i]], Float) , nv, xyzv, =
None)
if i =3D=3D imin and zaxis_min !=3D None and zaxis_min =3D=3D vc =
[i]:
# Don't send the "back" surface if it's below zaxis_min.
continue
else:
if color =3D=3D None :
pl3tree (nvb, xyzvb, (ones (len (nvb)) * colors =
[i]).astype ('b'),
split =3D 0, edges =3D edges)
else :
# N. B. Force edges to be on, otherwise the graph is =
empty.
pl3tree (nvb, xyzvb, "bg", split =3D 0, edges =3D 1)
if zaxis_max =3D=3D None or vc [imax - 1] < zaxis_max:
# send "front" surface if it's not beyond zaxis_max
if color =3D=3D None :
pl3tree (nv, xyzv, (ones (len (nv)) * colors [i]).astype =
('b'),
split =3D 0, edges =3D edges)
else :
pl3tree (nv, xyzv, "bg", split =3D 0, edges =3D 1)
def pl4cont (nverts, xyzverts, values, contours =3D 8, scale =3D "lin", =
clear =3D 1,
edges =3D 0, color =3D None, cmin =3D None, cmax =3D None,
caxis_min =3D None, caxis_max =3D None, split =3D 0) :
# pl4cont (nverts, xyzverts, values, contours =3D 8, scale =3D "lin", =
clear =3D 1,
# edges =3D 0, color =3D None, cmin =3D None, cmax =3D None,
# caxis_min =3D None, caxis_max =3D None, split =3D 0)
# Plot filled z contours on the specified surface. VALUES is
# a node-centered array the same length as SUM (NVERTS) whose
# contours will be drawn. NVERTS and
# XYZVERTS arrays specify the polygons for the surface being
# drawn. CONTOURS can be one of the following:
# N, an integer: Plot N contours (therefore, N+1 colored
# components of the surface)
# CVALS, a vector of floats: draw the contours at the
# specified levels.
# SCALE can be "lin", "log", or "normal" specifying the
# contour scale. (Only applicable if contours =3D N, of course).
# If CLEAR =3D=3D 1, clear the display list first.
# If EDGES =3D=3D 1, plot the edges.
# The algorithm is to apply slice2x repeatedly to the surface.
# If color =3D=3D None, then bytscl the palette into N + 1 colors
# and send each of the slices to pl3tree with the appropriate color.
# If color =3D=3D "bg", will plot only the edges.
# If CMIN is given, use it instead of the minimum c actually
# being plotted in the computation of contour levels. If CMAX is =
given,
# use it instead of the maximum c actually being plotted in the
# computation of contour levels. This is done so that a component
# of a larger graph will have the same colors at the same levels
# as every other component, rather than its levels being based
# on its own max and min, which may lie inside those of the
# rest of the graph.
# CAXIS_MIN and CAXIS_MAX represent axis limits on c as expressed
# by the user. If present, CAXIS_MIN will inhibit plotting of all
# lesser c values, and CAXIS_MAX will inhibit the plotting of all
# greater c values.
# =
------------------------------------------------------------------------
# 1. Get contour colors
if type (contours) =3D=3D IntType :
n =3D contours
if cmin !=3D None :
vcmin =3D cmin
minz =3D min (values)
else :
vcmin =3D min (values)
minz =3D vcmin
if cmax !=3D None :
vcmax =3D cmax
maxz =3D max (values)
else :
vcmax =3D max (values)
maxz =3D vcmax
if scale =3D=3D "lin" :
vc =3D vcmin + arange (1, n + 1, \
typecode =3D Float) * \
(vcmax - vcmin) / (n + 1)
elif scale =3D=3D "log" :
vc =3D vcmin + exp (arange (1, n + 1, \
typecode =3D Float) * \
log (vcmax - vcmin) / (n + 1))
elif scale =3D=3D "normal" :
zbar =3D add.reduce (values) / lzlin
zs =3D sqrt ( (add.reduce (values ** 2) - lzlin * zbar ** 2) =
/
(lzlin - 1))
z1 =3D zbar - 2. * zs
z2 =3D zbar + 2. * zs
diff =3D (z2 - z1) / (n - 1)
vc =3D z1 + arange (n) * diff
else :
raise _ContourError, "Incomprehensible scale parameter."
elif type (contours) =3D=3D ArrayType and contours.typecode () =
=3D=3D Float :
n =3D len (contours)
vc =3D sort (contours)
else :
raise _ContourError, "Incorrect contour specification."
if split =3D=3D 0 :
colors =3D (arange (n + 1, typecode =3D Float) * (199. / =
n)).astype ('b')
else :
colors =3D (arange (n + 1, typecode =3D Float) * (99. / =
n)).astype ('b')
# 2. Loop through slice2x calls
nv =3D array (nverts, copy =3D 1)
xyzv =3D array (xyzverts, copy =3D 1)
vals =3D array (values, copy =3D 1)
if clear =3D=3D 1 :
clear3 ( ) # Clear any previous plot or we're in trouble
# find imin--contours below this number need not be computed,
# and imax--contours at this level and above need not be computed.
imin =3D imax =3D 0
for i in range (n) :
if vc [i] <=3D minz :
imin =3D i + 1
if vc [i] >=3D maxz :
imax =3D i
break
if i =3D=3D n - 1 :
imax =3D n
# now make sure that the minimum and maximum contour levels =
computed
# are not outside the axis limits.
if caxis_min !=3D None and caxis_min > vc [imin] :
for i in range (imin, imax) :
if i + 1 < imax and caxis_min > vc [i + 1] :
imin =3D i + 1
else :
break
vc [imin] =3D caxis_min
if caxis_max !=3D None and caxis_max < vc [imax - 1] :
for i in range (imax - imin) :
if imax - 2 >=3D imin and caxis_max < vc [imax - 2] :
imax =3D imax - 1
else :
break
vc [imax - 1] =3D caxis_max
for i in range (n) :
if vc [i] <=3D minz :
continue
if vc [i] >=3D maxz :
break
[nv, xyzv, vals, nvb, xyzvb, d2] =3D \
slice2x (vc [i], nv, xyzv, vals)
if i =3D=3D imin and caxis_min !=3D None and caxis_min =3D=3D vc =
[i]:
# Don't send the "back" surface if it's below caxis_min.
continue
else:
if color =3D=3D None :
pl3tree (nvb, xyzvb, (ones (len (nvb)) * colors =
[i]).astype ('b'),
split =3D 0, edges =3D edges)
else :
# N. B. Force edges to be on, otherwise the graph is =
empty.
pl3tree (nvb, xyzvb, "bg", split =3D 0, edges =3D 1)
if caxis_max =3D=3D None or vc [imax - 1] < caxis_max:
# send "front" surface if it's not beyond caxis_max
if color =3D=3D None :
pl3tree (nv, xyzv, (ones (len (nv)) * colors [i]).astype =
('b'),
split =3D 0, edges =3D edges)
else :
pl3tree (nv, xyzv, "bg", split =3D 0, edges =3D 1)
def slice2x (plane, nverts, xyzverts, values =3D None) :
# slice2x (plane, nverts, xyzverts, values =3D None)
# Slice a polygon list, retaining only those polygons or
# parts of polygons on the positive side of PLANE, that is,
# the side where xyz(+)*PLANE(+:1:3)-PLANE(4) > 0.0.
# The NVERTS, VALUES, and XYZVERTS arrays serve as both
# input and output, and have the meanings of the return
# values from the slice3 function.
# Actually, since Python can't treat an argument as an output,
# this routine will return a sextuple of values (None for
# missing args).=20
# Note (ZCM 2/24/97) Reomved _slice2x as a global and added
# it as a final argument to slice2.
retval =3D slice2 (plane, nverts, xyzverts, values, 1)
retval =3D retval + [None] * (6 - len (retval))
return retval
_Pl3surfError =3D "Pl3surfError"
def pl3surf(nverts, xyzverts =3D None, values =3D None, cmin =3D None, =
cmax =3D None,
lim =3D None, edges =3D 0) :
# pl3surf (nverts, xyzverts)
# or pl3surf (nverts, xyzverts, values)
# Perform simple 3D rendering of an object created by slice3
# (possibly followed by slice2). NVERTS and XYZVERTS are polygon
# lists as returned by slice3, so XYZVERTS is sum(NVERTS)-by-3,
# where NVERTS is a list of the number of vertices in each polygon.
# If present, the VALUES should have the same length as NVERTS;
# they are used to color the polygon. If VALUES is not specified,
# the 3D lighting calculation set up using the light3 function
# will be carried out. Keywords cmin=3D and cmax=3D as for plf, pli,
# or plfp are also accepted. (If you do not supply VALUES, you
# probably want to use the ambient=3D keyword to light3 instead of
# cmin=3D here, but cmax=3D may still be useful.)
_draw3 =3D get_draw3_ ( )
if type (nverts) =3D=3D ListType :
list =3D nverts
nverts =3D list [0]
xyzverts =3D array (list [1], copy =3D 1)
values =3D list [2]
cmin =3D list [3]
cmax =3D list [4]
edges =3D list [6]
## Scale xyzverts to avoid loss of accuracy
minx =3D min (xyzverts [:, 0])
maxx =3D max (xyzverts [:, 0])
miny =3D min (xyzverts [:, 1])
maxy =3D max (xyzverts [:, 1])
minz =3D min (xyzverts [:, 2])
maxz =3D max (xyzverts [:, 2])
xyzverts [:, 0] =3D (xyzverts [:, 0] - minx) / (maxx - minx)
xyzverts [:, 1] =3D (xyzverts [:, 1] - miny) / (maxy - miny)
xyzverts [:, 2] =3D (xyzverts [:, 2] - minz) / (maxz - minz)
xyztmp =3D get3_xy (xyzverts, 1)
x =3D xyztmp [:, 0]
y =3D xyztmp [:, 1]
z =3D xyztmp [:, 2]
if values =3D=3D None :
# xyzverts [:, 0] =3D x
# xyzverts [:, 1] =3D y
# xyzverts [:, 2] =3D z
values =3D get3_light (xyztmp, nverts)
[list, vlist] =3D sort3d (z, nverts)
nverts =3D take (nverts, list)
values =3D take (values, list)
x =3D take (x, vlist)
y =3D take (y, vlist)
_square =3D get_square_ ( )
[_xfactor, _yfactor] =3D get_factors_ ()
xmax =3D max (x)
xmin =3D min (x)
ymax =3D max (y)
ymin =3D min (y)
xdif =3D xmax - xmin
ydif =3D ymax - ymin
if _xfactor !=3D 1. :
xmax =3D xmax + (_xfactor - 1.) * xdif /2.
xmin =3D xmin - (_xfactor - 1.) * xdif /2.
if _yfactor !=3D 1. :
ymax =3D ymax + (_yfactor - 1.) * ydif /2.
ymin =3D ymin - (_yfactor - 1.) * ydif /2.
if _square :
xdif =3D xmax - xmin
ydif =3D ymax - ymin
if xdif > ydif :
dif =3D (xdif - ydif) / 2.
ymin =3D ymin - dif
ymax =3D ymax + dif
elif ydif > xdif :
dif =3D (ydif - xdif) / 2.
xmin =3D xmin - dif
xmax =3D xmax + dif
plfp (values, y, x, nverts, cmin =3D cmin, cmax =3D cmax, legend =
=3D "",
edges =3D edges)
return [xmin, xmax, ymin, ymax]
nverts =3D array (nverts, Int)
xyzverts =3D array (xyzverts, Float )
if shape (xyzverts) [0] !=3D sum (nverts) or sum (less (nverts, 3)) =
or \
nverts.typecode () !=3D Int :
raise _Pl3surfError, "illegal or inconsistent polygon list"
if values !=3D None and len (values) !=3D len (nverts) :
raise _Pl3surfError, "illegal or inconsistent polygon color =
values"
if values !=3D None :
values =3D array (values, Float )
clear3 ( )
set3_object ( pl3surf, [nverts, xyzverts, values, cmin, cmax, lim, =
edges])
if (_draw3) :
# Plot the current list if _draw3 has been set.
call_idler ( )
if lim :
tmp =3D get3_xy (xyzverts, 1)
return max ( max (abs (tmp [:,0:2])))
else :
return None
# =
------------------------------------------------------------------------
_Pl3treeError =3D "Pl3treeError"
def pl3tree (nverts, xyzverts =3D None, values =3D None, plane =3D None,
cmin =3D None, cmax =3D None, split =3D 1, edges =3D 0) :
# pl3tree, nverts, xyzverts
# or pl3tree, nverts, xyzverts, values, plane
# Add the polygon list specified by NVERTS (number of vertices in
# each polygon) and XYZVERTS (3-by-sum(NVERTS) vertex coordinates)
# to the currently displayed b-tree. If VALUES is specified, it
# must have the same dimension as NVERTS, and represents the color
# of each polygon. (ZCM 7/18/97) Or, if VALUES =3D=3D "bg" =
("background")
# Then each polygon will be filled with the background color,
# giving just a wire frame. If VALUES is not specified, the polygons
# are assumed to form an isosurface which will be shaded by the
# current 3D lighting model; the isosurfaces are at the leaves of
# the b-tree, sliced by all of the planes. If PLANE is specified,
# the XYZVERTS must all lie in that plane, and that plane becomes
# a new slicing plane in the b-tree.
# Each leaf of the b-tree consists of a set of sliced isosurfaces.
# A node of the b-tree consists of some polygons in one of the
# planes, a b-tree or leaf entirely on one side of that plane, and
# a b-tree or leaf on the other side. The first plane you add
# becomes the root node, slicing any existing leaf in half. When
# you add an isosurface, it propagates down the tree, getting
# sliced at each node, until its pieces reach the existing leaves,
# to which they are added. When you add a plane, it also propagates
# down the tree, getting sliced at each node, until its pieces
# reach the leaves, which it slices, becoming the nodes closest to
# the leaves.
# This structure is relatively easy to plot, since from any
# viewpoint, a node can always be plotted in the order from one
# side, then the plane, then the other side.
# This routine assumes a "split palette"; the colors for the
# VALUES will be scaled to fit from color 0 to color 99, while
# the colors from the shading calculation will be scaled to fit
# from color 100 to color 199. (If VALUES is specified as a char
# array, however, it will be used without scaling.)
# You may specifiy a cmin=3D or cmax=3D keyword to affect the
# scaling; cmin is ignored if VALUES is not specified (use the
# ambient=3D keyword from light3 for that case).
# (ZCM 4/23/97) Add the split keyword. This will determine
# whether or not to split the palette (half to the isosurfaces
# for shading and the other half to plane sections for contouring).
# (ZCM 7/17/97) Add a calculation of the maximum and minimum
# of everything that is put into the tree. This cures distortion
# caused by loss of accuracy in orientation calculations.
# What is now put on the display list is pl3tree and [tree, minmax];
# both components are passed to _pl3tree to normalize results.
# avoid overhead of local variables for _pl3tree and _pl3leaf
# -- I don't know if this is such a big deal
_draw3 =3D get_draw3_ ()
if type (nverts) =3D=3D ListType :
_nverts =3D []
for i in range (len (nverts)) :
_nverts.append (nverts [i])
return _pl3tree (_nverts [0], nverts [1])
# We need copies of everything, or else arrays get clobbered.
nverts =3D array (nverts, Int)
xyzverts =3D array (xyzverts, Float )
if values =3D=3D "background" :
values =3D "bg"
elif values !=3D None and values !=3D "bg" :
values =3D array (values, values.typecode ())
if plane !=3D None :
plane =3D plane.astype (Float)
if shape (xyzverts) [0] !=3D sum (nverts) or sum (less (nverts, 3)) > =
0 or \
type (nverts [0]) !=3D IntType :
print "Dim1 of xyzverts ", shape (xyzverts) [0], " sum (nverts) =
",\
sum (nverts), " sum (less (nverts, 3)) ", sum (less (nverts, =
3)), \
" type (nverts [0]) ", `type (nverts [0])`
raise _Pl3treeError, "illegal or inconsistent polygon list."
if type (values) =3D=3D ArrayType and len (values) !=3D len (nverts) =
and \
len (values) !=3D sum (nverts) :
raise _Pl3treeError, "illegal or inconsistent polygon color =
values"
if type (values) =3D=3D ArrayType and len (values) =3D=3D sum =
(nverts) :
# We have vertex-centered values, which for Gist must be
# averaged over each cell
list =3D zeros (sum (nverts), Int)
array_set (list, cumsum (nverts) [0:-1], ones (len (nverts), Int))
tpc =3D values.typecode ()
values =3D (histogram (cumsum (list), values) / nverts).astype =
(tpc)
if plane !=3D None :
if (len (shape (plane)) !=3D 1 or shape (plane) [0] !=3D 4) :
raise _Pl3treeError, "illegal plane format, try plane3 =
function"
# Note: a leaf is going to be a list of lists.
leaf =3D [ [nverts, xyzverts, values, cmin, cmax, split, edges]]
## max and min of current leaf
minmax =3D array ( [min (xyzverts [:, 0]), max (xyzverts [:, 0]),
min (xyzverts [:, 1]), max (xyzverts [:, 1]),
min (xyzverts [:, 2]), max (xyzverts [:, 2])])
# retrieve current b-tree (if any) from 3D display list
_draw3_list =3D get_draw3_list_ ()
_draw3_n =3D get_draw3_n_ ()
try :
tree =3D _draw3_list [_draw3_n:]
except :
tree =3D []
if tree =3D=3D [] or tree [0] !=3D pl3tree :
tree =3D [plane, [], leaf, []]
else :
oldminmax =3D tree [1] [1]
tree =3D tree [1] [0]
## Find new minmax for whole tree
minmax =3D array ( [min (minmax [0], oldminmax [0]),
max (minmax [1], oldminmax [1]),
min (minmax [2], oldminmax [2]),
max (minmax [3], oldminmax [3]),
min (minmax [4], oldminmax [4]),
max (minmax [5], oldminmax [5])])
_pl3tree_add (leaf, plane, tree)
set_multiple_components (1)
tmp =3D has_multiple_components ()
clear3 ()
set_multiple_components (tmp)
# plist (tree)
set3_object (pl3tree, [tree, minmax])
if (_draw3) :
## Plot the current list
call_idler ( )
palette_dict =3D {
"earth.gp" :
[array ([0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 3, 4, 5, 5, 6, 7, 8,
8, 9, 10, 11, 11, 12, 13, 14, 15, 15, 16, 17, 18, 18, 19,
20, 21, 22, 22, 23, 24, 25, 26, 26, 27, 28, 29, 30, 31, =
31,
32, 33, 34, 35, 36, 36, 37, 38, 39, 40, 41, 41, 42, 43, =
44,
45, 46, 47, 48, 48, 48, 49, 49, 50, 50, 51, 51, 52, 52, =
53,
53, 54, 54, 55, 55, 56, 56, 57, 57, 58, 58, 59, 59, 60, =
61,
61, 62, 62, 63, 63, 64, 64, 65, 65, 66, 67, 67, 68, 68, =
69,
69, 70, 71, 73, 76, 78, 81, 83, 86, 88, 91, 94, 96, 99, =
101,
104, 106, 109, 111, 114, 117, 119, 121, 122, 124, 126, =
128,
129, 131, 133, 135, 136, 138, 140, 141, 143, 145, 147, =
149,
150, 152, 154, 156, 157, 159, 161, 163, 165, 166, 168, =
170,
172, 174, 175, 177, 179, 181, 183, 183, 184, 184, 185, =
185,
186, 186, 187, 187, 187, 188, 188, 189, 189, 190, 190, =
190,
191, 191, 192, 192, 193, 195, 196, 197, 198, 199, 201, =
202,
203, 204, 205, 207, 208, 209, 210, 211, 213, 214, 215, =
216,
217, 219, 220, 221, 222, 223, 225, 226, 227, 228, 229, =
231,
232, 233, 234, 235, 237, 238, 239, 240, 241, 243, 244, =
245,
246, 247, 249, 250, 251, 252, 253, 255], 'b'),
array ( [0, 0, 0, 0, 0, 0, 0, 0, 3, 6, 8, 11, 13, 16, 18, 21, 23, =
26,
28, 31, 33, 36, 38, 41, 43, 45, 48, 50, 52, 55, 57, 59, =
61,
64, 66, 68, 70, 72, 74, 77, 79, 81, 83, 85, 87, 89, 91, =
93,
95, 97, 99, 100, 102, 104, 106, 108, 109, 111, 113, 115,
116, 118, 120, 121, 123, 125, 126, 128, 128, 129, 129, =
130,
131, 131, 132, 133, 133, 134, 134, 135, 136, 136, 137, =
138,
138, 139, 140, 140, 141, 141, 142, 143, 143, 144, 145, =
145,
146, 146, 147, 148, 148, 149, 150, 150, 151, 151, 152, =
153,
153, 154, 155, 155, 156, 156, 157, 158, 158, 159, 160, =
160,
161, 161, 162, 163, 163, 164, 165, 165, 166, 166, 167, =
168,
168, 168, 169, 169, 170, 170, 171, 171, 172, 172, 172, =
173,
173, 174, 174, 175, 175, 175, 176, 176, 177, 177, 178, =
178,
179, 179, 179, 180, 180, 181, 181, 182, 182, 183, 183, =
182,
181, 181, 180, 179, 178, 177, 176, 175, 174, 173, 172, =
171,
170, 169, 168, 167, 166, 165, 164, 163, 163, 164, 164, =
165,
165, 166, 167, 167, 168, 169, 170, 171, 172, 173, 174, =
175,
176, 177, 178, 179, 181, 182, 184, 185, 187, 188, 190, =
192,
194, 196, 198, 200, 202, 204, 206, 208, 211, 213, 215, =
218,
221, 223, 226, 229, 232, 235, 238, 241, 244, 248, 251, =
255],
'b'),
array ( [0, 46, 58, 69, 81, 92, 104, 116, 116, 116, 116, 116, 117,
117, 117, 117, 117, 118, 118, 118, 118, 118, 119, 119, =
119,
119, 119, 120, 120, 120, 120, 120, 121, 121, 121, 121, =
121,
122, 122, 122, 122, 122, 123, 123, 123, 123, 123, 124, =
124,
124, 124, 124, 125, 125, 125, 125, 125, 126, 126, 126, =
126,
126, 127, 127, 127, 127, 127, 128, 126, 125, 124, 123, =
122,
120, 119, 118, 117, 115, 114, 113, 111, 110, 109, 108, =
106,
105, 104, 102, 101, 100, 98, 97, 96, 94, 93, 92, 90, 89, =
88,
86, 85, 84, 82, 81, 80, 78, 77, 76, 74, 73, 71, 70, 71, =
72,
72, 73, 73, 74, 75, 75, 76, 76, 77, 77, 78, 79, 79, 80, =
80,
81, 82, 82, 82, 83, 83, 83, 84, 84, 84, 85, 85, 85, 86, =
86,
86, 87, 87, 87, 88, 88, 88, 89, 89, 89, 90, 90, 90, 91, =
91,
91, 92, 92, 92, 93, 93, 93, 94, 94, 94, 95, 95, 95, 96, =
96,
97, 97, 97, 98, 98, 98, 99, 99, 99, 100, 100, 100, 101, =
101,
104, 106, 108, 111, 113, 116, 118, 121, 123, 126, 129, =
131,
134, 137, 139, 142, 145, 148, 150, 153, 156, 159, 162, =
165,
168, 170, 173, 176, 179, 182, 185, 189, 192, 195, 198, =
201,
204, 207, 211, 214, 217, 220, 224, 227, 230, 234, 237, =
241,
244, 248, 251, 255] , 'b')],
"gray.gp" : [
array ( [
0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 16, =
17,
18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, =
33,
34, 35, 36, 37, 38, 39, 41, 42, 43, 44, 45, 46, 47, 48, =
49,
50, 51, 52, 53, 54, 55, 57, 58, 59, 60, 61, 62, 63, 64, =
65,
66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 79, 80, =
81,
82, 83, 84, 85, 86, 87, 89, 90, 91, 92, 93, 94, 95, 96, =
97,
98, 99, 100, 101, 102, 103, 105, 106, 107, 108, 109, 110,
111, 112, 113, 114, 115, 116, 117, 118, 119, 121, 122, =
123,
124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, =
135,
137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, =
148,
149, 150, 151, 153, 154, 155, 156, 157, 158, 159, 160, =
161,
162, 163, 164, 165, 166, 167, 169, 170, 171, 172, 173, =
174,
175, 176, 177, 178, 179, 180, 181, 182, 183, 185, 186, =
187,
188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, =
199,
201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, =
212,
213, 214, 215, 217, 218, 219, 220, 221, 222, 223, 224, =
225,
226, 227, 228, 229, 230, 231, 233, 234, 235, 236, 237, =
238,
239, 240, 241, 242, 243, 244, 245, 246, 247, 249, 250, =
251,
252, 253, 254, 255
] , 'b'),
array ( [
0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 16, =
17,
18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, =
33,
34, 35, 36, 37, 38, 39, 41, 42, 43, 44, 45, 46, 47, 48, =
49,
50, 51, 52, 53, 54, 55, 57, 58, 59, 60, 61, 62, 63, 64, =
65,
66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 79, 80, =
81,
82, 83, 84, 85, 86, 87, 89, 90, 91, 92, 93, 94, 95, 96, =
97,
98, 99, 100, 101, 102, 103, 105, 106, 107, 108, 109, 110,
111, 112, 113, 114, 115, 116, 117, 118, 119, 121, 122, =
123,
124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, =
135,
137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, =
148,
149, 150, 151, 153, 154, 155, 156, 157, 158, 159, 160, =
161,
162, 163, 164, 165, 166, 167, 169, 170, 171, 172, 173, =
174,
175, 176, 177, 178, 179, 180, 181, 182, 183, 185, 186, =
187,
188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, =
199,
201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, =
212,
213, 214, 215, 217, 218, 219, 220, 221, 222, 223, 224, =
225,
226, 227, 228, 229, 230, 231, 233, 234, 235, 236, 237, =
238,
239, 240, 241, 242, 243, 244, 245, 246, 247, 249, 250, =
251,
252, 253, 254, 255
] , 'b'),
array ( [
0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 16, =
17,
18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, =
33,
34, 35, 36, 37, 38, 39, 41, 42, 43, 44, 45, 46, 47, 48, =
49,
50, 51, 52, 53, 54, 55, 57, 58, 59, 60, 61, 62, 63, 64, =
65,
66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 79, 80, =
81,
82, 83, 84, 85, 86, 87, 89, 90, 91, 92, 93, 94, 95, 96, =
97,
98, 99, 100, 101, 102, 103, 105, 106, 107, 108, 109, 110,
111, 112, 113, 114, 115, 116, 117, 118, 119, 121, 122, =
123,
124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, =
135,
137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, =
148,
149, 150, 151, 153, 154, 155, 156, 157, 158, 159, 160, =
161,
162, 163, 164, 165, 166, 167, 169, 170, 171, 172, 173, =
174,
175, 176, 177, 178, 179, 180, 181, 182, 183, 185, 186, =
187,
188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, =
199,
201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, =
212,
213, 214, 215, 217, 218, 219, 220, 221, 222, 223, 224, =
225,
226, 227, 228, 229, 230, 231, 233, 234, 235, 236, 237, =
238,
239, 240, 241, 242, 243, 244, 245, 246, 247, 249, 250, =
251,
252, 253, 254, 255
] , 'b')
],
"heat.gp" : [
array ( [
0, 1, 2, 4, 5, 7, 8, 10, 11, 13, 15, 17, 18, 20, 21, 23, =
24,
26, 27, 28, 30, 31, 33, 34, 36, 37, 39, 40, 42, 43, 46, =
47,
49, 50, 52, 53, 55, 56, 57, 59, 60, 62, 63, 65, 66, 68, =
69,
70, 72, 73, 76, 78, 79, 81, 82, 84, 85, 86, 88, 89, 92, =
94,
95, 97, 98, 99, 101, 102, 104, 105, 108, 110, 111, 113, =
114,
115, 117, 118, 120, 121, 123, 124, 126, 127, 128, 130, =
131,
133, 134, 136, 139, 140, 141, 143, 144, 146, 147, 149, =
150,
152, 153, 155, 156, 157, 159, 160, 162, 163, 165, 166, =
169,
170, 172, 173, 175, 176, 178, 179, 181, 182, 185, 186, =
188,
189, 191, 192, 194, 195, 197, 198, 201, 202, 204, 205, =
207,
208, 210, 211, 212, 214, 215, 217, 218, 220, 221, 223, =
224,
226, 227, 228, 231, 233, 234, 236, 237, 239, 240, 241, =
243,
244, 246, 247, 249, 250, 252, 253, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255
] , 'b'),
array ( [
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 5, 7, 9, =
11,
15, 17, 18, 20, 22, 24, 26, 28, 30, 32, 35, 37, 39, 41, =
43,
45, 47, 49, 51, 52, 54, 56, 58, 60, 62, 64, 66, 68, 69, =
71,
75, 77, 79, 81, 83, 85, 86, 88, 90, 92, 94, 96, 98, 100,
102, 103, 105, 107, 109, 111, 115, 117, 119, 120, 122, =
124,
126, 128, 130, 132, 136, 137, 139, 141, 143, 145, 147, =
149,
151, 153, 156, 158, 160, 162, 164, 166, 168, 170, 171, =
173,
175, 177, 179, 181, 183, 185, 187, 188, 190, 192, 196, =
198,
200, 202, 204, 205, 207, 209, 211, 213, 215, 217, 219, =
221,
222, 224, 226, 228, 230, 232, 236, 238, 239, 241, 243, =
245,
247, 249, 251, 253
] , 'b'),
array ( [
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
7, 11, 15, 19, 23, 27, 31, 35, 39, 43, 51, 54, 58, 62, =
66,
70, 74, 78, 82, 86, 90, 94, 98, 102, 105, 109, 113, 117,
121, 125, 133, 137, 141, 145, 149, 153, 156, 160, 164, =
168,
172, 176, 180, 184, 188, 192, 196, 200, 204, 207, 215, =
219,
223, 227, 231, 235, 239, 243, 247, 251
] , 'b')
],
"rainbow.gp" : [
array ( [
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 245, 240, 235, 229, 224, 219, =
213,
208, 202, 197, 192, 186, 181, 175, 170, 159, 154, 149, =
143,
138, 132, 127, 122, 116, 111, 106, 100, 95, 89, 84, 73, =
68,
63, 57, 52, 46, 41, 36, 30, 25, 19, 14, 9, 3, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 2, 7, 18, 24, 29, 34, 40, 45, 50, 56, 61, =
67,
72, 77, 83, 88, 93, 104, 110, 115, 120, 126, 131, 136, =
142,
147, 153, 158, 163, 169, 174, 180, 190, 196, 201, 206, =
212,
217, 223, 228, 233, 239, 244, 249, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255
] , 'b'),
array ( [
0, 0, 0, 0, 0, 0, 0, 0, 5, 11, 16, 22, 27, 32, 38, 43, =
48,
54, 59, 65, 70, 75, 81, 91, 97, 102, 108, 113, 118, 124,
129, 135, 140, 145, 151, 156, 161, 167, 178, 183, 188, =
194,
199, 204, 210, 215, 221, 226, 231, 237, 242, 247, 253, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 250, 239, 234, 228, 223, 218, 212, =
207,
201, 196, 191, 185, 180, 174, 169, 164, 153, 148, 142, =
137,
131, 126, 121, 115, 110, 105, 99, 94, 88, 83, 78, 67, 62,
56, 51, 45, 40, 35, 29, 24, 18, 13, 8, 2, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0
] , 'b'),
array ( [
42, 36, 31, 26, 20, 15, 10, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, =
1,
12, 17, 23, 28, 33, 39, 44, 49, 55, 60, 66, 71, 76, 82, =
87,
98, 103, 109, 114, 119, 125, 130, 135, 141, 146, 152, =
157,
162, 168, 173, 184, 189, 195, 200, 205, 211, 216, 222, =
227,
232, 238, 243, 248, 254, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, =
255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 254, 249, =
243,
233, 227, 222, 217, 211, 206, 201
] , 'b')
],
"stern.gp" : [
array ( [
0, 18, 36, 54, 72, 90, 108, 127, 145, 163, 199, 217, 235,
254, 249, 244, 239, 234, 229, 223, 218, 213, 208, 203, =
197,
192, 187, 182, 177, 172, 161, 156, 151, 146, 140, 135, =
130,
125, 120, 115, 109, 104, 99, 94, 89, 83, 78, 73, 68, 63, =
52,
47, 42, 37, 32, 26, 21, 16, 11, 6, 64, 65, 66, 67, 68, =
69,
70, 71, 72, 73, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, =
85,
86, 87, 88, 89, 90, 91, 92, 93, 94, 96, 97, 98, 99, 100,
101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, =
112,
113, 114, 115, 117, 118, 119, 120, 121, 122, 123, 124, =
125,
126, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, =
139,
140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, =
151,
152, 153, 154, 155, 156, 157, 158, 160, 161, 162, 163, =
164,
165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, =
176,
177, 178, 179, 181, 182, 183, 184, 185, 186, 187, 188, =
189,
190, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, =
203,
204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, =
215,
216, 217, 218, 219, 220, 221, 222, 224, 225, 226, 227, =
228,
229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, =
240,
241, 242, 243, 245, 246, 247, 248, 249, 250, 251, 252, =
253,
254
] , 'b'),
array ( [
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 15, 16, 17, =
18,
19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 32, 33, =
34,
35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, =
49,
50, 51, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 64, 65, =
66,
67, 68, 69, 70, 71, 72, 73, 75, 76, 77, 78, 79, 80, 81, =
82,
83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 96, 97, =
98,
99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, =
110,
111, 112, 113, 114, 115, 117, 118, 119, 120, 121, 122, =
123,
124, 125, 126, 128, 129, 130, 131, 132, 133, 134, 135, =
136,
137, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, =
149,
150, 151, 152, 153, 154, 155, 156, 157, 158, 160, 161, =
162,
163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, =
174,
175, 176, 177, 178, 179, 181, 182, 183, 184, 185, 186, =
187,
188, 189, 190, 192, 193, 194, 195, 196, 197, 198, 199, =
200,
201, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, =
213,
214, 215, 216, 217, 218, 219, 220, 221, 222, 224, 225, =
226,
227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, =
238,
239, 240, 241, 242, 243, 245, 246, 247, 248, 249, 250, =
251,
252, 253, 254
] , 'b'),
array ( [
0, 1, 3, 5, 7, 9, 11, 13, 15, 17, 21, 23, 25, 27, 29, 31,
33, 35, 37, 39, 41, 43, 45, 47, 49, 51, 53, 55, 57, 59, =
63,
65, 67, 69, 71, 73, 75, 77, 79, 81, 83, 85, 87, 89, 91, =
93,
95, 97, 99, 101, 105, 107, 109, 111, 113, 115, 117, 119,
121, 123, 127, 129, 131, 133, 135, 137, 139, 141, 143, =
145,
149, 151, 153, 155, 157, 159, 161, 163, 165, 167, 169, =
171,
173, 175, 177, 179, 181, 183, 185, 187, 191, 193, 195, =
197,
199, 201, 203, 205, 207, 209, 211, 213, 215, 217, 219, =
221,
223, 225, 227, 229, 233, 235, 237, 239, 241, 243, 245, =
247,
249, 251, 255, 251, 247, 243, 238, 234, 230, 226, 221, =
217,
209, 204, 200, 196, 192, 187, 183, 179, 175, 170, 166, =
162,
158, 153, 149, 145, 141, 136, 132, 128, 119, 115, 111, =
107,
102, 98, 94, 90, 85, 81, 77, 73, 68, 64, 60, 56, 51, 47, =
43,
39, 30, 26, 22, 17, 13, 9, 5, 0, 3, 7, 15, 19, 22, 26, =
30,
34, 38, 41, 45, 49, 57, 60, 64, 68, 72, 76, 79, 83, 87, =
91,
95, 98, 102, 106, 110, 114, 117, 121, 125, 129, 137, 140,
144, 148, 152, 156, 159, 163, 167, 171, 175, 178, 182, =
186,
190, 194, 197, 201, 205, 209, 216, 220, 224, 228, 232, =
235,
239, 243, 247, 251
] , 'b')
],
"yarg.gp" : [
array ( [
255, 254, 253, 252, 251, 250, 249, 248, 246, 245, 244, =
243,
242, 241, 240, 239, 238, 237, 236, 235, 234, 233, 232, =
230,
229, 228, 227, 226, 225, 224, 223, 222, 221, 220, 219, =
218,
217, 216, 214, 213, 212, 211, 210, 209, 208, 207, 206, =
205,
204, 203, 202, 201, 200, 198, 197, 196, 195, 194, 193, =
192,
191, 190, 189, 188, 187, 186, 185, 184, 182, 181, 180, =
179,
178, 177, 176, 175, 174, 173, 172, 171, 170, 169, 168, =
166,
165, 164, 163, 162, 161, 160, 159, 158, 157, 156, 155, =
154,
153, 152, 150, 149, 148, 147, 146, 145, 144, 143, 142, =
141,
140, 139, 138, 137, 136, 134, 133, 132, 131, 130, 129, =
128,
127, 126, 125, 124, 123, 122, 121, 120, 118, 117, 116, =
115,
114, 113, 112, 111, 110, 109, 108, 107, 106, 105, 104, =
102,
101, 100, 99, 98, 97, 96, 95, 94, 93, 92, 91, 90, 89, 88,
86, 85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, =
72,
70, 69, 68, 67, 66, 65, 64, 63, 62, 61, 60, 59, 58, 57, =
56,
54, 53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, =
40,
38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, =
24,
22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, =
6,
5, 4, 3, 2, 1, 0
] , 'b'),
array ( [
255, 254, 253, 252, 251, 250, 249, 248, 246, 245, 244, =
243,=20
242, 241, 240, 239, 238, 237, 236, 235, 234, 233, 232, =
230,=20
229, 228, 227, 226, 225, 224, 223, 222, 221, 220, 219, =
218,=20
217, 216, 214, 213, 212, 211, 210, 209, 208, 207, 206, =
205,=20
204, 203, 202, 201, 200, 198, 197, 196, 195, 194, 193, =
192,=20
291, 190, 189, 188, 187, 186, 185, 184, 182, 181, 180, =
179,=20
278, 177, 176, 175, 174, 173, 172, 171, 170, 169, 168, =
166,=20
265, 164, 163, 162, 161, 160, 159, 158, 157, 156, 155, =
154,=20
253, 152, 150, 149, 148, 147, 146, 145, 144, 143, 142, =
141,=20
240, 139, 138, 137, 136, 134, 133, 132, 131, 130, 129, =
128,=20
127, 126, 125, 124, 123, 122, 121, 120, 118, 117, 116, =
115,
114, 113, 112, 111, 110, 109, 108, 107, 106, 105, 104, =
102,
101, 100, 99, 98, 97, 96, 95, 94, 93, 92, 91, 90, 89, 88, =
86,
85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, 72, =
70,
69, 68, 67, 66, 65, 64, 63, 62, 61, 60, 59, 58, 57, 56, =
54,
53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, 40, =
38,
37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, =
22,
21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 6, =
5,
4, 3, 2, 1, 0
] , 'b'),
array ( [
255, 254, 253, 252, 251, 250, 249, 248, 246, 245, 244, =
243,
242, 241, 240, 239, 238, 237, 236, 235, 234, 233, 232, =
230,
229, 228, 227, 226, 225, 224, 223, 222, 221, 220, 219, =
218,
217, 216, 214, 213, 212, 211, 210, 209, 208, 207, 206, =
205,
204, 203, 202, 201, 200, 198, 197, 196, 195, 194, 193, =
192,
191, 190, 189, 188, 187, 186, 185, 184, 182, 181, 180, =
179,
178, 177, 176, 175, 174, 173, 172, 171, 170, 169, 168, =
166,
165, 164, 163, 162, 161, 160, 159, 158, 157, 156, 155, =
154,
153, 152, 150, 149, 148, 147, 146, 145, 144, 143, 142, =
141,
140, 139, 138, 137, 136, 134, 133, 132, 131, 130, 129, =
128,
127, 126, 125, 124, 123, 122, 121, 120, 118, 117, 116, =
115,
114, 113, 112, 111, 110, 109, 108, 107, 106, 105, 104, =
102,
101, 100, 99, 98, 97, 96, 95, 94, 93, 92, 91, 90, 89, 88,
86, 85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, =
72,
70, 69, 68, 67, 66, 65, 64, 63, 62, 61, 60, 59, 58, 57, =
56,
54, 53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, =
40,
38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, =
24,
22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8,
6, 5, 4, 3, 2, 1, 0
] , 'b')
]
}
def split_palette ( * name) :
# split_palette
# or split_palette ("palette_name.gp")
# split the current palette or the specified palette into two
# parts; colors 0 to 99 will be a compressed version of the
# original, while colors 100 to 199 will be a gray scale.
if len (name) > 0 :
dum =3D palette (name [0])
del dum
r =3D zeros (240, 'b')
g =3D zeros (240, 'b')
b =3D zeros (240, 'b')
dum =3D palette (r, g, b, query =3D 1)
del dum
try : # r may be all zeros, in which case the following will fail:
n =3D max (nonzero (r)) + 1 # (Figure out how many entries there =
are)
except :
n =3D 0
if n < 100 :
dum =3D palette ("earth.gp")
dum =3D palette (r, g, b, query =3D 1)
del dum
n =3D max (max (nonzero (r)), max (nonzero (g)),
max (nonzero (b))) + 1
newr =3D zeros (200, 'b')
newg =3D zeros (200, 'b')
newb =3D zeros (200, 'b')
newr [0:100] =3D interp (r [0:n].astype (Float), arange (n, typecode =
=3D Float ),
arange (100, typecode =3D Float ) * n / 100).astype ('b')
newg [0:100] =3D interp (g [0:n].astype (Float), arange (n, typecode =
=3D Float ),
arange (100, typecode =3D Float ) * n / 100).astype ('b')
newb [0:100] =3D interp (b [0:n].astype (Float), arange (n, typecode =
=3D Float ),
arange (100, typecode =3D Float ) * n / 100).astype ('b')
newr [100:200] =3D (arange (100, typecode =3D Int) * 255 / 99).astype =
('b')
newg [100:200] =3D (arange (100, typecode =3D Int) * 255 / 99).astype =
('b')
newb [100:200] =3D (arange (100, typecode =3D Int) * 255 / 99).astype =
('b')
palette (newr, newg, newb)
def split_bytscl (x, upper, cmin =3D None, cmax =3D None) :
# split_bytscl(x, 0)
# or split_bytscl(x, 1)
# as bytscl function, but scale to the lower half of a split
# palette (0-99, normally the color scale) if the second parameter
# is zero or nil, or the upper half (100-199, normally the gray
# scale) if the second parameter is non-zero.
x =3D bytscl (x, cmin =3D cmin, cmax =3D cmax, top =3D =
99).astype('b')
if upper :
x =3D x + 100
return x
def _pl3tree (tree, minmax) :
# tree is a 4-element list like this:
# [plane, back_tree, inplane_leaf, front_tree]
# plane=3D tree [0] is None if this is just a leaf
# in which case, only inplane_leaf is not None
# back_tree=3D tree [1] is the part behind plane
# inplane_leaf=3D tree [2] is the part in the plane itself
# front_tree=3D tree [3] is the part in front of plane
if tree =3D=3D None or tree =3D=3D [] :
return None
if tree [0] =3D=3D None or tree [0] =3D=3D [] :
# only the leaf is non-nil (but not a plane)
return _pl3leaf ( tree [2], 1, minmax)
# apply the 3D coordinate transform to two points along the
# normal of the splitting plane to judge which is in front
xyz =3D get3_xy (array ( [ [0., 0., 0.],
[tree [0] [0], tree [0] [1], tree [0] [2]]], Float), =
1)
[x, y, z] =3D [xyz [:, 0], xyz [:, 1], xyz [:, 2]]
# plot the parts in order toward the current viewpoint
if z [1] >=3D z [0] :
q1 =3D _pl3tree (tree [1], minmax)
q2 =3D _pl3leaf (tree [2], 0, minmax)
q3 =3D _pl3tree (tree [3], minmax)
else :
q1 =3D _pl3tree (tree [3], minmax)
q2 =3D _pl3leaf (tree [2], 0, minmax)
q3 =3D _pl3tree (tree [1], minmax)
if q1 !=3D None :
if q2 !=3D None and q3 =3D=3D None :
return [min (q2 [0], q1 [0]),
max (q2 [1], q1 [1]),
min (q2 [2], q1 [2]),
max (q2 [3], q1 [3])]
elif q2 =3D=3D None and q3 !=3D None :
return [min (q3 [0], q1 [0]),
max (q3 [1], q1 [1]),
min (q3 [2], q1 [2]),
max (q3 [3], q1 [3])]
elif q2 !=3D None and q3 !=3D None :
return [min (q3 [0], q2 [0], q1 [0]),
max (q3 [1], q2 [1], q1 [1]),
min (q3 [2], q2 [2], q1 [2]),
max (q3 [3], q2 [3], q1 [3])]
else :
return q1
elif q2 !=3D None :
if q3 =3D=3D None :
return q2
else :
return [min (q2 [0], q3 [0]),
max (q2 [1], q3 [1]),
min (q2 [2], q3 [2]),
max (q2 [3], q3 [3])]
elif q3 !=3D None :
return q3
else :
return None
## from lp import *
def _pl3leaf (leaf, not_plane, minmax) :
=20
# count number of polys, number of vertices
_nverts =3D _xyzverts =3D 0
if type (leaf) =3D=3D ListType and type (leaf [0]) =3D=3D ListType :
for i in range (len (leaf)) :
[_nverts, _xyzverts] =3D _pl3tree_count ( leaf [i], _nverts, =
_xyzverts )
else :
[_nverts, _xyzverts] =3D _pl3tree_count ( leaf , _nverts, =
_xyzverts)
# accumulate polys and vertices into a single polygon list
# The type of array required for palettes is "Py_GpColor",
# which translates to "PyArray_UBYTE", which is selected
# with a second argument of 'b' to the zeros() function.
## _values =3D zeros (_nverts, 'b') # See below
old_nverts =3D _nverts
_nverts =3D zeros (_nverts, Int)
_x =3D zeros (_xyzverts, Float )
_y =3D zeros (_xyzverts, Float )
if (not_plane) :
# Not just straight assignment; make _z a separate copy
_z =3D zeros (_xyzverts, Float )
else :
_z =3D None
_list =3D 1
_vlist =3D 1
if type (leaf) =3D=3D ListType and type (leaf [0]) =3D=3D ListType :
if leaf [0] [2] !=3D "bg" :
_values =3D zeros (old_nverts, 'b')
else :
_values =3D "bg"
for i in range (len (leaf)) :
[_list, _vlist, _edges] =3D _pl3tree_accum ( leaf [i] , =
not_plane,
_x, _y, _z, _list, _vlist, _values, _nverts, minmax)
else :
if leaf [2] !=3D "bg" :
_values =3D zeros (old_nverts, 'b')
else :
_values =3D "bg"
[_list, _vlist, _edges] =3D _pl3tree_accum ( leaf , not_plane,
_x, _y, _z, _list, _vlist, _values, _nverts, minmax)
# sort the single polygon list
if not_plane :
[_list, _vlist] =3D sort3d (_z, _nverts)
_nverts =3D take (_nverts, _list)
if _values !=3D "bg" :
_values =3D take (_values, _list)
_x =3D take (_x, _vlist)
_y =3D take (_y, _vlist)
_square =3D get_square_ ( )
[_xfactor, _yfactor] =3D get_factors_ ()
xmax =3D max (_x)
xmin =3D min (_x)
ymax =3D max (_y)
ymin =3D min (_y)
xdif =3D xmax - xmin
ydif =3D ymax - ymin
if _xfactor !=3D 1. :
xmax =3D xmax + (_xfactor - 1.) * xdif /2.
xmin =3D xmin - (_xfactor - 1.) * xdif /2.
if _yfactor !=3D 1. :
ymax =3D ymax + (_yfactor - 1.) * ydif /2.
ymin =3D ymin - (_yfactor - 1.) * ydif /2.
if _square :
xdif =3D xmax - xmin
ydif =3D ymax - ymin
if xdif > ydif :
dif =3D (xdif - ydif) / 2.
ymin =3D ymin - dif
ymax =3D ymax + dif
elif ydif > xdif :
dif =3D (ydif - xdif) / 2.
xmin =3D xmin - dif
xmax =3D xmax + dif
if _values =3D=3D "bg" :
_values =3D None
plfp (_values, _y, _x, _nverts, legend =3D "", edges =3D _edges)
return [xmin, xmax, ymin, ymax]
def _pl3tree_count (item, _nverts, _xyzverts) :
return [_nverts + len (item [0]), _xyzverts + len (ravel (item [1])) =
/ 3]
def _pl3tree_accum (item, not_plane, _x, _y, _z, _list, _vlist, _values,
_nverts, minmax) :
# (ZCM 2/24/97) Add _x, _y, _z , _list, _vlist, _nverts as parameters =
to
# avoid use of globals. return the new values of _list, _vlist.
# (ZCM 7/16/97) Return item [6] (whether to show edges)
# (ZCM 7/17/97) Add parameter minmax to normalize values
# N. B.:
# item [0] is nverts
# item [1] is xyzverts
# item [2] is values (if present)
# item [3] is cmin (if present)
# item [4] is cmax (if present)
# item [5] is split (1 =3D split the palette, 0 =3D do not split)
# item [6] is edges (1 =3D show edges, 0 =3D do not show edges)
# I have cleaned up what I think is extremely obscure Yorick code
# apparently designed to avoid some overhead.
# N. B. avoid splitting the palette if split is 0. (ZCM 4/23/97)
_xyzverts =3D array (item [1], copy =3D 1) # protect copy in tree
# Normalize copy (I'm only going to do this if it's an
# isosurface or is not a plane or has multiple components.=20
# There is a real problem here.
# You get bad distortion without doing this if one coordinate
# is many orders of magnitude larger than the others but the
# others have significant figues. You also get bad distortion
# by doing this in the case of a single plane section
# when one coordinate is insignificant with
# respect to the others and doesn't have significant digits.
# It is awfully hard to come up with a numerical criterion for this.)
if item [2] =3D=3D None or not_plane or has_multiple_components ():
minx =3D minmax [0]
maxx =3D minmax [1]
miny =3D minmax [2]
maxy =3D minmax [3]
minz =3D minmax [4]
maxz =3D minmax [5]
_xyzverts [:, 0] =3D (_xyzverts [:, 0] - minx) / (maxx - minx)
_xyzverts [:, 1] =3D (_xyzverts [:, 1] - miny) / (maxy - miny)
_xyzverts [:, 2] =3D (_xyzverts [:, 2] - minz) / (maxz - minz)
if item [2] =3D=3D None :
# this is an isosurface to be shaded (no values specified)
_xyzverts =3D get3_xy (_xyzverts, 1)
# accumulate nverts and values
incr =3D len (item [0])
_nverts [ _list - 1: _list - 1 + incr] =3D item [0]
if item [5] !=3D 0 :
_values [ _list - 1: _list - 1 + incr] =3D split_bytscl (
get3_light (_xyzverts, item [0]), 1, cmin =3D 0.0,
cmax =3D item [4]).astype ('b')
else : # no split
_values [ _list - 1: _list - 1 + incr] =3D bytscl (
get3_light (_xyzverts, item [0]), cmin =3D 0.0,
cmax =3D item [4]).astype ('b')
_list =3D _list + incr
# accumulate x, y, and z
incr =3D shape (_xyzverts) [0]
_x [_vlist - 1:_vlist - 1 + incr] =3D _xyzverts [:, 0]
_y [_vlist - 1:_vlist - 1 + incr] =3D _xyzverts [:, 1]
if not_plane :
_z [_vlist - 1:_vlist - 1 + incr] =3D _xyzverts [:, 2]
_vlist =3D _vlist + incr
else :
# this is to be pseudo-colored since values are given
if (not_plane) :
__xyz =3D get3_xy (_xyzverts, 1)
else :
__xyz =3D get3_xy (_xyzverts, 0)
# accumulate nverts and values
incr =3D len (item [0])
_nverts [ _list - 1: _list - 1 + incr] =3D item [0]
if item [2] !=3D "bg" :
if (item [2]).typecode () !=3D 'b' :
if item [5] !=3D 0 :
_values [ _list - 1: _list - 1 + incr] =3D split_bytscl (
item [2], 0, cmin =3D item [3], cmax =3D item =
[4]).astype ('b')
else :
_values [ _list - 1: _list - 1 + incr] =3D bytscl (
item [2], cmin =3D item [3], cmax =3D item [4]).astype =
('b')
else :
_values [ _list - 1: _list - 1 + incr] =3D item [2]
_list =3D _list + incr
# accumulate x, y, and z
incr =3D shape (__xyz) [0]
_x [_vlist - 1:_vlist - 1 + incr] =3D __xyz [:, 0]
_y [_vlist - 1:_vlist - 1 + incr] =3D __xyz [:, 1]
if not_plane :
_z [_vlist - 1:_vlist - 1 + incr] =3D __xyz [:, 2]
_vlist =3D _vlist + incr
return [_list, _vlist, item [6]]
def _pl3tree_add (leaf, plane, tree) :
if tree !=3D None and tree !=3D [] and \
not is_scalar (tree) and tree [0] !=3D None :
# tree has slicing plane, slice new leaf or plane and descend
[back, leaf] =3D _pl3tree_slice (tree [0], leaf)
if back :
if len (tree) >=3D 2 and tree [1] !=3D None and tree [1] !=3D =
[] :
_pl3tree_add (back, plane, tree [1])
else :
tree [1] =3D [None, [], back, []]
if (leaf) :
if len (tree) >=3D 4 and tree [3] !=3D None and tree [3] !=3D =
[] :
_pl3tree_add (leaf, plane, tree [3])
else :
tree [3] =3D [None, [], leaf, []]
elif plane !=3D None :
# tree is just a leaf, but this leaf has slicing plane
tree [0] =3D plane
tmp =3D tree [2]
tree [2] =3D leaf
leaf =3D tmp # swap new leaf with original leaf
[back, leaf] =3D _pl3tree_slice (plane, leaf)
if (back) :
tree [1] =3D [None, [], back, []]
if (leaf) :
tree [3] =3D [None, [], leaf, []]
else :
# tree is just a leaf and this leaf has no slicing plane
tree [2] =3D leaf + tree [2]
def _pl3tree_slice (plane, leaf) :
back =3D frnt =3D None
for ll in leaf :
# each item in the leaf list is itself a list
nvf =3D ll [0]
if nvf !=3D None :
nvb =3D array (nvf, copy =3D 1)
else :
nvb =3D None
xyzf =3D ll [1]
if xyzf !=3D None :
xyzb =3D array (xyzf, copy =3D 1)
else :
xyzb =3D None
valf =3D ll [2]
if valf !=3D None :
valb =3D array (valf, copy =3D 1)
else :
valb =3D None
if len (ll) > 4 :
ll4 =3D ll [4]
else :
ll4 =3D None
if len (ll) > 5 :
ll5 =3D ll [5]
else :
ll5 =3D 1
if len (ll) > 6 :
ll6 =3D ll [6]
else :
ll6 =3D 0
[nvf, xyzf, valf, nvb, xyzb, valb] =3D \
slice2x (plane, nvf, xyzf, valf)
if nvf !=3D None :
if frnt !=3D None :
frnt =3D [ [nvf, xyzf, valf, ll [3], ll4, ll5, ll6]] + frnt
else :
frnt =3D [ [nvf, xyzf, valf, ll [3], ll4, ll5, ll6]]
if nvb !=3D None :
if back !=3D None :
back =3D [ [nvb, xyzb, valb, ll [3], ll4, ll5, ll6]] + back
else :
back =3D [ [nvb, xyzb, valb, ll [3], ll4, ll5, ll6]]
return [back, frnt]
_Pl3tree_prtError =3D "Pl3tree_prtError"
def pl3tree_prt () :
_draw3_list =3D get_draw3_list ()
_draw3_n =3D get_draw3_n_ ()
if len (_draw3_list) >=3D _draw3_n :
tree =3D _draw3_list [_draw3_n:]
if tree =3D=3D None or tree =3D=3D [] or tree [0] !=3D pl3tree :
raise _Pl3tree_prtError, "<current 3D display not a pl3tree>"
else :
tree =3D tree [1]
_pl3tree_prt (tree, 0)
def _pl3tree_prt (tree, depth) :
if tree =3D=3D None or tree =3D=3D [] :
return
indent =3D (" " * (1 + 2 * depth)) [0:-1]
print indent + "+DEPTH=3D " + `depth`
if len (tree) !=3D 4 :
print indent + "***error - not a tree"
print indent + "plane=3D " + `tree [0]`
back =3D tree [1]
list =3D tree [2]
frnt =3D tree [3]
if back =3D=3D None or back =3D=3D [] :
print indent + "back =3D []"
else :
_pl3tree_prt (back, depth + 1)
for leaf in list :
print indent + "leaf length=3D " + `len (leaf)`
print indent + "npolys=3D " + `len (leaf [0])` + \
", nverts=3D " + `sum (leaf [0])`
print indent + "nverts=3D " + `shape (leaf [1]) [0]` + \
", nvals=3D " + `len (leaf [2])`
if frnt =3D=3D None or frnt =3D=3D [] :
print indent + "frnt =3D []"
else :
_pl3tree_prt (frnt, depth + 1)
print indent + "+DEPTH=3D " + `depth`
# =
------------------------------------------------------------------------
def _isosurface_slicer (m3, chunk, iso_index, _value) :
# Have to remember here that getv3 can either return an array of=20
# values, or a 2-list consisting of values and the corresponding cell
# indices, the latter in the case of an irregular grid.
# Note: (ZCM 2/24/97) I have fixed slicers to return the vertex
# information and what used to be the global _xyz3, or None. Hence
# returning the tuple [tmp, None].
tmp =3D getv3 (iso_index, m3, chunk)
if type (tmp) =3D=3D ListType :
tmp[0] =3D tmp [0] - _value
else :
tmp =3D tmp - _value
return [tmp, None]
def _plane_slicer (m3, chunk, normal, projection) :
# (ZCM 2/24/97) In all cases, return x as the last element of
# the tuple. This eliminates the global _xyz3.
x =3D xyz3(m3,chunk)
irregular =3D type (chunk) =3D=3D ListType and len (chunk) =3D=3D 2 \
or type (chunk) =3D=3D ArrayType and len (shape (chunk)) =3D=3D 1 =
\
and type (chunk [0]) =3D=3D IntType
if irregular :
# Need to return a list, the first component being the x's,
# the second being the relative cell list, and the third an offset
verts =3D m3 [1] [0]
cell_offset =3D 0
if type (verts) =3D=3D ListType :
totals =3D m3 [1] [3]
if type (chunk) =3D=3D ListType :
fin =3D chunk [0] [1]
else :
fin =3D chunk [-1]
for j in range (len (verts)) :
if fin <=3D totals [j] :
break
if j > 0 :
cell_offset =3D totals [j - 1]
if type (chunk) =3D=3D ListType :
clist =3D arange (0, chunk [0] [1] - chunk [0] [0], typecode =
=3D Int)
else :
clist =3D chunk - cell_offset
# In the irregular case we know x is ncells by 3 by something
return [ [x [:,0] * normal [0] + x [:,1] * normal [1] + \
x [:,2] * normal [2] - projection, clist, cell_offset], x]
elif len (shape (x)) =3D=3D 5 : # It's ncells X 3 X 2 X 2 X 2
return [x [:,0] * normal [0] + x [:,1] * normal [1] + \
x [:,2] * normal [2] - projection, x]
else : # It's 3 X ni X nj X nk
return [x [0] * normal [0] + x [1] * normal [1] + x [2] * normal =
[2] -\
projection, x]
def xyz3 (m3, chunk) :
# xyz3 (m3, chunk)
# return vertex coordinates for CHUNK of 3D mesh M3. The CHUNK
# may be a list of cell indices, in which case xyz3 returns a
# (dimsof(CHUNK))x3x2x2x2 list of vertex coordinates. CHUNK may
# also be a mesh-specific data structure used in the slice3
# routine, in which case xyz3 may return a 3x(ni)x(nj)x(nk)
# array of vertex coordinates. For meshes which are logically
# rectangular or consist of several rectangular patches, this
# is up to 8 times less data, with a concomitant performance
# advantage. Use xyz3 when writing slicing functions or coloring
# functions for slice3.
xyz =3D m3 [0] [0] (m3, chunk)
return xyz
def xyz3_rect (m3, chunk) :
m3 =3D m3 [1]
if len (shape (chunk)) !=3D 1 :
c =3D chunk
# The difference here is that our arrays are 0-based, while
# yorick's are 1-based; and the last element in a range is not
# included in the result array.
return m3 [1] [:,c [0, 0] - 1:1 + c [1, 0], c [0, 1] - 1:1 + c [1, =
1] ,
c [0, 2] - 1:1 + c [1, 2]]
else :
# Need to create an array of m3 [1] values the same size and shape
# as what to_corners3 returns.
# To avoid exceedingly arcane calculations attempting to
# go backwards to a cell list, this branch returns the list
# [<function values>, chunk]
# ???????????? ^^^^^^^^^^^^
# Then it is trivial for slice3 to find a list of cell
# numbers in which fi has a negative value.
dims =3D m3 [0]
indices =3D to_corners3 (chunk, dims [1] + 1, dims [2] + 1)
no_cells =3D shape (indices) [0]
indices =3D ravel (indices)
retval =3D zeros ( (no_cells, 3, 2, 2, 2), Float )
m30 =3D ravel (m3 [1] [0, ...])
retval [:, 0, ...] =3D reshape (take (m30, indices), (no_cells, 2, =
2, 2))
m31 =3D ravel (m3 [1] [1, ...])
retval [:, 1, ...] =3D reshape (take (m31, indices), (no_cells, 2, =
2, 2))
m32 =3D ravel (m3 [1] [2, ...])
retval [:, 2, ...] =3D reshape (take (m32, indices), (no_cells, 2, =
2, 2))
return retval
_xyz3Error =3D "xyz3Error"
def xyz3_irreg (m3, chunk) :
xyz =3D m3 [1] [1]
if type (chunk) =3D=3D ListType and len (chunk) =3D=3D 2 :
no_cells =3D chunk [0] [1] - chunk [0] [0]
if type (m3 [1] [0]) =3D=3D ListType :
totals =3D m3 [1] [3]
start =3D chunk [0] [0]
fin =3D chunk [0] [1]
for i in range (len (totals)) :
if fin <=3D totals [i] :
break
verts =3D m3 [1] [0] [i]
if i > 0 :
start =3D start - totals [i - 1]
fin =3D fin - totals [i - 1]
ns =3D verts [start:fin]
shp =3D shape (verts)
else :
ns =3D m3 [1] [0] [chunk [0] [0]:chunk [0] [1]] # This is a =
kXnv array
shp =3D shape (m3 [1] [0])
elif type (chunk) =3D=3D ArrayType and len (shape (chunk)) =3D=3D 1 =
and \
type (chunk [0]) =3D=3D IntType :
# chunk is an absolute cell list
no_cells =3D len (chunk)
if type (m3 [1] [0]) =3D=3D ListType :
totals =3D m3 [1] [3]
fin =3D max (chunk)
for i in range (len (totals)) :
if fin <=3D totals [i] :
break
verts =3D m3 [1] [0] [i]
if i > 0 :
start =3D totals [i - 1]
else :
start =3D 0
verts =3D m3 [1] [0] [i]
ns =3D take (verts, chunk - start)
shp =3D shape (verts)
else :
ns =3D take (m3 [1] [0], chunk)
shp =3D shape (m3 [1] [0])
else :
raise _xyz3Error, "chunk parameter has the wrong type."
if shp [1] =3D=3D 8 : # hex
retval =3D zeros ( (no_cells, 3, 2, 2, 2), Float)
retval [:, 0] =3D \
reshape (take (xyz [0], ravel (ns)), (no_cells, 2, 2, 2))
retval [:, 1] =3D \
reshape (take (xyz [1], ravel (ns)), (no_cells, 2, 2, 2))
retval [:, 2] =3D \
reshape (take (xyz [2], ravel (ns)), (no_cells, 2, 2, 2))
elif shp [1] =3D=3D 6 : # prism
retval =3D zeros ( (no_cells, 3, 3, 2), Float)
retval [:, 0] =3D \
reshape (take (xyz [0], ravel (ns)), (no_cells, 3, 2))
retval [:, 1] =3D \
reshape (take (xyz [1], ravel (ns)), (no_cells, 3, 2))
retval [:, 2] =3D \
reshape (take (xyz [2], ravel (ns)), (no_cells, 3, 2))
elif shp [1] =3D=3D 5 : # pyramid
retval =3D zeros ( (no_cells, 3, 5), Float)
retval [:, 0] =3D \
reshape (take (xyz [0], ravel (ns)), (no_cells, 5))
retval [:, 1] =3D \
reshape (take (xyz [1], ravel (ns)), (no_cells, 5))
retval [:, 2] =3D \
reshape (take (xyz [2], ravel (ns)), (no_cells, 5))
elif shp [1] =3D=3D 4 : # tet
retval =3D zeros ( (no_cells, 3, 4), Float)
retval [:, 0] =3D \
reshape (take (xyz [0], ravel (ns)), (no_cells, 4))
retval [:, 1] =3D \
reshape (take (xyz [1], ravel (ns)), (no_cells, 4))
retval [:, 2] =3D \
reshape (take (xyz [2], ravel (ns)), (no_cells, 4))
else :
raise _xyz3Error, "Funny number of cell faces: " + `shp [1]`
return retval
def xyz3_unif (m3, chunk) :
m3 =3D m3 [1]
n =3D m3 [1]
if len (chunk.shape) !=3D 1 :
c =3D chunk
i =3D c [0] - 1
dn =3D c [1] + 1 - i
xyz =3D zeros ( (3, dn [0], dn [1], dn [2]), Float)
else :
dims =3D m3 [0]
nj =3D dims [1]
nk =3D dims [2]
njnk =3D nj * nk
zchunk =3D chunk % nk
ychunk =3D chunk / nk % nj
xchunk =3D chunk / njnk
xyz =3D zeros ( (len (chunk), 3, 2, 2, 2), Float )
ijk0 =3D array ( [zeros ( (2, 2), Int ), ones ( (2, 2), Int )])
ijk1 =3D array ( [ [0, 0], [1, 1]], Int )
ijk1 =3D array ( [ijk1, ijk1] , Int )
ijk2 =3D array ( [ [0, 1], [0, 1]], Int )
ijk2 =3D array ( [ijk2, ijk2] , Int )
if len (n) =3D=3D 2: # we have dxdydz and x0y0z0
dxdydz =3D n [0]
x0y0z0 =3D n [1]
if len (shape (chunk)) !=3D 1:
# Convert the increment and size into array coordinates
# -- consecutive values
xx =3D arange (dn [0], typecode =3D Float ) * dxdydz [0] / (dn =
[0] - 1)
yy =3D arange (dn [1], typecode =3D Float ) * dxdydz [1] / (dn =
[1] - 1)
zz =3D arange (dn [2], typecode =3D Float ) * dxdydz [2] / (dn =
[2] - 1)
xyz [0] =3D x0y0z0 [0] + i [0] * dxdydz [0] + multiply.outer (
multiply.outer ( xx, ones (dn [1], Float )),
ones (dn [2], Float ))
xyz [1] =3D x0y0z0 [1] + i [1] * dxdydz [0] + \
multiply.outer ( ones (dn [0], Float ), \
multiply.outer ( yy, ones (dn [2], Float )))
xyz [2] =3D x0y0z0 [2] + i [2] * dxdydz [0] + \
multiply.outer ( ones (dn [0], Float ), \
multiply.outer ( ones (dn [1], Float ), zz))
else :
# -- nonconsecutive values
xyz [:, 0] =3D add.outer ( xchunk, ijk0) * dxdydz [0] + x0y0z0 =
[0]
xyz [:, 1] =3D add.outer ( ychunk, ijk1) * dxdydz [1] + x0y0z0 =
[1]
xyz [:, 2] =3D add.outer ( zchunk, ijk2) * dxdydz [2] + x0y0z0 =
[2]
elif type (n) =3D=3D ListType and len (n) =3D=3D 3: # We have three =
one-dimensional arrays.
xx =3D n [0]
yy =3D n [1]
zz =3D n [2]
n0 =3D len (xx)
n1 =3D len (yy)
n2 =3D len (zz)
if len (shape (chunk)) !=3D 1:
# Convert the increment and size into array coordinates
# -- consecutive values
xyz [0] =3D multiply.outer (
multiply.outer ( xx [i [0]:i [0] + n0], ones (n1, Float )), =
\
ones (n2, Float ))
xyz [1] =3D multiply.outer ( ones (n0, Float ), \
multiply.outer ( yy [i [1]: i[1] + n1], ones (n2, Float )))
xyz [2] =3D multiply.outer ( ones (n0, Float ), \
multiply.outer ( ones (n1, Float ), zz [i [2]: i[2] + n2]))
else :
# -- nonconsecutive values
xyz [:, 0] =3D reshape (take (xx, ravel (add.outer (xchunk, =
ijk0))),
(len (chunk), 2, 2, 2))
xyz [:, 1] =3D reshape (take (yy, ravel (add.outer (ychunk, =
ijk1))),
(len (chunk), 2, 2, 2))
xyz [:, 2] =3D reshape (take (zz, ravel (add.outer (zchunk, =
ijk2))),
(len (chunk), 2, 2, 2))
return xyz
def to_corners3 (list, nj, nk) :
# to_corners3(list, nj, nk)
# convert an array of cell indices in an (ni-1)-by-(NJ-1)-by-(NK-1)
# logically rectangular grid of cells into the list of
# len(LIST)-by-2-by-2-by-2 cell corner indices in the
# corresponding ni-by-NJ-by-NK array of vertices.
# Note that this computation in Yorick gives an absolute offset
# for each cell quantity in the grid. In Yorick it is legal to
# index a multidimensional array with an absolute offset. In
# Python it is not. However, an array can be flattened if
# necessary.
# Other changes from Yorick were necessitated by row-major
# order and 0-origin indices, and of course the lack of
# Yorick array facilities.
njnk =3D nj * nk
kk =3D list / (nk - 1)
list =3D list + kk + nk * (kk / (nj - 1))
adder =3D array ( [ [ [0, 1], [nk, nk + 1]],
[ [njnk, njnk + 1], [nk + njnk, nk + njnk + 1]]])
res =3D zeros ( (len (list), 2, 2, 2), Int)
for i in range (len(list)):
res [i] =3D adder + list [i]
return res
------=_NextPart_000_0011_01BD93B1.A9ACA3C0
Content-Type: application/octet-stream;
name="gist.py"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="gist.py"
# Copyright (c) 1996, 1997, The Regents of the University of California.
# All rights reserved. See LEGAL.LLNL for full text and disclaimer.
from Numeric import *
import sys, os # To be sure expand_path has posixpath and we have =
sys.path
from gistC import *
from help import help
from shapetest import *
# Parameters used by pltitle and xytitles
pltitle_height=3D 18;
pltitle_font=3D "helvetica";
# Parameters used by plmk and plmk_default
_plmk_count =3D 0
_plmk_msize =3D _plmk_color =3D _plmk_width =3D None
# delete the current graphics window, or graphics window N (0-7).
def winkill(*N):
if N:
window(N[0], display=3D"", hcp=3D"")
else:
window(display=3D"", hcp=3D"")
# Plot TITLE centered above the coordinate system for any of the
# standard Gist styles. You will need to customize this for other
# plot styles.
def pltitle(title):
vp =3D viewport()
xmidpt =3D (vp[0] + vp[1])/2.0
plt( title, xmidpt, vp[3] + 0.02,
font=3Dpltitle_font, justify=3D"CB", height=3Dpltitle_height)
# Set limits on the y-axis
def ylimits(ymin=3D'u',ymax=3D'u'): limits('u','u',ymin,ymax)
# Return the 1-origin zone index for the point clicked in
# for the default mesh, or for the mesh (X,Y) (region array IREG).
# Number of args can be 0, 2, or 3.
def moush(*arg):
narg =3D len(arg)
if narg =3D=3D 3: # (y, x, ireg)
xy =3D mouse (-1, 0, "<Click mouse in mesh>")
if xy =3D=3D None: return None
return mesh_loc (xy[1], xy[0], arg[0], arg[1], arg[2]);
elif narg =3D=3D 2: # (y, x)
xy =3D mouse (-1, 0, "<Click mouse in mesh>")
if xy =3D=3D None: return None
return mesh_loc (xy[1], xy[0], arg[0], arg[1]);
elif narg =3D=3D 0: # ()
xy =3D mouse (-1, 0, "<Click mouse in mesh>")
if xy =3D=3D None: return None
return mesh_loc (xy[1], xy[0]);
else:
print "Moush takes 0, 2, or 3 args: ( [ y, x [ , ireg ] ] )"
return None
# Create an encapsulated PostScript file. Requires Ghostscript and its
# associated ps2epsi utility.
def eps(name):
import os
name =3D name + ".ps"
window (hcp =3D name, dump =3D 1, legends =3D 0)
hcp ()
window (hcp=3D"")
os.system ("ps2epsi " + name)
os.system ("rm " + name)
def xytitles(xtitle =3D "", ytitle =3D "", delta =3D (0.,0.)):
vp =3D viewport()
xmidpt =3D (vp[0] + vp[1])/2.0
ymidpt =3D (vp[2] + vp[3])/2.0
if len(xtitle) > 0:
plt(xtitle, xmidpt, vp[2] - 0.050 + delta[1],
font=3Dpltitle_font, justify=3D"CT", height=3Dpltitle_height)
if len(ytitle) > 0:
plt(ytitle, vp[0] - 0.050 + delta[0], ymidpt,
font=3Dpltitle_font, justify=3D"CB", height=3Dpltitle_height, =
orient=3D1)
def _spanz(lb,ub,n):
if n < 3: raise ValueError, '3rd arg must be at least 3'
c =3D 0.5*(ub - lb)/(n - 1.0)
b =3D lb + c
a =3D (ub - c - b)/(n - 2.0)
return map(lambda x,A=3Da,B=3Db: A*x + B, arange(n-1))
# predefined markers: square, +, delta, circle, diamond, x, grad
_seq =3D _spanz(-pi,pi,37)
_plmk_markers =3D (
array([[-1,1,1,-1],[-1,-1,1,1]])*.007,
array([[-4,-1,-1,1,1,4,4,1,1,-1,-1,-4],
[-1,-1,-4,-4,-1,-1,1,1,4,4,1,1]])*.007/sqrt(7),
array([[-sqrt(3),sqrt(3),0],[-1,-1,2]])*.007/sqrt(.75*sqrt(3)),
array([cos(_seq),sin(_seq)])*.007/(pi/4.),
array([[-1,0,1,0],[0,-1,0,1]])*.007*sqrt(2),
array([[-1,-2.5,-1.5,0,1.5,2.5,1,2.5,1.5,0,-1.5,-2.5],
=
[0,-1.5,-2.5,-1,-2.5,-1.5,0,1.5,2.5,1,2.5,1.5]])*.007*sqrt(2)/sqrt(7),
array([[0,sqrt(3),-sqrt(3)],[-2,1,1]])*.007/sqrt(.75*sqrt(3))
)
del(_seq)
def =
plmk(y,x=3DNone,marker=3DNone,width=3DNone,color=3DNone,msize=3DNone):
global _plmk_count
global _plmk_color, _plmk_width, _plmk_msize
color_dict =3D { 'bg':-1, 'fg':-2, 'black':-3, 'white':-4, 'red':-5,
'green':-6, 'blue':-7, 'cyan':-8, 'magenta':-9, 'yellow':-10 }
z =3D None
if marker =3D=3D None:
marker =3D _plmk_markers[(_plmk_count)%7]
_plmk_count =3D _plmk_count + 1
elif type(marker) =3D=3D type(0):
marker =3D _plmk_markers[marker-1];
xm =3D marker[0]
ym =3D marker[1]
if not msize: msize =3D _plmk_msize;
if msize:
xm =3D xm * msize;
ym =3D ym * msize;
if not color: color =3D _plmk_color;
ecolor =3D color;
if type(color) =3D=3D type(""):
color =3D color_dict[color];
=20
if not width: width =3D _plmk_width;
if width >=3D 10:
if not color:
color =3D ecolor =3D -2 # magic number for "fg"
z =3D ones(1+len(y)) * color
z =3D z.astype(UnsignedInt8) # convert array to type <unsigned char>
width =3D None
n =3D ones(1+len(y));
n[0] =3D len(ym);
if not x: x =3D 1 + arange(len(y));
plfp( z, concatenate((ym,y)), concatenate((xm,x)),n,
edges=3D1, ewidth=3Dwidth, ecolor=3Decolor)
# Set default color, msize, and width values for plmk. Use
# width=3D10 to get solid fills. With no parameters, plmk_default
# restores the initial default values.
def plmk_default(color=3DNone, msize=3DNone, width=3DNone):
global _plmk_color, _plmk_width, _plmk_msize
if color: _plmk_color =3D color
if width: _plmk_width =3D width
if msize: _plmk_msize =3D msize
if not (color or width or msize):
_plmk_msize =3D _plmk_color =3D _plmk_width =3D None
from arrayfns import *
from types import *
def spann (zmin, zmax, n =3D 8, fudge =3D 0, force =3D 0) :
# return no more than N equally spaced "nice" numbers between
# ZMIN and ZMAX.
# HAR! I added a "force" parameter to force exactly n values
# to be returned.
dz =3D (zmax - zmin) / max (float (n), 0.)
if dz =3D=3D 0. :
dz =3D abs (zmin)
if (dz !=3D 0.) :
power =3D floor (log10 (dz) + 0.00001)
base =3D dz / (10. ** power)
if base > 5.00001 :
base =3D 1.0
power =3D power + 1.0
elif base > 2.00001 :
base =3D 5.0
else :
base =3D 2.0
# Round dz up to nearest "nice" number
dz =3D base * 10.0 ** power
zmin =3D ceil (zmin / dz - fudge)
zmax =3D floor (zmax / dz + fudge)
if (force =3D=3D 0) :
nz =3D int (zmax - zmin + 1.0)
else :
nz =3D n
if nz > 1 :
levs =3D span (zmin * dz, zmax * dz, nz)
else :
if nz < 1 :
if base < 1.5 :
base =3D 5.0
power =3D power - 1
elif base < 2.5 :
base =3D 1.0
else :
base =3D 2.0
dz =3D base * 10.0 ** power
zmin =3D ceil (zmin / dz + 0.001)
levs =3D array ( [zmin * dz])
else :
levs =3D array ( [-1.0, 1.0])
return (levs)
_ContourError =3D "ContourError"
def plfc (z, y, x, ireg, contours =3D 20, colors =3D None, region =3D 0,
triangle =3D None, scale =3D "lin") :
#
# plfc (z, y, x, ireg, contours =3D 20, colors =3D None, region =3D 0,
# triangle =3D None, scale =3D "lin")
#
# 1. Get contour colors
vcmin =3D min (ravel (z))
vcmax =3D max (ravel (z))
if type (contours) =3D=3D IntType :
n =3D contours
vc =3D zeros (n + 2, Float)
vc [0] =3D vcmin
vc [n + 1] =3D vcmax
if scale =3D=3D "lin" or scale =3D=3D None :
# This stuff is in lieu of the spann stuff in Yorick.
# Note that my calculation and spann give different
# results, and both are different from how plc plots
# the contour lines. Thus if you want filled contours
# with lines superimposed on top, the lines will not
# match the contours.
# The only way to get around this is to send both
# plc and plfc the array of levels that you want.
vc [1:n + 1] =3D vcmin + arange (1, n + 1, typecode =3D =
Float) * \
(vcmax - vcmin) / (n + 1)
# vc [1:n + 1] =3D spann (vcmin, vcmax, n, force =3D 1)
elif scale =3D=3D "log" :
vc [1:n + 1] =3D vcmin + exp (arange (1, n + 1, typecode =3D =
Float) * \
log (vcmax - vcmin) / (n + 1))
elif scale =3D=3D "normal" :
zlin =3D ravel (z)
lzlin =3D len (zlin)
zbar =3D add.reduce (zlin) / lzlin
zs =3D sqrt ( (add.reduce (zlin ** 2) - lzlin * zbar ** 2) /
(lzlin - 1))
z1 =3D zbar - 2. * zs
z2 =3D zbar + 2. * zs
diff =3D (z2 - z1) / (n - 1)
vc [1:n + 1] =3D z1 + arange (n) * diff
else :
raise _ContourError, "Incomprehensible scale parameter."
elif type (contours) =3D=3D ArrayType and contours.typecode () =
=3D=3D Float :
n =3D len (contours)
vc =3D zeros (n + 2, Float)
vc [0] =3D vcmin
vc [n + 1] =3D vcmax
vc [1:n + 1] =3D sort (contours)
else :
raise _ContourError, "Incorrect contour specification."
if colors =3D=3D None :
colors =3D (arange (n + 1, typecode =3D Float) * (199. / =
n)).astype ('b')
else :
colors =3D array (colors)
if len (colors) !=3D n + 1 :
raise "PLFC_Error", \
"colors must specify one more color than contours."
if colors.typecode !=3D 'b' :
colors =3D bytscl (colors)
if triangle =3D=3D None :
triangle =3D zeros (z.shape, Int16)
# Set mesh first
plmesh (y, x, ireg, triangle =3D triangle)
for i in range (n + 1) :
[nc, yc, xc] =3D contour (array ( [vc [i], vc [i + 1]]), z)
if (is_scalar(nc) and nc =3D=3D 0 or nc =3D=3D None) :
continue
plfp ( (ones (len (nc)) * colors [i]).astype ('b'),
yc, xc, nc, edges =3D 0)
------=_NextPart_000_0011_01BD93B1.A9ACA3C0--