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 /* file: read_mesh.c */
7 /* */
8 /* description: reading data of mesh and vectors in machine independent */
9 /* and native binary formats. */
10 /* */
11 /*--------------------------------------------------------------------------*/
12 /* */
13 /* authors: Alfred Schmidt */
14 /* Zentrum fuer Technomathematik */
15 /* Fachbereich 3 Mathematik/Informatik */
16 /* Univesitaet Bremen */
17 /* Bibliothekstr. 2 */
18 /* D-28359 Bremen, Germany */
19 /* */
20 /* Kunibert G. Siebert */
21 /* Institut fuer Mathematik */
22 /* Universitaet Augsburg */
23 /* Universitaetsstr. 14 */
24 /* D-86159 Augsburg, Germany */
25 /* */
26 /* http://www.mathematik.uni-freiburg.de/IAM/ALBERTA */
27 /* */
28 /* (c) by A. Schmidt and K.G. Siebert (1996-2003) */
29 /* */
30 /*--------------------------------------------------------------------------*/
31 /*--------------------------------------------------------------------------*/
32 /* to_do : number of dofs depending on node */
33 /* nach read_macro, damit dort pointer auf mel, v vorhanden sind!!! */
34 /* (fuers naechste Schreiben) */
35 /* error handling is terrible */
36 /*--------------------------------------------------------------------------*/
37
38 #include "alberta_intern.h"
39 #include "alberta.h"
40
41 #include <string.h>
42
43 static XDR *xdrp;
44 static FILE *file;
45
46 /*
47 WARNING:
48 XDR routines to read/write ALBERTA types must be changed if the
49 ALBERTA types change!
50
51 current state: REAL = double
52 U_CHAR = unsigned char
53 S_CHAR = signed char
54 DOF = int
55
56 Another WARNING! (D.K.)
57 XDR routines are not well documented in the "xdr" man page.
58 Do not change anything unless you know what you are doing!
59 */
60
61 typedef DOF_REAL_VEC DOF_VEC;
62
63 typedef enum dv_type_enum {
64 __DOF_REAL_VEC = 0,
65 __DOF_REAL_D_VEC,
66 __DOF_REAL_VEC_D,
67 __DOF_INT_VEC,
68 __DOF_SCHAR_VEC,
69 __DOF_UCHAR_VEC,
70 } DV_TYPE_ENUM;
71
72 static DOF_VEC *read_dof_vec_master(DV_TYPE_ENUM dv_type,
73 DOF_VEC *dv,
74 MESH *mesh, FE_SPACE *fe_space,
75 bool expect_next);
76
AI_fread(const void * ptr,size_t size,size_t nmemb,FILE * stream)77 size_t AI_fread(const void *ptr, size_t size, size_t nmemb, FILE *stream)
78 {
79 size_t nread;
80
81 nread = fwrite(ptr, size, nmemb, stream);
82
83 /* Now, what to do with RESULT? ;) */
84 return nread;
85 }
86
AI_xdr_int(XDR * xdr,void * ip)87 bool_t AI_xdr_int(XDR *xdr, void *ip)
88 {
89 FUNCNAME("AI_xdr_int");
90
91 TEST_EXIT(sizeof(int) == 4,
92 "sizeof(int) != 4, please start hacking!");
93 #if HAVE_XDR_INT32_T
94 return (xdr_int32_t(xdr,(int *)ip));
95 #elif HAVE_XDR_INT
96 return (xdr_int(xdr,(int *)ip));
97 #endif
98 }
99
AI_xdr_REAL(XDR * xdr,void * rp)100 bool_t AI_xdr_REAL(XDR *xdr, void *rp)
101 {
102 FUNCNAME("AI_xdr_REAL");
103
104 TEST_EXIT(sizeof(double) == sizeof(REAL),
105 "sizeof(double) != sizeof(REAL), please start hacking!");
106 return (xdr_double(xdr,(double *)rp));
107 }
108
AI_xdr_U_CHAR(XDR * xdr,void * ucp)109 bool_t AI_xdr_U_CHAR(XDR *xdr, void *ucp)
110 {
111 return (xdr_u_char(xdr,(u_char *)ucp));
112 }
113
AI_xdr_S_CHAR(XDR * xdr,void * cp)114 bool_t AI_xdr_S_CHAR(XDR *xdr, void *cp)
115 {
116 return (xdr_char(xdr,(char *)cp));
117 }
118
AI_xdr_DOF(XDR * xdr,void * dp)119 bool_t AI_xdr_DOF(XDR *xdr, void *dp)
120 {
121 FUNCNAME("AI_xdr_DOF");
122
123 TEST_EXIT(sizeof(DOF) == sizeof(int),
124 "sizeof(DOF) != sizeof(int), please start hacking!");
125 return (AI_xdr_int(xdr,(int *)dp));
126 }
127
128 #if 1
129 /* Needed to read old mesh-files. */
read_xdr_file(void * file,void * buffer,unsigned int size)130 static int read_xdr_file(void *file, void *buffer, unsigned int size)
131 {
132 return (int)fread(buffer, 1, (size_t)size, (FILE *)file);
133 }
134
write_xdr_file(void * file,void * buffer,unsigned int size)135 static int write_xdr_file(void *file, void *buffer, unsigned int size)
136 {
137 return fwrite(buffer, (size_t)size, 1, (FILE *)file) == 1 ? (int)size : 0;
138 }
139 #endif
140
AI_xdr_fopen(FILE * fp,enum xdr_op mode)141 XDR *AI_xdr_fopen(FILE *fp, enum xdr_op mode)
142 {
143 FUNCNAME("AI_xdr_open_file");
144 XDR *xdr;
145
146 if (!(xdr = MEM_ALLOC(1,XDR))) {
147 ERROR("can't allocate memory for xdr pointer.\n");
148 return NULL;
149 }
150
151 xdrstdio_create(xdr, file = fp, mode);
152 #if 0
153 /* only needed for old meshes, this function is now integrated into
154 * the version check
155 */
156 xdrrec_create(xdr, 65536, 65536, (caddr_t) file,
157 (int (*)())read_xdr_file, (int (*)())write_xdr_file);
158
159 xdr->x_op = mode;
160 xdr->x_public = (caddr_t)file;
161
162 if (mode == XDR_DECODE) {
163 xdrrec_skiprecord(xdr);
164 }
165 #endif
166
167 return xdr;
168 }
169
AI_xdr_open_file(const char * fn,enum xdr_op mode)170 XDR *AI_xdr_open_file(const char *fn, enum xdr_op mode)
171 {
172 if ((file = fopen(fn, (mode == XDR_DECODE) ? "r": "w"))) {
173 return AI_xdr_fopen(file, mode);
174 }
175 return NULL;
176 }
177
AI_xdr_close(XDR * xdr)178 bool AI_xdr_close(XDR *xdr)
179 {
180 FUNCNAME("AI_xdr_close");
181 if (!xdr) {
182 ERROR("NULL xdr pointer.\n");
183 return false;
184 }
185
186 //TODO: still needed? see quick hack concerning the retrying of
187 //reading the alberta version
188 #if 0
189 if (xdr->x_op == XDR_ENCODE) {
190 xdrrec_endofrecord(xdr, 1);
191 }
192 #endif
193
194 xdr_destroy(xdr);
195
196 MEM_FREE(xdr,1,XDR);
197
198 return true;
199 }
200
AI_xdr_close_file(XDR * xdr)201 bool AI_xdr_close_file(XDR *xdr)
202 {
203 FUNCNAME("AI_xdr_close_file");
204
205 if (!AI_xdr_close(xdr)) {
206 return false;
207 }
208
209 if (fclose(file)) {
210 ERROR("error closing file.\n");
211 }
212 return true;
213 }
214
215 /*--------------------------------------------------------------------------*/
216
217 /* void _AI_read_REAL() */
_AI_read_REAL(REAL * val)218 void _AI_read_REAL(REAL *val)
219 {
220 bool result;
221
222 if (xdrp)
223 result = AI_xdr_REAL(xdrp, val);
224 else
225 result = fread(val, sizeof(REAL), 1, file) == 1;
226 /* Now, what to do with RESULT? ;) */
227 (void)result;
228 }
229
_AI_read_int(int * val)230 void _AI_read_int(int *val)
231 {
232 bool result;
233
234 if (xdrp)
235 result = AI_xdr_int(xdrp, val);
236 else
237 result = fread(val, sizeof(int), 1, file) == 1;
238
239 /* Now, what to do with RESULT? ;) */
240 (void)result;
241 }
242
243 #if 0
244 /* NEVER use this. If you need 64 bits, then use xdr_int64_t() */
245 static void read_long(long *val)
246 {
247 bool result;
248
249 if (xdrp)
250 result = xdr_long(xdrp, val);
251 else
252 result = fread(val, sizeof(long int), 1, file) == 1;
253
254 /* Now, what to do with RESULT? ;) */
255 (void)result;
256 }
257 #endif
258
_AI_read_string(char * string,int strileng)259 void _AI_read_string(char *string, int strileng)
260 {
261 bool result;
262
263 if (xdrp) {
264 result = xdr_string(xdrp, &string, strileng+1);
265 } else {
266 result = fread(string, sizeof(char), strileng+1, file) == (strileng+1);
267 }
268
269 /* Now, what to do with RESULT? ;) */
270 (void)result;
271 }
272
_AI_read_var_string(char ** string)273 void _AI_read_var_string(char **string)
274 {
275 bool result;
276 int strileng;
277
278 _AI_read_int(&strileng);
279 if (strileng) {
280 *string = (char *)MEM_ALLOC(strileng+1, char);
281 }
282 if (xdrp) {
283 result = xdr_string(xdrp, string, strileng+1);
284 } else {
285 result = fread(*string, sizeof(char), strileng+1, file) == (strileng+1);
286 }
287
288 /* Now, what to do with RESULT? ;) */
289 (void)result;
290 }
291
_AI_read_vector(void * start,int n,size_t size,xdrproc_t xdrproc)292 void _AI_read_vector(void *start, int n, size_t size, xdrproc_t xdrproc)
293 {
294 bool result;
295
296 if (xdrp)
297 result = xdr_vector(xdrp, (char *)start, n, size, xdrproc);
298 else
299 result = fread(start, size, n, file) == n;
300
301 /* Now, what to do with RESULT? ;) */
302 (void)result;
303 }
304
_AI_read_U_CHAR(U_CHAR * val)305 void _AI_read_U_CHAR(U_CHAR *val)
306 {
307 bool result;
308
309 if (xdrp)
310 result = AI_xdr_U_CHAR(xdrp, val);
311 else
312 result = fread(val, sizeof(U_CHAR), 1, file) == 1;
313
314 /* Now, what to do with RESULT? ;) */
315 (void)result;
316 }
317
_AI_read_S_CHAR(S_CHAR * val)318 void _AI_read_S_CHAR(S_CHAR *val)
319 {
320 bool result;
321
322 if (xdrp)
323 result = AI_xdr_S_CHAR(xdrp, val);
324 else
325 result = fread(val, sizeof(S_CHAR), 1, file) == 1;
326
327 /* Now, what to do with RESULT? ;) */
328 (void)result;
329 }
330
331 /*--------------------------------------------------------------------------*/
332
333 static int n_vert_dofs;
334 static DOF **vert_dofs;
335
336 static int n_edge_dofs;
337 static DOF **edge_dofs;
338
339 static int n_face_dofs;
340 static DOF **face_dofs;
341
342 /*--------------------------------------------------------------------------*/
343
344 /****************************************************************************/
345 /* read_dofs(mesh, dof_ptr, type): Read a DOF vector of type "type" from */
346 /* the file. Allocate space for nonzero DOF_entries (which are overwritten */
347 /* with the indices from the file). */
348 /****************************************************************************/
349
read_dofs(MESH * mesh,DOF ** dof_ptr,int type)350 static void read_dofs(MESH *mesh, DOF **dof_ptr, int type)
351 {
352 int n_dof = mesh->n_dof[type];
353 DOF temp_dofs[n_dof];
354 DOF_ADMIN *admin;
355 int i, j, n, n0;
356
357 _AI_read_vector(temp_dofs, n_dof, sizeof(DOF), (xdrproc_t)AI_xdr_DOF);
358 *dof_ptr = AI_get_dof_memory(mesh, type);
359
360 for (i = 0; i < mesh->n_dof_admin; i++) {
361 admin = mesh->dof_admin[i];
362
363 n = admin->n_dof[type];
364 n0 = admin->n0_dof[type];
365
366 TEST_EXIT(n+n0 <= n_dof,
367 "dof_admin \"%s\": n=%d, n0=%d too large: ndof=%d\n",
368 admin->name, n, n0, n_dof);
369
370 for (j = 0; j < n; j++) {
371 (*dof_ptr)[n0+j] = temp_dofs[n0+j];
372 }
373 }
374 }
375
376 /*--------------------------------------------------------------------------*/
377
read_el_recursive(MESH * mesh,EL * parent,int level)378 static EL *read_el_recursive(MESH *mesh, EL *parent, int level)
379 {
380 FUNCNAME("read_el_recursive");
381 int dim = mesh->dim;
382 int i, j, n, node0;
383 EL *el;
384 U_CHAR uc, nc;
385
386 el = get_element(mesh);
387 mesh->n_hier_elements++;
388
389 #if ALBERTA_DEBUG
390 el->index = mesh->n_hier_elements;
391 #endif
392
393 _AI_read_U_CHAR(&uc);
394
395 if (dim > 1) {
396 _AI_read_U_CHAR(&nc);
397 if (nc) {
398 el->new_coord = get_real_d(mesh);
399
400 _AI_read_vector(el->new_coord, DIM_OF_WORLD, sizeof(REAL),
401 (xdrproc_t)AI_xdr_REAL);
402 }
403 else
404 el->new_coord = NULL;
405 }
406
407 if (mesh->n_dof[VERTEX] > 0) {
408 node0 = mesh->node[VERTEX];
409 for (i = 0; i < N_VERTICES(dim); i++) {
410 _AI_read_int(&j);
411
412 TEST_EXIT(j < n_vert_dofs,
413 "vert_dofs index too large: %d >= %d\n", j, n_vert_dofs);
414 el->dof[node0 + i] = vert_dofs[j];
415 }
416 }
417
418 if (dim > 1 && mesh->n_dof[EDGE] > 0) {
419 node0 = mesh->node[EDGE];
420 for (i = 0; i < N_EDGES(dim); i++) {
421 _AI_read_int(&j);
422
423 TEST_EXIT(j < n_edge_dofs,
424 "edge_dofs index too large: %d >= %d\n", j, n_edge_dofs);
425 if (j >= 0)
426 el->dof[node0 + i] = edge_dofs[j];
427 }
428 }
429
430 if (dim == 3 && (n = mesh->n_dof[FACE]) > 0) {
431 node0 = mesh->node[FACE];
432 for (i = 0; i < N_FACES_3D; i++) {
433 _AI_read_int(&j);
434
435 TEST_EXIT(j < n_face_dofs,
436 "face_dofs index too large: %d >= %d\n", j, n_face_dofs);
437 if (j >= 0)
438 el->dof[node0 + i] = face_dofs[j];
439 }
440 }
441
442 if ((n = mesh->n_dof[CENTER]) > 0) {
443 node0 = mesh->node[CENTER];
444
445 read_dofs(mesh, el->dof + node0, CENTER);
446 }
447
448 if (uc) {
449 el->child[0] = read_el_recursive(mesh, el, level+1);
450 el->child[1] = read_el_recursive(mesh, el, level+1);
451 }
452 else
453 mesh->n_elements++;
454
455 return el;
456 }
457
458 /*--------------------------------------------------------------------------*/
459
read_dof_admins(MESH * mesh)460 static void read_dof_admins(MESH *mesh)
461 {
462 FUNCNAME("read_dof_admins");
463 int i, n_dof_admin, iadmin, used_count;
464 int n_dof_el, n_dof[N_NODE_TYPES], n_node_el, node[N_NODE_TYPES];
465 int a_n_dof[N_NODE_TYPES];
466 U_CHAR flags;
467 char *name;
468 DOF_ADMIN *admin;
469
470
471 _AI_read_int(&n_dof_el);
472 _AI_read_vector(n_dof, N_NODE_TYPES, sizeof(int), (xdrproc_t)AI_xdr_int);
473 _AI_read_int(&n_node_el);
474 _AI_read_vector(node, N_NODE_TYPES, sizeof(int), (xdrproc_t)AI_xdr_int);
475 /* use data later for check */
476
477 _AI_read_int(&n_dof_admin);
478 for (iadmin = 0; iadmin < n_dof_admin; iadmin++) {
479 _AI_read_vector(a_n_dof, N_NODE_TYPES, sizeof(int), (xdrproc_t)AI_xdr_int);
480 _AI_read_int(&used_count);
481
482 _AI_read_var_string(&name);
483
484 _AI_read_U_CHAR(&flags);
485
486 admin = AI_get_dof_admin(mesh, name, a_n_dof);
487 admin->flags = (FLAGS)flags;
488
489 MEM_FREE(name, strlen(name)+1, char);
490
491 if (used_count > 0) {
492 enlarge_dof_lists(admin, used_count);
493 }
494 admin->used_count = used_count;
495 } /* end for (iadmin) */
496
497 for (i = 0; i < N_NODE_TYPES; i++) {
498 if (mesh->n_dof[i]) {
499 AI_get_dof_list(mesh, i);
500 }
501 }
502
503 AI_get_dof_ptr_list(mesh);
504
505 TEST(mesh->n_dof_el == n_dof_el,"wrong n_dof_el: %d %d\n",
506 mesh->n_dof_el, n_dof_el);
507 for (i = 0; i < N_NODE_TYPES; i++)
508 TEST(mesh->n_dof[i] == n_dof[i],"wrong n_dof[%d]: %d %d\n",
509 i, mesh->n_dof[i], n_dof[i]);
510 TEST(mesh->n_node_el == n_node_el,"wrong n_node_el: %d %d\n",
511 mesh->n_node_el, n_node_el);
512 for (i = 0; i < N_NODE_TYPES; i++)
513 TEST(mesh->node[i] == node[i],"wrong node[%d]: %d %d\n",
514 i, mesh->node[i], node[i]);
515 return;
516 }
517
518 /*--------------------------------------------------------------------------*/
519
520 static MESH *
read_mesh_master(REAL * timeptr,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int),MESH * master)521 read_mesh_master(REAL *timeptr,
522 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
523 MESH *master){
524 FUNCNAME("read_mesh_master");
525 MACRO_EL *mel;
526 int i, j, n;
527 REAL_D *v;
528 int neigh_i[N_NEIGH_MAX];
529 char *name, s[sizeof(ALBERTA_VERSION)];
530 int dim, iDIM_OF_WORLD, ne, nv;
531 REAL time;
532 REAL_D diam, bbox[2];
533 int n_vert, n_elements, n_hier_elements;
534 int n_edges;
535 int *vert_i, *mel_vertices = NULL;
536 static int funccount=0;
537 int n_faces, max_edge_neigh;
538 S_CHAR is_periodic;
539 static NODE_PROJECTION id_proj = { NULL };
540 MESH *mesh = NULL;
541 /*
542 #define ALBERTA_VERSION "ALBERTA: Version 1.2 "
543 #define ALBERTA_VERSION "ALBERTA: Version 1.00"
544 #define ALBERTA_VERSION "ALBERTA: Version 2.0 "
545 #define ALBERTA_VERSION "ALBERTA: Version 2.1 "
546 */
547 int major, minor;
548 char *vers;
549
550 _AI_read_string(s, sizeof(ALBERTA_VERSION)-1);
551
552 if (memcmp(s, ALBERTA_MAGIC, sizeof(ALBERTA_MAGIC)-1) != 0) {
553 /* try and read again, hoping the xdr pointer is still at the
554 * right position. ugly, but works.
555 */
556 WARNING("Invalid file id: \"%s\", expected \""ALBERTA_MAGIC"X.YZ\"\n", s);
557 MSG("Retrying in ALBERTA-1.2 compatibility mode ...\n");
558 AI_xdr_close(xdrp);
559
560 /* xdrrec_endofrecord(xdrp, 1); */
561
562 rewind( file /*stream*/);
563
564 if (!(xdrp = MEM_ALLOC(1,XDR))) {
565 ERROR("can't allocate memory for xdr pointer.\n");
566 return NULL;
567 }
568
569 xdrstdio_create(xdrp, file, XDR_DECODE);
570 xdrrec_create(xdrp, 65536, 65536, (void *)file,
571 (int (*)())read_xdr_file, (int (*)())write_xdr_file);
572
573 xdrp->x_op = XDR_DECODE;
574 xdrp->x_public = (caddr_t)file;
575
576 if (XDR_DECODE)
577 xdrrec_skiprecord(xdrp);
578 /*MSG("retry...\n");*/
579 _AI_read_string(s, sizeof(ALBERTA_VERSION)-1);
580 if (memcmp(s, ALBERTA_MAGIC, sizeof(ALBERTA_MAGIC)-1) != 0) {
581 ERROR("failed... AGAIN!\nabort...\n");
582 goto error_exit;
583 } else {
584 /*MSG("this time it works.\n");*/
585 }
586 }
587 vers = s + sizeof(ALBERTA_MAGIC)-1;
588 major = (int)(*vers++ - '0');
589 if (major < 1 || major > 2) {
590 ERROR("Unsupported ALBERTA major version: %d\n", major);
591 goto error_exit;
592 }
593 if (*vers++ != '.') {
594 ERROR("Invalid file id: \"%s\", expected \""ALBERTA_MAGIC"X.YZ\"\n", s);
595 ERROR("*** vers: '%c'\n", vers[-1]);
596 goto error_exit;
597 }
598 minor = (int)(*vers++ - '0');
599 if (minor < 0 || minor > 9) {
600 ERROR("Invalid file id: \"%s\", expected \""ALBERTA_MAGIC"X.YZ\"\n", s);
601 ERROR("*** vers: '%c'\n", vers[-1]);
602 goto error_exit;
603 }
604 if (*vers != ' ') {
605 if (major == 2 && minor == 0 && (int)*vers == 0) {
606 #if ALBERTA_DEBUG
607 /* V2.0 versions are one character too short */
608 WARNING("ALBERTA V2.0 bug detected, working around it.\n");
609 #endif
610 } else {
611 if (*vers < '0' || *vers > '9') {
612 ERROR("Invalid file id: \"%s\", expected \""ALBERTA_MAGIC"X.YZ\"\n", s);
613 ERROR("*** vers: '%c' (%d)\n", *vers, (int)*vers);
614 goto error_exit;
615 }
616 minor *= 10;
617 minor += (int)(*vers - '0');
618 }
619 }
620
621 #if 1
622 if (major < 2) {
623 #define SUP12 1
624 #if SUP12
625 if (major == 1 && minor == 2) {
626 return _AI_read_mesh_1_2(timeptr);
627 }
628 #endif
629 #if !SUP12
630 ERROR("ALBERTA V1 compatibility mode not yet implemented, sorry\n");
631 goto error_exit;
632 #endif
633 }
634 #endif
635 if (major >= 2) {
636 if (major > 2 || minor > 3) {
637 ERROR("Post-V2.3 ALBERTA file version detected (V%d.%d), aborting\n",
638 major, minor);
639 goto error_exit;
640 }
641 }
642
643 _AI_read_int(&dim);
644 if (dim > DIM_MAX) {
645 ERROR("dim == %d is greater than DIM_MAX == %d!\n", dim, DIM_MAX);
646 goto error_exit;
647 }
648
649 _AI_read_int(&iDIM_OF_WORLD);
650 if (iDIM_OF_WORLD != DIM_OF_WORLD) {
651 ERROR("wrong DIM_OF_WORLD %d. abort.\n", iDIM_OF_WORLD);
652 goto error_exit;
653 }
654
655 _AI_read_REAL(&time);
656 if (timeptr) *timeptr = time;
657
658 _AI_read_int(&i); /* length without terminating \0 */
659 if (i) {
660 name = MEM_ALLOC(i+1, char);
661 _AI_read_string(name, i);
662 }
663 else {
664 funccount++;
665 i=100;
666 name = MEM_ALLOC(i+1, char);
667 sprintf(name, "READ_MESH%d", funccount);
668 }
669
670 _AI_read_int(&n_vert);
671 if (dim > 1) {
672 _AI_read_int(&n_edges);
673 }
674
675 _AI_read_int(&n_elements);
676 _AI_read_int(&n_hier_elements);
677
678 if (dim == 3) {
679 _AI_read_int(&n_faces);
680 _AI_read_int(&max_edge_neigh);
681 }
682
683 mesh = GET_MESH(dim, name, NULL, NULL, NULL);
684
685 _AI_read_S_CHAR(&is_periodic);
686 mesh->is_periodic = is_periodic != 0;
687 if (mesh->is_periodic) {
688 _AI_read_int(&mesh->per_n_vertices);
689 if (dim > 1) {
690 _AI_read_int(&mesh->per_n_edges);
691 if (dim > 2) {
692 _AI_read_int(&mesh->per_n_faces);
693 }
694 }
695 _AI_read_int(&mesh->n_wall_trafos);
696 if (mesh->n_wall_trafos) {
697 int wt;
698
699 mesh->wall_trafos = MEM_ALLOC(mesh->n_wall_trafos, AFF_TRAFO *);
700 ((AFF_TRAFO **)mesh->wall_trafos)[0] =
701 MEM_ALLOC(mesh->n_wall_trafos, AFF_TRAFO);
702 for (wt = 0; wt < mesh->n_wall_trafos; wt++) {
703 ((AFF_TRAFO **)mesh->wall_trafos)[wt] = mesh->wall_trafos[0]+wt;
704 _AI_read_vector(mesh->wall_trafos[wt],
705 SQR(DIM_OF_WORLD) + DIM_OF_WORLD, sizeof(REAL),
706 (xdrproc_t)AI_xdr_REAL);
707 }
708 }
709 }
710
711 _AI_read_vector(mesh->bbox, 2*DIM_OF_WORLD, sizeof(REAL), (xdrproc_t)AI_xdr_REAL);
712 AXPBY_DOW(1.0, mesh->bbox[1], -1.0, mesh->bbox[0], mesh->diam);
713
714 read_dof_admins(mesh);
715
716 MEM_FREE(name, i+1, char);
717
718 _AI_read_int(&n_vert_dofs);
719
720 if (n_vert_dofs > 0) {
721 vert_dofs = MEM_ALLOC(n_vert_dofs, DOF *);
722 n = mesh->n_dof[VERTEX];
723 for (i = 0; i < n_vert_dofs; i++) {
724 vert_dofs[i] = _AI_get_dof(mesh, VERTEX, false);
725
726 _AI_read_vector(vert_dofs[i], n, sizeof(DOF), (xdrproc_t)AI_xdr_DOF);
727 }
728 }
729
730 if (dim > 1) {
731 _AI_read_int(&n_edge_dofs);
732
733 if (n_edge_dofs > 0) {
734 edge_dofs = MEM_ALLOC(n_edge_dofs, DOF *);
735
736 for (i = 0; i < n_edge_dofs; i++) {
737 read_dofs(mesh, edge_dofs + i, EDGE);
738 }
739 }
740 }
741
742 if (dim == 3) {
743 _AI_read_int(&n_face_dofs);
744
745 if (n_face_dofs > 0) {
746 face_dofs = MEM_ALLOC(n_face_dofs, DOF *);
747
748 for (i = 0; i < n_face_dofs; i++) {
749 read_dofs(mesh, face_dofs + i, FACE);
750 }
751 }
752 }
753
754 /* After we have read all DOFs we make sure that we have allocated
755 * the correct number of DOFs for each admin.
756 */
757 for (i = 0; i < mesh->n_dof_admin; i++) {
758 _AI_allocate_n_dofs(mesh->dof_admin[i], mesh->dof_admin[i]->used_count);
759 }
760
761 _AI_read_int(&ne);
762 _AI_read_int(&nv);
763
764 ((MESH_MEM_INFO *)mesh->mem_info)->count = nv;
765 v = ((MESH_MEM_INFO *)mesh->mem_info)->coords = MEM_ALLOC(nv, REAL_D);
766
767 for (i = 0; i < nv; i++) {
768 _AI_read_vector(v[i], DIM_OF_WORLD, sizeof(REAL), (xdrproc_t)AI_xdr_REAL);
769 }
770
771 for (j = 0; j < DIM_OF_WORLD; j++) {
772 bbox[0][j] = 1.E30;
773 bbox[1][j] = -1.E30;
774 }
775
776 for (i = 0; i < nv; i++) {
777 for (j = 0; j < DIM_OF_WORLD; j++) {
778 bbox[0][j] = MIN(bbox[0][j], v[i][j]);
779 bbox[1][j] = MAX(bbox[1][j], v[i][j]);
780 }
781 }
782
783 AXPBY_DOW(1.0, bbox[1], -1.0, bbox[0], diam);
784
785 mel = MEM_CALLOC(ne, MACRO_EL);
786
787 mesh->n_macro_el = ne;
788 mesh->macro_els = mel;
789
790 vert_i = mel_vertices = MEM_ALLOC(ne*N_VERTICES(dim), int);
791
792 for (n = 0; n < ne; n++) {
793 mel[n].index = n;
794
795 _AI_read_vector(vert_i, N_VERTICES(dim), sizeof(int), (xdrproc_t)AI_xdr_int);
796
797 for (i = 0; i < N_VERTICES(dim); i++) {
798 if (*vert_i >= 0 && *vert_i < nv) {
799 mel[n].coord[i] = v + *vert_i++;
800 } else {
801 #if 1
802 ERROR("Macro element vertex numbering is garbled");
803 goto error_exit;
804 #else
805 /* Error case???? */
806 mel[n].coord[i] = NULL;
807 #endif
808 }
809 }
810
811 if (major == 2 && minor == 0) {
812 S_CHAR dummy[N_VERTICES(dim)];
813 _AI_read_vector(dummy, N_VERTICES(dim), sizeof(S_CHAR),
814 (xdrproc_t)AI_xdr_S_CHAR);
815 }
816 _AI_read_vector(mel[n].wall_bound, N_WALLS(dim), sizeof(BNDRY_TYPE),
817 (xdrproc_t)AI_xdr_S_CHAR);
818 #if DIM_MAX > 2
819 if (dim == 3 && major == 2 && minor == 0) {
820 S_CHAR dummy[N_EDGES_3D];
821 _AI_read_vector(dummy, N_EDGES_3D, sizeof(S_CHAR), (xdrproc_t)AI_xdr_S_CHAR);
822 }
823 #endif
824
825 _AI_read_vector(neigh_i, N_NEIGH(dim), sizeof(int), (xdrproc_t)AI_xdr_int);
826
827 if (init_node_proj) {
828 mel[n].projection[0] = init_node_proj(mesh, mel + n, 0);
829 }
830 for (i = 0; i < N_NEIGH(dim); i++) {
831 if (init_node_proj && (dim > 1)) {
832 mel[n].projection[i+1] = init_node_proj(mesh, mel + n, i+1);
833 }
834
835 if ((neigh_i[i] >= 0) && (neigh_i[i] < ne)) {
836 mel[n].neigh[i] = mel + (neigh_i[i]);
837 } else {
838 mel[n].neigh[i] = NULL;
839 }
840 }
841
842 _AI_read_vector(mel[n].opp_vertex, N_NEIGH(dim), sizeof(S_CHAR),
843 (xdrproc_t)AI_xdr_S_CHAR);
844
845 /* V2.0 compatibility code. */
846 if (major == 2 && minor == 0) {
847 if (true) {
848 for (i = 0; i < N_NEIGH(dim); i++) {
849 _AI_read_vector(mel[n].neigh_vertices[i],
850 N_VERTICES(dim-1), sizeof(S_CHAR),
851 (xdrproc_t)AI_xdr_S_CHAR);
852 }
853 }
854 if (mesh->is_periodic) {
855 int wt_i[N_WALLS_MAX];
856 _AI_read_vector(wt_i, N_WALLS(dim), sizeof(int), (xdrproc_t)xdr_int);
857 for (i = 0; i < N_WALLS(dim); i++) {
858 if (wt_i[i] != -1) {
859 mel[n].wall_trafo[i] = mesh->wall_trafos[wt_i[i]];
860 } else {
861 mel[n].wall_trafo[i] = NULL;
862 }
863 }
864 }
865 }
866
867 #if DIM_MAX >= 3
868 if (dim == 3) {
869 _AI_read_U_CHAR(&(mel[n].el_type));
870 mel[n].orientation = AI_get_orientation_3d(&mel[n]);
871 }
872 #endif
873
874 if (major == 2 && minor >= 2) {
875 /* file-version V2.2 */
876 S_CHAR pr_status;
877 int mst_idx;
878
879 for (i = 0; i < 1+N_WALLS(dim); i++) {
880 _AI_read_S_CHAR(&pr_status);
881 if (init_node_proj == NULL) {
882 /* application specified projection binding has precedence. */
883 if (pr_status) {
884 mel[n].projection[i] = &id_proj;
885 }
886 }
887 }
888 _AI_read_int(&mst_idx);
889 _AI_read_S_CHAR(&mel[n].master.opp_vertex);
890 if (master) {
891 mel[n].master.macro_el = &master->macro_els[mst_idx];
892 }
893 }
894
895 /* >= V2.1 format */
896 if (major == 2 && minor >= 1) {
897 if (mesh->is_periodic) {
898 int wt_i[N_WALLS_MAX];
899
900 for (i = 0; i < N_NEIGH(dim); i++) {
901 _AI_read_vector(mel[n].neigh_vertices[i],
902 N_VERTICES(dim-1), sizeof(S_CHAR),
903 (xdrproc_t)AI_xdr_S_CHAR);
904 }
905 _AI_read_vector(wt_i, N_WALLS(dim), sizeof(int), (xdrproc_t)xdr_int);
906 for (i = 0; i < N_WALLS(dim); i++) {
907 if (wt_i[i] != -1) {
908 mel[n].wall_trafo[i] = mesh->wall_trafos[wt_i[i]];
909 } else {
910 mel[n].wall_trafo[i] = NULL;
911 }
912 }
913 } else {
914 /* Don't forget to set periodic info to invalid values (i.e. -1) */
915 memset(mel[n].neigh_vertices, -1, sizeof(mel[n].neigh_vertices));
916 }
917 }
918
919 mel[n].el = read_el_recursive(mesh, NULL, 0);
920 }
921
922 /* Generate the boundary classification for all lower-dimensional
923 * sub-simplices (vertices, in 3d: also edges).
924 */
925 _AI_fill_bound_info(mesh, mel_vertices, nv, ne, false);
926
927 /****************************************************************************/
928 /* The present mechanism only reads DOFs from the file which were used by */
929 /* at least one DOF_ADMIN. The missing element DOF pointers are filled by */
930 /* this routine (in memory.c). */
931 /****************************************************************************/
932 if (dim > 0) {
933 AI_fill_missing_dofs(mesh);
934 }
935
936 if (n_elements != mesh->n_elements) {
937 ERROR("n_elements != mesh->n_elements.\n");
938 goto error_exit;
939 }
940
941 if (n_hier_elements != mesh->n_hier_elements)
942 {
943 ERROR("n_hier_elements != mesh->n_hier_elements.\n");
944 goto error_exit;
945 }
946
947
948 if (mesh->n_dof[VERTEX])
949 {
950 if (n_vert != n_vert_dofs)
951 {
952 ERROR("n_vertices == %d != %d == n_vert_dofs.\n", n_vert, n_vert_dofs);
953 mesh->n_vertices = n_vert_dofs;
954 goto error_exit;
955 }
956 }
957 mesh->n_vertices = n_vert;
958
959 if (dim > 1)
960 mesh->n_edges = n_edges;
961
962 if (dim == 3) {
963 mesh->n_faces = n_faces;
964
965 mesh->max_edge_neigh = max_edge_neigh;
966 }
967
968 for (i=0; i < DIM_OF_WORLD; i++) {
969 if (ABS(mesh->diam[i]-diam[i]) > (mesh->diam[i]/10000.0)) {
970 WARNING("diam[%i] != mesh->diam[%i]. Parametric mesh?\n", i, i);
971 /* goto error_exit; */ /* nope, mesh could be parametric */
972 }
973 }
974
975 /* Read the magic cookie. */
976 _AI_read_int((int *)&mesh->cookie);
977
978 for (;;) {
979 _AI_read_string(s, 4);
980
981 /* file-version >= V2.2 */
982 if (major == 2 && minor >= 2 && memcmp(s, "LAPA", 4) == 0) {
983 S_CHAR not_all, use_reference;
984 DOF_REAL_D_VEC *coords;
985 DOF_UCHAR_VEC *touched_edges = NULL;
986 DOF_PTR_VEC *edge_pr;
987 /* PARAM_STRATEGY */ int strategy;
988 FLAGS flags;
989
990 _AI_read_S_CHAR(¬_all); /* not_all */
991 _AI_read_S_CHAR(&use_reference); /* use_reference_mesh */
992 _AI_read_int(&strategy);
993 coords = (DOF_REAL_D_VEC *)
994 read_dof_vec_master(__DOF_REAL_D_VEC, NULL, mesh, NULL, false);
995 if (not_all && coords->fe_space->bas_fcts->degree > 1) {
996 touched_edges = (DOF_UCHAR_VEC *)
997 read_dof_vec_master(__DOF_UCHAR_VEC, NULL, mesh, NULL, false);
998 }
999
1000 flags = strategy;
1001 if (mesh->is_periodic &&
1002 (coords->fe_space->admin->flags & ADM_PERIODIC) != 0) {
1003 flags |= PARAM_PERIODIC_COORDS;
1004 }
1005
1006 /* BUG: the "selective" feature will not really work. But so what. */
1007 use_lagrange_parametric(mesh,
1008 coords->fe_space->bas_fcts->degree,
1009 NULL /* i.e. not selective */,
1010 flags);
1011
1012 dof_copy_d(coords, get_lagrange_coords(mesh));
1013 free_dof_real_d_vec(coords);
1014
1015 edge_pr = get_lagrange_edge_projections(mesh);
1016 if (edge_pr != NULL) {
1017 static const NODE_PROJECTION dummy_projection = { NULL, };
1018 TEST_EXIT(touched_edges != NULL,
1019 "No touched-edges vector?\n");
1020 FOR_ALL_DOFS(touched_edges->fe_space->admin, {
1021 edge_pr->vec[dof] =
1022 touched_edges->vec[dof] ? (void *)&dummy_projection : NULL;
1023 });
1024 free_dof_uchar_vec(touched_edges);
1025
1026 /* Fill in the real projection pointers using a mesh loop. Is
1027 * there a better (more efficient) way?
1028 */
1029 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_MACRO_WALLS) {
1030 /* 1d: edge_projection refers to the "bulk" projection, if
1031 * only vertices (i.e. walls) are projected, then the
1032 * element stays affine.
1033 *
1034 * 2d: edge_projection refers to the bulk and wall
1035 * projections, where the edges are the walls.
1036 *
1037 * 3d: edge_projection is slightly more complicated.
1038 */
1039 switch (mesh->dim) {
1040 case 1: {
1041 int node_c = mesh->node[CENTER];
1042 int n0_edge_pr = edge_pr->fe_space->admin->n0_dof[CENTER];
1043 DOF edge_dof = el_info->el->dof[node_c][n0_edge_pr];
1044
1045 if (edge_pr->vec[edge_dof] != NULL) {
1046 edge_pr->vec[edge_dof] = (void *)wall_proj(el_info, -1);
1047 }
1048 break;
1049 }
1050 case 2: {
1051 int node_e = mesh->node[EDGE];
1052 int n0_edge_pr = edge_pr->fe_space->admin->n0_dof[EDGE];
1053
1054 int e;
1055 for (e = 0; e < N_EDGES_2D; n++) {
1056 DOF edge_dof = el_info->el->dof[node_e+e][n0_edge_pr];
1057
1058 if (edge_pr->vec[edge_dof] != NULL) {
1059 const NODE_PROJECTION *act_proj = wall_proj(el_info, e);
1060 if (act_proj == NULL) {
1061 act_proj = wall_proj(el_info, -1);
1062 }
1063 if (act_proj) {
1064 edge_pr->vec[edge_dof] = (void *)act_proj;
1065 }
1066 }
1067 }
1068 break;
1069 }
1070 case 3: {
1071 int node_e = mesh->node[EDGE];
1072 int n0_edge_pr = edge_pr->fe_space->admin->n0_dof[EDGE];
1073
1074 int e;
1075 for (e = 0; e < N_EDGES_3D; ++e) {
1076 DOF edge_pr_dof = el_info->el->dof[node_e+e][n0_edge_pr];
1077
1078 if (edge_pr->vec[edge_pr_dof] != NULL) {
1079 /* Look for a projection function that applies to edge e. */
1080 const NODE_PROJECTION *act_proj;
1081
1082 act_proj = wall_proj(el_info, face_of_edge_3d[e][0]);
1083 if (!act_proj) {
1084 act_proj = wall_proj(el_info, face_of_edge_3d[e][1]);
1085 }
1086 if (!act_proj) {
1087 act_proj = wall_proj(el_info, -1);
1088 }
1089
1090 if (act_proj != NULL) {
1091 /* replace by actual projection, but keep the
1092 * projected status if that was stored on disk.
1093 */
1094 edge_pr->vec[edge_pr_dof] = (void *)act_proj;
1095 }
1096 }
1097 }
1098 break;
1099 }
1100 } /* end switch (mesh->dim) */
1101 } TRAVERSE_NEXT();
1102 }
1103
1104 continue;
1105 }
1106
1107 /* file-version >= V2.3 */
1108 if (major == 2 && minor >= 2 && memcmp(s, "TMSH", 4) == 0) {
1109 int trace_id;
1110 MESH *slave;
1111 _AI_read_int((int *)&trace_id);
1112 /* This is a trace mesh, read it in */
1113 slave = read_mesh_master(NULL, init_node_proj, mesh);
1114 slave->trace_id = trace_id;
1115 ((MESH_MEM_INFO *)mesh->mem_info)->next_trace_id = trace_id+1;
1116
1117 continue;
1118 }
1119
1120 if (memcmp(s, "EOF.", 4) == 0) {
1121 break;
1122 }
1123
1124 ERROR("no FILE END MARK.\n");
1125 goto error_exit;
1126 }
1127
1128 if (master) {
1129 bind_submesh(master, mesh, NULL, NULL);
1130 }
1131
1132 return mesh;
1133
1134 error_exit:
1135 free_mesh(mesh);
1136 if (mel_vertices) {
1137 MEM_FREE(mel_vertices, nv, int);
1138 }
1139
1140 return NULL;
1141 }
1142
1143 MESH *
fread_mesh_xdr(FILE * fp,REAL * timeptr,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int),MESH * master)1144 fread_mesh_xdr(FILE *fp, REAL *timeptr,
1145 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
1146 MESH *master)
1147 {
1148 FUNCNAME("fread_mesh_xdr");
1149 MESH *mesh;
1150
1151 /* NOTE: xdrp is global */
1152 if ((xdrp = AI_xdr_fopen(fp, XDR_DECODE)) == NULL) {
1153 ERROR("Cannot convert file pointer to XDR handle\n");
1154 return NULL;
1155 }
1156 file = fp;
1157
1158 mesh = read_mesh_master(timeptr, init_node_proj, master);
1159
1160 AI_xdr_close(xdrp);
1161
1162 xdrp = NULL; /* do not close fp! */
1163 file = NULL;
1164
1165 return mesh;
1166 }
1167
read_mesh_xdr(const char * fn,REAL * timeptr,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int),MESH * master)1168 MESH *read_mesh_xdr(const char *fn, REAL *timeptr,
1169 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
1170 MESH *master)
1171 {
1172 FUNCNAME("read_mesh_xdr");
1173 MESH *mesh;
1174 FILE *fp;
1175
1176 if ((fp = fopen(fn, "rb")) == NULL) {
1177 ERROR("Cannot open file '%s'\n", fn);
1178 return NULL;
1179 }
1180
1181 mesh = fread_mesh_xdr(fp, timeptr, init_node_proj, master);
1182
1183 if (mesh) {
1184 MSG("File %s read.\n", fn);
1185 }
1186
1187 fclose(fp);
1188
1189 return mesh;
1190 }
1191
fread_mesh(FILE * fp,REAL * timeptr,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int),MESH * master)1192 MESH *fread_mesh(FILE *fp, REAL *timeptr,
1193 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
1194 MESH *master)
1195 {
1196 MESH *mesh;
1197 file = fp;
1198 xdrp = NULL;
1199
1200 mesh = read_mesh_master(timeptr, init_node_proj, master);
1201
1202 file = NULL; /* but do not close the file */
1203
1204 return mesh;
1205 }
1206
read_mesh(const char * fn,REAL * timeptr,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int),MESH * master)1207 MESH *read_mesh(const char *fn, REAL *timeptr,
1208 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
1209 MESH *master)
1210 {
1211 FUNCNAME("fread_mesh");
1212 MESH *mesh;
1213 FILE *fp;
1214
1215 if ((fp = fopen(fn, "rb")) == NULL) {
1216 ERROR("Cannot open file '%s'\n", fn);
1217 return NULL;
1218 }
1219
1220 mesh = fread_mesh(fp, timeptr, init_node_proj, master);
1221
1222 fclose(fp);
1223
1224 if (mesh) {
1225 MSG("File %s read.\n", fn);
1226 }
1227
1228 return mesh;
1229 }
1230
1231 /*--------------------------------------------------------------------------*/
1232 /* read DOF vectors of various types */
1233 /*--------------------------------------------------------------------------*/
1234
1235 typedef void (*ADD_TO_ADM_FCT)(DOF_VEC *dv, const DOF_ADMIN *admin);
1236 typedef DOF_VEC *(*GET_DV_FCT)(const char *name, const FE_SPACE *fesp);
1237
1238 static struct dv_read_struct {
1239 const char *name;
1240 ADD_TO_ADM_FCT add_to_admin;
1241 GET_DV_FCT get_dof_vec;
1242 xdrproc_t item_proc;
1243 size_t item_size;
1244 int stride;
1245 } dv_read_table[] = {
1246 /* [__DOF_REAL_VEC] = */ {
1247 "DOF_REAL_VEC ",
1248 (ADD_TO_ADM_FCT)add_dof_real_vec_to_admin,
1249 (GET_DV_FCT)get_dof_real_vec,
1250 (xdrproc_t)AI_xdr_REAL, sizeof(REAL), 1,
1251 },
1252 /* [__DOF_REAL_D_VEC] = */ {
1253 "DOF_REAL_D_VEC ",
1254 (ADD_TO_ADM_FCT)add_dof_real_d_vec_to_admin,
1255 (GET_DV_FCT)get_dof_real_d_vec,
1256 (xdrproc_t)AI_xdr_REAL, sizeof(REAL), DIM_OF_WORLD,
1257 },
1258 /* [__DOF_REAL_VEC_D] = */ {
1259 /* dummy, will redirect to __DOF_REAL_VEC or __DOF_REAL_VEC_D */
1260 "DOF_REAL_VEC_D ",
1261 },
1262 /* [__DOF_INT_VEC] = */ {
1263 "DOF_INT_VEC ",
1264 (ADD_TO_ADM_FCT)add_dof_int_vec_to_admin,
1265 (GET_DV_FCT)get_dof_int_vec,
1266 (xdrproc_t)AI_xdr_int, sizeof(int), 1,
1267 },
1268 /* [__DOF_SCHAR_VEC] = */ {
1269 "DOF_SCHAR_VEC ",
1270 (ADD_TO_ADM_FCT)add_dof_schar_vec_to_admin,
1271 (GET_DV_FCT)get_dof_schar_vec,
1272 (xdrproc_t)AI_xdr_S_CHAR, sizeof(S_CHAR), 1,
1273 },
1274 /* [__DOF_UCHAR_VEC] = */ {
1275 "DOF_UCHAR_VEC ",
1276 (ADD_TO_ADM_FCT)add_dof_uchar_vec_to_admin,
1277 (GET_DV_FCT)get_dof_uchar_vec,
1278 (xdrproc_t)AI_xdr_U_CHAR, sizeof(U_CHAR), 1,
1279 },
1280 };
1281
read_dof_vec_master(DV_TYPE_ENUM dv_type,DOF_VEC * dv,MESH * mesh,FE_SPACE * fe_space,bool expect_next)1282 static DOF_VEC *read_dof_vec_master(DV_TYPE_ENUM dv_type,
1283 DOF_VEC *dv,
1284 MESH *mesh, FE_SPACE *fe_space,
1285 bool expect_next)
1286 {
1287 FUNCNAME("read_dof_vec_master");
1288 struct dv_read_struct *rd_table;
1289 int i, last, stride, rdim;
1290 int n_dof[N_NODE_TYPES];
1291 unsigned int cookie;
1292 const DOF_ADMIN *admin = NULL;
1293 const BAS_FCTS *bas_fcts;
1294 char *name, *dv_name;
1295 char ss[17]; /* 17 = length of dofvectype with terminating \0 */
1296 char *s = ss;
1297 U_CHAR flags_and_stride = 0, flags;
1298 bool do_free_fe_space = false;
1299 int old_version = 0;
1300
1301 TEST_EXIT(mesh, "no mesh given\n");
1302
1303 if (dv_type == __DOF_REAL_VEC_D || dv_type == __DOF_REAL_D_VEC) {
1304 rdim = DIM_OF_WORLD;
1305 } else {
1306 rdim = 1;
1307 }
1308
1309 TEST_EXIT(fe_space == NULL || fe_space->rdim == rdim,
1310 "The range dimension %d of the given fe-space does not match "
1311 "the range dimension %d of the stored DOF-vector.\n",
1312 fe_space->rdim, rdim);
1313
1314 _AI_read_string(s, 16);
1315
1316
1317
1318 if (dv_type == __DOF_REAL_VEC_D) {
1319 if (memcmp(s, dv_read_table[__DOF_REAL_D_VEC].name, 12) != 0 &&
1320 memcmp(s, dv_read_table[__DOF_REAL_VEC].name, 12) != 0) {
1321 ERROR("invalid file id; %s\n", s);
1322 return NULL;
1323 }
1324 } else {
1325 if (memcmp(s, dv_read_table[dv_type].name, 12) != 0) {
1326
1327 /* if (memcmp(s, ALBERTA_MAGIC, sizeof(ALBERTA_MAGIC)-1) != 0) { */
1328 /* try and read again, hoping the xdr pointer is still at the right
1329 * position. ugly, but works.
1330 */
1331 WARNING("Invalid file id: \"%s\"\n", s);
1332 MSG("Retrying in ALBERTA-1.2 compatibility mode ...\n");
1333 AI_xdr_close(xdrp);
1334
1335 /* xdrrec_endofrecord(xdrp, 1); */
1336
1337 rewind( file /*stream*/);
1338
1339 if (!(xdrp = MEM_ALLOC(1,XDR))) {
1340 ERROR("can't allocate memory for xdr pointer.\n");
1341 return NULL;
1342 }
1343
1344 xdrstdio_create(xdrp, file, XDR_DECODE);
1345 xdrrec_create(xdrp, 65536, 65536, (void *)file,
1346 (int (*)())read_xdr_file, (int (*)())write_xdr_file);
1347
1348 xdrp->x_op = XDR_DECODE;
1349 xdrp->x_public = (caddr_t)file;
1350
1351 if (XDR_DECODE)
1352 xdrrec_skiprecord(xdrp);
1353 /*MSG("retry...\n");*/
1354 _AI_read_string(s, 16);
1355 if (memcmp(s, dv_read_table[dv_type].name, 12) != 0) {
1356 ERROR("failed... AGAIN!\nabort...\n");
1357 return NULL;
1358 } else {
1359 /*MSG("this time it works.\n");*/
1360 old_version = 1;
1361 }
1362 }
1363 }
1364
1365 _AI_read_int(&last); /* length of vector name */
1366 dv_name = MEM_ALLOC(last+1, char);
1367 _AI_read_string(dv_name, last);
1368
1369 if (!old_version) {
1370 _AI_read_U_CHAR(&flags_and_stride);
1371 } else {
1372 /* Goodness, preserve_coarse_dofs formerly was tied to the
1373 * mesh. How to handle this? Do it this way: as the flag was tied
1374 * to the mesh, the compat-read-mesh-code has already set the flag
1375 * in all DOF_ADMINs stored in the mesh file. So simply have a
1376 * look at the first admin (there always is one with DOFs at
1377 * vertices), if that one has the flags set, then all ADMINs must
1378 * have the flag set and we pretend that also this DOF_VECTOR
1379 * needs it.
1380 */
1381 TEST_EXIT(mesh->n_dof_admin != 0, "Not a single DOF_ADMIN in the mesh???\n");
1382 if (mesh->dof_admin[0]->flags & ADM_PRESERVE_COARSE_DOFS) {
1383 flags_and_stride = ADM_PRESERVE_COARSE_DOFS;
1384 }
1385 }
1386 stride = (flags_and_stride & 0x80) ? DIM_OF_WORLD : 1;
1387 flags = flags_and_stride & ADM_FLAGS_MASK;
1388
1389 if (memcmp(s, dv_read_table[__DOF_REAL_D_VEC].name, 12) == 0 &&
1390 stride == 1) {
1391 WARNING("Stride-mismatch, assuming pre-2.2 DOF_REAL_D_VEC\n");
1392 stride = DIM_OF_WORLD;
1393 }
1394
1395 if (dv_type == __DOF_REAL_VEC_D) {
1396 /* Switch to either DOF_REAL_VEC or DOF_REAL_D_VEC, based on the
1397 * "stride" vlaue.
1398 */
1399 dv_type = stride != 1 ? __DOF_REAL_D_VEC : __DOF_REAL_VEC;
1400 }
1401
1402 rd_table = &dv_read_table[dv_type];
1403
1404
1405 if (old_version) {
1406 _AI_read_vector(n_dof, mesh->dim + 1, sizeof(int), (xdrproc_t)AI_xdr_int);
1407 _AI_match_node_types(n_dof);
1408 for (i = mesh->dim + 1; i < N_NODE_TYPES; i++)
1409 n_dof[i] = 0;
1410 } else {
1411 _AI_read_vector(n_dof, N_NODE_TYPES, sizeof(int), (xdrproc_t)AI_xdr_int);
1412 }
1413
1414 _AI_read_int(&last); /* length of BAS_FCTS name */
1415
1416 if (last) {
1417 name = MEM_ALLOC(last+1, char);
1418 _AI_read_string(name, last);
1419
1420 if (fe_space && (bas_fcts = fe_space->bas_fcts)) {
1421 if (strcmp(bas_fcts->name, name) != 0) {
1422 ERROR("invalid name \"%s\" is not given fe_space->bas_fcts->name %s\n",
1423 name, bas_fcts->name);
1424 }
1425 } else { /* no given fe_space or no bas_fcts in given fe_space */
1426 bas_fcts = get_bas_fcts(mesh->dim, name);
1427 TEST_EXIT(bas_fcts != NULL, "cannot get bas_fcts <%s>\n", name);
1428
1429 if (fe_space) { /* use given fe_space */
1430 TEST_EXIT(bas_fcts->rdim <= fe_space->rdim,
1431 "Request for vector valued basis functions <%s> with "
1432 "scalar fe-space.\n", bas_fcts->name);
1433 fe_space->bas_fcts = bas_fcts;
1434 } else { /* create new fe_space */
1435 fe_space = (FE_SPACE *)get_fe_space(mesh, name, bas_fcts, rdim, flags);
1436 TEST_EXIT(fe_space != NULL,
1437 "cannot get fe_space for bas_fcts <%s>\n", name);
1438 do_free_fe_space = true;
1439 }
1440 }
1441 for (i = 0; i < N_NODE_TYPES; i++) {
1442 TEST_EXIT(n_dof[i] == bas_fcts->n_dof[i],
1443 "wrong n_dof in bas_fcts <%s>\n", name);
1444 }
1445 } else { /* no bas_fcts.name in file */
1446 if (fe_space) { /* use given fe_space */
1447 TEST_EXIT(admin = fe_space->admin, "no fe_space->admin");
1448 for (i = 0; i < N_NODE_TYPES; i++) {
1449 TEST_EXIT(n_dof[i] == admin->n_dof[i],
1450 "wrong n_dof in admin <%s>\n", NAME(admin));
1451 }
1452 } else { /* create new fe_space */
1453 fe_space = (FE_SPACE *)get_dof_space(mesh, NULL, n_dof, flags);
1454 TEST_EXIT(fe_space != NULL, "cannot get fe_space for given n_dof\n");
1455 TEST_EXIT(admin = fe_space->admin, "no admin in new fe_space\n");
1456 for (i = 0; i < N_NODE_TYPES; i++) {
1457 TEST_EXIT(n_dof[i] == admin->n_dof[i],
1458 "wrong n_dof in admin <%s>\n", NAME(admin));
1459 }
1460 do_free_fe_space = true;
1461 }
1462 }
1463 TEST_EXIT(fe_space,"still no fe_space\n");
1464 TEST_EXIT(admin = fe_space->admin, "still no admin\n");
1465
1466 dof_compress(mesh);
1467
1468 if (dv == NULL) {
1469 dv = dv_read_table[dv_type].get_dof_vec(dv_name, fe_space);
1470 }
1471 if (dv->name) {
1472 free((char *)dv->name);
1473 }
1474 dv->name = dv_name;
1475
1476 _AI_read_int(&last);
1477 if (last != admin->size_used) {
1478 ERROR("size of dof vector `%s' == %d "
1479 "does not fit to size_used == %d in admin `%s'\n",
1480 dv->name, last, admin->size_used, admin->name);
1481 ERROR_EXIT("cannot read incompatible data\n");
1482 }
1483
1484 if (last) {
1485 _AI_read_vector(
1486 dv->vec, last*rd_table->stride, rd_table->item_size, rd_table->item_proc);
1487 } else {
1488 ERROR("empty dof vector\n");
1489 dv->size = 0;
1490 dv->vec = NULL;
1491 }
1492
1493 /* Read the magic cookie. */
1494 if (!old_version) {
1495 _AI_read_int((int *)&cookie);
1496
1497 if (cookie != mesh->cookie) {
1498 WARNING("Mesh and DOF vector do not seem to match!\n");
1499 }
1500 }
1501 _AI_read_string(s, 4);
1502 if (memcmp(s, "EOF.", 4) != 0 &&
1503 (!expect_next || CHAIN_SINGLE(fe_space) || memcmp(s, "NEXT", 4) != 0)) {
1504 ERROR("no FILE END MARK.\n");
1505 }
1506
1507 if (do_free_fe_space) {
1508 free_fe_space(fe_space);
1509 }
1510
1511 return dv;
1512 }
1513
1514 /*--------------------------------------------------------------------------*/
1515
1516 static inline
fread_dof_vec_master(bool do_xdr,FILE * fp,MESH * mesh,FE_SPACE * fe_space,DV_TYPE_ENUM dv_type)1517 DOF_VEC *fread_dof_vec_master(bool do_xdr,
1518 FILE *fp, MESH *mesh, FE_SPACE *fe_space,
1519 DV_TYPE_ENUM dv_type)
1520 {
1521 FUNCNAME("fread_dof_vec_master");
1522 DOF_VEC *dv_chain;
1523 DOF_VEC *dv;
1524
1525 if (do_xdr) {
1526 if (!(xdrp = AI_xdr_fopen(fp, XDR_DECODE))) {
1527 ERROR("Cannot convert file handle to XDR handle\n");
1528 return NULL;
1529 }
1530 }
1531 file = fp;
1532
1533 dv = read_dof_vec_master(dv_type, NULL, mesh, fe_space, true /* next? */);
1534
1535 if (dv == NULL) {
1536 return NULL;
1537 }
1538
1539 fe_space = (FE_SPACE *)dv->fe_space;
1540
1541 /* more to read? */
1542 CHAIN_FOREACH(fe_space, dv->fe_space, FE_SPACE) {
1543 dv_chain = CHAIN_NEXT(dv, DOF_VEC);
1544 read_dof_vec_master(dv_type, dv_chain, mesh, fe_space,
1545 CHAIN_NEXT(fe_space, FE_SPACE) == dv->fe_space
1546 ? false : true);
1547 }
1548
1549 if (xdrp) {
1550 AI_xdr_close(xdrp);
1551 }
1552
1553 xdrp = NULL;
1554 file = NULL;
1555
1556 return dv;
1557 }
1558
1559 static inline
file_read_dof_vec_master(bool do_xdr,const char * fn,MESH * mesh,FE_SPACE * fe_space,DV_TYPE_ENUM dv_type)1560 DOF_VEC *file_read_dof_vec_master(bool do_xdr,
1561 const char *fn,
1562 MESH *mesh, FE_SPACE *fe_space,
1563 DV_TYPE_ENUM dv_type)
1564 {
1565 FUNCNAME("file_read_dof_vec_master");
1566 FILE *fp;
1567 DOF_VEC *dv;
1568
1569 if (!(fp = fopen(fn, "rb"))) {
1570 ERROR("Cannot open file '%s'\n",fn);
1571 return NULL;
1572 }
1573
1574 dv = fread_dof_vec_master(do_xdr, fp, mesh, fe_space, dv_type);
1575
1576 fclose(fp);
1577
1578 if (dv) {
1579 MSG("File '%s' read.\n",fn);
1580 }
1581
1582 return dv;
1583 }
1584
1585 /******************************************************************************/
1586
1587 #define DEFUN_READ_DOF_VEC(TYPE, type) \
1588 DOF_##TYPE * \
1589 read_dof_##type##_xdr(const char *fn, MESH *mesh, FE_SPACE *fe_space) \
1590 { \
1591 return (DOF_##TYPE *) \
1592 file_read_dof_vec_master(true, fn, mesh, fe_space, __DOF_##TYPE); \
1593 } \
1594 DOF_##TYPE * \
1595 read_dof_##type(const char *fn, MESH *mesh, FE_SPACE *fe_space) \
1596 { \
1597 return (DOF_##TYPE *) \
1598 file_read_dof_vec_master(false, fn, mesh, fe_space, __DOF_##TYPE); \
1599 } \
1600 DOF_##TYPE * \
1601 fread_dof_##type##_xdr(FILE *fp, MESH *mesh, FE_SPACE *fe_space) \
1602 { \
1603 return (DOF_##TYPE *) \
1604 fread_dof_vec_master(true, fp, mesh, fe_space, __DOF_##TYPE); \
1605 } \
1606 DOF_##TYPE * \
1607 fread_dof_##type(FILE *fp, MESH *mesh, FE_SPACE *fe_space) \
1608 { \
1609 return (DOF_##TYPE *) \
1610 fread_dof_vec_master(false, fp, mesh, fe_space, __DOF_##TYPE); \
1611 } \
1612 struct _AI_semicolon_dummy
1613
1614 DEFUN_READ_DOF_VEC(REAL_VEC, real_vec);
1615 DEFUN_READ_DOF_VEC(REAL_D_VEC, real_d_vec);
1616 DEFUN_READ_DOF_VEC(REAL_VEC_D, real_vec_d);
1617 DEFUN_READ_DOF_VEC(INT_VEC, int_vec);
1618 DEFUN_READ_DOF_VEC(SCHAR_VEC, schar_vec);
1619 DEFUN_READ_DOF_VEC(UCHAR_VEC, uchar_vec);
1620
1621 /*--------------------------------------------------------------------------*/
1622
1623 #if 0
1624 /* FORM
1625 len
1626 AFEM
1627 AMSH
1628 len
1629 data
1630 AVEC
1631 len
1632 data
1633 AVEC
1634 len
1635 data
1636
1637 where AVEC is one of DRV DRDV DUCV DSCV DINV
1638 */
1639
1640 /* Sequentially suck in an IFF file, store the data type in tag, and
1641 * return the size inf SIZE. When positioned at FORM, skip to ADTA,
1642 * when positioned at ADTA, skip to the first data record and read it
1643 * its tag into TAG and its size into SIZE. When positioned at a
1644 * data-record, read its tag into TAG and its size into SIZE.
1645 *
1646 * If "true" is returned and *size == -1 then we have an IO error. If
1647 * "true" is returned and *size >= 0 then it is no ALBERTA tag, but
1648 * still this may be a valid IFF-tag. If "false" is returned we have
1649 * hit one of the known tags and positioned at an ALBERTA data-record.
1650 */
1651 bool fread_iff(FILE *fp, char tag[4], unsigned int *size)
1652 {
1653 *size = -1;
1654
1655 if (fread(tag, 1, 4, fp) != 4) {
1656 return true;
1657 }
1658 if (is_iff_tag(tag, "FORM")) {
1659 if (fread(size, 4, 1, fp) != 1) {
1660 *size = -1;
1661 return true;
1662 }
1663 *size = ntohl(*size);
1664 if (fread(tag, 1, 4, fp) != 4) {
1665 *size = -1;
1666 return true;
1667 }
1668 if (!is_iff_tag(tag, IFF_TAG_ALBERTA)) {
1669 *size -= 4; /* amout to skip */
1670 return true;
1671 }
1672 if (fread(tag, 1, 4, fp) != 4) {
1673 *size = -1;
1674 return true;
1675 }
1676 }
1677 if (fread(size, 4, 1, fp) != 1) {
1678 *size = -1;
1679 return true;
1680 }
1681 *size = ntohl(*size);
1682
1683 if (is_iff_tag(tag, IFF_TAG_MESH) ||
1684 is_iff_tag(tag, IFF_TAG_REAL_VEC) ||
1685 is_iff_tag(tag, IFF_TAG_REAL_D_VEC) ||
1686 is_iff_tag(tag, IFF_TAG_INT_VEC) ||
1687 is_iff_tag(tag, IFF_TAG_UCHAR_VEC) ||
1688 is_iff_tag(tag, IFF_TAG_SCHAR_VEC)) {
1689
1690 return false; /* !!success */
1691 }
1692 return true;
1693 }
1694
1695 FILE *read_iff(const char *filename, char tag[4], unsigned int *size)
1696 {
1697 FILE *fp;
1698
1699 if ((fp = fopen(filename, "rb")) == NULL) {
1700 ERROR("Cannot open file '%s' for reading.\n", filename);
1701 return NULL;
1702 }
1703 if (fread_iff(fp, tag, size) == true) {
1704 fclose(fp);
1705 return NULL;
1706 }
1707 return fp;
1708 }
1709
1710 #endif
1711