1 | /***************************************
2 | $Header: /cvsroot/petscgraphics/illuminator.c,v 1.9 2004/08/17 15:05:56 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 | }