## Automatically adapted for scipy Oct 31, 2005 by # Copyright (c) 1996, 1997, The Regents of the University of California. # All rights reserved. See Legal.htm for full text and disclaimer. # The following is so I know about arrays: from scipy import * from numpy.core.umath import * from shapetest import * from surface import * from graftypes import * from plane import * from slice3 import * class Mesh3d ( Surface ) : """m = Mesh3d ( ) will create a mesh with coordinates and various other characteristics. The keywords are as follows: For all meshes: color_card = specifies which color card you wish to use, e. g., "rainbowhls" (the default), "random", etc. Although a characteristic of a Graph2d, it can be a Surface characteristic since 'link'ed surfaces can have different color cards. This is not implemented for Gist. Instead, Gist allows the option of a "split palette", which causes isosurfaces to be shaded and other surfaces to be colored by a subset of the current (or specified) palette. opt_3d = where is a string or a sequence of strings giving the 3d or 4d surface characteristics. A surface is colored by height in z if a 3d option is specified, and by the value of the function c if a 4d option is specified. With a wire grid option, the grid is colored; with a flat option, the quadrilaterals set off by grid points are colored; with a smooth option, the surface itself is colored by height; and with an iso option, the contour lines are colored. flat and iso options may be used together in any combination. wire grid options are independent of the other options. Legal arguments for set_3d_options are: 'wm'--monochrome wire grid; 'w3' and 'w4'--3d and 4d coloring of wire grid. 'f3' and 'f4'--flat 3d and 4d coloring options. 'i3' and 'i4'--3d and 4d isoline (contour line) options. 's3' and 's4'--3d and 4d smooth coloring options. For Gist plots, currently only 'wm', 'f3', and 'f4' are currently available. If others are specified, an intelligent guess will be made. mesh_type = in one of the wire modes, tells what form the wire grid takes: "x"--x lines only; "y": y lines only; "xy": both a lines and y lines. Not implemented for Gist plots. mask = controls hidden line removal. Allowed values are "none" : transparent graph; "min": simple masking; "max" : better masking; "sort": slowest but most sophisticated. Only the "none" and "sort" are available for Gist. z_c_switch = 0 or 1 : set to 1 means switch z and c in the plot. z_contours_scale, c_contours_scale = "lin" or "log" z_contours_array, c_contours_array = actual array of numbers to use for contours number_of_z_contours, number_of_c_contours = specifies how many contours to use; they will be computed automatically based on the data. The last six options are not available in Gist. For nonstructured meshes only: x = , y = , z = three vectors of equal lengths giving the coordinates of the nodes of a nonstructured mesh. c = , a vector of the same size as x, y, and z giving a data value at each of the points specified. c could also be an array of such vectors, when isosurfaces of more than one function are to be plotted. avs = 0 or 1: if 1, the input data represents a nonstructured mesh in a sort of AVS format, which will be explained in more detail below. The data will be translated into the Narcisse format prior to being sent to Narcisse. cell_descr = if present, this keyword signifies a nonstructured mesh submitted in the Narcisse format, as explained in the Narcisse manual. For nonstructured mesh data, one of these two keywords must be present. In the case avs = 1, at least one of the following four keywords (hex, tet, prism, or pyr) must be present also: hex = [ ] The entries in the list (in order) must be: an integer number n_zones, which is the number of hex cells in the mesh; and a matrix nz whose dimensions are n_zones by 8; nz [i][0], nz [i][1], ... nz [i][7] give the indices of the 8 vertices of the ith zone in canonical order (one side in the outward normal direction, then the corresponding vertices of the opposite side in the inward normal direction). tet = [ ] The list is the same format as for hex data. The matrix nz will be n_zones by 4, and each row gives the indices of the apex and then the base in inward normal order. prism = [ ] The list is the same format as for hex data. The matrix nz will be n_zones by 6, and each row gives the indices of one triangular side in the outward normal direction, then the corresponding vertices of the opposite side in the inward normal direction. pyr = [ ] The list is the same format as for hex data. The matrix nz will be n_zones by 5, and each row gives the indices of the apex and then the base in inward normal order. In the case of avs = 0, the following keywords are all required: cell_descr = , a vector of integers giving the description of the cells of the mesh, as specified in the Narcisse manual. no_cells = , the number of cells in the mesh. (Here cells refers to the number of faces of cells). For structured meshes only: x = , y = , z = three vectors giving the coordinates of the nodes of a structured rectangular mesh. c = , a three-dimensional array whose [i, j, k] element gives the associated data value at (x [i], y [j], z [k]). Could also be one less in each direction, giving a cell-centered quantity. c could also be an array of such arrays, when isosurfaces of more than one function are to be plotted. """ _MeshInitError = "MeshInitError" def type (self) : return Mesh3dType def __init__ ( self , *kwds , ** keywords ) : if len (kwds) == 1 : keywords = kwds[0] self.cell_dict = {} if keywords.has_key ("hex") : self.cell_dict ["hex"] = keywords ["hex"] if keywords.has_key ("tet") : self.cell_dict ["tet"] = keywords ["tet"] if keywords.has_key ("prism") : self.cell_dict ["prism"] = keywords ["prism"] if keywords.has_key ("pyr") : self.cell_dict ["pyr"] = keywords ["pyr"] if (keywords.has_key ("avs") and keywords ["avs"] == 0) or \ (not keywords.has_key ("avs") and not keywords.has_key ("cell_descr")) : # structured case works basically just like a Surface Surface.__init__ ( self, keywords ) self.structured = 1 else : # nonstructured case, at least initialize the generics self.structured = 0 self.generic_init ( keywords ) if keywords.has_key ("avs") : self.avs = keywords ["avs"] else : self.avs = None if not keywords.has_key ("x") or \ not keywords.has_key ("y") or \ not keywords.has_key ("z") or \ not keywords.has_key ("c") : raise self._MeshInitError , \ "x, y, z, and c keywords required for nonstructured mesh." self.x = keywords ["x"] self.y = keywords ["y"] self.z = keywords ["z"] self.c = keywords ["c"] if keywords.has_key ("cell_descr" ) : # already there if not keywords.has_key ("no_cells") : raise self._MeshInitError , \ "Keyword no_cells is required if avs is not specified." self.cell_descr = keywords ["cell_descr"] self.number_of_cells = keywords ["no_cells"] elif not keywords.has_key ("avs") : raise self._MeshInitError , \ "keywords must include avs = 1 or else cell_descr = values." # Conversion tables for AVS to Narcisse formats _cell_types = ["tet", "pyr", "prism", "hex"] _no_of_vertices = { "hex" : 8, "tet" : 4, "pyr" : 5, "prism" : 6 } _no_of_rect_faces = { "hex" : 6, "tet" : 0, "pyr" : 1, "prism" : 3 } _no_of_tri_faces = { "hex" : 0, "tet" : 4, "pyr" : 4, "prism" : 2 } # The next few items give the standard numbering of the faces of # the various cells. _hex_triangular_faces = [ ] _tet_triangular_faces = [ [0, 1, 2], [0, 3, 1], [0, 2, 3], [1, 3, 2]] _pyr_triangular_faces = [ [0, 1, 2], [0, 4, 1], [0, 3, 4], [0, 2, 3]] _prism_triangular_faces = [ [0, 1, 2], [3, 5, 4]] _hex_rectangular_faces = [ [0, 1, 2, 3], [1, 0, 4, 5], [2, 1, 5, 6], [3, 2, 6, 7], [0, 3, 7, 4], [4, 7, 6, 5]] _tet_rectangular_faces = [ ] _pyr_rectangular_faces = [ [1, 4, 3, 2]] _prism_rectangular_faces = [ [0, 3, 4, 1], [1, 4, 5, 2], [2, 5, 3, 0]] _tri_face_numbers = { "hex" : _hex_triangular_faces , "tet" : _tet_triangular_faces , "pyr" : _pyr_triangular_faces , "prism" : _prism_triangular_faces } _rect_face_numbers = { "hex" : _hex_rectangular_faces , "tet" : _tet_rectangular_faces , "pyr" : _pyr_rectangular_faces , "prism" : _prism_rectangular_faces } def get_verts_list (self) : """get_verts_list () returns a list of vertex arrays in the form wanted by mesh3. Note: if they are in avs order, they need to be put into Gist order (pyramids and tets are OK, only prisms and hexahedra need vertices permuted). """ verts = [] for k in self._cell_types : if self.cell_dict.has_key (k) : n_z = self.cell_dict[k][1] if n_z.shape [1] == 8 : verts.append ( take (n_z, array ( [0, 1, 3, 2, 4, 5, 7, 6]),axis=0)) elif n_z.shape [1] == 6 : verts.append ( take (n_z, array ( [3, 0, 4, 1, 5, 2]),axis=0)) else : verts.append (n_z) if len (verts) == 1 : verts = verts [0] return verts def create_Narcisse_format ( self ) : """create_Narcisse_format ( ) is an internal routine which takes input in (essentially) avs format and converts it into the format desired by Narcisse. """ if hasattr (self, "cell_descr") : # don't calculate it again. return if not self.cell_dict.has_key ("hex") \ and not self.cell_dict.has_key ("tet") \ and not self.cell_dict.has_key ("prism") \ and not self.cell_dict.has_key ("pyr") : raise _MeshInitError, \ "In avs mode, at least one of the keywords " + \ "(hex, tet, pyr, prism) is required." n_zones = { } nz = { } for k in self._cell_types : if self.cell_dict.has_key (k) : n_zones [k] = self.cell_dict [k][0] nz [k] = self.cell_dict[k][1] else : n_zones [k] = 0 nz [k] = [ ] self.cell_descr = zeros (30 * n_zones ["hex"] + 21 * n_zones ["pyr"] + 16 * n_zones ["tet"] + 23 * n_zones ["prism"], Int) self.number_of_cells = n_zones ["hex"] * 6 + n_zones ["pyr"] * 5 + \ n_zones ["tet"] * 4 + n_zones ["prism"] * 5 n = 0 start = 0 # Fill the start of the cell_descriptor array # The first part is simply a sequence of indices into the second # part of the self.cell_descr array; Each entry tells where the indices # for this particular face begin. cell_j = 0 for k in self._cell_types : nrf = self._no_of_rect_faces [k] ntf = self._no_of_tri_faces [k] for i in range (start, start + n_zones [k]) : for j in range (ntf) : # count triangular faces first n = n + 3 self.cell_descr [cell_j] = n cell_j = cell_j + 1 for j in range (nrf) : # count rectangles next n = n + 4 self.cell_descr [cell_j] = n cell_j = cell_j + 1 start = start + n_zones [k] # Next we go through each type of cell and enter the indices # of its faces, face-by-face, triangular sides first by convention. for k in self._cell_types : nrf = self._no_of_rect_faces [k] ntf = self._no_of_tri_faces [k] n_z = nz [k] for i in range (start, start + n_zones [k]) : n_z_ = n_z [i - start] for j in range (ntf) : # do triangular faces first ind = self._tri_face_numbers [k][j] for l in range (3) : self.cell_descr [cell_j + l] = n_z_ [ind [l]] cell_j = cell_j + 3 for j in range (nrf) : # do rectangular faces next ind = self._rect_face_numbers [k][j] for l in range (4) : self.cell_descr [cell_j + l] = n_z_ [ind [l]] cell_j = cell_j + 4 start = start + n_zones [k] def new ( self , ** keywords ) : """new (...keyword arguments...) allows you to reuse a previously defined mesh. """ del self.x, self.y, self.z, self.c, self.color_card, self.opt_3d, \ self.mask, self.z_c_switch, self.z_contours_scale, \ self.c_contours_scale, self.z_contours_array, \ self.c_contours_array, self.number_of_z_contours, \ self.number_of_c_contours, self.mesh_type if not self.structured : del self.cell_descr, self.number_of_cells self.__init__ ( keywords ) from plane import * class Slice : """ Class Slice is used to contain a geometric description of an isosurface or plane Slice of a Mesh3d object. It is created by the builtin function sslice, defined later in this module. Slice is created as follows: Slice (nv, xyzv, [, val [, plane [, iso [, opt_3d]]]]) where : nv is a one-dimensional integer array whose ith entry is the number of vertices of the ith face. xyzv is a two-dimensional array dimensioned sum (nv,axis=0) by 3. The first nv [0] triples are the vertices of face [0], the next nv [1] triples are the vertices of face [1], etc. val (if present) is an array the same length as nv whose ith entry specifies a color for face i. plane (if present) says that this is a plane Slice, and all the vertices xyzv lie in this plane. iso (if present) says that this is the isosurface for the given value. opt_3d (if present) specifies the 3d options used to color the slice. mask (if present) specifies how the surface is to be masked. contours (if present) is either an array of contour values or an integer giving the number of contours. scale (if present) is the contour scale: "lin", "log", or "normal". edges (if present) is 1 for showing edges, 0 for not. A Slice object or two Slice objects are created by a call to the function sslice (q. v.). The function sslice accepts either a mesh and a specification of how to slice it (isosurface or plane), or else a previously created Slice, a plane to slice it with, and whether you want to keep the resulting "front" Slice or both Slices. """ def type (self) : return Slice3dType def __init__ (self, nv, xyzv, val = None, plane = None, iso = None, smooth = None, opt_3d = ["wm", "f3"], mask = "max", contours = None, scale = "lin", edges = 0) : self.nv = nv self.xyzv = xyzv self.val = val self.plane = plane self.iso = iso self.smooth = smooth if type (opt_3d) != ListType : self.opt_3d = [opt_3d] else : self.opt_3d = opt_3d self.mask = mask self.contours = contours self.scale = scale self.edges = edges self.mesh_type = "xy" self.z_c_switch = 0 self.z_contours_scale = scale self.c_contours_scale = scale if type (self.contours) == IntType or self.contours is None : self.number_of_z_contours = self.contours self.number_of_c_contours = self.contours self.z_contours_array = None self.c_contours_array = None elif type (self.contours) == ArrayType : self.z_contours_array = self.contours self.c_contours_array = self.contours self.number_of_z_contours = len (self.contours) self.number_of_c_contours = len (self.contours) if edges == 0 : self.edges = "w3" in self.opt_3d or \ "w4" in self.opt_3d or \ "wm" in self.opt_3d def __del__ (self) : del self.nv del self.xyzv if self.val is not None: del self.val if self.plane is not None: del self.plane del self.mask if self.iso is not None: del self.iso if self.smooth is not None: del self.smooth del self.opt_3d if self.contours is not None: del self.contours del self.scale del self.edges def new (self, nv, xyzv, val = None, plane = None, iso = None, smooth = None, opt_3d = ["wm", "f3"], mask = "max", contours = None, scale = "lin", edges = 0) : del self.nv del self.xyzv del self.val del self.plane del self.mask del self.iso del self.smooth del self.opt_3d del self.contours del self.scale del self.edges self.__init__ (nv, xyzv, val, plane, iso, smooth, opt_3d, mask, contours, scale, edges) def set ( self , ** keywords ) : """ set (...keyword arguments...) allows you to set individual Slice characteristics. No error checking is done. """ for k in keywords.keys (): setattr (self, k, keywords [k]) if type (self.opt_3d) != ListType : self.opt_3d = [self.opt_3d] if self.edges == 0 : self.edges = "w3" in self.opt_3d or \ "w4" in self.opt_3d or \ "wm" in self.opt_3d _SliceError = "SliceError" def sslice ( v1, v2, varno = 1, nslices = 1, opt_3d = None) : """This function slices a Mesh3d or a previous Slice object, and returns a new Slice object, depending on how it is called. sslice (m, plane, varno = 1) returns a plane Slice through the specified mesh m. varno is the number of the variable to be used to color the Slice (1-based). sslice (m, val, varno = 1) returns an isosurface Slice of m. varno is the number of the variable with respect to which the Slice is taken (1-based) and val is its value on the isosurface. varno can be let default to 1 if there is only one function to worry about. sslice (s, plane, nslices = 1) slices the isosurface s and returns the part in "front" of the plane if nslices == 1, or returns a pair of Slices [front, back] if nslices == 2. If s is a Surface, then this operation will produce a Slice or pair of Slices from the given Surface. N. B. If you want just the "back" surface, you can achieve this by calling sslice with - plane instead of plane. The opt_3d argument can be used to overrule defaults, which are intelligently figured out, but may not always be what you want. """ if v1.type () == Mesh3dType : if type (v1.c) != ListType : funcs = [v1.c] else : funcs = v1.c if not hasattr (v2, "type") or v2.type () != PlaneType : try : val = float (v2) except : raise _SliceError, \ "second argument is not coercible to float." plane = None else : plane = v2 val = None if opt_3d is None : if plane is not None: opt_3d = ["wm", "f4"] else : opt_3d = ["f4"] # First we build a "mesh3" if v1 does not already have one. if v1.structured : #and not hasattr (v1, "m3") : setattr (v1, "m3", mesh3 (v1.x, v1.y, v1.z, funcs = funcs)) elif not v1.structured and v1.avs == 1 : #and not hasattr (v1, "m3") : if no_of_dims (v1.x) == 3 and \ no_of_dims (v1.y) == 3 and \ no_of_dims (v1.z) == 3 : setattr (v1, "m3", mesh3 (array ( [v1.x, v1.y, v1.z], Float), funcs = funcs)) elif no_of_dims (v1.x) == 1 and \ no_of_dims (v1.y) == 1 and \ no_of_dims (v1.z) == 1 : setattr (v1, "m3", mesh3 (v1.x, v1.y, v1.z, funcs = funcs, verts = v1.get_verts_list ())) else : raise _SliceError, \ "x, y, and z must have either 1 or 3 dimensions." elif not hasattr (v1, "m3") : raise _SliceError, \ "Option not implemented." # A Yorick 'mesh3' object now exists, slice it as specified. if plane is None : # An isosurface Slice [nv, xyzv, dum] = slice3 (v1.m3, varno, None, None, value = val) return Slice (nv, xyzv, iso = val, opt_3d = opt_3d) else : # a plane Slice [nv, xyzv, val] = slice3 (v1.m3, plane.rep (), None, None, varno ) return Slice (nv, xyzv, val, plane = plane, opt_3d = opt_3d) elif v1.type () == SurfaceType : plane = v2 if not hasattr (v2, "type") or v2.type () != PlaneType : raise _SliceError, \ "Second argument is not a Plane." if v1.z_c_switch : z = v1.c c = v1.z if "s3" in v1.opt_3d : scale = v1.c_contours_scale if v1.c_contours_array is not None: contours = v1.c_contours_array else : contours = v1.number_of_c_contours elif "s4" in v1.opt_3d : scale = v1.z_contours_scale if v1.z_contours_array is not None: contours = v1.z_contours_array else : contours = v1.number_of_z_contours else : z = v1.z c = v1.c if "s3" in v1.opt_3d : scale = v1.z_contours_scale if v1.z_contours_array is not None: contours = v1.z_contours_array else : contours = v1.number_of_z_contours elif "s4" in v1.opt_3d : scale = v1.c_contours_scale if v1.c_contours_array is not None: contours = v1.c_contours_array else : contours = v1.number_of_c_contours edges = "w3" in v1.opt_3d or "w4" in v1.opt_3d or "wm" in v1.opt_3d \ or "i3" in v1.opt_3d or "i4" in v1.opt_3d if "s4" in v1.opt_3d : [nv, xyzv, col] = slice3mesh (v1.x, v1.y, z, color = c, smooth = 1) elif "f3" in v1.opt_3d : [nv, xyzv, col] = slice3mesh (v1.x, v1.y, z, color = z, smooth = 0) contours = None scale = "lin" elif "f4" in v1.opt_3d : [nv, xyzv, col] = slice3mesh (v1.x, v1.y, z, color = c, smooth = 0) contours = None scale = "lin" else : [nv, xyzv, col] = slice3mesh (v1.x, v1.y, z) contours = None scale = "lin" if nslices == 1 : [nv, xyzv, col] = slice2 (plane.rep (), nv, xyzv, col) return Slice (nv, xyzv, col, opt_3d = v1.opt_3d, mask = v1.mask, contours = contours, scale = scale, edges = edges) elif nslices == 2 : [nv, xyzv, col, nb, xyzb, valb] = slice2x (plane.rep (), nv, xyzv, col) return [Slice (nv, xyzv, col, opt_3d = v1.opt_3d, mask = v1.mask, contours = contours, scale = scale, edges = edges), Slice (nb, xyzb, valb, opt_3d = v1.opt_3d, mask = v1.mask, contours = contours, scale = scale, edges = edges)] else : raise _SliceError, "Illegal number (" + `nslices` + \ ") of slices requested." elif v1.type () == Slice3dType : plane = v2 if not hasattr (v2, "type") or v2.type () != PlaneType : raise _SliceError, \ "Second argument is not a Plane." if nslices == 1 : [nv, xyzv, val] = slice2 (plane.rep (), v1.nv, v1.xyzv, v1.val) return Slice (nv, xyzv, val, plane = v1.plane, iso = v1.iso, opt_3d = v1.opt_3d) elif nslices == 2 : [nf, xyzf, valf, nb, xyzb, valb] = slice2x (plane.rep (), v1.nv, v1.xyzv, v1.val) return [Slice (nf, xyzf, valf, plane = v1.plane, iso = v1.iso, opt_3d = v1.opt_3d), Slice (nb, xyzb, valb, plane = v1.plane, iso = v1.iso, opt_3d = v1.opt_3d)] else : raise _SliceError, "Illegal number (" + `nslices` + \ ") of slices requested." else : raise _SliceError, "Can only slice a Mesh3d or another Slice."