1    | /***************************************
2    |   $Header$
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   | #include "structures.h"     /* For the isurface structure definition */
12   | #include <stdlib.h>      /* For malloc() */
13   | 
14   | /* Build with -DDEBUG for debugging output */
15   | #undef DPRINTF
16   | #ifdef DEBUG
17   | #define DPRINTF(fmt, args...) PetscPrintf (PETSC_COMM_WORLD, "%s: " fmt, __FUNCT__, args)
18   | #else
19   | #define DPRINTF(fmt, args...)
20   | #endif
21   | 
22   | #define MAX_TRIANGLES 10000000 /*+ Maximum number of generated triangles per
23   | 				node. +*/
24   | 
25   | 
26   | #undef __FUNCT__
27   | #define __FUNCT__ "ISurfCreate"
28   | 
29   | /*++++++++++++++++++++++++++++++++++++++
30   |   Constructor for ISurface data object.
31   | 
32   |   int SurfCreate Returns zero or an error code.
33   | 
34   |   Surf *newsurf Address in which to put the new ISurface object.
35   |   ++++++++++++++++++++++++++++++++++++++*/
36   | 
37   | int SurfCreate (ISurface *newsurf)
38   | {
39   |   struct isurface *thesurf;
40   | 
41   |   if (!(thesurf = (struct isurface *) malloc (sizeof (struct isurface))))
42   |     SETERRQ (PETSC_ERR_MEM, "Unable to allocate memory for isurface object");
43   | 
44   |   thesurf->vertices = NULL;
45   |   thesurf->vertisize = thesurf->num_triangles = 0;
46   |   *newsurf = (ISurface) thesurf;
47   |   return 0;
48   | }
49   | 
50   | 
51   | #undef __FUNCT__
52   | #define __FUNCT__ "ISurfDestroy"
53   | 
54   | /*++++++++++++++++++++++++++++++++++++++
55   |   Destructor for ISurface data object.
56   | 
57   |   int ISurfDestroy Returns zero or an error code.
58   | 
59   |   ISurface Surf ISurface object to destroy.
60   |   ++++++++++++++++++++++++++++++++++++++*/
61   | 
62   | int ISurfDestroy (ISurface Surf)
63   | {
64   |   struct isurface *thesurf = (struct isurface *) Surf;
65   | 
66   |   if (thesurf->vertices)
67   |     free (thesurf->vertices);
68   |   free (thesurf);
69   |   return 0;
70   | }
71   | 
72   | 
73   | #undef __FUNCT__
74   | #define __FUNCT__ "ISurfClear"
75   | 
76   | /*++++++++++++++++++++++++++++++++++++++
77   |   Clear out an isurface data object without freeing its memory.
78   | 
79   |   int SurfClear Returns zero or an error code.
80   | 
81   |   ISurface Surf ISurface object to clear out.
82   |   ++++++++++++++++++++++++++++++++++++++*/
83   | 
84   | int SurfClear (ISurface Surf)
85   | {
86   |   struct isurface *thesurf = (struct isurface *) Surf;
87   | 
88   |   return thesurf->num_triangles = 0;
89   | }
90   | 
91   | 
92   | #undef __FUNCT__
93   | #define __FUNCT__ "storetri"
94   | 
95   | /*++++++++++++++++++++++++++++++++++++++
96   |   This little inline routine just implements triangle storage.
97   | 
98   |   int storetri Returns 0 or an error code.
99   | 
100  |   PetscScalar x0 X-coordinate of corner 0.
101  | 
102  |   PetscScalar y0 Y-coordinate of corner 0.
103  | 
104  |   PetscScalar z0 Z-coordinate of corner 0.
105  | 
106  |   PetscScalar x1 X-coordinate of corner 1.
107  | 
108  |   PetscScalar y1 Y-coordinate of corner 1.
109  | 
110  |   PetscScalar z1 Z-coordinate of corner 1.
111  | 
112  |   PetscScalar x2 X-coordinate of corner 2.
113  | 
114  |   PetscScalar y2 Y-coordinate of corner 2.
115  | 
116  |   PetscScalar z2 Z-coordinate of corner 2.
117  | 
118  |   PetscScalar *color R,G,B,A quad for this triangle.
119  | 
120  |   struct isurface *thesurf ISurface object to put the triangles in.
121  |   ++++++++++++++++++++++++++++++++++++++*/
122  | 
123  | static inline int storetri
124  | (PetscScalar x0,PetscScalar y0,PetscScalar z0,
125  |  PetscScalar x1,PetscScalar y1,PetscScalar z1,
126  |  PetscScalar x2,PetscScalar y2,PetscScalar z2, PetscScalar *color,
127  |  struct isurface *thesurf) 
128  | {
129  |   if (thesurf->num_triangles >= thesurf->vertisize)
130  |     {
131  |       thesurf->vertices = (PetscScalar *) realloc
132  | 	(thesurf->vertices, 13*sizeof(double)*
133  | 	 (thesurf->vertisize = (thesurf->vertisize==0) ? 1000 :
134  | 	  ((thesurf->vertisize > MAX_TRIANGLES) ?
135  | 	   thesurf->vertisize : thesurf->vertisize*2)));
136  |       if (thesurf->num_triangles >= thesurf->vertisize)
137  | 	{
138  | 	  int rank, ierr;
139  | 	  ierr = MPI_Comm_rank (PETSC_COMM_WORLD, &rank); CHKERRQ (ierr);
140  | 	  printf ("CPU %d: Number of triangles exceeded %d, emptying array.\n"
141  | 		  "(The limit is MAX_TRIANGLES in illuminator.c line 21.)",
142  | 		  rank, MAX_TRIANGLES);
143  | 	  thesurf->num_triangles = 0;
144  | 	}
145  |       else
146  | 	DPRINTF ("Expanding vertices array to %d triangles, pointer at 0x%lx\n",
147  | 		 thesurf->vertisize, thesurf->vertices);
148  |     }
149  | 
150  |   thesurf->vertices[13*thesurf->num_triangles+0] = x0;
151  |   thesurf->vertices[13*thesurf->num_triangles+1] = y0;
152  |   thesurf->vertices[13*thesurf->num_triangles+2] = z0;
153  |   thesurf->vertices[13*thesurf->num_triangles+3] = x1;
154  |   thesurf->vertices[13*thesurf->num_triangles+4] = y1;
155  |   thesurf->vertices[13*thesurf->num_triangles+5] = z1;
156  |   thesurf->vertices[13*thesurf->num_triangles+6] = x2;
157  |   thesurf->vertices[13*thesurf->num_triangles+7] = y2;
158  |   thesurf->vertices[13*thesurf->num_triangles+8] = z2;
159  |   thesurf->vertices[13*thesurf->num_triangles+9] = color[0];
160  |   thesurf->vertices[13*thesurf->num_triangles+10] = color[1];
161  |   thesurf->vertices[13*thesurf->num_triangles+11] = color[2];
162  |   thesurf->vertices[13*thesurf->num_triangles+12] = color[3];
163  | 
164  |   thesurf->num_triangles++;
165  |   return 0;
166  | }
167  | 
168  | 
169  | #undef __FUNCT__
170  | #define __FUNCT__ "DrawTetWithPlane"
171  | 
172  | /* Simplifying macro for DrawTetWithPlane */
173  | #define COORD(c1, c2, index) ((index) * (c2) + (1.-(index)) * (c1))
174  | 
175  | /*++++++++++++++++++++++++++++++++++++++
176  |   This function calculates triangle vertices for an isoquant surface in a
177  |   linear tetrahedron, using the whichplane information supplied by the routine
178  |   calling this one, and "draws" them using storetri().  This is really an
179  |   internal function, not intended to be called by user programs.  It is used by
180  |   DrawTet() and DrawHex().
181  | 
182  |   int DrawTetWithPlane Returns 0 or an error code.
183  | 
184  |   PetscScalar x0 X-coordinate of vertex 0.
185  | 
186  |   PetscScalar y0 Y-coordinate of vertex 0.
187  | 
188  |   PetscScalar z0 Z-coordinate of vertex 0.
189  | 
190  |   PetscScalar f0 Function value at vertex 0.
191  | 
192  |   PetscScalar x1 X-coordinate of vertex 1.
193  | 
194  |   PetscScalar y1 Y-coordinate of vertex 1.
195  | 
196  |   PetscScalar z1 Z-coordinate of vertex 1.
197  | 
198  |   PetscScalar f1 Function value at vertex 1.
199  | 
200  |   PetscScalar x2 X-coordinate of vertex 2.
201  | 
202  |   PetscScalar y2 Y-coordinate of vertex 2.
203  | 
204  |   PetscScalar z2 Z-coordinate of vertex 2.
205  | 
206  |   PetscScalar f2 Function value at vertex 2.
207  | 
208  |   PetscScalar x3 X-coordinate of vertex 3.
209  | 
210  |   PetscScalar y3 Y-coordinate of vertex 3.
211  | 
212  |   PetscScalar z3 Z-coordinate of vertex 3.
213  | 
214  |   PetscScalar f3 Function value at vertex 3.
215  | 
216  |   PetscScalar isoquant Function value at which to draw triangle.
217  | 
218  |   PetscScalar edge0 Normalized intercept at edge 0, 0. at node 0, 1. at node 1.
219  | 
220  |   PetscScalar edge1 Normalized intercept at edge 1, 0. at node 1, 1. at node 2.
221  | 
222  |   PetscScalar edge3 Normalized intercept at edge 3, 0. at node 0, 1. at node 3.
223  | 
224  |   int whichplane Index of which edge intercept(s) is between zero and 1.
225  | 
226  |   PetscScalar *color R,G,B,A quad for this tetrahedron.
227  | 
228  |   struct isurface *thesurf ISurface object to put the triangles in.
229  |   ++++++++++++++++++++++++++++++++++++++*/
230  | 
231  | static inline int DrawTetWithPlane
232  | (PetscScalar x0,PetscScalar y0,PetscScalar z0,PetscScalar f0,
233  |  PetscScalar x1,PetscScalar y1,PetscScalar z1,PetscScalar f1,
234  |  PetscScalar x2,PetscScalar y2,PetscScalar z2,PetscScalar f2,
235  |  PetscScalar x3,PetscScalar y3,PetscScalar z3,PetscScalar f3,
236  |  PetscScalar isoquant, PetscScalar edge0,PetscScalar edge1,PetscScalar edge3,
237  |  int whichplane, PetscScalar *color, struct isurface *thesurf)
238  | {
239  |   PetscScalar edge2,edge4,edge5;
240  | 
241  |   switch (whichplane) {
242  |   case 7: { /* Triangles on edges 0,1,3; 3,1,5 */
243  |     edge5 = (isoquant - f2) / (f3 - f2);
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(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
250  | 		color, thesurf);
251  |       storetri (COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
252  | 		COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
253  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
254  | 		color, thesurf);
255  |     }
256  |     else {
257  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
258  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
259  | 		COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
260  | 		color, thesurf);
261  |       storetri (COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
262  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
263  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
264  | 		color, thesurf);
265  |     }
266  |     break;
267  |   }
268  | 
269  |   case 6: { /* Triangles on edges 1,2,4; 4,2,3 */
270  |     edge2 = (isoquant - f2) / (f0 - f2);
271  |     edge4 = (isoquant - f1) / (f3 - f1);
272  | 
273  |     /* Use edge 1 direction to guarantee right-handedness */
274  |     if (f2 > f1) {
275  |       storetri (COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
276  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
277  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
278  | 		color, thesurf);
279  |       storetri (COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
280  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
281  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
282  | 		color, thesurf);
283  |     }
284  |     else {
285  |       storetri (COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
286  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
287  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
288  | 		color, thesurf);
289  |       storetri (COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
290  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
291  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
292  | 		color, thesurf);
293  |     }
294  |     break;
295  |   }
296  | 
297  |   case 5: { /* Triangles on edges 0,2,3 */
298  |     edge2 = (isoquant - f2) / (f0 - f2);
299  | 
300  |     /* Use edge 0 direction to guarantee right-handedness */
301  |     if (f1 > f0)
302  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
303  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
304  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
305  | 		color, thesurf);
306  |     else
307  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
308  | 		COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
309  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
310  | 		color, thesurf);
311  |     break;
312  |   }
313  | 
314  |   case 4: { /* Triangles on edges 3,4,5 */
315  |     edge4 = (isoquant - f1) / (f3 - f1);
316  |     edge5 = (isoquant - f2) / (f3 - f2);
317  | 
318  |     /* Use edge 3 direction to guarantee right-handedness */
319  |     if (f3 > f0)
320  |       storetri (COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
321  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
322  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
323  | 		color, thesurf);
324  |     else
325  |       storetri (COORD(x0,x3,edge3), COORD(y0,y3,edge3), COORD(z0,z3,edge3),
326  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
327  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
328  | 		color, thesurf);
329  |     break;
330  |   }
331  | 
332  |   case 3: { /* Triangles on edges 0,1,4 */
333  |     edge4 = (isoquant - f1) / (f3 - f1);
334  | 
335  |     /* Use edge 0 direction to guarantee right-handedness */
336  |     if (f1 > f0)
337  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
338  | 		COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
339  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
340  | 		color, thesurf);
341  |     else
342  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
343  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
344  | 		COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
345  | 		color, thesurf);
346  |     break;
347  |   }
348  | 
349  |   case 2: { /* Triangles on edges 1,2,5 */
350  |     edge2 = (isoquant - f2) / (f0 - f2);
351  |     edge5 = (isoquant - f2) / (f3 - f2);
352  | 
353  |     /* Use edge 1 direction to guarantee right-handedness */
354  |     if (f2 > f1)
355  |       storetri (COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
356  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
357  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
358  | 		color, thesurf);
359  |     else
360  |       storetri (COORD(x1,x2,edge1), COORD(y1,y2,edge1), COORD(z1,z2,edge1),
361  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
362  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
363  | 		color, thesurf);
364  |     break;
365  |   }
366  | 
367  |   case 1: { /* Triangles on edges 0,2,4; 4,2,5 */
368  |     edge2 = (isoquant - f2) / (f0 - f2);
369  |     edge4 = (isoquant - f1) / (f3 - f1);
370  |     edge5 = (isoquant - f2) / (f3 - f2);
371  | 
372  |     /* Use edge 0 direction to guarantee right-handedness */
373  |     if (f1 > f0) {
374  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
375  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
376  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
377  | 		color, thesurf);
378  |       storetri (COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
379  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
380  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
381  | 		color, thesurf);
382  |     }
383  |     else {
384  |       storetri (COORD(x0,x1,edge0), COORD(y0,y1,edge0), COORD(z0,z1,edge0),
385  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
386  | 		COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
387  | 		color, thesurf);
388  |       storetri (COORD(x2,x0,edge2), COORD(y2,y0,edge2), COORD(z2,z0,edge2),
389  | 		COORD(x1,x3,edge4), COORD(y1,y3,edge4), COORD(z1,z3,edge4),
390  | 		COORD(x2,x3,edge5), COORD(y2,y3,edge5), COORD(z2,z3,edge5),
391  | 		color, thesurf);
392  |     }
393  |     break;
394  |   }
395  |   }
396  |   return 0;
397  | }
398  | #undef COORD
399  | 
400  | 
401  | #undef __FUNCT__
402  | #define __FUNCT__ "DrawTet"
403  | 
404  | /*++++++++++++++++++++++++++++++++++++++
405  |   This sets the edge and whichplane parameters and then passes everything to
406  |   DrawTetWithPlane to actually draw the triangle.  It is intended for use by
407  |   developers with distributed arrays based on tetrahedra, e.g. a finite element
408  |   mesh.
409  | 
410  |   int DrawTet Returns 0 or an error code.
411  | 
412  |   ISurface Surf ISurface object into which to draw this tetrahedron's isoquant.
413  | 
414  |   PetscScalar *coords Coordinates of tetrahedron corner points: x0, y0, z0, x1,
415  |   etc.
416  | 
417  |   PetscScalar *vals Function values at tetrahedron corners: f0, f1, f2, f3.
418  | 
419  |   PetscScalar isoquant Function value at which to draw triangle.
420  | 
421  |   PetscScalar *color R,G,B,A quad for this tetrahedron.
422  |   ++++++++++++++++++++++++++++++++++++++*/
423  | 
424  | int DrawTet
425  | (ISurface Surf, PetscScalar *coords, PetscScalar *vals, PetscScalar isoquant,
426  |  PetscScalar *color)
427  | {
428  |   PetscScalar edge0, edge1, edge3;
429  |   int whichplane, ierr=0;
430  |   struct isurface *thesurf = (struct isurface *) Surf;
431  | 
432  |   edge0 = (isoquant-vals[0]) / (vals[1]-vals[0]);
433  |   edge1 = (isoquant-vals[1]) / (vals[2]-vals[1]);
434  |   edge3 = (isoquant-vals[0]) / (vals[3]-vals[0]);
435  |   whichplane = (edge0>0. && edge0<1.) | ((edge1>0. && edge1<1.) << 1) |
436  |     ((edge3>0. && edge3<1.) << 2);
437  |   if (whichplane)
438  |     ierr = DrawTetWithPlane (coords[0],coords[1],coords[2],vals[0],
439  | 			     coords[3],coords[4],coords[5],vals[1],
440  | 			     coords[6],coords[7],coords[8],vals[2],
441  | 			     coords[9],coords[10],coords[11],vals[3],
442  | 			     isoquant, edge0,edge1,edge3, whichplane, color,
443  | 			     thesurf);
444  |   return ierr;
445  | }
446  | 
447  | 
448  | #undef __FUNCT__
449  | #define __FUNCT__ "DrawHex"
450  | 
451  | /*++++++++++++++++++++++++++++++++++++++
452  |   This divides a right regular hexahedron into tetrahedra, and loops over them
453  |   to generate triangles on each one.  It calculates edge and whichplane
454  |   parameters so it can use DrawTetWithPlane directly.
455  | 
456  |   int DrawHex Returns 0 or an error code.
457  | 
458  |   ISurface Surf ISurface object into which to draw this hexahedron's isoquant.
459  | 
460  |   PetscScalar *coords Coordinates of hexahedron corner points: xmin, xmax,
461  |   ymin, etc.
462  | 
463  |   PetscScalar *vals Function values at hexahedron corners: f0, f1, f2, etc.
464  | 
465  |   PetscScalar isoquant Function value at which to draw triangles.
466  | 
467  |   PetscScalar *color R,G,B,A quad for this hexahedron.
468  |   ++++++++++++++++++++++++++++++++++++++*/
469  | 
470  | int DrawHex
471  | (ISurface Surf, PetscScalar *coords, PetscScalar *vals, PetscScalar isoquant,
472  |  PetscScalar *color)
473  | {
474  |   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}};
475  |   int tet,node,ierr;
476  |   int c0,c1,c2,c3,whichplane;
477  |   PetscScalar edge0,edge1,edge3;
478  |   struct isurface *thesurf = (struct isurface *) Surf;
479  | 
480  |   for(tet=0; tet<6; tet++)
481  |     {
482  |       /* Within a tetrahedron, edges 0 through 5 connect corners:
483  | 	 0,1; 1,2; 2,0; 0,3; 1,3; 2,3 respectively */
484  |       c0 = tetras[tet][0];
485  |       c1 = tetras[tet][1];
486  |       c2 = tetras[tet][2];
487  |       c3 = tetras[tet][3];
488  |       edge0 = (isoquant-vals[c0]) / (vals[c1]-vals[c0]);
489  |       edge1 = (isoquant-vals[c1]) / (vals[c2]-vals[c1]);
490  |       edge3 = (isoquant-vals[c0]) / (vals[c3]-vals[c0]);
491  |       whichplane = (edge0>0. && edge0<1.) | ((edge1>0. && edge1<1.) << 1) |
492  | 	((edge3>0. && edge3<1.) << 2);
493  |       if (whichplane)
494  | 	{
495  | 	  ierr=DrawTetWithPlane
496  | 	    (coords[c0&1],coords[2+((c0&2)>>1)],coords[4+((c0&4)>>2)],vals[c0],
497  | 	     coords[c1&1],coords[2+((c1&2)>>1)],coords[4+((c1&4)>>2)],vals[c1],
498  | 	     coords[c2&1],coords[2+((c2&2)>>1)],coords[4+((c2&4)>>2)],vals[c2],
499  | 	     coords[c3&1],coords[2+((c3&2)>>1)],coords[4+((c3&4)>>2)],vals[c3],
500  | 	     isoquant, edge0,edge1,edge3, whichplane, color, thesurf);
501  | 	  CHKERRQ (ierr);
502  | 	}
503  |     }
504  | 
505  |   return 0;
506  | }
507  | 
508  | 
509  | #undef __FUNCT__
510  | #define __FUNCT__ "Draw3DBlock"
511  | 
512  | /*++++++++++++++++++++++++++++++++++++++
513  |   Calculate vertices of isoquant triangle(s) in a 3-D regular array of right
514  |   regular hexahedra.  This loops through a 3-D array and calls DrawHex to
515  |   calculate the triangulation of each hexahedral cell.
516  | 
517  |   int Draw3DBlock Returns 0 or an error code.
518  | 
519  |   ISurface Surf ISurface object into which to draw this block's isoquants.
520  | 
521  |   int xd Overall x-width of function value array.
522  | 
523  |   int yd Overall y-width of function value array.
524  | 
525  |   int zd Overall z-width of function value array.
526  | 
527  |   int xs X-index of the start of the array section we'd like to draw.
528  | 
529  |   int ys Y-index of the start of the array section we'd like to draw.
530  | 
531  |   int zs Z-index of the start of the array section we'd like to draw.
532  | 
533  |   int xm X-width of the array section we'd like to draw.
534  | 
535  |   int ym Y-width of the array section we'd like to draw.
536  | 
537  |   int zm Z-width of the array section we'd like to draw.
538  | 
539  |   PetscScalar *minmax Position of block corners: xmin, xmax, ymin, ymax, zmin,
540  |   zmax.
541  | 
542  |   PetscScalar *vals The array of function values at vertices.
543  | 
544  |   int skip Number of interlaced fields in this array.
545  | 
546  |   int n_quants Number of isoquant surfaces to draw (isoquant values).
547  | 
548  |   PetscScalar *isoquants Array of function values at which to draw triangles.
549  | 
550  |   PetscScalar *colors Array of color R,G,B,A quads for each isoquant.
551  |   ++++++++++++++++++++++++++++++++++++++*/
552  | 
553  | int Draw3DBlock
554  | (ISurface Surf,int xd,int yd,int zd, int xs,int ys,int zs,int xm,int ym,int zm,
555  |  PetscScalar *minmax, PetscScalar *vals, int skip, int n_quants,
556  |  PetscScalar *isoquants, PetscScalar *colors)
557  | {
558  |   int i,j,k,q, start, ierr;
559  |   PetscScalar boxcoords[6], funcs[8];
560  | 
561  |   /* Simple check */
562  |   if (xd<=xs+xm || yd<=ys+ym || zd<=zs+zm)
563  |     SETERRQ (PETSC_ERR_ARG_WRONGSTATE, "Array size mismatch");
564  | 
565  |   /* The big loop */
566  |   for (k=0; k<zm; k++)
567  |     {
568  |       boxcoords[4] = minmax[4] + (minmax[5]-minmax[4])*k/zm;
569  |       boxcoords[5] = minmax[4] + (minmax[5]-minmax[4])*(k+1)/zm;
570  |       for(j=0;j<ym;j++)
571  | 	{
572  | 	  boxcoords[2] = minmax[2] + (minmax[3]-minmax[2])*j/ym;
573  | 	  boxcoords[3] = minmax[2] + (minmax[3]-minmax[2])*(j+1)/ym;
574  | 	  for(i=0;i<xm;i++)
575  | 	    {
576  | 	      boxcoords[0] = minmax[0] + (minmax[1]-minmax[0])*i/xm;
577  | 	      boxcoords[1] = minmax[0] + (minmax[1]-minmax[0])*(i+1)/xm;
578  | 	      start = ((k+zs)*yd + j+ys)*xd + xs + i;
579  | 	      funcs[0] = vals[skip*start];
580  | 	      funcs[1] = vals[skip*(start+1)];
581  | 	      funcs[2] = vals[skip*(start+xd)];
582  | 	      funcs[3] = vals[skip*(start+xd+1)];
583  | 	      funcs[4] = vals[skip*(start+xd*yd)];
584  | 	      funcs[5] = vals[skip*(start+xd*yd+1)];
585  | 	      funcs[6] = vals[skip*(start+xd*yd+xd)];
586  | 	      funcs[7] = vals[skip*(start+xd*yd+xd+1)];
587  | 	      for (q=0; q<n_quants; q++)
588  | 		{
589  | 		  ierr = DrawHex (Surf, boxcoords, funcs, isoquants[q],
590  | 				  colors+4*q); CHKERRQ (ierr);
591  | 		}
592  | 	    }
593  | 	}
594  |     }
595  | 
596  |   return 0;
597  | }