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_xdr.c */
7 /* */
8 /* description: reading data of mesh and vectors in machine-independent */
9 /* binary format */
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.alberta-fem.de/ */
27 /* */
28 /*--------------------------------------------------------------------------*/
29 /* (c) by A. Schmidt and K.G. Siebert (1996-2004) */
30 /* */
31 /* This program is free software; you can redistribute it and/or modify */
32 /* it under the terms of the GNU General Public License as published by */
33 /* the Free Software Foundation; either version 2 of the License, or */
34 /* any later version. */
35 /* */
36 /* This program is distributed in the hope that it will be useful, */
37 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
38 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
39 /* GNU General Public License for more details. */
40 /*--------------------------------------------------------------------------*/
41 /*--------------------------------------------------------------------------*/
42 /* to_do : number of dofs depending on node */
43 /* nach read_macro, damit dort pointer auf mel, v vorhanden sind!!! */
44 /* (fuers naechste Schreiben) */
45 /* error handling is terrible */
46 /*--------------------------------------------------------------------------*/
47
48 #include <string.h>
49 #include <rpc/types.h>
50 #include <rpc/xdr.h>
51 #include "alberta_intern.h"
52 #include "alberta.h"
53
54 /*
55 XDR-Routinen zur Uebertragung von ALBERTA-Datentypen
56
57 muessen bei abgeaenderter Definition auch veraendert werden !!!
58
59 akt. Def.: REAL = double
60 U_CHAR = unsigned char
61 S_CHAR = signed char
62 DOF = int
63 */
64
65 #include "alberta_intern.h"
66
67 /*--------------------------------------------------------------------------*/
68
69 static DOF_ADMIN *admin = NULL;
70 static MESH *mesh = NULL;
71 static U_CHAR preserve_coarse_dofs;
72
73 static int n_vert_dofs;
74 static DOF **vert_dofs;
75
76 #if DIM_MAX > 1
77 static int n_edge_dofs;
78 static DOF **edge_dofs;
79 #endif
80
81
82 #if 0
83 #if DIM_MAX == 2
84 #define OLD_EDGE CENTER
85 #define OLD_CENTER EDGE
86 #endif
87
88 #if DIM_MAX == 3
89 #define OLD_EDGE CENTER
90 #define OLD_FACE EDGE
91 #define OLD_CENTER FACE
92 #endif
93 #endif
94
95
96 #if DIM_MAX == 3
97 static int n_face_dofs;
98 static DOF **face_dofs;
99 #endif
100
101 static EL *read_el_recursive(EL *el);
102 static void read_dof_admins_xdr(MESH *mesh, bool preserve_coarse);
103
104 /*--------------------------------------------------------------------------*/
105
_AI_read_mesh_1_2(REAL * timeptr)106 MESH *_AI_read_mesh_1_2(REAL *timeptr)
107 {
108 FUNCNAME("read_mesh_1_2");
109 MACRO_EL *mel;
110 int i, j, n;
111 REAL_D *v, x_min, x_max;
112 int neigh_i[N_NEIGH_MAX];
113 char *name, s[sizeof(ALBERTA_VERSION)];
114 int iDIM, iDIM_OF_WORLD, ne, nv;
115 REAL time, diam[DIM_OF_WORLD];
116 int n_vertices, n_elements, n_hier_elements;
117 #if DIM_MAX > 1
118 int n_edges;
119 #endif
120 int *vert_i;
121 int *mel_vertices;
122 static int funccount=0;
123 #if DIM_MAX==3
124 int n_faces, max_edge_neigh;
125 #endif
126
127 _AI_read_int(&iDIM);
128
129 #if 0
130 if (iDIM != DIM_MAX)
131 {
132 ERROR("wrong DIM %d. abort.\n", iDIM);
133 goto error_exit;
134 }
135 #else
136 if (iDIM > DIM_MAX) {
137 ERROR("dim == %d is greater than DIM_MAX == %d!\n", iDIM, DIM_MAX);
138 goto error_exit;
139 }
140 #endif
141
142 _AI_read_int(&iDIM_OF_WORLD);
143 if (iDIM_OF_WORLD != DIM_OF_WORLD)
144 {
145 ERROR("wrong DIM_OF_WORLD %d. abort.\n", iDIM_OF_WORLD);
146 goto error_exit;
147 }
148
149 _AI_read_REAL(&time);
150 if (timeptr) *timeptr = time;
151
152 _AI_read_int(&i); /* length without terminating \0 */
153 if(i) {
154 name = MEM_ALLOC(i+1, char);
155 _AI_read_string(name, i);
156 }
157 else {
158 funccount++;
159 i=100;
160 name = MEM_ALLOC(i+1, char);
161 sprintf(name, "READ_MESH%d", funccount);
162 }
163
164 _AI_read_int(&n_vertices);
165
166 #if DIM_MAX > 1
167 _AI_read_int(&n_edges);
168 #endif
169 _AI_read_int(&n_elements);
170 _AI_read_int(&n_hier_elements);
171
172 #if DIM_MAX == 3
173
174 _AI_read_int(&n_faces);
175 _AI_read_int(&max_edge_neigh);
176
177 #endif
178
179 _AI_read_vector((char *)diam, DIM_OF_WORLD, sizeof(REAL),
180 (xdrproc_t)AI_xdr_REAL);
181 _AI_read_U_CHAR(&preserve_coarse_dofs);
182
183 mesh = GET_MESH(iDIM, name, NULL, NULL, NULL);
184
185 read_dof_admins_xdr(mesh, preserve_coarse_dofs);
186
187 MEM_FREE(name, i+1, char);
188 /* mesh->parametric = NULL; */
189
190 _AI_read_int(&n_vert_dofs);
191
192 if (n_vert_dofs > 0)
193 {
194 vert_dofs = MEM_ALLOC(n_vert_dofs, DOF *);
195 n = mesh->n_dof[VERTEX];
196 for (i = 0; i < n_vert_dofs; i++)
197 {
198 vert_dofs[i] = get_dof(mesh, VERTEX);
199
200 _AI_read_vector((char *)vert_dofs[i], n, sizeof(DOF),
201 (xdrproc_t) AI_xdr_DOF);
202 }
203 }
204
205 #if DIM_MAX > 1
206 _AI_read_int(&n_edge_dofs);
207
208 if (n_edge_dofs > 0)
209 {
210 edge_dofs = MEM_ALLOC(n_edge_dofs, DOF *);
211 n = mesh->n_dof[EDGE];
212 for (i = 0; i < n_edge_dofs; i++)
213 {
214 edge_dofs[i] = get_dof(mesh, EDGE);
215
216 _AI_read_vector((char *)edge_dofs[i], n, sizeof(DOF),
217 (xdrproc_t) AI_xdr_DOF);
218 }
219 }
220 #endif
221
222 #if DIM_MAX==3
223
224 _AI_read_int(&n_face_dofs);
225
226 if (n_face_dofs > 0)
227 {
228 face_dofs = MEM_ALLOC(n_face_dofs, DOF *);
229 n = mesh->n_dof[FACE];
230 for (i = 0; i < n_face_dofs; i++)
231 {
232 face_dofs[i] = get_dof(mesh, FACE);
233
234 _AI_read_vector((char *)face_dofs[i], n, sizeof(DOF),
235 (xdrproc_t) AI_xdr_DOF);
236 }
237 }
238 #endif
239
240 _AI_read_int(&ne);
241 _AI_read_int(&nv);
242
243 ((MESH_MEM_INFO *)mesh->mem_info)->count = nv;
244 v = ((MESH_MEM_INFO *)mesh->mem_info)->coords = MEM_ALLOC(nv, REAL_D);
245
246 for (i = 0; i < nv; i++)
247 _AI_read_vector((char *)v[i], DIM_OF_WORLD, sizeof(REAL),
248 (xdrproc_t) AI_xdr_REAL);
249
250 for (j = 0; j < DIM_OF_WORLD; j++)
251 {
252 x_min[j] = 1.E30;
253 x_max[j] = -1.E30;
254 }
255
256 for (i = 0; i < nv; i++)
257 for (j = 0; j < DIM_OF_WORLD; j++)
258 {
259 x_min[j] = MIN(x_min[j], v[i][j]);
260 x_max[j] = MAX(x_max[j], v[i][j]);
261 }
262
263 for (j = 0; j < DIM_OF_WORLD; j++)
264 mesh->diam[j] = x_max[j] - x_min[j];
265
266
267 mel = MEM_CALLOC(ne, MACRO_EL);
268
269 mesh->n_macro_el = ne;
270 mesh->macro_els = mel;
271
272 vert_i = mel_vertices = MEM_ALLOC(ne*N_VERTICES(iDIM), int);
273
274 for (n = 0; n < ne; n++)
275 {
276 mel[n].index = n;
277 memset(mel[n].neigh_vertices, -1, sizeof(mel[n].neigh_vertices));
278
279 _AI_read_vector((char *)vert_i, N_VERTICES(iDIM), sizeof(int),
280 (xdrproc_t)xdr_int);
281
282 for (i = 0; i < N_VERTICES(iDIM); i++)
283 {
284 if ((*vert_i >= 0) && (*vert_i < nv))
285 mel[n].coord[i] = v + *vert_i++;
286 else
287 mel[n].coord[i] = NULL;
288 }
289
290 /* vertex- and edge_bound are set by _AI_fill_bound_info, therefore we
291 * just read dummies.
292 */
293 #if 1
294 if (iDIM == 1) {
295 _AI_read_vector(
296 (char *)mel[n].wall_bound, N_VERTICES(iDIM), sizeof(S_CHAR),
297 (xdrproc_t) AI_xdr_S_CHAR);
298 #if DIM_MAX > 1
299 } else {
300 S_CHAR dummy[N_VERTICES_MAX];
301 _AI_read_vector(
302 (char *)dummy, N_VERTICES(iDIM), sizeof(S_CHAR),
303 (xdrproc_t) AI_xdr_S_CHAR);
304 #endif
305 }
306
307 #if DIM_MAX >= 2
308 if (iDIM == 2) {
309 _AI_read_vector(
310 (char *)mel[n].wall_bound, N_EDGES(iDIM), sizeof(S_CHAR),
311 (xdrproc_t) AI_xdr_S_CHAR);
312 }
313 #endif
314
315 #if DIM_MAX >= 3
316 if (iDIM == 3) {
317 S_CHAR bound_sc[N_FACES(iDIM)+N_EDGES(iDIM)];
318 _AI_read_vector(bound_sc, (N_FACES(iDIM)+N_EDGES(iDIM)),
319 sizeof(S_CHAR), (xdrproc_t) AI_xdr_S_CHAR);
320 for (i = 0; i < N_FACES(iDIM); i++) {
321 mel[n].wall_bound[i] = bound_sc[i];
322 }
323 }
324 #endif
325 #endif
326 /******************************************************************************/
327
328 _AI_read_vector((char *)neigh_i, N_NEIGH(iDIM), sizeof(int),
329 (xdrproc_t) xdr_int);
330
331 for (i = 0; i < N_NEIGH(iDIM); i++)
332 {
333 if ((neigh_i[i] >= 0) && (neigh_i[i] < ne))
334 mel[n].neigh[i] = mel + (neigh_i[i]);
335 else
336 mel[n].neigh[i] = NULL;
337 }
338
339 _AI_read_vector((char *)mel[n].opp_vertex, N_NEIGH(iDIM), sizeof(U_CHAR),
340 (xdrproc_t) AI_xdr_U_CHAR);
341
342 #if DIM_MAX==3
343 _AI_read_U_CHAR(&(mel[n].el_type));
344
345 {
346 mel[n].orientation = AI_get_orientation_3d(&mel[n]);
347 }
348 #endif
349 mel[n].el = read_el_recursive(NULL);
350 }
351
352 /* Generate the boundary classification for all lower-dimensional
353 * sub-simplices (vertices, in 3d: also edges).
354 */
355 _AI_fill_bound_info(mesh, mel_vertices, nv, ne, false);
356
357 /****************************************************************************/
358 /* The present mechanism only reads DOFs from the file which were used by */
359 /* at least one DOF_ADMIN. The missing element DOF pointers are filled by */
360 /* this routine (in memory.c). */
361 /****************************************************************************/
362 if (iDIM > 0) {
363 AI_fill_missing_dofs(mesh);
364 }
365
366 if (n_elements != mesh->n_elements)
367 {
368 ERROR("n_elements != mesh->n_elements.\n");
369 goto error_exit;
370 }
371
372 if (n_hier_elements != mesh->n_hier_elements)
373 {
374 ERROR("n_hier_elements != mesh->n_hier_elements.\n");
375 goto error_exit;
376 }
377
378
379 if (mesh->n_dof[VERTEX])
380 {
381 if (n_vertices != n_vert_dofs)
382 {
383 ERROR("n_vertices != n_vert_dofs.\n");
384 mesh->n_vertices = n_vert_dofs;
385 goto error_exit;
386 }
387 }
388 mesh->n_vertices = n_vertices;
389
390 #if DIM_MAX > 1
391 mesh->n_edges = n_edges;
392 #endif
393
394 #if DIM_MAX == 3
395 mesh->n_faces = n_faces;
396
397 mesh->max_edge_neigh = max_edge_neigh;
398 #endif
399
400 for (i=0; i<DIM_OF_WORLD; i++)
401 {
402 if (ABS(mesh->diam[i]-diam[i]) > (mesh->diam[i]/10000.0))
403 {
404 ERROR("diam[%i] != mesh->diam[%i].\n",i,i);
405 goto error_exit;
406 }
407 }
408
409 name = s;
410 _AI_read_string(name, sizeof(s)-1);
411 if (strncmp(s, "EOF.", 4)) /* file end marker */
412 {
413 ERROR("no FILE END MARK.\n");
414 goto error_exit;
415 }
416
417 error_exit:
418
419 return mesh;
420 }
421
422 /*--------------------------------------------------------------------------*/
423 /****************************************************************************/
424 /* In ALBERTA_VERSION 1.2 the NODE_TYPES are defined as follows: */
425 /* */
426 /* #if DIM == 1 */
427 /* #define VERTEX 0 */
428 /* #define CENTER 1 */
429 /* #endif */
430 /* */
431 /* #if DIM == 2 */
432 /* #define VERTEX 0 */
433 /* #define EDGE 1 */
434 /* #define CENTER 2 */
435 /* #endif */
436 /* */
437 /* #if DIM == 3 */
438 /* #define VERTEX 0 */
439 /* #define EDGE 1 */
440 /* #define FACE 2 */
441 /* #define CENTER 3 */
442 /* #endif */
443 /* */
444 /* But in our new ALBERTA_VERSION the NODE_TYPES are defined as follows: */
445 /* */
446 /* enum node_types { */
447 /* VERTEX = 0, */
448 /* CENTER, */
449 /* EDGE, */
450 /* FACE, */
451 /* N_NODE_TYPES */
452 /* }; */
453 /* */
454 /* So there is something to be done to "match_node_types" */
455 /* */
456 /****************************************************************************/
_AI_match_node_types(int * node_vec)457 void _AI_match_node_types(int *node_vec)
458 {
459 #if DIM_MAX == 1
460 return;
461 #elif DIM_MAX == 2
462 int tmp = node_vec[CENTER];
463 node_vec[CENTER] = node_vec[EDGE];
464 node_vec[EDGE] = tmp;
465 #elif DIM_MAX == 3
466 int tmp = node_vec[CENTER];
467 node_vec[CENTER] = node_vec[FACE];
468 node_vec[FACE] = node_vec[EDGE];
469 node_vec[EDGE] = tmp;
470 #endif
471 }
472
473 static void
read_dof_admins_xdr(MESH * mesh,bool preserve_coarse_dofs)474 read_dof_admins_xdr(MESH *mesh, bool preserve_coarse_dofs)
475 {
476 FUNCNAME("read_dof_admins_xdr");
477 int i, n_dof_admin, iadmin, used_count;
478
479 int n_dof_el, n_dof[N_NODE_TYPES] = {0,};
480 int n_node_el, node[N_NODE_TYPES] = {0,};
481 int a_n_dof[N_NODE_TYPES] = {0,};
482 char *name;
483 int dim = mesh->dim;
484
485 _AI_read_int(&n_dof_el);
486 _AI_read_vector((char *)n_dof, dim+1, sizeof(int), (xdrproc_t) xdr_int);
487 _AI_match_node_types(n_dof); /* might be correct */
488
489 _AI_read_int(&n_node_el);
490 _AI_read_vector((char *)node, dim+1, sizeof(int), (xdrproc_t) xdr_int);
491 #if 0
492 _AI_match_node_types(node); /* this one stinks */
493 #else
494 /* mesh->node[] is the index of the first vertex, ... etc. DOF. As
495 * the ordering of the nodes has changed, there is more to do than
496 * do a simple reordering of the array. We actually only used to
497 * utilize the stored values for a consistency check. Therefore I
498 * choose to ignore the stored values altogether.
499 */
500 #endif
501 /* use data later for check */
502
503 _AI_read_int(&n_dof_admin);
504 for (iadmin = 0; iadmin < n_dof_admin; iadmin++)
505 {
506 _AI_read_vector((char *)a_n_dof, dim+1, sizeof(int),
507 (xdrproc_t) xdr_int);
508 _AI_match_node_types(a_n_dof); /* probably correct here ... */
509 _AI_read_int(&used_count);
510
511 _AI_read_int(&i); /* length without terminating \0 */
512 name = MEM_ALLOC(i+1, char);
513 _AI_read_string(name, i);
514
515 admin = AI_get_dof_admin(mesh, name, a_n_dof);
516 admin->flags = preserve_coarse_dofs ? ADM_PRESERVE_COARSE_DOFS : 0U;
517 name = NULL;
518
519 if (used_count > 0) enlarge_dof_lists(admin, used_count);
520 } /* end for (iadmin) */
521
522 for (i = 0; i < N_NODE_TYPES; i++) {
523 if (mesh->n_dof[i]) {
524 AI_get_dof_list(mesh, i);
525 }
526 }
527
528 AI_get_dof_ptr_list(mesh);
529
530 TEST(mesh->n_dof_el == n_dof_el,
531 "wrong n_dof_el: %d %d\n", mesh->n_dof_el, n_dof_el);
532 for (i=0; i<=dim; i++)
533 TEST(mesh->n_dof[i] == n_dof[i],
534 "wrong n_dof[%d]: %d %d\n", i, mesh->n_dof[i], n_dof[i]);
535 TEST(mesh->n_node_el == n_node_el,
536 "wrong n_node_el: %d %d\n", mesh->n_node_el, n_node_el);
537 #if 0
538 for (i=0; i<=dim; i++)
539 TEST(mesh->node[i] == node[i],
540 "wrong node[%d]: %d %d\n", i, mesh->node[i], node[i]);
541 #endif
542
543 return;
544 }
545
546 /*--------------------------------------------------------------------------*/
547
548
read_el_recursive(EL * parent)549 static EL *read_el_recursive(EL *parent)
550 {
551 FUNCNAME("read_el_recursive");
552 int i, j, n, node0;
553 EL *el;
554 U_CHAR uc;
555 #if DIM_MAX > 1
556 U_CHAR nc;
557 #endif
558
559 el = get_element(mesh);
560 mesh->n_hier_elements++;
561
562
563 #if ALBERTA_DEBUG
564 el->index = mesh->n_hier_elements;
565 #endif
566
567 _AI_read_U_CHAR(&uc);
568
569 #if DIM_MAX > 1
570 _AI_read_U_CHAR(&nc);
571 if (nc)
572 {
573 el->new_coord = get_real_d(mesh);
574
575 _AI_read_vector((char *)el->new_coord, DIM_OF_WORLD, sizeof(REAL),
576 (xdrproc_t) AI_xdr_REAL);
577 }
578 else
579 {
580 el->new_coord = NULL;
581 }
582 #endif
583
584 if (mesh->n_dof[VERTEX] > 0)
585 {
586 node0 = mesh->node[VERTEX];
587 for (i = 0; i < N_VERTICES(mesh->dim); i++)
588 {
589 _AI_read_int(&j);
590
591 TEST_EXIT(j < n_vert_dofs,
592 "vert_dofs index too large: %d >= %d\n", j, n_vert_dofs);
593 el->dof[node0 + i] = vert_dofs[j];
594 }
595 }
596
597 if ((!uc) || preserve_coarse_dofs)
598 {
599 #if DIM_MAX > 1
600 if (mesh->n_dof[EDGE] > 0)
601 {
602 node0 = mesh->node[EDGE];
603 for (i = 0; i < N_EDGES(mesh->dim); i++)
604 {
605 _AI_read_int(&j);
606
607 TEST_EXIT(j < n_edge_dofs,
608 "edge_dofs index too large: %d >= %d\n", j, n_edge_dofs);
609 el->dof[node0 + i] = edge_dofs[j];
610 }
611 }
612 #endif
613
614 #if DIM_MAX == 3
615 if ((n = mesh->n_dof[FACE]) > 0)
616 {
617 node0 = mesh->node[FACE];
618 for (i = 0; i < N_FACES(mesh->dim); i++)
619 {
620 _AI_read_int(&j);
621
622 TEST_EXIT(j < n_face_dofs,
623 "face_dofs index too large: %d >= %d\n", j, n_face_dofs);
624 el->dof[node0 + i] = face_dofs[j];
625 }
626 }
627 #endif
628
629 if ((n = mesh->n_dof[CENTER]) > 0)
630 {
631 node0 = mesh->node[CENTER];
632 el->dof[node0] = get_dof(mesh, CENTER);
633
634 _AI_read_vector((char *)el->dof[node0], n, sizeof(DOF),
635 (xdrproc_t) AI_xdr_DOF);
636 }
637 }
638
639 #if NEIGH_IN_EL
640 for (i = 0; i < N_NEIGH; i++)
641 {
642 el->neigh[i] = NULL;
643 el->opp_vertex[i] = 0;
644 }
645 #endif
646
647 if (uc)
648 {
649 el->child[0] = read_el_recursive(el);
650 el->child[1] = read_el_recursive(el);
651 }
652 else
653 mesh->n_elements++;
654
655 return(el);
656 }
657