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 | }