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:     write_mesh.c                                                   */
7 /*                                                                          */
8 /* description:  functions for writing meshes and DOF vectors in binary     */
9 /*               machine independent and native 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 #include "alberta_intern.h"
33 #include "alberta.h"
34 
35 static XDR  *xdrp;
36 static FILE *file;
37 
38 /*
39    WARNING:
40    XDR routines to read/write ALBERTA types must be changed if the
41    ALBERTA types change!
42 
43    current state:  REAL   = double
44                    U_CHAR = unsigned char
45                    S_CHAR = signed char
46                    DOF    = int
47 
48    Another WARNING! (D.K.)
49    XDR routines are not well documented in the "xdr" man page.
50    Do not change anything unless you know what you are doing!
51 */
52 
53 /*--------------------------------------------------------------------------*/
54 
55 typedef DOF_REAL_VEC DOF_VEC;
56 
57 static bool
58 write_dof_vec_master(const DOF_VEC *dv, const char *dofvectype,
59 		     const char term[]);
60 
61 #define DOF_REAL_VEC_NAME   "DOF_REAL_VEC    "
62 #define DOF_REAL_D_VEC_NAME "DOF_REAL_D_VEC  "
63 #define DOF_INT_VEC_NAME    "DOF_INT_VEC     "
64 #define DOF_SCHAR_VEC_NAME  "DOF_SCHAR_VEC   "
65 #define DOF_UCHAR_VEC_NAME  "DOF_UCHAR_VEC   "
66 
AI_fwrite(const void * ptr,size_t size,size_t nmemb,FILE * stream)67 size_t AI_fwrite(const void *ptr, size_t size, size_t nmemb, FILE *stream)
68 {
69   size_t nwritten;
70 
71   nwritten = fwrite(ptr, size, nmemb, stream);
72 
73   return nwritten;
74 }
75 
write_int(int val)76 static void write_int(int val)
77 {
78   bool result;
79 
80   if(xdrp)
81     result = xdr_int(xdrp, &val);
82   else
83     result = fwrite(&val, sizeof(int), 1, file) == 1;
84   /* Now, what to do with RESULT? ;) */
85   (void)result;
86 }
87 
88 #if 0
89 /* NEVER use this. If you need 64 bits, then use xdr_int64_t() */
90 static void write_long(long int val)
91 {
92   bool result;
93 
94   if(xdrp)
95     result = xdr_long(xdrp, &val);
96   else
97     result = fwrite(&val, sizeof(long int), 1, file) == 1;
98   /* Now, what to do with RESULT? ;) */
99 }
100 #endif
101 
write_REAL(REAL val)102 static void write_REAL(REAL val)
103 {
104   bool result;
105 
106   if(xdrp)
107     result = AI_xdr_REAL(xdrp, &val);
108   else
109     result = fwrite(&val, sizeof(REAL), 1, file) == 1;
110   /* Now, what to do with RESULT? ;) */
111   (void)result;
112 }
113 
write_string(const char * string,int write_length)114 static void write_string(const char *string, int write_length)
115 {
116   bool result;
117   int strileng = 0;
118 
119   if(string)
120     strileng = strlen(string);
121 
122   if(write_length)
123     write_int(strileng);
124 
125   if(strileng) {
126     if(xdrp)
127       result = xdr_string(xdrp, (char **)&string, strileng+1);
128     else
129       result = fwrite(string, sizeof(char), strileng+1, file) == (strileng+1);
130   }
131   /* Now, what to do with RESULT? ;) */
132   (void)result;
133 }
134 
write_vector(void * start,int n,size_t size,xdrproc_t xdrproc)135 static void write_vector(void *start, int n, size_t size, xdrproc_t xdrproc)
136 {
137   bool result;
138 
139   if(xdrp)
140     result = xdr_vector(xdrp, (char *)start, n, size, xdrproc);
141   else
142     result = fwrite(start, size, n, file) == n;
143   /* Now, what to do with RESULT? ;) */
144   (void)result;
145 }
146 
write_U_CHAR(U_CHAR val)147 static void write_U_CHAR(U_CHAR val)
148 {
149   bool result;
150 
151   if(xdrp)
152     result = AI_xdr_U_CHAR(xdrp, &val);
153   else
154     result = fwrite(&val, sizeof(U_CHAR), 1, file) == 1;
155   /* Now, what to do with RESULT? ;) */
156   (void)result;
157 }
158 
write_S_CHAR(S_CHAR val)159 static void write_S_CHAR(S_CHAR val)
160 {
161   bool result;
162 
163   if(xdrp)
164     result = AI_xdr_S_CHAR(xdrp, &val);
165   else
166     result = fwrite(&val, sizeof(S_CHAR), 1, file) == 1;
167   /* Now, what to do with RESULT? ;) */
168   (void)result;
169 }
170 
171 /*--------------------------------------------------------------------------*/
172 
write_mesh_master(MESH * mesh,REAL time)173 static bool write_mesh_master(MESH *mesh, REAL time)
174 {
175   FUNCNAME("write_mesh_master");
176   MACRO_EL  *mel, *meln;
177   DOF_ADMIN *admin, *face_admin, *vert_admin, *edge_admin;
178   DOF       *face_dofs = NULL, *vert_dofs = NULL, *edge_dofs = NULL;
179   DOF       **face_ptrs = NULL, **vert_ptrs = NULL, **edge_ptrs = NULL;
180   int       n_face_ptrs, n_vert_ptrs, n_edge_ptrs;
181   int       i, n, iadmin, dim;
182   int       neigh_i[N_NEIGH_MAX];
183   int       wt_i[N_WALLS_MAX];
184 
185   if (!mesh) {
186     ERROR("no mesh - no file created\n");
187     return true;
188   }
189   dim = mesh->dim;
190 
191   dof_compress(mesh);
192 
193   write_string(ALBERTA_VERSION, false); /* file marker */
194 
195   write_int(dim);
196 
197   write_int(DIM_OF_WORLD);
198 
199   write_REAL(time);
200 
201   write_string(mesh->name, true);
202 
203   write_int(mesh->n_vertices);
204 
205   if(dim > 1) {
206     write_int(mesh->n_edges);
207   }
208 
209   write_int(mesh->n_elements);
210   write_int(mesh->n_hier_elements);
211 
212   if(dim == 3) {
213     write_int(mesh->n_faces);
214     write_int(mesh->max_edge_neigh);
215   }
216 
217   write_S_CHAR((S_CHAR)mesh->is_periodic);
218   if (mesh->is_periodic) {
219     write_int(mesh->per_n_vertices);
220 
221     if(dim > 1) {
222       write_int(mesh->per_n_edges);
223       if(dim > 2) {
224 	write_int(mesh->per_n_faces);
225       }
226     }
227     write_int(mesh->n_wall_trafos);
228     for (i = 0; i < mesh->n_wall_trafos; i++) {
229       write_vector((REAL *)mesh->wall_trafos[i],
230 		   SQR(DIM_OF_WORLD) + DIM_OF_WORLD, sizeof(REAL),
231 		   (xdrproc_t)AI_xdr_REAL);
232     }
233   }
234 
235   write_vector(mesh->bbox,
236 	       2*DIM_OF_WORLD, sizeof(REAL), (xdrproc_t)AI_xdr_REAL);
237 
238   write_int(mesh->n_dof_el);
239   write_vector(mesh->n_dof, N_NODE_TYPES, sizeof(int), (xdrproc_t)xdr_int);
240   write_int(mesh->n_node_el);
241   write_vector(mesh->node, N_NODE_TYPES, sizeof(int), (xdrproc_t)xdr_int);
242 
243 
244 
245   write_int(mesh->n_dof_admin);
246   for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++) {
247     admin = mesh->dof_admin[iadmin];
248 
249     DEBUG_TEST_EXIT(admin, "dof admin no. %d not found!\n", iadmin);
250 
251     write_vector(admin->n_dof, N_NODE_TYPES, sizeof(int), (xdrproc_t)xdr_int);
252     write_int(admin->used_count);
253 
254     /* after dof_compress(), no more information is required */
255 
256     write_string(admin->name, true);
257 
258     write_U_CHAR(admin->flags & 0xff); /* preserve coarse dofs and periodic */
259   }
260 
261 
262   face_admin = NULL;
263   if (mesh->n_dof[FACE] > 0) {
264     for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++) {
265       if ((admin = mesh->dof_admin[iadmin]) != NULL) {
266 	if (admin->n_dof[FACE] > 0 &&
267 	    (face_admin == NULL ||
268 	     ((admin->flags & ADM_PRESERVE_COARSE_DOFS) != 0 &&
269 	      (face_admin->flags & ADM_PRESERVE_COARSE_DOFS) == 0)
270 	     ||
271 	     (((admin->flags & ADM_PRESERVE_COARSE_DOFS)
272 	       ==
273 	       (face_admin->flags & ADM_PRESERVE_COARSE_DOFS) &&
274 	       admin->n_dof[FACE] < face_admin->n_dof[FACE])))) {
275 	  face_admin = admin;
276 	}
277       }
278     }
279     DEBUG_TEST_EXIT(face_admin, "no admin with face dofs?\n");
280   }
281 
282   edge_admin = NULL;
283   if (mesh->n_dof[EDGE] > 0) {
284     if (face_admin && (face_admin->n_dof[EDGE] > 0)
285 	&& (face_admin->flags & ADM_PRESERVE_COARSE_DOFS)) {
286       edge_admin = face_admin;
287     } else {
288       for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++) {
289 	if ((admin = mesh->dof_admin[iadmin]) != NULL) {
290 	  if (admin->n_dof[EDGE] > 0 &&
291 	      (edge_admin == NULL ||
292 	       ((admin->flags & ADM_PRESERVE_COARSE_DOFS) != 0 &&
293 		(edge_admin->flags & ADM_PRESERVE_COARSE_DOFS) == 0)
294 	       ||
295 	       ((admin->flags & ADM_PRESERVE_COARSE_DOFS)
296 		==
297 		(edge_admin->flags & ADM_PRESERVE_COARSE_DOFS) &&
298 		admin->n_dof[EDGE] < edge_admin->n_dof[EDGE]))) {
299 	    edge_admin = admin;
300 	  }
301 	}
302       }
303       DEBUG_TEST_EXIT(edge_admin, "no admin with edge dofs?\n");
304     }
305   }
306 
307   vert_admin = NULL;
308   if (mesh->n_dof[VERTEX] > 0) {
309     if (face_admin && (face_admin->n_dof[VERTEX] > 0)) {
310       vert_admin = face_admin;
311     } else if (edge_admin && (edge_admin->n_dof[VERTEX] > 0)) {
312       vert_admin = edge_admin;
313     } else {
314       for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++) {
315 	if ((admin = mesh->dof_admin[iadmin])) {
316 	  if (admin->n_dof[VERTEX] > 0  &&
317 	      (vert_admin == NULL ||
318 	       admin->n_dof[VERTEX] < vert_admin->n_dof[VERTEX])) {
319 	    vert_admin = admin;
320 	  }
321 	}
322       }
323       DEBUG_TEST_EXIT(vert_admin, "no admin with vertex dofs?\n");
324     }
325   }
326 
327   /* On a periodic mesh we have duplicated DOF-pointers across
328    * periodic boundaries. So we first loop over the mesh and count the
329    * number of DOF pointers for each type of node (FACE/EDGE/VERTEX)
330    * and then loop again to do the actual dumping.
331    */
332   n_vert_ptrs = n_edge_ptrs = n_face_ptrs = 0;
333   TRAVERSE_FIRST(mesh, -1, CALL_EVERY_EL_PREORDER|FILL_NOTHING) {
334     EL *el = el_info->el;
335     int n0, node;
336     DOF dof;
337 
338     if (vert_admin) { /* shouldn't there be one, always? */
339       node = mesh->node[VERTEX];
340       n0   = vert_admin->n0_dof[VERTEX];
341       for (i = 0; i < N_VERTICES(dim); i++) {
342 	if ((dof = el->dof[node+i][n0]) >= DOF_UNUSED) {
343 	  el->dof[node+i][n0] = DOF_UNUSED - 1 -  dof;
344 	  n_vert_ptrs++;
345 	}
346       }
347     }
348 
349     if (edge_admin) {
350       node = mesh->node[EDGE];
351       n0   = edge_admin->n0_dof[EDGE];
352       for (i = 0; i < N_EDGES(dim); i++) {
353 	if ((dof = el->dof[node+i][n0]) >= DOF_UNUSED) {
354 	  el->dof[node+i][n0] = DOF_UNUSED - 1 -  dof;
355 	  n_edge_ptrs++;
356 	}
357       }
358     }
359 
360     if (face_admin) {
361       node = mesh->node[FACE];
362       n0   = face_admin->n0_dof[FACE];
363       for (i = 0; i < N_FACES(dim); i++) {
364 	if ((dof = el->dof[node+i][n0]) >= DOF_UNUSED) {
365 	  el->dof[node+i][n0] = DOF_UNUSED - 1 -  dof;
366 	  n_face_ptrs++;
367 	}
368       }
369     }
370 
371   } TRAVERSE_NEXT();
372 
373   /* Loop again over the mesh and enumerate the DOF-pointers, the
374    * original DOF values will we saved and restored later when dumping
375    * the elements to disk.
376    */
377   if (vert_admin) {
378     vert_ptrs = MEM_ALLOC(n_vert_ptrs, DOF *);
379     vert_dofs = MEM_ALLOC(n_vert_ptrs, DOF);
380   }
381   if (edge_admin) {
382     edge_ptrs = MEM_ALLOC(n_edge_ptrs, DOF *);
383     edge_dofs = MEM_ALLOC(n_edge_ptrs, DOF);
384   }
385   if (face_admin) {
386     face_ptrs = MEM_ALLOC(n_face_ptrs, DOF *);
387     face_dofs = MEM_ALLOC(n_face_ptrs, DOF);
388   }
389 
390   n_vert_ptrs = n_edge_ptrs = n_face_ptrs = 0;
391   TRAVERSE_FIRST(mesh, -1, CALL_EVERY_EL_PREORDER|FILL_NOTHING) {
392     DOF *dofptr;
393     EL *el = el_info->el;
394     int n0, node;
395 
396     if (vert_admin) { /* shouldn't there be one, always? */
397       node = mesh->node[VERTEX];
398       n0   = vert_admin->n0_dof[VERTEX];
399       for (i = 0; i < N_VERTICES(dim); i++) {
400 	if ((dofptr = el->dof[node+i])[n0] < DOF_UNUSED) {
401 	  dofptr[n0]               = DOF_UNUSED - 1 - dofptr[n0];
402 	  vert_ptrs[n_vert_ptrs++] = dofptr;
403 	}
404       }
405     }
406 
407     if (edge_admin) {
408       node = mesh->node[EDGE];
409       n0   = edge_admin->n0_dof[EDGE];
410       for (i = 0; i < N_EDGES(dim); i++) {
411 	if ((dofptr = el->dof[node+i])[n0] < DOF_UNUSED) {
412 	  dofptr[n0]               = DOF_UNUSED - 1 - dofptr[n0];
413 	  edge_ptrs[n_edge_ptrs++] = dofptr;
414 	}
415       }
416     }
417 
418     if (face_admin) {
419       node = mesh->node[FACE];
420       n0   = face_admin->n0_dof[FACE];
421       for (i = 0; i < N_FACES(dim); i++) {
422 	if ((dofptr = el->dof[node+i])[n0] < DOF_UNUSED) {
423 	  dofptr[n0]               = DOF_UNUSED - 1 - dofptr[n0];
424 	  face_ptrs[n_face_ptrs++] = dofptr;
425 	}
426       }
427     }
428 
429   } TRAVERSE_NEXT();
430 
431   write_int(n_vert_ptrs);
432   if (n_vert_ptrs) {
433     int n, n0;
434     n  = mesh->n_dof[VERTEX];
435     n0 = vert_admin->n0_dof[VERTEX];
436     for (i = 0; i < n_vert_ptrs; i++) {
437       write_vector(vert_ptrs[i], n, sizeof(DOF), (xdrproc_t)AI_xdr_DOF);
438       vert_dofs[i]     = vert_ptrs[i][n0]; /* undo this abuse later */
439       vert_ptrs[i][n0] = i;
440     }
441   }
442 
443   if(dim > 1) {
444     int n, n0;
445     write_int(n_edge_ptrs);
446     if (n_edge_ptrs) {
447       n  = mesh->n_dof[EDGE];
448       n0 = edge_admin->n0_dof[EDGE];
449       for (i = 0; i < n_edge_ptrs; i++) {
450 	write_vector(edge_ptrs[i], n, sizeof(DOF), (xdrproc_t)AI_xdr_DOF);
451 	edge_dofs[i]     = edge_ptrs[i][n0]; /* undo this abuse later */
452 	edge_ptrs[i][n0] = i;
453       }
454     }
455   }
456 
457   if(dim == 3) {
458     int n, n0;
459     write_int(n_face_ptrs);
460     if (n_face_ptrs) {
461       n  = mesh->n_dof[FACE];
462       n0 = face_admin->n0_dof[FACE];
463       for (i = 0; i < n_face_ptrs; i++) {
464 	write_vector(face_ptrs[i], n, sizeof(DOF), (xdrproc_t)AI_xdr_DOF);
465 	face_dofs[i]     = face_ptrs[i][n0]; /* undo this abuse later */
466 	face_ptrs[i][n0] = i;
467       }
468     }
469   }
470 
471   /*
472    * gather info about macro elements (vertices, ...)
473    */
474   {
475     typedef int  intNV[N_VERTICES_MAX];
476     intNV    *mcindex = MEM_ALLOC(mesh->n_macro_el, intNV);
477     REAL_D   *mccoord;
478     int    mccount, m;
479 
480     mccount = ((MESH_MEM_INFO *)(mesh->mem_info))->count;
481     mccoord  = ((MESH_MEM_INFO *)(mesh->mem_info))->coords;
482 
483     for (m = 0, mel = mesh->macro_els;
484 	 m < mesh->n_macro_el;
485 	 m++, mel = mesh->macro_els + m) {
486       for (i = 0; i < N_VERTICES(dim); i++) {
487 	mcindex[m][i] =
488 	  ((char *)(mel->coord[i])
489 	   -
490 	   (char *)mccoord)/(sizeof(char)*sizeof(REAL_D));
491       }
492       if (mel->index != m) {
493 	mel->index = m;
494       }
495     }
496 
497     write_int(mesh->n_macro_el);
498     write_int(mccount);                           /* number of macro coords */
499 
500     for (i = 0; i < mccount; i++) {
501       write_vector(mccoord[i], DIM_OF_WORLD, sizeof(REAL),
502 		   (xdrproc_t)AI_xdr_REAL);
503     }
504 
505     TRAVERSE_FIRST(mesh, -1, CALL_EVERY_EL_PREORDER|FILL_NOTHING) {
506       EL *el = el_info->el;
507       int node, n0;
508 
509       if (el_info->level == 0) {
510 	/* this is a macro element */
511 	MACRO_EL *mel = (MACRO_EL *)el_info->macro_el;
512 
513 	m = mel->index;
514 
515 	write_vector(mcindex[m], N_VERTICES(dim), sizeof(int),
516 		     (xdrproc_t)xdr_int);
517 	write_vector(mel->wall_bound, N_WALLS(dim), sizeof(BNDRY_TYPE),
518 		     (xdrproc_t)AI_xdr_S_CHAR);
519 
520 	for (i = 0; i < N_NEIGH(dim); i++) {
521 	  if ((meln = mel->neigh[i])) {
522 	    neigh_i[i] = meln->index;
523 	  } else {
524 	    neigh_i[i] = -1;
525 	  }
526 	}
527 
528 	write_vector(neigh_i, N_NEIGH(dim), sizeof(int), (xdrproc_t)xdr_int);
529 	write_vector(mel->opp_vertex, N_NEIGH(dim), sizeof(S_CHAR),
530 		     (xdrproc_t)AI_xdr_S_CHAR);
531 #if DIM_MAX > 2
532 	if(dim == 3) {
533 	  write_U_CHAR(mel->el_type);
534 	}
535 #endif
536 
537 	/* Dump the projection status of the element/faces */
538 	for (i = 0; i < 1+N_WALLS(dim); i++) {
539 	  write_S_CHAR(mel->projection[i] != NULL);
540 	}
541 
542 	/* dump the binding to the master mesh, if this is a trace mesh */
543 	if (mel->master.macro_el) {
544 	  write_int(mel->master.macro_el->index);
545 	  write_S_CHAR(mel->master.opp_vertex);
546 	} else {
547 	  write_int(-1);
548 	  write_S_CHAR(-1);
549 	}
550 
551 	/* dump the status of periodic boundaries */
552 
553 	if (mesh->is_periodic) {
554 	  for (i = 0; i < N_NEIGH(dim); i++) {
555 	    write_vector(mel->neigh_vertices[i], N_VERTICES(dim-1),
556 			 sizeof(S_CHAR), (xdrproc_t)AI_xdr_S_CHAR);
557 	  }
558 	  for (i = 0; i < N_WALLS(dim); i++) {
559 	    AFF_TRAFO *wt;
560 	    int j;
561 	    if ((wt = mel->wall_trafo[i]) != NULL) {
562 	      for (j = 0;
563 		   j < mesh->n_wall_trafos && mesh->wall_trafos[j] != wt; j++);
564 	      DEBUG_TEST_EXIT(j < mesh->n_wall_trafos,
565 			      "Iconsistent or unknown wall transformations.\n");
566 	      wt_i[i] = j;
567 	    } else {
568 	      wt_i[i] = -1;
569 	    }
570 	  }
571 	  write_vector(wt_i, N_WALLS(dim), sizeof(int), (xdrproc_t)xdr_int);
572 	}
573       }
574 
575       /* dump the element to disk (macro or not) */
576 
577       if (el->child[0]) {
578 	DEBUG_TEST_EXIT(el->child[1], "child 0 but no child 1\n");
579 	write_U_CHAR(true);
580       } else {
581 	write_U_CHAR(false);
582       }
583 
584       if(dim > 1) {
585 	if (el->new_coord) {
586 	  write_U_CHAR(true);
587 	  write_vector(el->new_coord, DIM_OF_WORLD, sizeof(REAL),
588 		       (xdrproc_t)AI_xdr_REAL);
589 	} else {
590 	  write_U_CHAR(false);
591 	}
592       }
593 
594       if (mesh->n_dof[VERTEX] > 0) {
595 	node = mesh->node[VERTEX];
596 	n0   = vert_admin->n0_dof[VERTEX];
597 	for (i = 0; i < N_VERTICES(dim); i++) {
598 	  write_int(el->dof[node + i][n0]);
599 	}
600       }
601 
602       if (dim > 1 && mesh->n_dof[EDGE] > 0) {
603 	node = mesh->node[EDGE];
604 	n0   = edge_admin->n0_dof[EDGE];
605 
606 	for (i = 0; i < N_EDGES(dim); i++) {
607 	  if(el->dof[node + i][n0] > DOF_UNUSED) {
608 	    write_int(el->dof[node + i][n0]);
609 	  } else {
610 	    write_int(DOF_UNUSED);
611 	  }
612 	}
613       }
614 
615       if (dim == 3 && mesh->n_dof[FACE] > 0) {
616 	node = mesh->node[FACE];
617 	n0   = face_admin->n0_dof[FACE];
618 
619 	for (i = 0; i < N_FACES_3D; i++) {
620 	  if(el->dof[node + i][n0] > DOF_UNUSED) {
621 	    write_int(el->dof[node + i][n0]);
622 	  } else {
623 	    write_int(DOF_UNUSED);
624 	  }
625 	}
626       }
627 
628       if ((n = mesh->n_dof[CENTER]) > 0) {
629 	node = mesh->node[CENTER];
630 
631 	write_vector(el->dof[node], n, sizeof(DOF), (xdrproc_t)AI_xdr_DOF);
632       }
633 
634     } TRAVERSE_NEXT();
635 
636     MEM_FREE(mcindex, mesh->n_macro_el, intNV);
637   }
638 
639   /* Write the magic cookie. */
640   write_int(mesh->cookie);
641 
642   if (dim == 3 && face_admin) {
643     /* It remains to undo the abuse of el->dof[][n0] */
644     int n0 = face_admin->n0_dof[FACE];
645     for (i = 0; i < n_face_ptrs; i++) {
646       face_ptrs[i][n0] = face_dofs[i];
647     }
648     MEM_FREE(face_ptrs, n_face_ptrs, DOF *);
649     MEM_FREE(face_dofs, n_face_ptrs, DOF);
650   }
651 
652   if (dim > 1 && edge_admin) {
653     /* It remains to undo the abuse of el->dof[][n0] */
654     int n0 = edge_admin->n0_dof[EDGE];
655     for (i = 0; i < n_edge_ptrs; i++) {
656       edge_ptrs[i][n0] = edge_dofs[i];
657     }
658     MEM_FREE(edge_ptrs, n_edge_ptrs, DOF *);
659     MEM_FREE(edge_dofs, n_edge_ptrs, DOF);
660   }
661 
662   if (vert_admin) {
663     /* It remains to undo the abuse of el->dof[][n0] */
664     int n0 = vert_admin->n0_dof[VERTEX];
665     for (i = 0; i < n_vert_ptrs; i++) {
666       vert_ptrs[i][n0] = vert_dofs[i];
667     }
668     MEM_FREE(vert_ptrs, n_vert_ptrs, DOF *);
669     MEM_FREE(vert_dofs, n_vert_ptrs, DOF);
670   }
671 
672   /* If this is a parametric mesh with the internal standard Lagrange
673    * parameterization, then we dump the co-ordinate and touched-edges
674    * information to disk here.
675    */
676   if (_AI_is_lagrange_parametric(mesh)) {
677     PARAMETRIC *parametric = mesh->parametric;
678     DOF_PTR_VEC *edge_projections;
679 
680     write_string("LAPA", false);
681     write_S_CHAR(parametric->not_all);
682     write_S_CHAR(parametric->use_reference_mesh);
683     write_int(_AI_lagrange_strategy(mesh));
684     write_dof_vec_master((DOF_VEC *)get_lagrange_coords(mesh),
685 			 DOF_REAL_D_VEC_NAME, "EOF.");
686     edge_projections = get_lagrange_edge_projections(mesh);
687     if (edge_projections != NULL) {
688       DOF_UCHAR_VEC *touched_edges =
689         get_dof_uchar_vec("touched edges",
690                           edge_projections->fe_space);
691       FOR_ALL_DOFS(edge_projections->fe_space->admin, {
692           touched_edges->vec[dof] = edge_projections->vec[dof] != NULL;
693         });
694       write_dof_vec_master((DOF_VEC *)touched_edges,
695                            DOF_UCHAR_VEC_NAME, "EOF.");
696       free_dof_uchar_vec(touched_edges);
697     }
698   }
699 
700   {
701     MESH_MEM_INFO *minfo = (MESH_MEM_INFO *)mesh->mem_info;
702 
703     /* If we have trace-meshes, then dump those to disk as well. */
704     for (i = 0; i < minfo->n_slaves; i++) {
705       write_string("TMSH", false);
706       write_int(minfo->slaves[i]->trace_id);
707       write_mesh_master(minfo->slaves[i], time);
708     }
709   }
710 
711   write_string("EOF.", false);                       /* file end marker */
712 
713   return false; /* Success status. */
714 }
715 
fwrite_mesh(MESH * mesh,FILE * fp,REAL time)716 bool fwrite_mesh(MESH *mesh, FILE *fp, REAL time)
717 {
718   bool result;
719 
720   file = fp;
721 
722   result = write_mesh_master(mesh, time);
723 
724   file = NULL;
725 
726   return result;
727 }
728 
write_mesh(MESH * mesh,const char * filename,REAL time)729 bool write_mesh(MESH *mesh, const char *filename, REAL time)
730 {
731   FUNCNAME("write_mesh");
732   FILE *fp;
733   bool result;
734 
735   if ((fp = fopen(filename,"wb")) == NULL) {
736     ERROR("Cannot open file '%s' for writing.\n", filename);
737     return true;
738   }
739 
740   result = fwrite_mesh(mesh, fp, time);
741 
742   fclose(fp);
743 
744   return result;
745 }
746 
fwrite_mesh_xdr(MESH * mesh,FILE * fp,REAL time)747 bool fwrite_mesh_xdr(MESH *mesh, FILE *fp, REAL time)
748 {
749   FUNCNAME("fwrite_mesh_xdr");
750   bool result;
751 
752   if (!(xdrp = AI_xdr_fopen(fp, XDR_ENCODE))) {
753     ERROR("Cannot convert file handle to XDR handle.\n");
754     return true;
755   }
756   file = fp;
757 
758   result = write_mesh_master(mesh, time);
759 
760   AI_xdr_close(xdrp);
761 
762   xdrp = NULL;
763   file = NULL;
764 
765   return result;
766 }
767 
write_mesh_xdr(MESH * mesh,const char * filename,REAL time)768 bool write_mesh_xdr(MESH *mesh, const char *filename, REAL time)
769 {
770   FUNCNAME("write_mesh_xdr");
771   FILE *fp;
772   bool result;
773 
774   if ((fp = fopen(filename,"wb")) == NULL) {
775     ERROR("Cannot open file '%s' for writing.\n", filename);
776     return true;
777   }
778 
779   result = fwrite_mesh_xdr(mesh, fp, time);
780 
781   fclose(fp);
782 
783   return result;
784 }
785 
786 /*--------------------------------------------------------------------------*/
787 /* write DOF vectors of various types                                       */
788 /*--------------------------------------------------------------------------*/
789 static bool
write_dof_vec_master(const DOF_VEC * dv,const char * dofvectype,const char term[])790 write_dof_vec_master(const DOF_VEC *dv, const char *dofvectype,
791 		     const char term[])
792 {
793   FUNCNAME("write_dof_vec_master");
794   U_CHAR          flags_and_stride;
795   const FE_SPACE  *fe_space;
796   const DOF_ADMIN *admin;
797   MESH            *mesh;
798   int              i, iadmin, last;
799 
800   if (!dv || !(fe_space = dv->fe_space)) {
801     ERROR("no %s or fe_space - no file created\n", dofvectype);
802     return(1);
803   }
804   if (!(admin = fe_space->admin) || !(mesh = admin->mesh)) {
805     ERROR("no dof_admin or dof_admin->mesh - no file created\n");
806     return true; /* crazy return convention: 1 means error, 0 means no error */
807   }
808 
809   dof_compress(mesh);
810 
811   iadmin = -1;
812   for (i = 0; i < mesh->n_dof_admin; i++) {
813     if (mesh->dof_admin[i] == admin) {
814       iadmin = i;
815       break;
816     }
817   }
818   if (iadmin < 0) {
819     ERROR("vec->admin not in mesh->dof_admin[] - no file created\n");
820     return true;
821   }
822 
823   last = admin->used_count;
824   DEBUG_TEST_EXIT(last <= dv->size,
825 		  "dof_vec->size %d < admin->size_used %d\n", dv->size, last);
826 
827   write_string(dofvectype, false);
828 
829   write_string(dv->name, true);
830 
831 #if (ADM_FLAGS_MASK & 0x80)
832 # error ADM_FLAGS_MASK overlaps the sign bit.
833 #endif
834 
835   /* We (mis-)use the sign bit to code the stride of DOF_REAL_VEC_Ds. */
836   flags_and_stride = admin->flags & ADM_FLAGS_MASK;
837   if (dv->reserved != 1) {
838     flags_and_stride |= 0x80;
839   }
840   write_U_CHAR(flags_and_stride);
841 
842   write_vector((void *)admin->n_dof, N_NODE_TYPES, sizeof(int),
843 	       (xdrproc_t)xdr_int);
844 
845   if (fe_space->bas_fcts)
846     write_string(fe_space->bas_fcts->name, true);
847   else
848     write_int(0);
849 
850 
851   write_int(last);
852 
853   if (last) {
854     if (!strncmp(dofvectype, "DOF_REAL_VEC    ", 12))
855       write_vector(dv->vec, last, sizeof(REAL), (xdrproc_t)AI_xdr_REAL);
856     else if (!strncmp(dofvectype, "DOF_REAL_D_VEC  ", 12))
857       write_vector(dv->vec, last*DIM_OF_WORLD, sizeof(REAL),
858 		   (xdrproc_t)AI_xdr_REAL);
859     else if (!strncmp(dofvectype, "DOF_INT_VEC     ", 12))
860       write_vector(dv->vec, last, sizeof(int), (xdrproc_t)xdr_int);
861     else if (!strncmp(dofvectype, "DOF_SCHAR_VEC   ", 12))
862       write_vector(dv->vec, last, sizeof(S_CHAR), (xdrproc_t)AI_xdr_S_CHAR);
863     else if (!strncmp(dofvectype, "DOF_UCHAR_VEC   ", 12))
864       write_vector(dv->vec, last, sizeof(U_CHAR), (xdrproc_t)AI_xdr_U_CHAR);
865     else
866       ERROR("Invalid file id '%s'.\n",dofvectype);
867   }
868 
869   /* Write the magic cookie. */
870   write_int(mesh->cookie);
871 
872   write_string(term, false);                         /* file end marker */
873 
874   return false; /* no error */
875 }
876 
877 static inline
fwrite_dof_vec_master(bool write_xdr,const DOF_VEC * dv,FILE * fp,const char * dofvectype)878 bool fwrite_dof_vec_master(bool write_xdr, const DOF_VEC *dv,
879 			   FILE *fp, const char *dofvectype)
880 {
881   FUNCNAME("fwrite_dof_vec_master");
882   const DOF_VEC *chain;
883   bool error, is_vec_d;
884 
885   if (write_xdr) {
886     if (!(xdrp = AI_xdr_fopen(fp, XDR_ENCODE))) {
887       ERROR("Cannot convert file handle to XDR handle.\n");
888       return true;
889     }
890   }
891   file = fp;
892 
893   is_vec_d = strcmp("DOF_REAL_VEC_D  ", dofvectype) == 0;
894 
895   chain = dv;
896   CHAIN_DO(chain, const DOF_VEC) {
897     if (is_vec_d) {
898       dofvectype = ((DOF_REAL_VEC_D *)chain)->stride == 1
899 	? "DOF_REAL_VEC    " : "DOF_REAL_D_VEC  ";
900     }
901     error =
902       write_dof_vec_master(chain, dofvectype,
903 			   CHAIN_NEXT(chain, DOF_VEC) == dv ? "EOF." : "NEXT");
904     if (error) {
905       break;
906     }
907   } CHAIN_WHILE(chain, const DOF_VEC);
908 
909   if(write_xdr) {
910     AI_xdr_close(xdrp);
911     xdrp = NULL;
912   }
913   file = NULL;
914 
915   return error;
916 }
917 
file_write_dof_vec_master(bool write_xdr,const DOF_VEC * dv,const char * fn,const char * dofvectype)918 static bool file_write_dof_vec_master(bool write_xdr, const DOF_VEC *dv,
919 				      const char *fn, const char *dofvectype)
920 {
921   FUNCNAME("file_write_dof_vec_master");
922   FILE *fp;
923   bool error;
924 
925   if ((fp = fopen(fn, "wb")) == NULL) {
926     ERROR("Cannot open file '%s' for writing.\n", fn);
927     return true;
928   }
929 
930   error = fwrite_dof_vec_master(write_xdr, dv, fp, dofvectype);
931 
932   fclose(fp);
933 
934   return error;
935 }
936 
937 /*--------------------------------------------------------------------------*/
938 
write_dof_real_vec_xdr(const DOF_REAL_VEC * dv,const char * fn)939 bool write_dof_real_vec_xdr(const DOF_REAL_VEC *dv, const char *fn)
940 {
941   return file_write_dof_vec_master(true, (const DOF_VEC *)dv,
942 				   fn, "DOF_REAL_VEC    ");
943 }
944 
write_dof_real_vec(const DOF_REAL_VEC * dv,const char * fn)945 bool write_dof_real_vec(const DOF_REAL_VEC *dv, const char *fn)
946 {
947   return file_write_dof_vec_master(false, (const DOF_VEC *)dv,
948 				   fn, "DOF_REAL_VEC    ");
949 }
950 
fwrite_dof_real_vec_xdr(const DOF_REAL_VEC * dv,FILE * fp)951 bool fwrite_dof_real_vec_xdr(const DOF_REAL_VEC *dv, FILE *fp)
952 {
953   return fwrite_dof_vec_master(true, (const DOF_VEC *)dv,
954 			       fp, "DOF_REAL_VEC    ");
955 }
956 
fwrite_dof_real_vec(const DOF_REAL_VEC * dv,FILE * fp)957 bool fwrite_dof_real_vec(const DOF_REAL_VEC *dv, FILE *fp)
958 {
959   return fwrite_dof_vec_master(false, (const DOF_VEC *)dv,
960 			       fp, "DOF_REAL_VEC    ");
961 }
962 
write_dof_real_vec_d_xdr(const DOF_REAL_VEC_D * dv,const char * fn)963 bool write_dof_real_vec_d_xdr(const DOF_REAL_VEC_D *dv, const char *fn)
964 {
965   return file_write_dof_vec_master(
966     true, (const DOF_VEC *)dv, fn, "DOF_REAL_VEC_D  ");
967 }
968 
write_dof_real_vec_d(const DOF_REAL_VEC_D * dv,const char * fn)969 bool write_dof_real_vec_d(const DOF_REAL_VEC_D *dv, const char *fn)
970 {
971   return file_write_dof_vec_master(
972     false, (const DOF_VEC *)dv, fn, "DOF_REAL_VEC_D  ");
973 }
974 
fwrite_dof_real_vec_d_xdr(const DOF_REAL_VEC_D * dv,FILE * fp)975 bool fwrite_dof_real_vec_d_xdr(const DOF_REAL_VEC_D *dv, FILE *fp)
976 {
977   return fwrite_dof_vec_master(
978     true, (const DOF_VEC *)dv, fp, "DOF_REAL_VEC_D  ");
979 }
980 
fwrite_dof_real_vec_d(const DOF_REAL_VEC_D * dv,FILE * fp)981 bool fwrite_dof_real_vec_d(const DOF_REAL_VEC_D *dv, FILE *fp)
982 {
983   return fwrite_dof_vec_master(
984     false, (const DOF_VEC *)dv, fp, "DOF_REAL_VEC_D  ");
985 }
986 
987 
988 
989 
write_dof_real_d_vec_xdr(const DOF_REAL_D_VEC * dv,const char * fn)990 bool write_dof_real_d_vec_xdr(const DOF_REAL_D_VEC *dv, const char *fn)
991 {
992   return file_write_dof_vec_master(true, (const DOF_VEC *)dv,
993 				   fn, "DOF_REAL_D_VEC  ");
994 }
995 
write_dof_real_d_vec(const DOF_REAL_D_VEC * dv,const char * fn)996 bool write_dof_real_d_vec(const DOF_REAL_D_VEC *dv, const char *fn)
997 {
998   return file_write_dof_vec_master(false, (const DOF_VEC *)dv,
999 				   fn, "DOF_REAL_D_VEC  ");
1000 }
1001 
fwrite_dof_real_d_vec_xdr(const DOF_REAL_D_VEC * dv,FILE * fp)1002 bool fwrite_dof_real_d_vec_xdr(const DOF_REAL_D_VEC *dv, FILE *fp)
1003 {
1004   return fwrite_dof_vec_master(true, (const DOF_VEC *)dv,
1005 			       fp, "DOF_REAL_D_VEC  ");
1006 }
1007 
fwrite_dof_real_d_vec(const DOF_REAL_D_VEC * dv,FILE * fp)1008 bool fwrite_dof_real_d_vec(const DOF_REAL_D_VEC *dv, FILE *fp)
1009 {
1010   return fwrite_dof_vec_master(false, (const DOF_VEC *)dv,
1011 			       fp, "DOF_REAL_D_VEC  ");
1012 }
1013 
write_dof_int_vec_xdr(const DOF_INT_VEC * dv,const char * fn)1014 bool write_dof_int_vec_xdr(const DOF_INT_VEC *dv, const char *fn)
1015 {
1016   return file_write_dof_vec_master(true, (const DOF_VEC *)dv,
1017 				   fn, "DOF_INT_VEC     ");
1018 }
1019 
write_dof_int_vec(const DOF_INT_VEC * dv,const char * fn)1020 bool write_dof_int_vec(const DOF_INT_VEC *dv, const char *fn)
1021 {
1022   return file_write_dof_vec_master(false, (const DOF_VEC *)dv,
1023 				   fn, "DOF_INT_VEC     ");
1024 }
1025 
fwrite_dof_int_vec_xdr(const DOF_INT_VEC * dv,FILE * fp)1026 bool fwrite_dof_int_vec_xdr(const DOF_INT_VEC *dv, FILE *fp)
1027 {
1028   return fwrite_dof_vec_master(true, (const DOF_VEC *)dv,
1029 			       fp, "DOF_INT_VEC     ");
1030 }
1031 
fwrite_dof_int_vec(const DOF_INT_VEC * dv,FILE * fp)1032 bool fwrite_dof_int_vec(const DOF_INT_VEC *dv, FILE *fp)
1033 {
1034   return fwrite_dof_vec_master(false, (const DOF_VEC *)dv,
1035 			       fp, "DOF_INT_VEC     ");
1036 }
1037 
write_dof_schar_vec_xdr(const DOF_SCHAR_VEC * dv,const char * fn)1038 bool write_dof_schar_vec_xdr(const DOF_SCHAR_VEC *dv, const char *fn)
1039 {
1040   return file_write_dof_vec_master(true, (const DOF_VEC *)dv,
1041 				   fn, "DOF_SCHAR_VEC   ");
1042 }
1043 
write_dof_schar_vec(const DOF_SCHAR_VEC * dv,const char * fn)1044 bool write_dof_schar_vec(const DOF_SCHAR_VEC *dv, const char *fn)
1045 {
1046   return file_write_dof_vec_master(false, (const DOF_VEC *)dv,
1047 				   fn, "DOF_SCHAR_VEC   ");
1048 }
1049 
fwrite_dof_schar_vec_xdr(const DOF_SCHAR_VEC * dv,FILE * fp)1050 bool fwrite_dof_schar_vec_xdr(const DOF_SCHAR_VEC *dv, FILE *fp)
1051 {
1052   return fwrite_dof_vec_master(true, (const DOF_VEC *)dv,
1053 			       fp, "DOF_SCHAR_VEC   ");
1054 }
1055 
fwrite_dof_schar_vec(const DOF_SCHAR_VEC * dv,FILE * fp)1056 bool fwrite_dof_schar_vec(const DOF_SCHAR_VEC *dv, FILE *fp)
1057 {
1058   return fwrite_dof_vec_master(false, (const DOF_VEC *)dv,
1059 			       fp, "DOF_SCHAR_VEC   ");
1060 }
1061 
write_dof_uchar_vec_xdr(const DOF_UCHAR_VEC * dv,const char * fn)1062 bool write_dof_uchar_vec_xdr(const DOF_UCHAR_VEC *dv, const char *fn)
1063 {
1064   return file_write_dof_vec_master(true, (const DOF_VEC *)dv,
1065 				   fn, "DOF_UCHAR_VEC   ");
1066 }
1067 
write_dof_uchar_vec(const DOF_UCHAR_VEC * dv,const char * fn)1068 bool write_dof_uchar_vec(const DOF_UCHAR_VEC *dv, const char *fn)
1069 {
1070   return file_write_dof_vec_master(false, (const DOF_VEC *)dv,
1071 				   fn, "DOF_UCHAR_VEC   ");
1072 }
1073 
fwrite_dof_uchar_vec_xdr(const DOF_UCHAR_VEC * dv,FILE * fp)1074 bool fwrite_dof_uchar_vec_xdr(const DOF_UCHAR_VEC *dv, FILE *fp)
1075 {
1076   return fwrite_dof_vec_master(true, (const DOF_VEC *)dv,
1077 			       fp, "DOF_UCHAR_VEC   ");
1078 }
1079 
fwrite_dof_uchar_vec(const DOF_UCHAR_VEC * dv,FILE * fp)1080 bool fwrite_dof_uchar_vec(const DOF_UCHAR_VEC *dv, FILE *fp)
1081 {
1082   return fwrite_dof_vec_master(false, (const DOF_VEC *)dv,
1083 			       fp, "DOF_UCHAR_VEC   ");
1084 }
1085 
1086 /*--------------------------------------------------------------------------*/
1087 /*  write_dof_matrix_pbm: print matrix structure as portable bitmap         */
1088 /*--------------------------------------------------------------------------*/
1089 
fwrite_dof_matrix_pbm(const DOF_MATRIX * matrix,FILE * file)1090 bool fwrite_dof_matrix_pbm(const DOF_MATRIX *matrix, FILE *file)
1091 {
1092   FUNCNAME("write_dof_matrix_pbm");
1093   int  i, j, jcol, size, matsize;
1094   MATRIX_ROW *row;
1095   char *pbm_row;
1096 
1097   TEST_EXIT(matrix->type == MATENT_REAL,
1098 	    "Only implemented for scalar matrices so far.\n");
1099 
1100   if (matrix->row_fe_space) {
1101     matsize = matrix->row_fe_space->admin->size_used;
1102   } else {
1103     matsize = matrix->size;
1104   }
1105 
1106   size = matsize + 1;
1107   pbm_row = MEM_CALLOC(size, char);
1108 
1109   fprintf(file, "P1\n");
1110   fprintf(file, "# ALBERTA output of DOF_MATRIX %s\n", matrix->name);
1111   fprintf(file, "%d %d\n", matsize, matsize);
1112 
1113   for (i=0; i < matsize; i++) {
1114     memset(pbm_row, '0', matsize);
1115     for (row = matrix->matrix_row[i]; row; row = row->next) {
1116 
1117       for (j=0; j<ROW_LENGTH; j++) {
1118 	jcol = row->col[j];
1119 
1120 	if (ENTRY_USED(jcol) && row->entry.real[j])
1121 	  pbm_row[jcol] = '1';
1122 
1123       }
1124     }
1125     fprintf(file, "%s\n", pbm_row);
1126   }
1127 
1128   MEM_FREE(pbm_row, size, char);
1129 
1130   return 0;
1131 }
1132 
write_dof_matrix_pbm(const DOF_MATRIX * matrix,const char * filename)1133 bool write_dof_matrix_pbm(const DOF_MATRIX *matrix,
1134 			  const char *filename)
1135 {
1136   FUNCNAME("write_dof_matrix_pbm");
1137   FILE *fp;
1138   bool result;
1139 
1140   if ((fp = fopen(filename,"w")) == NULL) {
1141     ERROR("cannot open file %s\n",filename);
1142     return true;
1143   }
1144 
1145   result = fwrite_dof_matrix_pbm(matrix, fp);
1146 
1147   fclose(fp);
1148 
1149   return result;
1150 }
1151 
1152 /*--------------------------------------------------------------------------*/
1153 
1154 #if 0
1155 /* FORM
1156    LENGTH
1157    AFEM
1158    AMSH
1159    LENGTH
1160    data
1161    AVEC
1162    LENGTH
1163    data
1164    AVEC
1165    LENGTH
1166    data
1167 
1168    where AVEC is one of DRV DRDV DUCV DSCV DINV
1169  */
1170 
1171 /* Write a new IFF-file. On return fp will be positioned at the FORM
1172  * record. If fp is already in IFF format, then just move to the FORM
1173  * record, otherwise generate a new FORM. On entry fp must be
1174  * positioned either at a FORM record which contains an ADTA id.
1175  */
1176 bool fwrite_iff(FILE *fp)
1177 {
1178   unsigned int size = htonl(4);
1179   char tag[4];
1180 
1181   if (fread(tag, 1, 4, fp) != 4 && !feof(fp)) {
1182     return true;
1183   }
1184   if (is_iff_tag(tag, "FORM")) {
1185     if (fseek(fp, 4, SEEK_CUR) != 0) {
1186       return true;
1187     }
1188     if (fread(tag, 1, 4, fp) != 4) {
1189       return true;
1190     }
1191     if (!is_iff_tag(tag, IFF_TAG_ALBERTA)) {
1192       return true;
1193     }
1194     if (fseek(fp, -12, SEEK_CUR) != 0) {
1195       return true;
1196     }
1197     return false;
1198   } else if (is_iff_tag(tag, IFF_TAG_ALBERTA)) {
1199     if (fseek(fp, -12, SEEK_CUR) != 0) {
1200       return true;
1201     }
1202     return fwrite_iff(fp);
1203   }
1204 
1205   /* ok, start a new FORM */
1206   if (fwrite("FORM", 1, 4, fp) != 4) {
1207     return true;
1208   }
1209   if (fwrite(&size, 4, 1, fp) != 1) {
1210     return true;
1211   }
1212   if (fwrite(IFF_TAG_ALBERTA, 1, 4, fp) != 4) {
1213     return true;
1214   }
1215   if (fseek(fp, -12, SEEK_CUR) != 0) {
1216     return true;
1217   }
1218 
1219   return false;
1220 }
1221 
1222 /* Open the named filename which must be an ALBERTA IFF file. Position
1223  * the file pointer on the FORM record. FILENAME will be created when
1224  * it does not exist. If APPEND, do not destroy the old data.
1225  */
1226 FILE *fopen_iff(const char *filename, bool append)
1227 {
1228   FILE *fp;
1229 
1230   if ((fp = fopen(filename, append ? "ab+" : "wb+")) == NULL) {
1231     ERROR("Cannot open file '%s' for writing.\n", filename);
1232     return NULL;
1233   }
1234   if (fseek(fp, 0, SEEK_SET) != 0) {
1235     fclose(fp);
1236     return NULL;
1237   }
1238   if (fwrite_iff(fp) == true) {
1239     fclose(fp);
1240     return NULL;
1241   }
1242   return fp;
1243 }
1244 
1245 /* Correct the size of the enclosing record, which must be a "FORM".
1246  * On entry the file must be positioned at the FORM record.
1247  */
1248 bool fterm_iff(FILE *fp)
1249 {
1250   char tag[4];
1251   unsigned int size;
1252 
1253   if (fread(tag, 1, 4, fp) != 4) {
1254     return true;
1255   }
1256   if (!is_iff_tag(tag, "FORM")) {
1257     return true;
1258   }
1259   if (fread(&size, 4, 1, fp) != 1) {
1260     return true;
1261   }
1262   if (fread(tag, 1, 4, fp) != 4) {
1263     return true;
1264   }
1265   if (!is_iff_tag(tag, IFF_TAG_ALBERTA)) {
1266     return true;
1267   }
1268 
1269   if (memcmp(tag, "FORM", 4) != 0) {
1270     return true;
1271   }
1272   size = htonl(size + 8); /* ADTA header + size */
1273   if (fwrite(&size, 4, 1, fp) != 1) {
1274     return true;
1275   }
1276 
1277   return false;
1278 }
1279 
1280 bool fclose_iff(FILE *fp)
1281 {
1282   bool result = fterm_iff(fp);
1283   fclose(fp);
1284   return result;
1285 }
1286 
1287 /* Write an AMSH record to an IFF file. On entry, the file must point
1288  * to the ADTA record. On return, the file will still be positioned at
1289  * the ADTA record. Obviously, all this can only work with seekable
1290  * files.
1291  */
1292 bool fwrite_mesh_iff(MESH *mesh, REAL time, FILE *fp)
1293 {
1294   long int pos, mshpos, mshend;
1295   char tag[4];
1296   unsigned int size, mshsize;
1297 
1298   if ((pos = ftell(fp)) == -1) {
1299     return true;
1300   }
1301 
1302   if (fread(tag, 1, 4, fp) != 4) {
1303     return true;
1304   }
1305   if (!is_iff_tag(tag, "FORM")) {
1306     return true;
1307   }
1308   if (fread(&size, 4, 1, fp) != 1) {
1309     return true;
1310   }
1311   size = ntohl(size);
1312   if (fread(tag, 1, 4, fp) != 4) {
1313     return true;
1314   }
1315   if (!is_iff_tag(tag, IFF_TAG_ALBERTA)) {
1316     return true;
1317   }
1318   if (fseek(fp, size-4, SEEK_CUR) != 0) {
1319     return true;
1320   }
1321   if ((mshpos = ftell(fp)) == -1) {
1322     return true;
1323   }
1324   if (fwrite(IFF_TAG_MESH, 1, 4, fp) != 4) {
1325     return true;
1326   }
1327   if (fwrite("SIZE", 1, 4, fp) != 4) { /* correct that later */
1328     return true;
1329   }
1330   fwrite_mesh_xdr(mesh, fp, time);
1331   if ((mshend = ftell(fp)) == -1) {
1332     return true;
1333   }
1334   mshsize = mshend - mshpos - 8;
1335 
1336   size += mshsize + 8;
1337   mshsize = htonl(mshsize);
1338 
1339   if (fseek(fp, mshpos+4, SEEK_SET) != 0) {
1340     return true;
1341   }
1342   if (fwrite(&mshsize, 4, 1, fp) != 1) { /* correct size */
1343     return true;
1344   }
1345 
1346   if (fseek(fp, pos+4, SEEK_SET) != 0) {
1347     return true;
1348   }
1349   if (fwrite(&size, 4, 1, fp) != 1) { /* correct size */
1350     return true;
1351   }
1352   if (fseek(fp, pos, SEEK_SET) != 0) {
1353     return true;
1354   }
1355 
1356   return false;
1357 }
1358 
1359 
1360 /* Write an AVEC record to an IFF file. On entry, the file must point
1361  * to the ADTA record. On return, the file will still be positioned at
1362  * the ADTA record. Obviously, all this can only work with seekable
1363  * files.
1364  */
1365 bool fwrite_dof_vec_master_iff(const DOF_VEC *vec,
1366 			       const char *ifftag, const char *type, FILE *fp)
1367 {
1368   long int pos, vecpos, vecend;
1369   char tag[4];
1370   unsigned int size, vecsize;
1371 
1372   if ((pos = ftell(fp)) == -1) {
1373     return true;
1374   }
1375 
1376   if ((pos = ftell(fp)) == -1) {
1377     return true;
1378   }
1379 
1380   if (fread(tag, 1, 4, fp) != 4) {
1381     return true;
1382   }
1383   if (!is_iff_tag(tag, "FORM")) {
1384     return true;
1385   }
1386   if (fread(&size, 4, 1, fp) != 1) {
1387     return true;
1388   }
1389   size = ntohl(size);
1390   if (fread(tag, 1, 4, fp) != 4) {
1391     return true;
1392   }
1393   if (!is_iff_tag(tag, IFF_TAG_ALBERTA)) {
1394     return true;
1395   }
1396   if (fseek(fp, size-4, SEEK_CUR) != 0) {
1397     return true;
1398   }
1399 
1400   if ((vecpos = ftell(fp)) == -1) {
1401     return true;
1402   }
1403   if (fwrite(ifftag, 1, 4, fp) != 4) {
1404     return true;
1405   }
1406   if (fwrite("SIZE", 1, 4, fp) != 4) { /* correct that later */
1407     return true;
1408   }
1409   if (fwrite_dof_vec_master(true, vec, fp, type)) {
1410     return true;
1411   }
1412   if ((vecend = ftell(fp)) == -1) {
1413     return true;
1414   }
1415   vecsize = vecend - vecpos - 8;
1416 
1417   size += vecsize + 8;
1418   vecsize = htonl(vecsize);
1419 
1420   if (fseek(fp, vecpos+4, SEEK_SET) != 0) {
1421     return true;
1422   }
1423   if (fwrite(&vecsize, 4, 1, fp) != 1) { /* correct size */
1424     return true;
1425   }
1426 
1427   if (fseek(fp, pos+4, SEEK_SET) != 0) {
1428     return true;
1429   }
1430   if (fwrite(&size, 4, 1, fp) != 1) { /* correct size */
1431     return true;
1432   }
1433   if (fseek(fp, pos, SEEK_SET) != 0) {
1434     return true;
1435   }
1436 
1437   return false;
1438 }
1439 
1440 bool fwrite_dof_real_vec_iff(const DOF_REAL_VEC *v, FILE *fp)
1441 {
1442   return fwrite_dof_vec_master_iff((const DOF_VEC *)v,
1443 				   IFF_TAG_REAL_VEC, "DOF_REAL_VEC    ", fp);
1444 }
1445 bool fwrite_dof_real_d_vec_iff(const DOF_REAL_D_VEC *v, FILE *fp)
1446 {
1447   return fwrite_dof_vec_master_iff((const DOF_VEC *)v,
1448 				   IFF_TAG_REAL_D_VEC, "DOF_REAL_D_VEC  ", fp);
1449 }
1450 bool fwrite_dof_int_vec_iff(const DOF_INT_VEC *v, FILE *fp)
1451 {
1452   return fwrite_dof_vec_master_iff((const DOF_VEC *)v,
1453 				   IFF_TAG_INT_VEC, "DOF_INT_VEC     ", fp);
1454 }
1455 bool fwrite_dof_schar_vec_iff(const DOF_SCHAR_VEC *v, FILE *fp)
1456 {
1457   return fwrite_dof_vec_master_iff((const DOF_VEC *)v,
1458 				   IFF_TAG_SCHAR_VEC, "DOF_SCHAR_VEC    ", fp);
1459 }
1460 bool fwrite_dof_uchar_vec_iff(const DOF_UCHAR_VEC *v, FILE *fp)
1461 {
1462   return fwrite_dof_vec_master_iff((const DOF_VEC *)v,
1463 				   IFF_TAG_UCHAR_VEC, "DOF_UCHAR_VEC    ", fp);
1464 }
1465 #endif
1466