The reader should be familiar with many of the functions in gistCmodule from CHAPTER 5: "Two-Dimensional Plotting Functions", page 27. Most of the gistCmodule functions discussed there, and many of the helper functions, are pretty close to literal translations of the equivalent functions in Gist, the main difference geing the superstructure built on top of them in order to handle PyObjects.. In addition to the plotting functions, a number of functions in gistCmodule are essential for maintenance of that module, and are discussed here. We also discuss briefly a few other functions which are not literal translations of Gist functions.

One of the primary challenges facing the developer of Python extensions is correct management of the reference counting for Python objects. Memory leaks will result if the programmer fails to decrement the reference count in temporary objects. On the other hand, decrementing the reference count too early can cause an object to go away that is referred to later, which can cause a segmentation fault when it is referenced. We have semi-automated the process in gistCmodule by maintaining a list of all PyObjects created in the process of running one of the module's functions. This list is declared as follows:

- #define ARRAY_LIST_SIZE 30
- static PyObject * PyArrayList [ARRAY_LIST_SIZE];
- static int array_list_length = 0;

- static int addToArrayList (PyObject * obj)
- static void clearArrayList ()
- static void removeFromArrayList (PyObject * obj)
- static void takeOffArrayList (PyObject * obj)

addToArrayList places obj on PyArrayList and returns 1 if successful. If obj is NULL or the list is full, returns 0. clearArrayList DECREF's everything on the list, and sets array_list_length to 0. This needs to be done prior to any error return. removeFromArrayList DECREF's obj (if it is on the list), removes it from the list, and compresses the list. takeOffArrayList removes obj from the list and compresses the list, but does not DECREF obj. This is done, for example, when obj is to be returned to the caller.

addToArrayList occurs throughout gistCmodule primarily in macros which create arrays. All of these macros use the TRY macro, which is defined as follows:

- #define TRY(e, m) do{if(!(e)){clearArrayList(); \
- clearFreeList(0);clearMemList();return m;}} while(0)
The idea behind TRY is that generally Python functions return 0 or NULL if an error occurred. In this case it is necessary to get rid of all temporary objects and memory which was allocated up to this point. clearArrayList was discussed above, clearFreeList and clearMemList are discussed later in the chapter.

- #define GET_ARR(ap,op,type,dim,cast) \
- TRY(addToArrayList((PyObject *)(ap=(PyArrayObject *)\
- PyArray_ContiguousFromObject(op,type,dim,dim))), \
- (cast)PyErr_NoMemory ())
This macro is the usual protocol for creating a contiguous array from a PyObject which has been sent as an argument to a function.

- #define NEW_ARR(ap,n,dims,type,cast) \
- TRY(addToArrayList((PyObject *)(ap=\
- (PyArrayObject *)PyArray_FromDims(n,dims,type))), \
- (cast)PyErr_NoMemory ())
This macro is used usually when creating an array whose dimensions are known and which is to be filled with computed data.

- #define RET_ARR(op,ndim,dim,type,data,cast)\
- TRY(addToArrayList(op=\
- PyArray_FromDimsAndData(ndim,dim,type,data)), \
- (cast)PyErr_NoMemory ())
This final macro is used when we have a block of data and we wish to create an array containing this data, usually as a return value from a function. In order to keep this object from being permanent, use the following macro:

- #define SET_OWN(op) \
- ( (PyArrayObject *) op)->flags |= OWN_DATA
This macro sets a flag in the PyObject which tells Python that it can be DECREF'ed.

ArrayObjects are defined as follows:

- typedef struct arrayobject {
- void * data ;
- int size ;
- char typecode ;
- } ArrayObject;
These objects are used primarily in the slice2 routines to store temporary results during the calculation. The final results are passed back in PyArrayObjects created by RET_ARR. Two lists of ArrayObjects are maintained by the slice2 suite: list 0 by slice2 itself, and list 1, which is used by _slice2_part, which is called by slice2. These lists are declared as follows:

