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