1 /*---------------------------------------------------------------------------*/
2 /* ALBERTA: an Adaptive multi Level finite element toolbox using */
3 /* Bisectioning refinement and Error control by Residual */
4 /* Techniques for scientific Applications */
5 /* */
6 /* www.alberta-fem.de */
7 /* */
8 /*---------------------------------------------------------------------------*/
9 /* */
10 /* file: macro.c */
11 /* */
12 /* description: reading of macro triangulations; file includes */
13 /* 1d/macro_1d.c, 2d/macro_2d.c, 3d/macro_3d.c which contain */
14 /* the dimension dependent parts. */
15 /* */
16 /*---------------------------------------------------------------------------*/
17 /* */
18 /* authors: Alfred Schmidt */
19 /* Zentrum fuer Technomathematik */
20 /* Fachbereich 3 Mathematik/Informatik */
21 /* Universitaet Bremen */
22 /* Bibliothekstr. 2 */
23 /* D-28359 Bremen, Germany */
24 /* */
25 /* Kunibert G. Siebert */
26 /* Institut fuer Mathematik */
27 /* Universitaet Augsburg */
28 /* Universitaetsstr. 14 */
29 /* D-86159 Augsburg, Germany */
30 /* */
31 /* Daniel Koester */
32 /* Institut fuer Mathematik */
33 /* Albert-Ludwigs-Universitaet Freiburg */
34 /* Hermann-Herder-Str. 10 */
35 /* D-79104 Freiburg */
36 /* */
37 /* Claus-Justus Heine */
38 /* Abteilung fuer Angewandte Mathematik */
39 /* Universitaet Freiburg */
40 /* Hermann-Herder-Strasse 10 */
41 /* 79104 Freiburg, Germany */
42 /* */
43 /* http://www.alberta-fem.de */
44 /* */
45 /* (c) by A. Schmidt and K.G. Siebert (1996-2003), */
46 /* C.-J. Heine (2006), D. Koester (2004-2007). */
47 /* */
48 /*---------------------------------------------------------------------------*/
49
50 #include "alberta_intern.h"
51 #include "alberta.h"
52
53 #define VERT_IND(dim, i, j) ((i)*N_VERTICES(dim)+(j))
54 #define NEIGH_IND(dim, i, j) ((i)*N_NEIGH(dim)+(j))
55
56 /*--------------------------------------------------------------------------*/
57 /* opp_vertex() checks whether the vertex/edge/face with vertices */
58 /* test[0],.., test[dim-1] is part of mel's boundary. It returns the */
59 /* opposite vertex if true else -1. */
60 /*--------------------------------------------------------------------------*/
61
opp_vertex(int dim,int * mel_vert,int * test,S_CHAR * mapping)62 static S_CHAR opp_vertex(int dim, int *mel_vert, int *test, S_CHAR *mapping)
63 {
64 int i, j, nv = 0, ov = 0;
65
66 for (i = 0; i < N_VERTICES(dim); i++) {
67 if (nv < i-1) return(-1);
68
69 for (j = 0; j < dim; j++) {
70 if (mel_vert[i] == test[j]) {
71 /*--------------------------------------------------------------------------*/
72 /* i is a common vertex */
73 /*--------------------------------------------------------------------------*/
74 if (mapping)
75 mapping[j] = i;
76 ov += i;
77 nv++;
78 break;
79 }
80 }
81
82 }
83 if (nv != dim) return(-1);
84 /*--------------------------------------------------------------------------*/
85 /* the opposite vertex is DIM! - (sum of indices of common vertices) */
86 /*--------------------------------------------------------------------------*/
87 if (dim == 1)
88 return(1-ov);
89 else if(dim == 2)
90 return(3-ov);
91 else
92 return(6-ov);
93 }
94
95 /*--------------------------------------------------------------------------*/
96 /* compute_neigh_fast() is an algorithm meant to speed up the task of */
97 /* computing neighbours. It does not use an N^2-algorithm. */
98 /* The idea is to link vertices to elements sharing them to make the */
99 /* search for neighbours more efficient - at the cost of some additional */
100 /* temporary memory usage. Daniel Koester */
101 /*--------------------------------------------------------------------------*/
102
compute_neigh_fast(MACRO_DATA * data)103 void compute_neigh_fast(MACRO_DATA *data)
104 {
105 FUNCNAME("compute_neigh_fast");
106 int dim = data->dim;
107 int i, j, k, l, wt, index, vertices[DIM_MAX] = { 0, }, info=0;
108 int neigh_found = false;
109
110 struct vert2elem {
111 struct vert2elem *next;
112 int mel;
113 };
114
115 typedef struct vert2elem VERT2ELEM;
116
117 VERT2ELEM *buffer, *buffer2;
118 VERT2ELEM **list = MEM_CALLOC(data->n_total_vertices, VERT2ELEM *);
119
120 if (!data->neigh)
121 data->neigh = MEM_ALLOC(data->n_macro_elements*N_NEIGH(dim), int);
122
123 if (!data->opp_vertex)
124 data->opp_vertex = MEM_ALLOC(data->n_macro_elements*N_NEIGH(dim), int);
125
126 /*--------------------------------------------------------------------------*/
127 /* first initialize elements (-2 as "undefined") */
128 /*--------------------------------------------------------------------------*/
129
130 for (i = 0; i < data->n_macro_elements; i++)
131 for (j = 0; j < N_NEIGH(dim); j++)
132 data->neigh[NEIGH_IND(dim, i, j)] = -2;
133
134 /*--------------------------------------------------------------------------*/
135 /* fill the array "list" of linked lists */
136 /*--------------------------------------------------------------------------*/
137
138 for(i = 0; i < data->n_macro_elements; i++) {
139 for(j = 0; j < N_VERTICES(dim); j++) {
140 buffer = list[(index=data->mel_vertices[VERT_IND(dim, i, j)])];
141 list[index] = MEM_ALLOC(1, VERT2ELEM);
142 list[index]->next = buffer;
143 list[index]->mel = i;
144 }
145 }
146
147 /*--------------------------------------------------------------------------*/
148 /* here comes the actual checking... */
149 /*--------------------------------------------------------------------------*/
150
151 for (i = 0; i < data->n_macro_elements; i++)
152 {
153 INFO(info, 4, "Current element %d\n", i);
154 INFO(info, 6, "with vertices: ");
155
156 for(j = 0; j < N_VERTICES(dim); j++)
157 PRINT_INFO(info, 6, "%d ", data->mel_vertices[VERT_IND(dim, i, j)]);
158 PRINT_INFO(info, 6, "\n");
159
160 for (j = 0; j < N_NEIGH(dim); j++)
161 {
162 if (data->neigh[NEIGH_IND(dim, i, j)] == -2)
163 {
164 INFO(info, 8, "looking for neighbour no %d\n", j);
165
166 /* If this face (j) has a wall-transformation attached to it,
167 * then we must actually check whether the transformed indices
168 * belong to the neighbour.
169 */
170 if (data->n_wall_vtx_trafos &&
171 (wt = data->el_wall_vtx_trafos[NEIGH_IND(dim, i, j)])) {
172 if (wt > 0) {
173 for (k = 0; k < N_VERTICES(dim-1); k++) {
174 index =
175 data->mel_vertices[VERT_IND(dim, i, (j+k+1)%N_VERTICES(dim))];
176 for (l = 0; l < N_VERTICES(dim-1); l++) {
177 if (data->wall_vtx_trafos[wt-1][l][0] == index) {
178 vertices[k] = data->wall_vtx_trafos[wt-1][l][1];
179 }
180 }
181 }
182 } else {
183 for (k = 0; k < N_VERTICES(dim-1); k++) {
184 index =
185 data->mel_vertices[VERT_IND(dim, i, (j+k+1)%N_VERTICES(dim))];
186 for (l = 0; l < N_VERTICES(dim-1); l++) {
187 if (data->wall_vtx_trafos[-wt-1][l][1] == index) {
188 vertices[k] = data->wall_vtx_trafos[-wt-1][l][0];
189 }
190 }
191 }
192 }
193 } else {
194 for (k = 0; k < N_VERTICES(dim-1); k++)
195 vertices[k] =
196 data->mel_vertices[VERT_IND(dim, i, (j+k+1)%N_VERTICES(dim))];
197 }
198
199 buffer = list[vertices[0]];
200
201 neigh_found = false;
202
203 while(buffer) {
204 if(buffer->mel != i) {
205 if ((l = opp_vertex(dim, data->mel_vertices+buffer->mel*N_VERTICES(dim), vertices, NULL)) != -1) {
206 TEST_EXIT(!neigh_found,
207 "Found two neighbours on wall %d of macro el %d!\n",
208 j, i);
209 data->neigh[NEIGH_IND(dim, i, j)] = buffer->mel;
210 data->neigh[NEIGH_IND(dim, buffer->mel, l)] = i;
211 data->opp_vertex[NEIGH_IND(dim, i, j)] = l;
212 data->opp_vertex[NEIGH_IND(dim, buffer->mel, l)] = j;
213 INFO(info, 8, "found element %d as neighbour...\n", buffer->mel);
214 neigh_found = true;
215 }
216 }
217
218 buffer = buffer->next;
219 }
220
221 if(!neigh_found) {
222 INFO(info, 8,
223 "no neighbour %d of element %d found: Assuming a boundary...\n",
224 j, i);
225
226 data->neigh[NEIGH_IND(dim, i, j)] = -1;
227 }
228 }
229 }
230 }
231
232 /*--------------------------------------------------------------------------*/
233 /* now is the time to clean up */
234 /*--------------------------------------------------------------------------*/
235
236
237 for(i = 0; i < data->n_total_vertices; i++) {
238 buffer = list[i];
239
240 while(buffer) {
241 buffer2 = buffer->next;
242 MEM_FREE(buffer, 1, VERT2ELEM);
243
244 buffer = buffer2;
245 }
246 }
247
248 MEM_FREE(list, data->n_total_vertices, VERT2ELEM *);
249
250 return;
251 }
252
253
254 /* Sets the boundary of all edges without neigbour to boundary type
255 * TYPE. if overwrite==true, then overwrite old boundary information,
256 * otherwise keep it, and only set type TYPE for previously
257 * unclassified boundary segments.
258 */
default_boundary(MACRO_DATA * data,U_CHAR type,bool overwrite)259 void default_boundary(MACRO_DATA *data, U_CHAR type, bool overwrite)
260 {
261 int i, dim = data->dim;
262
263 if(!data->boundary) {
264 data->boundary =
265 MEM_CALLOC(data->n_macro_elements*N_NEIGH(dim), BNDRY_TYPE);
266 }
267
268 for (i = 0; i < data->n_macro_elements * N_NEIGH(dim); i++) {
269 if (data->neigh[i] < 0 &&
270 (overwrite || data->boundary[i] == INTERIOR)) {
271 data->boundary[i] = type;
272 }
273 }
274 }
275
276 /* Transform the global representation of the wall-transformation to a
277 * local per-element representation (i.e. using the local numbering of
278 * the vertices of an element and its peridic neighbour).
279 */
_AI_compute_element_wall_transformations(MACRO_DATA * data)280 void _AI_compute_element_wall_transformations(MACRO_DATA *data)
281 {
282 int i, j, wtv, v, match, notmatched = 0, dim = data->dim;
283
284 memset(data->el_wall_vtx_trafos, 0,
285 sizeof(int)*N_WALLS(dim)*data->n_macro_elements);
286
287 for (i = 0; i < data->n_macro_elements; i++) {
288 for (j = 0; j < data->n_wall_vtx_trafos; j++) {
289 match = 0;
290 for (v = 0; v < N_VERTICES(dim); v++) {
291 for (wtv = 0; wtv < N_VERTICES(dim-1); wtv++) {
292 if (data->mel_vertices[i*N_VERTICES(dim)+v]
293 ==
294 data->wall_vtx_trafos[j][wtv][0]) {
295 break;
296 }
297 }
298 if (wtv < N_VERTICES(dim-1)) {
299 match++;
300 } else {
301 notmatched = v;
302 }
303 }
304 if (match == N_VERTICES(dim-1)) { /* wall-trafo j is for us */
305 data->el_wall_vtx_trafos[i*N_WALLS(dim)+notmatched] = j+1;
306 } else { /* test the inverse of this wall transformation */
307 match = 0;
308 for (v = 0; v < N_VERTICES(dim); v++) {
309 for (wtv = 0; wtv < N_VERTICES(dim-1); wtv++) {
310 if (data->mel_vertices[i*N_VERTICES(dim)+v]
311 ==
312 data->wall_vtx_trafos[j][wtv][1]) {
313 break;
314 }
315 }
316 if (wtv < N_VERTICES(dim-1)) {
317 match++;
318 } else {
319 notmatched = v;
320 }
321 }
322 if (match == N_VERTICES(dim-1)) { /* (wall-trafo j)^{-1} is for us */
323 data->el_wall_vtx_trafos[i*N_WALLS(dim)+notmatched] = -(j+1);
324 }
325 }
326 }
327 }
328 }
329
330 #include "macro_1d.c"
331 #if DIM_MAX > 1
332 #include "macro_2d.c"
333 #if DIM_MAX > 2
334 #include "macro_3d.c"
335 #endif
336 #endif
337
338
339 /*--------------------------------------------------------------------------*/
340 /* read data->neigh into mel[].neigh[] */
341 /* fill opp_vertex values and do a check on neighbour relations */
342 /*--------------------------------------------------------------------------*/
343
fill_neigh_info(MACRO_EL * mel,const MACRO_DATA * data)344 static void fill_neigh_info(MACRO_EL *mel, const MACRO_DATA *data)
345 {
346 FUNCNAME("check_neigh_info");
347 MACRO_EL *neigh;
348 int i, j, k, l, index, dim = data->dim;
349
350 for (i = 0; i < data->n_macro_elements; i++) {
351 for (j = 0; j < N_NEIGH(dim); j++) {
352 mel[i].neigh[j] =
353 ((index = data->neigh[NEIGH_IND(dim, i, j)]) >= 0) ? (mel+index) : NULL;
354 }
355 }
356
357 for (i = 0; i < data->n_macro_elements; i++) {
358 for (j = 0; j < N_NEIGH(dim); j++) {
359 for (l = 0; l < N_VERTICES(dim-1); l++) {
360 mel[i].neigh_vertices[j][l] = -1;
361 }
362 if ((neigh = mel[i].neigh[j])) {
363 if (data->n_wall_vtx_trafos) {
364 /* Strange things can happen with periodic meshes :) */
365 int vertices[N_VERTICES(DIM_MAX-1)] = { 0, };
366 int wt;
367
368 if ((wt = data->el_wall_vtx_trafos[NEIGH_IND(dim, i, j)])) {
369 if (wt > 0) {
370 for (k = 0; k < N_VERTICES(dim-1); k++) {
371 index = data->mel_vertices[VERT_IND(dim, i, (j+k+1)%(dim+1))];
372 for (l = 0; l < N_VERTICES(dim-1); l++) {
373 if (data->wall_vtx_trafos[wt-1][l][0] == index) {
374 vertices[k] = data->wall_vtx_trafos[wt-1][l][1];
375 }
376 }
377 }
378 } else {
379 for (k = 0; k < N_VERTICES(dim-1); k++) {
380 index = data->mel_vertices[VERT_IND(dim, i, (j+k+1)%(dim+1))];
381 for (l = 0; l < N_VERTICES(dim-1); l++) {
382 if (data->wall_vtx_trafos[-wt-1][l][1] == index) {
383 vertices[k] = data->wall_vtx_trafos[-wt-1][l][0];
384 }
385 }
386 }
387 }
388 } else {
389 for (k = 0; k < N_VERTICES(dim-1); k++) {
390 vertices[k] =
391 data->mel_vertices[VERT_IND(dim, i, (j+k+1)%(dim+1))];
392 }
393 }
394 k = opp_vertex(dim,
395 data->mel_vertices+neigh->index*N_VERTICES(dim),
396 vertices, wt ? mel[i].neigh_vertices[j] : NULL);
397 TEST_EXIT(k >= 0 && neigh->neigh[k] == mel+i,
398 "el %d is no neighbour of neighbour %d!\n",
399 mel[i].index, neigh->index);
400 } else {
401 for (k = 0; k < N_NEIGH(dim); k++)
402 if (neigh->neigh[k] == mel+i) break;
403
404 TEST_EXIT(k<N_NEIGH(dim), "el %d is no neighbour of neighbour %d!\n",
405 mel[i].index, neigh->index);
406
407 }
408 TEST_EXIT(data->opp_vertex == NULL ||
409 data->opp_vertex[i*N_NEIGH(dim)+j] == k,
410 "Inconsistent computations of opp_vertex!\n");
411 mel[i].opp_vertex[j] = k;
412 } else {
413 mel[i].opp_vertex[j] = -1;
414 }
415 }
416 }
417
418 return;
419 }
420
421
422 /*--------------------------------------------------------------------------*/
423 /* domain size */
424 /*--------------------------------------------------------------------------*/
425
calculate_size(MESH * mesh,const MACRO_DATA * data)426 static void calculate_size(MESH *mesh, const MACRO_DATA *data)
427 {
428 int i, j;
429
430 for (j = 0; j < DIM_OF_WORLD; j++)
431 {
432 mesh->bbox[0][j] = data->coords[0][j];
433 mesh->bbox[1][j] = data->coords[0][j];
434 }
435
436 for (i = 0; i < mesh->n_vertices; i++)
437 {
438 for (j = 0; j < DIM_OF_WORLD; j++)
439 {
440 mesh->bbox[0][j] = MIN(mesh->bbox[0][j], data->coords[i][j]);
441 mesh->bbox[1][j] = MAX(mesh->bbox[1][j], data->coords[i][j]);
442 }
443 }
444
445 AXPBY_DOW(1.0, mesh->bbox[1], -1.0, mesh->bbox[0], mesh->diam);
446
447 return;
448 }
449
450 /* read_comment(file): consume the current line starting at the first
451 * occurence of '#'
452 */
read_comment(FILE * file)453 static void read_comment(FILE *file)
454 {
455 int c;
456
457 for (;;) {
458 while (isspace(c = fgetc(file)));
459 if (c == '#') {
460 int waste = fscanf(file, "%*[^\r\n]"); /* <- Could this possibly work? */
461 (void)waste;
462 } else {
463 ungetc(c, file);
464 return;
465 }
466 }
467 }
468
469 /*--------------------------------------------------------------------------*/
470 /* read_indices() reads dim+1 indices from file into id[0-dim], */
471 /* returns true if dim+1 inputs arguments could be read successfully by */
472 /* fscanf(), else false */
473 /*--------------------------------------------------------------------------*/
474
read_indices(int dim,FILE * file,int id[])475 static int read_indices(int dim, FILE *file, int id[])
476 {
477 int i;
478
479 for (i = 0; i <= dim; i++) {
480 read_comment(file);
481 if (fscanf(file, "%d", id+i) != 1)
482 return false;
483 }
484 return true;
485 }
486
487 enum {
488 KEY_DIM = 0,
489 KEY_DOW,
490 KEY_NV,
491 KEY_NEL,
492 KEY_V_COORD,
493 KEY_EL_VERT,
494 KEY_EL_BND,
495 KEY_EL_NEIGH,
496 KEY_EL_TYPE,
497 KEY_N_WALL_VTX_TRAFOS,
498 KEY_WALL_VTX_TRAFOS,
499 KEY_N_WALL_TRAFOS,
500 KEY_WALL_TRAFOS,
501 KEY_EL_WALL_TRAFOS,
502 N_KEYS
503 };
504
505 #define N_MIN_KEYS 6
506 static const char *keys[N_KEYS] = {
507 "DIM", /* 0 */
508 "DIM_OF_WORLD", /* 1 */
509 "number of vertices", /* 2 */
510 "number of elements", /* 3 */
511 "vertex coordinates", /* 4 */
512 "element vertices", /* 5 */
513 "element boundaries", /* 6 */
514 "element neighbours", /* 7 */
515 "element type", /* 8 */
516 "number of wall vertex transformations", /* 9 */
517 "wall vertex transformations",/* 10 */
518 "number of wall transformations", /* 11 */
519 "wall transformations",/* 12 */
520 "element wall transformations", /* 13 */
521 };
522
get_key_no(const char * key)523 static int get_key_no(const char *key)
524 {
525 int i;
526
527 for (i = 0; i < N_KEYS; i++)
528 if (!strcmp(keys[i], key)) return(i);
529
530 return(-1);
531 }
532
read_key(const char * line)533 static const char *read_key(const char *line)
534 {
535 static char key[100];
536 char *k = key;
537
538 while(isspace(*line)) line++;
539 while((*k++ = *line++) != ':');
540 *--k = '\0';
541
542 return((const char *) key);
543 }
544
545 /*--------------------------------------------------------------------------*/
546 /* read_macro_data(): */
547 /* read macro triangulation from ascii file in ALBERTA format */
548 /* fills macro_data structure */
549 /* called by read_macro() */
550 /*--------------------------------------------------------------------------*/
551
read_macro_data(const char * filename)552 static MACRO_DATA *read_macro_data(const char *filename)
553 {
554 FUNCNAME("read_macro_data");
555 FILE *file;
556 MACRO_DATA *macro_data = NULL;
557 int dim, dow, nv, ne, nwt, nwvt, i, j, k, ind[N_VERTICES_MAX];
558 REAL dbl;
559 char name[128], line[256];
560 int line_no, n_keys, i_key, sort_key[N_KEYS];
561 int nv_key, ne_key, nwvt_key, nwt_key;
562 int key_def[N_KEYS] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
563 int blah;
564 const char *key;
565
566 TEST_EXIT(filename, "no file specified; filename NULL pointer\n");
567 TEST_EXIT(strlen(filename) < (unsigned int) 127,
568 "can only handle filenames up to 127 characters\n");
569
570 TEST_EXIT((file=fopen(filename, "r")), "cannot open file %s\n", filename);
571 strncpy(name, filename, 127);
572
573 /*--------------------------------------------------------------------------*/
574 /* looking for all keys in the macro file ... */
575 /*--------------------------------------------------------------------------*/
576
577 line_no = n_keys = 0;
578 while (fgets(line, sizeof(line), file))
579 {
580 line_no++;
581 if (*line == '#') continue;
582 if (!strchr(line, ':')) continue;
583 key = read_key(line);
584 i_key = get_key_no(key);
585 TEST_EXIT(i_key >= 0,
586 "file %s: must not contain key %s on line %d\n",
587 name, key, line_no);
588 TEST_EXIT(!key_def[i_key],
589 "file %s: key %s defined second time on line %d\n", name, key, line_no);
590
591 sort_key[n_keys++] = i_key;
592 key_def[i_key] = true;
593 }
594
595 for (i_key = 0; i_key < N_MIN_KEYS; i_key++)
596 {
597 for (j = 0; j < n_keys; j++)
598 if (sort_key[j] == i_key) break;
599 TEST_EXIT(j < n_keys,
600 "file %s: You do not have specified data for %s in %s\n",
601 name, keys[i_key]);
602
603 for (j = 0; j < n_keys; j++)
604 if (sort_key[j] == KEY_NV) break;
605 nv_key = j;
606 for (j = 0; j < n_keys; j++)
607 if (sort_key[j] == KEY_NEL) break;
608 ne_key = j;
609 for (j = 0; j < n_keys; j++)
610 if (sort_key[j] == KEY_N_WALL_VTX_TRAFOS) break;
611 nwvt_key = j;
612 for (j = 0; j < n_keys; j++)
613 if (sort_key[j] == KEY_N_WALL_TRAFOS) break;
614 nwt_key = j;
615
616 switch(i_key)
617 {
618 case KEY_DIM:
619 case KEY_DOW:
620 TEST_EXIT(sort_key[i_key] < 2,
621 "file %s: You have to specify DIM or DIM_OF_WORLD before all other data\n", name);
622 break;
623 case KEY_V_COORD:
624 TEST_EXIT(nv_key < i_key,
625 "file %s: Before reading data for %s, you have to specify the %s\n",
626 name, keys[KEY_V_COORD], keys[KEY_NV]);
627 break;
628 case KEY_EL_VERT:
629 TEST_EXIT(nv_key < i_key && ne_key < i_key,
630 "file %s: Before reading data for %s, "
631 "you have to specify the %s and %s\n",
632 name, keys[KEY_EL_VERT], keys[KEY_NEL], keys[KEY_NV]);
633 case KEY_EL_BND:
634 case KEY_EL_NEIGH:
635 case KEY_EL_TYPE:
636 TEST_EXIT(ne_key < i_key,
637 "file %s: Before reading data for %s, "
638 "you have to specify the %s\n",
639 name, keys[i_key], keys[KEY_NEL]);
640 break;
641 case KEY_N_WALL_VTX_TRAFOS:
642 break;
643 case KEY_WALL_VTX_TRAFOS:
644 TEST_EXIT(nwvt_key < i_key,
645 "file %s: Before reading data for %s, "
646 "you have to specify the %s\n",
647 name, keys[i_key], keys[KEY_N_WALL_VTX_TRAFOS]);
648 break;
649 case KEY_N_WALL_TRAFOS:
650 break;
651 case KEY_WALL_TRAFOS:
652 TEST_EXIT(nwt_key < i_key,
653 "file %s: Before reading data for %s, "
654 "you have to specify the %s\n",
655 name, keys[i_key], keys[KEY_N_WALL_TRAFOS]);
656 case KEY_EL_WALL_TRAFOS:
657 TEST_EXIT(ne_key < i_key,
658 "file %s: Before reading data for %s, "
659 "you have to specify the %s\n",
660 name, keys[i_key], keys[KEY_NEL]);
661 break;
662 }
663 }
664
665 for (i_key = 0; i_key < N_KEYS; i_key++)
666 key_def[i_key] = false;
667
668 /*--------------------------------------------------------------------------*/
669 /* and now, reading data ... */
670 /*--------------------------------------------------------------------------*/
671
672 rewind(file);
673
674 for (i_key = 0; i_key < n_keys; i_key++) {
675 read_comment(file);
676
677 switch(sort_key[i_key]) {
678 case KEY_DIM:
679 TEST_EXIT(fscanf(file, "%*s %d", &dim) == 1,
680 "file %s: cannot read DIM correctly\n", name);
681 TEST_EXIT(dim <= DIM_MAX,
682 "file %s: dimension = %d > DIM_MAX = %d\n",
683 name, dim, DIM_MAX);
684 key_def[KEY_DIM] = true;
685 break;
686 case KEY_DOW:
687 TEST_EXIT(fscanf(file, "%*s %d", &dow) == 1,
688 "file %s: cannot read DIM_OF_WORLD correctly\n", name);
689 TEST_EXIT(dow == DIM_OF_WORLD,
690 "file %s: dimension of world = %d != DIM_OF_WORLD = %d\n",
691 name, dow, DIM_OF_WORLD);
692
693 key_def[KEY_DOW] = true;
694 break;
695 case KEY_NV:
696 TEST_EXIT(fscanf(file, "%*s %*s %*s %d", &nv) == 1,
697 "file %s: cannot read number of vertices correctly\n", name);
698 TEST_EXIT(nv > 0,
699 "file %s: number of vertices = %d must be bigger than 0\n", name, nv);
700
701 key_def[KEY_NV] = true;
702
703 if(key_def[KEY_NEL])
704 macro_data = alloc_macro_data(dim, nv, ne);
705
706 break;
707 case KEY_NEL:
708 TEST_EXIT(fscanf(file, "%*s %*s %*s %d", &ne) == 1,
709 "file %s: cannot read number of elements correctly\n", name);
710 TEST_EXIT(ne > 0,
711 "file %s: number of elements = %d must be bigger than 0\n", name, ne);
712
713 key_def[KEY_NEL] = true;
714
715 if(key_def[KEY_NV])
716 macro_data = alloc_macro_data(dim, nv, ne);
717
718 break;
719 case KEY_V_COORD:
720 blah = fscanf(file, "%*s %*s");
721 (void)blah;
722 for (i = 0; i < nv; i++) {
723 for (j = 0; j < DIM_OF_WORLD; j++) {
724 read_comment(file);
725 TEST_EXIT(fscanf(file, "%lf", &dbl) == 1,
726 "file %s: error while reading coordinates, check file\n",
727 name);
728
729 macro_data->coords[i][j] = dbl;
730 }
731 }
732
733 key_def[KEY_V_COORD] = true;
734 break;
735 case KEY_EL_VERT:
736 blah = fscanf(file, "%*s %*s");
737 (void)blah;
738 /*--------------------------------------------------------------------------*/
739 /* global index of vertices for each single element */
740 /*--------------------------------------------------------------------------*/
741
742 for (i = 0; i < ne; i++) {
743 TEST_EXIT(read_indices(dim, file, ind),
744 "file %s: cannot read vertex indices of element %d\n",
745 name, i);
746
747 for (j = 0; j < N_VERTICES(dim); j++)
748 macro_data->mel_vertices[VERT_IND(dim, i, j)] = ind[j];
749 }
750
751 key_def[KEY_EL_VERT] = true;
752 break;
753 case KEY_EL_BND:
754 blah = fscanf(file, "%*s %*s");
755 (void)blah;
756
757 if(dim == 0)
758 ERROR_EXIT("Boundary types do not make sense in 0d!\n");
759 else {
760 /*--------------------------------------------------------------------------*/
761 /* read boundary type of each vertex/edge/face (in 1d/2d/3d) */
762 /*--------------------------------------------------------------------------*/
763 macro_data->boundary = MEM_ALLOC(ne*N_NEIGH(dim), BNDRY_TYPE);
764
765 for (i = 0; i < ne; i++) {
766 TEST_EXIT(read_indices(dim, file, ind),
767 "file %s: cannot read boundary types of element %d\n", name, i);
768
769 for(j = 0; j < N_NEIGH(dim); j++)
770 macro_data->boundary[NEIGH_IND(dim, i, j)] = (BNDRY_TYPE)ind[j];
771 }
772 }
773
774 key_def[KEY_EL_BND] = true;
775 break;
776 case KEY_EL_NEIGH:
777 blah = fscanf(file, "%*s %*s");
778 (void)blah;
779
780 if(dim == 0)
781 ERROR_EXIT("Neighbour indices do not make sense in 0d!\n");
782 else {
783 /*--------------------------------------------------------------------------*/
784 /* read neighbour indices: */
785 /*--------------------------------------------------------------------------*/
786 macro_data->neigh = MEM_ALLOC(ne*N_NEIGH(dim), int);
787
788 for (i = 0; i < ne; i++) {
789 TEST_EXIT(read_indices(dim, file, ind),
790 "file %s: cannot read neighbour info of element %d\n", name, i);
791
792 for(j = 0; j < N_NEIGH(dim); j++)
793 macro_data->neigh[NEIGH_IND(dim, i, j)] = ind[j];
794 }
795 }
796
797 key_def[KEY_EL_NEIGH] = true;
798 break;
799 case KEY_EL_TYPE:
800 blah = fscanf(file, "%*s %*s");
801 (void)blah;
802 /*--------------------------------------------------------------------------*/
803 /* el_type is handled just like bound and neigh above */
804 /*--------------------------------------------------------------------------*/
805
806 if(dim < 3) {
807 WARNING("File %s: element type only used in 3d; will ignore data for el_type\n", name);
808 }
809 #if DIM_MAX > 2
810 else {
811 macro_data->el_type = MEM_ALLOC(ne, U_CHAR);
812
813 for (i = 0; i < ne; i++) {
814 read_comment(file);
815 TEST_EXIT(fscanf(file, "%d", &j) == 1,
816 "file %s: cannot read el_type of element %d\n", name, i);
817
818 macro_data->el_type[i] = (U_CHAR) j;
819 }
820 }
821 #endif
822 key_def[KEY_EL_TYPE] = true;
823 break;
824 case KEY_N_WALL_VTX_TRAFOS:
825 TEST_EXIT(fscanf(file, "%*s %*s %*s %*s %*s %d", &nwvt) == 1,
826 "file %s: cannot read number of "
827 "wall vertex transformations correctly\n",
828 name);
829 TEST_EXIT(nwvt > 0,
830 "file %s: number of wall transformations = %d must be bigger than 0\n",
831 name, nwvt);
832
833 macro_data->n_wall_vtx_trafos = nwvt;
834 macro_data->wall_vtx_trafos = (int (*)[N_VERTICES(DIM_MAX-1)][2])
835 MEM_ALLOC(nwvt*sizeof(*macro_data->wall_vtx_trafos), char);
836 macro_data->el_wall_vtx_trafos = MEM_ALLOC(ne*N_WALLS(dim), int);
837
838 key_def[KEY_N_WALL_VTX_TRAFOS] = true;
839 break;
840 case KEY_WALL_VTX_TRAFOS:
841 blah = fscanf(file, "%*s %*s %*s");
842 (void)blah;
843 for (i = 0; i < nwvt; i++) {
844 TEST_EXIT(read_indices(dim*2-1, file, ind),
845 "file %s: cannot read wall transformation %d\n", name, i);
846 for (j = 0; j < N_VERTICES(dim-1); j++) {
847 macro_data->wall_vtx_trafos[i][j][0] = ind[2*j];
848 macro_data->wall_vtx_trafos[i][j][1] = ind[2*j+1];
849 }
850 }
851 key_def[KEY_WALL_VTX_TRAFOS] = true;
852 break;
853 case KEY_N_WALL_TRAFOS:
854 TEST_EXIT(fscanf(file, "%*s %*s %*s %*s %d", &nwt) == 1,
855 "file %s: cannot read number of "
856 "wall transformations correctly\n",
857 name);
858 TEST_EXIT(nwt > 0,
859 "file %s: number of wall transformations = %d must be bigger than 0\n",
860 name, nwt);
861
862 macro_data->n_wall_trafos = nwt;
863 macro_data->wall_trafos = MEM_ALLOC(nwt, AFF_TRAFO);
864 key_def[KEY_N_WALL_TRAFOS] = true;
865 break;
866 case KEY_WALL_TRAFOS:
867 blah = fscanf(file, "%*s %*s");
868 (void)blah;
869 for (i = 0; i < nwt; i++) {
870 for (j = 0; j < (DIM_OF_WORLD+1); j++) {
871 for (k = 0; k < (DIM_OF_WORLD+1); k++) {
872 read_comment(file);
873 TEST_EXIT(fscanf(file, "%lf", &dbl) == 1,
874 "file %s: error while reading wall transformation\n",
875 name);
876 if (j < DIM_OF_WORLD) {
877 if (k < DIM_OF_WORLD) {
878 macro_data->wall_trafos[i].M[j][k] = dbl;
879 } else {
880 macro_data->wall_trafos[i].t[j] = dbl;
881 }
882 }
883 /* discard the last line (projective co-ordinates) */
884 }
885 }
886 }
887 key_def[KEY_WALL_TRAFOS] = true;
888 break;
889 case KEY_EL_WALL_TRAFOS:
890 blah = fscanf(file, "%*s %*s %*s");
891 (void)blah;
892
893 if(dim == 0) {
894 ERROR_EXIT("Wall transformations do not make sense in 0d!\n");
895 } else {
896 macro_data->el_wall_trafos = MEM_ALLOC(ne*N_NEIGH(dim), int);
897
898 for (i = 0; i < ne; i++) {
899 TEST_EXIT(read_indices(dim, file, ind),
900 "file %s: "
901 "cannot read the wall transformations for element %d\n",
902 name, i);
903
904 for(j = 0; j < N_NEIGH(dim); j++) {
905 macro_data->el_wall_trafos[NEIGH_IND(dim, i, j)] = ind[j];
906 }
907 }
908 }
909 key_def[KEY_EL_WALL_TRAFOS] = true;
910 break;
911 break;
912 }
913 }
914 fclose(file);
915
916 return(macro_data);
917 }
918
919 /*--------------------------------------------------------------------------*/
920 /* read macro triangulation from file "filename" into macro_data data in */
921 /* native binary format */
922 /*--------------------------------------------------------------------------*/
923
read_macro_data_bin(const char * name)924 static MACRO_DATA *read_macro_data_bin(const char *name)
925 {
926 FUNCNAME("read_macro_data_bin");
927 FILE *file;
928 MACRO_DATA *macro_data;
929 int dim, i, length, nv, ne;
930 char *s;
931 char record_written;
932 int blah;
933
934 TEST_EXIT(file = fopen(name, "rb"), "cannot open file %s\n", name);
935
936 length = MAX(strlen(ALBERTA_VERSION)+1, 21);
937 s = MEM_ALLOC(length, char);
938
939 blah = fread(s, sizeof(char), length, file);
940 (void)blah;
941 TEST_EXIT(!strncmp(s, "ALBERTA", 6), "file %s: unknown file id:\"%s\"\n", name, s);
942
943 MEM_FREE(s, length, char);
944
945 blah = fread(&i, sizeof(int), 1, file);
946 (void)blah;
947 TEST_EXIT(i == sizeof(REAL), "file %s: wrong sizeof(REAL) %d\n", name, i);
948
949 blah = fread(&dim, sizeof(int), 1, file);
950 (void)blah;
951 TEST_EXIT(dim <= DIM_OF_WORLD,
952 "file %s: dimension = %d > DIM_MAX = %d\n",
953 name, dim, DIM_MAX);
954
955 blah = fread(&i, sizeof(int), 1, file);
956 (void)blah;
957 TEST_EXIT(i == DIM_OF_WORLD,
958 "file %s: dimension of world = %d != DIM_OF_WORLD = %d\n",
959 name, i, DIM_OF_WORLD);
960
961 blah = fread(&nv, sizeof(int), 1, file);
962 (void)blah;
963 TEST_EXIT(nv > 0,
964 "file %s: number of vertices = %d must be bigger than 0\n", name, nv);
965
966 blah = fread(&ne, sizeof(int), 1, file);
967 (void)blah;
968 TEST_EXIT(ne > 0,
969 "file %s: number of elements = %d must be bigger than 0\n", name, ne);
970
971 macro_data = alloc_macro_data(dim, nv, ne);
972
973 blah = fread(macro_data->coords, sizeof(REAL_D), nv, file);
974 blah = fread(
975 macro_data->mel_vertices, sizeof(int), N_VERTICES(dim) * ne, file);
976
977 blah = fread(&record_written, sizeof(char), 1, file);
978 if(record_written)
979 {
980 macro_data->boundary = MEM_ALLOC(ne*N_NEIGH(dim), BNDRY_TYPE);
981 blah = fread(
982 macro_data->boundary, sizeof(BNDRY_TYPE), N_NEIGH(dim) * ne, file);
983 }
984
985 blah = fread(&record_written, sizeof(char), 1, file);
986 if(record_written) {
987 macro_data->neigh = MEM_ALLOC(ne*N_NEIGH(dim), int);
988 blah = fread(macro_data->neigh, sizeof(int), N_NEIGH(dim) * ne, file);
989 }
990
991 #if DIM_MAX > 2
992 if(dim == 3) {
993 int nonsense;
994 (void)nonsense;
995 nonsense = fread(&record_written, sizeof(char), 1, file);
996 if (record_written) {
997 macro_data->el_type = MEM_ALLOC(ne, U_CHAR);
998 nonsense = fread(macro_data->el_type, sizeof(U_CHAR), ne, file);
999 }
1000 }
1001 #endif
1002
1003 s = MEM_ALLOC(5, char);
1004
1005 TEST_EXIT(fread(s, sizeof(char), 4, file) == 4,
1006 "file %s: problem while reading FILE END MARK\n", name);
1007
1008 TEST_EXIT(!strncmp(s, "EOF.", 4), "file %s: no FILE END MARK\n", name);
1009 MEM_FREE(s, 5, char);
1010
1011 fclose(file);
1012
1013 return(macro_data);
1014 }
1015
1016
1017 /*--------------------------------------------------------------------------*/
1018 /* Some routines needed for interaction with xdr-files */
1019 /* WARNING: These will need to be adapted if ALBERTA data types REAL, REAL_D*/
1020 /* , etc. change! */
1021 /*--------------------------------------------------------------------------*/
1022
1023 static int xdr_dim = 0;
1024
xdr_REAL(XDR * xdr,REAL * rp)1025 static bool_t xdr_REAL(XDR *xdr, REAL *rp)
1026 {
1027 return (xdr_double(xdr, rp));
1028 }
1029
ALBERTA_DEFUNUSED(static bool_t xdr_U_CHAR (XDR * xdr,U_CHAR * ucp))1030 ALBERTA_DEFUNUSED(static bool_t xdr_U_CHAR(XDR *xdr, U_CHAR *ucp))
1031 {
1032 return (xdr_u_char(xdr, ucp));
1033 }
1034
xdr_S_CHAR(XDR * xdr,S_CHAR * cp)1035 static bool_t xdr_S_CHAR(XDR *xdr, S_CHAR *cp)
1036 {
1037 return (xdr_char(xdr,(char *)cp));
1038 }
1039
xdr_REAL_D(XDR * xdr,REAL_D * dp)1040 static bool_t xdr_REAL_D(XDR *xdr, REAL_D *dp)
1041 {
1042 return (xdr_vector(xdr, (char *)dp, DIM_OF_WORLD, sizeof(REAL), (xdrproc_t) xdr_REAL));
1043 }
1044
1045 #if __APPLE_CC__
1046 # define _XDR_PTR_T_ void *
1047 #else
1048 # define _XDR_PTR_T_ char *
1049 #endif
1050
read_xdr_file(_XDR_PTR_T_ xdr_file,_XDR_PTR_T_ buffer,int size)1051 static int read_xdr_file(_XDR_PTR_T_ xdr_file, _XDR_PTR_T_ buffer, int size)
1052 {
1053 return fread(buffer, (size_t)size, 1, (FILE *)xdr_file) == 1 ? size : 0;
1054 }
1055
write_xdr_file(_XDR_PTR_T_ xdr_file,_XDR_PTR_T_ buffer,int size)1056 static int write_xdr_file(_XDR_PTR_T_ xdr_file, _XDR_PTR_T_ buffer, int size)
1057 {
1058 return fwrite(buffer, (size_t)size, 1, (FILE *)xdr_file) == 1 ? size : 0;
1059 }
1060
1061
xdr_open_file(const char * filename,enum xdr_op mode)1062 static XDR *xdr_open_file(const char *filename, enum xdr_op mode)
1063 {
1064 XDR *xdr;
1065 FILE *xdr_file;
1066
1067 if (!(xdr = MEM_ALLOC(1, XDR)))
1068 {
1069 ERROR("can't allocate memory for xdr pointer.\n");
1070
1071 return NULL;
1072 }
1073
1074 if ((xdr_file = fopen(filename, mode == XDR_DECODE ? "r": "w")))
1075 {
1076 xdrrec_create(xdr, 65536, 65536, (caddr_t) xdr_file,
1077 read_xdr_file, write_xdr_file);
1078
1079 xdr->x_op = mode;
1080 xdr->x_public = (caddr_t)xdr_file;
1081
1082 if (mode == XDR_DECODE)
1083 xdrrec_skiprecord(xdr);
1084
1085 return xdr;
1086 }
1087 else
1088 {
1089 ERROR("error opening xdr file.\n");
1090
1091 MEM_FREE(xdr, 1, XDR);
1092
1093 return NULL;
1094 }
1095 }
1096
1097
xdr_close_file(XDR * xdr)1098 static int xdr_close_file(XDR *xdr)
1099 {
1100 if (!xdr)
1101 {
1102 ERROR("NULL xdr pointer.\n");
1103 return 0;
1104 }
1105
1106 if (xdr->x_op == XDR_ENCODE)
1107 xdrrec_endofrecord(xdr, 1);
1108
1109 if (fclose((FILE *) xdr->x_public))
1110 ERROR("error closing file.\n");
1111
1112 xdr_destroy(xdr);
1113
1114 MEM_FREE(xdr, 1, XDR);
1115
1116 return 1;
1117 }
1118
1119 /*--------------------------------------------------------------------------*/
1120 /* read_macro_data_xdr(): */
1121 /* read macro triangulation from file in xdr-format */
1122 /* fills macro_data structure */
1123 /* called by ....? */
1124 /*--------------------------------------------------------------------------*/
1125
read_macro_data_xdr(const char * name)1126 static MACRO_DATA *read_macro_data_xdr(const char *name)
1127 {
1128 FUNCNAME("read_macro_data_xdr");
1129 XDR *xdrp;
1130 MACRO_DATA *macro_data;
1131 int length, dow, nv, ne, size;
1132 char *s;
1133 bool_t record_written;
1134 caddr_t array_loc;
1135
1136
1137 TEST_EXIT(name, "no file specified; filename NULL pointer\n");
1138
1139 if (!(xdrp = xdr_open_file(name, XDR_DECODE)))
1140 ERROR_EXIT("cannot open file %s\n", name);
1141
1142 length = MAX(strlen(ALBERTA_VERSION)+1, 21); /* length with terminating \0 */
1143 s = MEM_ALLOC(length, char);
1144
1145 TEST_EXIT(xdr_string(xdrp, &s, length), "file %s: could not read file id\n", name);
1146 TEST_EXIT(!strncmp(s, "ALBERTA", 6), "file %s: unknown file id: \"%s\"\n", name, s);
1147
1148 MEM_FREE(s, length, char);
1149
1150 TEST_EXIT(xdr_int(xdrp, &xdr_dim),
1151 "file %s: could not read dimension correctly\n", name);
1152 TEST_EXIT(xdr_dim <= DIM_MAX,
1153 "file %s: dimension = %d > DIM_MAX = %d\n",
1154 name, xdr_dim, DIM_MAX);
1155
1156
1157 TEST_EXIT(xdr_int(xdrp, &dow), "file %s: could not read dimension of world correctly\n", name);
1158 TEST_EXIT(dow == DIM_OF_WORLD,
1159 "file %s: dimension of world = %d != DIM_OF_WORLD = %d\n",
1160 name, dow, DIM_OF_WORLD);
1161
1162 TEST_EXIT(xdr_int(xdrp, &nv),
1163 "file %s: cannot read number of vertices correctly\n", name);
1164 TEST_EXIT(nv > 0,
1165 "file %s: number of vertices = %d must be bigger than 0\n", name, nv);
1166
1167 TEST_EXIT(xdr_int(xdrp, &ne),
1168 "file %s: cannot read number of elements correctly\n", name);
1169 TEST_EXIT(ne > 0,
1170 "file %s: number of elements = %d must be bigger than 0\n", name, ne);
1171
1172 macro_data = alloc_macro_data(xdr_dim, nv, ne);
1173
1174 array_loc=(caddr_t) macro_data->coords;
1175 TEST_EXIT(xdr_array(xdrp, &array_loc, (u_int *) &nv, (u_int) nv, sizeof(REAL_D), (xdrproc_t) xdr_REAL_D),
1176 "file %s: error while reading coordinates, check file\n", name);
1177
1178 array_loc=(caddr_t) macro_data->mel_vertices;
1179 size = ne * N_VERTICES(xdr_dim);
1180 TEST_EXIT(xdr_array(xdrp, &array_loc, (u_int *) &size, (u_int) size,
1181 sizeof(int), (xdrproc_t) xdr_int),
1182 "file %s: cannot read vertex indices\n", name);
1183
1184 TEST_EXIT(xdr_bool(xdrp, &record_written),
1185 "file %s: could not determine whether to allocate memory for boundaries\n", name);
1186 if(record_written)
1187 {
1188 macro_data->boundary = MEM_ALLOC(ne*N_NEIGH(xdr_dim), BNDRY_TYPE);
1189
1190 array_loc=(caddr_t) macro_data->boundary;
1191 size = ne * N_NEIGH(xdr_dim);
1192 TEST_EXIT(xdr_array(xdrp, &array_loc, (u_int *) &size, (u_int) size,
1193 sizeof(BNDRY_TYPE), (xdrproc_t) xdr_S_CHAR),
1194 "file %s: could not read boundary types\n", name);
1195 }
1196
1197 TEST_EXIT(xdr_bool(xdrp, &record_written),
1198 "file %s: could not determine whether to allocate memory for neighbours\n", name);
1199 if(record_written)
1200 {
1201 macro_data->neigh = MEM_ALLOC(ne*N_NEIGH(xdr_dim), int);
1202
1203 array_loc=(caddr_t) macro_data->neigh;
1204 size = ne * N_NEIGH(xdr_dim);
1205 TEST_EXIT(xdr_array(xdrp, &array_loc, (u_int *) &size, (u_int) size,
1206 sizeof(int), (xdrproc_t) xdr_int),
1207 "file %s: could not read neighbor info\n", name);
1208 }
1209
1210 #if DIM_MAX > 2
1211 if(xdr_dim == 3) {
1212 TEST_EXIT(xdr_bool(xdrp, &record_written),
1213 "file %s: could not determine whether to allocate memory for element types\n", name);
1214 if(record_written)
1215 {
1216 macro_data->el_type = MEM_ALLOC(ne, U_CHAR);
1217
1218 array_loc=(caddr_t) macro_data->el_type;
1219 TEST_EXIT(xdr_array(xdrp, &array_loc, (u_int *) &ne, (u_int) ne, sizeof(U_CHAR), (xdrproc_t) xdr_U_CHAR),
1220 "file %s: cannot read element types\n", name);
1221 }
1222 }
1223 #endif
1224
1225 xdr_close_file(xdrp);
1226
1227 return(macro_data);
1228 }
1229
1230 /*--------------------------------------------------------------------------*/
1231 /* supported file types for macro data files: */
1232 /*--------------------------------------------------------------------------*/
1233
1234 typedef enum {ascii_format, binary_format, xdr_format} macro_format;
1235
1236 /*--------------------------------------------------------------------------*/
1237 /* read macro triangulation from file "filename" */
1238 /*--------------------------------------------------------------------------*/
1239
read_macro_master(const char * filename,macro_format format)1240 static MACRO_DATA *read_macro_master(const char *filename,
1241 macro_format format)
1242 {
1243 FUNCNAME("read_macro_master");
1244 MACRO_DATA *macro_data = NULL;
1245 char filenew[1024];
1246
1247 TEST_EXIT(filename, "no file specified; filename NULL pointer\n");
1248
1249 switch (format) {
1250 case ascii_format:
1251 macro_data = read_macro_data(filename);
1252 break;
1253 case binary_format:
1254 macro_data = read_macro_data_bin(filename);
1255 break;
1256 case xdr_format:
1257 macro_data = read_macro_data_xdr(filename);
1258 }
1259
1260 if (macro_data->n_wall_vtx_trafos > 0)
1261 _AI_compute_element_wall_transformations(macro_data);
1262 if (!macro_data->neigh && macro_data->dim > 0)
1263 compute_neigh_fast(macro_data);
1264 if (!macro_data->boundary && macro_data->dim > 0) {
1265 default_boundary(macro_data, DIRICHLET /* default boundary type */, true);
1266 }
1267
1268 snprintf(filenew, 1024, "%s.new", filename);
1269 macro_test(macro_data, filenew);
1270
1271 return macro_data;
1272 }
1273
1274
cleanup_write_macro(MACRO_DATA * data,DOF_INT_VEC * dof_vert_ind,TRAVERSE_STACK * stack)1275 static void cleanup_write_macro(MACRO_DATA *data, DOF_INT_VEC *dof_vert_ind,
1276 TRAVERSE_STACK *stack)
1277 {
1278
1279 if (data) free_macro_data(data);
1280 free_dof_int_vec(dof_vert_ind);
1281 if (stack)
1282 free_traverse_stack(stack);
1283
1284 return;
1285 }
1286
1287
1288 /*--------------------------------------------------------------------------*/
1289 /* mesh2macro_data(): counterpart to macro_data2mesh below. This routine */
1290 /* converts the information stored in the leaf elements of mesh to the raw */
1291 /* data type MACRO_DATA. */
1292 /*--------------------------------------------------------------------------*/
1293
mesh2macro_data(MESH * mesh)1294 MACRO_DATA *mesh2macro_data(MESH *mesh)
1295 {
1296 FUNCNAME("mesh2macro_data");
1297 MACRO_DATA *data;
1298 FLAGS fill_flag;
1299 PARAMETRIC *parametric;
1300 const DOF_ADMIN *admin;
1301 const FE_SPACE *fe_space;
1302 const int n_dof[N_NODE_TYPES] = {1, 0, 0, 0};
1303 DOF_INT_VEC *dof_vert_ind;
1304 int dim = mesh->dim, n0, ne, nv, i, *vert_ind = NULL;
1305 #if DIM_MAX > 2
1306 U_CHAR write_el_type;
1307 #endif
1308 static const REAL_B vertex_bary[N_VERTICES_LIMIT] = {
1309 INIT_BARY_3D(1.0, 0.0, 0.0, 0.0),
1310 INIT_BARY_3D(0.0, 1.0, 0.0, 0.0),
1311 INIT_BARY_3D(0.0, 0.0, 1.0, 0.0),
1312 INIT_BARY_3D(0.0, 0.0, 0.0, 1.0)
1313 };
1314
1315 /* we get a real vertex-only admin; this way it is easy to
1316 * manipulate, e.g., per-vertex data, and if we re-generate a mesh
1317 * from this macro data then the new mesh will have the same vertex
1318 * ordering as the vertex admin.
1319 */
1320 dof_compress(mesh);
1321
1322 fe_space = get_dof_space(mesh, "mesh2macro_data", n_dof, ADM_FLAGS_DFLT);
1323 admin = fe_space->admin;
1324
1325 n0 = admin->n0_dof[VERTEX];
1326
1327 parametric = mesh->parametric;
1328
1329 dof_vert_ind = get_dof_int_vec("vertex indices", fe_space);
1330 GET_DOF_VEC(vert_ind, dof_vert_ind);
1331 FOR_ALL_DOFS(admin, vert_ind[dof] = -1);
1332
1333 data = alloc_macro_data(dim, mesh->n_vertices, mesh->n_elements);
1334
1335 nv = ne = 0;
1336 #if DIM_MAX > 2
1337 write_el_type = false;
1338 #endif
1339
1340 /*--------------------------------------------------------------------------*/
1341 /* The first pass counts elements and vertices, checks these against the */
1342 /* entries of mesh->n_elements, mesh->n_vertices, and fills data->coords. */
1343 /* A check on whether an element has nonzero el_type is also done. */
1344 /*--------------------------------------------------------------------------*/
1345 TRAVERSE_FIRST(mesh, -1,
1346 CALL_LEAF_EL|FILL_COORDS|FILL_NEIGH|FILL_NON_PERIODIC) {
1347
1348 if (parametric) {
1349 if (parametric->init_element(el_info, parametric)) {
1350 /* really parametric, otherwise el_info->coord already
1351 * contains the coordinates. */
1352 parametric->coord_to_world(el_info, NULL, N_VERTICES(dim),
1353 vertex_bary, (REAL_D *)el_info->coord);
1354 }
1355 }
1356
1357 for (i = 0; i < N_VERTICES(dim); i++) {
1358 if (vert_ind[el_info->el->dof[i][n0]] == -1) {
1359 /*--------------------------------------------------------------------------*/
1360 /* assign a global index to each vertex */
1361 /*--------------------------------------------------------------------------*/
1362 vert_ind[el_info->el->dof[i][n0]] = el_info->el->dof[i][n0];
1363
1364 COPY_DOW(el_info->coord[i], data->coords[el_info->el->dof[i][n0]]);
1365
1366 nv++;
1367
1368 if(nv > mesh->n_vertices)
1369 {
1370 cleanup_write_macro(data, dof_vert_ind, stack);
1371 free_fe_space(fe_space);
1372 ERROR("mesh %s: n_vertices (==%d) is too small! Writing aborted\n",
1373 mesh->name, mesh->n_vertices);
1374 return NULL;
1375 }
1376 }
1377 }
1378
1379 ne++;
1380
1381 if(ne > mesh->n_elements) {
1382 cleanup_write_macro(data, dof_vert_ind, stack);
1383 ERROR("mesh %s: n_elements (==%d) is too small! Writing aborted\n",
1384 mesh->name, mesh->n_elements);
1385 return(NULL);
1386 }
1387 #if DIM_MAX > 2
1388 if(dim == 3 && el_info->el_type) write_el_type = true;
1389 #endif
1390 } TRAVERSE_NEXT();
1391
1392 if(ne < mesh->n_elements) {
1393 cleanup_write_macro(data, dof_vert_ind, NULL);
1394 free_fe_space(fe_space);
1395 ERROR("mesh %s: n_elements (==%d) is too large: only %d leaf elements counted -- writing aborted\n", mesh->name, mesh->n_elements, ne);
1396 return(NULL);
1397 }
1398
1399 if(nv < mesh->n_vertices) {
1400 cleanup_write_macro(data, dof_vert_ind, NULL);
1401 free_fe_space(fe_space);
1402 ERROR("mesh %s: n_vertices (==%d) is too large: only %d vertices counted -- allocation of macro data aborted\n", mesh->name, mesh->n_vertices, nv);
1403 return(NULL);
1404 }
1405
1406 if(dim > 0) {
1407 data->boundary = MEM_ALLOC(ne*N_NEIGH(dim), BNDRY_TYPE);
1408 }
1409
1410 #if DIM_MAX > 2
1411 if(write_el_type) {
1412 data->el_type = MEM_ALLOC(ne, U_CHAR);
1413 }
1414 #endif
1415
1416 if (mesh->is_periodic) {
1417 data->el_wall_vtx_trafos = MEM_ALLOC(ne*N_WALLS(dim), int);
1418
1419 if (mesh->n_wall_trafos) {
1420 data->n_wall_trafos = mesh->n_wall_trafos / 2;
1421 data->wall_trafos = MEM_ALLOC(data->n_wall_trafos, AFF_TRAFO);
1422 TEST_EXIT(mesh->n_wall_trafos % 2 == 0,
1423 "odd number of wall-transformations????");
1424 for (i = 0; i < data->n_wall_trafos; i++) {
1425 data->wall_trafos[i] = *mesh->wall_trafos[2*i];
1426 }
1427 data->el_wall_trafos =
1428 MEM_CALLOC(N_WALLS(dim) * data->n_macro_elements, int);
1429 }
1430 }
1431
1432 #if ALBERTA_DEBUG
1433 char *mel_comment = MEM_CALLOC(ne*80, char);
1434 data->mel_comment = MEM_CALLOC(ne, char *);
1435 #endif
1436
1437 ne = 0;
1438
1439 /* The second pass assigns mel_vertices, boundary, and if necessary el_type
1440 *
1441 * We use FILL_NON_PERIODIC to recover boundary information
1442 * potentially attached to periodic walls. The neighbourhood
1443 * connectivity of the mesh is computed via compute_neigh_fast()
1444 * later, so FILL_NON_PERIODIC does not hurt here.
1445 */
1446 fill_flag =
1447 CALL_LEAF_EL|FILL_COORDS|FILL_MACRO_WALLS|FILL_NEIGH|FILL_NON_PERIODIC;
1448 TRAVERSE_FIRST(mesh, -1, fill_flag) {
1449
1450 for (i = 0; i < N_VERTICES(dim); i++)
1451 data->mel_vertices[VERT_IND(dim, ne, i)] =
1452 vert_ind[el_info->el->dof[i][n0]];
1453
1454 if(dim > 0) {
1455 for (i = 0; i < N_NEIGH(dim); i++) {
1456 data->boundary[NEIGH_IND(dim, ne, i)] = wall_bound(el_info, i);
1457 }
1458 }
1459
1460 #if DIM_MAX > 2
1461 if(write_el_type) data->el_type[ne] = el_info->el_type;
1462 #endif
1463
1464 if (mesh->is_periodic) {
1465 int w, v;
1466
1467
1468 for (w = 0; w < N_WALLS(dim); w++) {
1469 EL *neigh = el_info->neigh[w], *el = el_info->el;
1470
1471 if (neigh) {
1472 const int *wall_ind;
1473 const int *wall_ind_n;
1474 int ov = el_info->opp_vertex[w];
1475
1476 wall_ind =
1477 sorted_wall_vertices(dim, w, wall_orientation(dim, el, w));
1478 wall_ind_n =
1479 sorted_wall_vertices(dim, ov, wall_orientation(dim, neigh, ov));
1480
1481 /* If this is a periodic wall, then _ALL_ vertices must be
1482 * different from the neighbour vertices, so we just need to
1483 * check one.
1484 */
1485 if (el->dof[wall_ind[0]][n0] != neigh->dof[wall_ind_n[0]][n0]) {
1486 /* Now construct the combinatorical wall vertex transformations.
1487 */
1488 int nwt = data->n_wall_vtx_trafos++;
1489
1490 if (nwt % 100 == 0) {
1491 data->wall_vtx_trafos = (int (*)[N_VERTICES(DIM_MAX-1)][2])
1492 MEM_REALLOC(data->wall_vtx_trafos,
1493 nwt*sizeof(*data->wall_vtx_trafos),
1494 (nwt+100)*sizeof(*data->wall_vtx_trafos),
1495 char);
1496 }
1497 for (v = 0; v < N_VERTICES(dim-1); v++) {
1498 data->wall_vtx_trafos[nwt][v][0] =
1499 vert_ind[el->dof[wall_ind[v]][n0]];
1500 data->wall_vtx_trafos[nwt][v][1] =
1501 vert_ind[neigh->dof[wall_ind_n[v]][n0]];
1502 }
1503 data->el_wall_vtx_trafos[ne*N_WALLS(dim)+w] = nwt + 1;
1504
1505 if (data->el_wall_trafos) {
1506 int mwall = el_info->macro_wall[w];
1507 AFF_TRAFO *wt = el_info->macro_el->wall_trafo[mwall];
1508 int i;
1509
1510 for (i = 0; i < mesh->n_wall_trafos; i++) {
1511 if (wt == mesh->wall_trafos[i]) {
1512 break;
1513 }
1514 }
1515 DEBUG_TEST_EXIT(i < mesh->n_wall_trafos, "Unknown wall-trafo?!");
1516 if (i % 2 == 0) {
1517 data->el_wall_trafos[ne*N_WALLS(dim)+w] = +(i / 2 + 1);
1518 } else {
1519 data->el_wall_trafos[ne*N_WALLS(dim)+w] = -(i / 2 + 1);
1520 }
1521 }
1522
1523 } else {
1524 data->el_wall_vtx_trafos[ne*N_WALLS(dim)+w] = 0;
1525 }
1526 } else {
1527 data->el_wall_vtx_trafos[ne*N_WALLS(dim)+w] = 0;
1528 }
1529 }
1530 }
1531
1532 #if ALBERTA_DEBUG
1533 data->mel_comment[ne] = mel_comment + ne*80;
1534 snprintf(data->mel_comment[ne], 80, " Id: %d; Level: %d;",
1535 el_info->el->index, el_info->level);
1536 #endif
1537
1538 ++ne;
1539 } TRAVERSE_NEXT();
1540
1541 /*--------------------------------------------------------------------------*/
1542 /* Finally, we compute neighbour information. This seems to be the easiest */
1543 /* solution, since neighbour information in ALBERTA is only available as */
1544 /* pointers. */
1545 /*--------------------------------------------------------------------------*/
1546
1547 if(dim > 0) {
1548 compute_neigh_fast(data);
1549 }
1550
1551 /* Cleanup the wall-transformations (need only half of them) */
1552 if (mesh->is_periodic) {
1553 int neigh, oppv, w;
1554 int nwt = 0;
1555 int (*old_vtx_wts)[N_VERTICES(DIM_MAX-1)][2] = data->wall_vtx_trafos;
1556
1557 data->wall_vtx_trafos = (int (*)[N_VERTICES(DIM_MAX-1)][2])
1558 MEM_ALLOC(data->n_wall_vtx_trafos/2*sizeof(*data->wall_vtx_trafos), char);
1559 for (ne = 0; ne < data->n_macro_elements; ne++) {
1560 for (w = 0; w < N_WALLS(dim); w++) {
1561 if (data->el_wall_vtx_trafos[ne*N_WALLS(dim)+w] > 0) {
1562 memcpy(data->wall_vtx_trafos[nwt],
1563 old_vtx_wts[data->el_wall_vtx_trafos[ne*N_WALLS(dim)+w]-1],
1564 sizeof(*data->wall_vtx_trafos));
1565 data->el_wall_vtx_trafos[ne*N_WALLS(dim)+w] = nwt+1;
1566 neigh = data->neigh[ne*N_NEIGH(dim)+w];
1567 oppv = data->opp_vertex[ne*N_NEIGH(dim)+w];
1568 data->el_wall_vtx_trafos[neigh*N_WALLS(dim)+oppv] = -(nwt+1);
1569 ++nwt;
1570 }
1571 }
1572 }
1573 TEST_EXIT(data->n_wall_vtx_trafos == 2*nwt,
1574 "Wall transformation do not seem to be reflexive!\n");
1575 MEM_FREE(old_vtx_wts,
1576 (data->n_wall_vtx_trafos + 100 - 1) / 100 * 100, char);
1577 data->n_wall_vtx_trafos = nwt;
1578 }
1579
1580 cleanup_write_macro(NULL, dof_vert_ind, NULL);
1581
1582 free_fe_space(fe_space);
1583
1584 return data;
1585 }
1586
1587
1588 /*--------------------------------------------------------------------------*/
1589 /* write_macro() writes the current mesh (at the level of leaf elements) as */
1590 /* a macro triangulation to the specified file */
1591 /*--------------------------------------------------------------------------*/
1592
write_macro_master(MESH * mesh,const char * filename,macro_format format)1593 static int write_macro_master(MESH *mesh, const char *filename,
1594 macro_format format)
1595 {
1596 FUNCNAME("write_macro_master");
1597 int result = 0; /* make gcc happy */
1598 MACRO_DATA *data;
1599
1600 if (!filename)
1601 {
1602 ERROR("no filename specified, filename is NULL pointer\n");
1603 return(0);
1604 }
1605
1606 if (!mesh)
1607 {
1608 ERROR("no mesh specified, mesh is NULL pointer\n");
1609 return(0);
1610 }
1611
1612 if(!(data = mesh2macro_data(mesh))) {
1613 ERROR("Could not convert mesh to a macro data structure!\n");
1614 return(0);
1615 }
1616
1617 switch(format)
1618 {
1619 case ascii_format:
1620 result=write_macro_data(data, filename);
1621 break;
1622 case binary_format:
1623 result=write_macro_data_bin(data, filename);
1624 break;
1625 case xdr_format:
1626 result=write_macro_data_xdr(data, filename);
1627 }
1628
1629 free_macro_data(data);
1630
1631 return(result);
1632 }
1633
1634
1635 /*--------------------------------------------------------------------------*/
1636 /* These routines are available to the user: */
1637 /*--------------------------------------------------------------------------*/
1638
1639 /*--------------------------------------------------------------------------*/
1640 /* initialize and clear macro data structures */
1641 /*--------------------------------------------------------------------------*/
1642
alloc_macro_data(int dim,int nv,int ne)1643 MACRO_DATA *alloc_macro_data(int dim, int nv, int ne)
1644 {
1645 FUNCNAME("alloc_macro_data");
1646 MACRO_DATA *data = MEM_CALLOC(1, MACRO_DATA);
1647
1648 data->dim = dim;
1649 data->n_total_vertices = nv;
1650 data->n_macro_elements = ne;
1651
1652 data->coords = MEM_ALLOC(nv, REAL_D);
1653
1654 data->mel_vertices = MEM_ALLOC(ne*N_VERTICES(dim), int);
1655
1656 return data;
1657 }
1658
free_macro_data(MACRO_DATA * data)1659 void free_macro_data(MACRO_DATA *data)
1660 {
1661 int dim = data->dim;
1662 int ne = data->n_macro_elements;
1663 int nv = data->n_total_vertices;
1664
1665 MEM_FREE(data->coords, nv, REAL_D);
1666
1667 MEM_FREE(data->mel_vertices, ne*N_VERTICES(dim), int);
1668
1669 if (data->neigh) MEM_FREE(data->neigh, ne*N_NEIGH(dim), int);
1670 if (data->opp_vertex) MEM_FREE(data->opp_vertex, ne*N_NEIGH(dim), int);
1671 if (data->boundary) MEM_FREE(data->boundary, ne*N_NEIGH(dim), BNDRY_TYPE);
1672 #if DIM_MAX > 2
1673 if (dim == 3 && data->el_type) MEM_FREE(data->el_type, ne, U_CHAR);
1674 #endif
1675 if (data->wall_vtx_trafos) {
1676 MEM_FREE(data->wall_vtx_trafos,
1677 data->n_wall_vtx_trafos*sizeof(*data->wall_vtx_trafos),
1678 char);
1679 }
1680 if (data->el_wall_vtx_trafos) {
1681 MEM_FREE(data->el_wall_vtx_trafos, ne*N_WALLS(dim), int);
1682 }
1683 if (data->wall_trafos) {
1684 MEM_FREE(data->wall_trafos, data->n_wall_trafos, AFF_TRAFO);
1685 }
1686 if (data->el_wall_trafos) {
1687 MEM_FREE(data->el_wall_trafos, ne*N_WALLS(dim), int);
1688 }
1689
1690 #if ALBERTA_DEBUG
1691 if (data->mel_comment) {
1692 MEM_FREE(data->mel_comment[0], ne * 80, char);
1693 MEM_FREE(data->mel_comment, ne, char *);
1694 }
1695 #endif
1696
1697 MEM_FREE(data, 1, MACRO_DATA);
1698
1699 return;
1700 }
1701
read_macro(const char * filename)1702 MACRO_DATA *read_macro(const char *filename)
1703 {
1704 return read_macro_master(filename, ascii_format);
1705 }
1706
read_macro_bin(const char * filename)1707 MACRO_DATA *read_macro_bin(const char *filename)
1708 {
1709 return read_macro_master(filename, binary_format);
1710 }
1711
read_macro_xdr(const char * filename)1712 MACRO_DATA *read_macro_xdr(const char *filename)
1713 {
1714 return read_macro_master(filename, xdr_format);
1715 }
1716
1717
write_macro(MESH * mesh,const char * filename)1718 bool write_macro(MESH *mesh, const char *filename)
1719 {
1720 return(write_macro_master(mesh, filename, ascii_format));
1721 }
1722
write_macro_bin(MESH * mesh,const char * filename)1723 bool write_macro_bin(MESH *mesh, const char *filename)
1724 {
1725 return(write_macro_master(mesh, filename, binary_format));
1726 }
1727
write_macro_xdr(MESH * mesh,const char * filename)1728 bool write_macro_xdr(MESH *mesh, const char *filename)
1729 {
1730 return(write_macro_master(mesh, filename, xdr_format));
1731 }
1732
1733 /*--------------------------------------------------------------------------*/
1734 /* write raw macro triangulation in "data" to "filename" in standard ALBERTA*/
1735 /* key format */
1736 /*--------------------------------------------------------------------------*/
1737
write_macro_data(MACRO_DATA * data,const char * filename)1738 bool write_macro_data(MACRO_DATA *data, const char *filename)
1739 {
1740 FUNCNAME("write_macro_data");
1741 FILE *macro_file;
1742 int i, j, dim = data->dim;
1743
1744 if (!(macro_file = fopen(filename, "w")))
1745 {
1746 ERROR("could not open file %s for writing\n", filename);
1747 return false;
1748 }
1749
1750 fprintf(macro_file, "%s: %d\n", keys[KEY_DIM], dim);
1751 fprintf(macro_file, "%s: %d\n\n", keys[KEY_DOW], DIM_OF_WORLD);
1752
1753 fprintf(macro_file, "%s: %d\n", keys[KEY_NV], data->n_total_vertices);
1754 fprintf(macro_file, "%s: %d\n\n", keys[KEY_NEL], data->n_macro_elements);
1755
1756 fprintf(macro_file, "%s:\n", keys[KEY_V_COORD]);
1757 for(i = 0; i < data->n_total_vertices; i++)
1758 for (j = 0; j < DIM_OF_WORLD; j++)
1759 fprintf(macro_file, "%23.16e%s", data->coords[i][j],
1760 j < DIM_OF_WORLD-1 ? " " : "\n");
1761
1762 fprintf(macro_file, "\n%s:\n", keys[KEY_EL_VERT]);
1763 for(i = 0; i < data->n_macro_elements; i++) {
1764 for (j = 0; j < N_VERTICES(dim); j++) {
1765 fprintf(macro_file, " %5d", data->mel_vertices[VERT_IND(dim, i, j)]);
1766 }
1767 #if ALBERTA_DEBUG
1768 if (data->mel_comment) {
1769 fprintf(macro_file, " # %s", data->mel_comment[i]);
1770 }
1771 #endif
1772 fprintf(macro_file, "\n");
1773 }
1774
1775
1776 if(data->boundary) {
1777 fprintf(macro_file, "\n%s:\n", keys[KEY_EL_BND]);
1778 for(i = 0; i < data->n_macro_elements; i++)
1779 for (j = 0; j < N_NEIGH(dim); j++)
1780 fprintf(macro_file, "%4d%s", data->boundary[NEIGH_IND(dim, i, j)],
1781 j < N_NEIGH(dim)-1 ? " " : "\n");
1782 }
1783
1784 if(data->neigh) {
1785 fprintf(macro_file, "\n%s:\n", keys[KEY_EL_NEIGH]);
1786 for(i = 0; i < data->n_macro_elements; i++)
1787 for (j = 0; j < N_NEIGH(dim); j++)
1788 fprintf(macro_file, "%4d%s", data->neigh[NEIGH_IND(dim, i, j)],
1789 j < N_NEIGH(dim)-1 ? " " : "\n");
1790 }
1791
1792 #if DIM_MAX > 2
1793 if (dim == 3 && data->el_type)
1794 {
1795 fprintf(macro_file, "\n%s:\n", keys[KEY_EL_TYPE]);
1796 for(i = 0; i < data->n_macro_elements; i++)
1797 fprintf(macro_file, "%d%s", data->el_type[i], ((i+1)%20) ? " ": "\n");
1798 }
1799 #endif
1800
1801 if (data->n_wall_trafos) {
1802 int wt;
1803
1804 fprintf(macro_file, "\n%s: %d\n",
1805 keys[KEY_N_WALL_TRAFOS], data->n_wall_trafos);
1806 if (data->el_wall_trafos) {
1807 fprintf(macro_file, "\n%s:\n", keys[KEY_EL_WALL_TRAFOS]);
1808 for(i = 0; i < data->n_macro_elements; i++) {
1809 for (j = 0; j < N_WALLS(dim); j++) {
1810 fprintf(macro_file, "%4d%s",
1811 data->el_wall_trafos[NEIGH_IND(dim, i, j)],
1812 j < N_NEIGH(dim)-1 ? " " : "\n");
1813 }
1814 }
1815 }
1816 fprintf(macro_file, "\n%s:\n", keys[KEY_WALL_TRAFOS]);
1817 for (wt = 0; wt < data->n_wall_trafos; wt++) {
1818 fprintf(macro_file, "# wall transformation #%d\n", i);
1819 for (i = 0; i < DIM_OF_WORLD; i++) {
1820 for (j = 0; j < DIM_OF_WORLD; j++) {
1821 fprintf(macro_file, "%23.16e ", data->wall_trafos[wt].M[i][j]);
1822 }
1823 fprintf(macro_file, "%23.16e\n", data->wall_trafos[wt].t[i]);
1824 }
1825 fprintf(macro_file, "0 0 0 1\n");
1826 }
1827 }
1828
1829 if (data->n_wall_vtx_trafos) {
1830 fprintf(macro_file, "\n%s: %d\n",
1831 keys[KEY_N_WALL_VTX_TRAFOS], data->n_wall_vtx_trafos);
1832 fprintf(macro_file, "\n%s:\n", keys[KEY_WALL_VTX_TRAFOS]);
1833 for (i = 0; i < data->n_wall_vtx_trafos; i++) {
1834 fprintf(macro_file, "# wall vertex transformation #%d\n", i);
1835 for (j = 0; j < N_VERTICES(dim-1); j++) {
1836 fprintf(macro_file, "%4d %4d\n",
1837 data->wall_vtx_trafos[i][j][0], data->wall_vtx_trafos[i][j][1]);
1838 }
1839 }
1840 }
1841
1842 fprintf(macro_file, "\n");
1843
1844 fclose(macro_file);
1845
1846 INFO(2, 2, "wrote macro file %s\n", filename);
1847
1848 return true;
1849 }
1850
1851 /*--------------------------------------------------------------------------*/
1852 /* write raw macro triangulation in "data" to "filename" in native binary */
1853 /* format */
1854 /*--------------------------------------------------------------------------*/
1855
write_macro_data_bin(MACRO_DATA * data,const char * filename)1856 bool write_macro_data_bin(MACRO_DATA *data, const char *filename)
1857 {
1858 FUNCNAME("write_macro_data_bin");
1859 FILE *file;
1860 int i, dim = data->dim;
1861 char record_written=1;
1862 char record_not_written=0;
1863 int blah;
1864
1865 if(!data) {
1866 ERROR("no data - no file created\n");
1867 return 0;
1868 }
1869
1870 if (!(file = fopen(filename, "wb"))) {
1871 ERROR("cannot open file %s\n", filename);
1872 return 0;
1873 }
1874
1875 blah = fwrite(ALBERTA_VERSION, sizeof(char), strlen(ALBERTA_VERSION)+1, file);
1876
1877 i = sizeof(REAL);
1878 blah = fwrite(&i, sizeof(int), 1, file);
1879
1880 blah = fwrite(&dim, sizeof(int), 1, file);
1881
1882 i = DIM_OF_WORLD;
1883 blah = fwrite(&i, sizeof(int), 1, file);
1884
1885 blah = fwrite(&(data->n_total_vertices), sizeof(int), 1, file);
1886 blah = fwrite(&(data->n_macro_elements), sizeof(int), 1, file);
1887
1888 blah = fwrite(data->coords, sizeof(REAL_D), data->n_total_vertices, file);
1889 blah = fwrite(data->mel_vertices, sizeof(int),
1890 N_VERTICES(dim) * data->n_macro_elements, file);
1891
1892 if(data->boundary) {
1893 blah = fwrite(&record_written, sizeof(char), 1, file);
1894 blah = fwrite(data->boundary, sizeof(BNDRY_TYPE),
1895 N_NEIGH(dim) * data->n_macro_elements, file);
1896 } else {
1897 blah = fwrite(&record_not_written, sizeof(char), 1, file);
1898 }
1899
1900 if (data->neigh) {
1901 blah = fwrite(&record_written, sizeof(char), 1, file);
1902 blah = fwrite(data->neigh, sizeof(int),
1903 N_NEIGH(dim) * data->n_macro_elements, file);
1904 } else {
1905 blah = fwrite(&record_not_written, sizeof(char), 1, file);
1906 }
1907
1908 #if DIM_MAX > 2
1909 if (dim == 3 && data->el_type) {
1910 blah = fwrite(&record_written, sizeof(char), 1, file);
1911 blah = fwrite(data->el_type, sizeof(U_CHAR), data->n_macro_elements, file);
1912 } else {
1913 #endif
1914 blah = fwrite(&record_not_written, sizeof(char), 1, file);
1915 #if DIM_MAX > 2
1916 }
1917 #endif
1918
1919 blah = fwrite("EOF.", sizeof(char), 4, file);
1920 blah = fclose(file);
1921
1922 INFO(2, 2, "wrote macro binary-file %s\n", filename);
1923
1924 (void)blah;
1925
1926 return 1;
1927 }
1928
1929 /*--------------------------------------------------------------------------*/
1930 /* write raw macro triangulation in "data" to "filename" in xdr format */
1931 /*--------------------------------------------------------------------------*/
1932
write_macro_data_xdr(MACRO_DATA * data,const char * filename)1933 bool write_macro_data_xdr(MACRO_DATA *data, const char *filename)
1934 {
1935 FUNCNAME("write_macro_data_xdr");
1936 XDR *xdrp;
1937 int i, length;
1938 char *s;
1939 bool_t record_written=1;
1940 bool_t record_not_written=0;
1941
1942 caddr_t array_loc;
1943
1944 if(!data)
1945 {
1946 ERROR("no data - no file created\n");
1947 return(0);
1948 }
1949
1950 if (!(xdrp = xdr_open_file(filename, XDR_ENCODE)))
1951 {
1952 ERROR("cannot open file %s\n", filename);
1953 return(0);
1954 }
1955
1956 length = MAX(strlen(ALBERTA_VERSION) + 1, 5); /* length with terminating \0 */
1957 s=MEM_ALLOC(length, char);
1958 strcpy(s, ALBERTA_VERSION);
1959 xdr_string(xdrp, &s, length);
1960 MEM_FREE(s, length, char);
1961
1962 xdr_dim = data->dim;
1963
1964 xdr_int(xdrp, &xdr_dim);
1965
1966 i = DIM_OF_WORLD;
1967 xdr_int(xdrp, &i);
1968
1969 xdr_int(xdrp, &(data->n_total_vertices));
1970 xdr_int(xdrp, &(data->n_macro_elements));
1971
1972 array_loc=(caddr_t) data->coords;
1973 xdr_array(xdrp, &array_loc, (u_int *) &(data->n_total_vertices),
1974 (u_int) data->n_total_vertices, sizeof(REAL_D),
1975 (xdrproc_t) xdr_REAL_D);
1976
1977 array_loc=(caddr_t) data->mel_vertices;
1978 xdr_array(xdrp, &array_loc, (u_int *) &(data->n_macro_elements),
1979 (u_int) data->n_macro_elements * N_VERTICES(xdr_dim), sizeof(int),
1980 (xdrproc_t) xdr_int);
1981
1982 if(data->boundary) {
1983 xdr_bool(xdrp, &record_written);
1984 array_loc=(caddr_t) data->boundary;
1985 xdr_array(xdrp, &array_loc, (u_int *) &(data->n_macro_elements),
1986 (u_int) data->n_macro_elements * N_NEIGH(xdr_dim),
1987 sizeof(BNDRY_TYPE), (xdrproc_t) xdr_S_CHAR);
1988 }
1989 else xdr_bool(xdrp, &record_not_written);
1990
1991 if(data->neigh) {
1992 xdr_bool(xdrp, &record_written);
1993 array_loc=(caddr_t) data->neigh;
1994 xdr_array(xdrp, &array_loc, (u_int *) &(data->n_macro_elements),
1995 (u_int) data->n_macro_elements * N_NEIGH(xdr_dim), sizeof(int),
1996 (xdrproc_t) xdr_int);
1997 }
1998 else xdr_bool(xdrp, &record_not_written);
1999
2000 #if DIM_MAX > 2
2001 if (xdr_dim == 3 && data->el_type) {
2002 xdr_bool(xdrp, &record_written);
2003 array_loc=(caddr_t) data->el_type;
2004 xdr_array(xdrp, &array_loc, (u_int *) &(data->n_macro_elements),
2005 (u_int) data->n_macro_elements, sizeof(U_CHAR),
2006 (xdrproc_t) xdr_U_CHAR);
2007 }
2008 else
2009 #endif
2010 xdr_bool(xdrp, &record_not_written);
2011
2012 xdr_close_file(xdrp);
2013
2014 INFO(2, 2, "wrote macro xdr-file %s\n", filename);
2015
2016 return(1);
2017 }
2018
2019 static int find_matching_wall(MESH *mesh,
2020 int cur_mel, const AFF_TRAFO *wt, int wall,
2021 bool strict);
2022
2023 /* Transform the geomtric wall transformations specified by
2024 * init_wall_trafos() into topological wall transformations acting on
2025 * the vertex numbers of the elements. Check for consistency if the
2026 * macro triangulation specified combinatorical wall transformations
2027 * explicitly.
2028 *
2029 * This function is called _AFTER_ fill_neigh_info() and has to take
2030 * care of all additional neighbourhood relationships introduced by
2031 * periodic boundaries.
2032 *
2033 * This stuff looks a little bit ugly, but as the application has to
2034 * specifiy geometric wall-transformations anyway we can as well
2035 * derive the combinatorical transformations from the geometrical
2036 * ones.
2037 *
2038 * Return 0 if no face transformations apply, return 1 if at least one
2039 * is found and everything is consistent, return -1 if at least one is
2040 * found but the macro-triangulation is not compatible. If "strict"
2041 * bail out in this case.
2042 */
2043 static bool
init_wall_transformations(MESH * mesh,AFF_TRAFO * (* init_wall_trafos)(MESH * mesh,MACRO_EL * mel,int wall),bool strict)2044 init_wall_transformations(MESH *mesh,
2045 AFF_TRAFO *(*init_wall_trafos)(MESH *mesh,
2046 MACRO_EL *mel,
2047 int wall),
2048 bool strict)
2049 {
2050 FUNCNAME("init_wall_transformations");
2051 MACRO_EL *mel = mesh->macro_els;
2052 int dim = mesh->dim;
2053 int i, j, w;
2054 bool have_wall_trafos = false;
2055 AFF_TRAFO *cur_wt;
2056
2057 if (init_wall_trafos) {
2058 AFF_TRAFO **wall_trafos, *wt_neigh;
2059 int n_wall_trafos;
2060
2061 for (i = 0; i < mesh->n_macro_el; i++) {
2062 for (w = 0; w < N_WALLS(dim); w++) {
2063 cur_wt = mel[i].wall_trafo[w] = init_wall_trafos(mesh, mel + i, w);
2064 if (cur_wt) {
2065 if (mesh->wall_trafos == NULL) {
2066 mesh->wall_trafos =
2067 MEM_ALLOC(mesh->n_macro_el * N_WALLS(dim), AFF_TRAFO *);
2068 mesh->n_wall_trafos = 0;
2069 have_wall_trafos = true;
2070 }
2071 for (j = 0; j < mesh->n_wall_trafos; j++) {
2072 if (mesh->wall_trafos[j] == cur_wt) {
2073 break;
2074 }
2075 }
2076 if (j == mesh->n_wall_trafos) {
2077 ((AFF_TRAFO **)mesh->wall_trafos)[mesh->n_wall_trafos++] = cur_wt;
2078 }
2079 switch (find_matching_wall(mesh, i, cur_wt, w, strict)) {
2080 case 1:
2081 break; /* success */
2082 case -1: /* error */
2083 if (strict) {
2084 ERROR_EXIT(
2085 "Wall transformation identifies vertices in same element.\n");
2086 } else {
2087 WARNING(
2088 "Wall transformation identifies vertices in same element.\n");
2089 have_wall_trafos = false;
2090 }
2091 break;
2092 case 0: /* error */
2093 if (strict) {
2094 ERROR_EXIT(
2095 "Wall transformation does not seem to map walls to walls.\n");
2096 } else {
2097 WARNING(
2098 "Wall transformation does not seem to map walls to walls.\n");
2099 have_wall_trafos = false;
2100 }
2101 break;
2102 }
2103 }
2104 }
2105 }
2106
2107 TEST_EXIT(mesh->n_wall_trafos > 0, "No wall transformations apply???\n");
2108
2109 /* Make sure that the wall-transformations are grouped in pairs A,
2110 * A^{-1}
2111 */
2112
2113 /* The fancy allocation size is just to keep the memory allocation
2114 * book-keeping right; mesh_free() only frees mesh->wall_trafos.
2115 */
2116 wall_trafos = (AFF_TRAFO **)
2117 MEM_CALLOC(mesh->n_wall_trafos*(sizeof(AFF_TRAFO *) + sizeof(AFF_TRAFO)),
2118 char);
2119 n_wall_trafos = 0;
2120 if (have_wall_trafos) {
2121 /* Use the neighbourhood information */
2122 for (i = 0; i < mesh->n_macro_el; i++) {
2123 for (w = 0; w < N_WALLS(dim); w++) {
2124 if ((cur_wt = mel[i].wall_trafo[w]) != NULL) {
2125 for (j = 0; j < n_wall_trafos; j++) {
2126 if (wall_trafos[j] == cur_wt) {
2127 break;
2128 }
2129 }
2130 if (j == n_wall_trafos) {
2131 wt_neigh = mel[i].neigh[w]->wall_trafo[mel[i].opp_vertex[w]];
2132 wall_trafos[n_wall_trafos++] = cur_wt;
2133 wall_trafos[n_wall_trafos++] = wt_neigh;
2134 }
2135 }
2136 }
2137 }
2138 } else {
2139 /* Cannot use the neighbourhood information, invert the
2140 * transformation and do it the "hard way".
2141 */
2142 i = 0;
2143 do {
2144 cur_wt = mesh->wall_trafos[i];
2145 ((AFF_TRAFO **)mesh->wall_trafos)[i] = NULL;
2146 i++;
2147 for (j = 0; j < n_wall_trafos; j++) {
2148 if (wall_trafos[j] == cur_wt) {
2149 break;
2150 }
2151 }
2152 if (j == n_wall_trafos) {
2153 /* a new one, search for the inverse */
2154 AFF_TRAFO cur_inv, *inv_wt = NULL;
2155 REAL dist;
2156 int k;
2157
2158 INVAFF_DOW(cur_wt, &cur_inv);
2159 for (k = 0; k < mesh->n_wall_trafos; k++) {
2160 if (mesh->wall_trafos[k] == NULL) {
2161 continue;
2162 }
2163 inv_wt = mesh->wall_trafos[k];
2164 dist = MDST2_DOW((const REAL_D *)cur_inv.M,
2165 (const REAL_D *)inv_wt->M);
2166 dist += DST2_DOW(cur_inv.t, inv_wt->t);
2167 if (dist < 1e-24) {
2168 break;
2169 }
2170 }
2171 TEST_EXIT(k < mesh->n_wall_trafos,
2172 "Wall transformation without _distinct_ inverse. "
2173 "Involutions are not supported, sorry.\n");
2174 wall_trafos[n_wall_trafos++] = cur_wt;
2175 wall_trafos[n_wall_trafos++] = inv_wt;
2176 }
2177 } while (n_wall_trafos < mesh->n_wall_trafos);
2178 }
2179 /* strict or not strict: the following must hold. */
2180 TEST_EXIT(mesh->n_wall_trafos == n_wall_trafos,
2181 "Data inconsistency: mesh->n_wall_trafos = %d, counted %d.\n",
2182 mesh->n_wall_trafos, n_wall_trafos);
2183 MEM_FREE(mesh->wall_trafos, mesh->n_wall_trafos, AFF_TRAFO *);
2184 mesh->wall_trafos = wall_trafos;
2185 } else {
2186 have_wall_trafos = true;
2187 for (i = 0; i < mesh->n_macro_el; i++) {
2188 for (w = 0; w < N_WALLS(dim); w++) {
2189 if ((cur_wt = mel[i].wall_trafo[w]) != NULL) {
2190 switch (find_matching_wall(mesh, i, cur_wt, w, strict)) {
2191 case 1: /* success */
2192 break;
2193 case 0: /* error, not found */
2194 if (strict) {
2195 ERROR_EXIT(
2196 "Wall transformation does not seem to map walls to walls.\n");
2197 } else {
2198 WARNING(
2199 "Wall transformation does not seem to map walls to walls.\n");
2200 have_wall_trafos = false;
2201 }
2202 break;
2203 case -1: /* error, identification inside same elmenet */
2204 if (strict) {
2205 ERROR_EXIT(
2206 "Wall transformation identifies vertices in same element.\n");
2207 } else {
2208 WARNING(
2209 "Wall transformation identifies vertices in same element.\n");
2210 have_wall_trafos = false;
2211 }
2212 break;
2213 }
2214 } else {
2215 int wtno;
2216
2217 for (wtno = 0; wtno < mesh->n_wall_trafos; wtno++) {
2218 cur_wt = mesh->wall_trafos[wtno];
2219 switch (find_matching_wall(mesh, i, cur_wt, w, strict)) {
2220 case 1: /* success */
2221 mel[i].wall_trafo[w] = cur_wt;
2222 wtno = mesh->n_wall_trafos; /* break out of loop */
2223 break;
2224 case -1: /* error, identification inside same elmenet */
2225 if (strict) {
2226 ERROR_EXIT(
2227 "Wall transformation identifies vertices in same element.\n");
2228 } else {
2229 WARNING(
2230 "Wall transformation identifies vertices in same element.\n");
2231 have_wall_trafos = false;
2232 wtno = mesh->n_wall_trafos; /* break out of loop */
2233 }
2234 break;
2235 }
2236 }
2237 }
2238 }
2239 }
2240 }
2241
2242 return have_wall_trafos;
2243 }
2244
2245 /* Return:
2246 *
2247 * -1 if WT maps a vertex to another vertex in the same element (error case).
2248 * 0 if no matching wall in another element can be found.
2249 * 1 if a matching wall in another element can be found.
2250 */
find_matching_wall(MESH * mesh,int cur_mel,const AFF_TRAFO * wt,int wall,bool strict)2251 static int find_matching_wall(MESH *mesh,
2252 int cur_mel, const AFF_TRAFO *wt, int wall,
2253 bool strict)
2254 {
2255 FUNCNAME("find_matching_wall");
2256 MACRO_EL *mel = mesh->macro_els;
2257 int dim = mesh->dim;
2258 int j, k, l, wn, v;
2259 S_CHAR vmap[N_VERTICES(DIM_MAX-1)] = { 0, };
2260 REAL_D wall_img[N_VERTICES(DIM_MAX-1)] = { { 0.0, }, };
2261 REAL wall_vnrm[N_VERTICES(DIM_MAX-1)] = { 0.0, };
2262
2263 /* Compute the image of the wall-vertices under the action
2264 * of the _inverse_ wall transformation.
2265 */
2266 for (j = 0; j < N_VERTICES(dim-1); j++) {
2267 v = (wall + j + 1) % N_VERTICES(dim);
2268 AFFINV_DOW(wt, *mel[cur_mel].coord[v], wall_img[j]);
2269 wall_vnrm[j] = DIST_DOW(*mel[cur_mel].coord[wall], *mel[cur_mel].coord[v]);
2270 if (DIST_DOW(wall_img[j], *mel[cur_mel].coord[wall])
2271 <=
2272 1e-12*wall_vnrm[j]) {
2273 return -1; /* Maps to a vertex in the same element */
2274 }
2275 }
2276
2277 /* Now search for the matching wall */
2278 for (j = 0; j < mesh->n_macro_el; j++) {
2279 if (cur_mel == j) {
2280 continue; /* must not map to the same element */
2281 }
2282 for (wn = 0; wn < N_WALLS(dim); wn++) {
2283 if (mel[j].neigh[wn] && mel[j].neigh[wn] != mel+cur_mel) {
2284 continue;
2285 }
2286 for (l = 0; l < N_VERTICES(dim-1); l++) {
2287 vmap[l] = -1;
2288 }
2289 for (k = 0; k < N_VERTICES(dim-1); k++) {
2290 v = (wn + k + 1) % N_VERTICES(dim);
2291 for (l = 0; l < N_VERTICES(dim-1); l++) {
2292 if (vmap[l] == -1 &&
2293 DIST_DOW(wall_img[l], *mel[j].coord[v])
2294 <=
2295 1e-12*wall_vnrm[l]) { /* we have a match */
2296 vmap[l] = v;
2297 break;
2298 }
2299 }
2300 if (l == N_VERTICES(dim-1)) { /* does not match */
2301 break;
2302 }
2303 }
2304 if (k == N_VERTICES(dim-1)) { /* found it */
2305 break;
2306 }
2307 }
2308 if (wn < N_WALLS(dim)) { /* found a matching wall */
2309
2310 if (mel[cur_mel].neigh[wall] == NULL) {
2311 /* Add neighbourhood and opp-vertex specification. */
2312 mel[cur_mel].neigh[wall] = mel+j;
2313 mel[cur_mel].opp_vertex[wall] = wn;
2314 } else if (!(mel[cur_mel].neigh[wall] == mel+j &&
2315 mel[cur_mel].opp_vertex[wall] == wn)) {
2316 if (strict) {
2317 ERROR_EXIT("Wall transformation violates mesh connectivity.\n");
2318 } else {
2319 WARNING("Wall transformation violates mesh connectivity.\n");
2320 return false;
2321 }
2322 }
2323
2324 if (mel[j].neigh[wn] == NULL) {
2325 mel[j].neigh[wn] = mel+cur_mel;
2326 mel[j].opp_vertex[wn] = wall;
2327 } else if (!(mel[j].neigh[wn] == mel+cur_mel &&
2328 mel[j].opp_vertex[wn] == wall)) {
2329 if (strict) {
2330 ERROR_EXIT("Wall transformation violates mesh connectivity.\n");
2331 } else {
2332 WARNING("Wall transformation violates mesh connectivity.\n");
2333 return false;
2334 }
2335 }
2336
2337 if (mel[cur_mel].neigh_vertices[wall][0] != -1) {
2338 /* already have a periodic mapping for this wall, so
2339 * check for consistency
2340 */
2341 for (k = 0; k < N_VERTICES(dim-1); k++) {
2342 if (!(mel[cur_mel].neigh_vertices[wall][k] == vmap[k])) {
2343 if (strict) {
2344 ERROR_EXIT("Geometric and topological wall transformations "
2345 "do not match.\n");
2346 } else {
2347 WARNING("Geometric and topological wall transformations "
2348 "do not match.\n");
2349 return false;
2350 }
2351 }
2352 }
2353 } else {
2354 /* Fill the neigh_vertices information */
2355 for (k = 0; k < N_VERTICES(dim-1); k++) {
2356 mel[cur_mel].neigh_vertices[wall][k] = vmap[k];
2357 }
2358 }
2359 break;
2360 }
2361 }
2362
2363 /* Return "true" if we found a matching wall, otherwise "false". */
2364 return j < mesh->n_macro_el;
2365 }
2366
check_wall_transformations(MESH * mesh,bool strict)2367 static bool check_wall_transformations(MESH *mesh, bool strict)
2368 {
2369 FUNCNAME("check_wall_transformations");
2370 MACRO_EL *mel = mesh->macro_els;
2371 int dim = mesh->dim;
2372 int i, j, k, w, wn, v, nv, nnv, wnv;
2373 bool need_refine = false;
2374
2375 /* We do not allow wall transformations to map a vertex of a simplex
2376 * to another vertex on the same simplex; check this. Also check for
2377 * reflexivity.
2378 */
2379 for (i = 0; i < mesh->n_macro_el; i++) {
2380 for (w = 0; w < N_WALLS(dim); w++) {
2381 if (mel[i].neigh_vertices[w][0] != -1) {
2382 MACRO_EL *mel_n = mel[i].neigh[w];
2383 wn = mel[i].opp_vertex[w];
2384
2385 if ((mel_n = mel[i].neigh[w]) == NULL) {
2386 if (strict) {
2387 ERROR_EXIT("Wall transformation, but no neighour.\n");
2388 } else {
2389 WARNING("Wall transformation, but no neighour.\n");
2390 need_refine = true;
2391 continue;
2392 }
2393 }
2394 for (j = 0; j < N_VERTICES(dim-1); j++ ) {
2395 REAL_D *ncoord;
2396
2397 v = (w + j + 1) % N_VERTICES(dim);
2398 nv = mel[i].neigh_vertices[w][j];
2399 if (nv < wn) {
2400 wnv = nv + N_VERTICES(dim) - wn - 1;
2401 } else {
2402 wnv = nv - wn - 1;
2403 }
2404 nnv = mel_n->neigh_vertices[wn][wnv];
2405 TEST_EXIT(v == nnv,
2406 "Wall transformations are not inverse to each other.\n");
2407
2408 ncoord = mel_n->coord[nv];
2409 for (k = 0; k < N_VERTICES(dim); k++) {
2410 if (!(ncoord != mel[i].coord[k])) {
2411 if (strict) {
2412 ERROR_EXIT("Vertices must not be mapped to vertices of "
2413 "the same element.\n");
2414 } else {
2415 WARNING("Vertices must not be mapped to vertices of "
2416 "the same element.\n");
2417 need_refine = true;
2418 }
2419 }
2420 }
2421 }
2422 }
2423 }
2424 }
2425
2426 return !need_refine;
2427 }
2428
2429 /* Try to refine the mesh s.t. that the
2430 * no-vertex-maps-to-the-same-element property is fulfilled. We do so
2431 * by using global refinement ATM. More sophisticated algorithms could
2432 * be implemented.
2433 *
2434 * We only support the case when the periodic structure in induced by
2435 * geometric face-transformations, simply because the case of pure
2436 * combinatorical face-transformations is more complicated and I'm not
2437 * in the mood to implement that right now.
2438 */
2439 static void
try_resolve_periodic_walls(MESH * mesh,const MACRO_DATA * data,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int),AFF_TRAFO * (* init_wall_trafos)(MESH *,MACRO_EL *,int wall))2440 try_resolve_periodic_walls(MESH *mesh, const MACRO_DATA *data,
2441 NODE_PROJECTION *(*init_node_proj)(
2442 MESH *, MACRO_EL *, int),
2443 AFF_TRAFO *(*init_wall_trafos)(
2444 MESH *, MACRO_EL *, int wall))
2445 {
2446 FUNCNAME("try_resolve_periodic_walls");
2447 MACRO_DATA np_data[1]; /* not-periodic clone of "data" */
2448 MACRO_DATA *p_data; /* macro-data generated from globally refined mesh */
2449 int dim = mesh->dim;
2450 MESH *np_mesh, *p_mesh, swap;
2451 int i, j, el, nwt;
2452
2453 *np_data = *data;
2454 np_data->n_wall_trafos = 0;
2455 np_data->wall_trafos = NULL;
2456 np_data->el_wall_trafos = NULL;
2457 np_data->n_wall_vtx_trafos = 0;
2458 np_data->wall_vtx_trafos = NULL;
2459 np_data->el_wall_vtx_trafos = NULL;
2460
2461 /* Generate a new mesh and refine it globally s.t. each edge is
2462 * bisected at least once. If the issue can be resolved at all by
2463 * refinement, then this will hack it.
2464 */
2465 np_mesh =
2466 _AI_get_mesh(dim, "temporary periodic mesh", np_data, init_node_proj, NULL,
2467 true /* strict_periodic */);
2468 global_refine(np_mesh, np_mesh->dim, FILL_NOTHING);
2469 p_data = mesh2macro_data(np_mesh);
2470
2471 p_data->wall_trafos = MEM_ALLOC(mesh->n_wall_trafos/2, AFF_TRAFO);
2472 p_data->n_wall_trafos = mesh->n_wall_trafos/2;
2473 p_data->el_wall_trafos =
2474 MEM_CALLOC(p_data->n_macro_elements*N_WALLS(dim), int);
2475 for (i = 0; i < p_data->n_wall_trafos; i++) {
2476 p_data->wall_trafos[i] = *mesh->wall_trafos[2*i];
2477 }
2478
2479 /* Now we have to generate the missing periodic information. We
2480 * exploit the fact that mesh2macro_data() enumerates the macro
2481 * elements in the order of an ordinary leaf-element traversal.
2482 *
2483 * We also exploit the fact that the macro-elements of mesh and
2484 * np_mesh occur in the same order.
2485 *
2486 * The periodic mesh was carrying geometric wall-transformations,
2487 * and these are still available, we simply have to collect them.
2488 */
2489 el = 0;
2490 TRAVERSE_FIRST(np_mesh, -1, CALL_LEAF_EL|FILL_MACRO_WALLS) {
2491 MACRO_EL *p_mel;
2492 int mi;
2493 AFF_TRAFO *wt;
2494
2495 mi = el_info->macro_el->index;
2496 p_mel = &mesh->macro_els[mi];
2497
2498 for (i = 0; i < N_WALLS(dim); i++) {
2499 int mwall = el_info->macro_wall[i];
2500 if (mwall < 0) {
2501 continue;
2502 }
2503 if ((wt = p_mel->wall_trafo[mwall]) == NULL) {
2504 continue;
2505 }
2506 for (j = 0;
2507 j < mesh->n_wall_trafos && mesh->wall_trafos[j] != wt;
2508 j++);
2509 DEBUG_TEST_EXIT(j < mesh->n_wall_trafos,
2510 "Wall transformation not found.");
2511 if (j % 2 != 0) {
2512 j = -j/2-1;
2513 } else {
2514 j = j/2+1;
2515 }
2516 p_data->el_wall_trafos[NEIGH_IND(dim, el, i)] = j;
2517 }
2518
2519 el++;
2520
2521 } TRAVERSE_NEXT();
2522
2523 p_mesh = _AI_get_mesh(dim, mesh->name, p_data, NULL, NULL, true /* strict */);
2524
2525 /* It remains to copy over the node-projections from the old mesh,
2526 * and to adjust the pointers to the face transformations in case
2527 * the application came along with its own init_wall_trafos() hook.
2528 */
2529 el = 0;
2530 TRAVERSE_FIRST(np_mesh, -1, CALL_LEAF_EL|FILL_MACRO_WALLS) {
2531
2532 p_mesh->macro_els[el].projection[0] = el_info->macro_el->projection[0];
2533
2534 for (i = 0; i < N_WALLS(dim); i++) {
2535 int mwall = el_info->macro_wall[i];
2536 if (mwall < 0) {
2537 continue;
2538 }
2539 p_mesh->macro_els[el].projection[i+1] =
2540 el_info->macro_el->projection[mwall];
2541 }
2542
2543 el++;
2544 } TRAVERSE_NEXT();
2545
2546 if (init_wall_trafos) {
2547 for (el = 0; el < p_mesh->n_macro_el; el++) {
2548 MACRO_EL *mel = &p_mesh->macro_els[el];
2549
2550 for (i = 0; i < N_WALLS(dim); i++) {
2551 if ((nwt = p_data->el_wall_trafos[NEIGH_IND(dim, el, i)]) > 0) {
2552 nwt--;
2553 mel->wall_trafo[i] = mesh->wall_trafos[2*nwt];
2554 } else if (nwt < 0) {
2555 nwt = - nwt - 1;
2556 mel->wall_trafo[i] = mesh->wall_trafos[2*nwt+1];
2557 }
2558 }
2559 }
2560 memcpy((AFF_TRAFO **)p_mesh->wall_trafos, mesh->wall_trafos,
2561 mesh->n_wall_trafos * sizeof(AFF_TRAFO *));
2562 }
2563
2564 free_mesh(np_mesh);
2565 free_macro_data(p_data);
2566
2567 swap = *mesh;
2568 *mesh = *p_mesh;
2569 *p_mesh = swap;
2570
2571 free_mesh(p_mesh);
2572 }
2573
2574 /* Compute the boundary classification for lower-dimensional
2575 * sub-simplices (vertices and edges (3d)) on the macro level. If
2576 * COUNT == TRUE, then also update the vertex, edge, face numbers in
2577 * MESH, otherwise just fill in the boundary information.
2578 */
_AI_fill_bound_info(MESH * mesh,int * mel_vertices,int nv,int ne,bool count)2579 void _AI_fill_bound_info(MESH *mesh,
2580 int *mel_vertices, int nv, int ne, bool count)
2581 {
2582 FUNCNAME("_AI_fill_bound_info");
2583
2584 switch (mesh->dim) {
2585 case 1:
2586 fill_bound_info_1d(mesh, mel_vertices, nv, ne);
2587 /*fill_mel_orientation_1d(mesh->macro_els, mesh->n_macro_el);*/
2588 break;
2589 #if DIM_MAX > 1
2590 case 2:
2591 fill_bound_info_2d(mesh, mel_vertices, nv, ne);
2592 if (count) {
2593 count_edges_2d(mesh);
2594 }
2595 /*fill_mel_orientation_2d(mesh->macro_els, mesh->n_macro_el);*/
2596 break;
2597 #if DIM_MAX > 2
2598 case 3:
2599 fill_bound_info_3d(mesh, mel_vertices, nv, ne);
2600 fill_more_bound_info_3d(mesh, mel_vertices, ne, count);
2601 /*fill_mel_orientation_3d(mesh->macro_els, mesh->n_macro_el);*/
2602 break;
2603 #endif
2604 #endif
2605 default:
2606 ERROR_EXIT("Illegal dimension %d!\n", mesh->dim);
2607 }
2608 }
2609
init_node_projections(MESH * mesh,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int))2610 static void init_node_projections(MESH *mesh,
2611 NODE_PROJECTION *(*init_node_proj)(
2612 MESH *, MACRO_EL *, int))
2613 {
2614 int i;
2615 #if DIM_MAX > 1
2616 int j, dim = mesh->dim;
2617 #endif
2618 MACRO_EL *mel = mesh->macro_els;
2619
2620 /****************************************************************************
2621 * Call the user-defined new vertex projection assignment routine
2622 * "init_node_proj".
2623 ***************************************************************************/
2624 if (init_node_proj) {
2625 for (i = 0; i < mesh->n_macro_el; i++) {
2626 mel[i].projection[0] = init_node_proj(mesh, mel + i, 0);
2627
2628 #if DIM_MAX > 1
2629 if(dim == 2)
2630 for(j = 1; j < N_NEIGH_2D + 1; j++)
2631 mel[i].projection[j] = init_node_proj(mesh, mel + i, j);
2632 #if DIM_MAX > 2
2633 else if(dim == 3)
2634 for(j = 1; j < N_NEIGH_3D + 1; j++)
2635 mel[i].projection[j] = init_node_proj(mesh, mel + i, j);
2636 #endif
2637 #endif
2638
2639 #define REALLY_DONT_USE_THIS_CODE true
2640 #if !REALLY_DONT_USE_THIS_CODE
2641 /* No point in doing this at the moment, DK. */
2642 /* cH: I think we really should do this. Shouldn't we? FIXME! */
2643 /* cH: Update: there is no point in doing this _at_ _all_, the
2644 * code below breaks periodic meshes.
2645 */
2646 /****************************************************************************/
2647 /* If necessary, copy projections to neighbour elements. */
2648 /****************************************************************************/
2649 for (j = 0; j < N_NEIGH(dim); j++) {
2650 NODE_PROJECTION *tmp_proj;
2651 MACRO_EL *neigh;
2652
2653 if ((neigh = mel[i].neigh[j]) != NULL) {
2654 int k;
2655
2656 tmp_proj = mel[i].projection[0];
2657
2658 if (mel[i].projection[j+1])
2659 tmp_proj = mel[i].projection[j+1];
2660
2661 /* Search for the correct subsimplex on the neighbour element. */
2662 for (k = 0; k < N_NEIGH(dim); k++)
2663 if (neigh->neigh[k] == mel + i) break;
2664
2665 neigh->projection[k+1] = tmp_proj;
2666 }
2667 }
2668 #endif
2669 }
2670 }
2671 }
2672
2673 /***************************************************************************/
2674 /* macro_data2mesh(): */
2675 /* copy macro data to the MESH structure "mesh" provided: */
2676 /* 1) set most entries in "mesh" */
2677 /* 2) allocate macro elements and link them to "mesh" */
2678 /* 3) calculate macro element orientation for 3D */
2679 /* 4) calculate the mesh size for "mesh->bbox" */
2680 /* 5) Initialize slave meshes */
2681 /* */
2682 /* the entire MACRO_DATA structure can be freed after use! */
2683 /***************************************************************************/
2684
2685 void
macro_data2mesh(MESH * mesh,const MACRO_DATA * data,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int),AFF_TRAFO * (* init_wall_trafos)(MESH *,MACRO_EL *,int wall))2686 macro_data2mesh(MESH *mesh, const MACRO_DATA *data,
2687 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
2688 AFF_TRAFO *(*init_wall_trafos)(MESH *, MACRO_EL *, int wall))
2689 {
2690 _AI_macro_data2mesh(mesh, data, init_node_proj, init_wall_trafos, false);
2691 }
2692
2693 void
_AI_macro_data2mesh(MESH * mesh,const MACRO_DATA * data,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int),AFF_TRAFO * (* init_wall_trafos)(MESH *,MACRO_EL *,int wall),bool strict_periodic)2694 _AI_macro_data2mesh(MESH *mesh, const MACRO_DATA *data,
2695 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
2696 AFF_TRAFO *(*init_wall_trafos)(
2697 MESH *, MACRO_EL *, int wall),
2698 bool strict_periodic)
2699 {
2700 FUNCNAME("macro_data2mesh");
2701 int i, j, dim = data->dim;
2702 MACRO_EL *mel;
2703 REAL_D *newcoords;
2704
2705 TEST_EXIT(mesh, "no mesh, mesh is NULL pointer!\n");
2706
2707 mesh->dim = dim;
2708
2709 mesh->n_elements =
2710 mesh->n_hier_elements =
2711 mesh->n_macro_el = data->n_macro_elements;
2712 mesh->n_vertices = data->n_total_vertices;
2713
2714 mel = mesh->macro_els = MEM_CALLOC(data->n_macro_elements, MACRO_EL);
2715
2716 newcoords = MEM_ALLOC(data->n_total_vertices, REAL_D);
2717
2718 for(i = 0; i < data->n_total_vertices; i++)
2719 COPY_DOW(data->coords[i], newcoords[i]);
2720
2721 ((MESH_MEM_INFO *)mesh->mem_info)->count = data->n_total_vertices;
2722 ((MESH_MEM_INFO *)mesh->mem_info)->coords = newcoords;
2723
2724 for(i = 0; i < data->n_macro_elements; i++) {
2725 mel[i].el = get_element(mesh);
2726
2727 mel[i].index = i;
2728 #if ALBERTA_DEBUG
2729 mel[i].el->index = i;
2730 #endif
2731
2732 for(j = 0; j < N_VERTICES(dim); j++)
2733 mel[i].coord[j] = &newcoords[data->mel_vertices[VERT_IND(dim, i, j)]];
2734
2735 #if DIM_MAX >= 3
2736 if(dim == 3) {
2737 mel[i].el_type = data->el_type ? data->el_type[i] : 0;
2738 mel[i].orientation = AI_get_orientation_3d(mel + i);
2739 }
2740 #endif
2741 }
2742
2743 if (mesh->parametric)
2744 WARNING("mesh->bbox may not be computed correctly, "
2745 "problems with graphical output may occur\n");
2746 /* else */ /* should still be better than nothing */
2747 calculate_size(mesh, data);
2748
2749 if(dim > 0) {
2750 TEST_EXIT(data->neigh != NULL,
2751 "Neighbour information must be present!\n");
2752 TEST_EXIT(data->boundary != NULL,
2753 "Boundary information must be present!\n");
2754
2755 fill_neigh_info(mel, data);
2756
2757 /* Initialize periodic boundaries before filling the boundary
2758 * information. However, we need basic information about
2759 * wall-boundaries because init_wall_trafos() will probably need
2760 * such information to do its job. Therefore we fill the
2761 * wall-boundary information (without checking) at this point.
2762 *
2763 * We also compute the number of vertex orbits (per_n_vertices)
2764 * here.
2765 *
2766 * Note: the application may choose to override the
2767 * wall-transformations specified in the macro triangulation;
2768 * init_wall_trafos() takes precedence, at any case.
2769 *
2770 * For simplicity we store the wall-transformation _AND_ their
2771 * inverses.
2772 */
2773 if (data->n_wall_trafos && init_wall_trafos == NULL) {
2774 int wt;
2775
2776 mesh->n_wall_trafos = 2*data->n_wall_trafos;
2777 mesh->wall_trafos = MEM_ALLOC(mesh->n_wall_trafos, AFF_TRAFO *);
2778 ((AFF_TRAFO **)mesh->wall_trafos)[0] =
2779 MEM_ALLOC(mesh->n_wall_trafos, AFF_TRAFO);
2780 for (wt = 0; wt < data->n_wall_trafos; wt++) {
2781 ((AFF_TRAFO **)mesh->wall_trafos)[2*wt] = mesh->wall_trafos[0]+2*wt;
2782 *mesh->wall_trafos[2*wt] = data->wall_trafos[wt];
2783 ((AFF_TRAFO **)mesh->wall_trafos)[2*wt+1] = mesh->wall_trafos[0]+2*wt+1;
2784 INVAFF_DOW(&data->wall_trafos[wt], mesh->wall_trafos[2*wt+1]);
2785 }
2786 mesh->is_periodic = true;
2787 if (data->el_wall_trafos) {
2788 /* Even the mapping to the macro element walls was
2789 * specified. Use it. In this case the mapping of faces to
2790 * faces must be consistent and we will _NOT_ attempt to
2791 * resolve conflicts by global refinement.
2792 */
2793 for (i = 0; i < mesh->n_macro_el; i++) {
2794 for (j = 0; j < N_WALLS(dim); j++) {
2795 int wt;
2796 if ((wt = data->el_wall_trafos[NEIGH_IND(dim, i, j)]) > 0) {
2797 mel[i].wall_trafo[j] = mesh->wall_trafos[2*(wt - 1)];
2798 } else if (wt < 0) {
2799 mel[i].wall_trafo[j] = mesh->wall_trafos[-2*(wt + 1) + 1];
2800 }
2801 }
2802 }
2803 }
2804 }
2805 if (init_wall_trafos || mesh->n_wall_trafos > 0) {
2806 int (*wall_vtx_trafos)[N_VERTICES(DIM_MAX-1)][2];
2807 int nwt;
2808
2809 for (i = 0; i < mesh->n_macro_el; i++) {
2810 for (j = 0; j < N_NEIGH(dim); j++) {
2811 mel[i].wall_bound[j] = data->boundary[NEIGH_IND(dim, i, j)];
2812 }
2813 }
2814 mesh->is_periodic = true;
2815 if (!init_wall_transformations(mesh, init_wall_trafos, strict_periodic) ||
2816 !check_wall_transformations(mesh, strict_periodic)) {
2817 if (strict_periodic) {
2818 ERROR_EXIT("No non-trivial wall-transformation, or "
2819 "incompatible macro-triangulation.\n");
2820 } else {
2821 WARNING("Trying to resolve periodic "
2822 "boundaries by global refinement.\n");
2823 }
2824 try_resolve_periodic_walls(mesh, data,
2825 init_node_proj, init_wall_trafos);
2826 return;
2827 }
2828
2829 nwt = _AI_compute_macro_wall_trafos(mesh, &wall_vtx_trafos);
2830 mesh->per_n_vertices = mesh->n_vertices;
2831 (void)_AI_wall_trafo_vertex_orbits(dim, wall_vtx_trafos, nwt,
2832 NULL, &mesh->per_n_vertices);
2833 MEM_FREE(wall_vtx_trafos, nwt*N_VERTICES(DIM_MAX-1)*2, int);
2834 } else {
2835 mesh->is_periodic = data->n_wall_vtx_trafos > 0;
2836 if (mesh->is_periodic) {
2837 mesh->per_n_vertices = mesh->n_vertices;
2838 (void)_AI_wall_trafo_vertex_orbits(dim, data->wall_vtx_trafos,
2839 data->n_wall_vtx_trafos,
2840 0, &mesh->per_n_vertices);
2841 }
2842 }
2843
2844 for(i = 0; i < data->n_macro_elements; i++) {
2845 for(j = 0; j < N_WALLS(dim); j++) {
2846 mesh->macro_els[i].wall_bound[j] = data->boundary[NEIGH_IND(dim, i, j)];
2847 }
2848 }
2849
2850 _AI_fill_bound_info(mesh,
2851 data->mel_vertices, mesh->n_vertices, mesh->n_elements,
2852 true);
2853
2854 }
2855
2856 init_node_projections(mesh, init_node_proj);
2857
2858 }
2859
macro_test(MACRO_DATA * data,const char * new_name)2860 void macro_test(MACRO_DATA *data, const char *new_name)
2861 {
2862 FUNCNAME("macro_test");
2863
2864 switch(data->dim) {
2865 case 0:
2866 break;
2867 case 1:
2868 macro_test_1d(data, new_name);
2869 break;
2870 #if DIM_MAX > 1
2871 case 2:
2872 macro_test_2d(data, new_name);
2873 break;
2874 #endif
2875 #if DIM_MAX > 2
2876 case 3:
2877 macro_test_3d(data, new_name);
2878 break;
2879 #endif
2880 default:
2881 ERROR_EXIT("Illegal dim == %d!\n", data->dim);
2882 }
2883
2884 /* We do not allow mapping a vertex of a simplex to another vertex
2885 * on the same element. Otherwise we would have to cope with
2886 * pathological cases which can occur only with the unrefined macro
2887 * triangulation (almost only ...). The benefit is that we can use
2888 * the global numbering of the vertices (by means of a vertex
2889 * DOF_ADMIN) to determine things like orientation etc.
2890 */
2891 if (data->n_wall_vtx_trafos) {
2892 int el, wt, v, w, wv;
2893 int dim = data->dim;
2894
2895 for (el = 0; el < data->n_macro_elements; ++el) {
2896 for (w = 0; w < N_WALLS(dim); w++) {
2897 if ((wt = data->el_wall_vtx_trafos[N_WALLS(dim)*el+w]) != 0) {
2898 int img;
2899 if (wt > 0) {
2900 --wt;
2901 img = 1;
2902 } else {
2903 wt = -wt-1;
2904 img = 0;
2905 }
2906 for (wv = 0; wv < N_VERTICES(dim-1); wv++) {
2907 int vertex_image = data->wall_vtx_trafos[wt][wv][img];
2908 for (v = 0; v < N_VERTICES(dim); v++) {
2909 TEST_EXIT(vertex_image
2910 != data->mel_vertices[N_VERTICES(dim)*el+v],
2911 "ERROR: Unsupported feature in the context of periodic "
2912 "meshes: The walls of elements may not be mapped onto "
2913 "another wall on the same element; you have to refine "
2914 "your macro triangulation. Element nr: %d, "
2915 "wall trafo: %d, vertex (src/dst): %d/%d\n",
2916 el, wt,
2917 data->wall_vtx_trafos[wt][wv][1-img],
2918 vertex_image);
2919 }
2920 }
2921 }
2922 }
2923 }
2924 }
2925 }
2926