- #define MAX_NO_LISTS 2
- #define MAX_LIST_SIZE 30
- static ArrayObject * freeList [MAX_NO_LISTS] [MAX_LIST_SIZE];
- static int freeListLen [MAX_NO_LISTS] = {0, 0};

- static ArrayObject * allocateArray (int size, char tc,
- int nlist)
- static ArrayObject * copyArray (ArrayObject * a)
- static ArrayObject * arrayFromPointer (int size, char tc,
- void * data, int nlist)
- static void freeArray (ArrayObject * a, int n)
- static void clearFreeList (int n)
- static int addToFreeList (ArrayObject * x, int n)
- static void removeArrayOnly (ArrayObject * x, int n)
- static void removeFromFreeList (ArrayObject * x, int n)

allocateArray allocates an appropriate amount of space for size items of type tc. It then creates an ArrayObject containing this data and puts it on freeList [nlist]. copyArray makes and returns a copy of a. It does not add a to any freeList. arrayFromPointer creates an arrayObject whose data pointer points to data; it is assumed that the caller has supplied correct size and tc arguments. The resulting object is placed on freeList [nlist]. freeArray frees a's data and then a itself, and removes it from freeList [n] if it is there. clearFreeList frees everything on freeList [n] and sets the list length to 0. addToFreeList adds x to freeList [n], if it can. removeArrayOnly removes the array from freeList [n], then frees x without freeing its data. This would most likely be done when RET_ARR creates a PyArrayObject which points to x's data. removeFromFreeList frees x's data, then x itself, and removes x from freeList [n].

Occasionally in gistCmodule it is necessary to malloc a block of memory which is not contained inside some type of object. MemList is used to keep track of such memory:

- #define MEM_LIST_SIZE 15
- static void * PyMemList [MEM_LIST_SIZE];
- static int mem_list_length = 0;

- static int addToMemList (void * addr)
- static void clearMemList ()

The first function adds an address to MemList; the second frees everything on MemList and sets its length back to 0.

- # Set mesh first
- plmesh (y, x, ireg, triangle = triangle)
- [nc, yc, xc] = contour (level, z)

The calling sequence given above emphasizes that mesh parameters should be set by a call to plmesh prior to calling contour. plmesh arguments are explained in section "plmesh: Set Default Mesh" on page 29. If level is a scalar floating point number, then the the returned values are the points at that contour level. All such points lie on edges of the mesh. If a contour curve closes, the final point is the same as the initial point (i.e., that point is included twice in the returned list). If level is a sequence of two reals, then contour returns the points of a set of polygons which outline the regions between the two contour levels. The returned values are in the form required for arguments of plfp (see Section 5.1.8 "plfp: Plot a List of Filled Polygons" on page 39).These will include points on the mesh boundary which lie between the levels, in addition to the edge points for both levels. The polygons are closed, simply connected, and will not contain more than about 4000 points (larger polygons are split into pieces with a few points repeated where the pieces join).

The 2D filled contour plot routine plfc (see Section 5.1.7 "plfc: Plot filled contours" on page 37) operates by calling contour with pairs of adjacent contour levels, and then calling plfp with the output and a single color, inside a loop. contour needed to be programmed in C because it can be called many times to do a single filled contour plot, and the calculations take too long to be performed in interpreted code. contour calls lower level Gist routines that do most of the work.

The 3D graphics in Gist itself is still experimental, and virtually all the computational functions are written in Yorick, an interpreted language. Many of the PyGist 3D computations were translated into Python from Yorick, originally including the slice2 and slice2x functions, and their auxiliary, _slice2_part. When we impleoemnted contours and filled contours on surfaces (which is currently not implemented in Gist itself), we used slice2 and slice2x to compute the contours and the polygon lists enclosed within them. These computations were much too slow, so we rewrote slice2 in C (slice2x remains in Python; it just calls slice2 with a parameter set) and put them into the gistCmodule. The user interface to these functions has been discussed in a previous section (7.4.3 "slice2 and slice2x: Slicing Surfaces with planes"), but we discuss them here from the viewpoint of implementation. We also discuss the ``hidden'' function _slice2_part here for the first time.

