1    | /***************************************
2    |   $Header: /cvsroot/petscgraphics/illuminator.c,v 1.7 2003/05/06 23:15:25 hazelsct Exp $
3    | 
4    |   This is the illuminator.c main file.  It has all of the routines which
5    |   compute the triangulation in a distributed way.
6    | ***************************************/
7    | 
8    | 
9    | #include "config.h" /* esp. for inline */
10   | #include "illuminator.h" /* Just to make sure the interface is "right" */
11   | 
12   | #define MAX_TRIANGLES 2000000 /*+ Maximum number of generated triangles per
13   | 			      node. +*/
14   | 
15   | /*+ Number of triangles in this node.  Its value is initialized to zero, and
16   |   incremented as each triangle is added, then reset to zero when the
17   |   triangulation is displayed. +*/
18   | int num_triangles=0; 
19   | 
20   | /*+ Array of vertex corners, whose size is fixed by MAX_TRIANGLES.  For each
21   |   triangle, this array has the coordinates of the three nodes, and its R, G, B
22   |   and A color values, hence 13 PetscScalars for each triangle. +*/
23   | PetscScalar vertices[13*MAX_TRIANGLES];
24   | 
25   | 
26   | #undef __FUNCT__
27   | #define __FUNCT__ "storetri"
28   | 
29   | /*++++++++++++++++++++++++++++++++++++++
30   |   This little inline routine just implements triangle storage.  Maybe it will
31   |   be more sophisticated in the future.
32   | 
33   |   int storetri Returns 0 or an error code.
34   | 
35   |   PetscScalar x0 X-coordinate of corner 0.
36   | 
37   |   PetscScalar y0 Y-coordinate of corner 0.
38   | 
39   |   PetscScalar z0 Z-coordinate of corner 0.
40   | 
41   |   PetscScalar x1 X-coordinate of corner 1.
42   | 
43   |   PetscScalar y1 Y-coordinate of corner 1.
44   | 
45   |   PetscScalar z1 Z-coordinate of corner 1.
46   | 
47   |   PetscScalar x2 X-coordinate of corner 2.
48   | 
49   |   PetscScalar y2 Y-coordinate of corner 2.
50   | 
51   |   PetscScalar z2 Z-coordinate of corner 2.
52   | 
53   |   PetscScalar *color R,G,B,A quad for this triangle.
54   |   ++++++++++++++++++++++++++++++++++++++*/
55   | 
56   | static inline int storetri
57   | (PetscScalar x0,PetscScalar y0,PetscScalar z0,PetscScalar x1,PetscScalar y1,
58   |  PetscScalar z1,PetscScalar x2,PetscScalar y2,PetscScalar z2,
59   |  PetscScalar *color) 
60   | {
61   |   vertices[13*num_triangles+0] = x0;
62   |   vertices[13*num_triangles+1] = y0;
63   |   vertices[13*num_triangles+2] = z0;
64   |   vertices[13*num_triangles+3] = x1;
65   |   vertices[13*num_triangles+4] = y1;
66   |   vertices[13*num_triangles+5] = z1;
67   |   vertices[13*num_triangles+6] = x2;
68   |   vertices[13*num_triangles+7] = y2;
69   |   vertices[13*num_triangles+8] = z2;
70   |   vertices[13*num_triangles+9] = color[0];
71   |   vertices[13*num_triangles+10] = color[1];
72   |   vertices[13*num_triangles+11] = color[2];
73   |   vertices[13*num_triangles+12] = color[3];
74   | 
75   |   if (++num_triangles >= MAX_TRIANGLES)
76   |     return --num_triangles;
77   |   return 0;
78   | }
79   | 
80   | 
81   | #undef __FUNCT__
82   | #define __FUNCT__ "DrawTetWithPlane"
83   | 
84   | /* Simplifying macro for DrawTetWithPlane */
85   | #define COORD(c1, c2, index) ((index) * (c2) + (1.-(index)) * (c1))
86   | 
87   | /*++++++++++++++++++++++++++++++++++++++
88   |   This function calculates triangle vertices for an isoquant surface in a
89   |   linear tetrahedron, using the whichplane information supplied by the routine
90   |   calling this one, and "draws" them using storetri.  This is really an
91   |   internal function, not intended to be called by user programs.  It is used by
92   |   DrawTet and DrawHex.
93   | 
94   |   int DrawTetWithPlane Returns 0 or an error code.
95   | 
96   |   PetscScalar x0 X-coordinate of vertex 0.
97   | 
98   |   PetscScalar y0 Y-coordinate of vertex 0.
99   | 
100  |   PetscScalar z0 Z-coordinate of vertex 0.
101  | 
102  |   PetscScalar f0 Function value at vertex 0.
103  | 
104  |   PetscScalar x1 X-coordinate of vertex 1.
105  | 
106  |   PetscScalar y1 Y-coordinate of vertex 1.
107  | 
108  |   PetscScalar z1 Z-coordinate of vertex 1.
109  | 
110  |   PetscScalar f1 Function value at vertex 1.
111  | 
112  |   PetscScalar x2 X-coordinate of vertex 2.
113  | 
114  |   PetscScalar y2 Y-coordinate of vertex 2.
115  | 
116  |   PetscScalar z2 Z-coordinate of vertex 2.
117  | 
118  |   PetscScalar f2 Function value at vertex 2.
119  | 
120  |   PetscScalar x3 X-coordinate of vertex 3.
121  | 
122  |   PetscScalar y3 Y-coordinate of vertex 3.
123  | 
124  |   PetscScalar z3 Z-coordinate of vertex 3.
125  | 
126  |   PetscScalar f3 Function value at vertex 3.
127  | 
128  |   PetscScalar isoquant Function value at which to draw triangle.
129  | 
130  |   PetscScalar edge0 Normalized intercept at edge 0, 0. at node 0, 1. at node 1.
131  | 
132  |   PetscScalar edge1 Normalized intercept at edge 1, 0. at node 1, 1. at node 2.
133  | 
134  |   PetscScalar edge3 Normalized intercept at edge 3, 0. at node 0, 1. at node 3.
135  | 
136  |   int whichplane Index of which edge intercept(s) is between zero and 1.
137  | 
138  |   PetscScalar *color R,G,B,A quad for this tetrahedron.
139  |   ++++++++++++++++++++++++++++++++++++++*/
140  | 
141  | static inline int DrawTetWithPlane
142  | (PetscScalar x0,PetscScalar y0,PetscScalar z0,PetscScalar f0,
143  |  PetscScalar x1,PetscScalar y1,PetscScalar z1,PetscScalar f1,
144  |  PetscScalar x2,PetscScalar y2,PetscScalar z2,PetscScalar f2,
145  |  PetscScalar x3,PetscScalar y3,PetscScalar z3,PetscScalar f3,
146  |  PetscScalar isoquant, PetscScalar edge0,PetscScalar edge1,PetscScalar edge3,
147  |  int whichplane, PetscScalar *color)
148  | {
149  |   PetscScalar edge2,edge4,edge5;
150  | 
151  |   switch (whichplane) {
152  |   case 7: { /* Triangles on edges 0,1,3; 3,1,5 */
153  |     edge5 = (isoquant - f2) / (f3 - f2);
154  | 
155  |     /* Use edge 0 direction to guarantee right-handedness */
156  |     if (f1 > f0) {
157  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
158  | 		COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
159  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
160  | 		color);
161  |       storetri (COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
162  | 		COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
163  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
164  | 		color);
165  |     }
166  |     else {
167  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
168  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
169  | 		COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
170  | 		color);
171  |       storetri (COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
172  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
173  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
174  | 		color);
175  |     }
176  |     break;
177  |   }
178  | 
179  |   case 6: { /* Triangles on edges 1,2,4; 4,2,3 */
180  |     edge2 = (isoquant - f2) / (f0 - f2);
181  |     edge4 = (isoquant - f1) / (f3 - f1);
182  | 
183  |     /* Use edge 1 direction to guarantee right-handedness */
184  |     if (f2 > f1) {
185  |       storetri (COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
186  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
187  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
188  | 		color);
189  |       storetri (COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
190  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
191  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
192  | 		color);
193  |     }
194  |     else {
195  |       storetri (COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
196  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
197  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
198  | 		color);
199  |       storetri (COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
200  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
201  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
202  | 		color);
203  |     }
204  |     break;
205  |   }
206  | 
207  |   case 5: { /* Triangles on edges 0,2,3 */
208  |     edge2 = (isoquant - f2) / (f0 - f2);
209  | 
210  |     /* Use edge 0 direction to guarantee right-handedness */
211  |     if (f1 > f0)
212  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
213  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
214  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
215  | 		color);
216  |     else
217  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
218  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
219  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
220  | 		color);
221  |     break;
222  |   }
223  | 
224  |   case 4: { /* Triangles on edges 3,4,5 */
225  |     edge4 = (isoquant - f1) / (f3 - f1);
226  |     edge5 = (isoquant - f2) / (f3 - f2);
227  | 
228  |     /* Use edge 3 direction to guarantee right-handedness */
229  |     if (f3 > f0)
230  |       storetri (COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
231  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
232  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
233  | 		color);
234  |     else
235  |       storetri (COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
236  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
237  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
238  | 		color);
239  |     break;
240  |   }
241  | 
242  |   case 3: { /* Triangles on edges 0,1,4 */
243  |     edge4 = (isoquant - f1) / (f3 - f1);
244  | 
245  |     /* Use edge 0 direction to guarantee right-handedness */
246  |     if (f1 > f0)
247  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
248  | 		COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
249  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
250  | 		color);
251  |     else
252  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
253  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
254  | 		COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
255  | 		color);
256  |     break;
257  |   }
258  | 
259  |   case 2: { /* Triangles on edges 1,2,5 */
260  |     edge2 = (isoquant - f2) / (f0 - f2);
261  |     edge5 = (isoquant - f2) / (f3 - f2);
262  | 
263  |     /* Use edge 1 direction to guarantee right-handedness */
264  |     if (f2 > f1)
265  |       storetri (COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
266  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
267  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
268  | 		color);
269  |     else
270  |       storetri (COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
271  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
272  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
273  | 		color);
274  |     break;
275  |   }
276  | 
277  |   case 1: { /* Triangles on edges 0,2,4; 4,2,5 */
278  |     edge2 = (isoquant - f2) / (f0 - f2);
279  |     edge4 = (isoquant - f1) / (f3 - f1);
280  |     edge5 = (isoquant - f2) / (f3 - f2);
281  | 
282  |     /* Use edge 0 direction to guarantee right-handedness */
283  |     if (f1 > f0) {
284  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
285  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
286  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
287  | 		color);
288  |       storetri (COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
289  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
290  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
291  | 		color);
292  |     }
293  |     else {
294  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
295  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
296  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
297  | 		color);
298  |       storetri (COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
299  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
300  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
301  | 		color);
302  |     }
303  |     break;
304  |   }
305  |   }
306  |   return 0;
307  | }
308  | #undef COORD
309  | 
310  | 
311  | #undef __FUNCT__
312  | #define __FUNCT__ "DrawTet"
313  | 
314  | /*++++++++++++++++++++++++++++++++++++++
315  |   This sets the edge and whichplane parameters and then passes everything to
316  |   DrawTetWithPlane to actually draw the triangle.  It is intended for use by
317  |   developers with distributed arrays based on tetrahedra, e.g. a finite element
318  |   mesh.
319  | 
320  |   int DrawTet Returns 0 or an error code.
321  | 
322  |   PetscScalar *coords Coordinates of tetrahedron corner points: x0, y0, z0, x1, etc.
323  | 
324  |   PetscScalar *vals Function values at tetrahedron corners: f0, f1, f2, f3.
325  | 
326  |   PetscScalar isoquant Function value at which to draw triangle.
327  | 
328  |   PetscScalar *color R,G,B,A quad for this tetrahedron.
329  |   ++++++++++++++++++++++++++++++++++++++*/
330  | 
331  | int DrawTet
332  | (PetscScalar *coords, PetscScalar *vals, PetscScalar isoquant,
333  |  PetscScalar *color)
334  | {
335  |   PetscScalar edge0, edge1, edge3;
336  |   int whichplane, ierr=0;
337  | 
338  |   edge0 = (isoquant-vals[0]) / (vals[1]-vals[0]);
339  |   edge1 = (isoquant-vals[1]) / (vals[2]-vals[1]);
340  |   edge3 = (isoquant-vals[0]) / (vals[3]-vals[0]);
341  |   whichplane = (edge0>0. && edge0<1.) | ((edge1>0. && edge1<1.) << 1) |
342  |     ((edge3>0. && edge3<1.) << 2);
343  |   if (whichplane)
344  |     ierr = DrawTetWithPlane (coords[0],coords[1],coords[2],vals[0],
345  | 			     coords[3],coords[4],coords[5],vals[1],
346  | 			     coords[6],coords[7],coords[8],vals[2],
347  | 			     coords[9],coords[10],coords[11],vals[3],
348  | 			     isoquant, edge0,edge1,edge3, whichplane, color);
349  |   return ierr;
350  | }
351  | 
352  | 
353  | #undef __FUNCT__
354  | #define __FUNCT__ "DrawHex"
355  | 
356  | /*++++++++++++++++++++++++++++++++++++++
357  |   This divides a right regular hexahedron into tetrahedra, and loops over them
358  |   to generate triangles on each one.  It calculates edge and whichplane
359  |   parameters so it can use DrawTetWithPlane directly.
360  | 
361  |   int DrawHex Returns 0 or an error code.
362  | 
363  |   PetscScalar *coords Coordinates of hexahedron corner points: xmin, xmax, ymin, etc.
364  | 
365  |   PetscScalar *vals Function values at hexahedron corners: f0, f1, f2, etc.
366  | 
367  |   PetscScalar isoquant Function value at which to draw triangles.
368  | 
369  |   PetscScalar *color R,G,B,A quad for this hexahedron.
370  |   ++++++++++++++++++++++++++++++++++++++*/
371  | 
372  | int DrawHex
373  | (PetscScalar *coords, PetscScalar *vals, PetscScalar isoquant,
374  |  PetscScalar *color)
375  | {
376  |   int tetras[6][4] = {{0,1,2,4}, {1,2,4,5}, {2,4,5,6}, {1,3,2,5}, {2,5,3,6}, {3,6,5,7}};
377  |   int tet,node,ierr;
378  |   int c0,c1,c2,c3,whichplane;
379  |   PetscScalar edge0,edge1,edge3;
380  | 
381  |   for(tet=0; tet<6; tet++)
382  |     {
383  |       /* Within a tetrahedron, edges 0 through 5 connect corners:
384  | 	 0,1; 1,2; 2,0; 0,3; 1,3; 2,3 respectively */
385  |       c0 = tetras[tet][0];
386  |       c1 = tetras[tet][1];
387  |       c2 = tetras[tet][2];
388  |       c3 = tetras[tet][3];
389  |       edge0 = (isoquant-vals[c0]) / (vals[c1]-vals[c0]);
390  |       edge1 = (isoquant-vals[c1]) / (vals[c2]-vals[c1]);
391  |       edge3 = (isoquant-vals[c0]) / (vals[c3]-vals[c0]);
392  |       whichplane = (edge0>0. && edge0<1.) | ((edge1>0. && edge1<1.) << 1) |
393  | 	((edge3>0. && edge3<1.) << 2);
394  |       if (whichplane)
395  | 	{
396  | 	  ierr=DrawTetWithPlane
397  | 	    (coords[c0&1],coords[2+((c0&2)>>1)],coords[4+((c0&4)>>2)],vals[c0],
398  | 	     coords[c1&1],coords[2+((c1&2)>>1)],coords[4+((c1&4)>>2)],vals[c1],
399  | 	     coords[c2&1],coords[2+((c2&2)>>1)],coords[4+((c2&4)>>2)],vals[c2],
400  | 	     coords[c3&1],coords[2+((c3&2)>>1)],coords[4+((c3&4)>>2)],vals[c3],
401  | 	     isoquant, edge0,edge1,edge3, whichplane, color); CHKERRQ (ierr);
402  | 	}
403  |     }
404  | 
405  |   return 0;
406  | }
407  | 
408  | 
409  | #undef __FUNCT__
410  | #define __FUNCT__ "Draw3DBlock"
411  | 
412  | /*++++++++++++++++++++++++++++++++++++++
413  |   Calculate vertices of isoquant triangle(s) in a 3-D array of right regular
414  |   hexahedra.  This loops through a 3-D array and calls DrawHex to calculate the
415  |   triangulation of each hexahedral cell.
416  | 
417  |   int Draw3DBlock Returns 0 or an error code.
418  | 
419  |   int xd Overall x-width of function value array.
420  | 
421  |   int yd Overall y-width of function value array.
422  | 
423  |   int zd Overall z-width of function value array.
424  | 
425  |   int xs X-index of the start of the array section we'd like to draw.
426  | 
427  |   int ys Y-index of the start of the array section we'd like to draw.
428  | 
429  |   int zs Z-index of the start of the array section we'd like to draw.
430  | 
431  |   int xm X-width of the array section we'd like to draw.
432  | 
433  |   int ym Y-width of the array section we'd like to draw.
434  | 
435  |   int zm Z-width of the array section we'd like to draw.
436  | 
437  |   PetscScalar *minmax Position of block corners: xmin, xmax, ymin, ymax, zmin, zmax.
438  | 
439  |   PetscScalar *vals The array of function values at vertices.
440  | 
441  |   int skip Number of interlaced fields in this array.
442  | 
443  |   int n_quants Number of isoquant surfaces to draw (isoquant values).
444  | 
445  |   PetscScalar *isoquants Array of function values at which to draw triangles.
446  | 
447  |   PetscScalar *colors Array of color R,G,B,A quads for each isoquant.
448  |   ++++++++++++++++++++++++++++++++++++++*/
449  | 
450  | int Draw3DBlock
451  | (int xd,int yd,int zd, int xs,int ys,int zs, int xm,int ym,int zm,
452  |  PetscScalar *minmax, PetscScalar *vals, int skip, int n_quants,
453  |  PetscScalar *isoquants, PetscScalar *colors)
454  | {
455  |   int i,j,k,q, start, ierr;
456  |   PetscScalar boxcoords[6], funcs[8];
457  | 
458  |   /* Simple check */
459  |   if (xd<=xs+xm || yd<=ys+ym || zd<=zs+zm)
460  |     IllErrorHandler (PETSC_ERR_ARG_WRONGSTATE, "Array size mismatch");
461  | 
462  |   /* The big loop */
463  |   for (k=0; k<zm; k++)
464  |     {
465  |       boxcoords[4] = minmax[4] + (minmax[5]-minmax[4])*k/zm;
466  |       boxcoords[5] = minmax[4] + (minmax[5]-minmax[4])*(k+1)/zm;
467  |       for(j=0;j<ym;j++)
468  | 	{
469  | 	  boxcoords[2] = minmax[2] + (minmax[3]-minmax[2])*j/ym;
470  | 	  boxcoords[3] = minmax[2] + (minmax[3]-minmax[2])*(j+1)/ym;
471  | 	  for(i=0;i<xm;i++)
472  | 	    {
473  | 	      boxcoords[0] = minmax[0] + (minmax[1]-minmax[0])*i/xm;
474  | 	      boxcoords[1] = minmax[0] + (minmax[1]-minmax[0])*(i+1)/xm;
475  | 	      start = ((k+zs)*yd + j+ys)*xd + xs + i;
476  | 	      funcs[0] = vals[skip*start];
477  | 	      funcs[1] = vals[skip*(start+1)];
478  | 	      funcs[2] = vals[skip*(start+xd)];
479  | 	      funcs[3] = vals[skip*(start+xd+1)];
480  | 	      funcs[4] = vals[skip*(start+xd*yd)];
481  | 	      funcs[5] = vals[skip*(start+xd*yd+1)];
482  | 	      funcs[6] = vals[skip*(start+xd*yd+xd)];
483  | 	      funcs[7] = vals[skip*(start+xd*yd+xd+1)];
484  | 	      for (q=0; q<n_quants; q++)
485  | 		{
486  | 		  ierr = DrawHex (boxcoords, funcs, isoquants[q], colors+4*q);
487  | 		  CHKERRQ (ierr);
488  | 		}
489  | 	    }
490  | 	}
491  |     }
492  | 
493  |   return 0;
494  | }