1 /*****************************************************************************
2 *
3 * Elmer, A Finite Element Software for Multiphysical Problems
4 *
5 * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library (in file ../LGPL-2.1); if not, write
19 * to the Free Software Foundation, Inc., 51 Franklin Street,
20 * Fifth Floor, Boston, MA 02110-1301 USA
21 *
22 *****************************************************************************/
23
24 /*******************************************************************************
25 *
26 * MATC surface display routines.
27 *
28 *******************************************************************************
29 *
30 * Author: Juha Ruokolainen
31 *
32 * Address: CSC - IT Center for Science Ltd.
33 * Keilaranta 14, P.O. BOX 405
34 * 02101 Espoo, Finland
35 * Tel. +358 0 457 2723
36 * Telefax: +358 0 457 2302
37 * EMail: Juha.Ruokolainen@csc.fi
38 *
39 * Date: 30 May 1996
40 *
41 * Modified by:
42 *
43 * Date of modification:
44 *
45 ******************************************************************************/
46 /***********************************************************************
47 |
48 | Written ??-???-88 / JPR
49 | Last edited 16-FEB-89 / OP
50 |
51 * File: c3d
52 ^
53 | USING THE C3D (preliminary)
54 |
55 | Full usage will be included later!
56 |
57 ^**********************************************************************/
58
59
60 /*
61 * $Id: c3d.c,v 1.2 2005/05/27 12:26:14 vierinen Exp $
62 *
63 * $Log: c3d.c,v $
64 * Revision 1.2 2005/05/27 12:26:14 vierinen
65 * changed header install location
66 *
67 * Revision 1.1.1.1 2005/04/14 13:29:14 vierinen
68 * initial matc automake package
69 *
70 * Revision 1.2 1998/08/01 12:34:30 jpr
71 *
72 * Added Id, started Log.
73 *
74 *
75 */
76
77 #include "elmer/matc.h"
78
79 #define C3D_MASK 9
80 #define C3D_HALFMASK (1 << (C3D_MASK - 1))
81
82 /***********************************************************************
83 Typedefinitions for panel tree
84 ***********************************************************************/
85 struct node
86 {
87 int x, y, z, d;
88 };
89
90 struct element
91 {
92 struct node *node[4];
93 int d, z;
94 };
95
96 struct el_tree
97 {
98 struct el_tree *left, *right;
99 struct element *entry;
100 };
101
102 void C3D_Contour(double *, int, int);
103 void C3D_Levels(int);
104
105 void C3D_Show_Tri(int *, int *, int *d);
106 void C3D_Show_Elem(struct element *el);
107 void C3D_SelCol(int);
108
109 void C3D_Show_El_Tree(struct el_tree *head);
110 void C3D_Add_El_Tree(struct el_tree *head, struct el_tree *add);
111
112 void C3D_Pcalc(int, int, int, int, int, int, int *, int *, int *,int *);
113 void C3D_AreaFill(int, int, int *, int *);
114
c3d_gc3d(var)115 VARIABLE *c3d_gc3d(var) VARIABLE *var;
116 {
117 C3D_Contour(MATR(var), NROW(var), NCOL(var));
118 return (VARIABLE *)NULL;
119 }
120
c3d_gc3dlevels(var)121 VARIABLE *c3d_gc3dlevels(var) VARIABLE *var;
122 {
123 C3D_Levels((int)*MATR(var));
124 return (VARIABLE *)NULL;
125 }
126
127 /***********************************************************************
128 Globals needed for C3D All start with the prefix c3d_
129 ***********************************************************************/
130 static int c3d_clevels = 10,
131 c3d_perspective = FALSE;
132 #pragma omp threadprivate (c3d_clevels, c3d_perspective)
133
134 /***********************************************************************
135 void C3D_Contour(matrix, nr, nc) double *matrix; int nr, nc;
136 ? Main function to be used.
137 | This function displays a matrix seen from the current C3D-viewpoint
138 | and with number of levels selected.
139 |
140 | NOTICE! CALL FROM FORTRAN:
141 |
142 | CALL C3D_Contour( a , %val( DIMY ) , %val( DIMX ) )
143 |
144 | The a should be of type double precision (or REAL*8 in VAX/VMS)
145 | and the dimensions are in reverse order
146 ! No error messages are generated
147 & C3D...
148 ***********************************************************************/
C3D_Contour(double * matrix,int nr,int nc)149 void C3D_Contour(double *matrix, int nr, int nc)
150 {
151 double xmin, xmax, ymin, ymax, dmin, dmax;
152 double x, y, z, d, xs, ys, zs;
153
154 struct el_tree *tr, *trb, *element_el_tree = NULL;
155 struct element *el, *elements = NULL;
156 struct node *nodes = NULL;
157
158 GMATRIX gm;
159
160 static GMATRIX ident =
161 {
162 1, 0, 0, 0,
163 0, 1, 0, 0,
164 0, 0, 1, 0,
165 0, 0, 0, 1
166 };
167 #pragma omp threadprivate (ident)
168
169 int i, j, k, n;
170
171 nodes = (struct node *)ALLOCMEM(nr * nc * sizeof(struct node));
172
173 xmin = ymin = dmin = 1.0e20;
174 xmax = ymax = dmax = -1.0e20;
175
176 for(i = 0, n = 0; i < nr; i++)
177 for(j = 0; j < nc; j++, n++) {
178 z = matrix[n]; dmin = min(dmin, z); dmax = max(dmax, z);
179 }
180
181 for(i = 0, n = 0; i < nr; i++)
182 for(j = 0; j < nc; j++, n++)
183 {
184 x = 2 * (double)i / (double)nr - 1;
185 y = 2 * (double)j / (double)nc - 1;
186 d = (matrix[n] - dmin) / (dmax - dmin);
187 z = 2 * d - 1;
188 gra_mtrans(x, y, z, &xs, &ys, &zs);
189 xs *= (1 << 20);
190 ys *= (1 << 20);
191 zs *= (1 << 20);
192 nodes[n].x = xs;
193 nodes[n].y = ys;
194 nodes[n].z = zs;
195 nodes[n].d = (c3d_clevels * d + 1) * (1 << C3D_MASK);
196 xmin = min(xmin, xs);
197 xmax = max(xmax, xs);
198 ymin = min(ymin, ys);
199 ymax = max(ymax, ys);
200 }
201
202 for(i = 0, n = 0; i < nr; i++)
203 for(j = 0; j < nc; j++, n++)
204 {
205 nodes[n].x = 4095 * (nodes[n].x-xmin) / (xmax-xmin);
206 nodes[n].y = 4095 * (nodes[n].y-ymin) / (ymax-ymin);
207 }
208
209 elements = (struct element *)
210 ALLOCMEM((nr - 1) * (nc - 1) * sizeof(struct element));
211
212 trb = (struct el_tree *)
213 ALLOCMEM((nr -1) * (nc - 1) * sizeof(struct el_tree));
214
215 for(i = 0, n = 0; i < nr - 1; i++)
216 for(j = 0; j < nc - 1; j++, n++) {
217 tr = &trb[n];
218 el = tr -> entry = &elements[n];
219 el -> node[0] = &nodes[nc*i + j];
220 el -> node[1] = &nodes[nc*(i+1) + j];
221 el -> node[2] = &nodes[nc*(i+1) + j + 1];
222 el -> node[3] = &nodes[nc*i + j + 1];
223 el -> d = 0; el -> z = 0;
224 for(k = 0; k < 4; k++) {
225 el -> d += el -> node[k] -> d;
226 el -> z += el -> node[k] -> z;
227 }
228 el -> d = (el -> d + 2) >> 2;
229 tr -> left = tr -> right = NULL;
230 if (element_el_tree == NULL)
231 element_el_tree = tr;
232 else
233 C3D_Add_El_Tree(element_el_tree, tr);
234 }
235
236 GRA_GETMATRIX(gm);
237 GRA_SETMATRIX(ident);
238 GRA_WINDOW(0.0,4096.0,0.0,4096.0,-1.0,1.0);
239
240 C3D_Show_El_Tree(element_el_tree);
241
242 if (gra_state.pratio>0)
243 GRA_PERSPECTIVE(gra_state.pratio);
244 GRA_SETMATRIX(gm);
245 GRA_FLUSH();
246
247 FREEMEM(elements);
248 FREEMEM(trb);
249 FREEMEM(nodes);
250 }
251
C3D_Levels(levels)252 void C3D_Levels(levels) int levels;
253 /***********************************************************************
254 | Set the number of colorlevels to be used in coloring
255 | use 0 for singlecolor hidden surface plot
256 ! Number of levels must be between 0 - 13
257 ^**********************************************************************/
258 {
259 if (levels >= 0)
260 c3d_clevels = levels;
261 else
262 error("C3D_Levels: level parameter negative, not changed.\n");
263 }
264
C3D_Persp(tf)265 void C3D_Persp(tf) int tf;
266 /***********************************************************************
267 | Use perspective or orthogonal transformation ?
268 ! Number of levels must be between 0 - 13
269 ^**********************************************************************/
270 {
271 c3d_perspective = tf;
272 }
273
C3D_Add_El_Tree(head,add)274 void C3D_Add_El_Tree(head, add) struct el_tree *head, *add;
275 /***********************************************************************
276 | Add a single element into the current leaf of the element tree
277 | Only internal use
278 ^**********************************************************************/
279 {
280 if (add -> entry -> z > head -> entry -> z) {
281 if (head -> left == NULL) {
282 head -> left = add;
283 }
284 else {
285 C3D_Add_El_Tree(head -> left, add);
286 }
287 }
288 else if (add -> entry -> z < head -> entry -> z) {
289 if (head -> right == NULL) {
290 head -> right = add;
291 }
292 else {
293 C3D_Add_El_Tree(head -> right, add);
294 }
295 }
296 else {
297 add -> left = head -> left;
298 head -> left = add;
299 }
300 }
301
C3D_Show_El_Tree(head)302 void C3D_Show_El_Tree(head) struct el_tree *head;
303 /***********************************************************************
304 | Display the contents of the current leaf of the element tree
305 | Only internal use
306 ^**********************************************************************/
307 {
308 if (head == NULL) return;
309 C3D_Show_El_Tree(head->left);
310 C3D_Show_Elem(head->entry);
311 C3D_Show_El_Tree(head->right);
312 }
313
C3D_Free_El_Tree(head)314 void C3D_Free_El_Tree(head) struct el_tree *head;
315 /***********************************************************************
316 | Free the memory allocated to the current leaf of the element tree
317 | Only internal use
318 ^**********************************************************************/
319 {
320 if (head == NULL) return;
321 C3D_Free_El_Tree(head->left);
322 C3D_Free_El_Tree(head->right);
323 FREEMEM(head);
324 }
325
C3D_Convex_Test(x,y)326 int C3D_Convex_Test(x, y) int x[], y[];
327 /***********************************************************************
328 | Test four-node ploygon
329 | Only internal use
330 ^**********************************************************************/
331 {
332 int a1, a2, at1, at2, amax, aind;
333
334 amax = abs(y[0]*(x[2]-x[1])+y[1]*(x[0]-x[2])+y[2]*(x[1]-x[0]));
335 aind = 3;
336
337 at1 = abs(y[2]*(x[0]-x[3])+y[3]*(x[2]-x[0])+y[0]*(x[3]-x[2]));
338
339 a1 = amax + at1;
340
341 if (at1 > amax) {
342 amax = at1; aind = 1;
343 }
344
345 at1 = abs(y[1]*(x[3]-x[2])+y[2]*(x[1]-x[3])+y[3]*(x[2]-x[1]));
346
347 if (at1 > amax) {
348 amax = at1; aind = 0;
349 }
350
351 at2 = abs(y[3]*(x[1]-x[0])+y[0]*(x[3]-x[1])+y[1]*(x[0]-x[3]));
352
353 if (at2 > amax) {
354 aind = 2;
355 }
356
357 a2 = at1 + at2;
358
359 if (a1 == a2) return -1;
360
361 return aind;
362 }
363
C3D_Show_Elem(el)364 void C3D_Show_Elem(el) struct element *el;
365 /***********************************************************************
366 | Display a single element by coloring and transformation
367 | Only internal use
368 ^**********************************************************************/
369 {
370 int x[5], y[5], z[5], zz, xp, yp;
371 int xi[5], yi[5], i;
372 int d[5], zp, col[3];
373
374 Point p[5];
375
376 for(i = 0; i != 4; i++)
377 {
378 x[i] = el -> node[i] -> x;
379 y[i] = el -> node[i] -> y;
380 d[i] = el -> node[i] -> d;
381 }
382
383 if ((d[0]>>C3D_MASK) == (d[1]>>C3D_MASK) &&
384 (d[0]>>C3D_MASK) == (d[2]>>C3D_MASK) &&
385 (d[0]>>C3D_MASK) == (d[3]>>C3D_MASK))
386 {
387 C3D_SelCol(d[0]>>C3D_MASK);
388 C3D_AreaFill(4, TRUE, x, y);
389 return;
390 }
391
392 switch(C3D_Convex_Test(x,y))
393 {
394 case 0:
395 C3D_Show_Tri(x, y, d);
396 xi[0] = x[2]; xi[1] = x[3]; xi[2] = x[0];
397 yi[0] = y[2]; yi[1] = y[3]; yi[2] = y[0];
398 col[0] = d[2]; col[1] = d[3]; col[2] = d[0];
399 C3D_Show_Tri(xi, yi, col);
400 break;
401
402 case 1:
403 C3D_Show_Tri(&x[1], &y[1], &d[1]);
404 xi[0] = x[0]; xi[1] = x[1]; xi[2] = x[3];
405 yi[0] = y[0]; yi[1] = y[1]; yi[2] = y[3];
406 col[0] = d[0]; col[1] = d[1]; col[2] = d[3];
407 C3D_Show_Tri(xi, yi, col);
408 break;
409
410 case 2:
411 C3D_Show_Tri(x, y, d);
412 xi[0] = x[2]; xi[1] = x[3]; xi[2] = x[0];
413 yi[0] = y[2]; yi[1] = y[3]; yi[2] = y[0];
414 col[0] = d[2]; col[1] = d[3]; col[2] = d[0];
415 C3D_Show_Tri(xi, yi, col);
416 break;
417
418 case 3:
419 C3D_Show_Tri(&x[1], &y[1], &d[1]);
420 xi[0] = x[0]; xi[1] = x[1]; xi[2] = x[3];
421 yi[0] = y[0]; yi[1] = y[1]; yi[2] = y[3];
422 col[0] = d[0]; col[1] = d[1]; col[2] = d[3];
423 C3D_Show_Tri(xi, yi, col);
424 break;
425
426 default:
427 xp = 0; yp = 0;
428 for(i = 0; i != 4; i++)
429 {
430 xp += x[i]; yp += y[i];
431 }
432 xp = (xp + 2) >> 2; yp = (yp + 2) >> 2;
433 zp = el -> d;
434
435 xi[0] = x[0]; xi[1] = x[1]; xi[2] = xp;
436 yi[0] = y[0]; yi[1] = y[1]; yi[2] = yp;
437 col[0] = d[0]; col[1] = d[1]; col[2] = zp;
438 C3D_Show_Tri(xi, yi, col);
439
440 xi[0] = x[1]; xi[1] = x[2];
441 yi[0] = y[1]; yi[1] = y[2];
442 col[0] = d[1]; col[1] = d[2];
443 C3D_Show_Tri(xi, yi, col);
444
445 xi[0] = x[2]; xi[1] = x[3];
446 yi[0] = y[2]; yi[1] = y[3];
447 col[0] = d[2]; col[1] = d[3];
448 C3D_Show_Tri(xi, yi, col);
449
450 xi[0] = x[3]; xi[1] = x[0];
451 yi[0] = y[3]; yi[1] = y[0];
452 col[0] = d[3]; col[1] = d[0];
453 C3D_Show_Tri(xi, yi, col);
454 break;
455 }
456
457 p[0].x = (int)(x[0]+0.5); p[0].y = (int)(y[0]+0.5); p[0].z = 0.0;
458 p[1].x = (int)(x[1]+0.5); p[1].y = (int)(y[1]+0.5); p[1].z = 0.0;
459 p[2].x = (int)(x[2]+0.5); p[2].y = (int)(y[2]+0.5); p[2].z = 0.0;
460 p[3].x = (int)(x[3]+0.5); p[3].y = (int)(y[3]+0.5); p[3].z = 0.0;
461 p[4].x = (int)(x[0]+0.5); p[4].y = (int)(y[0]+0.5); p[4].z = 0.0;
462 GRA_COLOR(1);
463 GRA_POLYLINE(5, p);
464 }
465
C3D_Show_Tri(x,y,d)466 void C3D_Show_Tri(x, y, d) int x[3], y[3], d[3];
467 /***********************************************************************
468 | Only internal use
469 ^**********************************************************************/
470 {
471 int xx[128], yy[128], dd[128], px[7], py[7];
472 int i, j, k, n = 0;
473
474 if ((d[0] >> C3D_MASK) == (d[1] >> C3D_MASK) &&
475 (d[0] >> C3D_MASK) == (d[2] >> C3D_MASK))
476 {
477 C3D_SelCol(d[0] >> C3D_MASK);
478 C3D_AreaFill(3, FALSE, x, y);
479 return;
480 }
481
482 C3D_Pcalc(x[0],y[0],d[0],x[1],y[1],d[1],&i,&xx[n],&yy[n],&dd[n]); n += i;
483 C3D_Pcalc(x[1],y[1],d[1],x[2],y[2],d[2],&i,&xx[n],&yy[n],&dd[n]); n += i;
484 C3D_Pcalc(x[2],y[2],d[2],x[0],y[0],d[0],&i,&xx[n],&yy[n],&dd[n]); n += i;
485
486 for(i = 0; i < 2; i++)
487 {
488 xx[n+i] = xx[i];
489 yy[n+i] = yy[i];
490 dd[n+i] = dd[i];
491 }
492
493 for(i = 0, k = 0; i < n - 2; k = 0, i++)
494 {
495 px[k] = xx[i]; py[k++] = yy[i];
496 px[k] = xx[i+1]; py[k++] = yy[i+1];
497
498 if (dd[i] == dd[i+1]) {
499 i++; px[k] = xx[i+1]; py[k++] = yy[i+1];
500 }
501
502 for(j = n - 1; j > i; j--) {
503 if (dd[i] == dd[j]) {
504 if (dd[j-1] == dd[j]) {
505 px[k] = xx[j-1]; py[k++] = yy[j-1];
506 }
507 px[k] = xx[j]; py[k++] = yy[j];
508 px[k] = xx[j+1]; py[k++] = yy[j+1];
509 if (dd[j] == dd[j+1]) {
510 j++; px[k] = xx[j+1]; py[k++] = yy[j+1];
511 }
512 break;
513 }
514 }
515
516 if (k > 2) {
517 C3D_SelCol(dd[i]);
518 C3D_AreaFill(k, FALSE, px, py);
519 }
520
521 }
522 }
523
C3D_Pcalc(x0,y0,d0,x1,y1,d1,n,xx,yy,dd)524 void C3D_Pcalc(x0, y0, d0, x1, y1, d1, n, xx, yy, dd)
525 /**/ int x0, y0, d0, x1, y1, d1, *n, xx[], yy[], dd[];
526 /***********************************************************************
527 | Only internal use
528 ^**********************************************************************/
529 {
530 int x, y, deltax, deltay, deltad;
531 int i, d, ds;
532
533 *n = abs((d1 >> C3D_MASK) - (d0 >> C3D_MASK));
534
535 xx[0] = x0; yy[0] = y0; dd[0] = d0 >> C3D_MASK;
536 (*n)++;
537
538 if (*n != 1) {
539
540 deltad = (d1 >= d0 ? 1 : (-1));
541 ds = d0 & ((1 << C3D_MASK)-1);
542
543 if (d1 > d0) {
544 ds = (1 << C3D_MASK) - ds;
545 }
546
547 d = abs(d1 - d0);
548
549 if (x1 > x0) {
550 deltax = ((x1 - x0) << (C3D_MASK << 1)) / d;
551 deltax >>= C3D_MASK;
552 x = x0 + ((ds * deltax + C3D_HALFMASK) >> C3D_MASK);
553 }
554 else {
555 deltax = ((x0 - x1) << (C3D_MASK << 1)) / d;
556 deltax >>= C3D_MASK;
557 x = x0 - ((ds * deltax + C3D_HALFMASK) >> C3D_MASK);
558 deltax = -deltax;
559 }
560
561 if (y1 > y0) {
562 deltay = ((y1 - y0) << (C3D_MASK << 1)) / d;
563 deltay >>= C3D_MASK;
564 y = y0 + ((ds * deltay + C3D_HALFMASK) >> C3D_MASK);
565 }
566 else {
567 deltay = ((y0 - y1) << (C3D_MASK << 1)) / d;
568 deltay >>= C3D_MASK;
569 y = y0 - ((ds * deltay + C3D_HALFMASK) >> C3D_MASK);
570 deltay = -deltay;
571 }
572
573 for(i = 1; i != *n; i++) {
574 dd[i] = dd[i-1] + deltad;
575 xx[i] = x; yy[i] = y;
576 x = x + deltax; y = y + deltay;
577 }
578 }
579 }
580
C3D_SelCol(col)581 void C3D_SelCol(col) int col;
582 /***********************************************************************
583 | Only internal use
584 ^**********************************************************************/
585 {
586 GRA_COLOR(++col);
587 }
588
C3D_AreaFill(n,border,x,y)589 void C3D_AreaFill(n, border, x, y) int x[], y[], n, border;
590 /***********************************************************************
591 | Only internal use
592 ^**********************************************************************/
593 {
594 int i, j;
595 Point p[8];
596
597 while(n > 0 && x[n-1] == x[0] && y[n-1] == y[0]) n--;
598
599 for(i = 0; i < n; i++) {
600 p[i].x = (int)(x[i]+0.5); p[i].y = (int)(y[i]+0.5); p[i].z = 0.0;
601 }
602
603 for(i = 0; i < n - 1; i++)
604 if (p[i].x == p[i+1].x && p[i].y == p[i+1].y) {
605 for(j = i + 1; j < n - 1; j++) {
606 p[j].x = p[j+1].x; p[j].y = p[j+1].y;
607 }
608 n--;
609 }
610
611 if (n > 2) {
612 GRA_AREAFILL(n, p);
613 if (border)
614 {
615 p[n].x = p[0].x;
616 p[n].y = p[0].y;
617 p[n].z = 0;
618 GRA_COLOR(1);
619 GRA_POLYLINE(n+1, p);
620 }
621 }
622 }
623