- [nverts, xyzverts, values] = slice2 (plane, nv,
- xyzv, vals = None, _slice2x = 0)
- [nverts, xyzverts, values, nvertb, xyzvertb, valueb] =
- slice2x (plane, nv, xyzv, vals)
- static int _slice2_part (ArrayObject * xyzc,
- ArrayObject * keep, ArrayObject * next, ArrayObject * dp,
- ArrayObject * prev, ArrayObject * last,
- ArrayObject * valc, ArrayObject ** xyzc_new,
- ArrayObject ** nvertc, ArrayObject ** valc_new,
- int freexyzc, int freevalc)

The argument plane can be either a scalar or a plane3 (see "Creating a Plane" on page 61); nv is an array of integers, the i^{th} entry of which gives the number of vertices of the i^{th} polygonal cell; xyzv are the vertices of the coordinatesof the cells, with each consecutive nv [i] entries representing the vertices of the i^{th} cell; and vals being a set of values, one for each cell. These arguments are the same format as returned by slice3 and slice3mesh.

If plane is a plane3, then vals (if not None) is a cell-centered set of values expressing the color of each cell, and the outputs nverts, xyzverts, and values represent the polygons and their colors (if any) describing the portion of the sliced surface that is on the positive side of the plane. That's all you get with slice2. With slice2x, you get in addition nvertb, xyzvertb, and valueb, which describe the part of the surface on the negative side of the slicing plane. Warning: one of these specifications could be None, None, None if the entire surface lies on one side of the plane.

If plane is a scalar value, then vals must be present and must be node-centered. In this case, the outputs nverts, xyzverts, and values represent the polygons and their colors (if any) describing the portion of the sliced surface where vals on the vertices are greater than or equal to the scalar value plane. (This actually allows you to form an arbitrary two-dimensional slice of a surface.) With slice2x, you get in addition nvertb, xyzvertb, and valueb, which describe the part of the surface where vals on the vertices are less than the scalar value plane.

The optional parameter _slice2x, if 1, tells slice2 to return slices on both sides of the slicing surface or plane; if not present, or 0, then the slice on ``top'' is returned. slice2 works by deciding which polygons lie entirely ``above'' the slicing surface, which ones lie entirely ``below'' the slicing surface, and which ones are cut by the surface. If _slice2x is 0, then the ones ``below'' the surface are discarded. slice2 then calls _slice2_part with the polygons to be cut by the plane; once to get the cut polygons ``above'' the surface, then, if _slice2x is 1, a second time to get the cut polygons ``below'' the surface. The list of uncut and cut polygons ``above'' the surface is concatenated and returned (_slice2x == 0); the list of uncut and cut polygons ``below'' the surface is concatenated and returned also if _slice2x is 1.

In the case of a plane slice, suppose that the equation of the slicing plane is

Then a point (*x*1*, y*1*, z*1) is considered to be on the positive side of the plane if

*ax*1 + *by*1 + *cz*1 - *d* >= *_slice2_precision*

*ax*1 + *by*1 + *cz*1 - *d* < *_slice2_precision*

For a discussion of _slice2_precision, and how to get and set its value, see Section 8.12 "Controlling Points Close to the Slicing Plane: _slice2_precision" on page 92.

In the case of a slicing surface, vertex i is considered to be above the surface if

*vals* *[i]* - *plane* >= *_slice2_precision*

*vals* *[i]* - *plane* < *_slice2_precision*

For all intents and purposes, the user may assume that _slice2_precision is 0.0, as this is the default. However, we allow you to change this if you think you have good reason.

There is a conceptual difficulty for the case of a quad face all four of whose edges are cut by the slicing plane or surface. 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 we don't believe is worthwhile, because the X choice is unique, and we don't see why we shouldn't use it here. For contouring routines, we reject the X choice on aesthetic grounds, and perhaps that will prove to be the case here as well - but we 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 we doubt there will be an aesthetic problem. Anyway, our implementation of slice3, slice2, and _slice2_part 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 aesthetically as every zone triangulated the same way!

[Top] [Prev] [Next] [Bottom]

support@icf.llnl.gov